Using the AMBER forcefield with NAMD
In our lab the standard protocol for running MD simulations entails the NAMD with the CHARMM forcefield. So, every time we wish to push our lack with some other widespread and well-known forcefield, is mainly for comparison reasons.
Using the OPLS-AA is quite straightforward and described here Using the OPLS-AA forcefield
Using the AMBER forcefield is more tricky especially if you want to use it with the NAMD engine. AMBER11 is not free, but the AMBERTOOLS is. Assuming you have download it and correctly set it up in your local machine you can use the Xleap tool to set up your simulation system and prepare the parameter and topology files. If you feel lucky…, you can work directly in Norma where the AMBERTOOLs package is already installed
For those who like the NMG-famous “RTFM-way” here is the manual AmberTools
For the rest that are in favor of the DO-ER way here are the steps to follow:
- Run XLeap
/usr/local/amber11/exe/xleap -s -f /usr/local/amber11/dat/leap/cmd/leaprc.ff99SB
This will execute the start-up script for the ff99SB forcefield. The -s flag is to ignore any defined defaults which might override our selection
The XLeap has an known bug that needs the NumLock to be turned-off! This does not apply to the terminal version (tleap).
- Load the PDB file
DTRW = loadPDB "DTRW.pdb"
This will load the local pdb file. You can also built your extended peptide, but the result is that, there are so many steric clashes with the side-chains (especially if you have proline residues) that you need to perform a short minimization of the structure before continuing. To avoid all that, prepare you system with the old-good way and take the pdb file. You can use the pdb output from ribosome but the Leap cannot understand some of the protons (H) in order to write the parameter file. I use the “aligned.pdb”, which is the output from the moleman.
edit DTRW
to observe the structure. Here is the output of the program so far:
Welcome to LEaP!
Sourcing: /usr/local/amber11/dat/leap/cmd/leaprc.ff99SB
Log file: ./leap.log
Loading parameters: /usr/local/amber11/dat/leap/parm/parm99.dat
Reading title:
PARM99 for DNA,RNA,AA, organic molecules, TIP3P wat. Polariz.& LP incl.02/04/99
Loading parameters: /usr/local/amber11/dat/leap/parm/frcmod.ff99SB
Reading force field modification type file (frcmod)
Reading title:
Modification/update of parm99.dat (Hornak & Simmerling)
Loading library: /usr/local/amber11/dat/leap/lib/all_nucleic94.lib
Loading library: /usr/local/amber11/dat/leap/lib/all_amino94.lib
Loading library: /usr/local/amber11/dat/leap/lib/all_aminoct94.lib
Loading library: /usr/local/amber11/dat/leap/lib/all_aminont94.lib
Loading library: /usr/local/amber11/dat/leap/lib/ions94.lib
Loading library: /usr/local/amber11/dat/leap/lib/solvents.lib
> DTRW = loadPDB "DTRW.pdb"
Loading PDB file: ./DTRW.pdb
Added missing heavy atom: .R<CTRP 4>.A<OXT 25>
total atoms in file: 40
Leap added 37 missing atoms according to residue templates:
1 Heavy
36 H / lone pairs
>
If you see something in the spirit of “FATAL: Atom .R<ILE 26>.A<CD 21> does not have a type.” just look carefully in you PDB file and correct the “false” atom types to the atom types that leap reads [ “all_amino94.lib” file ].
If you wish to add capping ends (like ACE and/or NME) in a pre-existing PDB file, just edit your PDB and change the atom name in one of the atoms that belong to the free end and leap will take care of the rest (assuming you know which atoms need to be changed)
- Since I want to make a system comparable to the one I have already using the VMD, I will ionize and solvate my peptide using approximately the same number of atoms, i.e. 1 Na+, 1Cl- and 269 water residues (807 atoms):
> addions DTRW Na+ 1 Cl- 1 Adding 2 counter ions to "DTRW" using 1A grid Grid extends from solute vdw + 1.87 to 7.78 Resolution: 1.00 Angstrom. grid build: 0 sec (no solvent present) Calculating grid charges charges: 0 sec Placed Na+ in DTRW at (8.73, -1.41, -1.22). Placed Cl- in DTRW at (6.73, -5.41, -2.22). Done adding ions.
- and I will use a truncated octahedral water box
> solvateoct DTRW TIP3PBOX 2.6
Scaling up box by a factor of 1.257315 to meet diagonal cut criterion
Solute vdw bounding box: 21.781 12.827 8.739
Total bounding box for atom centers: 28.319 28.319 28.319
(box expansion for 'iso' is 171.1%)
Solvent unit box: 18.774 18.774 18.774
Volume: 12601.657 A^3 (oct)
Total mass 5499.378 amu, Density 0.725 g/cc
Added 270 residues.
The above line will put an 2.6Å buffer of TIP3P water around the peptide in each direction. I want approximately the same number of atoms so the 2.6 was selected because it adds 270 waters (compare with 269 of the VMD script).
- Save a PDB file of the ionized and solvated system (similar to the ionized.pdb)
> savePDB DTRW DTRW.amber.pdb Writing pdb file: DTRW.amber.pdb Converting N-terminal residue name to PDB format: NASP -> ASP Converting C-terminal residue name to PDB format: CTRP -> TRP
- built the prmtop and inpcrd files. The prmtop, contains a description of the molecular topology and the necessary force field parameters, and the inpcrd, contains a description of the atom coordinates and optionally velocities and current periodic box dimensions
> saveamberparm DTRW dtrw.prmtop dtrw.inpcrd Checking Unit. Building topology. Building atom parameters. Building bond parameters. Building angle parameters. Building proper torsion parameters. Building improper torsion parameters. total 19 improper torsions applied Building H-Bond parameters. Not Marking per-residue atom chain types. Marking per-residue atom chain types. (Residues lacking connect0/connect1 - these don't have chain types marked: res total affected CTRP 1 NASP 1 WAT 270 ) (no restraints)
Before continuing, observe the two systems prepared with AMBER's XLeap, and VMD:
As you see, the orthogonal frame is not the same, so we need to define different unit cell translations for the amber-generated truncated octahedron in the NAMD script. The simulation box of truncated octahedron generated by leap gives a triclinic box with lattice vectors:
| d | 0 | 0 |
| -d/3 | (2*sqrt2*d)/3 | 0 |
| -d/3 | -(sqrt2*d)/3 | -(sqrt6*d)/3 |
Since AMBER treats all periodic cells as triclinic, unit cell directions are 109° apart and have distance of half the diagonal across the cube: (cube edge dist) * (sqrt3 / 2),
where the cube edge dist is the total bounding box in the AMBER output (in our case 28.319). This gives a diagonal of 24.52Å. Leap adds on a little to avoid bad VdW contacts across the periodic boundaries, thus the value 25.39Å in the last line of the inpcrd file.
We use again the little trick where we tell NAMD that the cell edge is somewhat smaller to avoid significant volume reduction and large negative pressures at the beginning of the simulation when turning the barostat on. So, we consider d = 24.00 (and not 25.39, nor 24.52 ) in the previous matrix.
The new unit cell translations is not the only difference. There are still some things you need to consider with the new force-field: AMBER forcefield parameters. Here is how my heating-up and equilibration scripts look like:
The protonation state specification in AMBER is different. The correct residue names are HIP, HID, and HIE.