Dynamics: Description and Discussion There are four separate dynamics integrators available in CHARMM: (This discussion does not apply to multi-body dynamics, which has a separate set of integrators). See *note Mbond:(mbond.doc). Name Keyword Module Original Verlet ORIG dynamcv.src Leapfrog Verlet LEAP dynamc.src (default) Velocity Verlet VVER dynamvv.src 4-D L-F Verlet VER4 dynam4.src All methods are based on the Verlet scheme, and when used without any special features, provide identical trajectories for short simulations. All methods allow SHAKE. The ORIG integrator is a standard 3-step Verlet integrator with few frills. It allows: Langevin Dynamics (LANG) Thermodynamic Simulation Method (TSM) The LEAP integrator is similar to the ORIG integrator, but does provide increased accuracy (esp. for single precision version of CHARMM). It allows: Langevin dynamics (LANG) (with accurate temperatures printed) Constant Temperature and Pressure (CPT) (based on Berendsen's method) Accurate pressures with SHAKE High frequency correction to the total energy Parallel code Free energy equilibration indicator (deltaF*V) (with PERT) Thermodynamic Simulation Method (TSM) The VVER integrator also provides increase accuracy. It allows: Constant Temperature (NOSE) (Nose-Hoover method) Multiple Time Step (MTS) The VER4 integrator enables the energy embedding technique that entails placing a molecule into a higher spatial dimension [Crippen, G. M. & Havel,T.F. (1990) J.Chem.Inf.Comput.Sci. Vol 30, 222-227]. The possibility of surmounting energy barriers with these added degrees of freedom may lead to lower energy minima. Here, this is accomplished by molecular dynamics in four dimensions. Specifically, another cartesian coordinates was added to the usual X, Y, and Z coordinates in the LEAPfrog VERLet algorithm. In order to generate a dynamics trajectory, all requirements for evaluating the energy must be met. See *note Energy:(energy.doc)Needs. * Menu: * Syntax:: Syntax of the dynamics command * Description:: Description of the keywords and options * Recommended:: Recommended input options and values * Discussion:: Running dynamics * Output:: Output from a dynamics run * Trajectory:: Trajectory manipulation and I/O * Merge:: Merging or breaking up trajectory files into different size pieces. Resampling at a larger interval. Least squares fit to a reference. Recentering, or undoing image operations. * Reorient:: Reorienting a coordinate trajectory * RMSDyn:: Computes the RMS difference between two trajectories * Format:: formatting and unformatting a dynamics trajectory * CVELocity:: Constant velocity dynamics * TMD:: Targeted molecular dynamics * TSALlis:: Molecular dynamics in the Tsallis ensemble * CENT:: The reCENTering command * Monitor:(monitor.doc). Monitor dihedral transitions * CPT dynamics:(pressure.doc). CPT dynamics * 4-D dynamics:(fourd.doc). 4-D dynamics * Pressure:(pressure.doc)Pressure. The pressure command * MTS:(mts.doc). Multiple Time Scales Method * Nose:(nose.doc). Nose-Hoover Dynamics * MBOND:(mbond.doc)Dynamic. Multi-body Dynamics

Syntax for the Dynamics Command DYNAmics { [ LEAP ] [ VERLet ] } ! Dynamics with the leap-frog integrator { [ NEW ] [ LANGevin ] } { [ CPT ] [ NOLAngevin ] } other-specs cpt-spec { [ EULEr ] } { ORIG } ! Dynamics with the verlet integrator DYNAmics { VVERlet [ NOSE ] } ! Dynamics with the velocity verlet integrator { } other-specs DYNAmics { LEAPfrog VER4 four-dim-spec } ! Four dimensional dynamics { } other-specs DYNAmics { MBONd mbond-spec } other-specs ! Multibody dynamics other-specs::= [NSTEp integer] [nonbond-spec] [hbond-spec] [frequency-spec] {[TIMEstp real]} {STARt } [unit-spec] [temperature-spec] [options-spec] { AKMA real } {RESTart} hbond-spec::= see *note Hbonds:(hbonds.doc). nonbond-spec::= See *note Nbonds:(nbonds.doc). mbond-spec::= See *note Mbond:(mbond.doc)DynDesc. frequency-spec::= [INBFrq integer] [IEQFrq integer] [IHBFrq integer] [IHTFrq integer] [IPRFrq integer] [NPRInt integer] [NSAVC integer] [NSAVV integer] [NTRFrq integer] [ILBFrq integer] [ISVFRQ integer] [NSAVL integer] [IMGFrq integer] [IXTFrq integer] [NLAT integer] [NBLCkfep integer] unit-spec::= [IUNCrd integer] [IUNRea integer] [IUNVel integer] [IUNWri integer] [KUNIt integer] [CRAShu integer] [BACKup integer] [IUNLdm integer] [ILAT integer] [IBLCkfep integer] [ILAP integer] [ILAF integer] temperature-spec::= [FINAlt real] [FIRStt real] [TEMInc real] [TSTRuc real] [TWINDH real] [TWINDL real] [TBATh real] options-spec::= [IASOrs integer] [IASVel integer] [ICHEcw integer] [ISCAle integer] [ISCVel integer] [ISEEd integer] [SCALe real] [NDEGg integer] [RBUFfer real] [AVERage] [ECHEck real] [TOL real] cpt-spec::= See *note cpt:(pressure.doc) four-dim-spec::= [FIL4dimension] [SKBOnd] [SKANgle] [SKDIhedral] [SKVDerWaals] [SKELectrostatics] [K4DInitial real] [INC4Dforce integer] [DEC4Dforce integer] [MULTK4di real] [E4FILLcoordinates real] [FNLT4 real] [FSTT4 real] [TIN4 real] [IHT4 integer] [IEQ4 integer] [ICH4 integer] [TWH4 real] [TWL4 real]

Options common to minimization and dynamics The following table describes the keywords which apply to both minimization and dynamics. Keyword Default Purpose NSTEP 100 The number of steps to be taken. This is the number of dynamics steps which is also equal to the number of energy evaluations. INBFRQ 50 The frequency of regenerating the non-bonded list. The list is regenerated if the current step number modulo INBFRQ is zero and if INBFRQ is non-zero. Specifying zero prevents the non-bonded list from being regenerated at all. INBFRQ = -1 --> all lists are updated when necessary (heuristic test). IHBFRQ 50 The frequency of regenerating the hydrogen bond list. Analogous to INBFRQ ILBFRQ 50 The frequency for checking whether an atom is in the Langevin region, defined by RBUF, or not. IMGFrq 0 The frequency for the image update (only used if IMAGES or CRYSTAL is in use). The image update creates image atoms needed for the energy computation from the list of allowed symmetry transformations. Recommended value: 50, if a 2A buffer is used between CUTIM and CUTNB. IXTFrq 0 The frequency for the crystal update (only used of CRYSTAL is in use). The crystal update generates a new list of allowed symmetry transformations. This option is only required if the size or shape of the periodic box (i.e. CPT) can change during a simulation (or minimization). Recommended value: 1000 (if running CPT dynamics). non-bond- The specifications for generating the non-bonded list. -spec See *note Nbonds:(nbonds.doc). hbond- The specifications for generating the hydrogen bond list. -spec See *note Hbonds:(hbonds.doc). [ STRT ] STRT The dynamics is assumed to start from the input [ ] coordinates using an assignment of velocities given by [ ] IASVEL. No restart file is read. [ REST ] The dynamics is restarted by reading the restart file from unit IUNREA. TIMESTP 0.001 Time step for dynamics in picoseconds. The default value is 0.001 picoseconds. IUNREA -1 Fortran unit from which the dynamics restart file should be read. A value of -1 means don't read any file IUNWRI -1 Fortran unit on which the dynamics restart file for the present run is to be written. A value of -1 means don't read any file. Formatted output. IUNCRD -1 Fortran unit on which the coordinates of the dynamics run are to be saved. A value of -1 means no coordinates should be written. Unformatted output. IUNLDM -1 Fortran unit on which the biasing potentials, the histograms of the lambda variables of the dynamics run are to be saved. A value of -1 means no histograms should be written. Unformatted output (for details see node: output). IUNVEL -1 Fortran unit on which the velocities of the dynamics run are to be saved. -1 means don't write. Unformatted output. KUNIT -1 Fortran unit on which the total energy and some of its components along with the temperature during the run are written using formatted output. CRASHU -1 Fortran unit where a single DCL command file will be written. If the machine crashes before a restart file is written, this file won't be touched. If the crash occurs after a restart is written but before the run completes, this file will contain the line, "$ @CRASH". If the run completes, the file will contain the line, "$ @COMPLET". This allows for an automatic recovery system after crashes. NSAVC 10 The step frequency for writing coordinates. NSAVL 0 The step frequency for writing lambda histograms. NSAVV 10 The step frequency for writing velocities. NPRINT 10 The step frequency for storing on KUNIT as well as printing on unit 6, the energy data of the dynamics run. IPRFRQ 100 The step frequency for calculating averages and rms fluctuations of the major energy values. If this number is less than NTRFRQ and NTRFRQ is not equal to 0, square root of negative number errors will occur. ISVFRQ NSTEP The step frequency for writing a restart file. IHTFRQ 0 The step frequency for heating the molecule in increments of TEMINC degrees in the heating portion of a dynamics run. Zero means do no heating. IEQFRQ 0 The step frequency for assigning or scaling velocities to FINALT temperature during the equilibration stage of the dynamics run. NTRFRQ 0 The step frequency for stopping the rotation and translation of the molecule during dynamics. This operation is done automatically after any heating. FIRSTT 0.0 The initial temperature at which the velocities have to be assigned at to begin the dynamics run. Important only for the initial stage of a dynamics run. FINALT 298.0 The desired final (equilibrium) temperature for the system. Important for all stages except initiation. TEMINC 5.0 The temperature increment to be given to the system every IHTFRQ steps. Important in the heating stage. TSTRUC -999. The temperature at which the starting structure has been equilibrated. Used to assign velocities so that equal partition of energy will yield the correct equilibrated temperature. -999. is a default which causes the program to assign velocities at T=1.25*FIRSTT. TWINDH 10.0 The temperature deviation from FINALT to be allowed on the high temperature side.(+ve). i.e. high side of the temperature window. Useful during equilibration. TWINDL -10.0 The temperature deviation from FINALT to be allowed on the low temperature side.(-ve). i.e. low side of the temperature window. Useful during equilibration. TBATH FINALT The temperature of the heatbath in Langevin dynamics. When set to zero it allows one to do purely dissipative (quenched) dynamics. RBUF 0.0 Inner radius of the buffer, or Langevin, region sphere. All atoms with radial positions greater than RBUF angstroms are propagated by Langevin dynamics, if the dynamics keyword LANGevin has been specified. IASORS 0 The option for scaling or assigning of velocities during heating (every IHTFRQ steps) or equilibration (every IEQFRQ steps). This keyword does not control the initial assignment of velocities. .eq. 0 - scale velocities. (use ISCVEL option) .ne. 0 - assign velocities. (use IASVEL option) IASVEL 1 The option for the choice of method for the assignment of velocities during heating and equilibration when IASORS is nonzero. This option also controls the initial assignment of velocities (when not RESTart) regardless of the IASORS value. .eq. 0 - Use the comparison coordinate values in AKMA units (sorry) with the STRT option. If NTRFRQ is positive, then net trans/rot will be removed first. This option supresses other assignments of velocity. .gt. 0 - gaussian distribution of velocity. (+ve) .lt. 0 - uniform distribution of velocity. (-ve) kinetic energy of 3N velocity components are same. ISEED 314159 The seed for the random number generator used for assigning velocities. ISCVEL 0 The option for two ways of scaling velocities. .eq. 0 - single scale factor for all atoms .ne. 0 - a scale factor for each atom proportional to the kinetic energy average ratio between the system and along every degree of freedom for that atom. ICHECW 1 The option for checking to see if the average temperature of the system lies within the allotted temperature window (between FINALT+TWINDH and FINALT+TWINDL) every IEQFRQ steps. .eq. 0 - do not check i.e. assign or scale velocities. .ne. 0 - check window i.e. assign or scale velocities only if average temperature lies outside the window. ISCALE 0 This option is to allow the user to scale the velocities by a factor SCALE at the beginning of a restart run. This may be useful in changing the desired temperature. .eq. 0 no scaling done (usual input value) .ne. 0 scale velocities by SCALE. WARNING: Please use this option only when you are changing the temperature of the run. SCALE 1. Scale factor for the previous option. NDEGF computed Number of degrees of freedom to use in computing the temperature. If not specified on any call, the value is computed. This specification is not remembered between successive calls to dynamics. AVERAGE no When saving coordinates every NSAVC steps, this option will cause the average structure of the last NSAVC dynamics steps to be written instead of the final snapshot coordinate set. This option is primarily used for making smooth movies. ECHECK 20.0 The maximum amount the total energy may change on any step. TOL 1.0E-10 The shake tolerance (if SHAKE is in use). PCONst false Flag to indicate that constant pressure code will be used. PINTernal true Flag to indicate that the internal pressure will be coupled the reference pressure. PEXTernal false Flag to indicate that the external pressure will be coupled to the reference pressure. PCOUpling 0.0 The coupling decay time in picoseconds for the pressure. A good value for this is 5 ps. COMPress 0.0 The compressibility in atm**-1. A good value for proteins is 4.63e-5 PREFerence 1.0 The reference pressure in atmospheres. VOLUme computed The volume in Angstroms**3 to use for the pressure calculation denominator. This value is calculated if the CRYStal feature is use. TCONst false Flag to indicate that constant temperature code will be used. TCOUpling 0.0 The coupling decay time in picoseconds for the temperature. A good value for this is 5 ps. TREFerence FINALT The reference temperature for constant temperature simulations. MBOND Signifies that the dynamics run will be based on a multi-body simulation. If no bodies have been defined, this produces an error. Many of the standard dynamics options retain their meaning, in this case, but the following options are NOT supported: SHAKE, CONSTANT PRESSURE, NOSE, LEAPFROG, VER4, LANGEVIN. See *note Mbond:(mbond.doc)DynDesc for a description of the MBOND dynamics options. NLAT 0 The step frequency for writing instantaneous lambda temperature trajectories in lambda-dynamics. NBLCkfep 0 The step frequency for writing the energy decompostion trajectories in free energy calculation. ILAT -1 Fortran unit on which the histograms of the lambda temperature are to be saved. A value of -1 means no histograms should be written. IBLCkfep -1 Fortran unit on which the histograms of the energy decomposition are to be saved. A value of -1 means no histograms should be written. This file is used in post processing in TSM module. ILAP -1 Fortran unit on which the histograms of (Vi - Fi) are to be saved. A value of -1 means no histograms should be written. NSAVL is used for step frequency of printing (Vi - Fi) information. ILAF -1 Fortran unit on which the histograms of the restraining potential are to be saved. A value of -1 means no histograms should be written. NSAVL is used for step frequency of printing the restraining potential. TARGet false Flag to indicate that TMD code will be used. INRT 10 Number of step that one needs to get rid of artificial rotational motion in rigid TMD simulation. NTAR ntaom Number of atoms that are subject to TMD constraint. DINC 0.0 RMS increment in TMD simulation.

Recommended CHARMM input for dynamics. This section is intended only as a guide in setting up a dynamics simulation input file. Changes should be made as necessary according to personal tastes and project requirements. 1) For heating and early equilibration: DYNAMICS LEAP VERLET RESTART(*) NSTEP 20000 TIMESTEP 0.001(+) - IPRFRQ 1000 IHTFRQ 1000 IEQFRQ 5000 NTRFRQ 5000 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 100 NSAVC 100 NSAVV 0 INBFRQ 25 - hbond-spec nonbond-spec - FIRSTT 100.0 FINALT 300.0 TEMINC 100.0 - IASORS 1 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 10.0 TWINDL -10.0 (*) Except for first run, the use STRT in place of RESTART (+) If bonds to hydrogen atoms are SHAKEd 2) For late equilibration and analysis runs: DYNAMICS LEAP VERLET RESTART NSTEP 20000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 2000 IEQFRQ 5000(*) NTRFRQ 5000 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 - hbond-spec nonbond-spec - FIRSTT 100.0 FINALT 300.0 TEMINC 100.0 - IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 1 TWINDH 10.0 TWINDL -10.0 (*) Window checking should be disabled for the analysis run (i.e. IEQFRQ=0) if you want a real microcanonical ensemble. 3) For heating, equilibration and analysis runs using Langevin dynamics:(+) DYNA LEAP LANGEVIN STRT(*) NSTEP 20000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 - ILBFRQ 1000 RBUFFER 0.0 TBATH 300.0 - hbond-spec nonbond-spec - FIRSTT 300.0 FINALT 300.0 - IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0 (+) Note that the friction coefficients, in units of 1/ps, must first be initialized by filling the array FBETA with the SCALAR command SCALAR FBETA SET <real> <optional atom selection> (*) Except for first run, the use STRT in place of RESTART 4) For quenched molecular dynamics: For the first run (STRT), read velocities into the comparison coordinate set, or this should directly follow a former dynamics command. DYNA VERLET STRT(*) NSTEP 10000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 200 IEQFRQ 200 NTRFRQ 400 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 50 NSAVC 50 NSAVV 0 IHBFRQ 0 INBFRQ 25 - hbond-spec nonbond-spec - TSTRUC 300.0 FIRSTT 300.0 FINALT 0.0 TEMINC -30.0 - IASORS 0 IASVEL 0 ISCVEL 0 ICHECW 1 TWINDH 0.0 or equivalently with Langevin (dissipative) dynamics DYNA LANGEVIN STRT(*) NSTEP 10000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 4000 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 50 NSAVC 50 NSAVV 0 IHBFRQ 0 INBFRQ 25 - hbond-spec nonbond-spec - TSTRUC 300.0 FIRSTT 300.0 FINALT 300.0 - ILBFRQ 1000 RBUFFER 0.0 TBATH 0.0 - IASORS 1 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 (*) For first run, use RESTART otherwise The IASVEL 0 option causes the comparison coordinates to be used for the initial velocities (AKMA units). For subsequent runs the options IASORS 1 and IASVEL 1 may be used if random velocities are to be periodically assigned. 5) For constant temperature and/or pressure dynamics DYNA LEAP VERLET STRT(*) NSTEP 20000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 - PCONst PINTernal COMPress 4.63e-5 PCOUpling 5.0 PREFerence 1.0 - TCONst TCOUpling 5.0 TREFerence 300.0 - hbond-spec nonbond-spec - FIRSTT 300.0 FINALT 300.0 - IASORS 0 IASVEL 1 ISCVEL 0 ICHECW 0 TWINDH 0.0 TWINDL 0.0 6) For multi-body dynamics (assumes an atomistic equilibration has already been performed, substructures defined and modes generated): DYNA MBOND LOBA RESTART NSTEP 20000 TIMESTEP 0.001 - IPRFRQ 1000 IHTFRQ 0 IEQFRQ 0 NTRFRQ 0 - IUNREA 30 IUNWRI 31 IUNCRD 50 IUNVEL -1 KUNIT 70 - NPRINT 100 NSAVC 100 NSAVV 0 IHBFRQ 0 INBFRQ 25 - TCONst TCOUpling 5.0 TREFerence 300.0 - hbond-spec nonbond-spec - FIRSTT 300.0 FINALT 300.0 - MBPRlev -1

