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(smtbq,PairSMTBQ);
17 // clang-format on
18 #else
19 
20 #ifndef LMP_PAIR_SMTBQ_H
21 #define LMP_PAIR_SMTBQ_H
22 
23 #include "pair.h"
24 
25 namespace LAMMPS_NS {
26 
27 class PairSMTBQ : public Pair {
28  public:
29   PairSMTBQ(class LAMMPS *);
30   virtual ~PairSMTBQ();
31   virtual void compute(int, int);
32 
33   void settings(int, char **);
34   void coeff(int, char **);
35   void init_style();
36   double init_one(int, int);
37   double memory_usage();
38 
39  protected:
40   struct Param {
41     double sto, n0, ne, chi, dj;
42     double R, dzeta;    //Rayon slater
43     double cutsq;
44     double qform, masse;    // Charge formelle
45     char *nom;
46   };
47 
48   struct Intparam {
49     double a, p, ksi, q, r0, dc1, dc2;
50     double abuck, rhobuck, neig_cut, aOO, bOO, r1OO, r2OO;
51     char *typepot, *mode;
52     int intsm;
53   };
54 
55   double rmin, dr, ds;    // table parameter
56   int kmax;
57   bigint Qstep;        // Frequency of charge resolution
58   double precision;    // accuracy of convergence
59   int loopmax;         // max of iteration
60 
61   double cutmax;       // max cutoff for all elements
62   char *QEqMode;       // name of QEqMode
63   char *Bavard;        // Verbose parameter
64   char *writepot;      // write or not the electronegativity component
65   char *writeenerg;    // write or not the energy component
66   char *QInitMode;     // mode of initialization of charges
67   double zlim1QEq;     // z limit for QEq equilibration
68   double zlim2QEq;     // z limit for QEq equilibration
69   double QOxInit;      // Initial charge for oxygen atoms (if the charge is not specified)
70   int maxintparam;     // max # of interaction sets
71   int maxintsm;        // max # of interaction SM
72   double r1Coord, r2Coord;
73   Param *params;          // parameter set for an I atom
74   Intparam *intparams;    // parameter set for an I interaction
75 
76   int nmax, *nQEqall, *nQEqaall, *nQEqcall;
77   double *qf, *q1, *q2, Nevery, Neverypot;
78 
79   // Coulombian interaction
80   double *esm, **fafb, **dfafb, *fafbOxOxSurf, *dfafbOxOxSurf, *fafbTiOxSurf, *dfafbTiOxSurf;
81   double *potqn, *dpotqn, Vself, *Zsm, *coord, *fafbOxOxBB, *dfafbOxOxBB, *fafbTiOxBB, *dfafbTiOxBB;
82   int **intype, **coultype;
83   int *NCo;
84   double coordOxBulk, coordOxSurf, ROxSurf, coordOxBB, ROxBB;
85 
86   // Covalent interaction
87   double *ecov, *potmad, *potself, *potcov;    //, *chimet;
88   double **tabsmb, **dtabsmb, **tabsmr, **dtabsmr, *sbcov, *sbmet;
89   double ncov;
90 
91   // Neighbor Table
92   int nteam, cluster, *hybrid;
93   int *nvsm, **vsm, *flag_QEq;
94 
95   // Parallelisation
96   int me, nproc;
97   double *tab_comm;
98 
99   // GAMMAS function
100   double *fct;
101 
102   // HERE its routines
103   // =====================================
104   void allocate();
105   virtual void read_file(char *);
106 
107   void tabsm();
108   void tabqeq();
109 
110   void potqeq(int, int, double, double, double, double &, int, double &);
111   void pot_ES(int, int, double, double &);
112   void pot_ES2(int, int, double, double &);
113 
114   double self(Param *, double);
115   double qfo_self(Param *, double);
116 
117   virtual void repulsive(Intparam *, double, int, int, double &, int, double &);
118   virtual void rep_OO(Intparam *, double, double &, int, double &);
119   virtual void Attr_OO(Intparam *, double, double &, int, double &);
120   virtual void attractive(Intparam *, double, int, int, double, int, double);
121   void f_att(Intparam *, int, int, double, double &);
122   void pot_COV(Param *, int, double &);
123   double potmet(Intparam *, double, int, double, int, double);
124 
125   double fcoupure(double, double, double);
126   double fcoupured(double, double, double);
127 
128   double fcoup2(double, double, double);
129   double Intfcoup2(double, double, double);
130   double Primfcoup2(double, double, double);
131 
132   void groupBulkFromSlab_QEq();
133   void groupQEqAllParallel_QEq();
134   void groupQEqAll_QEq();
135   void groupSurface_QEq();
136   void QForce_charge(int);
137   void Charge();
138   void Init_charge(int *, int *, int *);
139   void CheckEnergyVSForce();
140 
141   // ===========================================
142   // Communication pack
143   int pack_forward_comm(int, int *, double *, int, int *);
144   void unpack_forward_comm(int, int, double *);
145   int pack_reverse_comm(int, int, double *);
146   void unpack_reverse_comm(int, int *, double *);
147   void forward(double *);
148   void reverse(double *);
149   void forward_int(int *);
150   void reverse_int(int *);
151 
152   int Tokenize(char *, char ***);
153 
min(const T & a,const T & b)154   template <class T> const T &min(const T &a, const T &b)
155   {
156     return !(b < a) ? a : b;    // or: return !comp(b,a)?a:b; for the comp version
157   }
158 
159   // =============================================
160   // Gammas function
161   void gammas(double &, double &, double &, double &, double &, double &, double &, double &,
162               double &, double &, double &, double &, double &, double &, double &);
163 
164   void css(double &, double, double, double, double, double, double &, double &, double &, double &,
165            double &, double &, double &, double &, double &);
166 
167   double coeffs(int, int, int);
168 
169   double binm(int, int);
170   void caintgs(double, int, double *);
171   void cbintgs(double, int, double *);
172 
173   // =====================================
174   // short range neighbor list
175 
176   void add_pages(int howmany = 1);
177   //  void Short_neigh();
178   int maxpage, pgsize, oneatom, **pages;
179   double cutmin;
180 };
181 
182 }    // namespace LAMMPS_NS
183 
184 #endif
185 #endif
186