University of Kurdistan, Sanandaj logo

Department of Chemistry, Faculty of Science, University of Kurdistan

Molecular dynamics simulations using AMBER

After setting up the protein and parametrizing hetero groups (cf. AMBER-setup and AMBER-Parm pages) we can start to set up MD simulations.

1.     Protonate and solvate the structure

Run tleap as the following

tleap -s -f tleap.in

where the tleap.in file is as:

source leaprc.ff14SB

source leaprc.gaff

source leaprc.water.tip3p

addPath /lunarc/nobackup/projects/teobio/Mehdi/glxii/Parms

loadAmberPrep zn.in

loadAmberPrep xoh.in

loadAmberPrep o2x.in

loadAmberPrep slg.in

loadAmberPrep gsh.in

loadAmberPrep gbp.in

loadAmberParams zn.dat

loadAmberParams xoh.dat

loadAmberParams o2x.dat

loadAmberParams slg.dat

loadAmberParams gsh.dat

loadAmberParams gbp.dat

x=loadpdb 01-tided_up

solvateOct x TIP3PBOX 10

saveamberparm x prmtop prmcrd

savepdb x pdbout

quit

 

In this step, note the following points:

·        Add "TER" between the chains and hetero groups coordinates in the final PDB file (01-tided_up).

·        If Cys-Cys crosslinks exist in the PDB structure, add extra bonds in the tleap.in file. (e.g., bond x.22.SG x.65.SG; note that the atom names must be capital; also note that tleap renumbers the residues, especially if there are several subunits).

·        You need to change and sort the atom names of the hetero groups (column 3) as they appear in the .in files. For example, we edited the lines containing GBP as follows.

HETATM 4083  C1  GBP B 464      27.140  47.038  60.878  1.00 43.76           C

HETATM 4084  O1  GBP B 464      27.407  47.659  61.940  1.00 45.46           O

HETATM 4085  O2  GBP B 464      26.979  47.636  59.785  1.00 45.08           O

HETATM 4086  N1  GBP B 464      28.149  44.845  60.180  1.00 42.86           N

HETATM 4087  C2  GBP B 464      26.999  45.486  60.865  1.00 42.41           C

HETATM 4088  C3  GBP B 464      26.858  44.817  62.259  1.00 39.40           C

HETATM 4089  C4  GBP B 464      23.765  40.788  63.655  1.00 24.26           C

HETATM 4090  S1  GBP B 464      24.903  40.411  62.262  1.00 27.99           S

HETATM 4091  C5  GBP B 464      25.645  45.211  63.151  1.00 34.52           C

HETATM 4092  C6  GBP B 464      24.354  44.388  63.029  1.00 32.05           C

HETATM 4093  O3  GBP B 464      23.381  44.874  62.439  1.00 32.66           O

HETATM 4094  N2  GBP B 464      24.355  43.170  63.578  1.00 26.25           N

HETATM 4095  C7  GBP B 464      23.215  42.228  63.508  1.00 24.01           C

HETATM 4096  C8  GBP B 464      22.119  42.465  64.618  1.00 21.59           C

HETATM 4097  O4  GBP B 464      21.071  41.810  64.582  1.00 20.46           O

HETATM 4098  N3  GBP B 464      22.477  43.285  65.617  1.00 20.07           N

HETATM 4099  C9  GBP B 464      21.677  43.493  66.841  1.00 21.16           C

HETATM 4100  C10 GBP B 464      21.340  44.944  67.176  1.00 21.22           C

HETATM 4101  O5  GBP B 464      21.720  45.837  66.388  1.00 25.56           O

HETATM 4102  O6  GBP B 464      20.693  45.180  68.211  1.00 18.65           O

HETATM 4103  N4  GBP B 464      24.505  40.113  59.521  1.00 35.22           N

HETATM 4104  O7  GBP B 464      25.829  40.195  59.532  1.00 37.31           O

HETATM 4105  C11 GBP B 464      23.951  40.271  60.735  1.00 31.76           C

HETATM 4106  O8  GBP B 464      22.741  40.219  60.938  1.00 33.64           O

HETATM 4107  C12 GBP B 464      23.856  39.902  58.227  1.00 35.58           C

HETATM 4108  C13 GBP B 464      22.626  40.498  57.912  1.00 36.70           C

HETATM 4109  C14 GBP B 464      24.500  39.089  57.282  1.00 36.16           C

HETATM 4110  C15 GBP B 464      22.032  40.300  56.668  1.00 36.01           C

HETATM 4111  C16 GBP B 464      23.909  38.894  56.039  1.00 36.35           C

HETATM 4112  C17 GBP B 464      22.689  39.488  55.748  1.00 36.44           C

HETATM 4113  Br1 GBP B 464      21.893  39.204  54.066  1.00 35.21          BR

 

·        You may want to add counter ions into your system to neutralize it or add an ionic strength to it. If so, add the following to your tleap.in file after solvateOct.

addIonsRand x Na+ 96 Cl- 110 4

 And also, add the following at the beginning of the tleap.in file to read the ion parameters.

