LAMMPS Angle Potential in Degrees

[Joshua Moore]

Hello,

How do people handle the fact that the LAMMPS angle potential is in degrees and not radians?

If I convert the angle distribution to degrees, will csg_inverse calculate the correct comparison?

Thanks.

Josh

[Christoph Junghans]

Before running the actual simulation in each iteration step,
csg_inverse will convert its internal format (.pot file) to a LAMMPS
table (name given in the xml file), it will do the conversion from
radians to degrees by default.
There is also a scaling factor for the non-bonded interactions in case
you need to convert from nm to Angstrom.

Christoph
>
> Thanks.
>
> Josh
>
> --
> You received this message because you are subscribed to the Google Groups
> "votca" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to votca+un...@googlegroups.com.
> To post to this group, send email to vo...@googlegroups.com.
> Visit this group at https://groups.google.com/group/votca.
> For more options, visit https://groups.google.com/d/optout.



--
Christoph Junghans
Web: http://www.compphys.de

[Joshua Moore]

Hi Christoph,

Thank you for the quick reply.  Maybe I am not setting something right in the settings.xml file, because the angle.pot file in step000 is still from 0 to 3.1415.  Also my target distribution function is from 0 to 3.1415 in radians, calculated with csg_stat.  Do I need to convert that first before feeding it through csg_inverse to LAMMPS?

  <bonded>
    <!-- name of the interaction -->
    <name>angle</name>
    <min>0.0</min>
    <max>3.1415</max>
    <step>0.01</step>
    <inverse>
      <kBT>0.59616108</kBT>
      <!-- target distribution -->
      <target>angle.dist.tgt</target>
      <lammps>
      <kBT>0.59616108</kBT>
        <table>table_a1.pot</table>
      </lammps>
    </inverse>
  </bonded>

[Christoph Junghans]

2017-02-13 13:43 GMT-07:00 Joshua Moore <jdm...@ncsu.edu>:
> Hi Christoph,
>
> Thank you for the quick reply. Maybe I am not setting something right in
> the settings.xml file, because the angle.pot file in step000 is still from 0
> to 3.1415. Also my target distribution function is from 0 to 3.1415 in
> radians, calculated with csg_stat. Do I need to convert that first before
> feeding it through csg_inverse to LAMMPS?

step_000 doesn't run any simulation and hence there is no
table_a1.pot, look in step_001.

Christoph

[Joshua Moore]

Thanks.  I see now in step001.  It does convert the table_a1.pot file as you said.

The issue seems to be that the LAMMPS potential file generated converts the min/max from radians to degrees, but the <step>0.01</step> isn't converted, so it tries to write out > 180 degrees in the LAMMPS file (18001 points) or up to 18001*0.01*180/3.1415 = 10314.12 degrees.  Is there a special way to specify this step size after the unit conversion?

Thanks for all your help with this.

Here is tail of the table_a1.pot LAMMPS file.

17989  1.03064e+04   7.9841910e-04   6.9526739e-04
17990  1.03069e+04   3.9962120e-04   6.9677605e-04
17991  1.03075e+04   0.0000000e+00   6.9814120e-04
17992  1.03081e+04  -4.0036230e-04   6.9936284e-04
17993  1.03087e+04  -8.0138340e-04   7.0044097e-04
17994  1.03092e+04  -1.2029811e-03   7.0137559e-04
17995  1.03098e+04  -1.6050732e-03   7.0216670e-04
17996  1.03104e+04  -2.0075775e-03   7.0281430e-04
17997  1.03109e+04  -2.4104117e-03   7.0331839e-04
17998  1.03115e+04  -2.8134936e-03   7.0367897e-04
17999  1.03121e+04  -3.2167410e-03   7.0389604e-04
18000  1.03127e+04  -3.6200716e-03   7.0396961e-04
18001  1.03132e+04  -4.0234033e-03   7.0389966e-04

