AMOEBA

From UC_Chemistry

Revision as of 15:51, 24 April 2009 David Rogers (Talk | contribs) (new page) ← Previous diff		Current revision David Rogers (Talk | contribs) (→Simulation Setup with AMBER - formatting)
Line 18:		Line 18:
-	Since I was having some difficulties with AMBER's leaprc, I had to pair theirs down to make my own.	+	Since I was having some difficulties with AMBER's leaprc, I had to pair theirs down to make my own.
	''diff leaprc.sys $AMBERHOME/dat/leap/gleap/leaprc.amoeba''		''diff leaprc.sys $AMBERHOME/dat/leap/gleap/leaprc.amoeba''

---

Current revision

We have recently been doing some simulations using AMOEBA as implemented in AMBER10. The attached package Media:AMOEBA.tgz has the files needed to set up this calculation.

Simulation Setup with AMBER

The basic process is to use sleap < anion.leaprc

anion.leaprc

source leaprc.sys

sys=loadpdb anion.pdb
set sys box {18.761 18.761 18.761}
saveamoebaparm sys anion.prmtop anion.crd

quit

anion.pdb has to have a TER after every molecule (including water) or else sleap freaks out.

Since I was having some difficulties with AMBER's leaprc, I had to pair theirs down to make my own.

diff leaprc.sys $AMBERHOME/dat/leap/gleap/leaprc.amoeba

2,5d1
< loadoff amoeba_amino.off
< loadoff amoeba_aminont.off
< loadoff amoeba_aminoct.off
< loadoff amoeba_watbox.off
7c3,6
< loadamoebaparams amoebapro.prm
---
> loadoff wat.off
> loadoff sys.off
>
> loadamoebaparams amoeba.prm
76d74
<

The wat.off and sys.off stuff is important because the last line uses the plain old amoeba.prm file from TINKER's distribution -- which has a different atom type numbering from amoebapro.prm.

diff $AMBERHOME/dat/leap/gleap/amoeba_wat.off wat.off

4,6c4,6
<  "O" "202" 0 1 131072 1 8 -0.834000
<  "H1" "203" 0 1 131072 2 1 0.417000
<  "H2" "203" 0 1 131072 3 1 0.417000
---
>  "O" "22" 0 1 131072 1 8 -0.834000
>  "H1" "23" 0 1 131072 2 1 0.417000
>  "H2" "23" 0 1 131072 3 1 0.417000
8,10c8,10
<  "O" "202" 0 -1 0.0
<  "H1" "203" 0 -1 0.0
<  "H2" "203" 0 -1 0.0
---
>  "O" "22" 0 -1 0.0
>  "H1" "23" 0 -1 0.0
>  "H2" "23" 0 -1 0.0

Now sys.off contains the ion parameters and was generated using sleap <sys.cmd. There was some wierd error generating topologies when I tried to save all 4 anions to a single sys.off file, so I just uncommented the section I needed and re-made sys.off before building each anion.

That should be about it.

Simulation Runs

This is what sander was designed to do, so all I should have to tell you is that the command to run is

sander -O -i $prmdir/$job.in -o $locdir/$out/$job.log \
       -p $prmdir/anion.top -c $in -r restrt

where $in is the restart file, and that the sander parameter file looks like

equ.in

60 ps Equilib. Run
 &cntrl
  nstlim = 48000, dt = 0.00125,
  ntpr = 800, ntwe = 800, ntwx = 8000,
  ntx = 5, irest = 0,
  jfastw = 4, cut = 9.0,
  ntt = 1, tempi = 298.15, temp0 = 298.15, tautp = 0.2,
  ntp = 1, pres0 = 1.0, taup = 0.5, comp = 44.6,
  iamoeba=1, ntb = 2, iwrap = 1, nscm = 1,
 /
 &ewald
  order=8, ew_coeff=0.4, nfft1=32,nfft2=32,nfft3=32, skinnb = 0.1,
 /
 &amoeba
  dipole_scf_tol=1e-4, do_vdw_taper=1, do_vdw_longrange=1,
 /

or

run.in

2000 ps MD Run
 &cntrl
  nstlim = 1600000, dt = 0.00125, t = 0.0,
  ntpr = 400, ntwe = 400, ntwx = 400,
  ntx = 5, irest = 1, ioutfm = 1,
  jfastw = 4, cut = 9.0,
  ntt = 2, tempi = 298.15, temp0 = 298.15, vrand=200, ig=-1,
  ntp = 1, pres0 = 1.0, taup = 1.5, comp = 44.6,
  iamoeba=1, ntb = 2, iwrap = 1, nscm = 1,
 /
 &ewald
  order=8, ew_coeff=0.4, nfft1=32,nfft2=32,nfft3=32, skinnb = 0.1,
 /
 &amoeba
  dipole_scf_tol=1e-4, do_vdw_taper=1, do_vdw_longrange=1,
 /

However you will of course have to mess with little details like the ntx=5,irest=1 or establishing and using a consistent file location/naming scheme. The zipped info. contains such artifacts from my setup.

Simulation Analysis

We want to be able to get energy differences when we change the topology around. This can be done using sander's groupfile specification. As a word of warning, you will need to add a net charge correction of -332.0716 0.5 Q^2 / pi V ew_coeff^2 kCal/mol to the energies computed in this way if there is a net charge of Q on your system. Other than that, the required info. is contained in the package's fe.in and calc_de.job files. The enbrk.awk file is a handy little tool for parsing the .inf output files from sander, but you should always check that your files look the same way mine did or else you won't be pulling out the correct energies! It should be pretty easy to understand what it is doing, so there is no excuse for not tweaking it.