The following is the full procedure for installing the AMBER force field port for GROMACS (AMBER-in-GROMACS, AMBER-with-GROMACS, AMBER-on-GROMACS, whatever you want to call it) developed by Eric Sorin at California State University, Long Beach, providing a bit more depth in the installation process (specifically for GROMACS 3.3.x) and a few modified GROMACS files.
As brief background, AMBER (Assisted Model Building and Energy Refinement) is one of THE dominant molecular mechanics/molecular dynamics (MM/MD) force fields used today in biochemical simulations. The motivation for this page (my installing AMBER for use in GROMACS) stems from the current Nanorex focus on Structural DNA Nanotechnology (SDN) modeling, for which we're working on a reduced model force field for large-structure energy minimizations and, importantly, integrating the GROMACS MM/MD package for use via our CAD interface. You can read more about this in the poster presented at FNANO08 this past April. As a force field validated for DNA simulations, AMBER meets our needs of performing atomistic simulations on DNA nanostructures. While NAMD is also a possibility for DNA simulations, GROMACS meets Nanorex's open source needs.
Needless to say, finding that Sorin and co-workers had ported AMBER to GROMACS was a wonderful discovery (and that this same port appears to be driving the Folding@Home project). The installation directions on the AMBER-on-GROMACS website are good, but they don't mention a few important steps that I spent no small amount of time trying to figure out (and I'm definitely not complaining. 99% of the work was the porting, which has been graciously handed to the GROMACS community on a silver platter). Below is the complete set of steps required to get, as the first example, the Dickerson.pdb sample file provided with the porting files to energy minimize and MD correctly. Issues related to DNA-based simulations not using the Dickerson.pdb file will be covered in another post.
QuteMol rendering of Dickerson.pdb (from Drew, H. R., Wing, R. M., Takano, T., Broka, C., Tanaka, S., Itakura, K. & Dickerson, R. E. (1981). Structure of a B-DNA dodecamer: conformation and dynamics. Proc. Natl. Acad. Sci. USA 78, 2179-2183. Protein Data Bank structure 1BNA).
Install/compile a custom installation of GROMACS
The AMBER port to GROMACS does not require modification of the GROMACS code but does require a few changes to force field text files. In the interest of remembering what was changed and what wasn't, I recommend a custom compilation of GROMACS or, at the very least, the installation of another copy of GROMACS you can modify the top directory of without risking changes to an already working GROMACS version.
If you don't know how to compile your own version of GROMACS, I recommend taking a look at Compiling Single-Precision And Double-Precision GROMACS 3.3.3 With OpenMPI 1.2.6 Under OSX 10.5 (Leopard).
My installation directory is /usr/local/gromacs333-amber (which I will refer to here as /ULG-A/ and all directory calls will be to this. Obviously, change all directory calls to match yours.
Download GROMACS-Ported AMBER force field files
I'm running a double-precision MPI version of GROMACS 3.3.3, meaning I'm using the GROMACS-ported AMBER force field files for GROMACS 3.3.1 (May 2006). And get the version "with pdfs." If you're going to be doing these calculations, you should know where the math come from!
Copy vdwradii.dat and aminoacids-NA.dat into /usr/local/gromacs-amber/share/gromacs/top
You will need to be logged in as root to do this (sudo cp FILENAME /ULG-A/top). The AMBER port file vdwradii.dat differs from the vdwradii.dat file in the /ULG-A/top directory of the GROMACS install by the addition of 5 lines:
??? P 0.15
??? LP1 0
??? LP2 0
LYSH MNZ1 0
LYSH MNZ2 0
With no removals or number changes.
The aminoacids.dat file is a little different. In order to use pdb2gmx to generate nucleic acid topology files, GROMACS requires the list of residue codes that reside in aminoacids-NA.dat. So far as I can tell, using aminoacids-NA.dat instead of aminoacids.dat does not change the handling of amino acids (why two files exist for a reason other than amino acid-centric organization is beyond my pay scale), so we'll be using aminoacids-NA.dat exclusively.
Move/delete/rename aminoacids.dat and rename aminoacids-NA.dat to aminoacids.dat
The aminoacids-NA.dat file includes 32 additional residue codes, accounting for the 3' (includes terminal H atom), 5' (included terminal H atom), in-strand (nucleic acid repeat unit) and molecule (individual hydrogen-terminated sugar and nucleic acid) topology information codes for DNA and RNA found in the .rtp files
Copy the ffamber files of interest into the /ULG-A/top directory
If you're doing this in Linux or OSX, you might as well move all of the files you want installed into a single directory and su cp * /ULG-A/top. Note that the .itp and .gro files in your ffamber_v3.3.1 (for the 3.3.x GROMACS) are the same for all of the AMBER force field flavors you can install. Further, all of the ffamber* files have unique names (94, 99, 03, etc.), so you can install all the force field files into /ULG-A/top and use whichever.
Modify /ULG-A/top/FF.dat to reflect the new force field files
The FF.dat file as installed by GROMACS looks like the following:
9
ffG43a1 GROMOS96 43a1 force field ffG43b1 GROMOS96 43b1 vacuum force field ffG43a2 GROMOS96 43a2 force field (improved alkane dihedrals) ffG45a3 GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) ffG53a5 GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656) ffG53a6 GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) ffoplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) ffencadv Encad all-atom force field, using scaled-down vacuum charges ffencads Encad all-atom force field, using full solvent charges
To include the AMBER force field files (so that these separate files can be called during pdb2gmx), we modify the FF.dat file by (1) incrementing the total number of force fields and (2) adding the force field code and text description to this file. The new FF.dat file will look like this if you added all of the ffamber force field files:
17
ffG43a1 GROMOS96 43a1 force field ffG43b1 GROMOS96 43b1 vacuum force field ffG43a2 GROMOS96 43a2 force field (improved alkane dihedrals) ffG45a3 GROMOS96 45a3 force field (Schuler JCC 2001 22 1205) ffG53a5 GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656) ffG53a6 GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656) ffoplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals) ffencadv Encad all-atom force field, using scaled-down vacuum charges ffencads Encad all-atom force field, using full solvent charges ffamber94 AMBER94 Cornell et al. (1995), JACS 117, 5179-5197 ffamber96 AMBER96 Kollman (1996), Acc. Chem. Res. 29, 461-469 ffamberGS AMBER99GS Garcia & Sanbonmatsu (2002), PNAS 99, 2782-2787 ffamberGSs AMBER99GSs Nymeyer & Garcia (2003) PNAS 100, 13934-13939 ffamber99 AMBER99 Wang et al. (2000), J. Comp. Chem. 21, 1049-1074 ffamber99p AMBER99p Sorin & Pande (2005). Biophys. J. 88(4), 2472-2493 ffamber99SB AMBER99sb Hornak et. al (2006). Proteins 65, 712-725 ffamber03 AMBER03 Duan et al. (2003), J. Comp. Chem. 24, 1999-2012
Add and increment numbers as appropriate for the force field you want to use. You can download the complete FF.dat at the following link and comment:
Download amber_FF.txt, delete "amber_", change ".txt" to ".dat", place in /ULG-A/top
source /UGL-A/bin/GMXRC
The GMXRC file contains path-dependent settings for your shell.
In theory, you're all done. At least, this is the end of the installation process according to the official AMBER-on-GROMACS list. The IMPORTANT NOTES section contains relevant information but does not complete the installation. The additional steps are provided below.
Modify the /ULG-A/top/spc.itp file
If you're running through a complete energy minimization calculation with the website installation followed, your first source of error (hopefully) comes in the form of…
Program grompp_amber, VERSION 3.3.3 Source code file: toppush.c, line: 1193 Fatal error: [ file "/usr/local/gromacs333_amber/share/gromacs/top/spc.itp", line 32 ]: Atom index (1) in bonds out of bounds (1-0). This probably means that you have inserted topology section "bonds" in a part belonging to a different molecule than you intended to. In that case move the "bonds" section to the right molecule.
The origin of this error is the lack of amberXX water parameters in the spc.itp file (nothing directly to do with the reported GROMACS error). The fix for this is straightforward, simply modifying the spc.itp file in the /ULG-A/top directory as follows below. The original spc.itp file below…
[ moleculetype ]
; molname nrexcl
SOL 2
[ atoms ]
; nr type resnr residue atom cgnr charge mass
#ifdef _FF_GROMACS
1 OW 1 SOL OW 1 -0.82
2 HW 1 SOL HW1 1 0.41
3 HW 1 SOL HW2 1 0.41
#endif
#ifdef _FF_GROMOS96
#ifdef HEAVY_H
1 OW 1 SOL OW 1 -0.82 9.95140
2 H 1 SOL HW1 1 0.41 4.03200
3 H 1 SOL HW2 1 0.41 4.03200
#else
1 OW 1 SOL OW 1 -0.82 15.99940
2 H 1 SOL HW1 1 0.41 1.00800
3 H 1 SOL HW2 1 0.41 1.00800
#endif
#endif
#ifdef _FF_OPLS
1 opls_116 1 SOL OW 1 -0.82
2 opls_117 1 SOL HW1 1 0.41
3 opls_117 1 SOL HW2 1 0.41
#endif
#ifdef FLEXIBLE
[ bonds ]
; i j funct length force.c.
1 2 1 0.1 345000 0.1 345000
1 3 1 0.1 345000 0.1 345000
[angles ]
; i j k funct angle force.c.
2 1 3 1 109.47 383 109.47 383
#else
[ settles ]
; OW funct doh dhh
1 1 0.1 0.16330
[ exclusions ]
1 2 3
2 1 3
3 1 2
#endif
is modified to this…
[ moleculetype ] ; molname nrexcl SOL 2 [ atoms ] ; nr type resnr residue atom cgnr charge mass #ifdef _FF_GROMACS 1 OW 1 SOL OW 1 -0.82 2 HW 1 SOL HW1 1 0.41 3 HW 1 SOL HW2 1 0.41 #endif #ifdef _FF_GROMOS96 #ifdef HEAVY_H 1 OW 1 SOL OW 1 -0.82 9.95140 2 H 1 SOL HW1 1 0.41 4.03200 3 H 1 SOL HW2 1 0.41 4.03200 #else 1 OW 1 SOL OW 1 -0.82 15.99940 2 H 1 SOL HW1 1 0.41 1.00800 3 H 1 SOL HW2 1 0.41 1.00800 #endif #endif #ifdef _FF_OPLS 1 opls_116 1 SOL OW 1 -0.82 2 opls_117 1 SOL HW1 1 0.41 3 opls_117 1 SOL HW2 1 0.41 #endif #ifdef _FF_AMBER94 ; also applies to FF_AMBER96, FF_AMBERGS, and FF_AMBERGSs 1 amber94_42 1 SOL OW 1 -0.82 15.99940 2 amber94_27 1 SOL HW1 1 0.41 1.00800 3 amber94_27 1 SOL HW2 1 0.41 1.00800 #endif #ifdef _FF_AMBER99 ; also applies to FF_AMBER99P, FF_AMBER99SB, FF_AMBER03 1 amber99_54 1 SOL OW 1 -0.82 15.99940 2 amber99_55 1 SOL HW1 1 0.41 1.00800 3 amber99_55 1 SOL HW2 1 0.41 1.00800 #endif #ifdef FLEXIBLE [ bonds ] ; i j funct length force.c. 1 2 1 0.1 345000 0.1 345000 1 3 1 0.1 345000 0.1 345000 [ angles ] ; i j k funct angle force.c. 2 1 3 1 109.47 383 109.47 383 #else [ settles ] ; OW funct doh dhh 1 1 0.1 0.16330 [ exclusions ] 1 2 3 2 1 3 3 1 2 #endif
This FFAMBER-included spc.itp file is downloadable in the following form and comment:
Download amber_spc.txt, delete "amber_" change ".txt" to ".itp", place in /ULG-A/top
At this point, pdb2gmx, editconf, genbox, grompp, and genion should all work without error. The next problem comes with the post-genion grompp, which can't find ion information…
Modify ions.itp
Similar to the spc.itp error, the error that comes up with the inclusion of counterions confounds the mind because everything seems to be in place.
Program grompp_amber, VERSION 3.3.3
Source code file: toppush.c, line: 1396
Fatal error:
No such moleculetype Na
With Na being whatever counterion you're trying. The AMBER-on-GROMACS website reports the following for AMBER ion inclusion:
Our AMBER ports include common ion definitions, which are listed in the ffamber*.rtp files (just below the TIP water models). This allows the AMBER ports to be used without modification or use of the GROMACS ions.itp file. At the moment these include Cl- , IB+, Na+, K+ , Rb+, Cs+, Li+ , Ca2+, Mg2+, Zn2+, Sr2+, and Ba2+. If your pdb file has ions present and pdb2gmx does not properly convert those ions, please check the atom and residue name of your ions and rename them if necessary to agree with the .rtp file. If you are using ion-related GROMACS tools, such as genion, you will need to enter the AMBER ion definition to the ions.itp file in the "top" directory of the GROMACS distribution.
The bases for the ion use are in the .rtp files, but the implementation, at least for GROMACS 3.3.3 and possibly earlier versions, is not as found in these files.
The proper format for both the FFAMBER_94 and FFAMBER_99 ion types in ions.itp is as follows:
#ifdef _FF_AMBER94
[ moleculetype ]
; molname nrexcl
Na 1
[ atoms ]
; id at type res nr residu name at name cgnr charge mass
1 amber94_31 1 Na Na 1 1 22.99000
#endif
#ifdef _FF_AMBER99
[ moleculetype ]
; molname nrexcl
Na 1
[ atoms ]
; id at type res nr residu name at name cgnr charge mass
1 amber99_31 1 Na Na 1 1 22.9900
#endif
This is the format used for the FF_GROMOS96 ions with one change. The atom type for the FF_AMBERXX parameters is NOT the atom label, but instead the AMBER force field label. In the case of Na, this is amber94_31 (for AMBER94) and amber99_31 (for AMBER99). Note that the mass must be specified (you have to scroll down the ions.itp file to see FF_GROMOS96. FF_GROMACS does not require mass specification).
For quick reference, here are the following [ atoms ] specifications for FF_AMBER94 and FF_AMBER_99.
#ifdef _FF_AMBER94 ; id at type res nr residu name at name cgnr charge mass 1 amber94_32 1 IB IB 1 1 131.00000 1 amber94_56 1 Li Li 1 1 6.94000 1 amber94_31 1 Na Na 1 1 22.99000 1 amber94_51 1 K K 1 1 39.10000 1 amber94_52 1 Rb Rb 1 1 85.47000 1 amber94_53 1 Cs Cs 1 1 132.90000 1 amber94_33 1 Mg Mg 1 2 24.30500 1 amber94_15 1 Ca Ca 1 2 40.08000 1 amber94_57 1 Zn Zn 1 2 65.40000 1 amber94_58 1 Sr Sr 1 2 87.62000 1 amber94_59 1 Ba Ba 1 2 137.33000 1 amber94_30 1 Cl Cl 1 -1 35.45000 #endif #ifdef _FF_AMBER99 ; id at type res nr residu name at name cgnr charge mass 1 amber99_32 1 IB IB 1 1 131.00000 1 amber99_56 1 Li Li 1 1 6.94000 1 amber99_31 1 Na Na 1 1 22.99000 1 amber99_51 1 K K 1 1 39.10000 1 amber99_52 1 Rb Rb 1 1 85.47000 1 amber99_53 1 Cs Cs 1 1 132.90000 1 amber99_33 1 Mg Mg 1 2 24.30500 1 amber99_15 1 Ca Ca 1 2 40.08000 1 amber99_57 1 Zn Zn 1 2 65.40000 1 amber99_58 1 Sr Sr 1 2 87.62000 1 amber99_59 1 Ba Ba 1 2 137.33000 1 amber99_30 1 Cl Cl 1 -1 35.45000 #endif
To save yourself plenty of trouble, download the ions.itp file with all of the ion parameters in the form of the following link:
Download amber_ions.txt, delete "amber_", change ".txt" to ".itp", place in /ULG-A/top
With the ions.itp file complete, your GROMACS energy minimizations and dynamics simulations of amino acid and nucleic acid structures should go without hitch.
My thanks to Alan Wilter S. da Silva, D.Sc. – CCPN Research Associate, Department of Biochemistry, University of Cambridge, for catching a formatting error in the first version of the amber_ions.itp file.
If you find problems, questions, concerns, incompatibilities, etc., please let me know by setting up a user account and posting a comment in this post or email me so I can post whatever you might find.
And don't forget to cite ffAMBER if you use it!
Sorin & Pande (2005). Biophys. J. 88(4), 2472-2493
chemistry.csulb.edu/ffamber
www.gromacs.org
chemistry.csulb.edu/esorin
chemistry.csulb.edu
amber.scripps.edu
www.nanorex.com
en.wikipedia.org/wiki/DNA_nanotechnology
www.somewhereville.com/rescv/fnano08_2008_poster.jpg
www.cs.duke.edu/~reif/FNANO
www.ks.uiuc.edu/Research/namd
en.wikipedia.org/wiki/Open_source
folding.stanford.edu
qutemol.sourceforge.net
www.pdb.org
www.pdb.org/pdb/explore/explore.do?structureId=1BNA
www.linux.org
www.apple.com/macosx
www.gromacs.org/documentation/reference/online/pdb2gmx.html
Quick Note: To test the Dickerson.pdb structure provided with the AMBER in GROMACS force field files, you can use the following set of commands:
pdb2gmx -nomerge -f Dickerson.pdb -o Dickerson_pdb2gmx.gro -p Dickerson_pdb2gmx.top
editconf -f Dickerson_pdb2gmx.gro -o Dickerson_editconf.gro -d 1.0 -bt triclinic
genbox -cp Dickerson_editconf.gro -cs -o Dickerson_genbox.gro -p Dickerson_pdb2gmx.top
grompp -f em -c Dickerson_genbox.gro -p Dickerson_pdb2gmx.top -o Dickerson_gromppem.tpr
genion -np 22 -norandom -pname Na -o Dickerson_genion.gro -s Dickerson_gromppem.tpr -p Dickerson_pdb2gmx.top
grompp -f Dickerson_em -c Dickerson_genion.gro -p Dickerson_pdb2gmx.top -o Dickerson_grompp2em.tpr
mdrun -s Dickerson_grompp2em.tpr -o Dickerson_md_em.trr -c Dickerson_md_em.pdb -v
with the following Dickerson_em.mdp file (save this as a file with the name Dickerson_em.mdp in the same directory as Dickerson.pdb):
cpp = /usr/bin/cpp
define = -DFLEXIBLE
integrator = steep
nsteps = 25000
emtol = 10.0
emstep = 0.01
nstcgsteep = 100
coulombtype = PME
rvdw = 1.0
rlist = 1.1
rcoulomb = 1.1
pme_order = 4
ewald_rtol = 1e-5
vdwtype = shift
ns_type = grid
nstlist = 10
Hi, your installation notes worked wonderfully. Could you share a pr.mdp file too?
Regrads.
I'm glad it works somewhere other than my office!
A boilerplate pr.mdp file is provided below. After a proper energy minimization, the PR run should do little more than "prime" the system for the MD simulation (I liken it to priming a cymbal with gentle taps before actually hitting it). Settings are, for lack of a better explanation why, an average of other pr.mdp files I found on the 'net while learning the GROMACS ropes (and there's a considerable bit of commented-out flags).
I assume your cpp is sitting in /usr/bin. The -traditional flag is commented on at Sorin's AMBER/GROMACS page (chemistry.csulb.edu/ffamber).
Ignore the asterisks…
***********************************************************************************
; PRE-PROCESSING CODE
; define = -DFLEXIBLE
define = -DPOSRES
title = AMBER
cpp = /usr/bin/cpp -traditional
; RUN CONTROL
integrator = md
dt = 0.001
nsteps = 10000
comm_mode = Linear ; for periodic systems
nstcomm = 1
; OUTPUT CONTROL
nstxout = 500
nstvout = 500
nstfout = 500
nstlog = 10
nstenergy = 10
nstxtcout = 500
xtc_precision = 1000
; if only the DNA group is to be written to .xtc and .edr, then specify DNA
xtc_grps = DNA SOL
energygrps = DNA SOL
; NEIGHBOR SEARCHING
nstlist = 10
ns_type = grid
pbc = xyz
rlist = 1.2
; ELECTROSTATICS
rcoulomb = 1.2
coulombtype = PME
; VDW
vdwtype = shift
rvdw_switch = 1.0
rvdw = 1.1
; EWALD
fourierspacing = 0.10
pme_order = 6
ewald_rtol = 1e-5
optimize_fft = yes
; TEMPERATURE COUPLING
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = DNA SOL
ref_t = 300 300
; PRESSURE COUPLING
Pcoupl = berendsen
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; VELOCITY GENERATION
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
; BONDS
constraints = all-bonds
; constraint_algorithm =
; unconstrained_start =
; shake_tol =
; lincs_order =
; lincs_iter =
; lincs_warnangle =
; morse =