Equilibration (part 2)

After the calculation is finished, you should have a couple more files in your directory.  Let's take a look at these results.

1. wt1mg_eq_v.rst <-- the restart file for our next calculation
2. wt1mg_eq_v.crd <-- saved snapshots
3. wt1mg_eq_v.out <-- output file containing various energy terms as well as temperature and what not.

Recalling that our goal for this part of equilibration was to raise the temperature of the system from 100K to 300K.  To examine the result, let's extract the temperature information from the output file and plot it on a graph.

On the command line, do the following:

$ grep TEMP wt1mg_eq_v.out | awk '{print $6, $9}' > temp.dat

This will extract the time and temperature information from the output file and save them in a new file called temp.dat.  You can then plot this data in your favorite graphing or spreadsheet program.  I used Excel to generate the following graph.  Remember to exclude the last two data items because they are the average and rms deviation, not the actual data points.

From this graph, we can see that the temperature of the system was gradually increased over the 10 ps simulation time.  It is not quite 300K yet, but for our purposes, it is close enough.

Next, we will equilibrate the system using pressure and temperature control to adjust the density of water to experimental values.

The input file is as follows:

Constant pressure constant temperature equilibration stage 2
  nstlim=5000, dt=0.002, ntx=5, irest=1, ntpr=500, ntwr=5000, ntwx=5000,

  temp0=300.0, ntt=1, tautp=2.0,
  ntb=2, ntp=1,

  ntc=2, ntf=2,

Notice that the new input file is almost identical to the previous one with exception of the highlighted items.  The first change is to reflect the fact that we are now continuing from a previous simulation and we will be using the velocity information from the restart file (ntx=5, irest=1).  We also switched on constant pressure (ntb=2) with isotropic position scaling (ntp=1) and gotten rid of unused parameters (tempi and ig).  For constant pressure simulatoins, we want to turn off the respa option.

Let's continue the fun:

$ sander -O -i eq_pt.in -p wt1mg.parm7 -c wt1mg_eq_v.rst -r wt1mg_eq_pt.rst -x wt1mg_eq_pt.crd -o wt1mg_eq_pt.out