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: Trung Dac Nguyen (ORNL)
17 References: Fennell and Gezelter, JCP 124, 234104 (2006)
18 ------------------------------------------------------------------------- */
19
20 #include "pair_coul_dsf.h"
21
22 #include <cmath>
23 #include <cstring>
24 #include "atom.h"
25 #include "comm.h"
26 #include "force.h"
27 #include "neighbor.h"
28 #include "neigh_list.h"
29 #include "memory.h"
30 #include "math_const.h"
31 #include "error.h"
32
33 using namespace LAMMPS_NS;
34 using namespace MathConst;
35
36 #define EWALD_F 1.12837917
37 #define EWALD_P 0.3275911
38 #define A1 0.254829592
39 #define A2 -0.284496736
40 #define A3 1.421413741
41 #define A4 -1.453152027
42 #define A5 1.061405429
43
44 /* ---------------------------------------------------------------------- */
45
PairCoulDSF(LAMMPS * lmp)46 PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) {}
47
48 /* ---------------------------------------------------------------------- */
49
~PairCoulDSF()50 PairCoulDSF::~PairCoulDSF()
51 {
52 if (copymode) return;
53
54 if (allocated) {
55 memory->destroy(setflag);
56 memory->destroy(cutsq);
57 }
58 }
59
60 /* ---------------------------------------------------------------------- */
61
compute(int eflag,int vflag)62 void PairCoulDSF::compute(int eflag, int vflag)
63 {
64 int i,j,ii,jj,inum,jnum;
65 double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,ecoul,fpair;
66 double r,rsq,forcecoul,factor_coul;
67 double prefactor,erfcc,erfcd,t;
68 int *ilist,*jlist,*numneigh,**firstneigh;
69
70 ecoul = 0.0;
71 ev_init(eflag,vflag);
72
73 double **x = atom->x;
74 double **f = atom->f;
75 double *q = atom->q;
76 int nlocal = atom->nlocal;
77 double *special_coul = force->special_coul;
78 int newton_pair = force->newton_pair;
79 double qqrd2e = force->qqrd2e;
80
81 inum = list->inum;
82 ilist = list->ilist;
83 numneigh = list->numneigh;
84 firstneigh = list->firstneigh;
85
86 // loop over neighbors of my atoms
87
88 for (ii = 0; ii < inum; ii++) {
89 i = ilist[ii];
90 qtmp = q[i];
91 xtmp = x[i][0];
92 ytmp = x[i][1];
93 ztmp = x[i][2];
94 jlist = firstneigh[i];
95 jnum = numneigh[i];
96
97 if (eflag) {
98 double e_self = -(e_shift/2.0 + alpha/MY_PIS) * qtmp*qtmp*qqrd2e;
99 ev_tally(i,i,nlocal,0,0.0,e_self,0.0,0.0,0.0,0.0);
100 }
101
102 for (jj = 0; jj < jnum; jj++) {
103 j = jlist[jj];
104 factor_coul = special_coul[sbmask(j)];
105 j &= NEIGHMASK;
106
107 delx = xtmp - x[j][0];
108 dely = ytmp - x[j][1];
109 delz = ztmp - x[j][2];
110 rsq = delx*delx + dely*dely + delz*delz;
111
112 if (rsq < cut_coulsq) {
113 r = sqrt(rsq);
114 prefactor = qqrd2e*qtmp*q[j]/r;
115 erfcd = exp(-alpha*alpha*rsq);
116 t = 1.0 / (1.0 + EWALD_P*alpha*r);
117 erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
118
119 forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +
120 r*f_shift) * r;
121 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
122 fpair = forcecoul / rsq;
123
124 f[i][0] += delx*fpair;
125 f[i][1] += dely*fpair;
126 f[i][2] += delz*fpair;
127 if (newton_pair || j < nlocal) {
128 f[j][0] -= delx*fpair;
129 f[j][1] -= dely*fpair;
130 f[j][2] -= delz*fpair;
131 }
132
133 if (eflag) {
134 ecoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
135 if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor;
136 } else ecoul = 0.0;
137
138 if (evflag) ev_tally(i,j,nlocal,newton_pair,
139 0.0,ecoul,fpair,delx,dely,delz);
140 }
141 }
142 }
143
144 if (vflag_fdotr) virial_fdotr_compute();
145 }
146
147 /* ----------------------------------------------------------------------
148 allocate all arrays
149 ------------------------------------------------------------------------- */
150
allocate()151 void PairCoulDSF::allocate()
152 {
153 allocated = 1;
154 int n = atom->ntypes;
155
156 memory->create(setflag,n+1,n+1,"pair:setflag");
157 for (int i = 1; i <= n; i++)
158 for (int j = i; j <= n; j++)
159 setflag[i][j] = 0;
160
161 memory->create(cutsq,n+1,n+1,"pair:cutsq");
162 }
163
164 /* ----------------------------------------------------------------------
165 global settings
166 ------------------------------------------------------------------------- */
167
settings(int narg,char ** arg)168 void PairCoulDSF::settings(int narg, char **arg)
169 {
170 if (narg != 2) error->all(FLERR,"Illegal pair_style command");
171
172 alpha = utils::numeric(FLERR,arg[0],false,lmp);
173 cut_coul = utils::numeric(FLERR,arg[1],false,lmp);
174 }
175
176 /* ----------------------------------------------------------------------
177 set coeffs for one or more type pairs
178 ------------------------------------------------------------------------- */
179
coeff(int narg,char ** arg)180 void PairCoulDSF::coeff(int narg, char **arg)
181 {
182 if (narg != 2) error->all(FLERR,"Incorrect args for pair coefficients");
183 if (!allocated) allocate();
184
185 int ilo,ihi,jlo,jhi;
186 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
187 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
188
189 int count = 0;
190 for (int i = ilo; i <= ihi; i++) {
191 for (int j = MAX(jlo,i); j <= jhi; j++) {
192 setflag[i][j] = 1;
193 count++;
194 }
195 }
196
197 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
198 }
199
200 /* ----------------------------------------------------------------------
201 init specific to this pair style
202 ------------------------------------------------------------------------- */
203
init_style()204 void PairCoulDSF::init_style()
205 {
206 if (!atom->q_flag)
207 error->all(FLERR,"Pair style coul/dsf requires atom attribute q");
208
209 neighbor->request(this,instance_me);
210
211 cut_coulsq = cut_coul * cut_coul;
212 double erfcc = erfc(alpha*cut_coul);
213 double erfcd = exp(-alpha*alpha*cut_coul*cut_coul);
214 f_shift = -(erfcc/cut_coulsq + 2.0/MY_PIS*alpha*erfcd/cut_coul);
215 e_shift = erfcc/cut_coul - f_shift*cut_coul;
216 }
217
218 /* ----------------------------------------------------------------------
219 init for one type pair i,j and corresponding j,i
220 ------------------------------------------------------------------------- */
221
init_one(int,int)222 double PairCoulDSF::init_one(int /*i*/, int /*j*/)
223 {
224 return cut_coul;
225 }
226
227 /* ----------------------------------------------------------------------
228 proc 0 writes to restart file
229 ------------------------------------------------------------------------- */
230
write_restart(FILE * fp)231 void PairCoulDSF::write_restart(FILE *fp)
232 {
233 write_restart_settings(fp);
234
235 int i,j;
236 for (i = 1; i <= atom->ntypes; i++)
237 for (j = i; j <= atom->ntypes; j++) {
238 fwrite(&setflag[i][j],sizeof(int),1,fp);
239 }
240 }
241
242 /* ----------------------------------------------------------------------
243 proc 0 reads from restart file, bcasts
244 ------------------------------------------------------------------------- */
245
read_restart(FILE * fp)246 void PairCoulDSF::read_restart(FILE *fp)
247 {
248 read_restart_settings(fp);
249 allocate();
250
251 int i,j;
252 int me = comm->me;
253 for (i = 1; i <= atom->ntypes; i++)
254 for (j = i; j <= atom->ntypes; j++) {
255 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
256 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
257 }
258 }
259
260 /* ----------------------------------------------------------------------
261 proc 0 writes to restart file
262 ------------------------------------------------------------------------- */
263
write_restart_settings(FILE * fp)264 void PairCoulDSF::write_restart_settings(FILE *fp)
265 {
266 fwrite(&alpha,sizeof(double),1,fp);
267 fwrite(&cut_coul,sizeof(double),1,fp);
268 fwrite(&offset_flag,sizeof(int),1,fp);
269 fwrite(&mix_flag,sizeof(int),1,fp);
270 }
271
272 /* ----------------------------------------------------------------------
273 proc 0 reads from restart file, bcasts
274 ------------------------------------------------------------------------- */
275
read_restart_settings(FILE * fp)276 void PairCoulDSF::read_restart_settings(FILE *fp)
277 {
278 if (comm->me == 0) {
279 utils::sfread(FLERR,&alpha,sizeof(double),1,fp,nullptr,error);
280 utils::sfread(FLERR,&cut_coul,sizeof(double),1,fp,nullptr,error);
281 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
282 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
283 }
284 MPI_Bcast(&alpha,1,MPI_DOUBLE,0,world);
285 MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world);
286 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
287 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
288 }
289
290 /* ---------------------------------------------------------------------- */
291
single(int i,int j,int,int,double rsq,double factor_coul,double,double & fforce)292 double PairCoulDSF::single(int i, int j, int /*itype*/, int /*jtype*/, double rsq,
293 double factor_coul, double /*factor_lj*/,
294 double &fforce)
295 {
296 double r,erfcc,erfcd,prefactor,t;
297 double forcecoul,phicoul;
298
299 forcecoul = phicoul = 0.0;
300 if (rsq < cut_coulsq) {
301 r = sqrt(rsq);
302 prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r;
303 erfcd = exp(-alpha*alpha*rsq);
304 t = 1.0 / (1.0 + EWALD_P*alpha*r);
305 erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd;
306
307 forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS*erfcd +
308 r*f_shift) * r;
309 if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
310
311 phicoul = prefactor * (erfcc - r*e_shift - rsq*f_shift);
312 if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor;
313 }
314
315 fforce = forcecoul / rsq;
316
317 return phicoul;
318 }
319
320 /* ---------------------------------------------------------------------- */
321
extract(const char * str,int & dim)322 void *PairCoulDSF::extract(const char *str, int &dim)
323 {
324 if (strcmp(str,"cut_coul") == 0) {
325 dim = 0;
326 return (void *) &cut_coul;
327 }
328 return nullptr;
329 }
330