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 authors: Arben Jusufi, Axel Kohlmeyer (Temple U.)
17 ------------------------------------------------------------------------- */
18
19 #include "pair_gauss_cut.h"
20
21 #include <cmath>
22 #include "atom.h"
23 #include "comm.h"
24 #include "force.h"
25 #include "neigh_list.h"
26 #include "memory.h"
27 #include "error.h"
28
29 #include "math_const.h"
30
31 using namespace LAMMPS_NS;
32 using namespace MathConst;
33
34 /* ---------------------------------------------------------------------- */
35
PairGaussCut(LAMMPS * lmp)36 PairGaussCut::PairGaussCut(LAMMPS *lmp) : Pair(lmp)
37 {
38 respa_enable = 0;
39 writedata = 1;
40 }
41
42 /* ---------------------------------------------------------------------- */
43
~PairGaussCut()44 PairGaussCut::~PairGaussCut()
45 {
46 if (allocated) {
47 memory->destroy(setflag);
48 memory->destroy(cutsq);
49
50 memory->destroy(cut);
51 memory->destroy(hgauss);
52 memory->destroy(sigmah);
53 memory->destroy(rmh);
54 memory->destroy(pgauss);
55 memory->destroy(offset);
56 }
57 }
58
59 /* ---------------------------------------------------------------------- */
60
compute(int eflag,int vflag)61 void PairGaussCut::compute(int eflag, int vflag)
62 {
63 int i,j,ii,jj,inum,jnum,itype,jtype;
64 double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
65 double rsq,r,rexp,ugauss,factor_lj;
66 int *ilist,*jlist,*numneigh,**firstneigh;
67
68 evdwl = 0.0;
69 ev_init(eflag,vflag);
70
71 double **x = atom->x;
72 double **f = atom->f;
73 int *type = atom->type;
74 int nlocal = atom->nlocal;
75 double *special_lj = force->special_lj;
76 int newton_pair = force->newton_pair;
77
78 inum = list->inum;
79 ilist = list->ilist;
80 numneigh = list->numneigh;
81 firstneigh = list->firstneigh;
82
83 // loop over neighbors of my atoms
84
85 for (ii = 0; ii < inum; ii++) {
86 i = ilist[ii];
87 xtmp = x[i][0];
88 ytmp = x[i][1];
89 ztmp = x[i][2];
90 itype = type[i];
91 jlist = firstneigh[i];
92 jnum = numneigh[i];
93
94 for (jj = 0; jj < jnum; jj++) {
95 j = jlist[jj];
96 factor_lj = special_lj[sbmask(j)];
97 j &= NEIGHMASK;
98
99 delx = xtmp - x[j][0];
100 dely = ytmp - x[j][1];
101 delz = ztmp - x[j][2];
102 rsq = delx*delx + dely*dely + delz*delz;
103 jtype = type[j];
104
105 if (rsq < cutsq[itype][jtype]) {
106 r = sqrt(rsq);
107 rexp = (r-rmh[itype][jtype])/sigmah[itype][jtype];
108 ugauss = pgauss[itype][jtype]*exp(-0.5*rexp*rexp);
109 fpair = factor_lj*rexp/r*ugauss/sigmah[itype][jtype];
110
111 f[i][0] += delx*fpair;
112 f[i][1] += dely*fpair;
113 f[i][2] += delz*fpair;
114 if (newton_pair || j < nlocal) {
115 f[j][0] -= delx*fpair;
116 f[j][1] -= dely*fpair;
117 f[j][2] -= delz*fpair;
118 }
119
120 if (eflag) {
121 evdwl = ugauss - offset[itype][jtype];
122 evdwl *= factor_lj;
123 }
124
125 if (evflag) ev_tally(i,j,nlocal,newton_pair,
126 evdwl,0.0,fpair,delx,dely,delz);
127 }
128 }
129 }
130
131 if (vflag_fdotr) virial_fdotr_compute();
132 }
133
134
135 /* ----------------------------------------------------------------------
136 allocate all arrays
137 ------------------------------------------------------------------------- */
138
allocate()139 void PairGaussCut::allocate()
140 {
141 allocated = 1;
142 int n = atom->ntypes;
143
144 memory->create(setflag,n+1,n+1,"pair:setflag");
145 for (int i = 1; i <= n; i++)
146 for (int j = i; j <= n; j++)
147 setflag[i][j] = 0;
148
149 memory->create(cutsq,n+1,n+1,"pair:cutsq");
150
151 memory->create(cut,n+1,n+1,"pair:cut");
152 memory->create(hgauss,n+1,n+1,"pair:hgauss");
153 memory->create(sigmah,n+1,n+1,"pair:sigmah");
154 memory->create(rmh,n+1,n+1,"pair:rmh");
155 memory->create(pgauss,n+1,n+1,"pair:pgauss");
156 memory->create(offset,n+1,n+1,"pair:offset");
157 }
158
159 /* ----------------------------------------------------------------------
160 global settings
161 ------------------------------------------------------------------------- */
162
settings(int narg,char ** arg)163 void PairGaussCut::settings(int narg, char **arg)
164 {
165 if (narg != 1) error->all(FLERR,"Illegal pair_style command");
166
167 cut_global = utils::numeric(FLERR,arg[0],false,lmp);
168
169 // reset cutoffs that have been explicitly set
170
171 if (allocated) {
172 int i,j;
173 for (i = 1; i <= atom->ntypes; i++)
174 for (j = i; j <= atom->ntypes; j++)
175 if (setflag[i][j]) cut[i][j] = cut_global;
176 }
177 }
178
179 /* ----------------------------------------------------------------------
180 set coeffs for one or more type pairs
181 ------------------------------------------------------------------------- */
182
coeff(int narg,char ** arg)183 void PairGaussCut::coeff(int narg, char **arg)
184 {
185 if (narg < 5 || narg > 6) error->all(FLERR,"Incorrect args for pair coefficients");
186 if (!allocated) allocate();
187
188 int ilo,ihi,jlo,jhi;
189 utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
190 utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
191
192 double hgauss_one = utils::numeric(FLERR,arg[2],false,lmp);
193 double rmh_one = utils::numeric(FLERR,arg[3],false,lmp);
194 double sigmah_one = utils::numeric(FLERR,arg[4],false,lmp);
195 if (sigmah_one <= 0.0)
196 error->all(FLERR,"Incorrect args for pair coefficients");
197
198
199 double cut_one = cut_global;
200 if (narg == 6) cut_one = utils::numeric(FLERR,arg[5],false,lmp);
201
202 int count = 0;
203 for (int i = ilo; i <= ihi; i++) {
204 for (int j = MAX(jlo,i); j <= jhi; j++) {
205 hgauss[i][j] = hgauss_one;
206 sigmah[i][j] = sigmah_one;
207 rmh[i][j] = rmh_one;
208 cut[i][j] = cut_one;
209 setflag[i][j] = 1;
210 count++;
211 }
212 }
213
214 if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
215 }
216
217 /* ----------------------------------------------------------------------
218 init for one type pair i,j and corresponding j,i
219 ------------------------------------------------------------------------- */
220
init_one(int i,int j)221 double PairGaussCut::init_one(int i, int j)
222 {
223 if (setflag[i][j] == 0) {
224 hgauss[i][j] = mix_energy(fabs(hgauss[i][i]), fabs(hgauss[j][j]),
225 fabs(sigmah[i][i]), fabs(sigmah[j][j]));
226
227 // If either of the particles is repulsive (ie, if hgauss > 0),
228 // then the interaction between both is repulsive.
229 double sign_hi = (hgauss[i][i] >= 0.0) ? 1.0 : -1.0;
230 double sign_hj = (hgauss[j][j] >= 0.0) ? 1.0 : -1.0;
231 hgauss[i][j] *= MAX(sign_hi, sign_hj);
232
233 sigmah[i][j] = mix_distance(sigmah[i][i], sigmah[j][j]);
234 rmh[i][j] = mix_distance(rmh[i][i], rmh[j][j]);
235 cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
236 }
237
238 pgauss[i][j] = hgauss[i][j] / sqrt(MY_2PI) / sigmah[i][j];
239
240 if (offset_flag) {
241 double rexp = (cut[i][j]-rmh[i][j])/sigmah[i][j];
242 offset[i][j] = pgauss[i][j] * exp(-0.5*rexp*rexp);
243 } else offset[i][j] = 0.0;
244
245 hgauss[j][i] = hgauss[i][j];
246 sigmah[j][i] = sigmah[i][j];
247 rmh[j][i] = rmh[i][j];
248 pgauss[j][i] = pgauss[i][j];
249 offset[j][i] = offset[i][j];
250 cut[j][i] = cut[i][j];
251
252 // compute I,J contribution to long-range tail correction
253 // count total # of atoms of type I and J via Allreduce
254
255 if (tail_flag) {
256 int *type = atom->type;
257 int nlocal = atom->nlocal;
258
259 double count[2],all[2];
260 count[0] = count[1] = 0.0;
261 for (int k = 0; k < nlocal; k++) {
262 if (type[k] == i) count[0] += 1.0;
263 if (type[k] == j) count[1] += 1.0;
264 }
265 MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world);
266 }
267
268 return cut[i][j];
269 }
270
271 /* ----------------------------------------------------------------------
272 proc 0 writes to restart file
273 ------------------------------------------------------------------------- */
274
write_restart(FILE * fp)275 void PairGaussCut::write_restart(FILE *fp)
276 {
277 write_restart_settings(fp);
278
279 int i,j;
280 for (i = 1; i <= atom->ntypes; i++)
281 for (j = i; j <= atom->ntypes; j++) {
282 fwrite(&setflag[i][j],sizeof(int),1,fp);
283 if (setflag[i][j]) {
284 fwrite(&hgauss[i][j],sizeof(double),1,fp);
285 fwrite(&rmh[i][j],sizeof(double),1,fp);
286 fwrite(&sigmah[i][j],sizeof(double),1,fp);
287 fwrite(&cut[i][j],sizeof(double),1,fp);
288 }
289 }
290 }
291
292 /* ----------------------------------------------------------------------
293 proc 0 reads from restart file, bcasts
294 ------------------------------------------------------------------------- */
295
read_restart(FILE * fp)296 void PairGaussCut::read_restart(FILE *fp)
297 {
298 read_restart_settings(fp);
299 allocate();
300
301 int i,j;
302 int me = comm->me;
303 for (i = 1; i <= atom->ntypes; i++)
304 for (j = i; j <= atom->ntypes; j++) {
305 if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
306 MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
307 if (setflag[i][j]) {
308 if (me == 0) {
309 utils::sfread(FLERR,&hgauss[i][j],sizeof(double),1,fp,nullptr,error);
310 utils::sfread(FLERR,&rmh[i][j],sizeof(double),1,fp,nullptr,error);
311 utils::sfread(FLERR,&sigmah[i][j],sizeof(double),1,fp,nullptr,error);
312 utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
313 }
314 MPI_Bcast(&hgauss[i][j],1,MPI_DOUBLE,0,world);
315 MPI_Bcast(&rmh[i][j],1,MPI_DOUBLE,0,world);
316 MPI_Bcast(&sigmah[i][j],1,MPI_DOUBLE,0,world);
317 MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
318 }
319 }
320 }
321
322 /* ----------------------------------------------------------------------
323 proc 0 writes to restart file
324 ------------------------------------------------------------------------- */
325
write_restart_settings(FILE * fp)326 void PairGaussCut::write_restart_settings(FILE *fp)
327 {
328 fwrite(&cut_global,sizeof(double),1,fp);
329 fwrite(&offset_flag,sizeof(int),1,fp);
330 fwrite(&mix_flag,sizeof(int),1,fp);
331 }
332
333 /* ----------------------------------------------------------------------
334 proc 0 reads from restart file, bcasts
335 ------------------------------------------------------------------------- */
336
read_restart_settings(FILE * fp)337 void PairGaussCut::read_restart_settings(FILE *fp)
338 {
339 int me = comm->me;
340 if (me == 0) {
341 utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
342 utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error);
343 utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
344 }
345 MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
346 MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
347 MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
348 }
349
350 /* ----------------------------------------------------------------------
351 proc 0 writes to data file
352 ------------------------------------------------------------------------- */
353
write_data(FILE * fp)354 void PairGaussCut::write_data(FILE *fp)
355 {
356 for (int i = 1; i <= atom->ntypes; i++)
357 fprintf(fp,"%d %g %g %g\n",i,hgauss[i][i],rmh[i][i],sigmah[i][i]);
358 }
359
360 /* ----------------------------------------------------------------------
361 proc 0 writes all pairs to data file
362 ------------------------------------------------------------------------- */
363
write_data_all(FILE * fp)364 void PairGaussCut::write_data_all(FILE *fp)
365 {
366 for (int i = 1; i <= atom->ntypes; i++)
367 for (int j = i; j <= atom->ntypes; j++)
368 fprintf(fp,"%d %d %g %g %g %g\n",i,j,hgauss[i][j],rmh[i][j],sigmah[i][j],cut[i][j]);
369 }
370
371 /* ---------------------------------------------------------------------- */
372
single(int,int,int itype,int jtype,double rsq,double,double factor_lj,double & fforce)373 double PairGaussCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
374 double /*factor_coul*/, double factor_lj,
375 double &fforce)
376 {
377 double r, rexp,ugauss,phigauss;
378
379 r=sqrt(rsq);
380 rexp = (r-rmh[itype][jtype])/sigmah[itype][jtype];
381 ugauss = pgauss[itype][jtype]*exp(-0.5*rexp*rexp);
382
383 fforce = factor_lj*rexp/r*ugauss/sigmah[itype][jtype];
384
385 phigauss = ugauss - offset[itype][jtype];
386 return factor_lj*phigauss;
387 }
388
389 /* ---------------------------------------------------------------------- */
memory_usage()390 double PairGaussCut::memory_usage()
391 {
392 const int n=atom->ntypes;
393
394 double bytes = Pair::memory_usage();
395
396 bytes += 7.0*((n+1.0)*(n+1.0) * sizeof(double) + (n+1.0)*sizeof(double *));
397 bytes += 1.0*((n+1.0)*(n+1.0) * sizeof(int) + (n+1.0)*sizeof(int *));
398
399 return bytes;
400 }
401