The Molecular Dynamics Module

Introduction

The MD module within ChemShell is based on routines from the DL_POLY package (version 2.11) and supports a range of direct dynamics integration algorithms which may be coupled with any ChemShell module capable of providing energy and gradients. The integration in performed within the constant NVE ensemble, and uses the quaternion and shake algorithms when rigid bodies and interatomic distance constraints are present. Periodic boundary conditions are supported.

The Tcl interface to the module is based on an object-oriented scheme, similar to that used (for example) by the Tk toolkit. A command, in this case dynamics is invoked once, creating the MD object. At this stage the internal memory required is allocated, any book-keeping is performed and the module which will be used to compute energies and gradients is initialised.

Once created, the name of the object created becomes a Tcl command which can be used to invoke the a set of methods. At the present time, no specific simulation protocols are explicitly programmed, the required simulation protocol is provided as part of the users script, the relevant methods being invoked as required.

The `dynamics` Command

Command Line Arguments

Argument Argument type Mandatory Default To specify

theory= string no mopac module used to compute energy and forces

coords= Fragment tag yes none initial configuration of the system

timestep= real no 0.0005 integration timestep (ps)

temperature= real no 293 simulation temperature (K)

mcstep= real no 0.2 Max step displacement (a.u.) for Monte Carlo

rigid_groups= Tcl List no " " specification of rigid group (quaternion definitions)

constraints= Tcl List no " " specification interatomic distances for SHAKE

shake_h= keyword no none SHAKE constraints for X-H distances ("input" retains the lengths from the input structure, "ideal" uses ideal values).

ensemble= keyword no NVE Choice of ensemble (NVE, NVT, NPT, NST, MC or FRIC)
FRIC performs friction dynamics. The (constant) friction is specified by the friction argument.

friction= real no 0.0 Sets the (constant) friction parameter for [-1.0 ... +1.0]. Use with ensemble FRIC.

frozen= Tcl List of integers and ranges no none List of frozen atoms

masses= Tcl List no none List of masses that differ from standard values. A nested Tcl List : each entry being a list that contains the number of the atom and the mass of the atom.

nosehoover= real no 0 0: If specified with a T=const. ensemble, uses the standard Berendsen thermostat.
1: If specified with a T=const. ensemble, uses a simple Nose-Hoover thermostat.
The characteristic time of the thermostat is specified by the taut argument. If the Nose-Hoover thermostat is enabled, the method "friction" is made available within the MD object. The total energy includes thermostat contributions.
>1: Nose-Hoover chain thermostat with the specified chain length. The characteristic time for the 1st thermostat, which couples to the system, is given by taut. The other thermostats have a characteristic time of taut/(degrees of freedom).

taut= real no 0.5 Tau(t) for Berendsen Thermostat (ps)

taup= real no 5.0 Tau(p) for Berendsen Barostat (ps)

compute_pressure= Boolean no no Whether to compute pressure and virial (for NVT simulation)

trajectory_type= integer no 0 Specification of data written to file dynamics.trj
0: Positions (Bohr)
1: Positions and Velocities (Ang/ps)
2: Positions, Velocities and Forces
10,11,12: as above, but frozen atoms are not printed.

trajectory_file= string no dynamics.trj file for trajectory output

list_option= Output keyword no medium Provide additional output quantities (none, medium, full = debug).

energy_unit= keyword no au Basic energy unit for energy output (au, kcal mol-1, ev).

reaction_coord= Tcl List no " " reaction coordinate= { type atoms value }
type: Type of constraint,
Type "1" is bond length, d(R1,R2). Two atoms, R1, R2 must be given. Value is the value of the bond length in a.u. (Bohr).
Type "2" for the difference of bond length, d(R1,R2)-d(R2,R3) atoms: atom numbers involved in the constraint. Three atoms, R1, R2, R3 must be given. Value: value of the constraint in a.u. (Bohr).
Type "3" is a torsion angle. Four atoms must be given. Val is the value in rad.
The resulting force may be obtained by get constraint_force in atomic untis (Hartree/Bohr or Hartree/rad, respectively). Depending on energy_unit the unit of the force changes.

Ensemble specification

The available ensembles are NVE (constant energy), NVT (constant volume and temperature) and NPT(constant pressure and temperature, with isotropic cell shape changes), NST (constant pressure and temperature, with anisotropic cell shape changes) and MC (monte carlo).

