MD

Description of Molecular Dynamics
Molecular dynamics simulates the motion of a system by numerically integrating Newton's second law of motion.The objective of a MD simulation is usually to determine the time correlation between events or to get a statistically-valid collection of structures (called an ensemble) that satisfies ergodicity. Currently, MD simulations tend to run from several nanoseconds to around a microsecond (although millisecond scale MD simulations are possible on advanced hardware). This is often too short to simulate biological processes such as the folding of complex polypeptides. Biophysicists often therefore look to use the results of MD simulations to make statistical arguments about the behavior of a structure instead of directly observing long scale behavior. This tutorial will focus on the practical issues of setting up a molecular dynamics simulation in CHARMM. We will begin with a discussion of the prerequisites for undertaking a successful simulation.

Prerequisites
In order to begin an MD simulation, you should have a structure that has been reasonably well minimized under the exact same conditions as you are planning to use for the dynamics. This means that the non-bonded set-up (including Ewald) and periodic boundary conditions should be exactly the same. However, in many cases it is undesirable to minimize the structure too much, as it may deform in an undesirable manner. Please see the Minimization page for more details.

If there are problems with the way that the model was built, they will likely manifest themselves in the dynamics run. In many cases, the potential energy of the system will explode to unrealistic levels. In other cases, the structure will twist into physically impossible conformations. You should keep an eye on the energy and position of your system during dynamics to make sure that this does not happen. Visualizing the MD trajectory is highly recommended as a check for possible problems. Also, pay careful attention to the warnings CHARMM prints out and resist the temptation to disable ECHECK (or set it to an unreasonably high value). Exceeding the energy change toleration (ECHECK) and deviations in SHAKE can indicate model building errors.

It is usually necessary to heat and equilibrate a structure at the desired temperature after minimization. This topic is discussed further below. For finicky structures, it might be necessary to begin the equilibration at a shorter time step (e.g. 0.5 fs instead of 1 fs) to keep within the energy change tolerance.

Choosing an ensemble
As mentioned above, one of the purposes of running molecular dynamics is to sample a collection of structures.

A note about the DYNAmics command
Molecular dynamics simulations in CHARMM are run via the DYNAmics command. This command takes a lot of options, and for that reason the first six characters of each subcommand name are significant (as opposed to the first four elsewhere in CHARMM). All of the key words listed on this page of the tutorial are options to the dynamics command. Some of the more important options that will be used in almost every MD run are:


 * NSTEP : specifies the number of steps to be run
 * TIMESTep : The time step in picoseconds, 0.001 is 1 fs.
 * NPRINT : specifies the frequency at which the energy should be printed out, e.g. NPRINT 100 prints the energy every one hunder steps.
 * NSAVC : frequency to write the coordinates to a trajectory file
 * IUNCRD : Unit number of the coordinate trajectory file (must by opened before DYNAmics are invoked)

A complete dynamics set-up might look like:

! open the trajectory file open unit 50 write unform name trajectory.trj dyna start nstep 1000000 timestep 0.001 nprint 1000 nsavc 100 iuncrd 50

This will run 1 ns of molecular dynamics, saving a coordinate trajectory frame every 50 steps (0.05 picoseconds).It is also possible to save the energy, temperature, and velocities to a trajectory, but this is less used in practice. The only option that you have not seen here is START, which tells CHARMM to start a new run (as opposed to continuing an old run). We'll talk about restarting molecular dynamics simulations below.

DYNAmics, like other CHARMM commands, will use the previously established non-bond configuration unless these are over-ridden in the DYNAmics command itself. The recommended practice is to set up your non-bond options in advanced so as not to clutter up the DYNA command (it has enough options already!).

Basic MD Input Script
Preparation After the topology file, parameter file, and coordinates have been read (see the Basic Input Script section of the tutorial if you are unfamiliar with reading in structures), the molecular dynamics (MD) simulation part of the input script can be started. This tutorial, however, will use a PDB that has already been minimized, solvated, and then minimized again. This thus assumes that the prior parts of the tutorial have been done. The files needed are available at this tutorials web site: http://www.charmmtutorial.org.

