Molecular Dynamics

Description of Molecular Dynamics
Molecular dynamics is one of CHARMM's primary features as it enables users to simulate molecular movement with Newton's second law of motion. This section of the tutorial will explain how to create a basic CHARMM input script for a molecular dynamic simulation.

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.