next up previous contents index
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 online-help 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.

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:

  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 -123

Integer random number seed

Arbitrary title
 100                mdlog.P
 71                 mdlog.Q
  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.

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

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.

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.

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Å.

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 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

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.

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.

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 H-C distance multiplied by the absolute value of const.

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:

     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:

     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).
     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.
     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 (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 non-dynamical geometry optimization can be carried out in a subdirectory (relax).
     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.


This keyword allows to carry out Tully-type 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
 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:

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

keyword needed to collect Cartesian non-adiabatic coupling vectors along the trajectory
$exopt 1

keyword needed to ensure dynamics starting in S1
$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 Tamm-Dancoff 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 S1-S0 energy gap drops below <real> eV. As default a switch to S0 is enforced if the S1 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.

next up previous contents index
Next: Keywords for Module Mpshift Up: Format of Keywords and Previous: Keywords for wave function   Contents   Index