Running Molecular Dynamics The theoretical basis for dynamical simulations is elementary physics. The force on a particle is equal to the negative gradient of the potential energy of the particle. CHARMM can solve this equation numerically for all atoms in the molecule. A simple second order predictor two step method due to Verlet is used for integration. The choice of the integration step size is very important. One must weigh the increased accuracy of using a small step size against the longer real time that can be simulated with a given amount of execution time when a larger step size is used. The time step may be entered in picoseconds (using the TIMESTP keyword). CHARMM provides information on the accuracy of the numerical solution. Since the system has no external forces, the total energy should be conserved. Numerical errors will result in some fluctuations in the total energy so a good test is to compare the fluctuations in total energy to the fluctuations in kinetic energy as these fluctuations are proportional to the heat capacity of the system. See the next node for a description of dynamics output. Because the force constants for the bonds and bond angles are fairly large, it is reasonable under certain circumstances to constrain their values during dynamics. Such constraints are applicable if the harmonic motions are weakly coupled to other motions. The advantage of such constraints is that the step size of the numerical integration may be increased without sacrificing accuracy as these terms have the largest gradients in macromolecules simulated at physiological temperatures. We use the SHAKE algorithm for applying the constraints, see *note shake:(cons.doc)SHAKE. SHAKE can be applied to just the bonds involved with hydrogens, all bonds, all bonds and the angles involving hydrogens, or all bonds and angles. A dynamics run has basically four parts; initialization, heating, equilibration, and the simulation itself. Initialization means providing an initial position and velocity for all the atoms. Heating is the process of increasing the kinetic energy of the system up to a final temperature at which the simulation will be conducted. Equilibration is the process where the kinetic energy and the potential energy of the system evenly distribute themselves throughout the system. Only when the average temperature of the system stabilizes can one collect the trajectory information for analysis. The initial coordinates of a simulation are obtained after applying the minimization algorithm to a complete coordinate set. One cannot start with a system with a large potential energy as it will quickly heat up to unreasonable temperatures. For initializing the velocities, the user can specify velocities from the comparison coordinates (IASVEL 0), a uniform distribution of kinetic energy along each coordinate with random sign of the motion along each axis (IASVEL -1) or a Gaussian distribution of velocities (IASVEL 1 the default). The temperature at which velocities are assigned is determined by FIRSTT and TSTRUC by the algorithm: Tassign = 2*(FIRSTT-TSTRUC) + TSTRUC. For a harmonic system equilibrated to TSTRUC equal partition of the energy will result in an equilibrated temperature of roughly FIRSTT. If TSTRUC is not specified 1.25*FIRSTT will be used for assignment. Velocities may also be passed to dynamics in the comparison coordinate set (as opposed to assignment). This allows the user considerable flexibility in setting up the initial conditions. The heating of system is performed gently by increasing the kinetic energy by a small amount periodically. The number of integration steps between heating applications, the final temperature, and the kinetic energy increment are all user specified. In addition, there is a choice in the method of increasing the kinetic energy of the system. One may scale existing velocities or reassign them. The velocities can be scaled by either one scale factor calculated for the kinetic energy of the system averaged over many time steps or by scale factors established for each atom base ed on the ratio of its time averaged kinetic energy with that of the system. If reassignment is chosen, the velocities can have either a uniform or Gaussian distribution. To equilibrate the structure, one can specify a window around the final temperature where velocity adjustments will be made. The choice of velocity adjustments is the same as described above for heating. For the actual run, CHARMM will output the position and velocities of all atoms at intervals specified by the user. The temperature window can be set larger so that any gross conformational changes which result in a different potential energy will cause the temperature to be maintained. At any time energy is added to the system, the angular momentum of the system will be reduced to zero and translational motion will be stopped. One can also request that these operations be performed at any time during the dynamics run. The use of a restart file is essential for running dynamics. The restart file contains information about the most recent coordinate sets necessary for the VERLET algorithm. In addition the values of the energy accumulators are stored. All other information (such as SHAKE, fixed atoms, harmonic constraints, friction coefficients) has to be regenerated before invoking a dynamics restart. When the run is initiated, a restart file must be written using the IUNWRI keyword. As the dynamics routine complete NCYCLE, see *note Output::, steps of dynamics, the Fortran unit specified by IUNWRI will be rewound and a restart file will be written. In case of crashes, one has restart files corresponding to various points in the run. The CRASHU variable may prove valuable. Successive runs of CHARMM to continue the dynamics run must read the previous restart file using the IUNREA keyword and write it out for the next part of the run. Restarts may be done to reset various options, or to break up a long run into several shorter runs. Restart files will only run with the version of CHARMM they are created with. There are many numbers giving the frequency of actions to be taken during dynamics such as updating the non-bonded list, heating the molecule etc. Some of these numbers are adjusted along with the number of steps to run so that numbers all have a common divisor. At the present time, there are combinations which result in errors. At some point an attempt may be made to catalog all the actions, and check for erroneous processing. If one is interested in simulating the motion of part of the system with the rest of the system remaining fixed, it is possible to fix atoms in place, see *note fix:(cons.doc)fixed atom. If this is done, there are several effect on the dynamics. First, since the system is now anchored in space, the center of mass motion and total angular velocity is never stopped. Second, the number of degrees of freedom used for calculating the temperature is set to the number of free atoms times 3 minus 6. Third, the coordinate and velocity trajectory files will contain the position of the fixed atoms only once, and all other records will hold just the moving atoms. This saves a great deal of disk space. Trajectory files can be merged, broken in smaller pieces, and sampled at different intervals. Likewise, said operations can be performed on coordinate trajectories while rotating the coordinates to match a given coordinate set. When the DYNAmics command exits, the main coordinate set contains the final coordinate positions from the last energy evaluation and the comparison coordinates will contain the final velocities In AKMA units. Finally, a brief discussion of the Langevin dynamics algorithm is presented. The Langevin dynamics algorithm presently in CHARMM was intented to be used primarily with the "Stochastic Boundary Molecular Dynamics" method and consequently has been restricted to an algorithm which is valid only for the case FBETA*TIMESTEP<<1.0. That is for cases where relatively small friction coefficients are used. Typically values of FBETA*TIMESTEP up to about 0.3 still produce a stable dynamics which also satisfy the fluctuation-dissipation theorem. The algorithm itself reduces to the Verlet algorithm when FBETA is zero and consequently may be used to do regular dynamics, actually it is the same routine which does both dynamics. In using Langevin dynamics care must be taken to first initialize the array FBETA by using the scalar commands e.g., CHARMM >SCALAR FBETA SET <real> <atom selection> Failure to do this just means you are doing regular dynamics so no warning is given by CHARMM.

