question on oplsaa_moltemplate.py

Hello,

I am wondering whether the stretch constant read from the OPLSAA.prm parameter file by the oplsaa_moltemplate.py script needs to be divided by 2 as below (x[2]/2).

—oplsaa_moltemplate.py file—

[snip]

g.write(“write_once(“In Settings”) {\n”)
index1=0
for x in bond:
for y in atom_lookup.get(x[0],[]):
for z in atom_lookup.get(x[1],[]):
g.write(“bond_coeff @bond:{}-{} harmonic {} {}\n”.format(y,z,x[2]/2,x[3]))
h.write("@bond:{0}-{1} @atom:{0} @atom:{1}\n".format(y,z))
g.write("} #(end of bond_coeffs)\n\n")

[snip]

I think the OPLSAA bond stretching potential is identical in form to the LAMMPS “harmonic” bond style and therefore does not need a division by 2.

Could you please comment on this?

Thanks,

Hello,

I am wondering whether the stretch constant read from the OPLSAA.prm
parameter file by the oplsaa_moltemplate.py script needs to be divided by 2
as below (x[2]/2).

---oplsaa_moltemplate.py file---
[snip]
g.write("write_once(\"In Settings\") {\n")
index1=0
for x in bond:
  for y in atom_lookup.get(x[0],[]):
    for z in atom_lookup.get(x[1],[]):
      g.write("bond_coeff @bond:{}-{} harmonic {}
{}\n".format(y,z,x[2]/2,x[3]))
      h.write("@bond:{0}-{1} @atom:{0} @atom:{1}\n".format(y,z))
g.write("} #(end of bond_coeffs)\n\n")
[snip]

I think the OPLSAA bond stretching potential is identical in form to the
LAMMPS "harmonic" bond style and therefore does not need a division by 2.

Could you please comment on this?

have you checked out this?
http://lammps.sandia.gov/doc/bond_harmonic.html

Yes. That form in the LAMMPS doc is the same as OPLSAA and therefore K should be what is in the oplsaa parameter file oplsaa.prm without dividing by 2. This is my understanding but wanted to hear from others to confirm.

Thanks,

Yes. That form in the LAMMPS doc is the same as OPLSAA and therefore K
should be what is in the oplsaa parameter file oplsaa.prm without dividing
by 2. This is my understanding but wanted to hear from others to confirm.

it has been a while since i wrote this, but you can check out this:
https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---part-2

the inputs there have been confirmed to be correct (multiple people
sent in bug reports for typos and little oversights), so you could
check if moltemplate reproduces those setups and would have an
indirect confirmation.

axel.

Yes. That form in the LAMMPS doc is the same as OPLSAA and therefore K
should be what is in the oplsaa parameter file oplsaa.prm without dividing
by 2. This is my understanding but wanted to hear from others to confirm.

Great discussion.
According to the supplemental information for the OPLSAA paper (free download)
http://pubs.acs.org/doi/abs/10.1021/ja9621760
OPLSAA uses this convention Ubond(r) = k*(r-r0)^2

As Axel pointed out, in LAMMPS, bond_style harmonic uses this convention.
I can only assume that the parameters in the "oplsaa.txt" file also
obey this convention. (That file is read by the oplsaa_moltemplate.py
script.)

Forgive me for posting code which I do not have the expertise to
evaluate. Disclaimer. I am not an expert using OPLSAA. Jason
Lambert graciously wrote the "oplsaa_moltemplate.py" script and
provided the "oplsaa.txt" file which it reads as an input file (after
some minor manual editing on the part of the user). I did some
googling, and I noticed that this file is very similar to the
"oplsaa.prm" file located at this web site:

http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm

I suspect it uses the same parameter convents.

However I just realized that the two files are not identical. (I am
emailing Jason now to ask about the difference.)

it has been a while since i wrote this, but you can check out this:
https://sites.google.com/site/akohlmey/software/topotools/topotools-tutorial---part-2

This discussion is useful to me. (Thanks for the warning about the
partial charges.)

the inputs there have been confirmed to be correct (multiple people
sent in bug reports for typos and little oversights), so you could
check if moltemplate reproduces those setups and would have an
indirect confirmation.

Yes. This is critical.

Thank you.
I am hoping people will post problems.
It's the only way to squash the bugs. I would be surprised if
everything works perfectly on the first try. (In addition to the atom
type numbers, I'm also worried about whether moltemplate constructs
data files with the correct number of improper interactions.)
Please post any bugs you find!

Thanks Valmor and Axel
...and Jason Lambert, of coarse!

Andrew
P.S.
There are a surprsingly large number of atom types (906) in the
"oplsaa.prm" file (and the "oplsaa.txt" file). I suspect this is
because of the large number of possible partial charges for each atom
type. (I say this because Jason's script assigns charge by atom type,
and it seems to generate neutral systems.)

Indeed, there was a serious bug in oplsaa_moltemplate.py causing
all of the bonds and all of the bond-angles to be too weak by a factor
of 2, and fixed another problem with the 1-4 interactions. Today I
posted a new version of moltemplate which contains a corrected version
of "oplsaa_moltemplate.py" which fixes this bug. If you plan to use
moltemplate to build systems using the OPLSAA force-field, please
update to the latest version of moltemplate (version 1.22 or later).

   Credit goes to Miguel Gonzales for digging through the details,
convincing me that there was a bug, and providing a code fix. Thanks
Miguel!

Andrew