Next: Keywords for Module Mpshift
Up: Format of Keywords and
Previous: Keywords for wave function
Contents
Index
Keywords for Module FROG
The ab initio molecular dynamics
(MD) program frog needs a command file
named mdmaster.
The interactive Mdprep
program manages the generation of
mdmaster and associated files. It is always a good idea to
let Mdprep check over mdmaster before starting
an MD run. Mdprep has onlinehelp for all menus.
In this implementation of ab initio MD, time
is divided into steps of equal duration Δt.
Every step, the energy and its gradient are calculated and
these are used by the frog to work out the new
coordinates for the next step along the dynamical trajectory.
Both the accuracy of the trajectory and the total computation
time thus depend crucially on the time step
chosen in Mdprep.
A bad choice of timestep will result in integration errors
and cause fluctuations and drift in the total energy.
As a general rule of thumb, a timestep Δt should be chosen
which is no longer than one tenth of the shortest
vibrational period of the system to be simulated.
Note that Mdprep will transform velocities so that the
total linear and angular momentum is zero.
(Actually, for the Leapfrog algorithm,
initial velocities are
Δt/2 before the starting time).
The following keywords are vital for frog:
 $nsteps 75

Number of MD time steps to be carried out.
$nsteps
is decreased by 1 every time frog is run
and JOBEX md stops when $nsteps
reaches 0.
 $natoms 9

Number of atoms in system.
 $current file=mdlog.aa

The file containing the current position, velocity,
time and timestep,
that is, the input configuration.
During an MD run the $current
information
is generally kept at the end of the $log
file.
 $log file=mdlog.ZZ

The file to which the trajectory should be logged,
i.e. the output:
t=time (a.u.);
atomic positions x,y,z (Bohr) and symbols at t;
timestep (au) Δt;
atomic symbols and velocities x,y,z (au) at
t  (Δt/2);
kinetic energy (H) interpolated at t, ab initio
potential
energy (H) calculated at t, and pressure recorded at
the
barrier surface (atomic units, 1 au = 29.421 TPa) during
the
corresponding timestep;
ab initio potential energy gradients x,y,z
(H/Bohr) at t.
This file can be
manipulated with LOG2? tools after the MD run
(Section 1.5).
 $turbomole file=control

Where to look for TURBOMOLE keywords $grad
etc.
 $md_status

The status of the MD run is a record of the action
carried out during the previous MD step,
along with the duration of that step.
The format matches that of $md_action
below.
Canonical dynamics is supported using the NoséHoover
thermostat. This option can be enabled in Mdprep or
by the following syntax:
$md_status
canonical T=500 t=100
from t= 25.0000000000 until t= 0.00000000000
Here, T specifies the temperature of the thermostat in K (500 K in the
example) and t specifies the thermostat relaxation time in a.u. (100
a.u. in the example). It is advisable to choose the thermostat
relaxation 210 times larger than the time step. Note that
userdefined actions are presently not supported in canonical
dynamics mode.
These are optional keywords:
 $seed 123

Integer random number seed
 $title

Arbitrary title
$log_history
100 mdlog.P
71 mdlog.Q
$ke_control
length 50
response 1
To determine the
trends in kinetic energy and total energy (average values and
overall drifts) it is necessary to read the history of energy
statistics over the recent MD steps. The number of MD steps
recorded so far in each log file are therefore kept in the
$log_history entry: this is updated by the program each
step. The length of records needed for reliable statistics and the
number of steps over which changes are made to kinetic energy
(response) are specified in $ke_control.
$barrier angstroms
type elps
limits 5.0 10.0 7.5
constant 2.0
thickness 1.0
temperature 300.0

 $barrier specifies a virtual cavity for simulating condensed
phases. The optional
flag, angstroms, can be used to indicate that
data will be entered in Ångstrøms rather than Bohr.
 type

can be one of orth, elps,
or none, for orthorhombic, ellipsoidal, or
no barrier (the default) respectively.
 limits

are the +x,y,z sizes of the cavity. In this
case, an ellipsoid with a major axis of 20Å along y, semimajor of
15Å on z, and minor of 10Å on x.
 constant

is the Hooke's Law force constant in atomic
units of force (H/Bohr) per length unit. Here, it is
2.0H/Bohr/Ångstrøm,
a bastard combination of units.
 thickness

is the effective limit to the restorative force
of the barrier. For this system, an atom at 5Å into the
barrier will feel the same force as at 1.0Å.
 temperature

denotes the temperature of the cavity
walls in Kelvin. If the system quasitemperature is below this
setpoint, particles will be accelerated on their return to the
interior. Alternately, they will be retarded if the system is too
warm. A temperature of 0.0K will turn off wall
temperature control, returning molecules to the system with the same
momentum as when they encountered the barrier.
$constraints angstroms
tolerance 0.05
adjpercyc 0.25
type H O 0.9 1.2
type F C 0.0 1.7
type H C 1.0 1.2
2 1 0.0
3 1 1.54
4 1 1.0
 $constraints

specifies and/or automatically
generates atomic distance constraints. The optional
flag, angstroms, can be used to indicate that
data will be entered in Ångstrøms rather than Bohr.
 tolerance

