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