Van der Waals Parameters for Cations in AMBER Format

The AMBER force field contains parameters for a variety of ions; however, if we are interested in studying a system containing a species not included in the standard parameter set we must append the force field with suitable quantities. The "state of the art" for monovalent and divalent cations is attributable to Aqvist [J. Aqvist, J. Phys. Chem. 94, 8021 (1990)]. We are charged with determining the parameters rij* (the internuclear separation of the ij pair at the potential minimum) and e, the potential well depth for the ion at this minimum value.

Given is the Lennard-Jones potential:

[1] U(r) = e (r*/ r)12 - 2e (r*/ r)6

which is commonly rewritten for the ij pair of ions as:

[2] U(rij) = (AiAj / rij12) - (BiBj / rij6)

where the A and B parameters are given in Aqvist. To find rij*, we take the derivative of [2] with respect to rij:

[3] dU(rij) / drij = -12AiAjrij-13 + 6BiBjrij-7

and set the righthand side of [3] to zero:

[4] -12AiAjrij-13 + 6BiBjrij-7 = 0

then solve for rij at this minimum which is what we call rij*:

[5] rij* = (2AiAj / BiBj)1/6

Using [5] in concert with Aqvist's paper, we can proceed to determine rij* for Mg2+, Ca2+, Sr2+, and Ba2+ (the monovalent species are already in the AMBER parameter set) where species i is one of the cations and species j is the oxygen in water. As an example, we will calculate rij* for Ca2+.

From Aqvist, we have ACa2+ = 264.1 and BCa2+ = 18.82. The water model used in AMBER is TIP3P [W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983)]; the cation interacts with the oxygen in water, thus we need the A and B parameters for the oxygen in TIP3P water (AO = 762.89 and BO = 24.39).

[6] rO-Ca2+* = (2 × 762.89 × 264.1 / 24.39 × 18.82)1/6 = 3.09437 Å

We are almost finished. In [6], rO-Ca2+* is the sum of rO* and rCa2+*:

[7] rO-Ca2+* = rO* + rCa2+*

Rearranging [7], and using rO* = 1.768 Å from Aqvist (note that this number includes the hydrogens implicitly!):

[8] rCa2+* = rO-Ca2+* - rO* = 3.09437 - 1.768 = 1.3264 Å

For e, we just compare [1] to [2]: It is clear that

[9] AiAj = eij (rij*)12

and

[10] BiBj = 2 e ij (rij*)6

Since i = j for each species, we will drop the subscripts in the following discussion. If we square [10]:

[11] B4 = 4 e2 (r*)12

we can solve [11] for (r*)12:

[12] (r*)12 = B4 / 4 e2

Rearrange [9] to

[13] (r*)12 = A2 / e

We can now equate [12] to [13] and solve for e:

[14] e = B4 / 4 A2

Returning to the Ca2+ example, recall from Aqvist we have ACa2+ = 264.1 and BCa2+ = 18.82. Plugging these quantities into [14]:

[15] e (Ca2+) = 18.824 / 4 × 264.12 = 0.44966 kcal/mol

Finally, we are ready to make our file that includes these data. All we need to do is name the file - let's choose frcmod_Ca.in - and set it up as follows (there are a total of 6 lines in the file, including the blank line):

\# First line; parameters for Ca++.
MASS
CA 40.08

NONB
CA 1.3264 0.44966 (adjusted, from Aqvist)

That's all. In tLEaP, we would do the following:

> loadamberparams frcmod_Ca.in
Loading parameters: ./frcmod_Ca.in
Reading force field mod type file (frcmod)
>

and so on. Questions and comments can be directed to Todd J. Minehardt (tjm@princeton.edu).