Multisite Divalent Cation Models
This page describes how to implement out multisite cation models as described in:
A. Saxena and D. Sept, Multisite ion models that
improve coordination and free energy calculations in molecular dynamics simulations, J. Chem. Theor. Comp.
9(8) 3538-42, 2013.
Note that these instructions describe the implementation in Amber, but
other implementations have been performed successfully and are discussed below.
Also note that these procedures and files differ slightly from the
supplemental materials published with the manuscript, however these should be
simpler and faster.
- First, replace the spherical ion with the multisite ion structure from either the calcium or magensium pdb files.
- Use xleap to generate *parm (topology) and *crd (structure) files for the
system. Addition of water and/or counterion can also be done at this step.
Also, make sure to load the *lib files (cal.lib and mag.lib) and *frcmod files (frcmod.mca and frcmod_mg.mmg). This must
be done so that the new atom types are recognized. N.B. that the frcmod and cal.lib files
are one of the changes mentioned above.
- As input, Amber uses the r and ε notation, but since our dummy
atoms have only a repulsive part of the potential and no attractive part, there
is no way to make this work mathematically. However, when xleap generates the
topology file, it converts this back to the A/B form that we discuss in the
paper. The parameters are correct for the central atom (XZ) and the A paramters
are correct for the dummy atoms, however the B values are non-zero and need to be
To do this, we provide a simple python script (zerobvalues.py). This
script will go through the topology file generated above, locate the dummay
atoms, determine the dummy atom type, and correct the LENNARD_JONES_B
Example of Calcium in Water
>tleap -s -f leaprc.ff99SB
> loadoff bondmiss_cal.lib
Loading library: ./bondmiss_cal.lib
> loadamberparams frcmod_ca.mca
Loading parameters: ./frcmod_ca.mca
Reading force field modification type file (frcmod)
\# parameters for multisite Ca++ from Saxena and Sept JCTC (2013)
> mol = loadpdb calcium_ms.pdb
Loading PDB file: ./calcium_ms.pdb
total atoms in file: 8
> solvatebox mol TIP3PBOX 10.0
Solute vdw bounding box: 2.712 2.628 2.800
Total bounding box for atom centers: 22.712 22.628 22.800
Solvent unit box: 18.774 18.774 18.774
Total vdw box size: 25.763 25.372 25.855 angstroms.
Volume: 16900.198 A^3
Total mass 6057.424 amu, Density 0.595 g/cc
Added 334 residues.
> addions2 mol Cl- 2
Adding 2 counter ions to "mol" using 1A grid
Grid extends from solute vdw + 2.47 to 8.24
Resolution: 1.00 Angstrom.
grid build: 0 sec
Calculating grid charges
charges: 0 sec
Placed Cl- in mol at (12.65, 7.32, -1.54).
Placed Cl- in mol at (-13.35, -2.68, 4.46).
Done adding ions.
> saveamberparm mol Ca_Model_bond.top Ca_Model_bond.crd
This is a simple script designed to zero out the Lennard-Jones
B values for dummy atoms in our multisite ion model.
Please input topology file: Ca_Model_bond.top
Modified parameter file will be Ca_Model_bond.top.mod
Topology file has 1012 atoms and 5 atom types
D1 atom at position 7 in atom list
D1 atom has type 2
Making substitution in LENNARD_JONES_BCOEF at positions: [2, 3, 5, 8, 12]
Implementation in Gromacs
This is under development.