Unlike bond-angles and dihedral-angles, there is no universal way to guess impropers.
Improper angle generation is typically only used to enforce planarity. For this reason, one would think they would be very simple. However due to unfortunate decisions made in the distant past, improper angle generation is highly force-field and software dependent. Topotools (if I’m not mistaken), follows the CHARMm, which are similar to conventions used by AMBER, and CVFF. On the other hand, the moltemplate molecule builder supports the the OPLSAA, AMBER, and CVFF conventions, but can be customized to use other rules (explained below. Disclaimer: I wrote moltemplate.)
If you are trying to design your own force-field, then you have to be careful to clarify: when you want a new improper interaction to be generated, and the order you want the atoms to appear in the list. (In this particular case, moltemplate may be useful.) However, if you are trying to use one of the existing force-fields, then great care must be taken to make sure you are using the correct conventions (which are often poorly documented by the people who originally created them). In that case, definitely take a look at some of the other lammps molecule builders. Some of those tools support a wider range of force-fields than topotools or moltemplate do.
So, in summary, the way you generate improper interactions depends on which force-field you decide to use. This will probably determine which molecule-builder you decide to use as well.
------------ DETAILS (feel free to skip) -------------
-
Firsty, there are multiple conflicting definitions for what an “improper angle” is. (Here is yet another one. This is the definition that many of the improper_styles in LAMMPS currently use.)
-
Order of atoms varies. Which atom is the “central” atom? As you noticed, different molecule-builders assume different atom-type ordering. Molecule-builders designed to work with AMBER, CHARMm, and CVFF, the “central” atom is usually 2nd atom in the list of atom-types, and for OPLS, it is the 3rd atom in the lis:
a) Again, if you are using topotools to generate impropers, it assumes that the 3rd atom in the list of atom types is the “central” atom.
b) If you are using moltemplate, then by default, it assumes that the 1st atom is the “central” atom, although this can be overridden. -
Atom-type and Atom-ID order may differ.
To add improper interactions to a LAMMPS simulation, the molecule-builder software will create an “Impropers” section and insert it into the DATA file that LAMMPS will read. Each line in that section contains a list of 4 atoms, and the order that the atoms should appear in that section depends entirely on the “improper_style” that you have selected (which may also depend on the force-field you are using). Specifying the atoms in a different order will change the forces acting on those 4 atoms. However, the order that atoms appear in these lines may differ from the order they need to appear in the “atom-type-list” discussed above. LAMMPS does not know or care about how the atoms in an improper interaction are bonded together, or which atom is the central atom. The molecule-building software you use must get these details right for the force-field you are interested in. -
Different force-fields have different improper-angle symmetry. For example:
a) If you use the OPLSAA force-field, then you don’t want to generate an improper angle between atoms 1,2,3,4 if there is already an existing improper interaction for atoms 4,2,3,1, because the force acting on these atoms should be identical, because for improper_style opls, the formula that is used to calculate these force is invariant when you swap these atoms.
b) However if you are using the CVFF force-field, you DO want to create separate improper interactions for both 1,2,3,4 AND 4,2,3,1, because the formula for the energy does change when swapping atoms 1 and 4, in this case. (… at least that’s my impression.)
— customization ----
In topotools you can customize all these rule by editing the source code in topotools as you are doing.
In moltemplate, you can select from one of 4 different rules by choosing from one of the alternate files ending in “imp.py” in the “src/nbody_alternate_symmetry/” subdirectory. (See attached. Note: you can select between these force-field conventions without editing any code. See the OPLSAA and AMBER examples to see how this is done.) To create your own rule-set, select one of those files, and modify it to match your needs. These are short files. Only 13 lines long, excluding comments.)
Cheers
Andrew
P.S. Thanks to Arun and Fu-Chang Sun and others for helping me fix these issues, and to Paul Saxe for explaining them to me.
nbody_Impropers.py (2.85 KB)
opls_imp.py (3.3 KB)
gaff_imp.py (3.41 KB)
class2_imp.py (1.57 KB)