Setup the system with counter ions in a box of water.

Next, we turn from "vacuum" simulations to consider how to set up and carry out molecular dynamics simulations in a box of solvent water, using periodic boundary conditions. This example is typical of many molecular dynamics simulations begin carried out with AMBER.

Use LEaP to set up the prmtop file

This requires a modified version of the LEaP input file we used.  The file,  build_pcy, calls the solvateBox command as well as addIons:

+-----------------------------------------------------------+
| 			build_pcy			    |
+-----------------------------------------------------------+
|addAtomTypes {                                             |
|{"SM" "S"  "sp3"}                                          |
|}                                                          |
|loadamberparams frcmod.pcy                                 |
|loadAmberPrep mem.in                                       |
|set MEM restype protein                                    |
|loadAmberPrep hicu.in                                      |
|set HIC restype protein                                    |
|plc = loadPdb 1plc.protein.pdb                             |
|bond plc.87.ND1 plc.37.CU                                  |
|bond plc.84.SG plc.37.CU                                   |
|bond plc.92.SD plc.37.CU                                   |
|addIons plc Na+ 0                                          |
|solvateBox plc TIP3PBOX 10.0 0.8                           |
|saveAmberParm plc prmtop.wat prmcrd.wat                    |
+-----------------------------------------------------------+

The above file is read into LEaP as follows.   ">" is the prompt that LEaP provides:

tleap -f leaprc.ff94

      > source build_pcy

      > quit

Run solvated molecular dynamics simulation.

We will first run a simple minimization in water to remove initial bad contacts:

 +--------------------------------------------------------+
 |			  min2.in			  |
 +--------------------------------------------------------+
 |							  |
 |  minimization run					  |
 | &cntrl						  |
 |    imin=1, maxcyc=500,                                 |
 |    ntpr=25, cut=10.0,                                  |
 | /     						  |
 |							  |
 +--------------------------------------------------------+

Here is the command to run:

We then perform a  short constant volume run using Langevin dynamics.  The constant volume simulation is a necessary step in equilibrating the system and preparing for subsequent constant pressure molecular dynamics.   

+-----------------------------------------------------------------+
|        			 md1.in				  |
+-----------------------------------------------------------------+
|								  |
|molecular dynamics run 					  |
| &cntrl							  |
|    imin=0, irest=0, ntx=1, tempi=0.0,	   (no restart MD)	  |
|    ntt=3, temp0=300.0, tautp=1.0,	   (temperature control)  |
|    ntp=0, taup=2.0, nrespa=1,            (constant volume )     |
|    ntb=1, ntc=2, ntf=2, 		   (SHAKE, periodic bc.)  |
|    nstlim=40000, dt=0.001,               (run for 0.040 nsec)   |
|    ntwe=100, ntwx=100, ntpr=100, 	   (output frequency)	  |
|    gamma_ln=10,                          (collision frequency)  |
| / 								  | 
+-----------------------------------------------------------------+

We use the min2.rst file created above in the minimzation step to begin the dynamics with the following command:

sander -O -i md1.in -c min2.rst -p prmtop.wat
-o md1.out -r md1.rst -x md1.mdcrd -e md1.ene &

The output will include the md1.out file giving information about the progress of the trajectory, along with md1.mdcrd and md1.ene files that contain the coordinates and energy information at every 100th step, respectively.  Many of the analysis programs in AMBER can use these sorts of files.

Next, we restart the dynamics and perform constant pressure or "production" dynamics (md2.in):

+-----------------------------------------------------------------+
|    				 md2.in				  |
+-----------------------------------------------------------------+
|								  |
|molecular dynamics run 					  |
| &cntrl							  |
|    imin=0, irest=1, ntx=7, tempi=0.0,	   (restart MD)	          |
|    ntt=3, temp0=300.0, tautp=1.0,	   (temperature control)  |
|    ntp=1, taup=2.0, nrespa=1,            (constant pressure )   |
|    ntb=2, ntc=2, ntf=2, 		   (SHAKE, periodic bc.)  |
|    nstlim=60000, dt=0.001,               (run for 0.060 nsec)   |
|    ntwe=1000, ntwx=1000, ntpr=1000, 	   (output frequency)	  |
|    gamma_ln=10,                          (collision frequency)  | 
| / 								  |
|								  |
+-----------------------------------------------------------------+

Command to use:

    sander -O -i md2.in -c md1.rst -p prmtop.wat
    -o md2.out -r md2.rst -x md2.mdcrd -e md2.ene &



Coordinate and energy information is now updated every 1000th step of the trajectory.  Actual "production" computations would include many, many more MD steps than is given here, however, these examples will get you started.  *Note that following this step, the final density (after 100ps) and the average density are quite different (1.0046 vs 0.9367, respectively), suggesting that the simulation is not yet equilibrated.  It is good practice to break up production runs into manageable portions for analysis and to verify that the simulation is providing useful data.