The system is a solvated protein where the boundaries of the system are determined using the crystal command. This command defines the repeating unit cell of the solvated protein, which is typically used to represent the protein in a liquid solution. Instead of setting up all the crystal and images information again, the solvated and minimized structure can be read in and the crystal can be easily rebuilt and setup with the transformation file being read:

Before the actual molecular dynamics input command line can be written, SHAKe is used to restrict hydrogen bond lengths (as in minimization).

Now the dynamics command will be written. Since MD simulations can take from a few minutes to a few months, depending on the system size, CHARMM offers a feature to help make longer runs more practical through restart files. Restart files created by CHARMM store the coordinates and velocities of the current systems atoms and allow the user to restart a run or make a longer run into several shorter runs. To restart a previous simulation, two restart files should be opened, one to be read and one to be written to. If the user is starting a simulation from scratch, then only the written file is needed so for this tutorial, therefore the OPEN UNIT 51 READ restart file line will be commented out for the first run. This can be done with a standard open read/write command that should resemble the following:

For MD runs, CHARMM can create a trajectory file that contains the coordinates of the system and therefore we have to tell CHARMM to open the trajectory then write to it.

CHARMM knows it will be performing a MD run with the "DYNA" command. After the "DYNA" command, "LEAP" or "CPT" should be specified. LEAP alone results in a simulation that is constant in number of particles, volume, and energy (NVE). LEAP should be used if the protein has not been solvated or if periodic boundary conditions are not to be used. If the protein has been solvated and periodic boundary conditions are being employed then "CPT" (constant pressure and temperature) should be used to simulate experimental conditions. Throughout the remainder of this tutorial, it will be assumed that the protein is solvated and PBC are applied. After "CPT", CHARMM should know if this is a new simulation or a continuation from a previous run. If the initial coordinates are based on a previous run, then the "RESTart" command is used. If this is a new run then STRT should be used. The "IUNRea" and"IUNWri" commands followed by a number define the read and write restart file, respectively (the commands are sent the unit number to the respective file). Similarly, "IUNCrd" defines the trajectory file. If the user wants to write out the total energy, its components, and temperature then they can use the "KUNIt" command followed by the unit number. This tutorial will not take advantage of the kunit feature so a value of -1 will be assigned to it.

Following the "STRT" command will be the "NSTEp" command, which is the number of dynamics steps that CHARMM will perform. For this tutorial, NSTEp will be set to 1000. For "TIMEstep", CHARMMs unit of time measurement is picoseconds (ps). In this tutorial a TIMEstep of 0.001 picoseconds (one femtosecond) will be requested using the command "TIMEstep 0.001". The number of steps times the timestep results in the total time of the simulation. This simulation will represent 1 ps (0.001 * 1000). Just in case the MD run goes out of control, CHARMM can terminate the simulation based on the "ECHEck" command value. ECHEck is the amount of total energy that can change between steps. The number should be in the hundreds for small systems and in the thousands for larger systems. If the simulation is to go haywire then the energy change will probably be much greater, but the maximum value ECHEck can be is 9999. This tutorial will set ECHEck to 500. Since "CPT" was used, the command "PCONs" follows and indicates to CHARMM that the user wants to run constant pressure dynamics. Another indicator, "PINTernal", informs CHARMM that internal pressure will join together with the reference pressure. To set reference pressure, use the keyword "PREF". The units are in standard atmospheres so in this tutorial the value of 1.0 will be used. The command "PMASs 20" sets the mass of the pressure piston to 2000 and "PGAMma 0.0" is the Langevin piston collision frequency. As PMASs and TMASs are dependent on system size the following algorithm is helpful for determining optimal values:

So far, the input script should resemble the following:

