Atom order of improper dihedral

Hi guys,

I found some different explanations for the atom orders of impropers in LAMMPS data file, and I feel confused.

According to this LAMMPS doc:
https://lammps.sandia.gov/doc/2001/data_format.html

Impropers
1 improper-type atom-1 atom-2 atom-3 atom-4 (atom-2 is central atom)

In another LAMMPS official doc:

https://lammps.sandia.gov/doc/improper_harmonic.html

If the 4 atoms in an improper quadruplet (listed in the data file read by the read_data command) are ordered I,J,K,L then X is the angle between the plane of I,J,K and the plane of J,K,L. Alternatively, you can think of atoms J,K,L as being in a plane, and atom I above the plane, and X as a measure of how far out-of-plane I is with respect to the other 3 atoms.
Note that defining 4 atoms to interact in this way, does not mean that bonds necessarily exist between I-J, J-K, or K-L, as they would in a linear dihedral. Normally, the bonds I-J, I-K, I-L would exist for an improper to be defined between the 4 atoms.

It’s not very clear which atom of the I-J-K-L is the central atom. But if I understand correctly, the last sentence indicates the first atom, I, is the central atom because it bonds to the other atoms.

But from this thread:
https://sourceforge.net/p/lammps/mailman/message/34565747/
Andrew and Arun suggested that atom-3 was the central atom (Let’s talk about OPLS ff and harmonic-style improper here).
The oplsaa in moltemplate also seems to take the atom-3 the central atom, according to the related source codes (oplsaa.lt, nbody_Impropers.py opls_imp.py); however, a case I carried out suggested that the moltemplate may actually take the atom-2 the central atom (could be a bug). Here is the case: I used the original file ethylene.lt given by moltemplate_2019-8-28 to generate a data file. The atoms are numbered as:
3H 5H
\ /
1C = 2C
/
4H 6H
The impropers generated by moltemplate.sh combined with oplsaa.lt turned out as following:
1 1 1 2 5 6
2 1 1 2 6 5
3 1 2 1 3 4
4 1 2 1 4 3
5 1 3 1 2 4
6 1 5 1 2 6

If I understand correctly, here moltemplate takes 2-atom the central atom. By the way, here moltemplate seems to generate redundant impropers though it may not affect the computation in LAMMPS.
I also made some data files for other molecules using moltemplate. For those cases, I feel moltemplate could sometimes take atom-2 the central atom and sometimes take the atom-3, but I’m not quite sure.

I hope Andrew can take a look to this issue. Any comments are appreciated!

Thanks,
Lingnan

three comments on this:

  1. what matters is the description for the individual improper dihedral. failing that, you would have to look at the comments in the source code. if i remember correctly, different force fields and thus different improper styles would follow different conventions.
  2. for as long as the improper interactions are symmetric, the order should not matter
  3. any moltemplate specific comments i have to defer to andrew

axel.

Thanks Axel

As Axel mentioned: what matters is the equation describing the energy
as a function of the 4 atoms. It is up to you to interpret which of
the 4 atoms is the "central" atom.

I found some different explanations for the atom orders of impropers in LAMMPS data file, and I feel confused.

According to this LAMMPS doc:
https://lammps.sandia.gov/doc/2001/data_format.html

Impropers
1 improper-type atom-1 atom-2 atom-3 atom-4 (atom-2 is central atom)

This is not true in general. (Sorry, I realize this is confusing.)
This documentation is for the 2001 version of LAMMPS. (Perhaps at
that time, LAMMPS only supported improper styles which were intended
to be used with the 2nd atom in the center.) But there are many
improper styles these days, and each one uses different equations
describing the force between these 4 atoms.

In another LAMMPS official doc:
https://lammps.sandia.gov/doc/improper_harmonic.html

If the 4 atoms in an improper quadruplet (listed in the data file read by the read_data command) are ordered I,J,K,L then X is the angle between the plane of I,J,K and the plane of J,K,L. Alternatively, you can think of atoms J,K,L as being in a plane, and atom I above the plane, and X as a measure of how far out-of-plane I is with respect to the other 3 atoms.
Note that defining 4 atoms to interact in this way, does not mean that bonds necessarily exist between I-J, J-K, or K-L, as they would in a linear dihedral. Normally, the bonds I-J, I-K, I-L would exist for an improper to be defined between the 4 atoms.

It's not very clear which atom of the I-J-K-L is the central atom. But if I understand correctly, the last sentence indicates the first atom, I, is the central atom because it bonds to the other atoms.

Unlike the 2001 documentation, the current LAMMPS documentation is
intentionally vague (which is wise).
The current LAMMPS documentation for each reveals only the formula for
the energy as a function of the atom coordinates (for that
dihedral_style). You can use that interaction any way you like.
There may not even be a central atom. The 4 atoms do not even need to
be bonded together (although leaving out bonds would probably be a bad
idea...)

But from this thread:
https://sourceforge.net/p/lammps/mailman/message/34565747/
Andrew and Arun suggested that atom-3 was the central atom (Let's talk about OPLS ff and harmonic-style improper here).

Here is another message from that same thread you found.
I discuss the huge variety of improper styles in this post:
https://sourceforge.net/p/lammps/mailman/message/34567561/
(See correction below.)

Different force fields (eg. OPLS, AMBER, CHARMM, GROMOS, COMPASS,
PCFF, CVFF, UFF, MMFF94, ...) use different improper styles.

Nobody agrees which atom is the central atom. Perhaps you are asking
a more specific question question?

"Which atom is the central atom for the OPLSAA force field?"
(Or, perhaps more specifically: Which atom is the central atom for the
OPLSAA force field as implemented in moltemplate?)

In that case, the 3rd atom (ie., atom#2, since they start at 0) is
used as the central atom. For OPLSAA-specific details, see:
https://github.com/jewettaij/moltemplate/blob/master/moltemplate/nbody_alt_symmetry/opls_imp.py
For AMBER/GAFF/GAFF2 in moltemplate, it is the 3rd atom as well:
https://github.com/jewettaij/moltemplate/blob/master/moltemplate/nbody_alt_symmetry/gaff_imp.py
For COMPASS in moltemplate, it is the 2nd atom:
https://github.com/jewettaij/moltemplate/blob/master/moltemplate/nbody_alt_symmetry/cenJsortIKL.py

But this is specific to moltemplate. In moltemplate, force field
parameters are saved in a format where the atoms in each improper
interaction are listed in the same order that they will eventually
appear in the DATA file which moltemplate will generate (and LAMMPS
will read). (See the "Impropers By Type" sections in "oplsaa.lt" or
"gaff.lt".) Other molecule-builder programs will use their own
conventions for atom-order. For example the EMC, AmberTools, and
PSFgen molecule builders could conceivably use an internal file
force-field file format which lists the atoms in a different order
(but I am too lazy to check whether this is true).

Correction: In
(https://sourceforge.net/p/lammps/mailman/message/34567561) I said:
"If you are using moltemplate, then by default, it assumes that the
1st atom is the "central" atom, although this can be overridden."

...I should have added that most force-fields in moltemplate will
override this choice. (See above.)

I hope this helps

Andrew

Andrew and Axel,

Thank you very much for your comments and detailed explanation! Appreciate it.

Best,
Lingnan