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