diff --git a/doc/src/Eqs/dihedral_table_cut.jpg b/doc/src/Eqs/dihedral_table_cut.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..c124184c1773c60fe2564d8cabce2f6c7d7883c0
Binary files /dev/null and b/doc/src/Eqs/dihedral_table_cut.jpg differ
diff --git a/doc/src/Eqs/dihedral_table_cut.tex b/doc/src/Eqs/dihedral_table_cut.tex
new file mode 100644
index 0000000000000000000000000000000000000000..3cc1d331a8978d76e668de3f83cf0d4264be3301
--- /dev/null
+++ b/doc/src/Eqs/dihedral_table_cut.tex
@@ -0,0 +1,11 @@
+\documentclass[12pt]{article}
+\pagestyle{empty}
+\begin{document}
+
+\begin{eqnarray*}
+ f(\theta) & = & K \qquad\qquad\qquad\qquad\qquad\qquad \theta < \theta_1 \\
+ f(\theta) & = & K \left(1-\frac{(\theta - \theta_1)^2}{(\theta_2 - \theta_1)^2}\right) \qquad \theta_1 < \theta < \theta_2
+\end{eqnarray*}
+
+\end{document}
+
diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index cdbdec9ee3ef6e7a8d13cbb28b53d9d41034c487..3dabdbeaa1e1091db35b798f19cf9f3b2b74848f 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -1212,7 +1212,8 @@ package"_Section_start.html#start_3.
"nharmonic (o)"_dihedral_nharmonic.html,
"quadratic (o)"_dihedral_quadratic.html,
"spherical (o)"_dihedral_spherical.html,
-"table (o)"_dihedral_table.html :tb(c=4,ea=c)
+"table (o)"_dihedral_table.html,
+"table/cut"_dihedral_table_cut.html :tb(c=4,ea=c)
:line
diff --git a/doc/src/dihedral_table_cut.txt b/doc/src/dihedral_table_cut.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1c83d4ffa09f7221ac8c3859594abc754b9b28c4
--- /dev/null
+++ b/doc/src/dihedral_table_cut.txt
@@ -0,0 +1,205 @@
+"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
+
+[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.
+
+[Restrictions:]
+
+This dihedral style can only be used if LAMMPS was built with the
+USER-MISC package. See the "Making LAMMPS"_Section_start.html#start_3
+section for more info on packages.
+
+[Related commands:]
+
+"dihedral_coeff"_dihedral_coeff.html, "dihedral_style table"_dihedral_table.html
+
+[Default:] none
+
+:link(dihedralcut-Salerno)
+[(Salerno)] Salerno, Bernstein, J Chem Theory Comput, --, ---- (2018).
diff --git a/doc/src/dihedrals.txt b/doc/src/dihedrals.txt
index 500a6a52bf69b7beb641c414f6dee354d6ca6b7b..a862bf50a0a9c357ffeed1129e3b67423259447c 100644
--- a/doc/src/dihedrals.txt
+++ b/doc/src/dihedrals.txt
@@ -19,6 +19,7 @@ Dihedral Styles :h1
dihedral_quadratic
dihedral_spherical
dihedral_table
+ dihedral_table_cut
dihedral_zero
dihedral_charmm
dihedral_class2
diff --git a/doc/src/lammps.book b/doc/src/lammps.book
index 1f4bf157f4829df2446970b8c9550f7098b765f8..0764c593f7e56cc251588d800a53d06a11e232b4 100644
--- a/doc/src/lammps.book
+++ b/doc/src/lammps.book
@@ -582,6 +582,7 @@ dihedral_opls.html
dihedral_quadratic.html
dihedral_spherical.html
dihedral_table.html
+dihedral_table_cut.html
dihedral_zero.html
lammps_commands_improper.html
diff --git a/src/.gitignore b/src/.gitignore
index 3e06a33580122b07bf88f6e5a0fb8bf992752073..3ae05079d211d2c043615cf93be2013449090506 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -324,6 +324,8 @@
/dihedral_spherical.h
/dihedral_table.cpp
/dihedral_table.h
+/dihedral_table_cut.cpp
+/dihedral_table_cut.h
/dump_atom_gz.cpp
/dump_atom_gz.h
/dump_xyz_gz.cpp
diff --git a/src/USER-MISC/README b/src/USER-MISC/README
index 7fc931b9625ed9fa6884bb5b1564692cf77564dd..222c06b4c426ff1adf271dcccf7438921eb5e5ff 100644
--- a/src/USER-MISC/README
+++ b/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
diff --git a/src/USER-MISC/dihedral_table_cut.cpp b/src/USER-MISC/dihedral_table_cut.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7345bcdad164c2df0cd598c8606b8f26cf1fc511
--- /dev/null
+++ b/src/USER-MISC/dihedral_table_cut.cpp
@@ -0,0 +1,1500 @@
+/* ----------------------------------------------------------------------
+ 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.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: K. Michael Salerno (NRL)
+ Based on tabulated dihedral (dihedral_table.cpp) by Andrew Jewett
+------------------------------------------------------------------------- */
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include <cassert>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include "dihedral_table_cut.h"
+#include "atom.h"
+#include "neighbor.h"
+#include "update.h"
+#include "domain.h"
+#include "comm.h"
+#include "force.h"
+#include "citeme.h"
+#include "math_const.h"
+#include "math_extra.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+using namespace MathConst;
+using namespace std;
+using namespace MathExtra;
+
+
+static const char cite_dihedral_tablecut[] =
+ "dihedral_style tablecut command:\n\n"
+ "@Article{Salerno17,\n"
+ " author = {K. M. Salerno and N. Bernstein},\n"
+ " title = {Persistence Length, End-to-End Distance, and Structure of Coarse-Grained Polymers},\n"
+ " journal = {J.~Chem.~Theory Comput.},\n"
+ " year = 2018,\n"
+ " DOI = 10.1021/acs.jctc.7b01229"
+ "}\n\n";
+
+/* ---------------------------------------------------------------------- */
+
+#define TOLERANCE 0.05
+#define SMALL 0.0000001
+
+// ------------------------------------------------------------------------
+// The following auxiliary functions were left out of the
+// DihedralTable class either because they require template parameters,
+// or because they have nothing to do with dihedral angles.
+// ------------------------------------------------------------------------
+
+// -------------------------------------------------------------------
+// --------- The function was stolen verbatim from the ---------
+// --------- GNU Scientific Library (GSL, version 1.15) ---------
+// -------------------------------------------------------------------
+
+/* Author: Gerard Jungman */
+/* for description of method see [Engeln-Mullges + Uhlig, p. 96]
+ *
+ * diag[0] offdiag[0] 0 ..... offdiag[N-1]
+ * offdiag[0] diag[1] offdiag[1] .....
+ * 0 offdiag[1] diag[2]
+ * 0 0 offdiag[2] .....
+ * ... ...
+ * offdiag[N-1] ...
+ *
+ */
+// -- (A non-symmetric version of this function is also available.) --
+
+enum { //GSL status return codes.
+ GSL_FAILURE = -1,
+ GSL_SUCCESS = 0,
+ GSL_ENOMEM = 8,
+ GSL_EZERODIV = 12,
+ GSL_EBADLEN = 19
+};
+
+
+static int solve_cyc_tridiag( const double diag[], size_t d_stride,
+ const double offdiag[], size_t o_stride,
+ const double b[], size_t b_stride,
+ double x[], size_t x_stride,
+ size_t N, bool warn)
+{
+ int status = GSL_SUCCESS;
+ double * delta = (double *) malloc (N * sizeof (double));
+ double * gamma = (double *) malloc (N * sizeof (double));
+ double * alpha = (double *) malloc (N * sizeof (double));
+ double * c = (double *) malloc (N * sizeof (double));
+ double * z = (double *) malloc (N * sizeof (double));
+
+ if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0) {
+ if (warn)
+ fprintf(stderr,"Internal Cyclic Spline Error: failed to allocate working space\n");
+
+ if (delta) free(delta);
+ if (gamma) free(gamma);
+ if (alpha) free(alpha);
+ if (c) free(c);
+ if (z) free(z);
+ return GSL_ENOMEM;
+ }
+ else
+ {
+ size_t i, j;
+ double sum = 0.0;
+
+ /* factor */
+
+ if (N == 1)
+ {
+ x[0] = b[0] / diag[0];
+ free(delta);
+ free(gamma);
+ free(alpha);
+ free(c);
+ free(z);
+ return GSL_SUCCESS;
+ }
+
+ alpha[0] = diag[0];
+ gamma[0] = offdiag[0] / alpha[0];
+ delta[0] = offdiag[o_stride * (N-1)] / alpha[0];
+
+ if (alpha[0] == 0) {
+ status = GSL_EZERODIV;
+ }
+
+ for (i = 1; i < N - 2; i++)
+ {
+ alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
+ gamma[i] = offdiag[o_stride * i] / alpha[i];
+ delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
+ if (alpha[i] == 0) {
+ status = GSL_EZERODIV;
+ }
+ }
+
+ for (i = 0; i < N - 2; i++)
+ {
+ sum += alpha[i] * delta[i] * delta[i];
+ }
+
+ alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
+
+ gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
+
+ alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2];
+
+ /* update */
+ z[0] = b[0];
+ for (i = 1; i < N - 1; i++)
+ {
+ z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
+ }
+ sum = 0.0;
+ for (i = 0; i < N - 2; i++)
+ {
+ sum += delta[i] * z[i];
+ }
+ z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
+ for (i = 0; i < N; i++)
+ {
+ c[i] = z[i] / alpha[i];
+ }
+
+ /* backsubstitution */
+ x[x_stride * (N - 1)] = c[N - 1];
+ x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
+ if (N >= 3)
+ {
+ for (i = N - 3, j = 0; j <= N - 3; j++, i--)
+ {
+ x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
+ }
+ }
+ }
+
+ free (z);
+ free (c);
+ free (alpha);
+ free (gamma);
+ free (delta);
+
+ if ((status == GSL_EZERODIV) && warn)
+ fprintf(stderr, "Internal Cyclic Spline Error: Matrix must be positive definite.\n");
+
+ return status;
+} //solve_cyc_tridiag()
+
+/* ----------------------------------------------------------------------
+ spline and splint routines modified from Numerical Recipes
+------------------------------------------------------------------------- */
+
+static int cyc_spline(double const *xa,
+ double const *ya,
+ int n,
+ double period,
+ double *y2a, bool warn)
+{
+
+ double *diag = new double[n];
+ double *offdiag = new double[n];
+ double *rhs = new double[n];
+ double xa_im1, xa_ip1;
+
+ // In the cyclic case, there are n equations with n unknows.
+ // The for loop sets up the equations we need to solve.
+ // Later we invoke the GSL tridiagonal matrix solver to solve them.
+
+ for(int i=0; i < n; i++) {
+
+ // I have to lookup xa[i+1] and xa[i-1]. This gets tricky because of
+ // periodic boundary conditions. We handle that now.
+ int im1 = i-1;
+ if (im1<0) {
+ im1 += n;
+ xa_im1 = xa[im1] - period;
+ }
+ else
+ xa_im1 = xa[im1];
+
+ int ip1 = i+1;
+ if (ip1>=n) {
+ ip1 -= n;
+ xa_ip1 = xa[ip1] + period;
+ }
+ else
+ xa_ip1 = xa[ip1];
+
+ // Recall that we want to find the y2a[] parameters (there are n of them).
+ // To solve for them, we have a linear equation with n unknowns
+ // (in the cyclic case that is). For details, the non-cyclic case is
+ // explained in equation 3.3.7 in Numerical Recipes in C, p. 115.
+ diag[i] = (xa_ip1 - xa_im1) / 3.0;
+ offdiag[i] = (xa_ip1 - xa[i]) / 6.0;
+ rhs[i] = ((ya[ip1] - ya[i]) / (xa_ip1 - xa[i])) -
+ ((ya[i] - ya[im1]) / (xa[i] - xa_im1));
+ }
+
+ // Because this matrix is tridiagonal (and cyclic), we can use the following
+ // cheap method to invert it.
+ if (solve_cyc_tridiag(diag, 1,
+ offdiag, 1,
+ rhs, 1,
+ y2a, 1,
+ n, warn) != GSL_SUCCESS) {
+ if (warn)
+ fprintf(stderr,"Error in inverting matrix for splines.\n");
+
+ delete [] diag;
+ delete [] offdiag;
+ delete [] rhs;
+ return 1;
+ }
+ delete [] diag;
+ delete [] offdiag;
+ delete [] rhs;
+ return 0;
+} // cyc_spline()
+
+/* ---------------------------------------------------------------------- */
+
+// cyc_splint(): Evaluates a spline at position x, with n control
+// points located at xa[], ya[], and with parameters y2a[]
+// The xa[] must be monotonically increasing and their
+// range should not exceed period (ie xa[n-1] < xa[0] + period).
+// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)]
+// "period" is typically 2*PI.
+static double cyc_splint(double const *xa,
+ double const *ya,
+ double const *y2a,
+ int n,
+ double period,
+ double x)
+{
+ int klo = -1;
+ int khi = n;
+ int k;
+ double xlo = xa[n-1] - period;
+ double xhi = xa[0] + period;
+ while (khi-klo > 1) {
+ k = (khi+klo) >> 1; //(k=(khi+klo)/2)
+ if (xa[k] > x) {
+ khi = k;
+ xhi = xa[k];
+ }
+ else {
+ klo = k;
+ xlo = xa[k];
+ }
+ }
+
+ if (khi == n) khi = 0;
+ if (klo ==-1) klo = n-1;
+
+ double h = xhi-xlo;
+ double a = (xhi-x) / h;
+ double b = (x-xlo) / h;
+ double y = a*ya[klo] + b*ya[khi] +
+ ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+
+ return y;
+
+} // cyc_splint()
+
+
+static double cyc_lin(double const *xa,
+ double const *ya,
+ int n,
+ double period,
+ double x)
+{
+ int klo = -1;
+ int khi = n;
+ int k;
+ double xlo = xa[n-1] - period;
+ double xhi = xa[0] + period;
+ while (khi-klo > 1) {
+ k = (khi+klo) >> 1; //(k=(khi+klo)/2)
+ if (xa[k] > x) {
+ khi = k;
+ xhi = xa[k];
+ }
+ else {
+ klo = k;
+ xlo = xa[k];
+ }
+ }
+
+ if (khi == n) khi = 0;
+ if (klo ==-1) klo = n-1;
+
+ double h = xhi-xlo;
+ double a = (xhi-x) / h;
+ double b = (x-xlo) / h;
+ double y = a*ya[klo] + b*ya[khi];
+
+ return y;
+
+} // cyc_lin()
+
+
+
+
+// cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
+// with n control points at xa[], ya[], with parameters y2a[].
+// The xa[] must be monotonically increasing and their
+// range should not exceed period (ie xa[n-1] < xa[0] + period).
+// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)]
+// "period" is typically 2*PI.
+static double cyc_splintD(double const *xa,
+ double const *ya,
+ double const *y2a,
+ int n,
+ double period,
+ double x)
+{
+ int klo = -1;
+ int khi = n; // (not n-1)
+ int k;
+ double xlo = xa[n-1] - period;
+ double xhi = xa[0] + period;
+ while (khi-klo > 1) {
+ k = (khi+klo) >> 1; //(k=(khi+klo)/2)
+ if (xa[k] > x) {
+ khi = k;
+ xhi = xa[k];
+ }
+ else {
+ klo = k;
+ xlo = xa[k];
+ }
+ }
+
+ if (khi == n) khi = 0;
+ if (klo ==-1) klo = n-1;
+
+ double yhi = ya[khi];
+ double ylo = ya[klo];
+ double h = xhi-xlo;
+ double g = yhi-ylo;
+ double a = (xhi-x) / h;
+ double b = (x-xlo) / h;
+ // Formula below taken from equation 3.3.5 of "numerical recipes in c"
+ // "yD" = the derivative of y
+ double yD = g/h - ( (3.0*a*a-1.0)*y2a[klo] - (3.0*b*b-1.0)*y2a[khi] ) * h/6.0;
+ // For rerefence: y = a*ylo + b*yhi +
+ // ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
+
+ return yD;
+
+} // cyc_splintD()
+
+
+/* ---------------------------------------------------------------------- */
+
+DihedralTableCut::DihedralTableCut(LAMMPS *lmp) : Dihedral(lmp)
+{
+ if (lmp->citeme) lmp->citeme->add(cite_dihedral_tablecut);
+ ntables = 0;
+ tables = NULL;
+ checkU_fname = checkF_fname = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+DihedralTableCut::~DihedralTableCut()
+{
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(setflag_d);
+ memory->destroy(setflag_aat);
+
+ memory->destroy(k1);
+ memory->destroy(k2);
+ memory->destroy(k3);
+ memory->destroy(phi1);
+ memory->destroy(phi2);
+ memory->destroy(phi3);
+
+ memory->destroy(aat_k);
+ memory->destroy(aat_theta0_1);
+ memory->destroy(aat_theta0_2);
+
+ for (int m = 0; m < ntables; m++) free_table(&tables[m]);
+ memory->sfree(tables);
+ memory->sfree(checkU_fname);
+ memory->sfree(checkF_fname);
+
+ memory->destroy(setflag);
+ memory->destroy(tabindex);
+
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralTableCut::compute(int eflag, int vflag)
+{
+
+ int i1,i2,i3,i4,i,j,k,n,type;
+ double edihedral;
+ double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm;
+ double fphi,fpphi;
+ double r1mag2,r1,r2mag2,r2,r3mag2,r3;
+ double sb1,rb1,sb2,rb2,sb3,rb3,c0,r12c1;
+ double r12c2,costh12,costh13,costh23,sc1,sc2,s1,s2,c;
+ double phi,sinphi,a11,a22,a33,a12,a13,a23,sx1,sx2;
+ double sx12,sy1,sy2,sy12,sz1,sz2,sz12;
+ double t1,t2,t3,t4;
+ double da1,da2;
+ double s12,sin2;
+ double dcosphidr[4][3],dphidr[4][3],dthetadr[2][4][3];
+ double fabcd[4][3];
+
+ edihedral = 0.0;
+ if (eflag || vflag) ev_setup(eflag,vflag);
+ else evflag = 0;
+
+ double **x = atom->x;
+ double **f = atom->f;
+ int **dihedrallist = neighbor->dihedrallist;
+ int ndihedrallist = neighbor->ndihedrallist;
+ int nlocal = atom->nlocal;
+ int newton_bond = force->newton_bond;
+
+ for (n = 0; n < ndihedrallist; n++) {
+ i1 = dihedrallist[n][0];
+ i2 = dihedrallist[n][1];
+ i3 = dihedrallist[n][2];
+ i4 = dihedrallist[n][3];
+ type = dihedrallist[n][4];
+
+ // 1st bond
+
+ vb1x = x[i1][0] - x[i2][0];
+ vb1y = x[i1][1] - x[i2][1];
+ vb1z = x[i1][2] - x[i2][2];
+
+ // 2nd bond
+
+ vb2x = x[i3][0] - x[i2][0];
+ vb2y = x[i3][1] - x[i2][1];
+ vb2z = x[i3][2] - x[i2][2];
+
+ vb2xm = -vb2x;
+ vb2ym = -vb2y;
+ vb2zm = -vb2z;
+
+ // 3rd bond
+
+ vb3x = x[i4][0] - x[i3][0];
+ vb3y = x[i4][1] - x[i3][1];
+ vb3z = x[i4][2] - x[i3][2];
+
+ // distances
+
+ r1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z;
+ r1 = sqrt(r1mag2);
+ r2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z;
+ r2 = sqrt(r2mag2);
+ r3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z;
+ r3 = sqrt(r3mag2);
+
+ sb1 = 1.0/r1mag2;
+ rb1 = 1.0/r1;
+ sb2 = 1.0/r2mag2;
+ rb2 = 1.0/r2;
+ sb3 = 1.0/r3mag2;
+ rb3 = 1.0/r3;
+
+ c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3;
+
+ // angles
+
+ r12c1 = rb1*rb2;
+ r12c2 = rb2*rb3;
+ costh12 = (vb1x*vb2x + vb1y*vb2y + vb1z*vb2z) * r12c1;
+ costh13 = c0;
+ costh23 = (vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z) * r12c2;
+
+ // cos and sin of 2 angles and final c
+
+ sin2 = MAX(1.0 - costh12*costh12,0.0);
+ sc1 = sqrt(sin2);
+ if (sc1 < SMALL) sc1 = SMALL;
+ sc1 = 1.0/sc1;
+
+ sin2 = MAX(1.0 - costh23*costh23,0.0);
+ sc2 = sqrt(sin2);
+ if (sc2 < SMALL) sc2 = SMALL;
+ sc2 = 1.0/sc2;
+
+ s1 = sc1 * sc1;
+ s2 = sc2 * sc2;
+ s12 = sc1 * sc2;
+ c = (c0 + costh12*costh23) * s12;
+
+ // error check
+
+ if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
+ int me;
+ MPI_Comm_rank(world,&me);
+ if (screen) {
+ char str[128];
+ sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT " "
+ TAGINT_FORMAT " " TAGINT_FORMAT,
+ me,update->ntimestep,
+ atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
+ error->warning(FLERR,str,0);
+ fprintf(screen," 1st atom: %d %g %g %g\n",
+ me,x[i1][0],x[i1][1],x[i1][2]);
+ fprintf(screen," 2nd atom: %d %g %g %g\n",
+ me,x[i2][0],x[i2][1],x[i2][2]);
+ fprintf(screen," 3rd atom: %d %g %g %g\n",
+ me,x[i3][0],x[i3][1],x[i3][2]);
+ fprintf(screen," 4th atom: %d %g %g %g\n",
+ me,x[i4][0],x[i4][1],x[i4][2]);
+ }
+ }
+
+ if (c > 1.0) c = 1.0;
+ if (c < -1.0) c = -1.0;
+ double phil = acos(c);
+ phi = acos(c);
+
+ sinphi = sqrt(1.0 - c*c);
+ sinphi = MAX(sinphi,SMALL);
+
+ // n123 = vb1 x vb2
+
+ double n123x = vb1y*vb2z - vb1z*vb2y;
+ double n123y = vb1z*vb2x - vb1x*vb2z;
+ double n123z = vb1x*vb2y - vb1y*vb2x;
+ double n123_dot_vb3 = n123x*vb3x + n123y*vb3y + n123z*vb3z;
+ if (n123_dot_vb3 > 0.0) {
+ phil = -phil;
+ phi = -phi;
+ sinphi = -sinphi;
+ }
+
+ a11 = -c*sb1*s1;
+ a22 = sb2 * (2.0*costh13*s12 - c*(s1+s2));
+ a33 = -c*sb3*s2;
+ a12 = r12c1 * (costh12*c*s1 + costh23*s12);
+ a13 = rb1*rb3*s12;
+ a23 = r12c2 * (-costh23*c*s2 - costh12*s12);
+
+ sx1 = a11*vb1x + a12*vb2x + a13*vb3x;
+ sx2 = a12*vb1x + a22*vb2x + a23*vb3x;
+ sx12 = a13*vb1x + a23*vb2x + a33*vb3x;
+ sy1 = a11*vb1y + a12*vb2y + a13*vb3y;
+ sy2 = a12*vb1y + a22*vb2y + a23*vb3y;
+ sy12 = a13*vb1y + a23*vb2y + a33*vb3y;
+ sz1 = a11*vb1z + a12*vb2z + a13*vb3z;
+ sz2 = a12*vb1z + a22*vb2z + a23*vb3z;
+ sz12 = a13*vb1z + a23*vb2z + a33*vb3z;
+
+ // set up d(cos(phi))/d(r) and dphi/dr arrays
+
+ dcosphidr[0][0] = -sx1;
+ dcosphidr[0][1] = -sy1;
+ dcosphidr[0][2] = -sz1;
+ dcosphidr[1][0] = sx2 + sx1;
+ dcosphidr[1][1] = sy2 + sy1;
+ dcosphidr[1][2] = sz2 + sz1;
+ dcosphidr[2][0] = sx12 - sx2;
+ dcosphidr[2][1] = sy12 - sy2;
+ dcosphidr[2][2] = sz12 - sz2;
+ dcosphidr[3][0] = -sx12;
+ dcosphidr[3][1] = -sy12;
+ dcosphidr[3][2] = -sz12;
+
+ for (i = 0; i < 4; i++)
+ for (j = 0; j < 3; j++)
+ dphidr[i][j] = -dcosphidr[i][j] / sinphi;
+
+
+ for (i = 0; i < 4; i++)
+ for (j = 0; j < 3; j++)
+ fabcd[i][j] = 0;
+ edihedral = 0;
+
+
+ // set up d(theta)/d(r) array
+ // dthetadr(i,j,k) = angle i, atom j, coordinate k
+
+ for (i = 0; i < 2; i++)
+ for (j = 0; j < 4; j++)
+ for (k = 0; k < 3; k++)
+ dthetadr[i][j][k] = 0.0;
+
+ t1 = costh12 / r1mag2;
+ t2 = costh23 / r2mag2;
+ t3 = costh12 / r2mag2;
+ t4 = costh23 / r3mag2;
+
+ // angle12
+
+ dthetadr[0][0][0] = sc1 * ((t1 * vb1x) - (vb2x * r12c1));
+ dthetadr[0][0][1] = sc1 * ((t1 * vb1y) - (vb2y * r12c1));
+ dthetadr[0][0][2] = sc1 * ((t1 * vb1z) - (vb2z * r12c1));
+
+ dthetadr[0][1][0] = sc1 * ((-t1 * vb1x) + (vb2x * r12c1) +
+ (-t3 * vb2x) + (vb1x * r12c1));
+ dthetadr[0][1][1] = sc1 * ((-t1 * vb1y) + (vb2y * r12c1) +
+ (-t3 * vb2y) + (vb1y * r12c1));
+ dthetadr[0][1][2] = sc1 * ((-t1 * vb1z) + (vb2z * r12c1) +
+ (-t3 * vb2z) + (vb1z * r12c1));
+
+ dthetadr[0][2][0] = sc1 * ((t3 * vb2x) - (vb1x * r12c1));
+ dthetadr[0][2][1] = sc1 * ((t3 * vb2y) - (vb1y * r12c1));
+ dthetadr[0][2][2] = sc1 * ((t3 * vb2z) - (vb1z * r12c1));
+
+ // angle23
+
+ dthetadr[1][1][0] = sc2 * ((t2 * vb2x) + (vb3x * r12c2));
+ dthetadr[1][1][1] = sc2 * ((t2 * vb2y) + (vb3y * r12c2));
+ dthetadr[1][1][2] = sc2 * ((t2 * vb2z) + (vb3z * r12c2));
+
+ dthetadr[1][2][0] = sc2 * ((-t2 * vb2x) - (vb3x * r12c2) +
+ (t4 * vb3x) + (vb2x * r12c2));
+ dthetadr[1][2][1] = sc2 * ((-t2 * vb2y) - (vb3y * r12c2) +
+ (t4 * vb3y) + (vb2y * r12c2));
+ dthetadr[1][2][2] = sc2 * ((-t2 * vb2z) - (vb3z * r12c2) +
+ (t4 * vb3z) + (vb2z * r12c2));
+
+ dthetadr[1][3][0] = -sc2 * ((t4 * vb3x) + (vb2x * r12c2));
+ dthetadr[1][3][1] = -sc2 * ((t4 * vb3y) + (vb2y * r12c2));
+ dthetadr[1][3][2] = -sc2 * ((t4 * vb3z) + (vb2z * r12c2));
+
+ // angle/angle/torsion cutoff
+
+ da1 = acos(costh12) - aat_theta0_1[type] ;
+ da2 = acos(costh23) - aat_theta0_1[type] ;
+ double dtheta = aat_theta0_2[type]-aat_theta0_1[type];
+
+ fphi = 0.0;
+ fpphi = 0.0;
+ if (phil < 0) phil +=MY_2PI;
+ uf_lookup(type, phil, fphi, fpphi);
+
+ double gt = aat_k[type];
+ double gtt = aat_k[type];
+ double gpt = 0;
+ double gptt = 0;
+
+ if ( acos(costh12) > aat_theta0_1[type]) {
+ gt *= 1-da1*da1/dtheta/dtheta;
+ gpt = -aat_k[type]*2*da1/dtheta/dtheta;
+ }
+
+ if ( acos(costh23) > aat_theta0_1[type]) {
+ gtt *= 1-da2*da2/dtheta/dtheta;
+ gptt = -aat_k[type]*2*da2/dtheta/dtheta;
+ }
+
+ if (eflag) edihedral = gt*gtt*fphi;
+
+ for (i = 0; i < 4; i++)
+ for (j = 0; j < 3; j++)
+ fabcd[i][j] -= - gt*gtt*fpphi*dphidr[i][j]
+ - gt*gptt*fphi*dthetadr[1][i][j] + gpt*gtt*fphi*dthetadr[0][i][j];
+
+ // apply force to each of 4 atoms
+
+ if (newton_bond || i1 < nlocal) {
+ f[i1][0] += fabcd[0][0];
+ f[i1][1] += fabcd[0][1];
+ f[i1][2] += fabcd[0][2];
+ }
+
+ if (newton_bond || i2 < nlocal) {
+ f[i2][0] += fabcd[1][0];
+ f[i2][1] += fabcd[1][1];
+ f[i2][2] += fabcd[1][2];
+ }
+
+ if (newton_bond || i3 < nlocal) {
+ f[i3][0] += fabcd[2][0];
+ f[i3][1] += fabcd[2][1];
+ f[i3][2] += fabcd[2][2];
+ }
+
+ if (newton_bond || i4 < nlocal) {
+ f[i4][0] += fabcd[3][0];
+ f[i4][1] += fabcd[3][1];
+ f[i4][2] += fabcd[3][2];
+ }
+
+ if (evflag)
+ ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,
+ fabcd[0],fabcd[2],fabcd[3],
+ vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z);
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralTableCut::allocate()
+{
+ allocated = 1;
+ int n = atom->ndihedraltypes;
+
+ memory->create(k1,n+1,"dihedral:k1");
+ memory->create(k2,n+1,"dihedral:k2");
+ memory->create(k3,n+1,"dihedral:k3");
+ memory->create(phi1,n+1,"dihedral:phi1");
+ memory->create(phi2,n+1,"dihedral:phi2");
+ memory->create(phi3,n+1,"dihedral:phi3");
+
+ memory->create(aat_k,n+1,"dihedral:aat_k");
+ memory->create(aat_theta0_1,n+1,"dihedral:aat_theta0_1");
+ memory->create(aat_theta0_2,n+1,"dihedral:aat_theta0_2");
+
+ memory->create(setflag,n+1,"dihedral:setflag");
+ memory->create(setflag_d,n+1,"dihedral:setflag_d");
+ memory->create(setflag_aat,n+1,"dihedral:setflag_aat");
+
+ memory->create(tabindex,n+1,"dihedral:tabindex");
+ //memory->create(phi0,n+1,"dihedral:phi0"); <-equilibrium angles not supported
+ memory->create(setflag,n+1,"dihedral:setflag");
+
+ for (int i = 1; i <= n; i++)
+ setflag[i] = setflag_d[i] = setflag_aat[i] = 0;
+}
+
+void DihedralTableCut::settings(int narg, char **arg)
+{
+ if (narg != 2) error->all(FLERR,"Illegal dihedral_style command");
+
+ if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
+ else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
+ else error->all(FLERR,"Unknown table style in dihedral style table_cut");
+
+ tablength = force->inumeric(FLERR,arg[1]);
+ if (tablength < 3)
+ error->all(FLERR,"Illegal number of dihedral table entries");
+ // delete old tables, since cannot just change settings
+
+ for (int m = 0; m < ntables; m++) free_table(&tables[m]);
+ memory->sfree(tables);
+
+ if (allocated) {
+ memory->destroy(setflag);
+ memory->destroy(tabindex);
+ }
+ allocated = 0;
+
+ ntables = 0;
+ tables = NULL;
+}
+
+/* ----------------------------------------------------------------------
+ set coeffs for one or more types
+ arg1 = "aat" -> AngleAngleTorsion coeffs
+ arg1 -> Dihedral coeffs
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::coeff(int narg, char **arg)
+{
+
+ if (narg != 7) error->all(FLERR,"Incorrect args for dihedral coefficients");
+ if (!allocated) allocate();
+ int ilo,ihi;
+ force->bounds(FLERR,arg[0],atom->ndihedraltypes,ilo,ihi);
+
+ int count = 0;
+
+
+ double k_one = force->numeric(FLERR,arg[2]);
+ double theta0_1_one = force->numeric(FLERR,arg[3]);
+ double theta0_2_one = force->numeric(FLERR,arg[4]);
+
+ // convert theta0's from degrees to radians
+
+ for (int i = ilo; i <= ihi; i++) {
+ aat_k[i] = k_one;
+ aat_theta0_1[i] = theta0_1_one/180.0 * MY_PI;
+ aat_theta0_2[i] = theta0_2_one/180.0 * MY_PI;
+ setflag_aat[i] = 1;
+ count++;
+ }
+
+ int me;
+ MPI_Comm_rank(world,&me);
+ tables = (Table *)
+ memory->srealloc(tables,(ntables+1)*sizeof(Table), "dihedral:tables");
+ Table *tb = &tables[ntables];
+ null_table(tb);
+ if (me == 0) read_table(tb,arg[5],arg[6]);
+ bcast_table(tb);
+
+ // --- check the angle data for range errors ---
+ // --- and resolve issues with periodicity ---
+
+ if (tb->ninput < 2) {
+ string err_msg;
+ err_msg = string("Invalid dihedral table length (")
+ + string(arg[5]) + string(").");
+ error->one(FLERR,err_msg.c_str());
+ }
+ else if ((tb->ninput == 2) && (tabstyle == SPLINE)) {
+ string err_msg;
+ err_msg = string("Invalid dihedral spline table length. (Try linear)\n (")
+ + string(arg[5]) + string(").");
+ error->one(FLERR,err_msg.c_str());
+ }
+
+ // check for monotonicity
+ for (int i=0; i < tb->ninput-1; i++) {
+ if (tb->phifile[i] >= tb->phifile[i+1]) {
+ stringstream i_str;
+ i_str << i+1;
+ string err_msg =
+ string("Dihedral table values are not increasing (") +
+ string(arg[5]) + string(", ")+i_str.str()+string("th entry)");
+ if (i==0)
+ err_msg += string("\n(This is probably a mistake with your table format.)\n");
+ error->all(FLERR,err_msg.c_str());
+ }
+ }
+
+ // check the range of angles
+ double philo = tb->phifile[0];
+ double phihi = tb->phifile[tb->ninput-1];
+ if (tb->use_degrees) {
+ if ((phihi - philo) >= 360) {
+ string err_msg;
+ err_msg = string("Dihedral table angle range must be < 360 degrees (")
+ +string(arg[5]) + string(").");
+ error->all(FLERR,err_msg.c_str());
+ }
+ }
+ else {
+ if ((phihi - philo) >= MY_2PI) {
+ string err_msg;
+ err_msg = string("Dihedral table angle range must be < 2*PI radians (")
+ + string(arg[5]) + string(").");
+ error->all(FLERR,err_msg.c_str());
+ }
+ }
+
+ // convert phi from degrees to radians
+ if (tb->use_degrees) {
+ for (int i=0; i < tb->ninput; i++) {
+ tb->phifile[i] *= MY_PI/180.0;
+ // I assume that if angles are in degrees, then the forces (f=dU/dphi)
+ // are specified with "phi" in degrees as well.
+ tb->ffile[i] *= 180.0/MY_PI;
+ }
+ }
+
+ // We want all the phi dihedral angles to lie in the range from 0 to 2*PI.
+ // But I don't want to restrict users to input their data in this range.
+ // We also want the angles to be sorted in increasing order.
+ // This messy code fixes these problems with the user's data:
+ {
+ double *phifile_tmp = new double [tb->ninput]; //temporary arrays
+ double *ffile_tmp = new double [tb->ninput]; //used for sorting
+ double *efile_tmp = new double [tb->ninput];
+
+ // After re-imaging, does the range of angles cross the 0 or 2*PI boundary?
+ // If so, find the discontinuity:
+ int i_discontinuity = tb->ninput;
+ for (int i=0; i < tb->ninput; i++) {
+ double phi = tb->phifile[i];
+ // Add a multiple of 2*PI to phi until it lies in the range [0, 2*PI).
+ phi -= MY_2PI * floor(phi/MY_2PI);
+ phifile_tmp[i] = phi;
+ efile_tmp[i] = tb->efile[i];
+ ffile_tmp[i] = tb->ffile[i];
+ if ((i>0) && (phifile_tmp[i] < phifile_tmp[i-1])) {
+ //There should only be at most one discontinuity, because we have
+ //insured that the data was sorted before imaging, and because the
+ //range of angle values does not exceed 2*PI.
+ i_discontinuity = i;
+ }
+ }
+
+ int I = 0;
+ for (int i = i_discontinuity; i < tb->ninput; i++) {
+ tb->phifile[I] = phifile_tmp[i];
+ tb->efile[I] = efile_tmp[i];
+ tb->ffile[I] = ffile_tmp[i];
+ I++;
+ }
+ for (int i = 0; i < i_discontinuity; i++) {
+ tb->phifile[I] = phifile_tmp[i];
+ tb->efile[I] = efile_tmp[i];
+ tb->ffile[I] = ffile_tmp[i];
+ I++;
+ }
+
+ // clean up temporary storage
+ delete[] phifile_tmp;
+ delete[] ffile_tmp;
+ delete[] efile_tmp;
+ }
+
+ // spline read-in and compute r,e,f vectors within table
+
+ spline_table(tb);
+ compute_table(tb);
+
+ // Optional: allow the user to print out the interpolated spline tables
+
+ if (me == 0) {
+ if (checkU_fname && (strlen(checkU_fname) != 0))
+ {
+ ofstream checkU_file;
+ checkU_file.open(checkU_fname, ios::out);
+ for (int i=0; i < tablength; i++) {
+ double phi = i*MY_2PI/tablength;
+ double u = tb->e[i];
+ if (tb->use_degrees)
+ phi *= 180.0/MY_PI;
+ checkU_file << phi << " " << u << "\n";
+ }
+ checkU_file.close();
+ }
+ if (checkF_fname && (strlen(checkF_fname) != 0))
+ {
+ ofstream checkF_file;
+ checkF_file.open(checkF_fname, ios::out);
+ for (int i=0; i < tablength; i++)
+ {
+ double phi = i*MY_2PI/tablength;
+ double f;
+ if ((tabstyle == SPLINE) && (tb->f_unspecified)) {
+ double dU_dphi =
+ // (If the user did not specify the forces now, AND the user
+ // selected the "spline" option, (as opposed to "linear")
+ // THEN the tb->f array is uninitialized, so there's
+ // no point to print out the contents of the tb->f[] array.
+ // Instead, later on, we will calculate the force using the
+ // -cyc_splintD() routine to calculate the derivative of the
+ // energy spline, using the energy data (tb->e[]).
+ // To be nice and report something, I do the same thing here.)
+ cyc_splintD(tb->phi, tb->e, tb->e2, tablength, MY_2PI,phi);
+ f = -dU_dphi;
+ }
+ else
+ // Otherwise we calculated the tb->f[] array. Report its contents.
+ f = tb->f[i];
+ if (tb->use_degrees) {
+ phi *= 180.0/MY_PI;
+ // If the user wants degree angle units, we should convert our
+ // internal force tables (in energy/radians) to (energy/degrees)
+ f *= MY_PI/180.0;
+ }
+ checkF_file << phi << " " << f << "\n";
+ }
+ checkF_file.close();
+ } // if (checkF_fname && (strlen(checkF_fname) != 0))
+ } // if (me == 0)
+
+ // store ptr to table in tabindex
+ count = 0;
+ for (int i = ilo; i <= ihi; i++)
+ {
+ tabindex[i] = ntables;
+ //phi0[i] = tb->phi0; <- equilibrium dihedral angles not supported
+ setflag[i] = 1;
+ count++;
+ }
+ ntables++;
+
+
+ if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
+
+ for (int i = ilo; i <= ihi; i++)
+ if (setflag_d[i] == 1 && setflag_aat[i] == 1 )
+ setflag[i] = 1;
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 writes out coeffs to restart file
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::write_restart(FILE *fp)
+{
+ fwrite(&tabstyle,sizeof(int),1,fp);
+ fwrite(&tablength,sizeof(int),1,fp);
+ fwrite(&k1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&k2[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&k3[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&phi1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&phi2[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&phi3[1],sizeof(double),atom->ndihedraltypes,fp);
+
+ fwrite(&aat_k[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fwrite(&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);
+}
+
+/* ----------------------------------------------------------------------
+ proc 0 reads coeffs from restart file, bcasts them
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::read_restart(FILE *fp)
+{
+ allocate();
+
+ if (comm->me == 0) {
+ fread(&tabstyle,sizeof(int),1,fp);
+ fread(&tablength,sizeof(int),1,fp);
+ fread(&k1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&k2[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&k3[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&phi1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&phi2[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&phi3[1],sizeof(double),atom->ndihedraltypes,fp);
+
+ fread(&aat_k[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&aat_theta0_1[1],sizeof(double),atom->ndihedraltypes,fp);
+ fread(&aat_theta0_2[1],sizeof(double),atom->ndihedraltypes,fp);
+ }
+
+ MPI_Bcast(&k1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&k2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&k3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&phi1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&phi2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&phi3[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+
+ MPI_Bcast(&aat_k[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&aat_theta0_1[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&aat_theta0_2[1],atom->ndihedraltypes,MPI_DOUBLE,0,world);
+ MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
+ MPI_Bcast(&tablength,1,MPI_INT,0,world);
+
+ allocate();
+ for (int i = 1; i <= atom->ndihedraltypes; i++) setflag[i] = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralTableCut::null_table(Table *tb)
+{
+ tb->phifile = tb->efile = tb->ffile = NULL;
+ tb->e2file = tb->f2file = NULL;
+ tb->phi = tb->e = tb->de = NULL;
+ tb->f = tb->df = tb->e2 = tb->f2 = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void DihedralTableCut::free_table(Table *tb)
+{
+ memory->destroy(tb->phifile);
+ memory->destroy(tb->efile);
+ memory->destroy(tb->ffile);
+ memory->destroy(tb->e2file);
+ memory->destroy(tb->f2file);
+
+ memory->destroy(tb->phi);
+ memory->destroy(tb->e);
+ memory->destroy(tb->de);
+ memory->destroy(tb->f);
+ memory->destroy(tb->df);
+ memory->destroy(tb->e2);
+ memory->destroy(tb->f2);
+}
+
+/* ----------------------------------------------------------------------
+ read table file, only called by proc 0
+------------------------------------------------------------------------- */
+
+static const int MAXLINE=2048;
+
+void DihedralTableCut::read_table(Table *tb, char *file, char *keyword)
+{
+ char line[MAXLINE];
+
+ // open file
+
+ FILE *fp = force->open_potential(file);
+ if (fp == NULL) {
+ string err_msg = string("Cannot open file ") + string(file);
+ error->one(FLERR,err_msg.c_str());
+ }
+
+ // loop until section found with matching keyword
+
+ while (1) {
+ if (fgets(line,MAXLINE,fp) == NULL) {
+ string err_msg=string("Did not find keyword \"")
+ +string(keyword)+string("\" in dihedral table file.");
+ error->one(FLERR, err_msg.c_str());
+ }
+ if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
+ if (line[0] == '#') continue; // comment
+ char *word = strtok(line," \t\n\r");
+ if (strcmp(word,keyword) == 0) break; // matching keyword
+ fgets(line,MAXLINE,fp); // no match, skip section
+ param_extract(tb,line);
+ fgets(line,MAXLINE,fp);
+ for (int i = 0; i < tb->ninput; i++)
+ fgets(line,MAXLINE,fp);
+ }
+
+ // read args on 2nd line of section
+ // allocate table arrays for file values
+
+ fgets(line,MAXLINE,fp);
+ param_extract(tb,line);
+ memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
+ memory->create(tb->efile,tb->ninput,"dihedral:efile");
+ memory->create(tb->ffile,tb->ninput,"dihedral:ffile");
+
+ // read a,e,f table values from file
+
+ int itmp;
+ for (int i = 0; i < tb->ninput; i++) {
+ // Read the next line. Make sure the file is long enough.
+ if (! fgets(line,MAXLINE,fp))
+ error->one(FLERR, "Dihedral table does not contain enough entries.");
+ // Skip blank lines and delete text following a '#' character
+ char *pe = strchr(line, '#');
+ if (pe != NULL) *pe = '\0'; //terminate string at '#' character
+ char *pc = line;
+ while ((*pc != '\0') && isspace(*pc))
+ pc++;
+ if (*pc != '\0') { //If line is not a blank line
+ stringstream line_ss(line);
+ if (tb->f_unspecified) {
+ //sscanf(line,"%d %lg %lg",
+ // &itmp,&tb->phifile[i],&tb->efile[i]);
+ line_ss >> itmp;
+ line_ss >> tb->phifile[i];
+ line_ss >> tb->efile[i];
+ }
+ else {
+ //sscanf(line,"%d %lg %lg %lg",
+ // &itmp,&tb->phifile[i],&tb->efile[i],&tb->ffile[i]);
+ line_ss >> itmp;
+ line_ss >> tb->phifile[i];
+ line_ss >> tb->efile[i];
+ line_ss >> tb->ffile[i];
+ }
+ if (! line_ss) {
+ stringstream err_msg;
+ err_msg << "Read error in table "<< keyword<<", near line "<<i+1<<"\n"
+ << " (Check to make sure the number of colums is correct.)";
+ if ((! tb->f_unspecified) && (i==0))
+ err_msg << "\n (This sometimes occurs if users forget to specify the \"NOF\" option.)\n";
+ error->one(FLERR, err_msg.str().c_str());
+ }
+ }
+ else //if it is a blank line, then skip it.
+ i--;
+ } //for (int i = 0; (i < tb->ninput) && fp; i++) {
+
+ fclose(fp);
+}
+
+/* ----------------------------------------------------------------------
+ build spline representation of e,f over entire range of read-in table
+ this function sets these values in e2file,f2file.
+ I also perform a crude check for force & energy consistency.
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::spline_table(Table *tb)
+{
+ memory->create(tb->e2file,tb->ninput,"dihedral:e2file");
+ memory->create(tb->f2file,tb->ninput,"dihedral:f2file");
+
+ if (cyc_spline(tb->phifile, tb->efile, tb->ninput,
+ MY_2PI,tb->e2file,comm->me == 0))
+ error->one(FLERR,"Error computing dihedral spline tables");
+
+ if (! tb->f_unspecified) {
+ if (cyc_spline(tb->phifile, tb->ffile, tb->ninput,
+ MY_2PI, tb->f2file, comm->me == 0))
+ error->one(FLERR,"Error computing dihedral spline tables");
+ }
+
+ // CHECK to help make sure the user calculated forces in a way
+ // which is grossly numerically consistent with the energy table.
+ if (! tb->f_unspecified) {
+ int num_disagreements = 0;
+ for (int i=0; i<tb->ninput; i++) {
+
+ // Calculate what the force should be at the control points
+ // by using linear interpolation of the derivatives of the energy:
+
+ double phi_i = tb->phifile[i];
+
+ // First deal with periodicity
+ double phi_im1, phi_ip1;
+ int im1 = i-1;
+ if (im1 < 0) {
+ im1 += tb->ninput;
+ phi_im1 = tb->phifile[im1] - MY_2PI;
+ }
+ else
+ phi_im1 = tb->phifile[im1];
+ int ip1 = i+1;
+ if (ip1 >= tb->ninput) {
+ ip1 -= tb->ninput;
+ phi_ip1 = tb->phifile[ip1] + MY_2PI;
+ }
+ else
+ phi_ip1 = tb->phifile[ip1];
+
+ // Now calculate the midpoints above and below phi_i = tb->phifile[i]
+ double phi_lo= 0.5*(phi_im1 + phi_i); //midpoint between phi_im1 and phi_i
+ double phi_hi= 0.5*(phi_i + phi_ip1); //midpoint between phi_i and phi_ip1
+
+ // Use a linear approximation to the derivative at these two midpoints
+ double dU_dphi_lo =
+ (tb->efile[i] - tb->efile[im1])
+ /
+ (phi_i - phi_im1);
+ double dU_dphi_hi =
+ (tb->efile[ip1] - tb->efile[i])
+ /
+ (phi_ip1 - phi_i);
+
+ // Now calculate the derivative at position
+ // phi_i (=tb->phifile[i]) using linear interpolation
+
+ double a = (phi_i - phi_lo) / (phi_hi - phi_lo);
+ double b = (phi_hi - phi_i) / (phi_hi - phi_lo);
+ double dU_dphi = a*dU_dphi_lo + b*dU_dphi_hi;
+ double f = -dU_dphi;
+ // Alternately, we could use spline interpolation instead:
+ // double f = - splintD(tb->phifile, tb->efile, tb->e2file,
+ // tb->ninput, MY_2PI, tb->phifile[i]);
+ // This is the way I originally did it, but I trust
+ // the ugly simple linear way above better.
+ // Recall this entire block of code doess not calculate
+ // anything important. It does not have to be perfect.
+ // We are only checking for stupid user errors here.
+
+ if ((f != 0.0) &&
+ (tb->ffile[i] != 0.0) &&
+ ((f/tb->ffile[i] < 0.5) || (f/tb->ffile[i] > 2.0))) {
+ num_disagreements++;
+ }
+ } // for (int i=0; i<tb->ninput; i++)
+
+ if ((num_disagreements > tb->ninput/2) && (num_disagreements > 2)) {
+ string msg("Dihedral table has inconsistent forces and energies. (Try \"NOF\".)\n");
+ error->all(FLERR,msg.c_str());
+ }
+
+ } // check for consistency if (! tb->f_unspecified)
+
+} // DihedralTable::spline_table()
+
+
+/* ----------------------------------------------------------------------
+ compute a,e,f vectors from splined values
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::compute_table(Table *tb)
+{
+ //delta = table spacing in dihedral angle for tablength (cyclic) bins
+ tb->delta = MY_2PI / tablength;
+ tb->invdelta = 1.0/tb->delta;
+ tb->deltasq6 = tb->delta*tb->delta / 6.0;
+
+ // N evenly spaced bins in dihedral angle from 0 to 2*PI
+ // phi,e,f = value at lower edge of bin
+ // de,df values = delta values of e,f (cyclic, in this case)
+ // phi,e,f,de,df are arrays containing "tablength" number of entries
+
+ memory->create(tb->phi,tablength,"dihedral:phi");
+ memory->create(tb->e,tablength,"dihedral:e");
+ memory->create(tb->de,tablength,"dihedral:de");
+ memory->create(tb->f,tablength,"dihedral:f");
+ memory->create(tb->df,tablength,"dihedral:df");
+ memory->create(tb->e2,tablength,"dihedral:e2");
+ memory->create(tb->f2,tablength,"dihedral:f2");
+
+ if (tabstyle == SPLINE) {
+ // Use cubic spline interpolation to calculate the entries in the
+ // internal table. (This is true regardless...even if tabstyle!=SPLINE.)
+ for (int i = 0; i < tablength; i++) {
+ double phi = i*tb->delta;
+ tb->phi[i] = phi;
+ tb->e[i]= cyc_splint(tb->phifile,tb->efile,tb->e2file,tb->ninput,MY_2PI,phi);
+ if (! tb->f_unspecified)
+ tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi);
+ }
+ } // if (tabstyle == SPLINE)
+ else if (tabstyle == LINEAR) {
+ if (! tb->f_unspecified) {
+ for (int i = 0; i < tablength; i++) {
+ double phi = i*tb->delta;
+ tb->phi[i] = phi;
+ tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
+ tb->f[i]= cyc_lin(tb->phifile,tb->ffile,tb->ninput,MY_2PI,phi);
+ }
+ }
+ else {
+ for (int i = 0; i < tablength; i++) {
+ double phi = i*tb->delta;
+ tb->phi[i] = phi;
+ tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
+ }
+ // In the linear case, if the user did not specify the forces, then we
+ // must generate the "f" array. Do this using linear interpolation
+ // of the e array (which itself was generated above)
+ for (int i = 0; i < tablength; i++) {
+ int im1 = i-1; if (im1 < 0) im1 += tablength;
+ int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
+ double dedx = (tb->e[ip1] - tb->e[im1]) / (2.0 * tb->delta);
+ // (This is the average of the linear slopes on either side of the node.
+ // Note that the nodes in the internal table are evenly spaced.)
+ tb->f[i] = -dedx;
+ }
+ }
+
+
+ // Fill the linear interpolation tables (de, df)
+ for (int i = 0; i < tablength; i++) {
+ int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
+ tb->de[i] = tb->e[ip1] - tb->e[i];
+ tb->df[i] = tb->f[ip1] - tb->f[i];
+ }
+ } // else if (tabstyle == LINEAR)
+
+
+
+ cyc_spline(tb->phi, tb->e, tablength, MY_2PI, tb->e2, comm->me == 0);
+ if (! tb->f_unspecified)
+ cyc_spline(tb->phi, tb->f, tablength, MY_2PI, tb->f2, comm->me == 0);
+}
+
+
+/* ----------------------------------------------------------------------
+ extract attributes from parameter line in table section
+ format of line: N value NOF DEGREES RADIANS
+ N is required, other params are optional
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::param_extract(Table *tb, char *line)
+{
+ //tb->theta0 = 180.0; <- equilibrium angles not supported
+ tb->ninput = 0;
+ tb->f_unspecified = false; //default
+ tb->use_degrees = true; //default
+
+ char *word = strtok(line," \t\n\r\f");
+ while (word) {
+ if (strcmp(word,"N") == 0) {
+ word = strtok(NULL," \t\n\r\f");
+ tb->ninput = atoi(word);
+ }
+ else if (strcmp(word,"NOF") == 0) {
+ tb->f_unspecified = true;
+ }
+ else if ((strcmp(word,"DEGREES") == 0) || (strcmp(word,"degrees") == 0)) {
+ tb->use_degrees = true;
+ }
+ else if ((strcmp(word,"RADIANS") == 0) || (strcmp(word,"radians") == 0)) {
+ tb->use_degrees = false;
+ }
+ else if (strcmp(word,"CHECKU") == 0) {
+ word = strtok(NULL," \t\n\r\f");
+ memory->sfree(checkU_fname);
+ memory->create(checkU_fname,strlen(word)+1,"dihedral_table:checkU");
+ strcpy(checkU_fname, word);
+ }
+ else if (strcmp(word,"CHECKF") == 0) {
+ word = strtok(NULL," \t\n\r\f");
+ memory->sfree(checkF_fname);
+ memory->create(checkF_fname,strlen(word)+1,"dihedral_table:checkF");
+ strcpy(checkF_fname, word);
+ }
+ // COMMENTING OUT: equilibrium angles are not supported
+ //else if (strcmp(word,"EQ") == 0) {
+ // word = strtok(NULL," \t\n\r\f");
+ // tb->theta0 = atof(word);
+ //}
+ else {
+ string err_msg("Invalid keyword in dihedral angle table parameters");
+ err_msg += string(" (") + string(word) + string(")");
+ error->one(FLERR,err_msg.c_str());
+ }
+ word = strtok(NULL," \t\n\r\f");
+ }
+
+ if (tb->ninput == 0)
+ error->one(FLERR,"Dihedral table parameters did not set N");
+
+} // DihedralTable::param_extract()
+
+
+/* ----------------------------------------------------------------------
+ broadcast read-in table info from proc 0 to other procs
+ this function communicates these values in Table:
+ ninput,phifile,efile,ffile,
+ f_unspecified,use_degrees
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::bcast_table(Table *tb)
+{
+ MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);
+
+ int me;
+ MPI_Comm_rank(world,&me);
+ if (me > 0) {
+ memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
+ memory->create(tb->efile,tb->ninput,"dihedral:efile");
+ memory->create(tb->ffile,tb->ninput,"dihedral:ffile");
+ }
+
+ MPI_Bcast(tb->phifile,tb->ninput,MPI_DOUBLE,0,world);
+ MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
+ MPI_Bcast(tb->ffile,tb->ninput,MPI_DOUBLE,0,world);
+
+ MPI_Bcast(&tb->f_unspecified,1,MPI_INT,0,world);
+ MPI_Bcast(&tb->use_degrees,1,MPI_INT,0,world);
+
+ // COMMENTING OUT: equilibrium angles are not supported
+ //MPI_Bcast(&tb->theta0,1,MPI_DOUBLE,0,world);
+}
+
+
+
+/* ----------------------------------------------------------------------
+ proc 0 writes to data file
+------------------------------------------------------------------------- */
+
+void DihedralTableCut::write_data(FILE *fp)
+{
+ for (int i = 1; i <= atom->ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g %g %g %g\n",i,
+ k1[i],phi1[i]*180.0/MY_PI,
+ k2[i],phi2[i]*180.0/MY_PI,
+ k3[i],phi3[i]*180.0/MY_PI);
+
+ fprintf(fp,"\nAngleAngleTorsion Coeffs\n\n");
+ for (int i = 1; i <= atom->ndihedraltypes; i++)
+ fprintf(fp,"%d %g %g %g\n",i,aat_k[i],
+ aat_theta0_1[i]*180.0/MY_PI,aat_theta0_2[i]*180.0/MY_PI);
+
+}
diff --git a/src/USER-MISC/dihedral_table_cut.h b/src/USER-MISC/dihedral_table_cut.h
new file mode 100644
index 0000000000000000000000000000000000000000..d01903c88b75d722aa2b5861a6cba4e00817ed87
--- /dev/null
+++ b/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.
+
+*/