Example (XMOLECULE): Liquid water (embedded QM/MM)

This examples propagates an ionized water cluster using the xmolecule electronic structure toolkit and the QM/MM embedding setup. The dynamics is modelled using fewest switches surface hopping and the electronic structure is calculated using Koopanns’ theorem.

See:

Examples/xpyder/embedding_xmol/

Preparation

Make sure you have installed XMOLECULE. You can find an output file for an equilibrated box of water molecule in the file water_prdc.gro. It has been generated with gromacs using the script gromacsWater.sh.

Based on the water structure, we generate initial conditions for the CDTK propagation. To that end, we run the following commands:

mkdir trj_gs
cp input_gs_embedding trj_gs
cd trj_gs
gmx2cdtk -g ../water_sample.gro -t ../water_sample.trr -n 1  --centerMethod first -s 18.9

This script cuts out water cluster with a radius of 18.9 a.u. around the first water (-n 1) molecule and prepares input files for the cdtk simulations. The initial atomic positions and velocities are taken from the files water_sample.gro and water_sample.trr

Details of the Inputfile input_gs_embedding

The file input_gs_embedding is the input file used for the calculation for an ionized water cluster:

$SYSTEM
   qchemistry = embedding
   xunit = bohr
$END
$embedding
   quantum = xmolecule
   region = 15
   forcefield = SPC/E
$end
$xmolecule
   HF=yes
   gto=6-31+G
   N=1
$end

$trajectory
   logn = 10
   dt = 0.5   fs
   tf = 80.0   fs
$end
$quantum
   type = gs
$end

$SYSTEM - qchemistry = xmolecule indicates that xmolecule is used for quantum chemistry engine.
$xmolecule
- gto = 6-31+G specifies the Gaussian basis set.
- HF = yes do Hartree-Fock calculations
- N = 1 sets the number of radial grid points to a minimum (the grid is not needed when employing GTO and HF)
$embedding
- quantum = xmolecule
- region = 15 specifies that the first 15 atoms are treated as QM region
- forcefield = SPC/E specifies the MM force field that is used for water
$trajectory
- logn = 10 write output files every 10th iteration. This reduces the io load.
- “dt” is the time step used for the calculation (in fs).
- “tf” is the time of the last time step (in fs).
quantum
- type = gs propagate on the ground state

Running ground-state dynamics

To run the first example, execute:

xpyder -i input_gs_embedding -d .

in the directory trj_gs.

Output data for ground-state dynamics

The folder trj_gs contains output files.

R.log contains the position of the atoms for each time step.
V.log contains the velocity of the atoms.
E.log contains the current potential energy, kinetic energy and total energy of the trajectory.
partial.log shows partial charges (Mulliken charges and MM charges) for different time steps.

Details of the inputfile for ionized-state dynamics

The file input_hole_embedding is the input file used for the calculation for an ionized water cluster. Compared to the file input_gs_embedding for the neutral ground state calculations, the following sections in the input file are modified:

$xmolecule
   HF=yes
   gto=6-31+G
   gs_occ=yes
   charge = 1
   N=1
$end

$quantum
   type = fssh
   nstates = 20
   istate = 2
   rescaling = nac
$end

The options charge = 1 and gs_occ = yes indicate to consider ionized configurations using Koopmanns’ theorem. The options type = fssh triggers fewest switches surface hopping. The option nstates = 20 gives the number of states considered. It MUST match the number of valence orbitals (here 5 water molecules = 20 valence orbitals) The option istate = 2 gives the index of starting state: this means that the hole is initially in HOMO-2. The option rescaling = nac specifies how rescaling of the velocity is done during a surface hop. “nac” means rescale velocities along the coupling vectors.

Running ionized-state dynamics

To run the first example, execute:

mkdir trj_hole
cp input_hole_embedding trj_hole
cd trj_hole
gmx2cdtk -g ../water_sample.gro -t ../water_sample.trr -n 1 --centerMethod first -s 18.9
xpyder -i input_hole_embedding -d .

Output data for ionized-state dynamics

The folder trj_hole contains the following output files.

R.log contains the position of the atoms for each time step.
V.log contains the velocity of the atoms.
E.log contains the current potential energy, kinetic energy and total energy of the trajectory.
Switch.log contains information about surface hops
S.log contains information about the current state index
C.log state coefficients
P.log state populations
V_ad.log adiabatic potential energies as a function of time

Details of the inputfile for excited-state dynamics

The file input_cis_embedding addresses the dynamics of a water cluster in the excited state. The excited state is described using configuration interaction singles. The following sections are modified compared to the earlier input files:

$xmolecule
  HF = yes
  CIS = yes
  nstates = 10
  gto = 6-31+G
$end

$quantum
   type = fssh
   nstates = 10
   istate  = 3
   rescaling = nac
   N = 1
$end

In the section $xmolecule, the line CIS = yes switches on configuration interaction singles calculation, nstates = 10 indicates that up to 10 states are calculated. In the section $quantum, the parameter “nstates” was set to 10 (the same number as in the section quantum), and “istate” is set to 3 (the initial state), which means the 3rd excited state (0 means ground state).

Running the excited-state trajectory

To run the second example, execute:

mkdir trj_cis
cp input_cis_embedding trj_cis
cd trj_cis
gmx2cdtk -g ../water_sample.gro -t ../water_sample.trr -n 1 --centerMethod first -s 18.9
xpyder -i input_cis_embedding -d .

Output data for the excited-state trajectory

The folder cis_embedding contains the same output files as for the earlier run for the ionized state.

Electrostatic embedding

All previous examples can be combined with electrostatic embedding. To do so, the additional option electrostatic = yes has to be added to the $embedding section. Furthermore, the option qmmm = yes has to be added to the $xmolecule section.