1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "pair_buck_coul_msm.h"
16 #include <cmath>
17 #include <cstring>
18 #include "atom.h"
19 #include "force.h"
20 #include "kspace.h"
21 #include "neigh_list.h"
22 #include "memory.h"
23 #include "error.h"
24 
25 using namespace LAMMPS_NS;
26 
27 
28 /* ---------------------------------------------------------------------- */
29 
PairBuckCoulMSM(LAMMPS * lmp)30 PairBuckCoulMSM::PairBuckCoulMSM(LAMMPS *lmp) : PairBuckCoulLong(lmp)
31 {
32   ewaldflag = pppmflag = 0;
33   msmflag = 1;
34   nmax = 0;
35   ftmp = nullptr;
36 }
37 
38 /* ---------------------------------------------------------------------- */
39 
~PairBuckCoulMSM()40 PairBuckCoulMSM::~PairBuckCoulMSM()
41 {
42   if (ftmp) memory->destroy(ftmp);
43 }
44 
45 /* ---------------------------------------------------------------------- */
46 
compute(int eflag,int vflag)47 void PairBuckCoulMSM::compute(int eflag, int vflag)
48 {
49   int i,j,ii,jj,inum,jnum,itype,jtype;
50   double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair,fcoul;
51   double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj;
52   double egamma,fgamma,prefactor;
53   double r,rexp;
54   int *ilist,*jlist,*numneigh,**firstneigh;
55   int eflag_old = eflag;
56 
57   if (force->kspace->scalar_pressure_flag && vflag) {
58     if (vflag > 2)
59       error->all(FLERR,"Must use 'kspace_modify pressure/scalar no' "
60         "to obtain per-atom virial with kspace_style MSM");
61 
62     if (atom->nmax > nmax) {
63       if (ftmp) memory->destroy(ftmp);
64       nmax = atom->nmax;
65       memory->create(ftmp,nmax,3,"pair:ftmp");
66     }
67     memset(&ftmp[0][0],0,nmax*3*sizeof(double));
68 
69     // must switch on global energy computation if not already on
70 
71     if (eflag == 0 || eflag == 2) {
72       eflag++;
73     }
74   }
75 
76   evdwl = ecoul = 0.0;
77   ev_init(eflag,vflag);
78 
79   double **x = atom->x;
80   double **f = atom->f;
81   double *q = atom->q;
82   int *type = atom->type;
83   int nlocal = atom->nlocal;
84   double *special_coul = force->special_coul;
85   double *special_lj = force->special_lj;
86   int newton_pair = force->newton_pair;
87   double qqrd2e = force->qqrd2e;
88 
89   inum = list->inum;
90   ilist = list->ilist;
91   numneigh = list->numneigh;
92   firstneigh = list->firstneigh;
93 
94   // loop over neighbors of my atoms
95 
96   for (ii = 0; ii < inum; ii++) {
97     i = ilist[ii];
98     qtmp = q[i];
99     xtmp = x[i][0];
100     ytmp = x[i][1];
101     ztmp = x[i][2];
102     itype = type[i];
103     jlist = firstneigh[i];
104     jnum = numneigh[i];
105 
106     for (jj = 0; jj < jnum; jj++) {
107       j = jlist[jj];
108       factor_lj = special_lj[sbmask(j)];
109       factor_coul = special_coul[sbmask(j)];
110       j &= NEIGHMASK;
111 
112       delx = xtmp - x[j][0];
113       dely = ytmp - x[j][1];
114       delz = ztmp - x[j][2];
115       rsq = delx*delx + dely*dely + delz*delz;
116       jtype = type[j];
117 
118       if (rsq < cutsq[itype][jtype]) {
119         r2inv = 1.0/rsq;
120         r = sqrt(rsq);
121 
122         if (rsq < cut_coulsq) {
123           prefactor = qqrd2e * qtmp*q[j]/r;
124           egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
125           fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
126           forcecoul = prefactor * fgamma;
127           if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
128         } else forcecoul = 0.0;
129 
130         if (rsq < cut_ljsq[itype][jtype]) {
131           r6inv = r2inv*r2inv*r2inv;
132           rexp = exp(-r*rhoinv[itype][jtype]);
133           forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
134         } else forcebuck = 0.0;
135 
136         if (!(force->kspace->scalar_pressure_flag && vflag)) {
137           fpair = (forcecoul + factor_lj*forcebuck) * r2inv;
138 
139           f[i][0] += delx*fpair;
140           f[i][1] += dely*fpair;
141           f[i][2] += delz*fpair;
142           if (newton_pair || j < nlocal) {
143             f[j][0] -= delx*fpair;
144             f[j][1] -= dely*fpair;
145             f[j][2] -= delz*fpair;
146           }
147         } else {
148           // separate Buck and Coulombic forces
149 
150           fpair = (factor_lj*forcebuck) * r2inv;
151 
152           f[i][0] += delx*fpair;
153           f[i][1] += dely*fpair;
154           f[i][2] += delz*fpair;
155           if (newton_pair || j < nlocal) {
156             f[j][0] -= delx*fpair;
157             f[j][1] -= dely*fpair;
158             f[j][2] -= delz*fpair;
159           }
160 
161           fcoul = (forcecoul) * r2inv;
162 
163           ftmp[i][0] += delx*fcoul;
164           ftmp[i][1] += dely*fcoul;
165           ftmp[i][2] += delz*fcoul;
166           if (newton_pair || j < nlocal) {
167             ftmp[j][0] -= delx*fcoul;
168             ftmp[j][1] -= dely*fcoul;
169             ftmp[j][2] -= delz*fcoul;
170           }
171         }
172 
173         if (eflag) {
174           if (rsq < cut_coulsq) {
175             ecoul = prefactor*egamma;
176             if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
177           } else ecoul = 0.0;
178           if (eflag_old && rsq < cut_ljsq[itype][jtype]) {
179             evdwl = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
180               offset[itype][jtype];
181             evdwl *= factor_lj;
182           } else evdwl = 0.0;
183         }
184 
185         if (evflag) ev_tally(i,j,nlocal,newton_pair,
186                              evdwl,ecoul,fpair,delx,dely,delz);
187       }
188     }
189   }
190 
191   if (vflag_fdotr) virial_fdotr_compute();
192 
193   if (force->kspace->scalar_pressure_flag && vflag) {
194     for (i = 0; i < 3; i++) virial[i] += force->pair->eng_coul/3.0;
195     for (i = 0; i < nmax; i++) {
196       f[i][0] += ftmp[i][0];
197       f[i][1] += ftmp[i][1];
198       f[i][2] += ftmp[i][2];
199     }
200   }
201 }
202 
203 /* ---------------------------------------------------------------------- */
204 
single(int i,int j,int itype,int jtype,double rsq,double factor_coul,double factor_lj,double & fforce)205 double PairBuckCoulMSM::single(int i, int j, int itype, int jtype,
206                                 double rsq,
207                                 double factor_coul, double factor_lj,
208                                 double &fforce)
209 {
210   double r2inv,r6inv,r,rexp,egamma,fgamma,prefactor;
211   double forcecoul,forcebuck,phicoul,phibuck;
212 
213   r2inv = 1.0/rsq;
214   if (rsq < cut_coulsq) {
215     r = sqrt(rsq);
216     prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
217     egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul);
218     fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul);
219     forcecoul = prefactor * fgamma;
220     if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
221   } else forcecoul = 0.0;
222   if (rsq < cut_ljsq[itype][jtype]) {
223     r6inv = r2inv*r2inv*r2inv;
224     r = sqrt(rsq);
225     rexp = exp(-r*rhoinv[itype][jtype]);
226     forcebuck = buck1[itype][jtype]*r*rexp - buck2[itype][jtype]*r6inv;
227   } else forcebuck = 0.0;
228   fforce = (forcecoul + factor_lj*forcebuck) * r2inv;
229 
230   double eng = 0.0;
231   if (rsq < cut_coulsq) {
232     phicoul = prefactor*egamma;
233     if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
234     eng += phicoul;
235   }
236   if (rsq < cut_ljsq[itype][jtype]) {
237     phibuck = a[itype][jtype]*rexp - c[itype][jtype]*r6inv -
238       offset[itype][jtype];
239     eng += factor_lj*phibuck;
240   }
241   return eng;
242 }
243 
244 /* ---------------------------------------------------------------------- */
245 
extract(const char * str,int & dim)246 void *PairBuckCoulMSM::extract(const char *str, int &dim)
247 {
248   dim = 0;
249   if (strcmp(str,"cut_coul") == 0) return (void *) &cut_coul;
250   dim = 2;
251   if (strcmp(str,"a") == 0) return (void *) a;
252   if (strcmp(str,"c") == 0) return (void *) c;
253 
254   return nullptr;
255 }
256