CHARMm Principles



4       Performing Molecular Dynamics

 One of the most important computational techniques in CHARMm is molecular dynamics simulations. In CHARMm, molecular dynamics simulations are performed using a classical mechanics approach in which Newton's equations of motion are integrated for all atoms in the system. As with energy evaluations, a defined PSF, a set of coordinates, and a set of parameters are required to initiate a molecular dynamics run.

Molecular dynamics can be used to generate a realistic picture of a structure's motions, perform conformational searching, do a time series analysis of structural and energetic properties, explore energy decay, and analyze solvent effects. In addition to exploring more global aspects of macromolecular structure than can be accomplished with minimization techniques, molecular dynamics can help with understanding critical aspects of protein function that involve both small-scale and large-scale atomic movements.

This chapter explains

Script examples are included throughout the chapter.


Understanding molecular dynamics

The essence of the molecular dynamics technique is in the numerical integration of Newton's second law relating the mass and acceleration of an atom in the system to the gradient of the potential energy field. Given acceleration, an approximate velocity for the atom can be computed for a given period of time and changes in atomic coordinates can be determined.

Several methods are available for numerical integration. In CHARMm, the Verlet method is implemented. Leap frog and velocity Verlet variants are also available. For detailed discussion of these methods and other numerical methods for integrating differential equations, see Haile (1992) or Allen and Tildesley (1987).

Basic steps

A typical molecular dynamics run involves the following basic steps:

1.   Preliminary preparation -- A molecular structure with all Cartesian coordinates defined is required for a dynamics simulation. After determining the internal coordinate values of the molecule, total energy as a function of the Cartesian coordinates is computed by evaluating the individual terms of the energy equation.

2.   Minimization -- All dynamics simulations begin with an initial structure that may be derived from experimental data. Energy minimization is performed on structures prior to dynamics to relax the conformation and remove steric overlap that produces bad contacts. In the absence of an experimental structure, a minimized ideal geometry can be used as a starting point.

3.   Heating -- A minimized structure represents the molecule at a temperature close to absolute zero. Heating is accomplished by initially assigning random velocities according to a Gaussian distribution appropriate for that low temperature and then running dynamics. The temperature is gradually increased by assigning greater random velocities to each atom at predetermined time intervals.

4.   Equilibration -- Equilibration is achieved by allowing the system to evolve spontaneously for a period of time and integrating the equations of motion until the average temperature and structure remain stable. This is facilitated by periodically reassigning velocities appropriate to the desired temperature. Generally, the procedure is continued until various statistical properties of the system become independent of time.

5.   Production -- In the final molecular dynamics simulation, CHARMm takes the equilibrated structure as its starting point. In a typical simulation, the trajectory traces the motions of the molecule through a period of at least 10 picoseconds. Just as with energy minimization, provision is made to update the nonbonded and hydrogen bonded lists periodically. Additional options are available, making the dynamics facility quite flexible.

6.   Quenching -- The logical opposite of heating, this optional step takes the molecule from the equilibrated temperature to zero. Quenching is a form of minimization, utilizing molecular dynamics to slowly remove all kinetic energy from the system.

Strictly speaking, minimization and heating are not necessary, provided the equilibration process is long enough. However, these steps can serve as a means to arrive at an equilibrated structure in an effective way.

A molecular dynamics run generates a dynamics trajectory consisting of a set of frames of coordinates and velocities that represent the trajectory of the atoms over time. Using trajectory data, you can compute the average structure and analyze fluctuations of geometric parameters, thermodynamics properties, and time-dependent processes of the molecule.

Preliminary analysis is possible using commands provided in the coordinate manipulation facility. Gross changes, as well as more detailed perturbations, can be monitored using correlation functions.

Because molecular dynamics runs often require considerable amounts of computer time, a restart facility is available that allows you to suspend the simulation and resume the calculation.

Temperature

The instantaneous kinetic temperature (T) at which each dynamics step is defined is directly proportional to the kinetic energy of the molecular system of N atoms and given by the equation:

Eq. 8     

Where:

mi and vi are the mass and velocity of the ith atom

kb is the Boltzmann constant

Average kinetic energy is proportional to the kinetic temperature:

Where:

Dynamics time step

Because the potential energy field changes as a molecular structure is adjusted, accelerations must also change. Consequently, the time step in molecular dynamics must be very small. The smaller the time step, the better the approximation.

The choice of the integration time step is very important. The time step should be chosen to be at least an order of magnitude smaller than the length of the time corresponding to the fastest motion in a molecule. Typically, this motion is the vibration of a bond containing a hydrogen (with a period less than 0.001 ps) that would require a time step of 0.0005 ps.

Another factor to be considered is the simulation time necessary to perform dynamics simulations. This time can range from ten to hundreds of picoseconds or even nanoseconds. The characteristic times for common molecular events are listed in the following table:
Event Time
Bond stretch   ~1 to 20 fs  
Elastic domain modes   100 fs to several ps  
Water reorientation   4 ps  
Inter-domain bending   10 ps to 100 ns  
Globular protein tumbling   1 to 10 ns  
Aromatic ring flipping   100 µs to several sec  
Allosteric shifts   2 µs to several sec  
Local denaturation   1 ms to several sec  

To allow larger time steps and shorter real-time dynamics simulations, the energy and force of certain high frequency motions such as bond stretching may be removed by constraining them to their initial or equilibrium values. This is a reasonable approximation, because it is assumed that the harmonic stretching of bonds is restricted by a very large force constant and, particularly for bonds containing hydrogens, is not significantly coupled to other motions of the system.

The method used in CHARMm to constrain these motions is the SHAKE algorithm. SHAKE may be applied to all bonds or only bonds containing hydrogens, and to all angles or only angles containing hydrogens. To maintain a high degree of accuracy during the dynamics simulation, as well as to provide a larger integration time step, apply SHAKE only to bonds containing hydrogens. A two-fold increase in the time step (to 0.001 ps) is then allowed.

Length of trajectory

The length of a trajectory needed to calculate a property depends on the time variation of the property under consideration. If the property is a slowly varying function, the dynamics integration should be extended to cover several periods.

If a property varies randomly about a mean value with a decay time of t, and the simulation is run for length T, the variance in the estimate will be proportional to the square root of t/T. When multiple independent samples are present (for example, each water molecule in a solvent simulation), averaging can be done over the samples to improve sampling.


Running dynamic simulation variants

Several variants of the standard molecular dynamics simulations are available in CHARMm and are described briefly in the following sections.

Langevin dynamics

The Langevin dynamics method (McCammon et al. 1976), (Levy et al. 1979) approximates a full molecular dynamics simulation of a system by eliminating unimportant or uninteresting degrees of freedom. The effects of the eliminated degrees of freedom are simulated by mean and stochastic forces. For example, instead of simulating hundreds of solvent molecules surrounding the solute molecules, the solvent can be ideally represented by a viscous fluid described in terms of dissipative and fluctuative equations.

Stochastic boundary molecular dynamics

The stochastic boundary molecular dynamics method uses a combination of both Langevin dynamics and Newtonian dynamics. With this method, the molecular system is partitioned into a reaction region where Newtonian dynamics simulation is run, a buffer region where Langevin dynamics is run, and a reservoir region.

In this way, atoms distant from the specific interactive sites in a large macromolecular system can be effectively eliminated from extensive analysis. This allows detailed studies of spatially localized portions of interacting molecular systems. Enzyme-substrate active site interactions can be effectively studied using this technique.

Constant pressure and temperature dynamics

CHARMm supports the use of constant pressure and constant temperature algorithms in molecular dynamics simulations. Only a few changes are needed to the standard CHARMm dynamics input file to run a constant pressure/temperature simulation. These changes include the specification of the isothermal compressibility (for solvated systems, a value of 0.0000463 atm.-1 is suggested), temperature and pressure coupling constants (5.0 ps for each), and a reference temperature and pressure (298.0 K and 1.0 atm., respectively).

For constant pressure dynamics, you must use the CHARMm crystal facility. Constant pressure scaling only works with cubic, orthorhombic, or triclinic unit cells. Shape matrix propagation and coordinates scaling for triclinic unit cells are done according to Brown and Clark (1991).