loadAmberParams frcmod.ionsjc_tip3p

 However, since it replaces any water molecule with a counter ion, the ions may end up inside the protein. This can be avoided by running tleap both without and with counter ions and saving pdbout files of both calculations. Then run changepdb, command fci (fix counter ions) on the pdbout file with counter ions (pdbout-ion), and read in the pdbout file without counter ions. The program overwrites the prmcrd file with counterions that are more than 4 Å from the protein (and also from each other).

 

changepdb <<EOF

pdbout-ion

fci

prmcrd-ion

pdbout

prmcrd

q

EOF

\mv prmtop-ion prmtop

 

2.     Treating metal centers

If there are metal centers in your protein, treat them as described in section 3 of the AMBER-Parm page)

 

3.     Run MD simulations

Run MD simulations with the following commands and the sander.in files listed below. The following uses the GPU-accelerated pmemd code of AMBER for the simulations. However, you can instead use the sander module of AmberTools if you have not purchased the AMBER program (replace pmemd.cuda with mpirun -bind-to core -np 18 sander.MPI). Where "18" is the number of CPU cores that I use in my simulations (you have to replace it with the number of your CPU cores), note that simulations with the sander module on CPUs are much slower than those with the pmemd.cuda on GPUs.

 

pmemd.cuda  -O -i sander.in1 -o sander.out1 -r mdrest1 -p prmtop -c prmcrd  -ref prmcrd

pmemd.cuda  -O -i sander.in2 -o sander.out2 -r mdrest2 -p prmtop -c mdrest1 -ref prmcrd

pmemd.cuda  -O -i sander.in3 -o sander.out3 -r mdrest3 -p prmtop -c mdrest2 -ref prmcrd

pmemd.cuda  -O -i sander.in4 -o sander.out4 -r mdrest4 -p prmtop -c mdrest3

pmemd.cuda  -O -i sander.in5 -o sander.out5 -r mdrest5 -p prmtop -c mdrest4 -x mdcrd5

 

The sander.in1 file

Restrained minimization

  &cntrl

   irest=0,ntx=1,

   imin=1,maxcyc=1000,drms=0.0001,ntmin=2,

   ntc=2,ntf=2,

   cut=8.0,

   ntpr=100,ntwx=0,ntwv=0,ntwe=0,

   ipol=0,igb=0,

   ntr=1,restraint_wt=100,

   restraintmask="!:WA= & !@H="

   nmropt=1

  &end

 

&wt type='END' &end

 DISANG=rst

 

The sander.in2 file

Restrained NVT equilibration

  &cntrl

   irest=0,ntx=1,ig=-1,

   nstlim=10000,dt=0.002,

   temp0=300.0,ntt=3,gamma_ln=2.0,

   ntc=2,ntf=2,

   cut=8.0,

   ntpr=500,ntwx=0,ntwv=0,ntwe=0,

   ntb=1,

   ipol=0,igb=0,

   ntr=1,restraint_wt=50,

   restraintmask="!:WA= & !@H="

   nmropt=1

  &end

 

 &wt type='END' &end

 DISANG=rst

 

The sander.in3 file

Restrained NPT equilibration

  &cntrl

   irest=1,ntx=5,

   nstlim=10000,dt=0.002,

   temp0=300.0,ntt=3,gamma_ln=2.0,

   ntc=2,ntf=2,

   cut=8.0,

   ntpr=500,ntwx=0,ntwv=0,ntwe=0,

   ntb=2,ntp=1,pres0=1.0,taup=1.0,

   ipol=0,igb=0,

   ntr=1,restraint_wt=50,

   restraintmask="!:WA= & !@H="

   nmropt=1

  &end

 

&wt type='END' &end

 DISANG=rst

 

The sander.in4 file

Equilibration, 1 ns NPT

  &cntrl

   irest=1,ntx=5,

   nstlim=500000,dt=0.002,

   temp0=300.0,ntt=3,gamma_ln=2.0,

   ntc=2,ntf=2,

   cut=8.0,

   ntpr=500,ntwx=0,ntwv=0,ntwe=0,

   ntb=2,ntp=1,pres0=1.0,taup=1.0,

   ipol=0,igb=0,

   ntr=0

   nmropt=1

  &end

 

&wt type='END' &end

 DISANG=rst

 

The sander.in5 file

Production 100 ns, sampling frequency 10 ps

  &cntrl

   irest=1,ntx=5,

   nstlim=50000000,dt=0.002,

   temp0=300.0,ntt=3,gamma_ln=2.0,

   ntc=2,ntf=2,

   cut=8.0,

   ntpr=5000,ntwx=5000,ntwv=0,ntwe=0,

   ntb=2,ntp=1,pres0=1.0,taup=1.0,

   ipol=0,igb=0,

   ntr=0

   nmropt=1

  &end

 

&wt type='END' &end

 DISANG=rst

 

The lines that are shown in boldface in the sander.in files are only necessary if you have a metal center in your protein (cf. section 2).