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(comb3,PairComb3); 17 // clang-format on 18 #else 19 20 #ifndef LMP_PAIR_COMB3_H 21 #define LMP_PAIR_COMB3_H 22 23 #include "pair.h" 24 25 namespace LAMMPS_NS { 26 27 class PairComb3 : public Pair { 28 public: 29 PairComb3(class LAMMPS *); 30 virtual ~PairComb3(); 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 virtual double combqeq(double *, int &); 38 double enegtot; 39 40 static constexpr int NPARAMS_PER_LINE = 74; 41 42 protected: 43 // general potential parameters 44 struct Param { 45 int ielement, jelement, kelement, powermint; 46 int ielementgp, jelementgp, kelementgp; //element group 47 int ang_flag, pcn_flag, rad_flag, tor_flag; //angle, coordination,radical, torsion flag 48 double lami, lambda, alfi, alpha1, alpha2, alpha3, beta; 49 double pcos6, pcos5, pcos4, pcos3, pcos2, pcos1, pcos0; 50 double gamma, powerm, powern, bigA, bigB1, bigB2, bigB3; 51 double bigd, bigr, cut, cutsq, c1, c2, c3, c4; 52 double p6p0, p6p1, p6p2, p6p3, p6p4, p6p5, p6p6; 53 double ptork1, ptork2; 54 double addrepr, addrep, vdwflag; 55 double QU, QL, DU, DL, Qo, dQ, aB, bB, nD, bD, qmin, qmax; 56 double chi, dj, dk, dl, dm, esm, cmn1, cmn2, pcmn1, pcmn2; 57 double coulcut, lcut, lcutsq; 58 double veps, vsig, pcna, pcnb, pcnc, pcnd, polz, curl, pcross; 59 double paaa, pbbb; 60 double curlcut1, curlcut2, curl0; 61 }; 62 63 // general setups 64 double PI, PI2, PI4, PIsq; // PIs 65 double cutmin; // min cutoff for all elements 66 double cutmax; // max cutoff for all elements 67 double precision; // tolerance for QEq convergence 68 Param *params; // parameter set for an I-J-K interaction 69 int debug_eng1, debug_eng2, debug_fq; // logic controlling debugging outputs 70 int pack_flag; 71 72 // Short range neighbor list 73 void Short_neigh(); 74 int pgsize, oneatom; 75 int *sht_num, **sht_first; 76 MyPage<int> *ipage; 77 78 // loop up tables and flags 79 int nmax, **intype; 80 int pol_flag, polar; 81 double *qf, **bbij, *charge, *NCo; 82 double *esm, **fafb, **dfafb, **ddfafb, **phin, **dphin, **erpaw; 83 double **vvdw, **vdvdw; 84 double **afb, **dafb; 85 double **dpl, bytes; 86 double *xcctmp, *xchtmp, *xcotmp; 87 88 // additional carbon parameters 89 int cflag; 90 int nsplpcn, nsplrad, nspltor; 91 int maxx, maxy, maxz, maxxc, maxyc, maxconj; 92 int maxxcn[4]; 93 double vmaxxcn[4], dvmaxxcn[4]; 94 int ntab; 95 double iin2[16][2], iin3[64][3]; 96 double brad[4], btor[4], bbtor, ptorr; 97 double fi_tor[3], fj_tor[3], fk_tor[3], fl_tor[3]; 98 double radtmp, fi_rad[3], fj_rad[3], fk_rad[3]; 99 100 double ccutoff[6], ch_a[7]; 101 102 //COMB3-v18 arrays for CHO 103 // We wanna dynamic arrays 104 // C angle arrays, size = ntab+1 105 double pang[20001]; 106 double dpang[20001]; 107 double ddpang[20001]; 108 109 //coordination spline arrays 110 double pcn_grid[4][5][5][5]; 111 double pcn_gridx[4][5][5][5]; 112 double pcn_gridy[4][5][5][5]; 113 double pcn_gridz[4][5][5][5]; 114 double pcn_cubs[4][4][4][4][64]; 115 116 //coordination spline arrays 117 double rad_grid[3][5][5][11]; 118 double rad_gridx[3][5][5][11]; 119 double rad_gridy[3][5][5][11]; 120 double rad_gridz[3][5][5][11]; 121 double rad_spl[3][4][4][10][64]; 122 123 //torsion spline arrays 124 double tor_grid[1][5][5][11]; 125 double tor_gridx[1][5][5][11]; 126 double tor_gridy[1][5][5][11]; 127 double tor_gridz[1][5][5][11]; 128 double tor_spl[1][4][4][10][64]; 129 130 // initialization functions 131 void allocate(); 132 void read_lib(); 133 void setup_params(); 134 virtual void read_file(char *); 135 136 // cutoff functions 137 double comb_fc(double, Param *); 138 double comb_fc_d(double, Param *); 139 double comb_fc_curl(double, Param *); 140 double comb_fc_curl_d(double, Param *); 141 double comb_fccc(double); 142 double comb_fccc_d(double); 143 double comb_fcch(double); 144 double comb_fcch_d(double); 145 double comb_fccch(double); 146 double comb_fccch_d(double); 147 double comb_fcsw(double); 148 149 // short range terms 150 void attractive(Param *, Param *, Param *, double, double, double, double, double, double, double, 151 double *, double *, double *, double *, double *, int, double); 152 virtual void comb_fa(double, Param *, Param *, double, double, double &, double &); 153 virtual void repulsive(Param *, Param *, double, double &, int, double &, double, double); 154 155 // bond order terms 156 double comb_bij(double, Param *, double, int, double); 157 double comb_gijk(double, Param *, double); 158 void comb_gijk_d(double, Param *, double, double &, double &); 159 double zeta(Param *, Param *, double, double, double *, double *, int, double); 160 void comb_bij_d(double, Param *, double, int, double &, double &, double &, double &, double &, 161 double &, double); 162 void coord(Param *, double, int, double &, double &, double &, double &, double &, double); 163 void comb_zetaterm_d(double, double, double, double, double, double *, double, double *, double, 164 double *, double *, double *, Param *, Param *, Param *, double); 165 void costheta_d(double *, double, double *, double, double *, double *, double *); 166 void force_zeta(Param *, Param *, double, double, double, double &, double &, double &, double &, 167 double &, double &, double &, double &, double &, double &, double &, double &, 168 double &, int, double &, double, double, int, int, int, double, double, double); 169 void cntri_int(int, double, double, double, int, int, int, double &, double &, double &, double &, 170 Param *); 171 172 // Legendre polynomials 173 void selfp6p(Param *, Param *, double, double &, double &); 174 double ep6p(Param *, Param *, double, double, double *, double *, double &); 175 void fp6p(Param *, Param *, double, double, double *, double *, double *, double *, double *); 176 177 // long range q-dependent terms 178 double self(Param *, double); 179 void tables(); 180 void potal_calc(double &, double &, double &); 181 void tri_point(double, int &, int &, int &, double &, double &, double &); 182 void vdwaals(int, int, int, int, double, double, double, double, double &, double &); 183 void direct(Param *, Param *, int, int, int, double, double, double, double, double, double, 184 double, double, double &, double &, int, int); 185 void field(Param *, Param *, double, double, double, double &, double &); 186 int heaviside(double); 187 double switching(double); 188 double switching_d(double); 189 double chicut1, chicut2; 190 191 // radical terms 192 double rad_init(double, Param *, int, double &, double); 193 void rad_calc(double, Param *, Param *, double, double, int, int, double, double); 194 void rad_int(int, double, double, double, int, int, int, double &, double &, double &, double &); 195 void rad_forceik(Param *, double, double *, double, double); 196 void rad_force(Param *, double, double *, double); 197 198 // torsion terms 199 double bbtor1(int, Param *, Param *, double, double, double, double *, double *, double *, 200 double); //modified by TAO 201 void tor_calc(double, Param *, Param *, double, double, int, int, double, double); 202 void tor_int(int, double, double, double, int, int, int, double &, double &, double &, double &); 203 void tor_force(int, Param *, Param *, double, double, double, double, double *, double *, 204 double *); //modified by TAO 205 206 // charge force terms 207 double qfo_self(Param *, double); 208 void qfo_short(Param *, Param *, double, double, double, double &, double &, int, int, int); 209 void qfo_direct(Param *, Param *, int, int, int, double, double, double, double, double, double &, 210 double &, double, double, int, int); 211 void qfo_field(Param *, Param *, double, double, double, double &, double &); 212 void qfo_dipole(double, int, int, int, int, double, double *, double, double, double, double &, 213 double &, int, int); 214 void qsolve(double *); 215 216 // dipole - polarization terms 217 double dipole_self(Param *, int); 218 void dipole_init(Param *, Param *, double, double *, double, int, int, int, double, double, 219 double, double, double, int, int); 220 void dipole_calc(Param *, Param *, double, double, double, double, double, int, int, int, double, 221 double, double, double, double, int, int, double &, double &, double *); 222 223 // communication functions 224 int pack_reverse_comm(int, int, double *); 225 void unpack_reverse_comm(int, int *, double *); 226 int pack_forward_comm(int, int *, double *, int, int *); 227 void unpack_forward_comm(int, int, double *); 228 }; 229 } // namespace LAMMPS_NS 230 231 #endif 232 #endif 233 234 /* ERROR/WARNING messages: 235 236 E: Illegal ... command 237 238 Self-explanatory. Check the input script syntax and compare to the 239 documentation for the command. You can use -echo screen as a 240 command-line option when running LAMMPS to see the offending line. 241 242 E: Incorrect args for pair coefficients 243 244 Self-explanatory. Check the input script or data file. 245 246 E: Pair style COMB3 requires atom IDs 247 248 This is a requirement to use the COMB3 potential. 249 250 E: Pair style COMB3 requires newton pair on 251 252 See the newton command. This is a restriction to use the COMB3 253 potential. 254 255 E: Pair style COMB3 requires atom attribute q 256 257 Self-explanatory. 258 259 E: All pair coeffs are not set 260 261 All pair coefficients must be set in the data file or by the 262 pair_coeff command before running a simulation. 263 264 E: Cannot open COMB3 lib.comb3 file 265 266 The COMB3 library file cannot be opened. Check that the path and name 267 are correct. 268 269 E: Cannot open COMB3 potential file %s 270 271 The specified COMB3 potential file cannot be opened. Check that the 272 path and name are correct. 273 274 E: Incorrect format in COMB3 potential file 275 276 Incorrect number of words per line in the potential file. 277 278 E: Illegal COMB3 parameter 279 280 One or more of the coefficients defined in the potential file is 281 invalid. 282 283 E: Potential file has duplicate entry 284 285 The potential file has more than one entry for the same element. 286 287 E: Potential file is missing an entry 288 289 The potential file does not have a needed entry. 290 291 E: Neighbor list overflow, boost neigh_modify one 292 293 There are too many neighbors of a single atom. Use the neigh_modify 294 command to increase the max number of neighbors allowed for one atom. 295 You may also want to boost the page size. 296 297 E: Error in vdw spline: inner radius > outer radius 298 299 A pre-tabulated spline is invalid. Likely a problem with the 300 potential parameters. 301 302 */ 303