Pressure and pressure statistics can be evaluated for static coordinate frames. The external isotropic pressure and tensor are calculated if a volume is present and a temperature is given. Average pressure and fluctuations are also computed.

Example: Molecular dynamics using the Berendsen CPT algorithm

This example sets up a crystal of the alanine tetrapeptide and performs a constant temperature and pressure dynamics simulation.


* Molecular dynamics using the Berendsen CPT algorithm.
*
stream datadir.def
Read rtf card...
<read topology file>
.
.
.
Read sequence card
<Read sequence>
.
.
.
Generate ALA4 setup
<read coordinates>
.
.
.
Read coor card
* OPTIMISED COORDINATES FOR THE FULL UNIT CELL OF CRYSTALLINE
* ALANINE (4 MOLECULES). THE LATTICE IS OPTIMISED WITH-
* A = 5.59967,
* B = 12.19617 AND C = 5.40430 ANGSTROMS.
<Read coordinates>
.
.
.
! Save the coordinates.
Coor copy comp
! Define the crystal.
Crystal Define orthorhombic 5.59967 12.19617 5.40430 90.0 -90.0 90.0
! Read in an existing transformation file.
Crystal Read card
<read crystal file>
.
.
.
! Calculate an energy with an update.
Energy imgfrq 10 inbfrq 10 ihbfrq 0 cutim 999.0
! Do 100 steps of constant temperature dynamics.
Dynamics CPT TConstant TCouple 0.4 TReference 300.0 -
                  strt nstep 100 timestep 0.001 -
                  iprfrq 100 ihtfrq 0 ieqfrq 0 ntrfrq 0  -
                  nprint 1 nsavc 0 nsavv 0 -
                  firstt 300.0 finalt 300.0
! Do a further 100 steps of constant pressure and temperature
! dynamics and test restart
open write card unit 1 name @9xtlala3.rst
Dynamics CPT TConstant TCouple 0.4 TReference 300.0 -
             PConstant Compressibility 0.00005 PCouple 0.1 -
             PReference 1.0 -
                  strt nstep 50 timestep 0.001 -
                  iprfrq 100 ihtfrq 0 ieqfrq 0 ntrfrq 0  -
                  nprint 1 nsavc 0 nsavv 0 iunwri 1 -
                  firstt 300.0 finalt 300.0
! now test restart; dynamics should use crystal data from the 
! restart file
!
open read card unit 1 name @9xtlala3.rst
Crystal Define orthorhombic 6. 12. 5. 90.0 90.0 90.0
Dynamics rest CPT TConstant TCouple 0.4 TReference 300.0 -
             PConstant Compressibility 0.00005 PCouple 0.1 -
             PReference 1.0 -
                  strt nstep 50 timestep 0.001 -
                  iprfrq 100 ihtfrq 0 ieqfrq 0 ntrfrq 0  -
                  nprint 1 nsavc 0 nsavv 0 iunrea 1 -
                  firstt 300.0 finalt 300.0
! Compare the coordinates.
Coordinates RMS Mass
Print Coordinates
Stop

Quenched molecular dynamics

 Quenched molecular dynamics can be run in at least two different ways, as illustrated in the following short examples.

Using Langevin dynamics

The following CHARMm script fragment could be used for quenched Langevin dynamics:


dyna langevin nstep 1000 TBATH 1000. npri 100 
dyna langevin nstep 1000 TBATH 500. npri 100 
dyna langevin nstep 1000 TBATH 500. npri 100
dyna langevin nstep 1000 TBATH 250. npri 100
dyna langevin nstep 1000 TBATH 125. npri 100
dyna langevin nstep 1000 TBATH 0. npri 100
Using Newtonian dynamics

The following CHARMm script fragment could be used for quenched Newtonian dynamics with velocity scaling:


dyna nstep 5000 FIRStt 1000. FINAlt 0. TEMInc -20. IHTFrq 100 -
	npri 100


Monitoring dynamics simulation output

 Dynamics simulations output a large amount of data:

Certain terms and values written to the output file should be monitored closely. After the molecular system has been equilibrated, a short trajectory of free dynamics should be run to determine the stability of the energy and the temperature.

