Molecular Dynamics

From CHARMM tutorial
Jump to navigationJump to search

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:

 * Molecular Dynamics
 *
 bomlev -2
 
 ! Read in Topology and  Parameter files
OPEN UNIT 1 CARD READ NAME top_all27_prot_na.rtf
READ RTF CARD UNIT 1
CLOSE UNIT 1

OPEN UNIT 1 CARD READ NAME par_all27_prot_na.prm
READ PARA CARD UNIT 1
CLOSE UNIT 1

! Read protein structure from PSF file
OPEN READ FORMATTED UNIT 27 NAME 1yjp-full-solv-min.psf
READ PSF CARD UNIT 27
CLOSE UNIT 27

OPEN READ UNIT 10 CARD NAME 1yjp-full-solv-min.crd
READ COOR CARD UNIT 10
CLOSE UNIT 10


! @XDIM IS THE SIZE OF THE BOX FROM THE SOLVATION SCRIPT
CRYStal DEFIne CUBIc @XDIM @XDIM @XDIM 90.0 90.0 90.0
CRYStal BUILD NOPER 0

OPEN UNIT 1 READ CARD NAME 1yjp_cube.xtl
CRYStal READ CARD UNIT 1
CLOSE UNIT 1

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:

!RESTART FILE THAT WILL BE WRITTEN

OPEN UNIT 41 WRITE CARD NAME 1yjp-md.res

!RESTART FILE THAT WILL BE READ

OPEN UNIT 51 READ CARD NAME dyn.rea

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.

!TRAJECTORY FILE NAME
OPEN UNIT 31 WRITE FILE NAME 1yjp-md.dcd

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:

! CALC PMASS = TOTAL MOLECULAR MASS / 50! CALC TMASS = PMASS * 10
SCALAR MASS STAT CALC PMASS = INT ( ?STOT 0.0 ) CALC TMASS = @PMASS * 10

So far, the input script should resemble the following:

DYNA CPT STRT -
  TIMEstep 0.001 -     !timestep in picoseconds
  NSTEp 1000 -         !number of steps and energy evaluations
  ECHEck 500. -        !max variation of energy step-to-step
  PCONs -              !use constant pressure
  PINT -               !use internal virial for calc pressure
  PREF 1.0 -           !reference pressure in atm
  PMASs 20 -           !piston mass
  PGAMma 0.0 -         !Langevin piston collision frequency

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* ps2). 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:

DYNA CPT STRT -
  TIMEstep 0.001 -     !timestep in picoseconds
  NSTEp 1000 -         !number of steps and energy evaluations
  ECHEck 500. -        !max variation of energy step-to-step
  PCONs -              !use constant pressure
  PINT -               !use internal virial for calc pressure
  PREF 1.0 -           !reference pressure in atm
  PMASs 20 -           !piston mass
  PGAMma 0.0 -         !Langevin piston collision frequency
  HOOVer -             !Hoover temperature control
  REFT 310.15 -        !Hoover reference temperature
  TMASs 190 -          !mass of thermal piston
  INBFrq -1 -          !lists updated when necessary (heuristic test)
  ATOM -               !electrostatics calculation (default) is done on atom-atom basis
  VATOm -              !needed for nonbond list generation
  CUTNb 16.0 -         !distance cutoff for neighbor list
  CTOFnb 12.0 -        !max distance for pair to consider in energy
  CDIE -               !use constant dielectric/r
  EPS 1. -             !potential for constant ditlectric
  CTONnb 8. -          !Distance at which smoothing function reduces pair's contribution
  VSWItch -            !Use switching function for van der Waals to avoid cutoff discontinuity
  WMIN 1.0 -           !warning cutoff for minimum atom-atom distance
  CUTIm 16.0 -         !cutoff for including image atom
  IMGFrq -1 -          !frequency to check for atoms in Langevin region

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:

