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.

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

$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, 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.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 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

$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 normalfont 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:

$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. 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 (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 optimisation 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.