Conversion of LJ parameters for ions

[rkingsbury]

Hi @htz1992213 , I'm trying to understand how to convert literature Lennard Jones parameters into LAMMPS format, and I'm concerned that the ones for ions stored in /data/ions may not have been converted correctly.

For example, consider Na+ from the Joung-Cheatham paper here

Here is what we have in mdgo for Na+ for the SPCE model:
https://github.com/htz1992213/mdgo/blob/c843a9643a67992cf5069e7be668b194e86ef889/mdgo/data/ion/joung_cheatham/spce/Na%2B.lmp#L16-L18

Here is the relevant excerpt from the paper:

The LAMMPS docs say that

Note that sigma is defined in the LJ formula as the zero-crossing distance for the potential, not as the energy minimum at 2^(1/6) sigma

For whatever reason, it's conventional to report sigma as r_min/2 in many papers I've seen. The potential should equal zero when sigma = r. So if I understand correctly, to convert rmin/2 to the sigma that LAMMPS wants, you should do

(r_min / 2 * 2) / 2^(1/6)

If I perform this conversion on the Na+ value from the paper, I get 2.159, vs. 2.238 in mdgo. Am I missing something? Also, would there be benefit to storing more decimal places for the epsilon parameter?

This is incredibly tricky and I'm really surprised that this subtle point about how LAMMPS defines LJ parameters is not documented more clearly (or that the research community doesn't just report sigma instead of r_min/2). Please let me know what you think, because as I add new parameters I want to make sure I'm converting them correctly.

[rkingsbury]

@kdfong this sounds like it's related to one of the subtle differences between codes you were telling me about. Do you think I'm understanding this correctly?

[rkingsbury]

As an additional point of reference, I found this example file from moltemplate that seems to confirm my understanding and my number for Na+:

# Note: Monovalent ion parameters for Ewald and SPC/E water
#       are from Joung & Cheatham JPCB (2008)   (Table 5)
#       (The widths of the ions, expressed in terms of Rmin/2, are:
#        1.212 and 2.711, for Na+ and Cl-, respectively)
# Note: They use  U(r) =   epsilon*((sigma/r)^12 - 2*(sigma/r)^6)
#            not  U(r) = 4*epsilon*((sigma/r)^12  -  (sigma/r)^6)
#       ...but this should not effect the epsilon value.
#       (LAMMPS uses later convention, in which case Rmin/2 = sigma/2^(5/6))

[htz1992213]

Hey @rkingsbury, I think that is a very interesting topic!!! Based on my understanding of the "Simulation conditions" part of the JC paper, I derived the conversion relation as follows:

So here, i should be the ion, and j is the water. They used 1.7699 Å as the Rmin/2 value of TIP4P water to convert from the JJ to JC. But it seems that I mistakenly used that 1.7699 Å value for converting ion parameters of other water models. So I believe the correct way here should be to use different Rmin/2 or σ values of the corresponding water model to convert ion parameters.
I think the simplified conversion using just Rmin/2 = σ/2^(5/6) as specified in moltempalte is not the most rigorous way. It could possibly cause a non-trivial difference if the Rmin/2 or σ value of a certain ion is too far way from that of water O. The importance of choosing a correct mixing rule is also highlighted in the JC paper:

Because the combining rules are different between different parametrizations and codes, care needs to be taken to use the correct mixing rule. The incorrect choice can lead to drastically different results.

[htz1992213]

@rkingsbury Have opened PR #14, all values now conform to the above JJ form (opls, geometric mixing). I believe moltemplate is more specialized in Amber FF, and amber utilizes arithmetic mixing, so that is why they used the simplified conversion in the example. What do you think would be a more suitable mixing rule here as the default? I checked reference papers, Amber also uses geometric mixing for ε, so I think the above form should be more suitable than the simplified form, right?

[rkingsbury]

Thanks @htz1992213 , this sounds more complicated than I realized and it sounds like geometric mixing is the way to go. . Let me read the paper you show above and I'll reply here again. My main, immediate interest in this is figuring out how to properly convert values from this paper so that they can be used with TIP4P-FB water. I'd also like to see if I can find an example somewhere that confirms this approach.

[rkingsbury]

OK, I've done a little more reading about this and I think we might be talking about slightly different things (or at least, there is a bigger philosophical question to address here).

If I understand correctly, the ion values presently stored in mdgo are intended to be pre-mixed with the respective water parameters, and if that's the case, then I agree with your logic on how to convert them. The mixing rules depend on the way the ion parameters were fit.

However, I'm not sure we should be storing pre-mixed parameters in mdgo. Because whether you need single atom parameters or pre-mixed parameters depends on how you set up the simulation. Wouldn't it be clearer and more flexible to just store the single atom parameters and rely either on 1) the user, 2) code we will write later, or 3) LAMMPS default settings to mix the parameters?

For example, in my runs, I have not been explicitly setting the pair_style for every pair type. Rather, I just provide the single atom parameters in my data file like this:

Pair Coeffs

1   0.000000  0.0000
2   0.749280  3.1655
3   0.170000  2.9600
4   0.250000  3.5500
5   0.170000  2.9600
6   0.170000  2.9600
7   0.070000  3.5500
8   0.070000  3.5500
9   0.070000  3.5500
10  0.070000  3.5500
11  0.070000  3.5500
12  0.070000  3.5500
13  0.030000  2.4200
14  0.030000  2.4200
15  0.030000  2.4200
16  0.030000  2.4200
17  0.030000  2.4200
18  0.025454  2.5836

And then rely on LAMMPS to mix them. It appears that LAMMPS uses geometric mixing by default (see pair_modify command), and from the documentation of the TIP4P model in LAMMPS it says:

For atom type pairs I,J and I != J, the epsilon and sigma coefficients and cutoff distance for all of the lj/cut pair styles can be mixed. The default mix value is geometric. See the “pair_modify” command for details.

So if we wanted to only store single atom parameters and externalize the mixing somewhere else, I think the conversion formula I posted at the top of this issue is correct. In addition, if we only store pre-mixed ion parameters between the ion and water O, it's not clear to me how the user or LAMMPS could easily calculate the correct pair parameters for ion-ion interactions (e.g., in a binary, ternary, or higher system with many ions).

One final complication is that some of the parameter sets use a combination of arithmetic and geometric mixing rules, which appears equivalent to the non-default arithmetic style supported by LAMMPS pair_modify. From the JC paper, just above the text you quoted:

The AMBER,(80, 81) CHARMM,(82) and NAMD(83) programs, by default, apply the Lorentz−Berthelot(84-87) combining rules which use the arithmetic mean for the combined Rmin and the geometric mean for the combined ε. Other programs, such as BOSS(88) and GROMOS,(89) use geometric mean combining rules.

The LB rules are also used in the newer parameter set that I want to add:

Whatever we decide to do, we should definitely document in the respective data files whether the parameters are pre-mixed or single atom, because if someone sets up a simulation like I have and does not realize the parameters are pre-mixed, LAMMPS is going to mix the parameters anyway and lead to unexpected results.

Given all of this, I would advocate only storing single atom values in mdgo and implementing mixing elsewhere. This is surprisingly complex to figure out though and I may be missing something. What do you think?

[rkingsbury]

Per discussion earlier today, I will work on updating the ion parameters so that what is stored in mdgo are the single atom parameters (i.e., sigma_i not sigma_ij). We'll then add some code to the input construction parts of mdgo that implements mixing rules for the parameters when appropriate based on the specific context.