Contents of a dynamics output Note: This description of the output of a command is not normally going to be part of the documentation of commands. The dynamics output is sufficiently confusing that this description is necessary. The first part of CHARMM's output after a dynamics command lists all of the options that apply to that part of the run. Then, any information about velocity assignments (temperature changes) follows. Any time the velocities are changed in an anisotropic way, the motion of and about the center of mass will be stopped. This results in a printout both before and after this operation of the "DETAILS ABOUT CENTRE OF MASS". Its position and velocity are output followed by the components of the angular momentum. The last line gives the translation kinetic energy of the system, and thus one should expect a drop in the total energy and temperature of the system afterwards. Non-bonded interaction and hydrogen bond updates will appear intermittently and are cleared labeled. Every NPRINT steps, the total energy and various contributions will be printed. This output is preceded by a title which gives the correspondence of numbers to energy names. After IPRFRQ steps will appear the averages and RMS fluctuations. After the second such printout of averages and RMS fluctuations, the averages and RMS fluctuations for the run upto the last turning of the molecule will be given. This gives you longer range statistics. Such a calculation will not be done if IPRFRQ equals NTRFRQ. The ratio of total energy to kinetic energy fluctuations is an excellent measure of the accuracy of the run. After the averages are printed, a least squares fit of the total energy against the step number will be made to look for drift in the energy. Two such values are printed, one for the last IPRFRQ steps, and one to the previous turn. Next, the initial energy for the statistics, both short range and long, are printed. Finally, the correlation coefficient of the energy versus step is given for both ranges. A value close to zero indicates no systematic drift; a magnitude near 1 means you have a real problem with the dynamics. This process of printout continues until the end of the run is reached. Just before the last energy is printed will appear a message about the writing of coordinates and velocities to their respective files. Output of the lambda dynamics and post-processing (a) Output The output of the lambda dynamics, i.e. the histograms and the biasing potentials on the lambda variables, is written in a separate file from the coordinate file. Parallel to a regular coordinate file of the dynamic run, the lambda dynamic output file will automatically include a header and an integer array. They are used to provide the information on the values of NSTEP, NSAVL, NPRIV etc. To name a title for the output file (in complying with the CHARMM file requirement), the command LDTItle (similar to TITLE command) can be used. E.g. LDTItle * mte: Methanol to ethane * output for lambda dynamics * will write out a title before any other output data. The information on biasing potentials will also be written out. It takes a similar form as they were read in (see BLOCK.DOC), i.e. INTEGER : total No. of biasing potentials. I J CLASS REF CFORCE NPOWER : the format for individual one. INTEGER : the total no. of blocks. To specify the output fortran unit and the writing frequency, keywords IUNLDM and NSAVL are used. They are treated in the same fashion as IUNCRD and NSAVC. There is no separate restart file for the lambda dynamics. The information necessary for restarting a lambda dynamics is included at the end of a regular dynamic restart file. Thus, to restart the lambda dynamics is exactly same as restarting a regular dynamics run except you have to specify IUNLDM and NSAVL. E.g !input title for lambda i/o LDTITLE * This is a test * output for lambdas * open unit 11 writ form name output_file open unit 12 read form name input_file open unit 15 writ file name histogram dyna rest leap time 0.001 - nstep 10 nprint 1 iprfrq 10 - iunrea 12 iunwri 11 iuncrd -1 nsavc 1 IUNLdm 15 NSAVL 5 - first 300. - inbfrq 40 nbxmod 5 atom cdie shif vatom vdist vshif - cutnb 8. ctofnb 7.5 ctonnb 6.5 eps 1. e14fac 0.4 wmin 1.5 - cutim 8. imgfrq 40 Caution: the file is an unformatted output. However, the order of the output is very similar to a regular output file: (1) header, icntrl (automatically written) (2) title (3) total no. of biasing potentials (4) form of each biasing potential ( total = Nbias) (5) total no. of blocks (6) lambda**2 ( total = No. of blocks) (b) Post-processing WHAM ?. END

