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 ------------------------------------------------------------------------- */
18 
19 #include "pair_dpd_ext_tstat.h"
20 
21 #include "atom.h"
22 #include "comm.h"
23 #include "error.h"
24 #include "force.h"
25 #include "neigh_list.h"
26 #include "random_mars.h"
27 #include "update.h"
28 
29 #include <cmath>
30 
31 using namespace LAMMPS_NS;
32 
33 #define EPSILON 1.0e-10
34 
35 /* ---------------------------------------------------------------------- */
36 
PairDPDExtTstat(LAMMPS * lmp)37 PairDPDExtTstat::PairDPDExtTstat(LAMMPS *lmp) : PairDPDExt(lmp)
38 {
39   single_enable = 0;
40   writedata = 1;
41 }
42 
43 /* ---------------------------------------------------------------------- */
44 
compute(int eflag,int vflag)45 void PairDPDExtTstat::compute(int eflag, int vflag)
46 {
47   int i,j,ii,jj,inum,jnum,itype,jtype;
48   double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
49   double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
50   double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd;
51   double P[3][3];
52   int *ilist,*jlist,*numneigh,**firstneigh;
53 
54   ev_init(eflag,vflag);
55 
56   // adjust sigma if target T is changing
57 
58   if (t_start != t_stop) {
59     double delta = update->ntimestep - update->beginstep;
60     if (delta != 0.0) delta /= update->endstep - update->beginstep;
61     temperature = t_start + delta * (t_stop-t_start);
62     double boltz = force->boltz;
63     for (i = 1; i <= atom->ntypes; i++) {
64       for (j = i; j <= atom->ntypes; j++) {
65         sigma[i][j] = sigma[j][i] = sqrt(2.0*boltz*temperature*gamma[i][j]);
66         sigmaT[i][j] = sigmaT[j][i] = sqrt(2.0*boltz*temperature*gammaT[i][j]);
67       }
68     }
69   }
70 
71   double **x = atom->x;
72   double **v = atom->v;
73   double **f = atom->f;
74   int *type = atom->type;
75   int nlocal = atom->nlocal;
76   double *special_lj = force->special_lj;
77   int newton_pair = force->newton_pair;
78   double dtinvsqrt = 1.0/sqrt(update->dt);
79 
80   inum = list->inum;
81   ilist = list->ilist;
82   numneigh = list->numneigh;
83   firstneigh = list->firstneigh;
84 
85   // loop over neighbors of my atoms
86 
87   for (ii = 0; ii < inum; ii++) {
88     i = ilist[ii];
89     xtmp = x[i][0];
90     ytmp = x[i][1];
91     ztmp = x[i][2];
92     vxtmp = v[i][0];
93     vytmp = v[i][1];
94     vztmp = v[i][2];
95     itype = type[i];
96     jlist = firstneigh[i];
97     jnum = numneigh[i];
98 
99     for (jj = 0; jj < jnum; jj++) {
100       j = jlist[jj];
101       factor_dpd = special_lj[sbmask(j)];
102       j &= NEIGHMASK;
103 
104       delx = xtmp - x[j][0];
105       dely = ytmp - x[j][1];
106       delz = ztmp - x[j][2];
107       rsq = delx*delx + dely*dely + delz*delz;
108       jtype = type[j];
109 
110       if (rsq < cutsq[itype][jtype]) {
111         r = sqrt(rsq);
112         if (r < EPSILON) continue;     // r can be 0.0 in DPD systems
113         rinv = 1.0/r;
114         delvx = vxtmp - v[j][0];
115         delvy = vytmp - v[j][1];
116         delvz = vztmp - v[j][2];
117         dot = delx*delvx + dely*delvy + delz*delvz;
118 
119         P[0][0] = 1.0 - delx*delx*rinv*rinv;
120         P[0][1] =     - delx*dely*rinv*rinv;
121         P[0][2] =     - delx*delz*rinv*rinv;
122 
123         P[1][0] = P[0][1];
124         P[1][1] = 1.0 - dely*dely*rinv*rinv;
125         P[1][2] =     - dely*delz*rinv*rinv;
126 
127         P[2][0] = P[0][2];
128         P[2][1] = P[1][2];
129         P[2][2] = 1.0 - delz*delz*rinv*rinv;
130 
131         wd = 1.0 - r/cut[itype][jtype];
132         wdPar = pow(wd,ws[itype][jtype]);
133         wdPerp = pow(wd,wsT[itype][jtype]);
134 
135         randnum = random->gaussian();
136         randnumx = random->gaussian();
137         randnumy = random->gaussian();
138         randnumz = random->gaussian();
139 
140         // drag force - parallel
141         fpair = -gamma[itype][jtype]*wdPar*wdPar*dot*rinv;
142 
143         // random force - parallel
144         fpair += sigma[itype][jtype]*wdPar*randnum*dtinvsqrt;
145 
146         fpairx = fpair*rinv*delx;
147         fpairy = fpair*rinv*dely;
148         fpairz = fpair*rinv*delz;
149 
150         // drag force - perpendicular
151         fpairx -= gammaT[itype][jtype]*wdPerp*wdPerp*
152                   (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
153         fpairy -= gammaT[itype][jtype]*wdPerp*wdPerp*
154                   (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
155         fpairz -= gammaT[itype][jtype]*wdPerp*wdPerp*
156                   (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
157 
158         // random force - perpendicular
159         fpairx += sigmaT[itype][jtype]*wdPerp*
160                   (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
161         fpairy += sigmaT[itype][jtype]*wdPerp*
162                   (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
163         fpairz += sigmaT[itype][jtype]*wdPerp*
164                   (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
165 
166         fpairx *= factor_dpd;
167         fpairy *= factor_dpd;
168         fpairz *= factor_dpd;
169 
170         f[i][0] += fpairx;
171         f[i][1] += fpairy;
172         f[i][2] += fpairz;
173         if (newton_pair || j < nlocal) {
174           f[j][0] -= fpairx;
175           f[j][1] -= fpairy;
176           f[j][2] -= fpairz;
177         }
178 
179         if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
180                       0.0,0.0,
181                       fpairx, fpairy, fpairz,
182                       delx,dely,delz);
183       }
184     }
185   }
186 
187   if (vflag_fdotr) virial_fdotr_compute();
188 }
189 
190 /* ----------------------------------------------------------------------
191    global settings
192 ------------------------------------------------------------------------- */
193 
settings(int narg,char ** arg)194 void PairDPDExtTstat::settings(int narg, char **arg)
195 {
196   if (narg != 4) error->all(FLERR,"Illegal pair_style command");
197 
198   t_start = utils::numeric(FLERR,arg[0],false,lmp);
199   t_stop = utils::numeric(FLERR,arg[1],false,lmp);
200   cut_global = utils::numeric(FLERR,arg[2],false,lmp);
201   seed = utils::inumeric(FLERR,arg[3],false,lmp);
202 
203   temperature = t_start;
204 
205   // initialize Marsaglia RNG with processor-unique seed
206 
207   if (seed <= 0) error->all(FLERR,"Illegal pair_style command");
208   delete random;
209   random = new RanMars(lmp,seed + comm->me);
210 
211   // reset cutoffs that have been explicitly set
212 
213   if (allocated) {
214     int i,j;
215     for (i = 1; i <= atom->ntypes; i++)
216       for (j = i; j <= atom->ntypes; j++)
217         if (setflag[i][j]) cut[i][j] = cut_global;
218   }
219 }
220 
221 /* ----------------------------------------------------------------------
222    set coeffs for one or more type pairs
223 ------------------------------------------------------------------------- */
224 
coeff(int narg,char ** arg)225 void PairDPDExtTstat::coeff(int narg, char **arg)
226 {
227   if (narg < 6 || narg > 7)
228     error->all(FLERR,"Incorrect args for pair coefficients");
229   if (!allocated) allocate();
230 
231   int ilo,ihi,jlo,jhi;
232   utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
233   utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
234 
235   double a0_one = 0.0;
236   double gamma_one = utils::numeric(FLERR,arg[2],false,lmp);
237   double gammaT_one = utils::numeric(FLERR,arg[3],false,lmp);
238   double ws_one = utils::numeric(FLERR,arg[4],false,lmp);
239   double wsT_one = utils::numeric(FLERR,arg[5],false,lmp);
240 
241   double cut_one = cut_global;
242   if (narg == 7) cut_one = utils::numeric(FLERR,arg[6],false,lmp);
243 
244   int count = 0;
245   for (int i = ilo; i <= ihi; i++) {
246     for (int j = MAX(jlo,i); j <= jhi; j++) {
247       a0[i][j] = a0_one;
248       gamma[i][j] = gamma_one;
249       gammaT[i][j] = gammaT_one;
250       ws[i][j] = ws_one;
251       wsT[i][j] = wsT_one;
252       cut[i][j] = cut_one;
253       setflag[i][j] = 1;
254       count++;
255     }
256   }
257 
258   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
259 }
260 
261 /* ----------------------------------------------------------------------
262    proc 0 writes to restart file
263 ------------------------------------------------------------------------- */
264 
write_restart(FILE * fp)265 void PairDPDExtTstat::write_restart(FILE *fp)
266 {
267   write_restart_settings(fp);
268 
269   int i,j;
270   for (i = 1; i <= atom->ntypes; i++)
271     for (j = i; j <= atom->ntypes; j++) {
272       fwrite(&setflag[i][j],sizeof(int),1,fp);
273       if (setflag[i][j]) {
274         fwrite(&gamma[i][j],sizeof(double),1,fp);
275         fwrite(&gammaT[i][j],sizeof(double),1,fp);
276         fwrite(&ws[i][j],sizeof(double),1,fp);
277         fwrite(&wsT[i][j],sizeof(double),1,fp);
278         fwrite(&cut[i][j],sizeof(double),1,fp);
279       }
280     }
281 }
282 
283 /* ----------------------------------------------------------------------
284    proc 0 reads from restart file, bcasts
285 ------------------------------------------------------------------------- */
286 
read_restart(FILE * fp)287 void PairDPDExtTstat::read_restart(FILE *fp)
288 {
289   read_restart_settings(fp);
290 
291   allocate();
292 
293   int i,j;
294   int me = comm->me;
295   for (i = 1; i <= atom->ntypes; i++)
296     for (j = i; j <= atom->ntypes; j++) {
297       if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error);
298       MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
299       if (setflag[i][j]) {
300         if (me == 0) {
301           utils::sfread(FLERR,&gamma[i][j],sizeof(double),1,fp,nullptr,error);
302           utils::sfread(FLERR,&gammaT[i][j],sizeof(double),1,fp,nullptr,error);
303           utils::sfread(FLERR,&ws[i][j],sizeof(double),1,fp,nullptr,error);
304           utils::sfread(FLERR,&wsT[i][j],sizeof(double),1,fp,nullptr,error);
305           utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error);
306         }
307         MPI_Bcast(&gamma[i][j],1,MPI_DOUBLE,0,world);
308         MPI_Bcast(&gammaT[i][j],1,MPI_DOUBLE,0,world);
309         MPI_Bcast(&ws[i][j],1,MPI_DOUBLE,0,world);
310         MPI_Bcast(&wsT[i][j],1,MPI_DOUBLE,0,world);
311         MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
312       }
313     }
314 }
315 
316 /* ----------------------------------------------------------------------
317    proc 0 writes to restart file
318 ------------------------------------------------------------------------- */
319 
write_restart_settings(FILE * fp)320 void PairDPDExtTstat::write_restart_settings(FILE *fp)
321 {
322   fwrite(&t_start,sizeof(double),1,fp);
323   fwrite(&t_stop,sizeof(double),1,fp);
324   fwrite(&cut_global,sizeof(double),1,fp);
325   fwrite(&seed,sizeof(int),1,fp);
326   fwrite(&mix_flag,sizeof(int),1,fp);
327 }
328 
329 /* ----------------------------------------------------------------------
330    proc 0 reads from restart file, bcasts
331 ------------------------------------------------------------------------- */
332 
read_restart_settings(FILE * fp)333 void PairDPDExtTstat::read_restart_settings(FILE *fp)
334 {
335   if (comm->me == 0) {
336     utils::sfread(FLERR,&t_start,sizeof(double),1,fp,nullptr,error);
337     utils::sfread(FLERR,&t_stop,sizeof(double),1,fp,nullptr,error);
338     utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error);
339     utils::sfread(FLERR,&seed,sizeof(int),1,fp,nullptr,error);
340     utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error);
341   }
342   MPI_Bcast(&t_start,1,MPI_DOUBLE,0,world);
343   MPI_Bcast(&t_stop,1,MPI_DOUBLE,0,world);
344   MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
345   MPI_Bcast(&seed,1,MPI_INT,0,world);
346   MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
347 
348   temperature = t_start;
349 
350   // initialize Marsaglia RNG with processor-unique seed
351   // same seed that pair_style command initially specified
352 
353   if (random) delete random;
354   random = new RanMars(lmp,seed + comm->me);
355 }
356 
357 /* ----------------------------------------------------------------------
358    proc 0 writes to data file
359 ------------------------------------------------------------------------- */
360 
write_data(FILE * fp)361 void PairDPDExtTstat::write_data(FILE *fp)
362 {
363   for (int i = 1; i <= atom->ntypes; i++)
364     fprintf(fp,"%d %g %g %g %g\n",i,gamma[i][i],gammaT[i][i],ws[i][i],wsT[i][i]);
365 }
366 
367 /* ----------------------------------------------------------------------
368    proc 0 writes all pairs to data file
369 ------------------------------------------------------------------------- */
370 
write_data_all(FILE * fp)371 void PairDPDExtTstat::write_data_all(FILE *fp)
372 {
373   for (int i = 1; i <= atom->ntypes; i++)
374     for (int j = i; j <= atom->ntypes; j++)
375       fprintf(fp,"%d %d %g %g %g %g %g\n",i,j,
376         gamma[i][j],gammaT[i][j],ws[i][j],wsT[i][j],cut[i][j]);
377 }
378