Rigid Group specification

The rigid body definitions are provided as a nested Tcl List . Each element of the main list defines a rigid group, and takes the form of a list of integers, specifying the atoms in that group.

Constraint specification

The constraints input for the SHAKE procedure is provided similarly as a list of lists, each sub-list in this case should contain 3 elements: the indices of the two atoms involved and the constrained bond length (in a.u.).

Pairlist Updates

If you are using a QM/MM or MM energy expression, you will probably need to do a periodic update of the pairlist. This can be requested using the update method of the dynamics object. The frequency of this will depend on the nature of your system. For the special case where you have requested to use all pairs, or to skip the pairlist, you may not need to run the update, but for periodic systems with an MM calculation present, it is essential for correct simulation. In the QM/MM case, if you are using a QM/MM pairlist this will be updated at the same time as the MM pairlist (only one update call is needed).

Methods available within the MD object

configure - modify simulation parameters
initvel - initialise random velocities
force - evaluate molecular forces
step - Take MD step
update - Request MM or QM/MM pairlist update
mctest - Test step (Monte Carlo only)
output - print data (debugging use only)
printe - print step number, kinetic, potential, total energy, temperature, pressure volume, virial, and friction.
get - Return a variable from the dynamics, takes a single keyword, which can be temperature, input_temperature, pressure, input_pressure, total_energy, kinetic_energy, potential_energy, friction, constraint_force or time
trajectory - Output the current configuration to the trajectory file
destroy - free memory and destroy object
dump - write restart files
load - read restart files
fcap - limit maximum force (during warm-up phase)

The configure method takes additional arguments of the same type as the dynamics command itself, except that only a subset of arguments are accepted. Options currently available include temperature pressure timestep taup taut ensemble energy_unit friction and trajectory_file. You cannot, for example, change the definition of constraints or the theory specification. It is used as follows

#
dynamics dyn1 coords=c \
   theory=dl_poly : { mm_defs=ff } \
   temperature = 300 \
   timestep = 0.001 
#
#  Some MD steps here
#

dyn1 configure temperature = 400 ensemble=NVT

#
# Second phase of MD here
#

Examples

A Box of Water

#
fragment c persistent
#
c_create coords=c {
coordinates angstrom
O     .6737543037        -2.533871039         4.341299547    
H     .3135705220        -2.431626305         3.414038166    
H     .9946264752E-01    -3.172722493         4.853215483    
# 99 other water molecules here 
cell angstrom
       15.0267420457         .0000000000         .0000000000
         .0000000000       15.0267420457         .0000000000
         .0000000000         .0000000000       15.0267420457
}
#
# All H-H's must be connected so that the topology generator
# adds in the bond terms
#
for {set i 1} {$i <= 300} {incr i 3} {
  set h1 [ expr $i + 1 ]
  set h2 [ expr $i + 2 ]
  add_connection coords=c i = $h1 j = $h2
}
#
# provide force-field parameters. For vdw use the 6-12 expression
# input in kcal/mol
#
read_input ff {
vdw   o o    628.841   633295.2
vdw   o h   .000000      0.000
vdw   h h   .000000      0.000
bond  o h   100.0        1.000
bond  h h   1000.0       1.63288
charge o  -0.82
charge h   0.41
}
#
# Create the dynamics object
#
dynamics dyn1 coords=c \
   theory=dl_poly : { mm_defs=ff } \
   temperature = 300 \
   timestep = 0.001 

# initialise random velocities
dyn1 initvel
#
set count 0
while {$count < 100} {

    # Update MM pairlist every 10 steps
    if { [ expr $count % 10 ] == 0 } { dyn1 update }
    dyn1 force

    # output the coordinates to the trajectory every 50 steps
    if { [ expr $count % 50 ] == 0 } { dyn1 trajectory }

    dyn1 printe
    dyn1 step
    incr count
}
dyn1 destroy
delete_object c
#
times
#

Using quaternions rigid water molecules

To employ rigid water molecules, the previous input can be modified by replacing the dynamics invocation as follows:

for {set i 1} {$i <= 300} {incr i 3} { 
  lappend rgd [list $i [expr $i + 1] [ expr $i + 2] ]
}
dynamics dyn1 coords=c \
   theory=dl_poly : { mm_defs=ff } \
   temperature = 300 \
   timestep = 0.001 \
   rigid_groups = $rgd

Using SHAKE constraints