Reading and writing trajectory frames with direct commands This facility allows the creation or manipulation of trajectory files The main uses of this facility are; 1) creating artificial trajectory files from coordinate frames 2) reading an existing trajectory frame by frame for analysis that requires access to a variety of CHARMM commands 3) modifying an existing trajectory (copy with changes) such as minimizing each frame or other operations. [Syntax TRAJectory command] =================================================================== There are four commands that comprise this facility. 1) Initializing trajectory I/O TRAJectory {read-spec} {write-spec} read-spec:== [FIRST unit] [NUNITint] [SKIP int] [BEGIN INT] [STOP INT] write-spec:== [IWRIte unit] [NWRIte int] [NFILE int] [EXPAnd] [NOTHer int] [DELTa real] [SKIP int] FIRStunit (IREAd) - first unit to read from (default: do not read) NUNIts (NREAd) - number of units to read from (default:1) SKIP - skip value for both reading and writing (default:1) IWRIte - first unit to write to (default: do not write) NWRIte - number of units to write to (default:1) NFILe - number of frames on each output file (default: total) EXPAnd - flag to free fixed atoms in copying (only if reading) NOTHer - number of frames in previous files (if not reading) (d:0) DELTa - output delta value (if not reading) (default:0.001) 2) Reading a frame TRAJectory READ [COMP] The reading command does not have any specifiers other than whether the comparison or main coordinates will be used. 3) Writing a frame TRAJectory WRITe [COMP] The writing command does not have any specifiers other than whether the comparison or main coordinates will be used. 4) Query a trajectory file TRAJectory QUERy UNIT integer The query command rewinds an open trajectory file and then reads the header information from this trajectory file. It prints a summary and sets the following command line substitution parameters: 'NFILE' - Number of frames in the trajectory file 'START' - Step number for the first frame 'SKIP' - Frequency at which frames were saved (NSTEP=NFILE*SKIP when not using restart files) 'NSTEP' - Total number of steps from the simulation 'NDEGF' - Number of degrees of freedom from the simulation (Can be use to get the temperature with velocity files). 'DELTA' - The dynamics step length (in picoseconds). This command, again, rewinds the trajectory file upon completion. =================================================================== There are three modes of operation; 1) Create a new trajectory. The IWRIte and NFILe keywords must be used. The default values for the others are listed above. If several files will be made in different CHARMM runs that will be linked together later, the NOTHer keyword value should be increased by NFILe on each subsequent run. EXAMPLE: Create a "movie" trajectory that involves the rotation of a single sidechain (residue 21). COOR AXIS SELE ATOM * 21 CA END SELE ATOM * 21 CB OPEN WRITE UNIT 22 FILE NAME TYR21.ROT TRAJECTORY IWRITE 22 NWRITE 1 NFILE 360 SKIP 1 * trajectory showing the rotation of sidechain 21 * SET 1 1 LABEL LOOP COOR ROTATE AXIS PHI 1.0 SELE ATOM * 21 * .AND. .NOT. ( TYPE C - .OR. TYPE N .OR. TYPE H ) END TRAJ WRITE INCR 1 BY 1 IF 1 LT 360.5 GOTO LOOP STOP =================================================================== 2) Reading an existing trajectory The FIRSTU (or IREAD) keyword must be used. The default NFILe value is 1 and the remaining values if not specified will be read from the trajectory file. EXAMPLE: find the structure with the lowest energy and save it. OPEN READ UNIT 22 FILE NAME DYN1.TRJ OPEN READ UNIT 23 FILE NAME DYN2.TRJ TRAJECTORY FIRSTU 22 NUNIT 2 SKIP 100 SET 1 1 SET 9 9999.0 CALC NTOT = ?NFILE * 2 LABEL LOOP TRAJ READ UPDATE ... ! depending on how much your atoms move, ! you may leave this outside the loop GETE IF 9 LT ?ENER GOTO NEXT SET 8 @1 COOR COPY COMP SET 9 ?ENER LABEL NEXT INCR 1 BY 1 IF 1 LT @NTOT GOTO LOOP OPEN WRITE CARD UNIT 12 NAME LOWE.CRD WRITE COOR COMP CARD UNIT 12 * structure with the lowest energy * frame number @8 with energy @9 * STOP =================================================================== 3) Copying from one trajectory to another. The operation of this command works in the same mode as the MERGE command, except a variety of CHARMM commands can be executed between reading and writing of frames. EXAMPLES: Create a new trajectory where every frame is minimized for 200 steps. OPEN READ UNIT 22 FILE NAME DYN.TRJ OPEN WRITE UNIT 32 FILE NAME DYN.MIN TRAJECTORY NUNIT 22 SKIP 100 IWRITE 32 * minimized trajectory * SET 1 1 LABEL LOOP TRAJ READ MINI ABNR NSTEP 200 TRAJ WRITE INCR 1 BY 1 IF 1 LT ?NFILE GOTO LOOP STOP