This tutorial will apply Nosé-Hoover MD and more information regarding this method can be found in the CHARMM documentation (http://www.charmm.org/documentation/c34b1/nose.html). To tell CHARMM the users intent, use the word "HOOVer" followed by the reference temperature (in kelvin) and the mass of the thermal piston (kcal * mol-1 * ps 2). Reference temperature can be set by "REFT" and the mass of the thermal piston can be set with "TMASs".

CHARMM needs to know the frequency to regenerate the non-bonded list, which is given by the INBFrq parameter. -1 is used then the list will only be updated when deemed necessary.

To inform CHARMM that it atom based van der Waal cutoffs will be used the "ATOM" and "VATOm" keywords are added. Additionally, the non-bond values have to be set: cutnb, ctonb and ctonnb.

Now the dynamics command should resemble the following:

Next the Ewald parameters need to be set. Additionally, for analysis purposes, "IPRFrq" is the step frequency for when CHARMM should calculate the averages and RMS fluctuations of major energy values. This should be an integer. This tutorial will set IPRFrq to 5000.

IHTFrq controls how often the temperature should be increased and this tutorial will set it as 0 since the user does not necessarily have to have a steady increase in temperature. If the user were to start a simulation from crystal coordinates, then the user will need to use this command to gradually increase the temperature to avoid simulation artifacts. Likewise, heating may be used to bring a low energy structure, such as one resulting from heavy minimization, up to a desired level.

While IHTFrq is active during the heating state, IEQFrq controls the step frequency of increasing or decreasing the velocities of the atoms, which then create a change in heat during the equilibrium stage (which will be explained later in this tutorial). This portion of the tutorial will set "IEQFrq" to zero.

NTRFrq is a CHARMM implemented check that makes sure no outside unnatural forces are acting upon the molecule causing it to rotate and translate. The number after the "NTRFrq" command is the step frequency in which CHARMM stops all rotation and translation of the molecule.

The input files DYNA command should now look like the following:

For diagnostic purposes, it is a good idea to print out the energy data to the output file every so often, which is controlled with NPRInt. 1000 is a general value for nprint, but since a shorter simulation is being requested NPRInt will be set to 100. Saving the coordinates frequently is another good idea and can be done with the command "NSAVc". Like NPRInt, the value of NSAVc will be set to 100, however, for longer simulations coordinate saving frequency should be increased. If the user wants to save the atoms velocities, which can be very helpful in analysis but can take massive amounts of space (depending upon the system size), the user can use the command "NSAVv" followed by the step frequency for saving velocities. The "IHBFrq" controls the frequency at which to regenerate the hydrogen bond list and will be set to 0 for this tutorial because modern force fields do not contain explicit hydrogen bond potentials. The last two lines should be the heating lines.

To specify the temperatures, CHARMM has created several commands depending on the stage (the stages are: initialization, heating, and equilibrium). The initial temperature can be specified with the "FIRStt" command. "FINAlt" is the desired equilibrium temperature for the system. The temperature increments can be specified with "TEMInc" although the value will not matter for the heating script since "IHTFrq" is set to 0. To tell CHARMM the temperature at which the structure was equilibrated, "TSTRuct" is used. TBATh is the temperature of the heat bath for Langevin dynamics. This tutorial will use values that run the simulation at approximately room temperature.

To initialize the velocities by assignment, "IASOrs" should be set to 1. IASOrs should be set to 0 if you are continuing this simulation from a restart file since reassigning the velocities is not desired. IASVel assigns velocities but only when these velocities are reinitialized. ISCVel controls the scaling of velocities, which is not used in this tutorial, so it should be set to 0. If the user is doing an NVE calculation rather than CPT, temperature will vary and therefore there should be a check to make sure the temperature is staying within the window. If the user is performing a CPT (Constant Pressure and Temperature) calculation, use the following command to ignore the temperature window check:

because the Hoover thermostat controls it. If an NVE calculation is being performed, set ICHEcw to 1 and TWINdh represents the highest temperature allowed and TWINdl represents the lowest temperature allowed.

Your finished DYNA command should resemble the following:

With the rest of the input file writing out the finished coordinates.