[Christoph Junghans]

2017-02-13 14:06 GMT-07:00 Joshua Moore <jdm...@ncsu.edu>:
> Thanks. I see now in step001. It does convert the table_a1.pot file as you
> said.
>
> The issue seems to be that the LAMMPS potential file generated converts the
> min/max from radians to degrees, but the <step>0.01</step> isn't converted,
> so it tries to write out > 180 degrees in the LAMMPS file (18001 points) or
> up to 18001*0.01*180/3.1415 = 10314.12 degrees. Is there a special way to
> specify this step size after the unit conversion?

Sorry, I just realized there is a bug the potential conversion for
lammps. The conversion factor is applied twice.
3.1415*(180/3.1415)**2=10313.5.

I fixed that now:
<https://github.com/votca/csg/commit/cb92bcd6b84a60b05698b5e409fc9d125959f03e>

That said, the support for lammps as a simulation backend is still
pretty rudimentary in VOTCA.

Christoph

[Joshua Moore]

Hi Christoph,

Thanks again for your quick message.

It seems to run ok after the update.

However, I think the forces are not correct in the LAMMPS table files.  I might be mistaken.

I'm attaching the input files needed to generate the run with (toVotca.tar.gz)

csg-inverse --options settings.xml

Below is a plot of what the Bond and Angle force from the generated table files (table_b1.pot and table_a1.pot) after step001.  I'm taking a simple derivative as -deltaU/delta_r or -deltaU/delta_theta, and they don't seem to match up with what VOTCA is outputting for the forces in the file.  I believe for the angle, that LAMMPS is expecting -dU/dtheta.

Consequently, the temperature blows up nearly immediately, and then you see warnings in LAMMPS like "WARNING: 87 of 91 force values in table are inconsistent with -dE/dr.  Should only be flagged at inflection points (../bond_table.cpp:380)".  For some reason the angle potential is not triggering this warning, but I believe it should because it also does not appear to be consistent.  (Note: checking the lammps source, it appears this check is not in the angle table code, so that's probably the reason.)

The energies and forces from the pair coefficients appear to be correct.

Any ideas about what I might be doing wrong, or a potential fix?  I'm also not getting perfect agreement to the potential of mean force (e.g., -kTln(Pr))for the step001 angle and bond potentials either, but this might partially be due to the shifting you are doing.

Thanks again for your help.

Josh

Bonds:

Angles:

[Christoph Junghans]

2017-02-13 18:24 GMT-07:00 Joshua Moore <jdm...@ncsu.edu>:

Hi Christoph,

Thanks again for your quick message.

It seems to run ok after the update.

However, I think the forces are not correct in the LAMMPS table files.  I might be mistaken.

I'm attaching the input files needed to generate the run with (toVotca.tar.gz)

csg-inverse --options settings.xml

Below is a plot of what the Bond and Angle force from the generated table files (table_b1.pot and table_a1.pot) after step001.  I'm taking a simple derivative as -deltaU/delta_r or -deltaU/delta_theta, and they don't seem to match up with what VOTCA is outputting for the forces in the file.  I believe for the angle, that LAMMPS is expecting -dU/dtheta.

We had a discussion about LAMMPS' table format a while ago here:

<https://groups.google.com/d/msg/votca/xGJnFSrFXYQ/Thg8kIyoYBoJ>

and the code was written by Frank here:

<https://groups.google.com/d/msg/votca/Jxsdz6v9QKw/OUVfYan8CwAJ>

I think, for bonds the 4th column is  -r*deltaU/delta_r, while for angle and dihedral it is just the negative derivative, but I guess the code was never really tested.

Have a look at:

< https://github.com/votca/csg/commit/6db1e9f22f4db0f0967e2ed525af31878cf73eaf>

Let me know if it fixes it. Otherwise, send me a patch for table_to_tab.pl.

Christoph

To unsubscribe from this group and stop receiving emails from it, send an email to votca+unsubscribe@googlegroups.com.

To post to this group, send email to vo...@googlegroups.com.
Visit this group at https://groups.google.com/group/votca.
For more options, visit https://groups.google.com/d/optout.

[Joshua Moore]

Hi Christoph,

Thanks for pointing me to these threads and the code.  I was trying to find it.  I didn't realize that's what you changed above for the angle.

I'll take a look and see if I can get it worked out.  It looks like it's all there already just needs tweaked a bit.  I think the bond force is in force units without the r and the angle issue might be because of the derivative with respect to theta in degrees.  That might be why you had that unit conversion in there previously for the angle but then it got double counted.  

It looks like it can be worked out easily.  I'll keep you posted.

Thanks.

Josh

You received this message because you are subscribed to a topic in the Google Groups "votca" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/votca/cmXf_FMEvKs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to votca+unsubscribe@googlegroups.com.

[Joshua Moore]

Hi Christoph,

I am still testing, but it appears to be correct if you remove the $r[$i] multiplying the force in table_to_tab.pl for bonds and angles.  The unit conversions seem to be correct after all.  I haven't tested dihedrals yet though.  I learned LAMMPS can handle dihedrals with degrees or radians but angles only with degrees in the table files but requires a term in the header. 

So change this line  in table_to_tab.pl for bonds in the LAMMPS section

printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]*$r[$i]);

