I was trying to use the GAFF force field in LAMMPS, and I used the amber2lammp.py tool to convert the Amber data file to a LAMMPS data file.
One thing I noticed about the dihedral parameters is that in GAFF, the dihedral interaction is defined with a "1/2" factor in the constant term (see
The only complete documentation I found to explain the format of the
gaff.dat file was here:
http://ambermd.org/formats.html
According to this documentation, the dihedral energy U(phi) is:
(PK/IDIVF) * (1 + cos(PN*phi - PHASE))
PN = the last column of gaff.dat (dihedrals section)
PHASE = the 2nd-to-last column of gaff.dat
PK = the 3rd-to-last column of gaff.dat
IDIVF = the 4th-to-last column of gaff.dat
If I am using the wrong formula (interpreting gaff.dat the wrong way),
please let me know.
Hello Wei. Amber2lammps.py is correct.
Hi Arun
I'm curious to know if the dihedral coeffs generated using
moltemplate's gaff.lt look correct to you. I'm using some confusing
documentation to generate them (see above) and I'd love to know if I
made a mistake. (I am working on fixing improper atom order, by the
way. Thanks for that bug report.)
GAFF), and in the meantime, the dihedral (harmonic style) in LAMMPS doesn't have the "1/2". So intuitively, for one same dihedral, the value showing in a Amber data file should be 2 times larger than that in a LAMMPS data file. However, the LAMMPS data file generated by the amb2lmp tool gives the same value as in the Amber .top file.
Is it correct? I also checked the parameter values in the "gaff.lt" file in moltemplate (but didn't use it), it should also yield the same value as the amb2lmp tool.
--- regarding moltemplate ---
The "gaff.lt" file which comes with moltemplate uses "dihedral_style
fourier" instead of "dihedral_style harmonic". The original gaff.dat
is converted to gaff.lt using the formula above, along with the
documentation for dihedral_style fourier:
Btw, is there some way to compute the energy of one specific dihedral in LAMMPS. I noticed the "compute pe/atom" command, but one atom is always involved in more than one dihedrals. How to deal with that?
I don't know how to calculate this during the simulation without
editing the LAMMPS code. However it's not difficult to do this after
you run the simulation (post-processing).
After running the simulation, assign the dihedral you care about to
a different type number (by modifying the data file, or using the set
command).
Then use the LAMMPS' "rerun" command to read a trajectory that you
created earlier, and recalculate the dihedral energy at each frame in
the trajectory
When you do this, set all the coeffs for other dihedral interactions
to zero (leaving the coeffs for that dihedral alone). Then print the
total dihedral energy using something like:
thermo_style custom step edihed
Cheers
Andrew