Merges or breaks up a trajectory into different numbers of files Frequently, one generates a trajectory into small files to minimize the CPU time of one job. However, so many files are usually hard to manage so it is desirable to merge said files into larger units. This command provides that capacity. In addition, it is possible to break up the trajectory into smaller pieces and to sample the trajectory less frequently than originally generated. Another option is to optionally rotate the structure at each frame to least squares fix a reference structure. [Syntax MERGE dynamics trajectories] MERGE [ COOR ] [FIRSTU unit-number] [NUNIT integer] [SKIP integer] [ VEL ] [OUTPutu unit-number] [NFILE integer] [ DRAW ] [BEGIN integer] [STOP integer] [first-atom-selection] [ XFLUct ] [ UNFOld ] [ ORIEnt [MASS] [WEIGht] [NOROt] [PRINT] second-atom-selection ] [ RECEnter second-atom-selection] Keyword table Option Default Purpose [COOR] COOR Specification of the type of trajectory file. COOR is [VEL ] coordinates; VEL is velocities. [DRAW] Make a CHARMM movie (do not write any files, just display) FIRSTU 51 The first unit of the trajectory to be read. NUNIT 1 The number of units to be read starting with FIRSTU SKIP 1 Only those coordinate whose dynamics step number modulo SKIP will be reoriented and written out. OUTPUTU 61 The first unit number of the output trajectory NFILE The number of coordinate sets written to each output merged file. If left out, this will be set to the number of coordinates in the first input file times the number of input files. WARNING: This default will generate a bad trajectory file if SKIP is not set to the interval actually present in the trajectories. Further, if you set its value to be larger than the number of coordinates that are actually written in any output file, you will have problems. The error that is generated results from the control array in the beginning specifying that there are more coordinates than actually exist in the file. EOF errors will result when the trajectory is read. BEGIN First step number to start reading from STOP Last step number to read first-atom-sel Selection of atoms to include in the output file. RECEnter Will re-center atoms based on the existing IMAGE transformations ("coor ... rece ..") thus HAS to be preceded by a normal image setup (read image, image byresidue ..) for the atoms (usually solvent) that are to be transformed as if the center of the primary box coincided with the center of geometry of the atoms in the second selection. In short: The second selection defines the origin of your lattice and the solvent molecules are put as close as possible to the solute, even if things drifted slightly out of the box during the simulation. Useful for calculation of solvation properties. Does not work with XFLUct or UNFOld. The possibly large amounts of output reporting on all transformation operations being performed may be suppressed by setting PRNLev 4, or PRNLev 1 to get rid of the SELECTE IMAGES BEING CENTERED messages as well. NOTE: Uses SAME second-selection as ORIENt, so if both RECEnter and ORIEnt are specified there should only be one second selection, which will first be used to define the recentering, then for the orienting. UNFOld removes the effects of image centering (not the same as RECEnter), ie will let a particle continue out to the right if that is what it was doing when PBC moved it back into the primary box during the simulation. XFLUct removes the effects of the box size/shape changes from constant pressure simulations. This allows an accurate calculation of transport properties (diffusion constants,..) from CPT trajectories. ORIEnt Flag to specify best fit rotation and translations. MASS Use mass weighting in best fit. WEIGht Use weighting array for best fit weights. NOROt Only translate in the best fit. PRINT Print the details of best fit second-atom-sel Selection of atoms to use in the best fit. NOTE: Uses SAME second-selection as RECEnter (see above) The title of the output trajectory will be copied from the input trajectory. NOTE: If the input trajectory is from a CRYStal simulation, the CRYStal setup has to be invoked also before the MERGe operation if CRYStal data is to be written out to the resulting trajectory.