to this

printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]);

It looks like you might have changed it for angles last night already? You can probably recombine bonds, angles, dihedrals, as this line is the same for all of them, although haven't tested dihedrals yet.

The LAMMPS manual also "Currently" says the force is in force units for the bond tables. It may very well have been different before.

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

Thanks. I am getting good agreement now to ibi for the hexane example with LAMMPS for pair potentials, bonds, and angles, although still running it longer. The temperature appears stable. I've run each iteration up to 250,000 MD steps so far.

Josh

On Tue, Feb 14, 2017 at 12:44 AM, Christoph Junghans <jung...@votca.org> wrote:

You received this message because you are subscribed to a topic in the Google Groups "votca" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/votca/cmXf_FMEvKs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to votca+unsubscribe@googlegroups.com.

[Christoph Junghans]

2017-02-14 17:58 GMT-07:00 Joshua Moore <jdm...@ncsu.edu>:
>
> Hi Christoph,
>

> I am still testing, but it appears to be correct if you remove the $r[$i] multiplying the force in table_to_tab.pl for bonds and angles. The unit conversions seem to be correct after all. I haven't tested dihedrals yet though. I learned LAMMPS can handle dihedrals with degrees or radians but angles only with degrees in the table files but requires a term in the header.
>
> So change this line in table_to_tab.pl for bonds in the LAMMPS section
>
> printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]*$r[$i]);
>
> to this
>
> printf(OUTFILE "%i %12.5e %15.7e %15.7e\n",$i+1,$r[$i],$pot[$i],-$pot_deriv[$i]);
>
> It looks like you might have changed it for angles last night already? You can probably recombine bonds, angles, dihedrals, as this line is the same for all of them, although haven't tested dihedrals yet.
>
> The LAMMPS manual also "Currently" says the force is in force units for the bond tables. It may very well have been different before.
>
> http://lammps.sandia.gov/doc/bond_table.html
>
> Thanks. I am getting good agreement now to ibi for the hexane example with LAMMPS for pair potentials, bonds, and angles, although still running it longer. The temperature appears stable. I've run each iteration up to 250,000 MD steps so far.

I fixed the forces for bond now:
<https://github.com/votca/csg/commit/567a37802a748bf3a73345f39fb02d1a5903aabe>

Christoph

>>> Angles:

>> You received this message because you are subscribed to a topic in the Google Groups "votca" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/votca/cmXf_FMEvKs/unsubscribe.

>> To unsubscribe from this group and all its topics, send an email to votca+un...@googlegroups.com.