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: Martin Svoboda (ICPF, UJEP), Martin Lísal (ICPF, UJEP)
17   based on pair style dpd by: Kurt Smith (U Pittsburgh)
18 ------------------------------------------------------------------------- */
19 
20 #include "pair_dpd_ext.h"
21 
22 #include "atom.h"
23 #include "comm.h"
24 #include "error.h"
25 #include "force.h"
26 #include "memory.h"
27 #include "neigh_list.h"
28 #include "neighbor.h"
29 #include "random_mars.h"
30 #include "update.h"
31 
32 #include <cmath>
33 
34 using namespace LAMMPS_NS;
35 
36 #define EPSILON 1.0e-10
37 
38 /* ---------------------------------------------------------------------- */
39 
PairDPDExt(LAMMPS * lmp)40 PairDPDExt::PairDPDExt(LAMMPS *lmp) : Pair(lmp)
41 {
42   writedata = 1;
43   random = nullptr;
44 }
45 
46 /* ---------------------------------------------------------------------- */
47 
~PairDPDExt()48 PairDPDExt::~PairDPDExt()
49 {
50   if (allocated) {
51     memory->destroy(setflag);
52     memory->destroy(cutsq);
53 
54     memory->destroy(cut);
55     memory->destroy(a0);
56     memory->destroy(gamma);
57     memory->destroy(gammaT);
58     memory->destroy(sigma);
59     memory->destroy(sigmaT);
60     memory->destroy(ws);
61     memory->destroy(wsT);
62   }
63 
64   if (random) delete random;
65 }
66 
67 /* ---------------------------------------------------------------------- */
68 
compute(int eflag,int vflag)69 void PairDPDExt::compute(int eflag, int vflag)
70 {
71   int i,j,ii,jj,inum,jnum,itype,jtype;
72   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpairx,fpairy,fpairz,fpair;
73   double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
74   double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd;
75   double P[3][3];
76   int *ilist,*jlist,*numneigh,**firstneigh;
77 
78   evdwl = 0.0;
79   ev_init(eflag,vflag);
80 
81   double **x = atom->x;
82   double **v = atom->v;
83   double **f = atom->f;
84   int *type = atom->type;
85   int nlocal = atom->nlocal;
86   double *special_lj = force->special_lj;
87   int newton_pair = force->newton_pair;
88   double dtinvsqrt = 1.0/sqrt(update->dt);
89 
90   inum = list->inum;
91   ilist = list->ilist;
92   numneigh = list->numneigh;
93   firstneigh = list->firstneigh;
94 
95   // loop over neighbors of my atoms
96 
97   for (ii = 0; ii < inum; ii++) {
98     i = ilist[ii];
99     xtmp = x[i][0];
100     ytmp = x[i][1];
101     ztmp = x[i][2];
102     vxtmp = v[i][0];
103     vytmp = v[i][1];
104     vztmp = v[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_dpd = special_lj[sbmask(j)];
112       j &= NEIGHMASK;
113 
114       delx = xtmp - x[j][0];
115       dely = ytmp - x[j][1];
116       delz = ztmp - x[j][2];
117       rsq = delx*delx + dely*dely + delz*delz;
118       jtype = type[j];
119 
120       if (rsq < cutsq[itype][jtype]) {
121         r = sqrt(rsq);
122         if (r < EPSILON) continue;     // r can be 0.0 in DPD systems
123         rinv = 1.0/r;
124         delvx = vxtmp - v[j][0];
125         delvy = vytmp - v[j][1];
126         delvz = vztmp - v[j][2];
127         dot = delx*delvx + dely*delvy + delz*delvz;
128 
129         P[0][0] = 1.0 - delx*delx*rinv*rinv;
130         P[0][1] =     - delx*dely*rinv*rinv;
131         P[0][2] =     - delx*delz*rinv*rinv;
132 
133         P[1][0] = P[0][1];
134         P[1][1] = 1.0 - dely*dely*rinv*rinv;
135         P[1][2] =     - dely*delz*rinv*rinv;
136 
137         P[2][0] = P[0][2];
138         P[2][1] = P[1][2];
139         P[2][2] = 1.0 - delz*delz*rinv*rinv;
140 
141         wd = 1.0 - r/cut[itype][jtype];
142         wdPar = pow(wd,ws[itype][jtype]);
143         wdPerp = pow(wd,wsT[itype][jtype]);
144 
145         randnum = random->gaussian();
146         randnumx = random->gaussian();
147         randnumy = random->gaussian();
148         randnumz = random->gaussian();
149 
150         // conservative force
151         fpair = a0[itype][jtype]*wd;
152 
153         // drag force - parallel
154         fpair -= gamma[itype][jtype]*wdPar*wdPar*dot*rinv;
155 
156         // random force - parallel
157         fpair += sigma[itype][jtype]*wdPar*randnum*dtinvsqrt;
158 
159         fpairx = fpair*rinv*delx;
160         fpairy = fpair*rinv*dely;
161         fpairz = fpair*rinv*delz;
162 
163         // drag force - perpendicular
164         fpairx -= gammaT[itype][jtype]*wdPerp*wdPerp*
165                   (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
166         fpairy -= gammaT[itype][jtype]*wdPerp*wdPerp*
167                   (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
168         fpairz -= gammaT[itype][jtype]*wdPerp*wdPerp*
169                   (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
170 
171         // random force - perpendicular
172         fpairx += sigmaT[itype][jtype]*wdPerp*
173                   (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
174         fpairy += sigmaT[itype][jtype]*wdPerp*
175                   (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
176         fpairz += sigmaT[itype][jtype]*wdPerp*
177                   (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
178 
179         fpairx *= factor_dpd;
180         fpairy *= factor_dpd;
181         fpairz *= factor_dpd;
182 
183         f[i][0] += fpairx;
184         f[i][1] += fpairy;
185         f[i][2] += fpairz;
186         if (newton_pair || j < nlocal) {
187           f[j][0] -= fpairx;
188           f[j][1] -= fpairy;
189           f[j][2] -= fpairz;
190         }
191 
192         if (eflag) {
193           // unshifted eng of conservative term:
194           // evdwl = -a0[itype][jtype]*r * (1.0-0.5*r/cut[itype][jtype]);
195           // eng shifted to 0.0 at cutoff
196           evdwl = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd*wd;
197           evdwl *= factor_dpd;
198         }
199         if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
200                       evdwl,0.0,
201                       fpairx, fpairy, fpairz,
202                       delx,dely,delz);
203       }
204     }
205   }
206 
207   if (vflag_fdotr) virial_fdotr_compute();
208 }
209 
210 /* ----------------------------------------------------------------------
211    allocate all arrays
212 ------------------------------------------------------------------------- */
213 
allocate()214 void PairDPDExt::allocate()
215 {
216   int i,j;
217   allocated = 1;
218   int n = atom->ntypes;
219 
220   memory->create(setflag,n+1,n+1,"pair:setflag");
221   for (i = 1; i <= n; i++)
222     for (j = i; j <= n; j++)
223       setflag[i][j] = 0;
224 
225   memory->create(cutsq,n+1,n+1,"pair:cutsq");
226 
227   memory->create(cut,n+1,n+1,"pair:cut");
228   memory->create(a0,n+1,n+1,"pair:a0");
229   memory->create(gamma,n+1,n+1,"pair:gamma");
230   memory->create(gammaT,n+1,n+1,"pair:gammaT");
231   memory->create(sigma,n+1,n+1,"pair:sigma");
232   memory->create(sigmaT,n+1,n+1,"pair:sigmaT");
233   memory->create(ws,n+1,n+1,"pair:ws");
234   memory->create(wsT,n+1,n+1,"pair:wsT");
235   for (i = 0; i <= atom->ntypes; i++)
236   {
237     for (j = 0; j <= atom->ntypes; j++)
238     {
239       sigma[i][j] = gamma[i][j] =sigmaT[i][j] = gammaT[i][j] = 0.0;
240       ws[i][j] = wsT[i][j] = 1.0;
241     }
242   }
243 }
244 
245 /* ----------------------------------------------------------------------
246    global settings
247 ------------------------------------------------------------------------- */
248 
settings(int narg,char ** arg)249 void PairDPDExt::settings(int narg, char **arg)
250 {
251   if (narg != 3) error->all(FLERR,"Illegal pair_style command");
252 
253   temperature = utils::numeric(FLERR,arg[0],false,lmp);
254   cut_global = utils::numeric(FLERR,arg[1],false,lmp);
255   seed = utils::inumeric(FLERR,arg[2],false,lmp);
256 
257   // initialize Marsaglia RNG with processor-unique seed
258 
259   if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
260   delete random;
261   random = new RanMars(lmp,seed + comm->me);
262 
263   // reset cutoffs that have been explicitly set
264 
265   if (allocated) {
266     int i,j;
267     for (i = 1; i <= atom->ntypes; i++)
268       for (j = i; j <= atom->ntypes; j++)
269         if (setflag[i][j]) {
270           cut[i][j] = cut_global;
271           cutsq[i][j] = cut_global*cut_global;
272         }
273   }
274 }
275 
276 /* ----------------------------------------------------------------------
277    set coeffs for one or more type pairs
278 ------------------------------------------------------------------------- */
279 
coeff(int narg,char ** arg)280 void PairDPDExt::coeff(int narg, char **arg)
281 {
282   if (narg < 7 || narg > 8)
283     error->all(FLERR,"Incorrect args for pair coefficients");
284   if (!allocated) allocate();
285 
286   int ilo,ihi,jlo,jhi;
287   utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
288   utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
289 
290   double a0_one = utils::numeric(FLERR,arg[2],false,lmp);
291   double gamma_one = utils::numeric(FLERR,arg[3],false,lmp);
292   double gammaT_one = utils::numeric(FLERR,arg[4],false,lmp);
293   double ws_one = utils::numeric(FLERR,arg[5],false,lmp);
294   double wsT_one = utils::numeric(FLERR,arg[6],false,lmp);
295 
296   double cut_one = cut_global;
297   if (narg == 8) cut_one = utils::numeric(FLERR,arg[7],false,lmp);
298 
299   int count = 0;
300   for (int i = ilo; i <= ihi; i++) {
301     for (int j = MAX(jlo,i); j <= jhi; j++) {
302       a0[i][j] = a0_one;
303       gamma[i][j] = gamma_one;
304       gammaT[i][j] = gammaT_one;
305       ws[i][j] = ws_one;
306       wsT[i][j] = wsT_one;
307       cut[i][j] = cut_one;
308       cutsq[i][j] = cut_one*cut_one;
309       setflag[i][j] = 1;
310       count++;
311     }
312   }
313 
314   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
315 }
316 
317 /* ----------------------------------------------------------------------
318    init specific to this pair style
319 ------------------------------------------------------------------------- */
320 
init_style()321 void PairDPDExt::init_style()
322 {
323   if (comm->ghost_velocity == 0)
324     error->all(FLERR,"Pair dpd requires ghost atoms store velocity");
325 
326   // if newton off, forces between atoms ij will be double computed
327   // using different random numbers
328 
329   if (force->newton_pair == 0 && comm->me == 0) error->warning(FLERR,
330       "Pair dpd needs newton pair on for momentum conservation");
331 
332   neighbor->request(this,instance_me);
333 }
334 
335 /* ----------------------------------------------------------------------
336    init for one type pair i,j and corresponding j,i
337 ------------------------------------------------------------------------- */
338 
init_one(int i,int j)339 double PairDPDExt::init_one(int i, int j)
340 {
341   if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
342 
343   sigma[i][j] = sqrt(2.0*force->boltz*temperature*gamma[i][j]);
344   sigmaT[i][j] = sqrt(2.0*force->boltz*temperature*gammaT[i][j]);
345 
346   cut[j][i] = cut[i][j];
347   cutsq[j][i] = cutsq[i][j];
348   a0[j][i] = a0[i][j];
349   gamma[j][i] = gamma[i][j];
350   gammaT[j][i] = gammaT[i][j];
351   sigma[j][i] = sigma[i][j];
352   sigmaT[j][i] = sigmaT[i][j];
353   ws[j][i] = ws[i][j];
354   wsT[j][i] = wsT[i][j];
355 
356   return cut[i][j];
357 }
358 
359 /* ----------------------------------------------------------------------
360    proc 0 writes to restart file
361 ------------------------------------------------------------------------- */
362 
write_restart(FILE * fp)363 void PairDPDExt::write_restart(FILE *fp)
364 {
365   write_restart_settings(fp);
366 
367   int i,j;
368   for (i = 1; i <= atom->ntypes; i++)
369     for (j = i; j <= atom->ntypes; j++) {
370       fwrite(&setflag[i][j],sizeof(int),1,fp);
371       if (setflag[i][j]) {
372         fwrite(&a0[i][j],sizeof(double),1,fp);
373         fwrite(&gamma[i][j],sizeof(double),1,fp);
374         fwrite(&gammaT[i][j],sizeof(double),1,fp);
375         fwrite(&ws[i][j],sizeof(double),1,fp);
376         fwrite(&wsT[i][j],sizeof(double),1,fp);
377         fwrite(&cut[i][j],sizeof(double),1,fp);
378       }
379     }
380 }
381 
382 /* ----------------------------------------------------------------------
383    proc 0 reads from restart file, bcasts
384 ------------------------------------------------------------------------- */
385 
read_restart(FILE * fp)386 void PairDPDExt::read_restart(FILE *fp)
387 {
388   read_restart_settings(fp);
389 
390   allocate();
391 
392   int i,j;
393   int me = comm->me;
394   for (i = 1; i <= atom->ntypes; i++)
395     for (j = i; j <= atom->ntypes; j++) {
396       if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
397       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
398       if (setflag[i][j]) {
399         if (me == 0) {
400           utils::sfread(FLERR,&a0[i][j],sizeof(double),1,fp,nullptr,error);
401           utils::sfread(FLERR,&gamma[i][j],sizeof(double),1,fp,nullptr,error);
402           utils::sfread(FLERR,&gammaT[i][j],sizeof(double),1,fp,nullptr,error);
403           utils::sfread(FLERR,&ws[i][j],sizeof(double),1,fp,nullptr,error);
404           utils::sfread(FLERR,&wsT[i][j],sizeof(double),1,fp,nullptr,error);
405           utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
406         }
407         MPI_Bcast(&a0[i][j],1,MPI_DOUBLE,0,world);
408         MPI_Bcast(&gamma[i][j],1,MPI_DOUBLE,0,world);
409         MPI_Bcast(&gammaT[i][j],1,MPI_DOUBLE,0,world);
410         MPI_Bcast(&ws[i][j],1,MPI_DOUBLE,0,world);
411         MPI_Bcast(&wsT[i][j],1,MPI_DOUBLE,0,world);
412         MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
413       }
414     }
415 }
416 
417 /* ----------------------------------------------------------------------
418    proc 0 writes to restart file
419 ------------------------------------------------------------------------- */
420 
write_restart_settings(FILE * fp)421 void PairDPDExt::write_restart_settings(FILE *fp)
422 {
423   fwrite(&temperature,sizeof(double),1,fp);
424   fwrite(&cut_global,sizeof(double),1,fp);
425   fwrite(&seed,sizeof(int),1,fp);
426   fwrite(&mix_flag,sizeof(int),1,fp);
427 }
428 
429 /* ----------------------------------------------------------------------
430    proc 0 reads from restart file, bcasts
431 ------------------------------------------------------------------------- */
432 
read_restart_settings(FILE * fp)433 void PairDPDExt::read_restart_settings(FILE *fp)
434 {
435   if (comm->me == 0) {
436     utils::sfread(FLERR,&temperature,sizeof(double),1,fp,nullptr,error);
437     utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
438     utils::sfread(FLERR,&seed,sizeof(int),1,fp,nullptr,error);
439     utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
440   }
441   MPI_Bcast(&temperature,1,MPI_DOUBLE,0,world);
442   MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
443   MPI_Bcast(&seed,1,MPI_INT,0,world);
444   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
445 
446   // initialize Marsaglia RNG with processor-unique seed
447   // same seed that pair_style command initially specified
448 
449   if (random) delete random;
450   random = new RanMars(lmp,seed + comm->me);
451 }
452 
453 /* ----------------------------------------------------------------------
454    proc 0 writes to data file
455 ------------------------------------------------------------------------- */
456 
write_data(FILE * fp)457 void PairDPDExt::write_data(FILE *fp)
458 {
459   for (int i = 1; i <= atom->ntypes; i++)
460     fprintf(fp,"%d %g %g %g %g %g\n",i,
461       a0[i][i],gamma[i][i],gammaT[i][i],ws[i][i],wsT[i][i]);
462 }
463 
464 /* ----------------------------------------------------------------------
465    proc 0 writes all pairs to data file
466 ------------------------------------------------------------------------- */
467 
write_data_all(FILE * fp)468 void PairDPDExt::write_data_all(FILE *fp)
469 {
470   for (int i = 1; i <= atom->ntypes; i++)
471     for (int j = i; j <= atom->ntypes; j++)
472       fprintf(fp,"%d %d %g %g %g %g %g %g\n",i,j,
473         a0[i][j],gamma[i][j],gammaT[i][j],ws[i][j],wsT[i][j],cut[i][j]);
474 }
475 
476 /* ---------------------------------------------------------------------- */
477 
single(int,int,int itype,int jtype,double rsq,double,double factor_dpd,double & fforce)478 double PairDPDExt::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
479                        double /*factor_coul*/, double factor_dpd, double &fforce)
480 {
481   double r,rinv,wd,phi;
482 
483   r = sqrt(rsq);
484   if (r < EPSILON) {
485     fforce = 0.0;
486     return 0.0;
487   }
488 
489   rinv = 1.0/r;
490   wd = 1.0 - r/cut[itype][jtype];
491   fforce = a0[itype][jtype]*wd * factor_dpd*rinv;
492 
493   phi = 0.5*a0[itype][jtype]*cut[itype][jtype] * wd*wd;
494   return factor_dpd*phi;
495 }
496