Because the system has no external forces, conservation of energy must be maintained. Check this by comparing the starting total energy term to the last total energy term. A ratio of the RMS fluctuation of the total energy to the RMS fluctuation of the kinetic energy for the free dynamics trajectory is an excellent measure of the accuracy of the run. A value of 0.001 or less is expected.

The average temperature and its fluctuation, printed at the end of the trajectory output, should also be noted. A temperature near the one specified by FINALT with a relatively small fluctuation is required. If a constant temperature cannot be maintained, the system may have heated too rapidly or may not have completely equilibrated.

Coordinate trajectory format

The coordinate trajectories produced by CHARMm are always in binary format. The principal reason for this is that binary files require much less disk space than their text file counterparts Also, reading binary data is much faster than reading text data. Because trajectories are used extensively during analysis, the speed with which they may be read is an important factor.

Binary formats, however, are machine dependent and cannot be easily moved from one machine to another. To overcome this limitation, CHARMm provides a way to format dynamics trajectories in a machine-independent manner. The resulting coordinate trajectory files can be easily copied from one machine to another. CHARMm is then used to return the files to the standard binary version to be used for analysis or graphic display. For this procedure to work properly, you must have CHARMm installed on both machines.

The format of the binary trajectory (DCD) file is illustrated in the following example:


HDR,ICNTRL ! character*4 HDR, integer ICNTRL(20)NTITL,(TITLE(J),J=1,NTITL) ! integer NTITL, CHARACTER*80 TITLE(*)
NATOM  ! integer NATOM



        FREEAT(I),I=1,NFREAT)      ! INTEGER FREEAT(*)



XTLABC  ! REAL*8 XTLABC(6)
(X(I),I=1,NATOM) ! real X(NATOM)
(Y(I),I=1,NATOM) ! real Y(NATOM)
(Z(I),I=1,NATOM) ! real Z(NATOM)
For all subsequent coordinate sets, only movable atoms are written. 


XTLABC
(X(FREEAT(I)),I=1,NFREAT)
(Y(FREEAT(I)),I=1,NFREAT)
(Z(FREEAT(I)),I=1,NFREAT)
Where HDR = `CORD' or `VELD' for coordinates and velocities, respectively: 


ICNTRL(1)=NFILE ! number of frames in trajectory file
ICNTRL(2)=NPRIV ! number of steps in previous run
ICNTRL(3)=NSAVC ! frequency of saving
ICNTRL(4)=NSTEP ! total number of steps 
NFILE=NSTEP/NSAVC
ICNTRL(8)=NDEGF ! number of degrees of freedom
ICNTRL(9)=NATOM-NFREAT ! number of fixed atoms
ICNTRL(10)=DELTA4 ! coded time step
ICNTRL(11)=stoi(XTLTYP,ALPHABET) ! coded crystallographic 
! group (or zero)
ICNTRL(20)=VERNUM ! version number
All other entries (12-19) are zero.


Using time sets and correlation functions

 After running a dynamics simulation on a molecule, you may want to extract data on the extent and nature of the changes that have occurred in the system. CHARMm can calculate a time series for an extensive set of properties of the total system trajectory or for some of its components.

The properties for which a time series can be calculated are averages over a set of individual properties of the molecule. Some of the individual properties encompass calculations on the position of an atom, the energy or geometry of internal coordinates, the vector joining two atoms, the scalar product of two such vectors, the scalar product of the position vectors of two atoms, fluctuations of the positions of two atoms, the scalar product of the functions of two atoms, or the velocity/kinetic energy of an atom.

The time series can also be evaluated for certain global properties like the radius of gyration and the number density of the molecule. The ability to average vectors or internal coordinates over a set of atoms is helpful in defining and studying the behavior of local regions. If the requested correlation is a cross-correlation between two different sets, then both sets are evaluated simultaneously and stored.

A time series can be manipulated, saved, plotted, or used to generate correlation functions. Examples of experimental quantities that can be calculated using correlation functions are frictional coefficients, IR line widths, fluorescence depolarization rates, spectral densities (the Fourier transform of the correlation function), and NMR relaxation times.

Alternatively, a covariance matrix can be computed for a collection of time series. This option computes a full matrix to use in entropy calculations or in other applications. For this option, only one TRAJectory command is allowed.

