25.2.30 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 online-help for all menus.

In this implementation of ab initio MD, time is divided into steps of equal duration ${\Delta}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 ${\Delta}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 ${\Delta}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) ${\Delta}t$;
atomic symbols and velocities x,y,z (au) at $t-({\Delta}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.7).

$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 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 2-10 times larger than the time step. Note that user-defined actions are presently not supported in canonical dynamics mode.

These are optional keywords:

$seed: 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 () are specified in $ke_control.

$barrier angstroms
  type        elps
  limits      5.0 10.0 7.5
  constant    2.0
  springlen   1.0
  temperature 300.0

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

type: can be one of , , or , 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, semi-major 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.0 H/Bohr/Ångstrøm, a bastard combination of units.
springlen: 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 quasi-temperature 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 of 0.0 K 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, , 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 +/- of the specified constraint. Additionally, the RMS deviation of all constrained distances must be below 2/3 of .
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 atom to each atom that is closer than and apply . The first line above examines each atom and looks for any atoms within 1.2 Å. The shortest distance, if any, is then fixed at 0.9 Å. Similarly, the second line binds each to the closest within 1.7 Å, but since =0.0, that distance is fixed at the current value. The third line attaches atoms to the appropriate nearby , but at the current average distance multiplied by the absolute value of .

Explicitly specified constraints are listed by atom index and supercede auto-generated 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.

User-defined 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.0 a.u., velocities are to be scaled every step to keep average total energy constant. Then, from 2500.0 a.u., gradual cooling at the default rate (annealing) is to occur until the time 3000.0 a.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. quasi-temperature) 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 minimum-energy structures. Temperature may be lowered gradually by a small factor each step (; default factor 0.905 over 100 steps) or lowered rapidly by reversing all uphill motion (; default factor -0.8 each step). The cooling factors may be changed from the default using . Another option allows the quenching part of the run to be logged to a separate file. Alternatively, a standard non-dynamical geometry optimization can be carried out in a subdirectory ().

$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 Tully-type fewest switches surface hopping (SH) [422]. This option is only available in combination with TDDFT. For the TDDFT surface hopping see Tapavicza et al. 2007 [423]; for the current implementation see Tapavicza et al. 2011 [424]. 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:

$currentstate
  1

$currentstate

keyword needed to ensure dynamics starting in S$_1$

$nacme

needed to compute non-adiabatic couplings; this keyword requires the use of weight derivatives in section dft

The excitations, excited state gradients and nonadiabatic coupling vectors are generated from a TDDFT calculation and stores by egrad. The gradients and the nonadiabatic vectors are stored under the data group and the excitations under . These data groups are automatically generated in an MD simulation.

$sh_coeffs file=sh_coeffs

collects amplitudes of the adiabatic states along the trajectory . This keyword is automatically generated in the first MD step, with the entire electronic population on the current active state. However, the user may specify the initial coefficients to be something different by manually creating the file.

$population_el

stores the electronic populations at each time step. This is generated by frog and the user need not specify.

$prob_fssh

stores the probability of hopping at each time step. This is generated by frog and the user need not specify.

$active

stores the active state for surface hopping at each time step. This is generated by frog and the user need not specify if is specified.

Special caution has to be taken to control problems related to conical intersections [425, 426]. 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 Tamm-Dancoff approximation (TDA) to TDDFT (see 8). 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$ TDDFT-TDA 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.

25.2.30 Keywords for Module frog

25.2.30 Keywords for Module `frog`