Reorienting a coordinate trajectory If one is interested in reorienting every set of coordinates found in a dynamics trajectory with respect to some reference structure, one can use the ORIEnt option in conjunction with the MERGe command. The process of reorienting a coordinate trajectory works as follows: A series of files containing the trajectory are assigned to successive units prior to a CHARMM run. The coordinates stored therein are presumed to have been written every NSAVC steps. CHARMM will read each coordinate, select some periodically, reorient them, and write them to successive units where each output file will have a user specified number of coordinates. The following table lists the options involved: Option Default Purpose ORIE .false. Specify that a least squares RMS fit will be done. MASS .false. Use a mass weighting in the fit WEIGH .false. Use the weighting array (wmain) in the fit NOROt .false. Just shift the centers to best fit. PRINt .false. Print what happened to each coordinate set. atom-selection all Select which atom to use in the fit. If atoms were fixed during the dynamics, the new trajectory produced will not have fixed atoms because the rotations applied to each coordinate set will be different thereby yielding different coordinates for the fixed atoms. Fixing the coordinates leads to a large space reductions, so the reorientation process will therefore result in potentially much larger trajectory files. See *note fix: (cons.doc)Fixed Atom.

Computes the RMS difference between two trajectory files and make a matrix of results. Large files should be reduced with the MERGe command before processing this command. RMSDynmics ORIEnt [MASS] [WEIGht] [NOROt] [RMS] [atom-selection] FIRSTU unit-number [NUNIT integer] BEGIN integer [SKIP integer] STOP integer [IWRIte unit-number] [SECU unit-number] [BEG2 integer] [SKP2 integer] [STP2 integer] ( [IREAd unit-number] [JREAd unit-number]) ([IMAGes]) [MATRix] [PQUNit unit-number [PQSEed integer] [IOPT int] [MAXFn int] [NSIG int] ] IWRIte int - Unit for the output matrix. FIRSTu int - Unit number for first file containing trajectory 1 NUNIt int - Number of files for trajectory 1 BEGIn int - Starting step number STOP int - Ending step number SKIP int - Number of steps to skip (default 1) NB! BEGIN, SKIP and STOP have to be specified to allow proper memory allocation! The TRAJ QUERY command can retrieve these values The trajectory(ies) are read into memory before calculations begin; Memory usage may be reduced by decreasing the number of selected atoms, and by reducing the number coordinate sets (frames) used - increase SKIP SECU int NUN2 int BEG2 int Specifications for trajectory 2. Defaults to same values SKP2 int as used for trajectory 1. If SECU is not given, or same as STP2 int FIRStu only one trajectory will be analyzed. (IREAd int - unit number of the first trajectory file.) If only IREAd is specified, or if IREAd and JREAd are the same the RMSDs will be between frames in same file (JREAd int - unit number of the second trajectory file.) NB! IREAD and JREAD are obsolete as of c30 (but still supported) (IMAGes - Use image atoms for the analysis (not implemented)) ORIEnt - Do best fit of structures MASS - Use a mass weighting in best fit. WEIGht - Use the weighting array in best fit. NOROt - Best fit without letting the structures rotate. RMS - Do RMS fit between structures, otherwise, align structures with the axis. MATRix - output just the RMSDvalues in matrix format NOTE: You have to specify correct BEGIN and STOP values so that the correct amount of memory can be allocated PQUNit int - unit number for output of (P,Q)-values in a 2D-projection of the RSMD-map according to Levitt, M. J.Mol.Biol. (1983) 168, 621-657. CHARMM variable PQRES is set to the final value of the target fcn. This can be slow to converge. PRNL 6 (or greater) outputs the value of the target fcn for each iteration, allowing judgement of the convergence. If iterations stop due to maximum number of function evaluations reached, you can increase this with the MAXFn <integer> keyword, and see if it results in a qualitative change of the p{p,q}-pattern. If not all is probably OK. The RMSD-values can be printed also when this option is on. This turns on the MATRix flag so BEGIN and STOP have to be correct. Uses IMSL-routine ZXMIN; IOPT,MAXFN and NSIG can be specified to tune the behavior of ZXMIN. PQSEed int - seed for random generation of initial guesses for the (P,Q). prnl 6 gives a little more information about initial guesses etc. atom-selection - Atoms to use in the fitting procedure.

