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 author: Aidan Thompson (SNL)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_box_relax.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "compute.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "force.h"
27 #include "kspace.h"
28 #include "math_extra.h"
29 #include "modify.h"
30 #include "update.h"
31
32 #include <cmath>
33 #include <cstring>
34
35 using namespace LAMMPS_NS;
36 using namespace FixConst;
37
38 enum{NONE,XYZ,XY,YZ,XZ};
39 enum{ISO,ANISO,TRICLINIC};
40
41 #define MAX_LIFO_DEPTH 2 // 3 box0 arrays in *.h dimensioned to this
42
43 /* ---------------------------------------------------------------------- */
44
FixBoxRelax(LAMMPS * lmp,int narg,char ** arg)45 FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
46 Fix(lmp, narg, arg),
47 id_temp(nullptr), id_press(nullptr), tflag(0), pflag(0)
48 {
49 if (narg < 5) error->all(FLERR,"Illegal fix box/relax command");
50
51 scalar_flag = 1;
52 extscalar = 1;
53 global_freq = 1;
54 no_change_box = 1;
55
56 // default values
57
58 pcouple = NONE;
59 allremap = 1;
60 vmax = 0.0001;
61 deviatoric_flag = 0;
62 nreset_h0 = 0;
63
64 p_target[0] = p_target[1] = p_target[2] =
65 p_target[3] = p_target[4] = p_target[5] = 0.0;
66 p_flag[0] = p_flag[1] = p_flag[2] =
67 p_flag[3] = p_flag[4] = p_flag[5] = 0;
68
69 // turn on tilt factor scaling, whenever applicable
70
71 dimension = domain->dimension;
72
73 scaleyz = scalexz = scalexy = 0;
74 if (domain->yperiodic && domain->xy != 0.0) scalexy = 1;
75 if (domain->zperiodic && dimension == 3) {
76 if (domain->yz != 0.0) scaleyz = 1;
77 if (domain->xz != 0.0) scalexz = 1;
78 }
79
80 // set fixed-point to default = center of cell
81
82 fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]);
83 fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]);
84 fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]);
85
86 // process keywords
87
88 int iarg = 3;
89
90 while (iarg < narg) {
91 if (strcmp(arg[iarg],"iso") == 0) {
92 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
93 pcouple = XYZ;
94 p_target[0] = p_target[1] = p_target[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
95 p_flag[0] = p_flag[1] = p_flag[2] = 1;
96 if (dimension == 2) {
97 p_target[2] = 0.0;
98 p_flag[2] = 0;
99 }
100 iarg += 2;
101 } else if (strcmp(arg[iarg],"aniso") == 0) {
102 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
103 pcouple = NONE;
104 p_target[0] = p_target[1] = p_target[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
105 p_flag[0] = p_flag[1] = p_flag[2] = 1;
106 if (dimension == 2) {
107 p_target[2] = 0.0;
108 p_flag[2] = 0;
109 }
110 iarg += 2;
111 } else if (strcmp(arg[iarg],"tri") == 0) {
112 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
113 pcouple = NONE;
114 scalexy = scalexz = scaleyz = 0;
115 p_target[0] = p_target[1] = p_target[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
116 p_flag[0] = p_flag[1] = p_flag[2] = 1;
117 p_target[3] = p_target[4] = p_target[5] = 0.0;
118 p_flag[3] = p_flag[4] = p_flag[5] = 1;
119 if (dimension == 2) {
120 p_target[2] = p_target[3] = p_target[4] = 0.0;
121 p_flag[2] = p_flag[3] = p_flag[4] = 0;
122 }
123 iarg += 2;
124
125 } else if (strcmp(arg[iarg],"x") == 0) {
126 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
127 p_target[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
128 p_flag[0] = 1;
129 deviatoric_flag = 1;
130 iarg += 2;
131 } else if (strcmp(arg[iarg],"y") == 0) {
132 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
133 p_target[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
134 p_flag[1] = 1;
135 deviatoric_flag = 1;
136 iarg += 2;
137 } else if (strcmp(arg[iarg],"z") == 0) {
138 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
139 p_target[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
140 p_flag[2] = 1;
141 deviatoric_flag = 1;
142 iarg += 2;
143 if (dimension == 2)
144 error->all(FLERR,"Invalid fix box/relax command for a 2d simulation");
145
146 } else if (strcmp(arg[iarg],"yz") == 0) {
147 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
148 p_target[3] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
149 p_flag[3] = 1;
150 deviatoric_flag = 1;
151 scaleyz = 0;
152 iarg += 2;
153 if (dimension == 2)
154 error->all(FLERR,"Invalid fix box/relax command for a 2d simulation");
155 } else if (strcmp(arg[iarg],"xz") == 0) {
156 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
157 p_target[4] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
158 p_flag[4] = 1;
159 deviatoric_flag = 1;
160 scalexz = 0;
161 iarg += 2;
162 if (dimension == 2)
163 error->all(FLERR,"Invalid fix box/relax command for a 2d simulation");
164 } else if (strcmp(arg[iarg],"xy") == 0) {
165 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
166 p_target[5] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
167 p_flag[5] = 1;
168 deviatoric_flag = 1;
169 scalexy = 0;
170 iarg += 2;
171
172 } else if (strcmp(arg[iarg],"couple") == 0) {
173 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
174 if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
175 else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
176 else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
177 else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
178 else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
179 else error->all(FLERR,"Illegal fix box/relax command");
180 iarg += 2;
181
182 } else if (strcmp(arg[iarg],"dilate") == 0) {
183 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
184 if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
185 else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
186 else error->all(FLERR,"Illegal fix box/relax command");
187 iarg += 2;
188 } else if (strcmp(arg[iarg],"vmax") == 0) {
189 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
190 vmax = utils::numeric(FLERR,arg[iarg+1],false,lmp);
191 iarg += 2;
192 } else if (strcmp(arg[iarg],"nreset") == 0) {
193 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
194 nreset_h0 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
195 if (nreset_h0 < 0) error->all(FLERR,"Illegal fix box/relax command");
196 iarg += 2;
197 } else if (strcmp(arg[iarg],"scalexy") == 0) {
198 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
199 if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
200 else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
201 else error->all(FLERR,"Illegal fix box/relax command");
202 iarg += 2;
203 } else if (strcmp(arg[iarg],"scalexz") == 0) {
204 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
205 if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
206 else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
207 else error->all(FLERR,"Illegal fix box/relax command");
208 iarg += 2;
209 } else if (strcmp(arg[iarg],"scaleyz") == 0) {
210 if (iarg+2 > narg) error->all(FLERR,"Illegal fix box/relax command");
211 if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
212 else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
213 else error->all(FLERR,"Illegal fix box/relax command");
214 iarg += 2;
215 } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
216 if (iarg+4 > narg) error->all(FLERR,"Illegal fix box/relax command");
217 fixedpoint[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
218 fixedpoint[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
219 fixedpoint[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
220 iarg += 4;
221 } else error->all(FLERR,"Illegal fix box/relax command");
222 }
223
224 if (p_flag[0]) box_change |= BOX_CHANGE_X;
225 if (p_flag[1]) box_change |= BOX_CHANGE_Y;
226 if (p_flag[2]) box_change |= BOX_CHANGE_Z;
227 if (p_flag[3]) box_change |= BOX_CHANGE_YZ;
228 if (p_flag[4]) box_change |= BOX_CHANGE_XZ;
229 if (p_flag[5]) box_change |= BOX_CHANGE_XY;
230
231 if (allremap == 0) restart_pbc = 1;
232
233 // error checks
234
235 if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4]))
236 error->all(FLERR,"Invalid fix box/relax command for a 2d simulation");
237 if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
238 error->all(FLERR,"Invalid fix box/relax command for a 2d simulation");
239
240 if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
241 error->all(FLERR,"Invalid fix box/relax command pressure settings");
242 if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
243 error->all(FLERR,"Invalid fix box/relax command pressure settings");
244 if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
245 error->all(FLERR,"Invalid fix box/relax command pressure settings");
246 if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
247 error->all(FLERR,"Invalid fix box/relax command pressure settings");
248 if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
249 error->all(FLERR,"Invalid fix box/relax command pressure settings");
250
251 // require periodicity in tensile dimension
252
253 if (p_flag[0] && domain->xperiodic == 0)
254 error->all(FLERR,"Cannot use fix box/relax on a non-periodic dimension");
255 if (p_flag[1] && domain->yperiodic == 0)
256 error->all(FLERR,"Cannot use fix box/relax on a non-periodic dimension");
257 if (p_flag[2] && domain->zperiodic == 0)
258 error->all(FLERR,"Cannot use fix box/relax on a non-periodic dimension");
259
260 // require periodicity in 2nd dim of off-diagonal tilt component
261
262 if (p_flag[3] && domain->zperiodic == 0)
263 error->all(FLERR,
264 "Cannot use fix box/relax on a 2nd non-periodic dimension");
265 if (p_flag[4] && domain->zperiodic == 0)
266 error->all(FLERR,
267 "Cannot use fix box/relax on a 2nd non-periodic dimension");
268 if (p_flag[5] && domain->yperiodic == 0)
269 error->all(FLERR,
270 "Cannot use fix box/relax on a 2nd non-periodic dimension");
271
272 if (scaleyz == 1 && domain->zperiodic == 0)
273 error->all(FLERR,"Cannot use fix box/relax "
274 "with tilt factor scaling on a 2nd non-periodic dimension");
275 if (scalexz == 1 && domain->zperiodic == 0)
276 error->all(FLERR,"Cannot use fix box/relax "
277 "with tilt factor scaling on a 2nd non-periodic dimension");
278 if (scalexy == 1 && domain->yperiodic == 0)
279 error->all(FLERR,"Cannot use fix box/relax "
280 "with tilt factor scaling on a 2nd non-periodic dimension");
281
282 if (p_flag[3] && scaleyz == 1)
283 error->all(FLERR,"Cannot use fix box/relax with "
284 "both relaxation and scaling on a tilt factor");
285 if (p_flag[4] && scalexz == 1)
286 error->all(FLERR,"Cannot use fix box/relax with "
287 "both relaxation and scaling on a tilt factor");
288 if (p_flag[5] && scalexy == 1)
289 error->all(FLERR,"Cannot use fix box/relax with "
290 "both relaxation and scaling on a tilt factor");
291
292 if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5]))
293 error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in "
294 "fix box/relax with non-triclinic box");
295
296 if (pcouple == XYZ && dimension == 3 &&
297 (p_target[0] != p_target[1] || p_target[0] != p_target[2]))
298 error->all(FLERR,"Invalid fix box/relax pressure settings");
299 if (pcouple == XYZ && dimension == 2 && p_target[0] != p_target[1])
300 error->all(FLERR,"Invalid fix box/relax pressure settings");
301 if (pcouple == XY && p_target[0] != p_target[1])
302 error->all(FLERR,"Invalid fix box/relax pressure settings");
303 if (pcouple == YZ && p_target[1] != p_target[2])
304 error->all(FLERR,"Invalid fix box/relax pressure settings");
305 if (pcouple == XZ && p_target[0] != p_target[2])
306 error->all(FLERR,"Invalid fix box/relax pressure settings");
307
308 if (vmax <= 0.0) error->all(FLERR,"Illegal fix box/relax command");
309
310 // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof
311 // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof
312 // else pstyle = ANISO -> 3 dof
313
314 if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC;
315 else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
316 else pstyle = ANISO;
317
318 // create a new compute temp style
319 // id = fix-ID + temp
320 // compute group = all since pressure is always global (group all)
321 // and thus its KE/temperature contribution should use group all
322
323 id_temp = utils::strdup(std::string(id) + "_temp");
324 modify->add_compute(fmt::format("{} all temp",id_temp));
325 tflag = 1;
326
327 // create a new compute pressure style (virial only)
328 // id = fix-ID + press, compute group = all
329 // pass id_temp as 4th arg to pressure constructor
330
331 id_press = utils::strdup(std::string(id) + "_press");
332 modify->add_compute(fmt::format("{} all pressure {} virial",id_press, id_temp));
333 pflag = 1;
334
335 dimension = domain->dimension;
336 nrigid = 0;
337 rfix = nullptr;
338
339 current_lifo = 0;
340 }
341
342 /* ---------------------------------------------------------------------- */
343
~FixBoxRelax()344 FixBoxRelax::~FixBoxRelax()
345 {
346 delete [] rfix;
347
348 // delete temperature and pressure if fix created them
349
350 if (tflag) modify->delete_compute(id_temp);
351 if (pflag) modify->delete_compute(id_press);
352 delete [] id_temp;
353 delete [] id_press;
354 }
355
356 /* ---------------------------------------------------------------------- */
357
setmask()358 int FixBoxRelax::setmask()
359 {
360 int mask = 0;
361 mask |= MIN_ENERGY;
362 return mask;
363 }
364
365 /* ---------------------------------------------------------------------- */
366
init()367 void FixBoxRelax::init()
368 {
369 // set temperature and pressure ptrs
370
371 int icompute = modify->find_compute(id_temp);
372 if (icompute < 0)
373 error->all(FLERR,"Temperature ID for fix box/relax does not exist");
374 temperature = modify->compute[icompute];
375
376 icompute = modify->find_compute(id_press);
377 if (icompute < 0)
378 error->all(FLERR,"Pressure ID for fix box/relax does not exist");
379 pressure = modify->compute[icompute];
380
381 pv2e = 1.0 / force->nktv2p;
382
383 if (force->kspace) kspace_flag = 1;
384 else kspace_flag = 0;
385
386 // detect if any rigid fixes exist so rigid bodies move when box is remapped
387 // rfix[] = indices to each fix rigid
388
389 delete [] rfix;
390 nrigid = 0;
391 rfix = nullptr;
392
393 for (int i = 0; i < modify->nfix; i++)
394 if (modify->fix[i]->rigid_flag) nrigid++;
395 if (nrigid) {
396 rfix = new int[nrigid];
397 nrigid = 0;
398 for (int i = 0; i < modify->nfix; i++)
399 if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
400 }
401
402 // initial box dimensions
403
404 xprdinit = domain->xprd;
405 yprdinit = domain->yprd;
406 zprdinit = domain->zprd;
407 if (dimension == 2) zprdinit = 1.0;
408 vol0 = xprdinit * yprdinit * zprdinit;
409 h0[0] = domain->h[0];
410 h0[1] = domain->h[1];
411 h0[2] = domain->h[2];
412 h0[3] = domain->h[3];
413 h0[4] = domain->h[4];
414 h0[5] = domain->h[5];
415
416 // hydrostatic target pressure and deviatoric target stress
417
418 compute_press_target();
419 if (deviatoric_flag) compute_sigma();
420 }
421
422 /* ----------------------------------------------------------------------
423 compute energy and force due to extra degrees of freedom
424 ------------------------------------------------------------------------- */
425
min_energy(double * fextra)426 double FixBoxRelax::min_energy(double *fextra)
427 {
428 double eng,scale,scalex,scaley,scalez,scalevol;
429
430 temperature->compute_scalar();
431 if (pstyle == ISO) pressure->compute_scalar();
432 else {
433 temperature->compute_vector();
434 pressure->compute_vector();
435 }
436 couple();
437
438 // trigger virial computation on every iteration of minimizer
439
440 pressure->addstep(update->ntimestep+1);
441
442 // compute energy, forces for each extra degree of freedom
443 // returned eng = PV must be in units of energy
444 // returned fextra must likewise be in units of energy
445
446 if (pstyle == ISO) {
447 scale = domain->xprd/xprdinit;
448 if (dimension == 3) {
449 eng = pv2e * p_target[0] * (scale*scale*scale-1.0)*vol0;
450 fextra[0] = pv2e * (p_current[0] - p_target[0])*3.0*scale*scale*vol0;
451 } else {
452 eng = pv2e * p_target[0] * (scale*scale-1.0)*vol0;
453 fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*vol0;
454 }
455
456 } else {
457 fextra[0] = fextra[1] = fextra[2] = 0.0;
458 scalex = scaley = scalez = 1.0;
459 if (p_flag[0]) scalex = domain->xprd/xprdinit;
460 if (p_flag[1]) scaley = domain->yprd/yprdinit;
461 if (p_flag[2]) scalez = domain->zprd/zprdinit;
462 scalevol = scalex*scaley*scalez;
463 eng = pv2e * p_hydro * (scalevol-1.0)*vol0;
464 if (p_flag[0])
465 fextra[0] = pv2e * (p_current[0] - p_hydro)*scaley*scalez*vol0;
466 if (p_flag[1])
467 fextra[1] = pv2e * (p_current[1] - p_hydro)*scalex*scalez*vol0;
468 if (p_flag[2])
469 fextra[2] = pv2e * (p_current[2] - p_hydro)*scalex*scaley*vol0;
470
471 if (pstyle == TRICLINIC) {
472 fextra[3] = fextra[4] = fextra[5] = 0.0;
473 if (p_flag[3])
474 fextra[3] = pv2e*p_current[3]*scaley*yprdinit*scalex*xprdinit*yprdinit;
475 if (p_flag[4])
476 fextra[4] = pv2e*p_current[4]*scalex*xprdinit*scaley*yprdinit*xprdinit;
477 if (p_flag[5])
478 fextra[5] = pv2e*p_current[5]*scalex*xprdinit*scalez*zprdinit*xprdinit;
479 }
480
481 if (deviatoric_flag) {
482 compute_deviatoric();
483 if (p_flag[0]) fextra[0] -= fdev[0]*xprdinit;
484 if (p_flag[1]) fextra[1] -= fdev[1]*yprdinit;
485 if (p_flag[2]) fextra[2] -= fdev[2]*zprdinit;
486 if (pstyle == TRICLINIC) {
487 if (p_flag[3]) fextra[3] -= fdev[3]*yprdinit;
488 if (p_flag[4]) fextra[4] -= fdev[4]*xprdinit;
489 if (p_flag[5]) fextra[5] -= fdev[5]*xprdinit;
490 }
491
492 eng += compute_strain_energy();
493 }
494 }
495
496 return eng;
497 }
498
499 /* ----------------------------------------------------------------------
500 store extra dof values for minimization linesearch starting point
501 boxlo0,boxhi0 = box dimensions
502 box values are pushed onto a LIFO stack so nested calls can be made
503 values are popped by calling min_step(0.0)
504 ------------------------------------------------------------------------- */
505
min_store()506 void FixBoxRelax::min_store()
507 {
508 for (int i = 0; i < 3; i++) {
509 boxlo0[current_lifo][i] = domain->boxlo[i];
510 boxhi0[current_lifo][i] = domain->boxhi[i];
511 }
512 if (pstyle == TRICLINIC) {
513 boxtilt0[current_lifo][0] = domain->yz;
514 boxtilt0[current_lifo][1] = domain->xz;
515 boxtilt0[current_lifo][2] = domain->xy;
516 }
517 }
518
519 /* ----------------------------------------------------------------------
520 clear the LIFO stack for min_store
521 ------------------------------------------------------------------------- */
522
min_clearstore()523 void FixBoxRelax::min_clearstore()
524 {
525 current_lifo = 0;
526 }
527
528 /* ----------------------------------------------------------------------
529 push the LIFO stack for min_store
530 ------------------------------------------------------------------------- */
531
min_pushstore()532 void FixBoxRelax::min_pushstore()
533 {
534 if (current_lifo >= MAX_LIFO_DEPTH) {
535 error->all(FLERR,"Attempt to push beyond stack limit in fix box/relax");
536 return;
537 }
538 current_lifo++;
539 }
540
541
542 /* ----------------------------------------------------------------------
543 pop the LIFO stack for min_store
544 ------------------------------------------------------------------------- */
545
min_popstore()546 void FixBoxRelax::min_popstore()
547 {
548 if (current_lifo <= 0) {
549 error->all(FLERR,"Attempt to pop empty stack in fix box/relax");
550 return;
551 }
552 current_lifo--;
553 }
554
555 /* ----------------------------------------------------------------------
556 check if time to reset reference state. If so, do so.
557 ------------------------------------------------------------------------- */
558
min_reset_ref()559 int FixBoxRelax::min_reset_ref()
560 {
561 int itmp = 0;
562
563 // if nreset_h0 > 0, reset reference box
564 // every nreset_h0 timesteps
565 // only needed for deviatoric external stress
566
567 if (deviatoric_flag && nreset_h0 > 0) {
568 int delta = update->ntimestep - update->beginstep;
569 if (delta % nreset_h0 == 0) {
570 compute_sigma();
571 itmp = 1;
572 }
573 }
574 return itmp;
575 }
576
577 /* ----------------------------------------------------------------------
578 change the box dimensions by fraction ds = alpha*hextra
579 ------------------------------------------------------------------------- */
580
min_step(double alpha,double * hextra)581 void FixBoxRelax::min_step(double alpha, double *hextra)
582 {
583 if (pstyle == ISO) {
584 ds[0] = ds[1] = ds[2] = alpha*hextra[0];
585 } else {
586 ds[0] = ds[1] = ds[2] = 0.0;
587 if (p_flag[0]) ds[0] = alpha*hextra[0];
588 if (p_flag[1]) ds[1] = alpha*hextra[1];
589 if (p_flag[2]) ds[2] = alpha*hextra[2];
590 if (pstyle == TRICLINIC) {
591 ds[3] = ds[4] = ds[5] = 0.0;
592 if (p_flag[3]) ds[3] = alpha*hextra[3];
593 if (p_flag[4]) ds[4] = alpha*hextra[4];
594 if (p_flag[5]) ds[5] = alpha*hextra[5];
595 }
596 }
597 remap();
598 if (kspace_flag) force->kspace->setup();
599 }
600
601 /* ----------------------------------------------------------------------
602 max allowed step size along hextra
603 ------------------------------------------------------------------------- */
604
max_alpha(double * hextra)605 double FixBoxRelax::max_alpha(double *hextra)
606 {
607 double alpha = 1.0;
608 if (pstyle == ISO) alpha = vmax/fabs(hextra[0]);
609 else {
610 if (p_flag[0]) alpha = MIN(alpha,vmax/fabs(hextra[0]));
611 if (p_flag[1]) alpha = MIN(alpha,vmax/fabs(hextra[1]));
612 if (p_flag[2]) alpha = MIN(alpha,vmax/fabs(hextra[2]));
613 if (pstyle == TRICLINIC) {
614 if (p_flag[3]) alpha = MIN(alpha,vmax/fabs(hextra[3]));
615 if (p_flag[4]) alpha = MIN(alpha,vmax/fabs(hextra[4]));
616 if (p_flag[5]) alpha = MIN(alpha,vmax/fabs(hextra[5]));
617 }
618 }
619 return alpha;
620 }
621
622 /* ----------------------------------------------------------------------
623 return number of degrees of freedom added by this fix
624 ------------------------------------------------------------------------- */
625
min_dof()626 int FixBoxRelax::min_dof()
627 {
628 if (pstyle == ISO) return 1;
629 if (pstyle == TRICLINIC) return 6;
630 return 3;
631 }
632
633 /* ----------------------------------------------------------------------
634 dilate the box and owned/ghost atoms around center of box
635 ------------------------------------------------------------------------- */
636
remap()637 void FixBoxRelax::remap()
638 {
639 int i,n;
640
641 // rescale simulation box from linesearch starting point
642 // scale atom coords for all atoms or only for fix group atoms
643
644 double **x = atom->x;
645 int *mask = atom->mask;
646 n = atom->nlocal + atom->nghost;
647
648 // convert pertinent atoms and rigid bodies to lamda coords
649
650 if (allremap) domain->x2lamda(n);
651 else {
652 for (i = 0; i < n; i++)
653 if (mask[i] & groupbit)
654 domain->x2lamda(x[i],x[i]);
655 }
656
657 if (nrigid)
658 for (i = 0; i < nrigid; i++)
659 modify->fix[rfix[i]]->deform(0);
660
661 // reset global and local box to new size/shape
662
663 for (i = 0; i < 3; i++)
664 if (p_flag[i]) {
665 double currentBoxLo0 = boxlo0[current_lifo][i];
666 double currentBoxHi0 = boxhi0[current_lifo][i];
667 domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0 - fixedpoint[i])/domain->h[i]*ds[i]*h0[i];
668 domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0 - fixedpoint[i])/domain->h[i]*ds[i]*h0[i];
669 if (domain->boxlo[i] >= domain->boxhi[i])
670 error->all(FLERR,"Fix box/relax generated negative box length");
671 }
672
673 // scale tilt factors with cell, if set
674
675 if (scaleyz) domain->yz = (domain->boxhi[2] - domain->boxlo[2])*h0[3]/h0[2];
676 if (scalexz) domain->xz = (domain->boxhi[2] - domain->boxlo[2])*h0[4]/h0[2];
677 if (scalexy) domain->xy = (domain->boxhi[1] - domain->boxlo[1])*h0[5]/h0[1];
678
679 if (pstyle == TRICLINIC) {
680 if (p_flag[3]) domain->yz = boxtilt0[current_lifo][0]+ds[3]*yprdinit;
681 if (p_flag[4]) domain->xz = boxtilt0[current_lifo][1]+ds[4]*xprdinit;
682 if (p_flag[5]) domain->xy = boxtilt0[current_lifo][2]+ds[5]*xprdinit;
683 }
684
685 domain->set_global_box();
686 domain->set_local_box();
687
688 // convert pertinent atoms and rigid bodies back to box coords
689
690 if (allremap) domain->lamda2x(n);
691 else {
692 for (i = 0; i < n; i++)
693 if (mask[i] & groupbit)
694 domain->lamda2x(x[i],x[i]);
695 }
696
697 if (nrigid)
698 for (i = 0; i < nrigid; i++)
699 modify->fix[rfix[i]]->deform(1);
700 }
701
702 /* ---------------------------------------------------------------------- */
703
couple()704 void FixBoxRelax::couple()
705 {
706 double *tensor = pressure->vector;
707
708 if (pstyle == ISO)
709 p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
710 else if (pcouple == XYZ) {
711 double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
712 p_current[0] = p_current[1] = p_current[2] = ave;
713 } else if (pcouple == XY) {
714 double ave = 0.5 * (tensor[0] + tensor[1]);
715 p_current[0] = p_current[1] = ave;
716 p_current[2] = tensor[2];
717 } else if (pcouple == YZ) {
718 double ave = 0.5 * (tensor[1] + tensor[2]);
719 p_current[1] = p_current[2] = ave;
720 p_current[0] = tensor[0];
721 } else if (pcouple == XZ) {
722 double ave = 0.5 * (tensor[0] + tensor[2]);
723 p_current[0] = p_current[2] = ave;
724 p_current[1] = tensor[1];
725 } else {
726 p_current[0] = tensor[0];
727 p_current[1] = tensor[1];
728 p_current[2] = tensor[2];
729 }
730
731 if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2]))
732 error->all(FLERR,"Non-numeric pressure - simulation unstable");
733
734 // switch order from xy-xz-yz to Voigt ordering
735
736 if (pstyle == TRICLINIC) {
737 p_current[3] = tensor[5];
738 p_current[4] = tensor[4];
739 p_current[5] = tensor[3];
740
741 if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5]))
742 error->all(FLERR,"Non-numeric pressure - simulation unstable");
743 }
744 }
745
746 /* ---------------------------------------------------------------------- */
747
modify_param(int narg,char ** arg)748 int FixBoxRelax::modify_param(int narg, char **arg)
749 {
750 if (strcmp(arg[0],"temp") == 0) {
751 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
752 if (tflag) {
753 modify->delete_compute(id_temp);
754 tflag = 0;
755 }
756 delete [] id_temp;
757 id_temp = utils::strdup(arg[1]);
758
759 int icompute = modify->find_compute(arg[1]);
760 if (icompute < 0)
761 error->all(FLERR,"Could not find fix_modify temperature ID");
762 temperature = modify->compute[icompute];
763
764 if (temperature->tempflag == 0)
765 error->all(FLERR,
766 "Fix_modify temperature ID does not compute temperature");
767 if (temperature->igroup != 0 && comm->me == 0)
768 error->warning(FLERR,"Temperature for fix modify is not for group all");
769
770 // reset id_temp of pressure to new temperature ID
771
772 icompute = modify->find_compute(id_press);
773 if (icompute < 0)
774 error->all(FLERR,"Pressure ID for fix modify does not exist");
775 modify->compute[icompute]->reset_extra_compute_fix(id_temp);
776
777 return 2;
778
779 } else if (strcmp(arg[0],"press") == 0) {
780 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
781 if (pflag) {
782 modify->delete_compute(id_press);
783 pflag = 0;
784 }
785 delete [] id_press;
786 id_press = utils::strdup(arg[1]);
787
788 int icompute = modify->find_compute(arg[1]);
789 if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
790 pressure = modify->compute[icompute];
791
792 if (pressure->pressflag == 0)
793 error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
794 return 2;
795 }
796 return 0;
797 }
798
799 /* ----------------------------------------------------------------------
800 compute sigma tensor (needed whenever reference box is reset)
801 -----------------------------------------------------------------------*/
802
compute_sigma()803 void FixBoxRelax::compute_sigma()
804 {
805 double pdeviatoric[3][3];
806 double tmp1[3][3],sigma_tensor[3][3],h_invtmp[3][3];
807
808 // reset reference box dimensions
809
810 xprdinit = domain->xprd;
811 yprdinit = domain->yprd;
812 zprdinit = domain->zprd;
813 if (dimension == 2) zprdinit = 1.0;
814 vol0 = xprdinit * yprdinit * zprdinit;
815
816 h0_inv[0] = domain->h_inv[0];
817 h0_inv[1] = domain->h_inv[1];
818 h0_inv[2] = domain->h_inv[2];
819 h0_inv[3] = domain->h_inv[3];
820 h0_inv[4] = domain->h_inv[4];
821 h0_inv[5] = domain->h_inv[5];
822
823 h_invtmp[0][0] = h0_inv[0];
824 h_invtmp[1][1] = h0_inv[1];
825 h_invtmp[2][2] = h0_inv[2];
826 h_invtmp[1][2] = h0_inv[3];
827 h_invtmp[0][2] = h0_inv[4];
828 h_invtmp[0][1] = h0_inv[5];
829 h_invtmp[2][1] = 0.0;
830 h_invtmp[2][0] = 0.0;
831 h_invtmp[1][0] = 0.0;
832
833 // compute target deviatoric stress tensor pdevmod
834
835 pdeviatoric[0][0] = pdeviatoric[1][1] = pdeviatoric[2][2] = 0.0;
836 if (p_flag[0]) pdeviatoric[0][0] = p_target[0] - p_hydro;
837 if (p_flag[1]) pdeviatoric[1][1] = p_target[1] - p_hydro;
838 if (p_flag[2]) pdeviatoric[2][2] = p_target[2] - p_hydro;
839 pdeviatoric[1][2] = pdeviatoric[2][1] = p_target[3];
840 pdeviatoric[0][2] = pdeviatoric[2][0] = p_target[4];
841 pdeviatoric[0][1] = pdeviatoric[1][0] = p_target[5];
842
843 // Modify to account for off-diagonal terms
844 // These equations come from the stationarity relation:
845 // Pdev,sys = Pdev,targ*hinv^t*hdiag
846 // where:
847 // Pdev,sys is the system deviatoric stress tensor,
848 // Pdev,targ = pdeviatoric, effective target deviatoric stress
849 // hinv^t is the transpose of the inverse h tensor
850 // hdiag is the diagonal part of the h tensor
851
852 pdeviatoric[1][1] -= pdeviatoric[1][2]*h0_inv[3]*h0[1];
853 pdeviatoric[0][1] -= pdeviatoric[0][2]*h0_inv[3]*h0[1];
854 pdeviatoric[1][0] = pdeviatoric[0][1];
855 pdeviatoric[0][0] -= pdeviatoric[0][1]*h0_inv[5]*h0[0] +
856 pdeviatoric[0][2]*h0_inv[4]*h0[0];
857
858 // compute symmetric sigma tensor
859
860 MathExtra::times3(h_invtmp,pdeviatoric,tmp1);
861 MathExtra::times3_transpose(tmp1,h_invtmp,sigma_tensor);
862 MathExtra::scalar_times3(vol0,sigma_tensor);
863
864 sigma[0] = sigma_tensor[0][0];
865 sigma[1] = sigma_tensor[1][1];
866 sigma[2] = sigma_tensor[2][2];
867 sigma[3] = sigma_tensor[1][2];
868 sigma[4] = sigma_tensor[0][2];
869 sigma[5] = sigma_tensor[0][1];
870 }
871
872 /* ----------------------------------------------------------------------
873 compute strain energy
874 -----------------------------------------------------------------------*/
875
compute_strain_energy()876 double FixBoxRelax::compute_strain_energy()
877 {
878 // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units
879
880 double* h = domain->h;
881 double d0,d1,d2;
882
883 if (dimension == 3) {
884 d0 =
885 sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) +
886 sigma[5]*( h[1]*h[5]+h[3]*h[4]) +
887 sigma[4]*( h[2]*h[4]);
888 d1 =
889 sigma[5]*( h[5]*h[1]+h[4]*h[3]) +
890 sigma[1]*( h[1]*h[1]+h[3]*h[3]) +
891 sigma[3]*( h[2]*h[3]);
892 d2 =
893 sigma[4]*( h[4]*h[2]) +
894 sigma[3]*( h[3]*h[2]) +
895 sigma[2]*( h[2]*h[2]);
896 } else {
897 d0 = sigma[0]*(h[0]*h[0]+h[5]*h[5]) + sigma[5]*h[1]*h[5];
898 d1 = sigma[5]*h[5]*h[1] + sigma[1]*h[1]*h[1];
899 d2 = 0.0;
900 }
901
902 double energy = 0.5*(d0+d1+d2)*pv2e;
903 return energy;
904 }
905
906 /* ----------------------------------------------------------------------
907 compute deviatoric barostat force = h*sigma*h^t
908 -----------------------------------------------------------------------*/
909
compute_deviatoric()910 void FixBoxRelax::compute_deviatoric()
911 {
912 double* h = domain->h;
913
914 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ]
915 // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ]
916 // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ]
917
918 if (dimension == 3) {
919 fdev[0] = pv2e*(h[0]*sigma[0]+h[5]*sigma[5]+h[4]*sigma[4]);
920 fdev[1] = pv2e*(h[1]*sigma[1]+h[3]*sigma[3]);
921 fdev[2] = pv2e*(h[2]*sigma[2]);
922 fdev[3] = pv2e*(h[1]*sigma[3]+h[3]*sigma[2]);
923 fdev[4] = pv2e*(h[0]*sigma[4]+h[5]*sigma[3]+h[4]*sigma[2]);
924 fdev[5] = pv2e*(h[0]*sigma[5]+h[5]*sigma[1]+h[4]*sigma[3]);
925 } else {
926 fdev[0] = pv2e*(h[0]*sigma[0]+h[5]*sigma[5]);
927 fdev[1] = pv2e*(h[1]*sigma[1]);
928 fdev[5] = pv2e*(h[0]*sigma[5]+h[5]*sigma[1]);
929 }
930 }
931
932 /* ----------------------------------------------------------------------
933 compute hydrostatic target pressure
934 -----------------------------------------------------------------------*/
935
compute_press_target()936 void FixBoxRelax::compute_press_target()
937 {
938 pflagsum = p_flag[0] + p_flag[1] + p_flag[2];
939
940 p_hydro = 0.0;
941 for (int i = 0; i < 3; i++)
942 if (p_flag[i]) p_hydro += p_target[i];
943 if (pflagsum) p_hydro /= pflagsum;
944
945 for (int i = 0; i < 3; i++) {
946 if (p_flag[i] && fabs(p_hydro - p_target[i]) > 1.0e-6) deviatoric_flag = 1;
947 }
948
949 if (pstyle == TRICLINIC) {
950 for (int i = 3; i < 6; i++)
951 if (p_flag[i] && fabs(p_target[i]) > 1.0e-6) deviatoric_flag = 1;
952 }
953 }
954
955 /* ----------------------------------------------------------------------
956 compute PV and strain energy for access to the user
957 ---------------------------------------------------------------------- */
958
compute_scalar()959 double FixBoxRelax::compute_scalar()
960 {
961 double ftmp[6] = {0.0,0.0,0.0,0.0,0.0,0.0};
962 if (update->ntimestep == 0) return 0.0;
963 return min_energy(ftmp);
964 }
965