new dihedral table/cut command

doc/src/dihedral_table_cut.txt

@@ -0,0 +1,230 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+dihedral_style table/cut command :h3
+dihedral_style table/cut/omp command :h3
+
+[Syntax:]
+
+dihedral_style table/cut style Ntable :pre
+
+style = {linear} or {spline} = method of interpolation
+Ntable = size of the internal lookup table :ul
+
+[Examples:]
+
+dihedral_style table/cut spline 400
+dihedral_style table/cut linear 1000
+dihedral_coeff 1 aat 1.0 177 180 file.table DIH_TABLE1
+dihedral_coeff 2 aat 0.5 170 180 file.table DIH_TABLE2 :pre
+
+[Description:]
+
+The {table/cut} dihedral style creates interpolation tables of length
+{Ntable} from dihedral potential and derivative values listed in a
+file(s) as a function of the dihedral angle "phi".  In addition, an
+analytic cutoff that is quadratic in the bond-angle (theta) is applied
+in order to regularize the dihedral interaction.  The dihedral table
+files are read by the "dihedral_coeff"_dihedral_coeff.html command.
+
+The interpolation tables are created by fitting cubic splines to the
+file values and interpolating energy and derivative values at each of
+{Ntable} dihedral angles. During a simulation, these tables are used
+to interpolate energy and force values on individual atoms as
+needed. The interpolation is done in one of 2 styles: {linear} or
+{spline}.
+
+For the {linear} style, the dihedral angle (phi) is used to find 2
+surrounding table values from which an energy or its derivative is
+computed by linear interpolation.
+
+For the {spline} style, cubic spline coefficients are computed and
+stored at each of the {Ntable} evenly-spaced values in the
+interpolated table.  For a given dihedral angle (phi), the appropriate
+coefficients are chosen from this list, and a cubic polynomial is used
+to compute the energy and the derivative at this angle.
+
+The following coefficients must be defined for each dihedral type via
+the "dihedral_coeff"_dihedral_coeff.html command as in the example
+above.
+
+style (aat)
+cutoff prefactor
+cutoff angle1
+cutoff angle2
+filename
+keyword :ul
+
+The cutoff dihedral style uses a tabulated dihedral interaction with a 
+cutoff function
+
+:c,image(Eqs/dihedral_table_cut.jpg)
+
+The cutoff specifies an prefactor to the cutoff function.  While this value
+would ordinarily equal 1 there may be situations where the value should change.
+
+The cutoff angle1 specifies the angle (in degrees) below which the dihedral
+interaction is unmodified, i.e. the cutoff function is 1.
+
+The cutoff function is applied between angle1 and angle2, which is the angle at
+which the cutoff function drops to zero.  The value of zero effectively "turns
+off" the dihedral interaction.
+
+The filename specifies a file containing tabulated energy and
+derivative values. The keyword specifies a section of the file.  The
+format of this file is described below.
+
+:line
+
+The format of a tabulated file is as follows (without the
+parenthesized comments).  It can begin with one or more comment
+or blank lines.
+
+# Table of the potential and its negative derivative  :pre
+
+DIH_TABLE1                   (keyword is the first text on line)
+N 30 DEGREES                 (N, NOF, DEGREES, RADIANS, CHECKU/F)
+                             (blank line)
+1 -168.0 -1.40351172223 0.0423346818422
+2 -156.0 -1.70447981034 0.00811786522531
+3 -144.0 -1.62956100432 -0.0184129719987
+...
+30 180.0 -0.707106781187 0.0719306095245 :pre
+
+# Example 2: table of the potential. Forces omitted :pre
+
+DIH_TABLE2
+N 30 NOF CHECKU testU.dat CHECKF testF.dat :pre
+
+1 -168.0 -1.40351172223
+2 -156.0 -1.70447981034
+3 -144.0 -1.62956100432
+...
+30 180.0 -0.707106781187 :pre
+
+A section begins with a non-blank line whose 1st character is not a
+"#"; blank lines or lines starting with "#" can be used as comments
+between sections. The first line begins with a keyword which
+identifies the section. The line can contain additional text, but the
+initial text must match the argument specified in the
+"dihedral_coeff"_dihedral_coeff.html command. The next line lists (in
+any order) one or more parameters for the table. Each parameter is a
+keyword followed by one or more numeric values.
+
+Following a blank line, the next N lines list the tabulated values. On
+each line, the 1st value is the index from 1 to N, the 2nd value is
+the angle value, the 3rd value is the energy (in energy units), and
+the 4th is -dE/d(phi) also in energy units). The 3rd term is the
+energy of the 4-atom configuration for the specified angle.  The 4th
+term (when present) is the negative derivative of the energy with
+respect to the angle (in degrees, or radians depending on whether the
+user selected DEGREES or RADIANS).  Thus the units of the last term
+are still energy, not force. The dihedral angle values must increase
+from one line to the next.
+
+Dihedral table splines are cyclic.  There is no discontinuity at 180
+degrees (or at any other angle).  Although in the examples above, the
+angles range from -180 to 180 degrees, in general, the first angle in
+the list can have any value (positive, zero, or negative).  However
+the {range} of angles represented in the table must be {strictly} less
+than 360 degrees (2pi radians) to avoid angle overlap.  (You may not
+supply entries in the table for both 180 and -180, for example.)  If
+the user's table covers only a narrow range of dihedral angles,
+strange numerical behavior can occur in the large remaining gap.
+
+[Parameters:]
+
+The parameter "N" is required and its value is the number of table
+entries that follow. Note that this may be different than the N
+specified in the "dihedral_style table"_dihedral_style.html command.
+Let {Ntable} is the number of table entries requested dihedral_style
+command, and let {Nfile} be the parameter following "N" in the
+tabulated file ("30" in the sparse example above).  What LAMMPS does
+is a preliminary interpolation by creating splines using the {Nfile}
+tabulated values as nodal points.  It uses these to interpolate as
+needed to generate energy and derivative values at {Ntable} different
+points (which are evenly spaced over a 360 degree range, even if the
+angles in the file are not).  The resulting tables of length {Ntable}
+are then used as described above, when computing energy and force for
+individual dihedral angles and their atoms.  This means that if you
+want the interpolation tables of length {Ntable} to match exactly what
+is in the tabulated file (with effectively nopreliminary
+interpolation), you should set {Ntable} = {Nfile}.  To insure the
+nodal points in the user's file are aligned with the interpolated
+table entries, the angles in the table should be integer multiples of
+360/{Ntable} degrees, or 2*PI/{Ntable} radians (depending on your
+choice of angle units).
+
+The optional "NOF" keyword allows the user to omit the forces
+(negative energy derivatives) from the table file (normally located in
+the 4th column).  In their place, forces will be calculated
+automatically by differentiating the potential energy function
+indicated by the 3rd column of the table (using either linear or
+spline interpolation).
+
+The optional "DEGREES" keyword allows the user to specify angles in
+degrees instead of radians (default).
+
+The optional "RADIANS" keyword allows the user to specify angles in
+radians instead of degrees.  (Note: This changes the way the forces
+are scaled in the 4th column of the data file.)
+
+The optional "CHECKU" keyword is followed by a filename.  This allows
+the user to save all of the the {Ntable} different entries in the
+interpolated energy table to a file to make sure that the interpolated
+function agrees with the user's expectations.  (Note: You can
+temporarily increase the {Ntable} parameter to a high value for this
+purpose.  "{Ntable}" is explained above.)
+
+The optional "CHECKF" keyword is analogous to the "CHECKU" keyword.
+It is followed by a filename, and it allows the user to check the
+interpolated force table.  This option is available even if the user
+selected the "NOF" option.
+
+Note that one file can contain many sections, each with a tabulated
+potential. LAMMPS reads the file section by section until it finds one
+that matches the specified keyword.
+
+:line
+
+Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
+functionally the same as the corresponding style without the suffix.
+They have been optimized to run faster, depending on your available
+hardware, as discussed in "Section_accelerate"_Section_accelerate.html
+of the manual.  The accelerated styles take the same arguments and
+should produce the same results, except for round-off and precision
+issues.
+
+These accelerated styles are part of the GPU, USER-INTEL, KOKKOS,
+USER-OMP and OPT packages, respectively.  They are only enabled if
+LAMMPS was built with those packages.  See the "Making
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+You can specify the accelerated styles explicitly in your input script
+by including their suffix, or you can use the "-suffix command-line
+switch"_Section_start.html#start_6 when you invoke LAMMPS, or you can
+use the "suffix"_suffix.html command in your input script.
+
+See "Section_accelerate"_Section_accelerate.html of the manual for
+more instructions on how to use the accelerated styles effectively.
+
+[Restrictions:]
+
+This dihedral style can only be used if LAMMPS was built with the
+USER-MISC package.  See the "Making LAMMPS"_Section_start.html#2_3
+section for more info on packages.
+
+[Related commands:]
+
+"dihedral_coeff"_dihedral_coeff.html
+
+[Default:] none
+:line
+
+:link(dihedralcut-Salerno)
+[(Salerno)] Salerno, Bernstein, J Chem Theory Comput, --, ---- (2018).

