Amber Workshop - Tutorial 2

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2004

In the previous section we saw the problems associated with running Molecular Dynamics in vacuo, namely that for biomolecular systems, which exist in a solvated environment, it does not accurately represent the system. One solution to this is to use explicit solvent, which is the subject of section 5. Using explicit solvent, however, can be expensive computationally and thus it is sometimes better to include solvent effects within the force field equations. This approach is termed implicit solvent. An excellent implicit solvent method currently implemented in Sander v8.0 is the Born solvation model [V. Tsui & D.A. Case, Biopolymers (Nucl. Acid. Sci.) 56, 275-291 (2001)].

Use of the Born solvation model is controlled by the IGB flag. The default value for this flag is 0 which corresponds to no generalised Born term being used. Alternative values are IGB=1 which corresponds to the "standard" pairwise generalised Born model using the default radii set up by LEaP (the method we will be using here), IGB=2 a modified GB model developed by A. Onufriev, D. Bashford and D.A.Case and IGB=5 which is essentially the same as IGB=2 but with alternative parameters. For further details see the AMBER v8.0 manual.

In this example we shall re-run our in vacuo simulations using the standard pairwise Born solvation model.

The first stage is to minimise our initial system. We use the same prmtop and inpcrd files as previously (polyAT_vac.prmtop, polyAT_vac.inpcrd). This time around though we will modify our input file to turn on the generalised Born method (IGB=1). Here is what the input file looks like (polyAT_gb_init_min.in):

The first thing you should notice is that this takes considerably longer than the in vacuo minimisation. On a 3 GHz Pentium 4 machine it takes about 47.8 seconds (over 10 times longer than the vacuum simulation). Herein lies the problem with including solvent in simulations. As we will see though it is often an expense which cannot be avoided.

Lets take a quick look at the output file produced during the minimisation (polyAT_gb_init_min.out). Again you will see that the energy has dropped between the first and last steps:

In section 3.2 we used sander to run molecular dynamics of our in vacuo DNA 10-mer. We noted that the results were very susceptible to the long range cut off, with no cut off giving us an unravelling of the DNA and subsequent crash of the sander MD program. This instability was attributed to the in vacuo simulation not representing reality correctly. In this section we will re-run the two simulations using the Born implicit solvent models. The two input files we require are:

We are using the same settings as before, namely we turn off minimisation (IMIN=0). We disable periodicity (NTB=0) but this time set IGB=1 since we want to use the Born implicit solvent model. We will again write information to the output file and trajectory coordinates file every 100 steps (NTPR=100, NTWX=100). Two different values for the cutoff will be used, one run will be with a cutoff of 12 angstroms (CUT = 12.0) and one run will be without a cutoff (CUT = 999). For temperature regulation we will use the Langevin thermostat (NTT=3, GAMMA_LN=1) to maintain the temperature of our system at 300 K. We shall also set the initial and final temperatures to 300 K (TEMPI=300.0, TEMP0=300.0) which will mean our system's temperature should remain around 300 K. Again we will run the two examples for a total of 100,000 steps each (NSTLIM=100000) with a 1 fs time step (DT=0.001) giving simulation lengths of 100 ps (100,000 x 1fs). The two input files are shown below:

polyAT_gb_md1_12Acut.in	polyAT_gb_md1_nocut.in
10-mer DNA MD Generalised Born, 12 angstrom cut off &cntrl imin = 0, ntb = 0, igb = 1, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 12.0 /	10-mer DNA MD Generalised Born, infinite cut off &cntrl imin = 0, ntb = 0, igb = 1, ntpr = 100, ntwx = 100, ntt = 3, gamma_ln = 1.0, tempi = 300.0, temp0 = 300.0 nstlim = 100000, dt = 0.001, cut = 999 /

We now run the two jobs, remembering to set our starting structures to the minimised structure from 4.1 above (polyAT_gb_init_min.rst), using the commands:

Note: These will take considerably longer to run than in vacuo MD did, so it is probably a good idea to just test running these calculations on your machine and then end them prematurely and simply use the output files given below. The 12 angstrom cutoff run takes about 8940s (2 hours 29 mins ) on a single 3 GHz Pentium 4 machine. You can follow the progress of the job by following the output file with the following command:

Warning: the mdcrd files are over 15MB in size - you can download gzip compressed (5.3MB) versions here (polyAT_gb_md1_12Acut.mdcrd.gz, polyAT_gb_md1_nocut.mdcrd.gz). You will need to unzip them before use. e.g gunzip polyAT_gb_md1_12Acut.mdcrd.gz

The first thing you should notice about the two generalised Born simulations is that they both completed successfully. The inclusion of solvent in the model has stabilised the DNA. Lets compare the potential energies of the two runs, 12 angstrom cutoff and no cutoff. We will use the process_mdout.perl perl script.

The first thing you should notice is that the plots are now very similar. The cut off value has made much less of a difference than it did in vacuo. Lets see what difference there is in the RMSd:

Notice how both simulations are much more stable than the vacuum case. You should also observe how similar the RMSd fits are. Thus we see that the use of a cut off here has far less influence on the results. It does, however, make a big difference to the time required for the simulation, especially when we use explicit solvent as we will find in the next section. On a single processor 3 GHz Pentium 4 the no cutoff simulation took 11,540 s while the 12 angstrom cutoff simulation took 8931 s, an increase in efficiency of 22.6 %.

Note: the RMSd is still quite large, however, typically we would want our RMSd to be less than 1.5 - 2 angstoms. We will see in a minute that the use of explicit solvent and more importantly periodic boundary techniques can improve the simulation still further.

Finally, before we move to the next stage of this tutorial, lets take a look at a video of the 12 angstrom cutoff trajectory.

Run vmd (based on vmd 1.8.1):

vmd

Open our 12 angstrom cutoff trajectory file. Select:

File -> New Molecule

Browse for the polyAT_vac.prmtop file and then select Parm7.

Then hit 'Load'.

Now select the trajectory to load, in this case the 12 angstrom cutoff generalised Born trajectory (polyAT_gb_md1_12Acut.mdcrd).

Our simulation is in implicit solvent but is not a periodic boundary simulation so select CRD and NOT CRDBOX. When you hit "Load" again you should see all of the frames loaded into the main molecule window. 1,000 frames in total.

You should have a go at playing the video of the trajectory. Notice how the DNA holds its secondary structure much better than the in vacuo simulation. There is also considerably less movement in the extremities. This is the effect of the solvent. You should especially notice the huge change in the structure of the first two residue over the last 20 ps of the simulation. The direction of the chains coil has begun to invert. It is possible that this is a transformation from B-DNA to A-DNA. A much longer simulation would be required to check this, however.

Below is a series of snap shots of the structure vs time:

0 ps	10 ps	20 ps
30 ps	40 ps	50 ps
60 ps	70 ps	80 ps
90 ps	100 ps

Now we have seen what difference implicit solvent makes we can now move on to looking at the issues involved with minimising and equilibrating simulations in explicit solvent.