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