A typical protocol for using time series and correlation functions in CHARMm takes the following basic form:

1.   Read coordinate or velocity sets of a dynamics trajectory. Evaluate the specified property for each set, then store the data.

2.   Manipulate the time series (Qt) to evaluate certain functions of the time series.

3.   Integrate the resultant time series (Ft) subset to obtain correlation functions.

4.   Manipulate the correlation function (Ct) to determine some function of the correlation function (Gt).

5.   Analyze and plot either the time series (Ft) or the correlation functions (Ct) during any of these processes.

When the correlation function is calculated, a number of operations can be performed on it. The integral, the integral of the square, and the logarithm of the correlation function can be evaluated, and the Fourier transform of the correlation function can be taken to obtain the spectral density.


Example: Running enkephalin dynamics

A very simple molecular dynamics simulation is performed in this example. Following the simulation, a few of the CHARMm dynamic analysis commands are used to compute average atomic positions and fluctuations.

The following files are used or generated in this example:

ENKPA(BC).RST

ENKPA(BC).DCD

ENKPA(BC).DVL

ENKPA(BC).ENE

ENKPH.CRD

ENKPE.CRD

ENKPD.CRD

ENKPAVG.CRD

ENKPDIF.CRD

AMINO.BIN

PARM.BIN

ENKPMIN.CRD 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
*...
*...
* This CHARMm command file generates a dynamics trajectory
* for met-enkephalin, from which the average structure and
* fluctuations are computed.
*...
*
!
! Open and read an all hydrogen amino acid topology file
! and parameter file
OPEN READ UNIT 02 FILE NAME "$CHM_DATA/AMINOH.BIN"
READ RTF UNIT 02 FILE
CLOSE UNIT 01
OPEN READ UNIT 02 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 02 FILE
CLOSE UNIT 02
! Read the sequence for met-enkephalin
READ SEQUENCE CARD
* Pentapeptide sequence for met-enkephalin
*
  5
TYR GLY GLY PHE MET
! Generate the PSF with a segment identifier of ENKP.
! Set up the internal coordinate tables.
GENERATE ENKP SETUP
! Read in minimized coordinates 
OPEN READ UNIT 04 CARD NAME "ENKMIN.CRD"
READ COOR UNIT 04 CARD
CLOSE UNIT 04
! Shake all bonds with hydrogens
SHAKE BONH
! Open the files necessary for the dynamics trajectory:
!   A restart file
!   A coordinates file (the coordinate trajectory)
!   A velocities file (the velocity trajectory)
!   An energy file!
OPEN WRITE UNIT 31 CARD NAME ENKPA.RST
OPEN WRITE UNIT 32 FILE NAME ENKPA.DCD
OPEN WRITE UNIT 33 FILE NAME ENKPA.DVL
OPEN WRITE UNIT 34 CARD NAME ENKPA.ENE
DYNA STRT VERLET NSTEP 3000 TIME 0.001 RDIE -
   VSWITCH SWITCH -
   IPRFRQ 100 IHTFRQ 50 IEQFRQ 0 INBFRQ -1 IHBFRQ 0 -
   IUNREA -1 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT 34 -
   NPRINT 50 NSAVC 50 NSAVV 50 -
   FIRSTT 0.0 FINALT 300.0 TEMINC 5 -
   TWINDH 10.0 TWINDL -10.0 -
   IASORS 1 IASVEL 1 ICHECW 0
OPEN WRITE UNIT 41 CARD NAME ENKPH.CRD
WRITE COOROORDINATES UNIT 41 CARD
* COORDINATES AFTER HEATING
*
CLOSE UNIT 31
CLOSE UNIT 32
CLOSE UNIT 33
CLOSE UNIT 34
CLOSE UNIT 41
OPEN READ UNIT 30 CARD NAME ENKPA.RST
OPEN WRITE UNIT 31 CARD NAME ENKPB.RST
OPEN WRITE UNIT 32 FILE NAME ENKPB.DCD
OPEN WRITE UNIT 33 FILE NAME ENKPB.DVL
OPEN WRITE UNIT 34 CARD NAME ENKPB.ENE
DYNA RESTART VERLET NSTEP 2000 TIME 0.001 RDIE -
   VSWITHC SWITCH -
   IPRFRQ 100 IHTFRQ 0 IEQFRQ 200 INBFRQ -1 IHBFRQ 0 -
   IUNREA 30 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT 34 -
   NPRINT 50 NSAVC 50 NSAVV 50 -
   FIRSTT 300.0 FINALT 300.0 -
   TWINDH 10.0 TWINDL -10.0 -
   IASORS 0 ISCVEL 0 ICHECW 1
