Dihedral coefficients in loplsaa.lt of moltempalte may be problematic

LAMMPS LAMMPS Mailing List Mirror
1

Hi guys,

I compared the dihedral coefficients in the loplsaa.lt of the latest moltempalte (2019-8-28) with those of the original papers (Sui, Pluhackova, Böckmann, J.Chem.Theory.Comp (2012), 8(4), 1459 and Pluhackova,…,Böckmann, J.Phys.Chem.B (2015), 119(49), 15287), and I found some problems.

The dihedral style that Böckmann’s group use is Ryckaert-Bellemans (RB) function:
V = C0 + C1 cos(ϕ) + C2 cos(ϕ)^2 + C3 cos(ϕ)^3

In the moltemplate version (or more accurately, the Tinker version), the dihedral style is a Fourier function, which is called “opls” in LAMMPS:
V = 1/2 * [F1(1+cos(ϕ))+F2(1−cos(2ϕ))+F3(1+cos(3ϕ))+F4(1−cos(4ϕ))]

There is conversion relation to translate a Fourier function to the RB function:
C0 = F2 + (F1+F3)/2.0
C1 = 0.5 * (-F1+3F3)
C2 = -1.0F2 + 4*F4
C3 = -2 * F3

see Eq. (33) of this Gromacs doc:
http://manual.gromacs.org/current/reference-manual/functions/bonded-interactions.html

However, using the above relation to translate a RB function to a Fourier function could sometimes lead to mistake. If you use the above relation to convert the parameters for RB (C0~C3) to “opls” parameters (F1~F4), you’ll find you don’t really need C0, C1~C3 are sufficient to yield F1~F3, and F4 = 0. This is how Tinker convert the parameters to “opls” form. Then when you use these resulted F1~F4 to calculate back to C0’~C3’, you’ll find the C0’ don’t equal to the original C0. (I’ve tested the parameters from loplsaa.lt, trying to convert them back to the orignal RB form, and it turned out that I obtained correct C1~C3 but failed to get correct C0). Böckmann’s group also performed a conversion in their first lopls paper (Siu et al. 2012, see Table 2). If you use the above relation to convert their Fourier dihedral parameters (middle part of that table), you can’t get their correct C0 either.

Actually, given a RB function with C0,C1,C2,C3, the corresonding Fourier function should be:
V = 1/2 * [F0(1+cos(0)+F1(1+cos(ϕ))+F2(1−cos(2ϕ))+F3(1+cos(3ϕ))+F4(1−cos(4ϕ)]
= F0 + 1/2 * [F1(1+cos(ϕ))+F2(1−cos(2ϕ))+F3(1+cos(3ϕ))+F4(1−cos(4ϕ)]

I guess this may be what Böckmann’s group actually used for conversion. The constant term, F0, is essential to account for the C0. Using this form to translate their Fourier parameters, I successfully reproduced the C0~C4.

Therefore, I think we can’t use “opls” dihedral-style in LAMMPS to convert the RB dihedral function (like Tinker did). Instead, an option is to use “multi/harmonic” dihedral-style in LAMMPS which is the same form as the RB function.

I hope anyone can correct me if I’m wrong. Any comments or ideas are welcome!

Best,
Lingnan

2

I'm forwarding (BCC) this to Sebastian Echeverri, who sent me the
LOPLSAA parameters in TINKER format, in case he has any comments.
(Sebastion, don't feel obligated to reply to this message unless you
want to.)

Hi Lingnan

    Thank you for looking into this! I am really grateful when users
take the time to report problems back to us when there are issues.
There have been issues with moltemplate force fields in the past, and
user feedback is critical. Without people like you, moltemplate would
be unusable.
    I agree that, as a diagnostic, users will want to verify that the
energies calculated by LAMMPS match the energies calculated by other
programs. (...and they are wise to want to do this) But there are
some serious drawbacks to changing the dihedral style to enable us to
independently modify the C0 offset parameter for the LOPLSAA force
field. So I think I will postpone working on this suggestion for now.
Below I go into the reasons why:

1) (As you know) The C0 parameter does not effect the behavior of
the molecule. It is just an energy offset which does not effect the
forces that act on the molecule. It would not surprise me if we got
the C0 dihedral parameters wrong. (I'm mentioning this for the
benefit of anyone who might stumble upon this conversation via
google.)

2) GPU support: As you also realize, there is no way to
independently choose the C0 offset when using dihedral_style opls
(which is what we currently use). As you pointed out, we would have
to switch to dihedral_style multi/harmonic or (preferably!)
dihedral_style fourier, or something else. That would require either
converting the rest of the OPLS dihedral_coeffs in the original
"oplsaa.lt" file to use one of these more general dihedral styles, or
using some kind of hybrid dihedral style (eg "dihedral_style hybrid
opls fourier": opls for the portions of the molecule that use ordinary
OPLSAA parameters, and fourier or multi/harmonic for the LOPLSAA
parameters.) Unfortunately GPU acceleration has not yet been
implemented for these dihedral styles.

(Incidentally, if anybody stumbles into this conversation and wants
to implement GPU support for additional dihedral styles, please give
priority to dihedral styles that are really general like
"dihedral_style fourier" style. dihedral styles like "multi/harmonic"
have comparatively limited niche uses.)

3) In the future, I will stop using TINKER as a source for my
force-field parameters.
     Hopefully sometime next year, I hope to get a
