Computing Binding Enthalpy Values

(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 Andrew T. Fenley & Michael K. Gilson 2014

Section 2: Preparing the system

This section describes the steps required to prepare the systems you built in Section 1 for simulation in AMBER.

There are three general setup stages that come before the main production simulations (these are common with many other types of simulations):

Minimize system to relax bad contacts (jump here)
Heat the system to the target temperature in an NVT ensemble (jump here)
Equilibrate the system at the target temperature in an NPT ensemble (jump here)

Minimizing the system

The first stage involves minimizing the solute(s) and explicit water for each system. We will perform minimization stage in two steps: the first step restrains the atomic positions of the solute and only relaxes the water, and the second step releases the restraint and allows all atoms to relax. (The pure water simulation only needs one step of minimization since no solutes are present.)

The input files are shown below:

b2_min.in

B2: initial minimization solvent &cntrl imin = 1, maxcyc = 10000, ncyc = 500, ntb = 1, ntr = 1, cut = 10.0 / Hold the Guest (B2) Fixed 100.0 RES 1 END END

To run this simulation, we use the command:

pmemd.cuda -O -i b2_min.in -p b2.prmtop -c b2_iso.rst7 -o b2_min.out -r b2_min.rst7 -ref b2_iso.rst7

(If the computer you are using does not contain a GPU with pmemd.cuda installed, then use either pmemd or pmemd.MPI for all simulations.)

We will now minimize the other simulations containing a solute:

cb7_min.in

CB7: initial minimization solvent &cntrl imin = 1, maxcyc = 10000, ncyc = 500, ntb = 1, ntr = 1, cut = 10.0 / Hold the Host (CB7) Fixed 100.0 RES 1 END END

To run this simulation, we use the command:

pmemd.cuda -O -i cb7_min.in -p cb7.prmtop -c cb7_iso.rst7 -o cb7_min.out -r cb7_min.rst7 -ref cb7_iso.rst7

cb7_b2_min.in

CB7 and B2: initial minimization solvent &cntrl imin = 1, maxcyc = 10000, ncyc = 500, ntb = 1, ntr = 1, cut = 10.0 / Hold the Host (CB7) and Guest (B2) Fixed 100.0 RES 1 2 END END

To run this simulation, we use the command:

pmemd.cuda -O -i cb7_b2_min.in -p cb7_b2.prmtop -c cb7_b2_iso.rst7 -o cb7_b2_min.out -r cb7_b2_min.rst7 -ref cb7_b2_iso.rst7

For minimization of the pure water box and the second step of minimization for cb7, b2, and cb7_b2, we will remove any restraints applied to the system.

all_min.in

All atoms minimization &cntrl imin = 1, maxcyc = 10000, ncyc = 500, ntb = 1, ntr = 0, cut = 10.0 /

To run this simulation, we use the command:

pmemd.cuda -O -i all_min.in -p water.prmtop -c water_iso.rst7 -o water_min.out -r water_min.rst7

We will also complete the second step of minimization for the other systems:

pmemd.cuda -O -i all_min.in -p b2.prmtop -c b2_min.rst7 -o b2_min2.out -r b2_min2.rst7

pmemd.cuda -O -i all_min.in -p cb7.prmtop -c cb7_min.rst7 -o cb7_min2.out -r cb7_min2.rst7

pmemd.cuda -O -i all_min.in -p cb7_b2.prmtop -c cb7_b2_min.rst7 -o cb7_b2_min2.out -r cb7_b2_min2.rst7

All systems should be fully minimized at this point, and the following restart files will now begin the NVT equilibration stage: b2_min2.rst7, cb7_min2.rst7, cb7_b2_min2.rst7, and water_min.rst7. We recommend checking these structures in a visualization package to make sure the minimization process did not cause any issues with the systems.

Heating the system

In this step, we will heat the system to 300 K while keeping the volume of the box fixed.

nvt_md.in

Heating Equilibration - constant volume &cntrl imin = 0, irest = 0, ntx = 1, ntb = 1, cut = 10.0, ntr = 0, ntc = 2, ntf = 2, tempi = 0.0, temp0 = 300.0, ntt = 3, ig = -1, ioutfm = 1, iwrap = 1, gamma_ln = 1.0, nstlim = 1000000, dt = 0.002 ntpr = 1000, ntwx = 1000, ntwr = 1000 /

To run the heating simulations, you can use the following commands:

pmemd.cuda -O -i nvt_md.in -c water_min.rst7 -p water.prmtop -o water_nvt.out -r water_nvt.rst7 -x water_nvt.nc

pmemd.cuda -O -i nvt_md.in -c b2_min2.rst7 -p b2.prmtop -o b2_nvt.out -r b2_nvt.rst7 -x b2_nvt.nc

pmemd.cuda -O -i nvt_md.in -c cb7_min2.rst7 -p cb7.prmtop -o cb7_nvt.out -r cb7_nvt.rst7 -x cb7_nvt.nc

pmemd.cuda -O -i nvt_md.in -c cb7_b2_min2.rst7 -p cb7_b2.prmtop -o cb7_b2_nvt.out -r cb7_b2_nvt.rst7 -x cb7_b2_nvt.nc

After the simulations finish (each took ~15 minutes on a GTX 680 GPU), you should have the following files:

Output files
water_nvt.out
b2_nvt.out
cb7_nvt.out
cb7_b2_nvt.out

Trajectory files
water_nvt.nc
b2_nvt.nc
cb7_nvt.nc
cb7_b2_nvt.nc

Restart files
water_nvt.rst7
b2_nvt.rst7
cb7_nvt.rst7
cb7_b2_nvt.rst7

Equilibrating the system

Now that the systems are properly thermalized to 300 K, we need to set the pressure of the system to 1 bar. Therefore, we will now be running NPT simulations, where the box volume fluctuates, and the density of the system should reach a converged value (close to the density of water at standard conditions) for each simulation.

npt_md.in

Volume Equilibration MD - constant pressure &cntrl imin = 0, irest = 1, ntx = 5, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 10.0, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 5.0, ig = -1, ioutfm = 1, iwrap = 1, nstlim = 2500000, dt = 0.002, ntpr = 1000, ntwx = 1000, ntwr = 1000 / &ewald nfft1=32, nfft2=32, nfft3=32, order=4 /

Note, we are manually setting nfft1, nfft2, nfft3 all equal to 32 which is about the size of the box lengths for this small system. In general, one should be careful about these values.

To run the NPT simulations, you can use the following commands:

pmemd.cuda -O -i npt_md.in -c water_nvt.rst7 -p water.prmtop -o water_npt.out -r water_npt.rst7 -x water_npt.nc

pmemd.cuda -O -i npt_md.in -c b2_nvt.rst7 -p b2.prmtop -o b2_npt.out -r b2_npt.rst7 -x b2_npt.nc

pmemd.cuda -O -i npt_md.in -c cb7_nvt.rst7 -p cb7.prmtop -o cb7_npt.out -r cb7_npt.rst7 -x cb7_npt.nc

pmemd.cuda -O -i npt_md.in -c cb7_b2_nvt.rst7 -p cb7_b2.prmtop -o cb7_b2_npt.out -r cb7_b2_npt.rst7 -x cb7_b2_npt.nc

After the simulations finish (each took ~50 minutes on a GTX 680 GPU), you should have the following files:

Output files
water_npt.out
b2_npt.out
cb7_npt.out
cb7_b2_npt.out

Trajectory files (Note, these files are too large to provide online copies.)
water_npt.nc
b2_npt.nc
cb7_npt.nc
cb7_b2_npt.nc

Restart files
water_npt.rst7
b2_npt.rst7
cb7_npt.rst7
cb7_b2_npt.rst7

Now that the densities of the simulations have stabilized, we will run an additional NPT simulation that saves VOLUME information much more frequently. We will then use the average VOLUME as a way to assign the final box sizes of the systems during the production simulations. Note, we have also reduced the 'cut' value from 10.0 to 9.0 to match what will be used in the production simulation stage.

npt_md2.in

Volume Equilibration MD - constant pressure &cntrl imin = 0, irest = 1, ntx = 5, ntb = 2, pres0 = 1.0, ntp = 1, taup = 2.0, cut = 9.0, ntr = 0, ntc = 2, ntf = 2, tempi = 300.0, temp0 = 300.0, ntt = 3, gamma_ln = 5.0, ig = -1, ioutfm = 1, iwrap = 1, nstlim = 20000000, dt = 0.002, ntpr = 50, ntwx = 10000, ntwr = 10000 / &ewald nfft1=32, nfft2=32, nfft3=32, order=4 /

To run the second set of NPT simulations, you can use the following commands:

pmemd.cuda -O -i npt_md2.in -c water_npt.rst7 -p water.prmtop -o water_npt2.out -r water_npt2.rst7 -x water_npt2.nc

pmemd.cuda -O -i npt_md2.in -c b2_npt.rst7 -p b2.prmtop -o b2_npt2.out -r b2_npt2.rst7 -x b2_npt2.nc

pmemd.cuda -O -i npt_md2.in -c cb7_npt.rst7 -p cb7.prmtop -o cb7_npt2.out -r cb7_npt2.rst7 -x cb7_npt2.nc

pmemd.cuda -O -i npt_md2.in -c cb7_b2_npt.rst7 -p cb7_b2.prmtop -o cb7_b2_npt2.out -r cb7_b2_npt2.rst7 -x cb7_b2_npt2.nc

After the simulations finish (each took ~400 minutes on a GTX 680 GPU), you should have the following files:

Output files (Note, these files are too large to provide online copies.)
water_npt2.out
b2_npt2.out
cb7_npt2.out
cb7_b2_npt2.out

Trajectory files (Note, these files are too large to provide online copies.)
water_npt2.nc
b2_npt2.nc
cb7_npt2.nc
cb7_b2_npt2.nc

Restart files
water_npt2.rst7
b2_npt2.rst7
cb7_npt2.rst7
cb7_b2_npt2.rst7

We will now find the average box edge length by taking the cube-root of the average VOLUME for each simulation. This can be done via a simple combination of standard Unix/Linux tools: grep, head and awk.

grep VOL water_npt2.out | head -n -2 | awk '{sum+=$9}END{printf "%10.7f\n",(sum/NR)^(1/3)}'
grep VOL b2_npt2.out | head -n -2 | awk '{sum+=$9}END{printf "%10.7f\n",(sum/NR)^(1/3)}'
grep VOL cb7_npt2.out | head -n -2 | awk '{sum+=$9}END{printf "%10.7f\n",(sum/NR)^(1/3)}'
grep VOL cb7_b2_npt2.out | head -n -2 | awk '{sum+=$9}END{printf "%10.7f\n",(sum/NR)^(1/3)}'

Assuming everything has worked properly up until this point, the average box edges for each simulation should be nearly identical to the following values:

Water:      35.7337823
B2:         35.8011978
CB7:        35.9742065
CB7 and B2: 36.0325779

Finally, we need to modify the latest restart files such that the box lengths are set to the equilibrated average values above. Modify the *npt2.rst7 files and save them as *prod.rst7 (e.g. b2_prod.rst7).

Production Box Edges

==> b2_prod.rst7 <== 35.8011978 35.8011978 35.8011978 90.0000000 90.0000000 90.0000000 ==> cb7_prod.rst7 <== 35.9742065 35.9742065 35.9742065 90.0000000 90.0000000 90.0000000 ==> cb7_b2_prod.rst7 <== 36.0325779 36.0325779 36.0325779 90.0000000 90.0000000 90.0000000 ==> water_prod.rst7 <== 35.7337823 35.7337823 35.7337823 90.0000000 90.0000000 90.0000000

Everything is now ready to start production simulation runs. Now is also a good time to check the currently completed simulation trajectories to make sure there are no issues with the systems.

Click here to go to section 3

Click here to go back to section 1