OPEN WRITE UNIT 41 CARD NAME ENKPE.CRD
WRITE COOR UNIT 41 CARD
* COORDINATES AFTER EQUILIBRATION
*
CLOSE UNIT 30
CLOSE UNIT 31
CLOSE UNIT 32
CLOSE UNIT 33
CLOSE UNIT 34
CLOSE UNIT 41
OPEN READ UNIT 30 CARD NAME ENKPB.RST
OPEN WRIT UNIT 31 CARD NAME ENKPC.RST
OPEN WRIT UNIT 32 FILE NAME ENKPC.DCD
OPEN WRIT UNIT 33 FILE NAME ENKPC.DVL
OPEN WRIT UNIT 34 CARD NAME ENKPC.ENE
DYNA RESTART VERLET NSTEP 5000 TIME 0.001 RDIE -
   VSWITCH SWITCH -
   IPRFRQ 100 IHTFRQ 0  IEQFRQ 0 INBFRQ -1 IHBFRQ 0 -
   IUNREA 30 IUNWRI 31 IUNCRD 32 IUNVEL 33 KUNIT 34 -
   NPRINT 50 NSAVC 50 NSAVV 50 -
   FIRSTT 300.0 FINALT 300.0 -
   TWINDH 10.0 TWINDL -10.0 -
   ICHECW 0
OPEN WRITE UNIT 41 CARD NAME ENKPD.CRD
WRITE COOR UNIT 41 CARD
* COORDINATES AFTER DYNAMICS
*
CLOSE UNIT 30
CLOSE UNIT 31
CLOSE UNIT 32
CLOSE UNIT 33
CLOSE UNIT 34
CLOSE UNIT 41
COOR INIT
! Calculate the average structure and isotropic
! fluctuations from the free dynamics trajectory
OPEN READ UNIT 51 FILE NAME ENKPC.DCD
COOR DYNAMIC FIRSTU 51
CLOSE UNIT 51
OPEN WRITE UNIT 17 CARD NAME ENKPAVG.CRD
WRITE COORDINATES CARD UNIT 17
* Average structure and isotropic fluctuations
*
CLOSE UNIT 17
! Print out the average positions and fluctuations
! of the backbone atoms
PRINT COORDINATES SELE ATOM * * N .OR. ATOM * * CA .OR. -
               ATOM * * C END
! Compare the minimized coordinates with the
! average coordinates; store the average
! coordinates in the comparison set
COOR COPY COMP
COOR INIT
OPEN READ UNIT 21 CARD NAME enkmin.crd
READ COOR UNIT 21 CARD
CLOSE UNIT 21
COOR DIFF
OPEN WRITE UNIT 22 CARD NAME ENKPDIF.CRD
WRITE COOR UNIT 22 CARD
* Difference in atomic positions between dynamics
* average coordinates and minimized coordinates
*
CLOSE UNIT 22
PRINT COORDINATES SELE ATOM * * N .OR. ATOM * * CA .OR. -
               ATOM * * C END
! Restore the average coordinates to the main set
COOR INIT
COOR COPY
! Fill the internal coordinate arrays with the average
! ic values from the free dynamics trajectory
IC FILL
OPEN READ UNIT 51 FILE NAME ENKPC.DCD
IC DYNA AVERAGE FIRSTU 51
CLOSE UNIT 51
PRINT IC
! Fill the internal coordinate arrays with the
! fluctuations to the ic averages from the free
! dynamics trajectory
OPEN READ UNIT 51 FILE NAME ENKPC.DCD
IC DYNA FLUCTUATIONS FIRSTU 51
CLOSE UNIT 51
PRINT IC
STOP


