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 "fix_press_berendsen.h"
16
17 #include "atom.h"
18 #include "comm.h"
19 #include "compute.h"
20 #include "domain.h"
21 #include "error.h"
22 #include "fix_deform.h"
23 #include "force.h"
24 #include "kspace.h"
25 #include "modify.h"
26 #include "update.h"
27
28 #include <cmath>
29 #include <cstring>
30
31 using namespace LAMMPS_NS;
32 using namespace FixConst;
33
34 enum{NOBIAS,BIAS};
35 enum{NONE,XYZ,XY,YZ,XZ};
36 enum{ISO,ANISO};
37
38 /* ---------------------------------------------------------------------- */
39
FixPressBerendsen(LAMMPS * lmp,int narg,char ** arg)40 FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) :
41 Fix(lmp, narg, arg),
42 id_temp(nullptr), id_press(nullptr), tflag(0), pflag(0)
43 {
44 if (narg < 5) error->all(FLERR,"Illegal fix press/berendsen command");
45
46 // Berendsen barostat applied every step
47
48 nevery = 1;
49
50 // default values
51
52 pcouple = NONE;
53 bulkmodulus = 10.0;
54 allremap = 1;
55
56 for (int i = 0; i < 3; i++) {
57 p_start[i] = p_stop[i] = p_period[i] = 0.0;
58 p_flag[i] = 0;
59 p_period[i] = 0.0;
60 }
61
62 // process keywords
63
64 dimension = domain->dimension;
65
66 int iarg = 3;
67
68 while (iarg < narg) {
69 if (strcmp(arg[iarg],"iso") == 0) {
70 if (iarg+4 > narg)
71 error->all(FLERR,"Illegal fix press/berendsen command");
72 pcouple = XYZ;
73 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
74 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
75 p_period[0] = p_period[1] = p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
76 p_flag[0] = p_flag[1] = p_flag[2] = 1;
77 if (dimension == 2) {
78 p_start[2] = p_stop[2] = p_period[2] = 0.0;
79 p_flag[2] = 0;
80 }
81 iarg += 4;
82 } else if (strcmp(arg[iarg],"aniso") == 0) {
83 if (iarg+4 > narg)
84 error->all(FLERR,"Illegal fix press/berendsen command");
85 pcouple = NONE;
86 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
87 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
88 p_period[0] = p_period[1] = p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
89 p_flag[0] = p_flag[1] = p_flag[2] = 1;
90 if (dimension == 2) {
91 p_start[2] = p_stop[2] = p_period[2] = 0.0;
92 p_flag[2] = 0;
93 }
94 iarg += 4;
95
96 } else if (strcmp(arg[iarg],"x") == 0) {
97 if (iarg+4 > narg)
98 error->all(FLERR,"Illegal fix press/berendsen command");
99 p_start[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
100 p_stop[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
101 p_period[0] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
102 p_flag[0] = 1;
103 iarg += 4;
104 } else if (strcmp(arg[iarg],"y") == 0) {
105 if (iarg+4 > narg)
106 error->all(FLERR,"Illegal fix press/berendsen command");
107 p_start[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
108 p_stop[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
109 p_period[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
110 p_flag[1] = 1;
111 iarg += 4;
112 } else if (strcmp(arg[iarg],"z") == 0) {
113 if (iarg+4 > narg)
114 error->all(FLERR,"Illegal fix press/berendsen command");
115 p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
116 p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
117 p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
118 p_flag[2] = 1;
119 iarg += 4;
120 if (dimension == 2)
121 error->all(FLERR,"Invalid fix press/berendsen for a 2d simulation");
122
123 } else if (strcmp(arg[iarg],"couple") == 0) {
124 if (iarg+2 > narg)
125 error->all(FLERR,"Illegal fix press/berendsen command");
126 if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
127 else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
128 else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
129 else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
130 else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
131 else error->all(FLERR,"Illegal fix press/berendsen command");
132 iarg += 2;
133
134 } else if (strcmp(arg[iarg],"modulus") == 0) {
135 if (iarg+2 > narg)
136 error->all(FLERR,"Illegal fix press/berendsen command");
137 bulkmodulus = utils::numeric(FLERR,arg[iarg+1],false,lmp);
138 if (bulkmodulus <= 0.0)
139 error->all(FLERR,"Illegal fix press/berendsen command");
140 iarg += 2;
141 } else if (strcmp(arg[iarg],"dilate") == 0) {
142 if (iarg+2 > narg)
143 error->all(FLERR,"Illegal fix press/berendsen command");
144 if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
145 else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
146 else error->all(FLERR,"Illegal fix press/berendsen command");
147 iarg += 2;
148 } else error->all(FLERR,"Illegal fix press/berendsen command");
149 }
150
151 if (allremap == 0) restart_pbc = 1;
152
153 // error checks
154
155 if (dimension == 2 && p_flag[2])
156 error->all(FLERR,"Invalid fix press/berendsen for a 2d simulation");
157 if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
158 error->all(FLERR,"Invalid fix press/berendsen for a 2d simulation");
159
160 if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
161 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
162 if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
163 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
164 if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
165 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
166 if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
167 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
168 if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
169 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
170
171 if (p_flag[0] && domain->xperiodic == 0)
172 error->all(FLERR,
173 "Cannot use fix press/berendsen on a non-periodic dimension");
174 if (p_flag[1] && domain->yperiodic == 0)
175 error->all(FLERR,
176 "Cannot use fix press/berendsen on a non-periodic dimension");
177 if (p_flag[2] && domain->zperiodic == 0)
178 error->all(FLERR,
179 "Cannot use fix press/berendsen on a non-periodic dimension");
180
181 if (pcouple == XYZ && dimension == 3 &&
182 (p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
183 p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] ||
184 p_period[0] != p_period[1] || p_period[0] != p_period[2]))
185 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
186 if (pcouple == XYZ && dimension == 2 &&
187 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
188 p_period[0] != p_period[1]))
189 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
190 if (pcouple == XY &&
191 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
192 p_period[0] != p_period[1]))
193 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
194 if (pcouple == YZ &&
195 (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] ||
196 p_period[1] != p_period[2]))
197 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
198 if (pcouple == XZ &&
199 (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] ||
200 p_period[0] != p_period[2]))
201 error->all(FLERR,"Invalid fix press/berendsen pressure settings");
202
203 if ((p_flag[0] && p_period[0] <= 0.0) ||
204 (p_flag[1] && p_period[1] <= 0.0) ||
205 (p_flag[2] && p_period[2] <= 0.0))
206 error->all(FLERR,"Fix press/berendsen damping parameters must be > 0.0");
207
208 if (p_flag[0]) box_change |= BOX_CHANGE_X;
209 if (p_flag[1]) box_change |= BOX_CHANGE_Y;
210 if (p_flag[2]) box_change |= BOX_CHANGE_Z;
211
212 // pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof
213 // else pstyle = ANISO -> 3 dof
214
215 if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
216 else pstyle = ANISO;
217
218 // create a new compute temp style
219 // id = fix-ID + temp
220 // compute group = all since pressure is always global (group all)
221 // and thus its KE/temperature contribution should use group all
222
223 id_temp = utils::strdup(std::string(id) + "_temp");
224 modify->add_compute(fmt::format("{} all temp",id_temp));
225 tflag = 1;
226
227 // create a new compute pressure style
228 // id = fix-ID + press, compute group = all
229 // pass id_temp as 4th arg to pressure constructor
230
231 id_press = utils::strdup(std::string(id) + "_press");
232 modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
233 pflag = 1;
234
235 nrigid = 0;
236 rfix = nullptr;
237 }
238
239 /* ---------------------------------------------------------------------- */
240
~FixPressBerendsen()241 FixPressBerendsen::~FixPressBerendsen()
242 {
243 delete [] rfix;
244
245 // delete temperature and pressure if fix created them
246
247 if (tflag) modify->delete_compute(id_temp);
248 if (pflag) modify->delete_compute(id_press);
249 delete [] id_temp;
250 delete [] id_press;
251 }
252
253 /* ---------------------------------------------------------------------- */
254
setmask()255 int FixPressBerendsen::setmask()
256 {
257 int mask = 0;
258 mask |= END_OF_STEP;
259 return mask;
260 }
261
262 /* ---------------------------------------------------------------------- */
263
init()264 void FixPressBerendsen::init()
265 {
266 if (domain->triclinic)
267 error->all(FLERR,"Cannot use fix press/berendsen with triclinic box");
268
269 // insure no conflict with fix deform
270
271 for (int i = 0; i < modify->nfix; i++)
272 if (strcmp(modify->fix[i]->style,"deform") == 0) {
273 int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
274 if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
275 (p_flag[2] && dimflag[2]))
276 error->all(FLERR,"Cannot use fix press/berendsen and "
277 "fix deform on same component of stress tensor");
278 }
279
280 // set temperature and pressure ptrs
281
282 int icompute = modify->find_compute(id_temp);
283 if (icompute < 0)
284 error->all(FLERR,"Temperature ID for fix press/berendsen does not exist");
285 temperature = modify->compute[icompute];
286
287 if (temperature->tempbias) which = BIAS;
288 else which = NOBIAS;
289
290 icompute = modify->find_compute(id_press);
291 if (icompute < 0)
292 error->all(FLERR,"Pressure ID for fix press/berendsen does not exist");
293 pressure = modify->compute[icompute];
294
295 // Kspace setting
296
297 if (force->kspace) kspace_flag = 1;
298 else kspace_flag = 0;
299
300 // detect if any rigid fixes exist so rigid bodies move when box is remapped
301 // rfix[] = indices to each fix rigid
302
303 delete [] rfix;
304 nrigid = 0;
305 rfix = nullptr;
306
307 for (int i = 0; i < modify->nfix; i++)
308 if (modify->fix[i]->rigid_flag) nrigid++;
309 if (nrigid) {
310 rfix = new int[nrigid];
311 nrigid = 0;
312 for (int i = 0; i < modify->nfix; i++)
313 if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
314 }
315 }
316
317 /* ----------------------------------------------------------------------
318 compute T,P before integrator starts
319 ------------------------------------------------------------------------- */
320
setup(int)321 void FixPressBerendsen::setup(int /*vflag*/)
322 {
323 // trigger virial computation on next timestep
324
325 pressure->addstep(update->ntimestep+1);
326 }
327
328 /* ---------------------------------------------------------------------- */
329
end_of_step()330 void FixPressBerendsen::end_of_step()
331 {
332 // compute new T,P
333
334 if (pstyle == ISO) {
335 temperature->compute_scalar();
336 pressure->compute_scalar();
337 } else {
338 temperature->compute_vector();
339 pressure->compute_vector();
340 }
341 couple();
342
343 double delta = update->ntimestep - update->beginstep;
344 if (delta != 0.0) delta /= update->endstep - update->beginstep;
345
346 for (int i = 0; i < 3; i++) {
347 if (p_flag[i]) {
348 p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
349 dilation[i] =
350 pow(1.0 - update->dt/p_period[i] *
351 (p_target[i]-p_current[i])/bulkmodulus,1.0/3.0);
352 }
353 }
354
355 // remap simulation box and atoms
356 // redo KSpace coeffs since volume has changed
357
358 remap();
359 if (kspace_flag) force->kspace->setup();
360
361 // trigger virial computation on next timestep
362
363 pressure->addstep(update->ntimestep+1);
364 }
365
366 /* ---------------------------------------------------------------------- */
367
couple()368 void FixPressBerendsen::couple()
369 {
370 double *tensor = pressure->vector;
371
372 if (pstyle == ISO)
373 p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
374 else if (pcouple == XYZ) {
375 double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
376 p_current[0] = p_current[1] = p_current[2] = ave;
377 } else if (pcouple == XY) {
378 double ave = 0.5 * (tensor[0] + tensor[1]);
379 p_current[0] = p_current[1] = ave;
380 p_current[2] = tensor[2];
381 } else if (pcouple == YZ) {
382 double ave = 0.5 * (tensor[1] + tensor[2]);
383 p_current[1] = p_current[2] = ave;
384 p_current[0] = tensor[0];
385 } else if (pcouple == XZ) {
386 double ave = 0.5 * (tensor[0] + tensor[2]);
387 p_current[0] = p_current[2] = ave;
388 p_current[1] = tensor[1];
389 } else {
390 p_current[0] = tensor[0];
391 p_current[1] = tensor[1];
392 p_current[2] = tensor[2];
393 }
394 }
395
396 /* ----------------------------------------------------------------------
397 change box size
398 remap all atoms or fix group atoms depending on allremap flag
399 if rigid bodies exist, scale rigid body centers-of-mass
400 ------------------------------------------------------------------------- */
401
remap()402 void FixPressBerendsen::remap()
403 {
404 int i;
405 double oldlo,oldhi,ctr;
406
407 double **x = atom->x;
408 int *mask = atom->mask;
409 int nlocal = atom->nlocal;
410
411 // convert pertinent atoms and rigid bodies to lamda coords
412
413 if (allremap) domain->x2lamda(nlocal);
414 else {
415 for (i = 0; i < nlocal; i++)
416 if (mask[i] & groupbit)
417 domain->x2lamda(x[i],x[i]);
418 }
419
420 if (nrigid)
421 for (i = 0; i < nrigid; i++)
422 modify->fix[rfix[i]]->deform(0);
423
424 // reset global and local box to new size/shape
425
426 for (i = 0; i < 3; i++) {
427 if (p_flag[i]) {
428 oldlo = domain->boxlo[i];
429 oldhi = domain->boxhi[i];
430 ctr = 0.5 * (oldlo + oldhi);
431 domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
432 domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
433 }
434 }
435
436 domain->set_global_box();
437 domain->set_local_box();
438
439 // convert pertinent atoms and rigid bodies back to box coords
440
441 if (allremap) domain->lamda2x(nlocal);
442 else {
443 for (i = 0; i < nlocal; i++)
444 if (mask[i] & groupbit)
445 domain->lamda2x(x[i],x[i]);
446 }
447
448 if (nrigid)
449 for (i = 0; i < nrigid; i++)
450 modify->fix[rfix[i]]->deform(1);
451 }
452
453 /* ---------------------------------------------------------------------- */
454
modify_param(int narg,char ** arg)455 int FixPressBerendsen::modify_param(int narg, char **arg)
456 {
457 if (strcmp(arg[0],"temp") == 0) {
458 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
459 if (tflag) {
460 modify->delete_compute(id_temp);
461 tflag = 0;
462 }
463 delete [] id_temp;
464 id_temp = utils::strdup(arg[1]);
465
466 int icompute = modify->find_compute(arg[1]);
467 if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID");
468 temperature = modify->compute[icompute];
469
470 if (temperature->tempflag == 0)
471 error->all(FLERR,"Fix_modify temperature ID does not compute temperature");
472 if (temperature->igroup != 0 && comm->me == 0)
473 error->warning(FLERR,"Temperature for NPT is not for group all");
474
475 // reset id_temp of pressure to new temperature ID
476
477 icompute = modify->find_compute(id_press);
478 if (icompute < 0)
479 error->all(FLERR,"Pressure ID for fix press/berendsen does not exist");
480 modify->compute[icompute]->reset_extra_compute_fix(id_temp);
481
482 return 2;
483
484 } else if (strcmp(arg[0],"press") == 0) {
485 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
486 if (pflag) {
487 modify->delete_compute(id_press);
488 pflag = 0;
489 }
490 delete [] id_press;
491 id_press = utils::strdup(arg[1]);
492
493 int icompute = modify->find_compute(arg[1]);
494 if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
495 pressure = modify->compute[icompute];
496
497 if (pressure->pressflag == 0)
498 error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
499 return 2;
500 }
501 return 0;
502 }
503