Hello Andrew.
I had a look at gaff.lt and gaff.dat. The atoms are re-ordered for impropers compared to gaff.dat for my system.
Moltemplate:
@improper:n-c-c3-hn @atom:n @atom:c @atom:c3 @atom:hn
14 13 15 36 (atom ids)
Gaff.dat
c -c3-n -hn 1.1 180. 2. Junmei et al.1999 (n is the central atom)
13 15 14 36 (atom ids)
Why is that so ?
As such, you do not need to re-order the atoms in gaff.dat for the impropers (if you are using improper style cvff as moltemplate is using)
Consider the case when you have an improper IJKL, and and you use the improper style cvff.
Angle computed between the planes , IJK and JKL is the same when the central atom is K or when the central atom is I (as it is in lammps for improper cvff).
The improper style does not require the bond connectivity information to determine the angle between the planes.
In my opinion the order should be retained as is appears in gaff.dat
Geometrically the angle between the planes are not the same (13 15 36 15 14 36 (if you use moltemplate) ).
Can you please explain why the atom is being re-arranged in moltemplate for my system.
Thanks
Dear Arun
I'm not certain, but I suspect you have found a bug.
I changed the order in response to an email from Paul Saxe who
pointed out that in AMBER/GAFF, the third atom ("K") is assumed to be
the central atom, not the first atom ("I"). I thought that the change
I made at the time was enough, but I may be printing out the atoms in
the DATA file in the wrong order. (More specifically, in the
"Impropers" section of the data file). I was confused and I was under
the impression that this would have no effect on the angle between the
planes (up to a sign change). But I realize now after today's email
correspondence that I was (or am) confused.
Fortunately, if you have used moltemplate to generate your LAMMPS
files, I don't think this will cause huge problems because, either
way, the force-field will constrain the same 4 atoms to be coplanar.
However, if there is a bug, the flexibility of the molecule may be
effected. In that case, it really needs to be fixed. I'll think
about how to do this, but I might have wait a couple weeks before I
can implement anything.
From his email to me, Paul Saxe seems has been thinking about these
kinds of issues longer than I have. I'm not trying to shrug off my
own responsibilities with moltemplate. Paul's company sells MedeA
licenses, and this software is likely to have ironed out these kinds
of bugs. (And I certainly appreciate him helping me fix this issue in
moltemplate.)
Anyway I will fix this issue in moltemplate soon.
Feeling embarrassed.
Thanks Arun and Paul
Andrew
The issue is that moltemplate does not currently allow the
atom-order to printed in the DATA file to be different from the atom
order used for atom-type pattern matching (which generates the list of
improper interactions automatically). This would have to be
generalized to make it work with all the various force-field
conventions, and allow moltemplate users to select between them.
Incidentally, this bug in moltemplate (reported by Arun Srikant
Sridhar) was fixed in 2014-12-19. I announced it back then, but I
forgot to mention it on this thread.
I was searching google for gaff.dat, and I accidentally noticed that
this old thread has been left dangling, unresolved. Please forgive
the duplicate post.
All improper interactions (and other interactions) in GAFF and
OPLSAA should now be accurate.
Thanks again to Arun for reporting the bug.
Andrew
P.S. More generally, if you are attempting to write your own
force-field, or port an existing forcefield into moltemplate format,
you can now customize the atom order (and symmetry conventions), by
copying the tools/moltemplate/src/nbody_alternate_symmetry/gaff_imp.py
file to a new file in the same directory (for example, named
"new_ff.py"). Edit that file to specify the central atom, and make
other changes. Later, in the LT file containing your force-field
parameters, refer to this file ("new_ff.py") using this syntax:
write_once("Data Impropers By Type (new_ff.py)")
@improper:X-o-c-o @atom:* @atom:o @atom:c @atom:o
@improper:X-X-c-o @atom:* @atom:* @atom:c @atom:o
@improper:X-X-ca-ha @atom:* @atom:* @atom:ca @atom:ha
@improper:X-X-n-hn @atom:* @atom:* @atom:n @atom:hn
:
}
You can do the same thing with dihedrals. (For example, unlike most
force-fields, the formula for the energy of class2 dihedrals is not
symmetric with respect to atom-order-reversal. Energy(I,J,K,L) !=
Energy(L,K,J,I), so we would want moltemplate to generate two
different dihedral interactions for the same for atoms. Too see how
this is implemented in moltemplate, compare these two files:
"nbody_Dihedrals.py" with "nbody_alternate_symmetry/class2_dih.py", in
the "tools/moltemplate/src/" directory.)
Cheers