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