CHARMM-to-MOLTEMPLATE converter implemented. Then we will obtain all
of the OPLSAA parameters directly from the original files distributed
by the Jorgensen group which are in CHARMM (.INP, .PRM) format. That
way, the atom type names in moltemplate will match their types in the
original CHARMM files, and users will be able to take advantage of
PSFGEN and VMD/AUTOPSF to handle the atom-typing and partial charge
assignment. (Generally this will make it easier for moltemplate and
other tools like VMD to work well together.) When we do that, we can
start to worry about getting the C0 offsets right.
    As you realized, when converting this to moltemplate format, I
used the tinkerparm2lt.py script (the same script that I normally use
to convert TINKER files). So even if there wasn't a problem with GPU
compatibility, I'm a little reluctant to spend much more time on
anything that relates to maintaining the TINKER->MOLTEMPLATE
converter.

--- offer ---
   If you send me a version of "loplsaa.lt" using dihedral_style
multi/harmonic (or dihdral_style hybrid opls multi/harmonic), I will
mention it in the README files for the LOPLSAA examples, and make it
available for download. That way we can see if molecules using
LOPLSAA (and OPLSAA) have the correct energy, and users can feel
assured that they are using LOPLSAA correctly. (I am not certain what
we will discover if we try that, but I'd be interested to find out.)

   Thanks again for reporting the issue. I hope you don't feel
slighted. Even if I don't implement your suggestions, I really
appreciate being made aware of the issue, because I will likely get
asked this question again (and I'll keep my eyes open in case new GPU
compatible dihedral styles become available in the future).

Andrew

   P.S. If the LOPLSAA parameters are available somewhere in CHARMM
format, please let me know where to find them.

3

Hi Andrew,

Thanks for your detailed reply!

It’s very wise to stop using TINKER as a source and start developing the tool that converts parameters directly from the original files distributed by force field authors. I look forward to that converter.

Unfortunately, today I found another problem in that loplsaa.lt file, and unlike the dihedral_style, this has significant impact on the simulation.

In pair_coeff section, the sigma values of some newly added atoms are seemed to be 1/10 of what they are supposed to be. For example, the sigma of atom:81LL_b013_a013_d013LL_i013 (carbon in “Alkane -CH2-”) here is 0.35 A, but it should 3.5 A or within this order of magnitude. These parameters, though, are in accordance with those in the LOPLS paper (see Table 1, Pluhackova et al., J.Phys.Chem.B (2015)). However, I thought there must be a typo for the unit of sigma in that table (I haven’t confirmed with the authors), because I also checked the force field files the authors released (in Gromacs format, see https://www.biomemphys.nat.fau.de/downloads/), the unit should be nm rather than angstrom.

Because of these seemingly wrong sigma values, I kept receiving error “Out of range atoms” when I was trying to equilibrate a system using a data file generated by that force field file. I think that error must be caused by some unrealistic forces. When I change those sigma values, the errors disappeared and the results are normal.

Right now I’m working on a ester molecular system. I’ve done some changes to the files loplsaa.lt and oplsaa.lt. I’ll send you those files when I complete more tests and make sure there are no more problems. This won’t take long. Very glad I can do some contribution to the community.

Thanks again for your effort and dedication on the moltemplate.

Best,
Lingnan

4

Unfortunately, today I found another problem in that loplsaa.lt file, and unlike the dihedral_style, this has significant impact on the simulation.

In pair_coeff section, the sigma values of some newly added atoms are seemed to be 1/10 of what they are supposed to be. For example, the sigma of atom:81LL_b013_a013_d013LL_i013 (carbon in “Alkane -CH2-”) here is 0.35 A, but it should 3.5 A or within this order of magnitude.

Well, that is just plain embarrassing.
Thank you very much for catching this!
The “loplsaa.lt” file was generated by running “tinkerparm2lt.py” on a file named “loplsaa.prm”, which was in tinker format. (That file also has these mistakes. No excuses though. I should have looked at the final “loplsaa.lt” file carefully, before posting it.)

These parameters, though, are in accordance with those in the LOPLS paper (see Table 1, Pluhackova et al., J.Phys.Chem.B (2015)).

Yeah. After a quick glance at the file, the other pair_coeff commands seem okay. It looks like the epsilon values in the “pair_coeff” commands are correct (although I did not check all of them). (At least we remembered to convert from kJ/mole to kCal/mole.)

However, I thought there must be a typo for the unit of sigma in that table (I haven’t confirmed with the authors), because I also checked the force field files the authors released (in Gromacs format, see https://www.biomemphys.nat.fau.de/downloads/), the unit should be nm rather than angstrom.

Because of these seemingly wrong sigma values, I kept receiving error “Out of range atoms” when I was trying to equilibrate a system using a data file generated by that force field file. I think that error must be caused by some unrealistic forces. When I change those sigma values, the errors disappeared and the results are normal.

I’m sorry for making you suffer so.

Right now I’m working on a ester molecular system. I’ve done some changes to the files loplsaa.lt and oplsaa.lt. I’ll send you those files when I complete more tests and make sure there are no more problems. This won’t take long. Very glad I can do some contribution to the community.

Thank you.
Actually, I just made the modifications you suggested to the “loplsaa.lt” and “loplsaa_ext.prm” files and uploaded them to github.
If you would like credit for contributing, and if you use github, then please make a small change to the “loplsaa.lt” file (for example, add a blank line at the end of the file), and issue a pull-request. I will approve the pull request, so that you will appear on the list of contributors. (If you have time, please do this, because it helps me keep track of people who improved moltemplate.)

And I am still curious to hear the outcome of your tests.

If you spot any more errors, let me know. (There might be. You can probably tell that we were not very careful.)

Thanks again for your effort and dedication on the moltemplate.

No, thank you. Without people reporting bugs, moltemplate would suck. A lot.

Great job!

Andrew

5

Hi Andrew, Lingnan

I am glad to see that the implementation of the lopls is being used and that it is being critically revised. As Andrew mentioned, it is vital to have the input of the users to fix all the small bugs that can be in the software.

Thank you both for fixing the “loplsaa.lt” and “loplsaa_ext.prm” files

Cheers,
Sebastián