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