Format or unformat a dynamics trajectory DYNAmics FORMat FIRStunit <unit> NUNIt <int> BEGIn <int> SKIP <int> STOP <int> OUTPut <unit> OFFSet <int> SCALe <int> MODE <FORTRAN-FORMAT> DYNAmics UNFOrmat INPUt <unit> OUTPut <unit> These commands allow to convert binary trajectory files into a machine independent yet compact format and to convert them back into binary files. The defaults for OFFSet, SCALe and MODE are: OFFSet=600, SCALE=10000, and MODE=12Z6. The trajectory is converted into positive integers according to the formula <integer>=INT(<real>+OFFSET)*SCALE). The user has to make sure that all coordinates of the trajectory are within OFFSET angstroms. The precision may be increased by choosing a larger SCALE and FORTRAN-format, e.g. MODE=11Z7 OFFSET=100000. ("Z" is the hexadecimal format and is available on most machines.)

CONSTANT VELOCITY A constant velocity method has been developed for use with DYNA (right now, it only works with LEAP [in charmm] and LOBATTO [in MBO(N)D] integrators). The main purpose of this facility is to run simulations similar to atomic force microscopy. The constant velocity method, therefore, is used in conjunction with the NOE facility used to apply a 'spring' between two atoms. A constant velocity for an atom is entered via CVEL in CHARMM syntax: CVELocity <real> <sele first atom> <sele second atom> where <real> = constant velocity in Angst/ps; the constant velocity vector and direction is defined from <sele first atom> to <sele second atom>. the position of the <sele second atom>, typically a dummy atom, is moved to the position of <sele first atom> + 0.0001 Angstr. along the vector (because charmm does not like duplicate coordinates); <sele second atom> then traverses along the vector at the constant velocity rate. The second atom is not really needed, but it is helpful in analyzing the vector visually before running dynamics. Note: If you want to apply a spring between the constant velocity atom and the first atom in the vector, you must use (currently) the NOE facility in charmm. --- Here are the relavent syntax from a sample input file (typical usage). --- *afm.inp *Simulated Atomic Force Microscopy *Continually loops over 10ps segments of dynamics (NVT'ish) * ...lots of typical charmm stufff... !--------DEFINITIONS !Two atoms, one is the near the end of myosin, the other is a dummy atom ! to be cvel'ed define tip SELE atom dumm 1 dumm END define pp SELE atom hc 835 ca END !Actin binding region define actb sele segid hc .and. (resi 405:415 .or. resi 529:550 .or. - resi 626:647) end !----FORCES set f 4. !spring constant; See Grubmueller Science 1996, 271, 997 set com 100 !force used to pin actin binding site set max 80 !tot number of dyn runs--arbitrary right now !##CVEL <Angst/ps> <first_single_atom_selec> <second_single_atom_selec> !## These two atoms define the pulling vector; the first selection !## is the pull point, and the second selection is the atom that moves at !## constant velocity along the pull point. Currently, the 'spring' !## between these two atoms is defined using the NOE facility below. cvel @{cv} SELE pp END SELE tip END !set up spring between atoms in cvel noe reset assign SELE pp END SELE tip END - kmin 0.0 rmin 0.0 kmax @f rmax .00001 fmax 1000 PRINT ANAL end label skip !----Pin protein cons harm sele actb end force @{com} DYNAMICS MBOND (or LEAP) (re)START - dynamics equilibration or constant temperature method. !lots of loops over the above stop References: 1. Grubmueller Science 1996, 271, 997. 2. "The Evaluation Of Multi-Body Dynamics For Studying Ligand-Protein Interactions. Using MBO(N)D To Probe The Unbinding Pathways Of Cbz-Val-Phe-Phe-Val-Cbz From The Active Site Of Hiv-1 Protease" Chin, D. N.; Haney, D. N.; Delak, K.; Chun, H. M.; Padilla, C, In Rational Drug Design; Parrill, A., Reddy, R. Eds.; ACS Washington, 1998, in press. MODIFIED CODE IN CHARMM 26?? ------------------ a build/sgi/newmk/charmm.mk 9 blocks a build/sgi/newmk/dynamc.mk 25 blocks a build/sgi/newmk/mbond.mk 17 blocks a source/fcm/newfcm/cveloci.fcm 2 blocks a source/dynamc/newsrc/cveloci.src 9 blocks a source/dynamc/newsrc/dcntrl.src 110 blocks a source/dynamc/newsrc/dynamc.src 104 blocks a source/charmm/newsrc/charmm_main.src 47 blocks a source/charmm/newsrc/iniall.src 57 blocks a source/moldyn/newsrc/compin.f 24 blocks a source/moldyn/newsrc/delta_v.f 19 blocks a source/moldyn/newsrc/engmom.f 21 blocks a source/moldyn/newsrc/engmom_ke.f 9 blocks a source/moldyn/newsrc/mbdyna.f 58 blocks a source/moldyn/newsrc/ydot.f 74 blocks a source/moldyn/newsrc/CHARMM.INC a source/mbond/newsrc/mbback.src 52 blocks a source/mbond/newsrc/mbdyn.src 40 blocks

Implementation of Targeted Molecular Dynamics (TMD) This is only for TMD using holonomic RMS constraint. For doing Targeted Molecular Dynamics (TMD), one needs to define a moving coordinate and a target coordinate. You slowly pull the moving structure towards the target structure by gradually decreasing the RMS distance between two. The 'pulling' speed is defined by user. Several new keywords for charmm dynamic command (DYNAmic) have been added. To invoke TMD, very similar to Berendsen's constant temperature method, one needs to start with the keyword TARGet, then followed by a few others defined detailed parameters such as the number of atom that subject to TMD constraint (NTAR), number of dynamic steps for stopping artificial rotation (INRT), and RMS increment (DINC).

Molecular Dynamics in the Tsallis (Generalized) Ensemble Molecular dynamics that yield averages for a Tsallis (generalized) ensemble rather than a canonical one can now be performed. At present, this method is implemented only for the leapfrog Verlet integrator (dynamc.src); the TSALLIS keyword must be included in pref.dat when compiling. The method is invoked by adding to the DYNAmics command the keywords: TSALlis QTSAllis real EMIN real where QTSAllis is the Tsallis q and EMIN is the estimated minimum energy of the system. The default value of QTSAllis is 1, in which case the method reduces to standard (canonical) dynamics. Values of q larger than 1 effectively correspond to a smoothed potential in which the positions of the extrema are preserved. Estimates for EMIN should err lower than any possible energy of the system encountered during the simulation. It is important to note that the scale factor for the forces involves a temperature. The temperature employed in the Tsallis transformation corresponds to the one used to assign the velocities during heating and equilibration (TSTRUC initially and then FIRSTT + int*TEMINC). Thus, it is important to set FIRSTT, FINALT and TEMINC appropriately even if one is running Langevin dynamics at temperature TBATh (i.e., for equilibrium dynamics, set FIRSTT = FINALT = TBATh). Reference: Andricioaei, I. and Straub, J. E. (1997) On Monte Carlo and molecular dynamics methods inspired by Tsallis statistics: Methodology, optimization, and application to atomic clusters. J. Chem. Phys. 107, 9117-9124.

Description of the CENT Command [ CENT NCRES int ] Description: The reCENTering command allows to recenter the system at the geometric center of the first NCRES residues in the psf file. This keyword is useful when modelling a protein/water system using the periodic boundary conditions to prevent the protein from driffting outside of the primary unit cell. It can be replaced by the IMAGe keyword when the solute is a small organic molecule. Syntaxis: The Keyword CENT, which is specified in the DYNAmics command line, turns on the recentering option for the system at the start of a dynamics calculation (dcntrl.src) and at each update of the nonbonded list (heuristic.src) NCRES int - The first N (int) residues in the psf file, based on which the system will be centered.

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology

FDA/CBER/OVRR Biophysics Laboratory

Modified, updated and generalized by C.L. Brooks, III

The Scripps Research Institute