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