Example: Calculating correlations in enkephalin dynamics

 The trajectory from the enkephalin dynamics simulation of the previous example is used here to compute internal correlated motion of the pentapeptide.

The following files are used in this example:

AMINOH.BIN

PARM.BIN

ENKPC.DCD 


*...
* Copyright (c) 1994
* Molecular Simulations Inc.
* All Rights Reserved
*
! Open and read an all hydrogen amino acid topology file
! and parameter file
OPEN READ UNIT 02 FILE NAME "$CHM_DATA/AMINOH.BIN"
READ RTF UNIT 02 FILE
CLOSE UNIT 02
OPEN READ UNIT 02 FILE NAME "$CHM_DATA/PARM.BIN"
READ PARA UNIT 02 FILE
CLOSE UNIT 02
! Read the sequence for met-enkephalin
READ SEQUENCE CARD
* Pentapeptide sequence for met-enkephalin
*
 5
TYR GLY GLY PHE MET
! Generate the PSF with a segment identifier of ENKP.
! Set up the internal coordinate tables.
GENERATE ENKP SETUP
! Open the enkephalin dynamics file
OPEN READ UNIT 51 FILE NAME "$GUIDE0601/ENKPC.DCD"
! Enter the correlation facility; specify a maximum of
! 500 dynamics coordinate files and a maximum of 4
! time series (1 each for the dihedrals). The maximum
! number of atoms required is 20 (4 for each of the four
! dihedral time series atoms plus 1 for the codes value).
CORREL MAXTIM 500 MAXSERIES 10 MAXATOMS 20
! Time series for correlation of the phi and psi
! dihedrals of the second and third residues
ENTER PHI2 DIHE ENKP 1 C ENKP 2 N ENKP 2 CA ENKP 2 C GEOMETRY
ENTER PSI2 DIHE ENKP 2 N ENKP 2 CA ENKP 2 C ENKP 3 N GEOMETRY
ENTER PHI3 DIHE ENKP 2 C ENKP 3 N ENKP 3 CA ENKP 3 C GEOMETRY
ENTER PSI3 DIHE ENKP 3 N ENKP 3 CA ENKP 3 C ENKP 4 N GEOMETRY
! trajectory
TRAJECTORY FIRSTU 51 NUNIT 1 BEGIN 5050 STOP 10000 SKIP 50
! Report on statistics
SHOW ALL
! Apply correlation functions. First, use both the FFT
! and DIRECT methods to compute the correlation between
! phi and psi for the second residue (to check for
! accuracy). Then, apply the correlation function to
! all other combinations of phi and psi.
CORFUN PHI2 PSI2 FFT
CORFUN PHI2 PSI2 DIRECT
! Autocorrelation of the dihedrals
CORFUN PHI2 PHI2
CORFUN PSI2 PSI2
CORFUN PHI3 PHI3
CORFUN PSI3 PSI3
! Cross-correlation of the dihedrals
CORFUN PHI2 PSI3
CORFUN PSI2 PHI3
CORFUN PHI2 PHI3
CORFUN PSI2 PSI3
CORFUN PHI2 PSI2
CORFUN PHI3 PSI3
! Exit correlation facility
END
STOP


References

 Allen, M. P. and Tildesley, D. J., Computer Simulations of Liquids, Oxford University Press, Oxford, 1987.

Berendsen, H. J. C., Postma, J. P. M., van Gunsteren, A., DiNola, W. F., and Haak, J. R., J. Chem. Phys., 81(8) 3684 (1984).

Brooks III, C. L., Brunger, A. T., and Karplus, M., Biopolymers, 24 843 (1985).

Brown, D. and Clark, J. H. R., Comp. Phys. Comm., 62 360 (1991).

Haile, J. M., Molecular Dynamic Simulations: Elementary Methods, John Wiley & Sons, New York, 1992.

Levy, R. M., McCammon, J. A., and Karplus, M., Chem. Phys. Lett., 64 4 (1979).

McCammon, J. A., Gelin, B. R., Karplus, M., and Wolynes, P. G., Nature (London), 262 325-326 (1976).

CHARMm Principles



Last updated September 18, 1998 at 04:02PM PDT.
Copyright © 1997, Molecular Simulations, Inc. All rights reserved.