dihedral_style/opls parameters for lammps

Hello everyone,

I have been having a hard time figuring out where to get the OPLS dihedral parameters for lammps. I see in the OPLS dihedral form used by lammps, one specifies the V1 V2 V3 and V4 in the dihedral_coeff command. However, in all the OPLS parameter sets I have looked at, the one given in tinker being what I think is close to the latest version of OPLS (http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm), I only see a mention of what I believe to be V1 V2 V3 and f1 f2 f3 where the latter are the phase angle offsets. The Watkins paper cited in the lammps manual only provides parameters for perfluouroalkanes, so I was wondering where I need to go to obtain these parameters for other molecules, or is there some conversion between the V1 V2 V3 f1 f2 f3 to the V1 V2 V3 V4 sets that I am missing.

thanks,

conor

Potential parameters are inputs to LAMMPS. We only provide a few

potentials in the potentials dir. So generally, you have to find
them in the literature or write them out from other software.

Steve

I have been having a hard time figuring out where to get the OPLS dihedral
parameters for lammps. I see in the OPLS dihedral form used by lammps, one
specifies the V1 V2 V3 and V4 in the dihedral_coeff command. However, in all
the OPLS parameter sets I have looked at, the one given in tinker being what
I think is close to the latest version of OPLS
(http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm), I only see a
mention of what I believe to be V1 V2 V3 and f1 f2 f3 where the latter are
the phase angle offsets.

...

or is there some
conversion between the V1 V2 V3 f1 f2 f3 to the V1 V2 V3 V4 sets that I am
missing.

You are asking for details about the format of the "oplsaa.prm" file,
which is is not part of LAMMPS. It is a file distributed with the
TINKER software, and used by MOLTEMPLATE.
Since I wrote moltemplate, I will do my best to answer.
If you discover that my answer is incorrect, please let me know.

However, this file is distributed by the Ponder lab, and they are the
authoritative source for information about this file.

