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 #ifdef PAIR_CLASS
15 // clang-format off
16 PairStyle(comb,PairComb);
17 // clang-format on
18 #else
19 
20 #ifndef LMP_PAIR_COMB_H
21 #define LMP_PAIR_COMB_H
22 
23 #include "pair.h"
24 
25 namespace LAMMPS_NS {
26 
27 class PairComb : public Pair {
28  public:
29   PairComb(class LAMMPS *);
30   virtual ~PairComb();
31   virtual void compute(int, int);
32   void settings(int, char **);
33   void coeff(int, char **);
34   void init_style();
35   double init_one(int, int);
36   double memory_usage();
37 
38   virtual double yasu_char(double *, int &);
39   double enegtot;
40 
41   static constexpr int NPARAMS_PER_LINE = 49;
42 
43  protected:
44   struct Param {
45     double lam11, lam12, lam21, lam22;
46     double c, d, h;
47     double gamma, powerm;
48     double powern, beta;
49     double biga1, biga2, bigb1, bigb2;
50     double bigd, bigr;
51     double cut, cutsq;
52     double c1, c2, c3, c4;
53     double plp1, plp3, plp6, a123, aconf;
54     double rlm1, rlm2;
55     double romiga, romigb, romigc, romigd, addrep;
56     double QU1, QL1, DU1, DL1, Qo1, dQ1, aB1, bB1, nD1, bD1;
57     double QU2, QL2, DU2, DL2, Qo2, dQ2, aB2, bB2, nD2, bD2;
58     double chi, dj, dk, dl, dm, esm1, esm2, cmn1, cmn2, cml1, cml2;
59     double coulcut, lcut, lcutsq, hfocor;
60     int ielement, jelement, kelement;
61     int powermint;
62   };
63 
64   double cutmax;    // max cutoff for all elements
65   double precision;
66   Param *params;    // parameter set for an I-J-K interaction
67 
68   int nmax;
69   double *qf;
70 
71   double *esm, **fafb, **dfafb, **ddfafb, **phin, **dphin, **erpaw;
72   double *charge;
73   int **intype, *typeno;
74   int *NCo, cor_flag, cuo_flag, cuo_flag1, cuo_flag2;
75   double **bbij;
76 
77   int pgsize;                   // size of neighbor page
78   int oneatom;                  // max # of neighbors for one atom
79   int *sht_num, **sht_first;    // short-range neighbor list
80   MyPage<int> *ipage;           // neighbor list pages
81   double cutmin;
82 
83   void allocate();
84   virtual void read_file(char *);
85   void setup_params();
86   virtual void repulsive(Param *, double, double &, int, double &, double, double);
87   double zeta(Param *, double, double, double *, double *);
88   void force_zeta(Param *, int, int, int, double, double, double, double, double &, double &,
89                   double &);
90   void attractive(Param *, double, double, double, double *, double *, double *, double *,
91                   double *);
92   double elp(Param *, double, double, double *, double *);
93   void flp(Param *, double, double, double *, double *, double *, double *, double *);
94   double comb_fc(double, Param *);
95   double comb_fc_d(double, Param *);
96   double comb_fc2(double);
97   double comb_fc2_d(double);
98   double comb_fc3(double);
99   double comb_fc3_d(double);
100   virtual double comb_fa(double, Param *, double, double);
101   virtual double comb_fa_d(double, Param *, double, double);
102   double comb_bij(double, Param *);
103   double comb_bij_d(double, Param *);
104 
comb_gijk(const double costheta,const Param * const param)105   inline double comb_gijk(const double costheta, const Param *const param) const
106   {
107     const double comb_c = param->c * param->c;
108     const double comb_d = param->d * param->d;
109     const double hcth = param->h - costheta;
110 
111     return param->gamma * (1.0 + comb_c / comb_d - comb_c / (comb_d + hcth * hcth));
112   }
113 
comb_gijk_d(const double costheta,const Param * const param)114   inline double comb_gijk_d(const double costheta, const Param *const param) const
115   {
116     const double comb_c = param->c * param->c;
117     const double comb_d = param->d * param->d;
118     const double hcth = param->h - costheta;
119     const double numerator = -2.0 * comb_c * hcth;
120     const double denominator = 1.0 / (comb_d + hcth * hcth);
121     return param->gamma * numerator * denominator * denominator;
122   }
123 
124   void comb_zetaterm_d(double, double *, double, double *, double, double *, double *, double *,
125                        Param *);
126   void costheta_d(double *, double, double *, double, double *, double *, double *);
127   double self(Param *, double, double);
128   void sm_table();
129   void potal_calc(double &, double &, double &);
130   void tri_point(double, int &, int &, int &, double &, double &, double &, int &);
131   void direct(int, int, int, int, double, double, double, double, double, double, double, double,
132               double, double &, double &);
133   void field(Param *, double, double, double, double &, double &);
134   double qfo_self(Param *, double, double);
135   void qfo_short(Param *, int, int, double, double, double, double &, double &);
136   void qfo_direct(int, int, int, int, double, double, double, double, double, double &);
137   void qfo_field(Param *, double, double, double, double &, double &);
138   void qsolve(double *);
139   void Over_cor(Param *, double, int, double &, double &);
140   int pack_reverse_comm(int, int, double *);
141   void unpack_reverse_comm(int, int *, double *);
142   int pack_forward_comm(int, int *, double *, int, int *);
143   void unpack_forward_comm(int, int, double *);
144 
145   void Short_neigh();
146 };
147 
148 }    // namespace LAMMPS_NS
149 
150 #endif
151 #endif
152 
153 /* ERROR/WARNING messages:
154 
155 E: Illegal ... command
156 
157 Self-explanatory.  Check the input script syntax and compare to the
158 documentation for the command.  You can use -echo screen as a
159 command-line option when running LAMMPS to see the offending line.
160 
161 E: Incorrect args for pair coefficients
162 
163 Self-explanatory.  Check the input script or data file.
164 
165 E: Pair style COMB requires atom IDs
166 
167 This is a requirement to use the AIREBO potential.
168 
169 E: Pair style COMB requires newton pair on
170 
171 See the newton command.  This is a restriction to use the COMB
172 potential.
173 
174 E: Pair style COMB requires atom attribute q
175 
176 Self-explanatory.
177 
178 E: All pair coeffs are not set
179 
180 All pair coefficients must be set in the data file or by the
181 pair_coeff command before running a simulation.
182 
183 E: Cannot open COMB potential file %s
184 
185 The specified COMB potential file cannot be opened.  Check that the
186 path and name are correct.
187 
188 E: Incorrect format in COMB potential file
189 
190 Incorrect number of words per line in the potential file.
191 
192 E: Illegal COMB parameter
193 
194 One or more of the coefficients defined in the potential file is
195 invalid.
196 
197 E: Potential file has duplicate entry
198 
199 The potential file has more than one entry for the same element.
200 
201 E: Potential file is missing an entry
202 
203 The potential file does not have a needed entry.
204 
205 W: Pair COMB charge %.10f with force %.10f hit min barrier
206 
207 Something is possibly wrong with your model.
208 
209 W: Pair COMB charge %.10f with force %.10f hit max barrier
210 
211 Something is possibly wrong with your model.
212 
213 E: Neighbor list overflow, boost neigh_modify one
214 
215 There are too many neighbors of a single atom.  Use the neigh_modify
216 command to increase the max number of neighbors allowed for one atom.
217 You may also want to boost the page size.
218 
219 */
220