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