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 <[email protected]> wrote: > > > 2017-02-13 18:24 GMT-07:00 Joshua Moore <[email protected]>: > >> 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/6db1e9f22f4db0f0967e2ed525af31 > 878cf73eaf> > Let me know if it fixes it. Otherwise, send me a patch for table_to_tab.pl > . > > Christoph > > >> 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: >> >> >> >> On Monday, February 13, 2017 at 6:11:08 PM UTC-5, Christoph Junghans >> wrote: >> >>> 2017-02-13 14:06 GMT-07:00 Joshua Moore <[email protected]>: >>> > 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/cb92bcd6b84a60b05698b5e >>> 409fc9d125959f03e> >>> >>> That said, the support for lammps as a simulation backend is still >>> pretty rudimentary in VOTCA. >>> >>> Christoph >>> >>> > >>> > 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 >>> > >>> > >>> > On Monday, February 13, 2017 at 3:59:10 PM UTC-5, Christoph Junghans >>> wrote: >>> >> >>> >> 2017-02-13 13:43 GMT-07:00 Joshua Moore <[email protected]>: >>> >> > 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 >>> >> > >>> >> > <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> >>> >> > >>> >> > >>> >> > >>> >> > On Monday, February 13, 2017 at 3:31:56 PM UTC-5, Christoph >>> Junghans >>> >> > wrote: >>> >> >> >>> >> >> 2017-02-13 12:56 GMT-07:00 Joshua Moore <[email protected]>: >>> >> >> > 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? >>> >> >> 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 [email protected]. >>> >> >> > To post to this group, send email to [email protected]. >>> >> >> > 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 >>> >> > >>> >> > -- >>> >> > 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 [email protected]. >>> >> > To post to this group, send email to [email protected]. >>> >> > 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 >>> > >>> > -- >>> > 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 [email protected]. >>> > To post to this group, send email to [email protected]. >>> > 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 >>> >> -- >> 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 [email protected]. >> To post to this group, send email to [email protected]. >> 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 > > -- > 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 > [email protected]. > To post to this group, send email to [email protected]. > > Visit this group at https://groups.google.com/group/votca. > For more options, visit https://groups.google.com/d/optout. > -- 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 [email protected]. To post to this group, send email to [email protected]. Visit this group at https://groups.google.com/group/votca. For more options, visit https://groups.google.com/d/optout.