is the convergence criterion for application
of constraints. All distances must be within +/
tolerance of the specified constraint. Additionally,
the RMS deviation of all constrained distances must be below 2/3 of
tolerance.
 adjpercyc

is the fraction of the total distance
correction to be applied on each constraint iteration.
 type X A const rmax

commands frog to
find the closest A atom to each atom X
that is closer than rmax and
apply const. The first type
line above examines each H atom and looks for any
O atoms within 1.2Å. The shortest distance,
if any, is then fixed at 0.9Å. Similarly, the second type
line binds each F to the closest C
within 1.7Å, but since const=0.0, that distance is fixed
at the current value. The third type line attaches H
atoms to the appropriate nearby C, but at the current
average HC distance multiplied by the absolute value of
const.
Explicitly specified constraints are listed by atom index and
supercede autogenerated constraints. A positive third number fixes
the constraint at that value, while zero fixes the constraint at the
current distance, and a negative number unsets the constraint.
The output of frog contains the full list of constrained atom pairs
and their current constraints in explicit format.
Userdefined instructions allow the user to tell frog to change some aspect
of the MD run at some point in time t=real number.
The same format is used for $md_status above.
Here is an example:
$md_action
fix total energy from t=2000.0
anneal from t=2500.0
free from t=3000.0
In this example, starting from the time 2000.0a.u.,
velocities are to be scaled every step to keep average total energy constant.
Then, from 2500.0a.u., gradual cooling at the default rate (annealing) is to
occur
until the time 3000.0a.u., when free Newtonian dynamics will resume.
Here are all the possible instructions:
$md_action
fix temperature from t=<real>
fix total energy from t=<real>


These commands cause velocities
to be scaled so as to
keep the average kinetic energy
(i.e. quasitemperature)
or the average total energy approximately constant.
This is only possible once enough information about run
history is available to give reliable statistics.
(Keywords $log_history, $ke_control).
$md_action
set temperature at t=<real> to x=<real> K
set total energy at t=<real> to x=<real> H
set kinetic energy at t=<real> to x=<real> H
set position file=<filename> at t=<real>
set velocity file=<filename> at t=<real>
set velocity at t=<real> random
set velocity at t=<real> zero


At some time during the ab initio MD run the user can specify a
new value for one of the dynamical variables.
The old value is discarded.
Single values are given by x=real number.
Vectors must be read in frog format from file=file.
$md_action
anneal from t=<real>
anneal from t=<real> x=<real>
quench from t=<real>
quench from t=<real> x=<real> file=<file>
relax at t=<real>


In Simulated Annealing MD,
the temperature of a run
is lowered so as to find minimumenergy structures.
Temperature may be lowered gradually by a small factor
each step (anneal; default factor 0.905 over 100 steps)
or lowered rapidly by reversing all uphill motion
(quench; default factor 0.8 each step).
The cooling factors may be changed from the default using x=.
Another option allows the quenching part of the run to be logged to a
separate file.
Alternatively, a standard nondynamical geometry
optimization can be carried out in a subdirectory (relax).
$md_action
free from t=<real>


Finally, this instruction turns off any previous action
and resumes free dynamics.
This is the default status of an MD run.
$surface_hopping


This keyword allows to carry out Tullytype fewest switches surface hopping (SH) [159].
This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [160]; for the current implementation see Tapavicza et al. 2011 [161].
In the current implementation the surface hopping algorithm only allows switches between the first excited singlet state and the ground state. However, total energies of higher excited states can be computed during the MD simulation.
The proper functioning of SH has only been tested for the option
$md_action
fix total energy from t= 0.00000000000
To carry out SH dynamics simulations, the keyword $surface_hopping
has to be added
to the control and mdmaster file.
In addition several keywords are required in the control file:
 $nacme

needed to compute nonadiabatic couplings; this keyword requires the use of weight derivatives
in section dft
 $nac

keyword needed to collect Cartesian nonadiabatic coupling vectors along the trajectory
 $exopt 1

keyword needed to ensure dynamics starting in S_{1}
 $ex_energies file=ex_energies

collects excitation energies along the trajectory
 $integral_ex file=integral_ex

collects time integration of excitation energies along the trajectory
 $sh_coeffs file=sh_coeffs

collects amplitudes of the adiabatic states along the trajectory
 $nac_matrix file=nac_matrix

collects NAC elements along the trajectory
Special caution has to be taken to control problems related to conical intersections [162,163].
At geometries where conical intersections between the ground and excited state are present, DFT often exhibits singlet instabilities, which leads to imaginary excitation energies in linear response TDDFT; in this case the MD run is terminated.
This problem can be circumvented by the use of the TammDancoff approximation (TDA) to TDDFT (see 7).
In addition an optional keyword for the md_master file can be used:
 $gap_threshold <real>

enforces a switch to the ground state in case the S_{1}S_{0} energy gap drops below <real>
eV. As default a switch to S_{0} is enforced if the S_{1} TDDFTTDA excitation energy becomes negative.
Often times if a switch is enforced due to a negative TDA excitation energy the potential energy surface is discontinuous and limited numerical precision of the nuclear forces may lead to a loss of total energy conservation. In this case the nuclear velocities are rescaled to obtain a conserved total energy.
Next: Keywords for Module Mpshift
Up: Format of Keywords and
Previous: Keywords for wave function
Contents
Index
TURBOMOLE