Implementing OPLS improper for alcene molecules

Dear all,

I am confused about how to implement OPLS improper dihedrals for alcene molecules.

I am using (L-)OPLS force field and I extracted the parameters from gromacs files. In these files it is indicated that impropers should be implemented with an energy of the form:
E = A [1-cos(2x)] with a stiffness A of 14.995 kcal/mol and x the dihedral angle.
This energy form can be implemented in Lammps using the cvff improper style.
But it is not clear from Gromacs files which improper angle x should refer to.

The documentation for cvff improper style indicates that, when specifying an improper dihedrals by a sequence of four atomes I-J-K-L, the improper angle x corresponds to the angle between I-J-K and J-K-L planes and I being the central atom.

Taking the example of ethene molecule

1H 2H
\ /
3C = 4C
/ \
5H 6H

I wrote:
- in the data file read with read_data command:

Impropers

1 3 1 4 5
1 4 3 2 6

- in main input file:

improper_style cvff
improper_coeff 1 14.995 -1 2
improper_coeff 2 14.995 -1 2

Is the stiffness value 14.995 Kcal/mol correct for this convention of the improper angle ?
I guess choosing an incorrect improper angle would lead to incorrect out-of-plane fluctuations for the atoms...

Thanks all,

Ludovic

Ordering of atoms for impropers is tricky.
I’m CCing Andrew Jewett who as I recall has looked into
this in detail. Maybe he will have ideas for you.

Steve

I am using (L-)OPLS force field and I extracted the parameters from
gromacs files. In these files it is indicated that impropers should be
implemented with an energy of the form:
E = A [1-cos(2x)] with a stiffness A of 14.995 kcal/mol and x the
dihedral angle.
This energy form can be implemented in Lammps using the cvff improper
style.

The implementation of OPLSAA in moltemplate is similar, but it uses improper style harmonic

improper_style harmonic
improper_coeff 1 15.0 180.0

This is compatible with your implementation because the 2x in your formula cancels with the (1/2) when you taylor-expand the cos() function.

But it is not clear from Gromacs files which improper angle x should
refer to.

https://github.com/jewettaij/moltemplate/tree/master/examples/all_atom/force_field_OPLSAA/ethylene+benzene

The documentation for cvff improper style indicates that, when
specifying an improper dihedrals by a sequence of four atomes I-J-K-L,
the improper angle x corresponds to the angle between I-J-K and J-K-L

Yes
(It varies depending on the improper style, but what you said is true for these two improper styles, at least:
https://lammps.sandia.gov/doc/improper_harmonic.html
https://lammps.sandia.gov/doc/improper_cvff.html)

I being the central atom.

This depends on the convention adopted by the people who implemented OPLSAA in GROMACS or MOLTEMPLATE. If my understanding is correct, for OPLSAA, it could be either the 2nd or 3rd atom (J or K), but not the first. In the moltemplate version of OPLSAA, the central atom is the 3rd atom (K). (Details here.)

The moltemplate implementation also creates 3 different improper interactions for each carbon atom (6 total for ethylene), because there are 3 different non-equivalent angles, depending on the order of the 3 surrounding atoms. I made this choice because it would seem arbitrary to pick one of these angles instead of the other 2, so moltemplate generates all 3. I thought this was the safest thing to do. If I generate 3 impropers per carbon atom, then perhaps I should divide the strength of each improper interaction by 3 to compensate. However this may not be correct. I may be making the improper interactions more stiff than they should be.

For example, I noticed that data file only contains 2 impropers (one for each carbon atom). AmberTools apparently also only creates 1 improper per carbon atom in ethylene (and it seems to choose which atom order to use based on its own arbitrary internal logic which is not documented anywhere). I myself am not certain whether I should be generating 1 improper interaction or 3 (per carbon atom). I had someone check this a couple years ago to verify I am generating impropers correctly, but perhaps we got it wrong. I’d be happy to have someone else verify this again. The only way to be sure is to compare the files created for GROMACS or the data files created by MOLTEMPLATE with the files created by the official software (BOSS or perhaps NAMD).

I don’t know if this helps you, but in any case, I have attached a small LAMMPS data file containing an ethylene molecule generated by moltemplate. I also included the (automatically generated) .in files that contain the force-field and charge information for that molecule. (The charge is not included in the data file, but in a separate file “system.in.charges”) Feel free to let me know if we messed up our implementation.

Impropers are always an incredible headache, because each force field implements them in completely different (somewhat arbitrary) ways, and some implementations are more symmetric than others.

Andrew

system.data (770 Bytes)

system.in.init (255 Bytes)

system.in.settings (279 Bytes)

system.in.charges (52 Bytes)