src/USER-MISC/README

@@ -37,6 +37,7 @@ dihedral_style nharmonic, Loukas Peristeras, loukas.peristeras at scienomics.com
 dihedral_style quadratic, Loukas Peristeras, loukas.peristeras at scienomics.com, 27 Oct 12
 dihedral_style spherical, Andrew Jewett, jewett.aij@gmail.com, 15 Jul 16
 dihedral_style table, Andrew Jewett, jewett.aij@gmail.com, 10 Jan 12
+dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18
 fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
 fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
 fix bond/react, Jacob Gissinger (CU Boulder), info at disarmmd.org, 24 Feb 2018

src/USER-MISC/dihedral_table_cut.h

@@ -0,0 +1,173 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef DIHEDRAL_CLASS
+
+DihedralStyle(table/cut,DihedralTableCut)
+
+#else
+
+#ifndef LMP_DIHEDRAL_TABLE_CUT_H
+#define LMP_DIHEDRAL_TABLE_CUT_H
+
+#include <stdio.h>
+#include "dihedral.h"
+
+namespace LAMMPS_NS {
+
+class DihedralTableCut : public Dihedral {
+ public:
+  DihedralTableCut(class LAMMPS *);
+  virtual ~DihedralTableCut();
+  virtual void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void write_restart(FILE *);
+  void read_restart(FILE *);
+  void write_data(FILE *);
+  double single(int type, int i1, int i2, int i3, int i4);
+
+ protected:
+  double *k1,*k2,*k3;
+  double *phi1,*phi2,*phi3;
+  double *aat_k,*aat_theta0_1,*aat_theta0_2;
+  int *setflag_d;
+  int *setflag_aat;
+
+  void allocate();
+
+  int tabstyle,tablength;
+  // double *phi0;       <- equilibrium angles not supported
+  char *checkU_fname;
+  char *checkF_fname;
+
+  struct Table {
+    int ninput;
+    //double phi0;      <-equilibrium angles not supported
+    int f_unspecified; // boolean (but MPI does not like type "bool")
+    int use_degrees;   // boolean (but MPI does not like type "bool")
+    double *phifile,*efile,*ffile;
+    double *e2file,*f2file;
+    double delta,invdelta,deltasq6;
+    double *phi,*e,*de,*f,*df,*e2,*f2;
+  };
+
+  int ntables;
+  Table *tables;
+  int *tabindex;
+
+  void null_table(Table *);
+  void free_table(Table *);
+  void read_table(Table *, char *, char *);
+  void bcast_table(Table *);
+  void spline_table(Table *);
+  void compute_table(Table *);
+
+  void param_extract(Table *, char *);
+
+  // --------------------------------------------
+  // ------------ inline functions --------------
+  // --------------------------------------------
+
+  // -----------------------------------------------------------
+  //   uf_lookup()
+  //   quickly calculate the potential u and force f at angle x,
+  //   using the internal tables tb->e and tb->f (evenly spaced)
+  // -----------------------------------------------------------
+  enum{LINEAR,SPLINE};
+
+  inline void uf_lookup(int type, double x, double &u, double &f)
+  {
+    Table *tb = &tables[tabindex[type]];
+    double x_over_delta = x*tb->invdelta;
+    int i = static_cast<int> (x_over_delta);
+    double a;
+    double b = x_over_delta - i;
+    // Apply periodic boundary conditions to indices i and i+1
+    if (i >= tablength) i -= tablength;
+    int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
+
+    switch(tabstyle) {
+      case LINEAR:
+        u = tb->e[i] + b * tb->de[i];
+        f = -(tb->f[i] + b * tb->df[i]); //<--works even if tb->f_unspecified==true
+        break;
+      case SPLINE:
+        a = 1.0 - b;
+        u = a * tb->e[i] + b * tb->e[ip1] +
+          ((a*a*a-a)*tb->e2[i] + (b*b*b-b)*tb->e2[ip1]) *
+          tb->deltasq6;
+        if (tb->f_unspecified)
+          //Formula below taken from equation3.3.5 of "numerical recipes in c"
+          //"f"=-derivative of e with respect to x (or "phi" in this case)
+          f = -((tb->e[i]-tb->e[ip1])*tb->invdelta +
+            ((3.0*a*a-1.0)*tb->e2[i]+(1.0-3.0*b*b)*tb->e2[ip1])*tb->delta/6.0);
+        else
+          f = -(a * tb->f[i] + b * tb->f[ip1] +
+            ((a*a*a-a)*tb->f2[i] + (b*b*b-b)*tb->f2[ip1]) *
+            tb->deltasq6);
+        break;
+    } // switch(tabstyle)
+  } // uf_lookup()
+
+
+  // ----------------------------------------------------------
+  //    u_lookup()
+  //  quickly calculate the potential u at angle x using tb->e
+  //-----------------------------------------------------------
+
+  inline void u_lookup(int type, double x, double &u)
+  {
+    Table *tb = &tables[tabindex[type]];
+    int N = tablength;
+
+    //  i = static_cast<int> ((x - tb->lo) * tb->invdelta); <-general version
+    double x_over_delta = x*tb->invdelta;
+    int    i = static_cast<int> (x_over_delta);
+    double b = x_over_delta - i;
+
+    // Apply periodic boundary conditions to indices i and i+1
+    if (i >= N) i -= N;
+    int ip1 = i+1; if (ip1 >= N) ip1 -= N;
+
+    if (tabstyle == LINEAR) {
+      u = tb->e[i] + b * tb->de[i];
+    }
+    else if (tabstyle == SPLINE) {
+      double a = 1.0 - b;
+      u = a * tb->e[i] + b * tb->e[ip1] +
+        ((a*a*a-a)*tb->e2[i] + (b*b*b-b)*tb->e2[ip1]) *
+        tb->deltasq6;
+    }
+  } // u_lookup()
+
+
+};
+
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+W: Dihedral problem: %d %ld %d %d %d %d
+
+Conformation of the 4 listed dihedral atoms is extreme; you may want
+to check your simulation geometry.
+
+E: Incorrect args for dihedral coefficients
+
+Self-explanatory.  Check the input script or data file.
+
+*/