Center atom for improper dihedral coefficients in oplsaa.it

[HUANGZR86]

Dear developer,

I am using the OPLS-AA parameters from the file oplsaa.lt. May I ask for the improper dihedral X_X_003_004, which one should be the center atom?

Best regards,
Huang Ziru

[jewettaij]

Hi Huang

For OPLSAA in moltemplate, the 3rd atom is the center atom.
(To be clear, this applies to the the 4 atoms listed on each line of the "Impropers" section of the "system.data" file generated by moltemplate.)

Details

Your example:

    improper_coeff @improper:X_X_003_004  10.5 180.0  # (moltemplate)
    # <==> "imptors 0 0 3 4  21.000 180.0 2" in oplsaa.prm (tinker)
    # <==> "define improper_O_C_X_Y 180.0 43.93200 2" in oplsaa.ff (gromacs)
   :
    @improper:X_X_003_004 @atom:* @atom:* @atom:*_b*_a*_d*_i003* @atom:*_b*_a*_d*_i004*

...In this case, the central atom corresponds to any atom that belongs to the "i003" group. (These groups are defined in the "replace" commands in the "oplsaa.lt" file. Multiple atoms belong to this group, as shown here.)

To create the "oplsaa.lt" file used by moltemplate, I converted the "oplsaa.prm" file which comes with TINKER into moltemplate (.lt) format. I think that TINKER also uses the 3rd atom. (But I could be wrong. Please let me know if this is not the case.)

Years ago, several moltemplate users and I spent days trying to confirm that this is the correct order. But this was impossible to confirm because some molecule builders (like BOSS and Materials Studio) do not document how the atom order is determined or where the central atom is. For this reason, we could not always reproduce the atom order created by these programs.

If this sounds like a serious problem. However perhaps it is not a problem: Fortunately, all of the improper interactions in OPLSAA are even functions((U(φ)=U(-φ)), and U(φ) has a minima at φ=0, which only occurs when the 4 atoms are coplanar (in the same plane). Changing the order of the 4 atoms might have the effect of strengthening or weakening this harmonic force that keeps them in the same plane (but probably not by more than a factor of 2). This will have hardly any effect on the flexibility of the molecule, as long as all 4 atoms remain in the same plane (or nearly in the same plane). It turns out that the flexibility of many molecules is mostly determined by the strength of the dihedral interactions (and also the 3-body "angle" interactions). So even if the improper interactions have the wrong atom order, as long as all 4 atoms remain nearly coplanar, then the mechanical behavior of the molecule should not be strongly affected.

That said, I could be wrong. If anyone thinks we made a mistake in our implementation of OPLSAA, then please let me know.

If you are curious how moltemplate works, see the "Impropers by Type" section of the "oplsaa.lt" file.

write_once("Data Impropers By Type (opls_imp.py)") {
  ...
}

Notice the "(opls_imp.py)". This tells moltemplate that the "opls_imp.py" file defines the central atom and the bond topology and symmetry for these kinds of interactions. Here is the relevant excerpt of the contents of the that file:

#    To find 4-body "improper" interactions,
#    (by default, most of the time), we would use this subgraph:
#           0
#           *                  1st bond connects atoms 2 and 0
#           |              =>  2nd bond connects atoms 2 and 1
#         _.*._                3rd bond connects atoms 2 and 3
#       *'  2  `*
#      1         3
#
# In OPLS (as implemented in TINKER .PRM files), the central atom
# is the third atom ("2").
bond_pattern = Ugraph([(2,0), (2,1), (2,3)])

Atom "2" is the central atom. (However atom "2" refers to the 3rd atom because the numbers begin at 0, not 1.)

Hope this helps.
-A