DYNA CPT STRT -
  TIMEstep 0.001 -     !timestep in picoseconds
  NSTEp 1000 -         !number of steps and energy evaluations
  ECHEck 500. -        !max variation of energy step-to-step
  PCONs -              !use constant pressure
  PINT -               !use internal virial for calc pressure
  PREF 1.0 -           !reference pressure in atm
  PMASs 20 -           !piston mass
  PGAMma 0.0 -         !Langevin piston collision frequency
  HOOVer reft -        !Hoover temperature control
  REFT 310.15 -        !Hoover reference temperature
  TMASs 190 -          !mass of thermal piston
  INBFrq -1 -          !lists updated when necessary (heuristic test)
  ATOM -               !electrostatics calculation (default) is done on atom-atom basis
  VATOm -              !needed for nonbond list generation
  CUTNb 16.0 -         !distance cutoff for neighbor list
  CTOFnb 12.0 -        !max distance for pair to consider in energy
  CDIE -               !use constant dielectric/r
  EPS 1. -             !potential for constant ditlectric
  CTONnb 8. -          !Distance at which smoothing function reduces pair's contribution
  VSWItch -            !Use switching function for van der Waals to avoid cutoff discontinuity
  WMIN 1.0 -           !warning cutoff for minimum atom-atom distance
  CUTIm 16.0 -         !cutoff for including image atom
  IMGFrq -1 -          !frequency to check for atoms in Langevin region
  EWALd -              !use Ewald summation
  PMEW -               !use Particle Mesh Ewald method
  FFTX 48 -            !number of FFT grid points
  FFTY 48 -            !number of FFT grid points
  FFTZ 48 -            !number of FFT grid points
  KAPPa .34 -          !width of the Gaussian distribution
  SPLIne -             !cubic spline interpolation
  ORDEr 6 -            !Order of the B-spline interpolation
  IPRFrq 10 -          !frequency for avg and rms energy
  IEQFrq 0  -          !frquency to assign or scale velocities to match FINALT temperature
  NTRFrq 10 -          !frequency for stopping the rotation and translation of the system during dynamics
  IUNWri 41 -          !unit to write restart file
  IUNCrd 31 -          !unit to write coordinated (unformatted)
  IUNRea -1 -          !unit to read in restart file
  KUNIt -1 -           !unit to write temperature and total energy

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:

ICHEcw 0 TWINdh 0.0 TWINdl 0.0

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:

DYNA CPT STRT -
  TIMEstep 0.001 -     !timestep in picoseconds
  NSTEp 1000 -         !number of steps and energy evaluations
  ECHEck 500. -        !max variation of energy step-to-step
  PCONs -              !use constant pressure
  PINT -               !use internal virial for calc pressure
  PREF 1.0 -           !reference pressure in atm
  PMASs 20 -           !piston mass
  PGAMma 0.0 -         !Langevin piston collision frequency
  HOOVer reft -        !Hoover temperature control
  REFT 310.15 -        !Hoover reference temperature
  TMASs 190 -          !mass of thermal piston
  INBFrq -1 -          !lists updated when necessary (heuristic test)
  ATOM -               !electrostatics calculation (default) is done on atom-atom basis
  VATOm -              !needed for nonbond list generation
  CUTNb 16.0 -         !distance cutoff for neighbor list
  CTOFnb 12.0 -        !max distance for pair to consider in energy
  CDIE -               !use constant dielectric/r
  EPS 1. -             !potential for constant ditlectric
  CTONnb 8. -          !Distance at which smoothing function reduces pair's contribution
  VSWItch -            !Use switching function for van der Waals to avoid cutoff discontinuity
  WMIN 1.0 -           !warning cutoff for minimum atom-atom distance
  CUTIm 16.0 -         !cutoff for including image atom
  IMGFrq -1 -          !frequency to check for atoms in Langevin region
  EWALd -              !use Ewald summation
  PMEW -               !use Particle Mesh Ewald method
  FFTX 48 -            !number of FFT grid points
  FFTY 48 -            !number of FFT grid points
  FFTZ 48 -            !number of FFT grid points
  KAPPa .34 -          !width of the Gaussian distribution
  SPLIne -             !cubic spline interpolation
  ORDEr 6 -            !Order of the B-spline interpolation
  IPRFrq 10 -          !frequency for avg and rms energy
  IEQFrq 0  -          !frquency to assign or scale velocities to match FINALT temperature
  NTRFrq 10 -          !frequency for stopping the rotation and translation of the system during dynamics
  IUNWri 41 -          !unit to write restart file
  IUNCrd 31 -          !unit to write coordinated (unformatted)
  IUNRea -1 -          !unit to read in restart file
  KUNIt -1 -           !unit to write temperature and total energy
  nprint 100 -         !step frequency for writing in kunit and printing energy on unit 6
  nsavc 100 -          !frequency for writing coordinates
  nsavv 0 -            !frequency for writing velocities
  ihbfrq 0 -           !frequency to regenerate hydrogen bond list
  firstt 310.15 -      !initial temperature of the system
  finalt 310.15 -      !final temperature of the system
  ihtfrq 0 -           !frequency to increase temperature by TEMINC
  teminc 0.0 -         !value to increase temperature by
  tbath 310.15 -       !temperature of heath bath
  iasors 1 -           !assign (NOT scale) velocities during heating/equil
  iasvel 1 -           !use gaussian distribution of velocities
  iscvel 0 -           !single scale factor
  ichecw 0 -           !checks to make sure avg temperature lies within a window
  twindh 0.0 -         !Highest deviation allowed of FINALT on the high side
  twindl 0.0           !Highest deviation allowed of FINALT on the low side

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