Anyway, in the 2001 JPCA paper cited by the LAMMPS manual,
(http://courses.chem.psu.edu/chem408/reading/MM_topics/jorgensen_perfluoro_jpca_2001.pdf)
the formula for the torsion energy is:

0.5*(V1*(1+cos(x+f1) + V2*(1-cos(2x+f2)) + V3*(1+cos(3x+f3)) +
V4*(1-cos(4x+f4)))

There is no general way to convert this formula to the formula used by
LAMMPS' "dihedral_style opls", for arbitrary values of f1,f2,f3,f4

http://lammps.sandia.gov/doc/dihedral_opls.html

For dihedral style opls, the energy is:
0.5*(K1*(1+cos(x)) + K2*(1+cos(2x)) + K3*(1+cos(3x)) + K4*(1+cos(4x)))

(Note there is a sign change in the two formulas.: "1-cos(2x+f2)"
instead of "1+cos(2x)".)

However, it turns out that, for every torsion interaction in OPLSAA:
f1=0
f2=180 # (this negates the sign change)
f3=0
f4=180
(See the torsion section of that "oplsaa.prm" file.) Consequently,
the "dihedral_style_opls" formula is general enough, and the
conversion you are looking for is:

K1 = V1
K2 = V2
K3 = V3
K4 = V4

(This is what Jason Lambert's "oplsaa_moltemplate.py" script does when
it converts the "oplsaa.prm" file into moltemplate format
(oplsaa.lt).)

Finally, I am having a hard time interpreting the parameters in the torsional part of
the oplsaa file. For instance,

torsion 13 3 20 13 4.669 0.0 1 5.124 180.0 2 0.000 0.0 3

I am assuming 4.669,5.124,and 0.000 correspond to a V1 V2 V3. Is this true?

Yes.

What are the 1 2 3 referring to along with the 0.0 180.0 0.0 (I am assuming these are a f1,f2,f3).

Yes.

(The "1", "2" and "3", are the integer frequencies in the Fourier
expansion. They never vary in the "oplsaa.prm" file. They are always
1,2,3.)

Unfortunately, the "oplsaa.prm" file does not include V4 parameters.
However, they are not needed if you are using the original, basic
OPLSAA force-field. In that case V4=0. (Apparently V4 is non-zero
for the perfluouroalkanes discussed in the 2001 paper, but not in the
original 1996 OPLSAA paper. See below.)

The Watkins paper cited in the lammps manual only
provides parameters for perfluouroalkanes, so I was wondering where I need
to go to obtain these parameters for other molecules,

The original OPLS paper is from 1996, and it is not only for perfluouroalkanes
http://pubs.acs.org/doi/abs/10.1021/ja9621760

It uses the same formula for the dihedral energy.

Incidentally, you do not have to use "dihedral_style opls" to simulate
molecules using the OPLS force-fields in LAMMPS. If you want a more
general formula, I recommend "dihedral_style fourier".
http://lammps.sandia.gov/doc/dihedral_fourier.html
That dihedral style will allow you to customize the f1,f2,f3,f4,...
offsets, and the frequencies as well (1,2,3,4,...)

If you need more detail than this, the TINKER mailing list is probably
the best source of information.

Andrew

Thank you so much for that thorough response. I apologize for bringing a tinker question to the lammps email list. I honestly thought lammps was doing something different because of that mysterious V4. Specifically about lammps, the dihedral_opls function I am looking at is (http://lammps.sandia.gov/doc/dihedral_opls.html):

0.5k1(1+cosx) + 0.5k2(1-cos2x)+0.5k3(1+cos3x)+0.5k4(1-cos4x)

The second and fourth term in the user manual already has the sign flipped, unlike the one you wrote. How would this still affect the sign flip you mentioned by the 180 degrees f2?

conor

Thank you so much for that thorough response. I apologize for bringing a
tinker question to the lammps email list.

No, that's okay.
I think by raising this issue, you might have found a bug in Jason's
"oplsaa_moltemplate.py" script (See below.)
If so, this could be be effecting many people. Let's fix it.

I honestly thought lammps was
doing something different because of that mysterious V4. Specifically about
lammps, the dihedral_opls function I am looking at is
(http://lammps.sandia.gov/doc/dihedral_opls.html):

0.5*k1*(1+cosx) + 0.5*k2*(1-cos2x)+0.5*k3*(1+cos3x)+0.5*k4*(1-cos4x)

You are right and I was wrong.
You kept pointing out there was a problem with the f1,f2,f3
parameters, and, for some reason, I kept not seeing it. (I kept
unconciously misreading the LAMMPS documentation the way I wanted it
to be.)

Thank you for your persistence in pointing this out to me, despite my
dismissals.

I have made some changes to "oplsaa_moltemplate.py" and I will test
them this weekend and release the update in the next few days.

Andrew

Awesome. Glad I could help out! I thought there was some mysterious trigonometric identity everyone else was in on that I had no idea about. So does this mean the convenient V1 = K1, V2 = K2, V3 = K3, V4 =K4 and ignore the f1,f2,f3 protocol will no longer be correct?

conor

You should see the tutorial on Axel Kohlmeyer's home page where he does
several small hydrocarbons using OPLSAA:

https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---part-2

You can see the parameters he uses (in.step2b and in.step2e_2) are
identical to 1-1-1-1 parameters in the oplsaal.prm file that comes with
Tinker5.

The lammps documentation for dihedral_style opls is correct for the
equations but the numeric example is total nonsense and worse than
useless. The values listed of 1 90.0 90.0 90.0 would give 6*90 kcal of
torsional energy for the V2 term in a normal single bond. A normal sp3-sp3
single bond has six torsions and each one would contribute
90*(1+cos2(phi)).

kevin

The lammps documentation for dihedral_style opls is correct for the
equations but the numeric example is total nonsense and worse than
useless. The values listed of 1 90.0 90.0 90.0 would give 690 kcal of
torsional energy for the V2 term in a normal single bond. A normal sp3-sp3
single bond has six torsions and each one would contribute
90
(1+cos2(phi)).

Please suggest a good example for the doc page and we’ll replace it.
I don’t use OPLS.

Steve

The lammps documentation for dihedral_style opls is correct for the
equations but the numeric example is total nonsense and worse than
useless. The values listed of 1 90.0 90.0 90.0 would give 6*90 kcal of
torsional energy for the V2 term in a normal single bond. A normal sp3-sp3
single bond has six torsions and each one would contribute
90*(1+cos2(phi)).

Please suggest a good example for the doc page and we'll replace it.
I don't use OPLS.

here are the numbers for the hydrocarbon inputs in the topotools
tutorial (in real units):

dihedral_coeff 1 1.740 -0.157 0.279 0.00 # CT-CT-CT-CT
dihedral_coeff 2 0.000 0.000 0.366 0.000 # CT-CT-CT-HC
dihedral_coeff 3 0.000 0.000 0.318 0.000 # HC-CT-CT-HC

Thank you so much for pointing out that tutorial. I had skimmed it previously, but missed the link to the input script. Just for clarification for ethane, the H is type 46 and the C is type 13 I believe. The torsion angle listed in the oplsaa.prm file is

torsion 46 13 13 46 0.000 0.0 1 0.000 180.0 2 0.300 0.0 3

What am I missing for the choice of the 0.318 value for the K3?

conor

Thank you so much for pointing out that tutorial. I had skimmed it
previously, but missed the link to the input script. Just for clarification
for ethane, the H is type 46 and the C is type 13 I believe. The torsion
angle listed in the oplsaa.prm file is

torsion 46 13 13 46 0.000 0.0 1 0.000 180.0 2 0.300 0.0 3

What am I missing for the choice of the 0.318 value for the K3?

that you should ask the tinker folks. i looked up the parameters in
the referenced publication.

axel.

Axel has already provided good number for oplsaa. I looked at the other
dihederal_style commands and I think they are all not correct for the
standard organic force fields (oplsaa, mmff95, amber, charmm, mm2, mm3,
class2). All of these force fields use some variant of the sum of a series
of K*cos(n*phi), where K is the energy of the torsional barrier and
usually given in kcal. Typical examples for some of these force fields
are:

CT-CT-CT-CT - backbone carbons in butane
          K1 K2 K3 K4
oplsaa: 1.74 -.157 .279 0
amber: 0 .18 .03 0
charmm: 0 .20 .03 0
mm3: .185 .170 .520 0
mmff94: .130 .681 .332 0

None of these force fields use the K4 term. You can't use the numbers
interchangeably because the 1,4 interaction depends on the vdw + charge +
torsion contributions, but the numerical values are all similar. This
applies to dihedral_styles charmm, class2, harmonic and opls. I don't know
anything about the helix style. It could be that the numeric values in the
documentation do apply to some coarse grained force fields and it would be
nice if someone would comment on this.

kevin

Hello developers,

I want to use the fix rigid/nve command to simulate the rigid body, which consists of many point particles (on spherical surface) and a dipole paricle (in the center of the spherical body). I hope the orientation of the dipole particle will be update with the rotation of the whole rigid body. However, I do not know whether the fix rigid/nve command could perform what i want.

Could anyone tell me and which is the subfunction in Fix_rigid.cpp for the aim.

Thanks in advane.

------------------发送

Best wishes from Suzhou, China

Wade

I believe fix rigid/nve should update the dipole. By inheriting

this bit of code from fix rigid:

if (eflags[i] & DIPOLE) {
MathExtra::quat_to_mat(quat[ibody],p);
MathExtra::matvec(p,dorient[i],mu[i]);
MathExtra::snormalize3(mu[i][3],mu[i],mu[i]);
}

You can check by dumping the dipole moments of one or more
atoms in rigid bodies, and verifying they are chaning as the rigid
body rotates.

Steve

I looked at the other
dihederal_style commands and I think they are all not correct for the
standard organic force fields (oplsaa, mmff95, amber, charmm, mm2, mm3,
class2).

If you want to suggest better coeff values for the doc page examples

for each of those dihedral styles, great. But you need to

give them to us in the format that matches the coeffs those
styles take. None of them are K1,K2,K3,K4.

Steve