(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 2015

Section 4: Running Minimization and MD (in implicit solvent)

Updated for AMBER 15

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 in this tutorial. Using explicit solvent, however, can be computationally expensive 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 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 generalized Born term being used. Alternative values are IGB=1 which corresponds to the "standard" pairwise generalized 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 manual.

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

4.1) Relaxing the System Prior to MD

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

So, now to run sander for minimization.

Input files: polyAT_gb_init_min.in, polyAT_vac.rst7, polyAT_vac.prmtop
Output files: polyAT_gb_init_min.out, polyAT_gb_init_min.ncrst

The first thing you should notice is that this takes considerably longer than the in vacuo minimization. On a 1.4 GHz i5 machine it takes about 22.8 seconds (about 10 times longer than the vacuum simulation). Herein lies the problem with including solvent in simulations. As we will see, it is often an expense that cannot be avoided.

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

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
      1      -2.2247E+03     1.9139E+01     9.3100E+01     C2'       634

   NSTEP       ENERGY          RMS            GMAX         NAME    NUMBER
    500      -3.3168E+03     8.7962E-01     6.8734E+00     N1        618

It is worth noting that the change in energy is significantly less than the in vacuo minimization. The starting energy in the Generalized Born case (-2224.7 kcal/mol) is significantly lower than in the in vacuo case (-1412.3 kcal/mol). This suggests that simply solvating our initial structure prior to minimization has significantly stabilized it. This stabilization will be even more apparent when we run molecular dynamics on the system.

Again we can create PDB files for the start (polyAT_vac.rst7) and final structures (polyAT_gb_init_min.ncrst). [polyAT_gb_init_min.pdb].

---

---

4.2) Running MD with generalized Born solvation

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 unraveling 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 minimization (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

As mentioned in the previous section when running production simulations with ntt=2 or 3 you should always change the random seed (ig) between restarts. You can achieve this automatically by setting ig=-1 in the ctrl namelist. We do not do this here in the tutorial for reproducibility but it is generally recommended that you do this in your calculations.

We now run the two jobs, remembering to set our starting structures to the minimized structure from 4.1 above (polyAT_gb_init_min.ncrst), 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 3654s (1 hours 1 min ) on a single 1.4 GHz Intel Core i5 machine. You can follow the progress of the job by following the output file with the following command:

Input files: polyAT_gb_md1_12Acut.in, polyAT_gb_md1_nocut.in, polyAT_gb_init_min.ncrst, polyAT_vac.prmtop

Output files: polyAT_gb_md1_12Acut.out, polyAT_gb_md1_nocut.out, polyAT_gb_md1_12Acut.ncrst, polyAT_gb_md1_nocut.ncrst, polyAT_gb_md1_12Acut.nc, polyAT_gb_md1_nocut.nc

S

4.3) Analyzing the results

The first thing you should notice about the two generalized Born simulations is that they both completed successfully. The inclusion of solvent in the model has stabilized the DNA. Lets compare the potential energies of the two runs, 12 angstrom cutoff and no cutoff. We will use the process_mdout.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:

Create the cpptraj input files:

polyAT_gb_md1_12Acut.calc_rms, polyAT_gb_md1_nocut.calc_rms.

Run the two cpptraj jobs with the following commands:

And then plot the RMSd fits:

---

---

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 cutoff 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 1.4 GHz Intel Core i5, the no cutoff simulation took 4,731 seconds while the 12 angstrom cutoff simulation took 3654 seconds, an increase in efficiency of 22.7 %. However, it is important to realize that the use of a cutoff here, due to the fact that there is no calculation of long range electrostatics as in an explicit solvent calculation, is still introducing additional approximations into the calculation. Thus if you can afford it, you should consider running implicit solvent GB simulations without a cutoff.

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.