1 /* -*- c++ -*- ----------------------------------------------------------
2    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3    https://www.lammps.org/, Sandia National Laboratories
4    Steve Plimpton, sjplimp@sandia.gov
5 
6    Copyright (2003) Sandia Corporation.  Under the terms of Contract
7    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8    certain rights in this software.  This software is distributed under
9    the GNU General Public License.
10 
11    See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13 
14 /* ----------------------------------------------------------------------
15    Contributing author: Andrew Jewett (jewett.aij at g mail)
16 ------------------------------------------------------------------------- */
17 
18 #ifdef DIHEDRAL_CLASS
19 // clang-format off
20 DihedralStyle(table,DihedralTable);
21 // clang-format on
22 #else
23 
24 #ifndef LMP_DIHEDRAL_TABLE_H
25 #define LMP_DIHEDRAL_TABLE_H
26 #include "dihedral.h"
27 
28 namespace LAMMPS_NS {
29 
30 class DihedralTable : public Dihedral {
31  public:
32   DihedralTable(class LAMMPS *);
33   virtual ~DihedralTable();
34   virtual void compute(int, int);
35   void settings(int, char **);
36   virtual void coeff(int, char **);
37   void write_restart(FILE *);
38   void read_restart(FILE *);
39   void write_restart_settings(FILE *);
40   void read_restart_settings(FILE *);
41   double single(int type, int i1, int i2, int i3, int i4);
42 
43  protected:
44   int tabstyle, tablength;
45   std::string checkU_fname;
46   std::string checkF_fname;
47 
48   struct Table {
49     int ninput;
50     int f_unspecified;    // boolean (but MPI does not like type "bool")
51     int use_degrees;      // boolean (but MPI does not like type "bool")
52     double *phifile, *efile, *ffile;
53     double *e2file, *f2file;
54     double delta, invdelta, deltasq6;
55     double *phi, *e, *de, *f, *df, *e2, *f2;
56   };
57 
58   int ntables;
59   Table *tables;
60   int *tabindex;
61 
62   virtual void allocate();
63   void null_table(Table *);
64   void free_table(Table *);
65   void read_table(Table *, char *, char *);
66   void bcast_table(Table *);
67   void spline_table(Table *);
68   void compute_table(Table *);
69 
70   void param_extract(Table *, char *);
71 
72   // --------------------------------------------
73   // ------------ inline functions --------------
74   // --------------------------------------------
75 
76   // -----------------------------------------------------------
77   //   uf_lookup()
78   //   quickly calculate the potential u and force f at angle x,
79   //   using the internal tables tb->e and tb->f (evenly spaced)
80   // -----------------------------------------------------------
81   enum { LINEAR, SPLINE };
82 
uf_lookup(int type,double x,double & u,double & f)83   inline void uf_lookup(int type, double x, double &u, double &f)
84   {
85     Table *tb = &tables[tabindex[type]];
86     double x_over_delta = x * tb->invdelta;
87     int i = static_cast<int>(x_over_delta);
88     double a;
89     double b = x_over_delta - i;
90     // Apply periodic boundary conditions to indices i and i+1
91     if (i >= tablength) i -= tablength;
92     int ip1 = i + 1;
93     if (ip1 >= tablength) ip1 -= tablength;
94 
95     switch (tabstyle) {
96       case LINEAR:
97         u = tb->e[i] + b * tb->de[i];
98         f = tb->f[i] + b * tb->df[i];    //<--works even if tb->f_unspecified==true
99         break;
100       case SPLINE:
101         a = 1.0 - b;
102         u = a * tb->e[i] + b * tb->e[ip1] +
103             ((a * a * a - a) * tb->e2[i] + (b * b * b - b) * tb->e2[ip1]) * tb->deltasq6;
104         if (tb->f_unspecified)
105           //Formula below taken from equation3.3.5 of "numerical recipes in c"
106           //"f"=-derivative of e with respect to x (or "phi" in this case)
107           f = (tb->e[i] - tb->e[ip1]) * tb->invdelta +
108               ((3.0 * a * a - 1.0) * tb->e2[i] + (1.0 - 3.0 * b * b) * tb->e2[ip1]) * tb->delta /
109                   6.0;
110         else
111           f = a * tb->f[i] + b * tb->f[ip1] +
112               ((a * a * a - a) * tb->f2[i] + (b * b * b - b) * tb->f2[ip1]) * tb->deltasq6;
113         break;
114     }    // switch(tabstyle)
115   }      // uf_lookup()
116 
117   // ----------------------------------------------------------
118   //    u_lookup()
119   //  quickly calculate the potential u at angle x using tb->e
120   //-----------------------------------------------------------
121 
u_lookup(int type,double x,double & u)122   inline void u_lookup(int type, double x, double &u)
123   {
124     Table *tb = &tables[tabindex[type]];
125     int N = tablength;
126 
127     //  i = static_cast<int> ((x - tb->lo) * tb->invdelta); <-general version
128     double x_over_delta = x * tb->invdelta;
129     int i = static_cast<int>(x_over_delta);
130     double b = x_over_delta - i;
131 
132     // Apply periodic boundary conditions to indices i and i+1
133     if (i >= N) i -= N;
134     int ip1 = i + 1;
135     if (ip1 >= N) ip1 -= N;
136 
137     if (tabstyle == LINEAR) {
138       u = tb->e[i] + b * tb->de[i];
139     } else if (tabstyle == SPLINE) {
140       double a = 1.0 - b;
141       u = a * tb->e[i] + b * tb->e[ip1] +
142           ((a * a * a - a) * tb->e2[i] + (b * b * b - b) * tb->e2[ip1]) * tb->deltasq6;
143     }
144   }    // u_lookup()
145 
146 };    //class DihedralTable
147 
148 }    // namespace LAMMPS_NS
149 
150 #endif    //#ifndef LMP_DIHEDRAL_TABLE_H
151 #endif    //#ifdef DIHEDRAL_CLASS ... #else
152