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:
17 ------------------------------------------------------------------------- */
18
19 #include "fix_npt_cauchy.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "compute.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "fix_deform.h"
27 #include "fix_store.h"
28 #include "force.h"
29 #include "group.h"
30 #include "irregular.h"
31 #include "kspace.h"
32 #include "math_extra.h"
33 #include "memory.h"
34 #include "modify.h"
35 #include "neighbor.h"
36 #include "respa.h"
37 #include "update.h"
38
39 #include <cmath>
40 #include <cstring>
41
42 using namespace LAMMPS_NS;
43 using namespace FixConst;
44
45 #define DELTAFLIP 0.1
46 #define TILTMAX 1.5
47
48 enum{NOBIAS,BIAS};
49 enum{NONE,XYZ,XY,YZ,XZ};
50 enum{ISO,ANISO,TRICLINIC};
51
52 /* ----------------------------------------------------------------------
53 NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
54 ---------------------------------------------------------------------- */
55
FixNPTCauchy(LAMMPS * lmp,int narg,char ** arg)56 FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) :
57 Fix(lmp, narg, arg),
58 rfix(nullptr), id_dilate(nullptr), irregular(nullptr),
59 id_temp(nullptr), id_press(nullptr),
60 eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr),
61 eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr),
62 etap_mass(nullptr), id_store(nullptr),init_store(nullptr)
63 {
64 if (narg < 4) error->all(FLERR,"Illegal fix npt/cauchy command");
65
66 dynamic_group_allow = 1;
67 ecouple_flag = 1;
68 time_integrate = 1;
69 scalar_flag = 1;
70 vector_flag = 1;
71 global_freq = 1;
72 extscalar = 1;
73 extvector = 0;
74
75 // for CauchyStat
76
77 alpha=0.001;
78 initRUN = 0;
79 restartPK = 0;
80 restart_global = 1;
81 restart_stored = 0;
82
83 // default values
84
85 pcouple = NONE;
86 drag = 0.0;
87 allremap = 1;
88 mtchain = mpchain = 3;
89 nc_tchain = nc_pchain = 1;
90 mtk_flag = 1;
91 deviatoric_flag = 0;
92 nreset_h0 = 0;
93 eta_mass_flag = 1;
94 omega_mass_flag = 0;
95 etap_mass_flag = 0;
96 flipflag = 1;
97 dipole_flag = 0;
98 dlm_flag = 0;
99
100 tcomputeflag = 0;
101 pcomputeflag = 0;
102
103 // turn on tilt factor scaling, whenever applicable
104
105 dimension = domain->dimension;
106
107 scaleyz = scalexz = scalexy = 0;
108 if (domain->yperiodic && domain->xy != 0.0) scalexy = 1;
109 if (domain->zperiodic && dimension == 3) {
110 if (domain->yz != 0.0) scaleyz = 1;
111 if (domain->xz != 0.0) scalexz = 1;
112 }
113
114 // set fixed-point to default = center of cell
115
116 fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]);
117 fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]);
118 fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]);
119
120 // used by FixNVTSllod to preserve non-default value
121
122 mtchain_default_flag = 1;
123
124 tstat_flag = 0;
125 double t_period = 0.0;
126
127 double p_period[6];
128 for (int i = 0; i < 6; i++) {
129 p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0;
130 p_flag[i] = 0;
131 }
132
133 // process keywords
134
135 int iarg = 3;
136
137 while (iarg < narg) {
138 if (strcmp(arg[iarg],"temp") == 0) {
139 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
140 tstat_flag = 1;
141 t_start = utils::numeric(FLERR,arg[iarg+1],false,lmp);
142 t_target = t_start;
143 t_stop = utils::numeric(FLERR,arg[iarg+2],false,lmp);
144 t_period = utils::numeric(FLERR,arg[iarg+3],false,lmp);
145 if (t_start <= 0.0 || t_stop <= 0.0)
146 error->all(FLERR,
147 "Target temperature for fix npt/cauchy cannot be 0.0");
148 iarg += 4;
149
150 } else if (strcmp(arg[iarg],"iso") == 0) {
151 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
152 pcouple = XYZ;
153 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
154 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
155 p_period[0] = p_period[1] = p_period[2] =
156 utils::numeric(FLERR,arg[iarg+3],false,lmp);
157 p_flag[0] = p_flag[1] = p_flag[2] = 1;
158 if (dimension == 2) {
159 p_start[2] = p_stop[2] = p_period[2] = 0.0;
160 p_flag[2] = 0;
161 }
162 iarg += 4;
163 } else if (strcmp(arg[iarg],"aniso") == 0) {
164 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
165 pcouple = NONE;
166 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
167 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
168 p_period[0] = p_period[1] = p_period[2] =
169 utils::numeric(FLERR,arg[iarg+3],false,lmp);
170 p_flag[0] = p_flag[1] = p_flag[2] = 1;
171 if (dimension == 2) {
172 p_start[2] = p_stop[2] = p_period[2] = 0.0;
173 p_flag[2] = 0;
174 }
175 iarg += 4;
176 } else if (strcmp(arg[iarg],"tri") == 0) {
177 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
178 pcouple = NONE;
179 scalexy = scalexz = scaleyz = 0;
180 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
181 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
182 p_period[0] = p_period[1] = p_period[2] =
183 utils::numeric(FLERR,arg[iarg+3],false,lmp);
184 p_flag[0] = p_flag[1] = p_flag[2] = 1;
185 p_start[3] = p_start[4] = p_start[5] = 0.0;
186 p_stop[3] = p_stop[4] = p_stop[5] = 0.0;
187 p_period[3] = p_period[4] = p_period[5] =
188 utils::numeric(FLERR,arg[iarg+3],false,lmp);
189 p_flag[3] = p_flag[4] = p_flag[5] = 1;
190 if (dimension == 2) {
191 p_start[2] = p_stop[2] = p_period[2] = 0.0;
192 p_flag[2] = 0;
193 p_start[3] = p_stop[3] = p_period[3] = 0.0;
194 p_flag[3] = 0;
195 p_start[4] = p_stop[4] = p_period[4] = 0.0;
196 p_flag[4] = 0;
197 }
198 iarg += 4;
199 } else if (strcmp(arg[iarg],"x") == 0) {
200 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
201 p_start[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
202 p_stop[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
203 p_period[0] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
204 p_flag[0] = 1;
205 deviatoric_flag = 1;
206 iarg += 4;
207 } else if (strcmp(arg[iarg],"y") == 0) {
208 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
209 p_start[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
210 p_stop[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
211 p_period[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
212 p_flag[1] = 1;
213 deviatoric_flag = 1;
214 iarg += 4;
215 } else if (strcmp(arg[iarg],"z") == 0) {
216 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
217 p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
218 p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
219 p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
220 p_flag[2] = 1;
221 deviatoric_flag = 1;
222 iarg += 4;
223 if (dimension == 2)
224 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
225
226 } else if (strcmp(arg[iarg],"yz") == 0) {
227 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
228 p_start[3] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
229 p_stop[3] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
230 p_period[3] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
231 p_flag[3] = 1;
232 deviatoric_flag = 1;
233 scaleyz = 0;
234 iarg += 4;
235 if (dimension == 2)
236 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
237 } else if (strcmp(arg[iarg],"xz") == 0) {
238 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
239 p_start[4] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
240 p_stop[4] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
241 p_period[4] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
242 p_flag[4] = 1;
243 deviatoric_flag = 1;
244 scalexz = 0;
245 iarg += 4;
246 if (dimension == 2)
247 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
248 } else if (strcmp(arg[iarg],"xy") == 0) {
249 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
250 p_start[5] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
251 p_stop[5] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
252 p_period[5] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
253 p_flag[5] = 1;
254 deviatoric_flag = 1;
255 scalexy = 0;
256 iarg += 4;
257
258 } else if (strcmp(arg[iarg],"couple") == 0) {
259 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
260 if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
261 else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
262 else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
263 else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
264 else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
265 else error->all(FLERR,"Illegal fix npt/cauchy command");
266 iarg += 2;
267
268 } else if (strcmp(arg[iarg],"drag") == 0) {
269 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
270 drag = utils::numeric(FLERR,arg[iarg+1],false,lmp);
271 if (drag < 0.0) error->all(FLERR,"Illegal fix npt/cauchy command");
272 iarg += 2;
273 } else if (strcmp(arg[iarg],"dilate") == 0) {
274 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
275 if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
276 else {
277 allremap = 0;
278 delete [] id_dilate;
279 id_dilate = utils::strdup(arg[iarg+1]);
280 int idilate = group->find(id_dilate);
281 if (idilate == -1)
282 error->all(FLERR,"Fix npt/cauchy dilate group ID does not exist");
283 }
284 iarg += 2;
285
286 } else if (strcmp(arg[iarg],"tchain") == 0) {
287 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
288 mtchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
289 // used by FixNVTSllod to preserve non-default value
290 mtchain_default_flag = 0;
291 if (mtchain < 1) error->all(FLERR,"Illegal fix npt/cauchy command");
292 iarg += 2;
293 } else if (strcmp(arg[iarg],"pchain") == 0) {
294 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
295 mpchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
296 if (mpchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command");
297 iarg += 2;
298 } else if (strcmp(arg[iarg],"mtk") == 0) {
299 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
300 if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
301 else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
302 else error->all(FLERR,"Illegal fix npt/cauchy command");
303 iarg += 2;
304 } else if (strcmp(arg[iarg],"tloop") == 0) {
305 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
306 nc_tchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
307 if (nc_tchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command");
308 iarg += 2;
309 } else if (strcmp(arg[iarg],"ploop") == 0) {
310 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
311 nc_pchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
312 if (nc_pchain < 0) error->all(FLERR,"Illegal fix npt/cauchy command");
313 iarg += 2;
314 } else if (strcmp(arg[iarg],"nreset") == 0) {
315 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
316 nreset_h0 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
317 if (nreset_h0 < 0) error->all(FLERR,"Illegal fix npt/cauchy command");
318 iarg += 2;
319 } else if (strcmp(arg[iarg],"scalexy") == 0) {
320 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
321 if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
322 else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
323 else error->all(FLERR,"Illegal fix npt/cauchy command");
324 iarg += 2;
325 } else if (strcmp(arg[iarg],"scalexz") == 0) {
326 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
327 if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
328 else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
329 else error->all(FLERR,"Illegal fix npt/cauchy command");
330 iarg += 2;
331 } else if (strcmp(arg[iarg],"scaleyz") == 0) {
332 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
333 if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
334 else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
335 else error->all(FLERR,"Illegal fix npt/cauchy command");
336 iarg += 2;
337 } else if (strcmp(arg[iarg],"flip") == 0) {
338 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
339 if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1;
340 else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
341 else error->all(FLERR,"Illegal fix npt/cauchy command");
342 iarg += 2;
343 } else if (strcmp(arg[iarg],"update") == 0) {
344 if (iarg+2 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
345 if (strcmp(arg[iarg+1],"dipole") == 0) dipole_flag = 1;
346 else if (strcmp(arg[iarg+1],"dipole/dlm") == 0) {
347 dipole_flag = 1;
348 dlm_flag = 1;
349 } else error->all(FLERR,"Illegal fix npt/cauchy command");
350 iarg += 2;
351 } else if (strcmp(arg[iarg],"alpha") == 0) {
352 alpha = utils::numeric(FLERR,arg[iarg+1],false,lmp);
353 iarg += 2;
354 } else if (strcmp(arg[iarg],"continue") == 0) {
355 if (strcmp(arg[iarg+1],"yes") != 0 && strcmp(arg[iarg+1],"no") != 0)
356 error->all(FLERR,"Illegal cauchystat continue value. "
357 "Must be 'yes' or 'no'");
358 restartPK = !strcmp(arg[iarg+1],"yes");
359 iarg += 2;
360 } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
361 if (iarg+4 > narg) error->all(FLERR,"Illegal fix npt/cauchy command");
362 fixedpoint[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
363 fixedpoint[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
364 fixedpoint[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
365 iarg += 4;
366
367 // disc keyword is also parsed in fix/nh/sphere
368
369 } else if (strcmp(arg[iarg],"disc") == 0) {
370 iarg++;
371
372 // keywords erate, strain, and ext are also parsed in fix/nh/uef
373
374 } else if (strcmp(arg[iarg],"erate") == 0) {
375 iarg += 3;
376 } else if (strcmp(arg[iarg],"strain") == 0) {
377 iarg += 3;
378 } else if (strcmp(arg[iarg],"ext") == 0) {
379 iarg += 2;
380
381 } else error->all(FLERR,"Illegal fix npt/cauchy command");
382 }
383
384 // error checks
385
386 if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4]))
387 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
388 if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
389 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
390 if (dimension == 2 && (scalexz == 1 || scaleyz == 1 ))
391 error->all(FLERR,"Invalid fix npt/cauchy command for a 2d simulation");
392
393 if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
394 error->all(FLERR,"Invalid fix npt/cauchy command pressure settings");
395 if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
396 error->all(FLERR,"Invalid fix npt/cauchy command pressure settings");
397 if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
398 error->all(FLERR,"Invalid fix npt/cauchy command pressure settings");
399 if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
400 error->all(FLERR,"Invalid fix npt/cauchy command pressure settings");
401 if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
402 error->all(FLERR,"Invalid fix npt/cauchy command pressure settings");
403
404 // require periodicity in tensile dimension
405
406 if (p_flag[0] && domain->xperiodic == 0)
407 error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension");
408 if (p_flag[1] && domain->yperiodic == 0)
409 error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension");
410 if (p_flag[2] && domain->zperiodic == 0)
411 error->all(FLERR,"Cannot use fix npt/cauchy on a non-periodic dimension");
412
413 // require periodicity in 2nd dim of off-diagonal tilt component
414
415 if (p_flag[3] && domain->zperiodic == 0)
416 error->all(FLERR,
417 "Cannot use fix npt/cauchy on a 2nd non-periodic dimension");
418 if (p_flag[4] && domain->zperiodic == 0)
419 error->all(FLERR,
420 "Cannot use fix npt/cauchy on a 2nd non-periodic dimension");
421 if (p_flag[5] && domain->yperiodic == 0)
422 error->all(FLERR,
423 "Cannot use fix npt/cauchy on a 2nd non-periodic dimension");
424
425 if (scaleyz == 1 && domain->zperiodic == 0)
426 error->all(FLERR,"Cannot use fix npt/cauchy "
427 "with yz scaling when z is non-periodic dimension");
428 if (scalexz == 1 && domain->zperiodic == 0)
429 error->all(FLERR,"Cannot use fix npt/cauchy "
430 "with xz scaling when z is non-periodic dimension");
431 if (scalexy == 1 && domain->yperiodic == 0)
432 error->all(FLERR,"Cannot use fix npt/cauchy "
433 "with xy scaling when y is non-periodic dimension");
434
435 if (p_flag[3] && scaleyz == 1)
436 error->all(FLERR,"Cannot use fix npt/cauchy with "
437 "both yz dynamics and yz scaling");
438 if (p_flag[4] && scalexz == 1)
439 error->all(FLERR,"Cannot use fix npt/cauchy with "
440 "both xz dynamics and xz scaling");
441 if (p_flag[5] && scalexy == 1)
442 error->all(FLERR,"Cannot use fix npt/cauchy with "
443 "both xy dynamics and xy scaling");
444
445 if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5]))
446 error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in "
447 "fix npt/cauchy with non-triclinic box");
448
449 if (pcouple == XYZ && dimension == 3 &&
450 (p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
451 p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] ||
452 p_period[0] != p_period[1] || p_period[0] != p_period[2]))
453 error->all(FLERR,"Invalid fix npt/cauchy pressure settings");
454 if (pcouple == XYZ && dimension == 2 &&
455 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
456 p_period[0] != p_period[1]))
457 error->all(FLERR,"Invalid fix npt/cauchy pressure settings");
458 if (pcouple == XY &&
459 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
460 p_period[0] != p_period[1]))
461 error->all(FLERR,"Invalid fix npt/cauchy pressure settings");
462 if (pcouple == YZ &&
463 (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] ||
464 p_period[1] != p_period[2]))
465 error->all(FLERR,"Invalid fix npt/cauchy pressure settings");
466 if (pcouple == XZ &&
467 (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] ||
468 p_period[0] != p_period[2]))
469 error->all(FLERR,"Invalid fix npt/cauchy pressure settings");
470
471 if (dipole_flag) {
472 if (!atom->sphere_flag)
473 error->all(FLERR,"Using update dipole flag requires atom style sphere");
474 if (!atom->mu_flag)
475 error->all(FLERR,"Using update dipole flag requires atom attribute mu");
476 }
477
478 if ((tstat_flag && t_period <= 0.0) ||
479 (p_flag[0] && p_period[0] <= 0.0) ||
480 (p_flag[1] && p_period[1] <= 0.0) ||
481 (p_flag[2] && p_period[2] <= 0.0) ||
482 (p_flag[3] && p_period[3] <= 0.0) ||
483 (p_flag[4] && p_period[4] <= 0.0) ||
484 (p_flag[5] && p_period[5] <= 0.0))
485 error->all(FLERR,"Fix npt/cauchy damping parameters must be > 0.0");
486
487 // set pstat_flag and box change and restart_pbc variables
488
489 pre_exchange_flag = 0;
490 pstat_flag = 0;
491 pstyle = ISO;
492
493 for (int i = 0; i < 6; i++)
494 if (p_flag[i]) pstat_flag = 1;
495
496 if (pstat_flag) {
497 if (p_flag[0]) box_change |= BOX_CHANGE_X;
498 if (p_flag[1]) box_change |= BOX_CHANGE_Y;
499 if (p_flag[2]) box_change |= BOX_CHANGE_Z;
500 if (p_flag[3]) box_change |= BOX_CHANGE_YZ;
501 if (p_flag[4]) box_change |= BOX_CHANGE_XZ;
502 if (p_flag[5]) box_change |= BOX_CHANGE_XY;
503 no_change_box = 1;
504 if (allremap == 0) restart_pbc = 1;
505
506 // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof
507 // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof
508 // else pstyle = ANISO -> 3 dof
509
510 if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC;
511 else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
512 else pstyle = ANISO;
513
514 // pre_exchange only required if flips can occur due to shape changes
515
516 if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
517 pre_exchange_flag = 1;
518 if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 ||
519 domain->xy != 0.0))
520 pre_exchange_flag = 1;
521 }
522
523 // convert input periods to frequencies
524
525 t_freq = 0.0;
526 p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0;
527
528 if (tstat_flag) t_freq = 1.0 / t_period;
529 if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
530 if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
531 if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
532 if (p_flag[3]) p_freq[3] = 1.0 / p_period[3];
533 if (p_flag[4]) p_freq[4] = 1.0 / p_period[4];
534 if (p_flag[5]) p_freq[5] = 1.0 / p_period[5];
535
536 // Nose/Hoover temp and pressure init
537
538 size_vector = 0;
539
540 if (tstat_flag) {
541 int ich;
542 eta = new double[mtchain];
543
544 // add one extra dummy thermostat, set to zero
545
546 eta_dot = new double[mtchain+1];
547 eta_dot[mtchain] = 0.0;
548 eta_dotdot = new double[mtchain];
549 for (ich = 0; ich < mtchain; ich++) {
550 eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
551 }
552 eta_mass = new double[mtchain];
553 size_vector += 2*2*mtchain;
554 }
555
556 if (pstat_flag) {
557 omega[0] = omega[1] = omega[2] = 0.0;
558 omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
559 omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0;
560 omega[3] = omega[4] = omega[5] = 0.0;
561 omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0;
562 omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0;
563 if (pstyle == ISO) size_vector += 2*2*1;
564 else if (pstyle == ANISO) size_vector += 2*2*3;
565 else if (pstyle == TRICLINIC) size_vector += 2*2*6;
566
567 if (mpchain) {
568 int ich;
569 etap = new double[mpchain];
570
571 // add one extra dummy thermostat, set to zero
572
573 etap_dot = new double[mpchain+1];
574 etap_dot[mpchain] = 0.0;
575 etap_dotdot = new double[mpchain];
576 for (ich = 0; ich < mpchain; ich++) {
577 etap[ich] = etap_dot[ich] =
578 etap_dotdot[ich] = 0.0;
579 }
580 etap_mass = new double[mpchain];
581 size_vector += 2*2*mpchain;
582 }
583
584 if (deviatoric_flag) size_vector += 1;
585 }
586
587 nrigid = 0;
588 rfix = nullptr;
589
590 if (pre_exchange_flag) irregular = new Irregular(lmp);
591 else irregular = nullptr;
592
593 // initialize vol0,t0 to zero to signal uninitialized
594 // values then assigned in init(), if necessary
595
596 vol0 = t0 = 0.0;
597
598 if (!tstat_flag)
599 error->all(FLERR,"Temperature control must be used with fix npt/cauchy");
600 if (!pstat_flag)
601 error->all(FLERR,"Pressure control must be used with fix npt/cauchy");
602
603 // create a new compute temp style
604 // id = fix-ID + temp
605 // compute group = all since pressure is always global (group all)
606 // and thus its KE/temperature contribution should use group all
607
608 id_temp = utils::strdup(std::string(id) + "_temp");
609 modify->add_compute(fmt::format("{} all temp",id_temp));
610 tcomputeflag = 1;
611
612 // create a new compute pressure style
613 // id = fix-ID + press, compute group = all
614 // pass id_temp as 4th arg to pressure constructor
615
616 id_press = utils::strdup(std::string(id) + "_press");
617 modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
618 pcomputeflag = 1;
619 }
620
621 /* ---------------------------------------------------------------------- */
622
post_constructor()623 void FixNPTCauchy::post_constructor()
624 {
625 CauchyStat_init();
626 }
627
628 /* ---------------------------------------------------------------------- */
629
~FixNPTCauchy()630 FixNPTCauchy::~FixNPTCauchy()
631 {
632 if (copymode) return;
633
634 delete [] id_dilate;
635 delete [] rfix;
636
637 delete [] id_store;
638 delete irregular;
639
640 // delete temperature and pressure if fix created them
641
642 if (tcomputeflag) modify->delete_compute(id_temp);
643 delete [] id_temp;
644
645 if (tstat_flag) {
646 delete [] eta;
647 delete [] eta_dot;
648 delete [] eta_dotdot;
649 delete [] eta_mass;
650 }
651
652 if (pstat_flag) {
653 if (pcomputeflag) modify->delete_compute(id_press);
654 delete [] id_press;
655 if (mpchain) {
656 delete [] etap;
657 delete [] etap_dot;
658 delete [] etap_dotdot;
659 delete [] etap_mass;
660 }
661 }
662 }
663
664 /* ---------------------------------------------------------------------- */
665
setmask()666 int FixNPTCauchy::setmask()
667 {
668 int mask = 0;
669 mask |= INITIAL_INTEGRATE;
670 mask |= FINAL_INTEGRATE;
671 mask |= INITIAL_INTEGRATE_RESPA;
672 mask |= PRE_FORCE_RESPA;
673 mask |= FINAL_INTEGRATE_RESPA;
674 if (pre_exchange_flag) mask |= PRE_EXCHANGE;
675 return mask;
676 }
677
678 /* ---------------------------------------------------------------------- */
679
init()680 void FixNPTCauchy::init()
681 {
682 // recheck that dilate group has not been deleted
683
684 if (allremap == 0) {
685 int idilate = group->find(id_dilate);
686 if (idilate == -1)
687 error->all(FLERR,"Fix npt/cauchy dilate group ID does not exist");
688 dilate_group_bit = group->bitmask[idilate];
689 }
690
691 // ensure no conflict with fix deform
692
693 if (pstat_flag)
694 for (int i = 0; i < modify->nfix; i++)
695 if (strcmp(modify->fix[i]->style,"deform") == 0) {
696 int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
697 if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
698 (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) ||
699 (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5]))
700 error->all(FLERR,"Cannot use fix npt/cauchy and fix deform "
701 "on same component of stress tensor");
702 }
703
704 // set temperature and pressure ptrs
705
706 int icompute = modify->find_compute(id_temp);
707 if (icompute < 0)
708 error->all(FLERR,"Temperature ID for fix npt/cauchy does not exist");
709 temperature = modify->compute[icompute];
710
711 if (temperature->tempbias) which = BIAS;
712 else which = NOBIAS;
713
714 if (pstat_flag) {
715 icompute = modify->find_compute(id_press);
716 if (icompute < 0)
717 error->all(FLERR,"Pressure ID for fix npt/cauchy does not exist");
718 pressure = modify->compute[icompute];
719 }
720
721 // set timesteps and frequencies
722
723 dtv = update->dt;
724 dtf = 0.5 * update->dt * force->ftm2v;
725 dthalf = 0.5 * update->dt;
726 dt4 = 0.25 * update->dt;
727 dt8 = 0.125 * update->dt;
728 dto = dthalf;
729
730 p_freq_max = 0.0;
731 if (pstat_flag) {
732 p_freq_max = MAX(p_freq[0],p_freq[1]);
733 p_freq_max = MAX(p_freq_max,p_freq[2]);
734 if (pstyle == TRICLINIC) {
735 p_freq_max = MAX(p_freq_max,p_freq[3]);
736 p_freq_max = MAX(p_freq_max,p_freq[4]);
737 p_freq_max = MAX(p_freq_max,p_freq[5]);
738 }
739 pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
740 }
741
742 if (tstat_flag)
743 tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
744
745 // tally the number of dimensions that are barostatted
746 // set initial volume and reference cell, if not already done
747
748 if (pstat_flag) {
749 pdim = p_flag[0] + p_flag[1] + p_flag[2];
750 if (vol0 == 0.0) {
751 if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
752 else vol0 = domain->xprd * domain->yprd;
753 h0_inv[0] = domain->h_inv[0];
754 h0_inv[1] = domain->h_inv[1];
755 h0_inv[2] = domain->h_inv[2];
756 h0_inv[3] = domain->h_inv[3];
757 h0_inv[4] = domain->h_inv[4];
758 h0_inv[5] = domain->h_inv[5];
759 }
760 }
761
762 boltz = force->boltz;
763 nktv2p = force->nktv2p;
764
765 if (force->kspace) kspace_flag = 1;
766 else kspace_flag = 0;
767
768 if (utils::strmatch(update->integrate_style,"^respa")) {
769 nlevels_respa = ((Respa *) update->integrate)->nlevels;
770 step_respa = ((Respa *) update->integrate)->step;
771 dto = 0.5*step_respa[0];
772 }
773
774 // detect if any rigid fixes exist so rigid bodies move when box is remapped
775 // rfix[] = indices to each fix rigid
776
777 delete [] rfix;
778 nrigid = 0;
779 rfix = nullptr;
780
781 for (int i = 0; i < modify->nfix; i++)
782 if (modify->fix[i]->rigid_flag) nrigid++;
783 if (nrigid) {
784 rfix = new int[nrigid];
785 nrigid = 0;
786 for (int i = 0; i < modify->nfix; i++)
787 if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
788 }
789 }
790
791 /* ----------------------------------------------------------------------
792 compute T,P before integrator starts
793 ------------------------------------------------------------------------- */
794
setup(int)795 void FixNPTCauchy::setup(int /*vflag*/)
796 {
797 // tdof needed by compute_temp_target()
798
799 t_current = temperature->compute_scalar();
800 tdof = temperature->dof;
801
802 // t_target is needed by NVT and NPT in compute_scalar()
803 // If no thermostat or using fix nphug,
804 // t_target must be defined by other means.
805
806 if (tstat_flag && strstr(style,"nphug") == nullptr) {
807 compute_temp_target();
808 } else if (pstat_flag) {
809
810 // t0 = reference temperature for masses
811 // cannot be done in init() b/c temperature cannot be called there
812 // is b/c Modify::init() inits computes after fixes due to dof dependence
813 // guesstimate a unit-dependent t0 if actual T = 0.0
814 // if it was read in from a restart file, leave it be
815
816 if (t0 == 0.0) {
817 t0 = temperature->compute_scalar();
818 if (t0 == 0.0) {
819 if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
820 else t0 = 300.0;
821 }
822 }
823 t_target = t0;
824 }
825
826 if (pstat_flag) compute_press_target();
827
828 if (pstat_flag) {
829 if (pstyle == ISO) pressure->compute_scalar();
830 else pressure->compute_vector();
831 couple();
832 pressure->addstep(update->ntimestep+1);
833 }
834
835 // masses and initial forces on thermostat variables
836
837 if (tstat_flag) {
838 eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
839 for (int ich = 1; ich < mtchain; ich++)
840 eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
841 for (int ich = 1; ich < mtchain; ich++) {
842 eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1] -
843 boltz * t_target) / eta_mass[ich];
844 }
845 }
846
847 // masses and initial forces on barostat variables
848
849 if (pstat_flag) {
850 double kt = boltz * t_target;
851 double nkt = (atom->natoms + 1) * kt;
852
853 for (int i = 0; i < 3; i++)
854 if (p_flag[i])
855 omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
856
857 if (pstyle == TRICLINIC) {
858 for (int i = 3; i < 6; i++)
859 if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
860 }
861
862 // masses and initial forces on barostat thermostat variables
863
864 if (mpchain) {
865 etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
866 for (int ich = 1; ich < mpchain; ich++)
867 etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
868 for (int ich = 1; ich < mpchain; ich++)
869 etap_dotdot[ich] =
870 (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
871 boltz * t_target) / etap_mass[ich];
872 }
873 }
874 }
875
876 /* ----------------------------------------------------------------------
877 1st half of Verlet update
878 ------------------------------------------------------------------------- */
879
initial_integrate(int)880 void FixNPTCauchy::initial_integrate(int /*vflag*/)
881 {
882 // update eta_press_dot
883
884 if (pstat_flag && mpchain) nhc_press_integrate();
885
886 // update eta_dot
887
888 if (tstat_flag) {
889 compute_temp_target();
890 nhc_temp_integrate();
891 }
892
893 // need to recompute pressure to account for change in KE
894 // t_current is up-to-date, but compute_temperature is not
895 // compute appropriately coupled elements of mvv_current
896
897 if (pstat_flag) {
898 if (pstyle == ISO) {
899 temperature->compute_scalar();
900 pressure->compute_scalar();
901 } else {
902 temperature->compute_vector();
903 pressure->compute_vector();
904 }
905 couple();
906 pressure->addstep(update->ntimestep+1);
907 }
908
909 if (pstat_flag) {
910 compute_press_target();
911 nh_omega_dot();
912 nh_v_press();
913 }
914
915 nve_v();
916
917 // remap simulation box by 1/2 step
918
919 if (pstat_flag) remap();
920
921 nve_x();
922
923 // remap simulation box by 1/2 step
924 // redo KSpace coeffs since volume has changed
925
926 if (pstat_flag) {
927 remap();
928 if (kspace_flag) force->kspace->setup();
929 }
930 }
931
932 /* ----------------------------------------------------------------------
933 2nd half of Verlet update
934 ------------------------------------------------------------------------- */
935
final_integrate()936 void FixNPTCauchy::final_integrate()
937 {
938 nve_v();
939
940 // re-compute temp before nh_v_press()
941 // only needed for temperature computes with BIAS on reneighboring steps:
942 // b/c some biases store per-atom values (e.g. temp/profile)
943 // per-atom values are invalid if reneigh/comm occurred
944 // since temp->compute() in initial_integrate()
945
946 if (which == BIAS && neighbor->ago == 0)
947 t_current = temperature->compute_scalar();
948
949 if (pstat_flag) nh_v_press();
950
951 // compute new T,P after velocities rescaled by nh_v_press()
952 // compute appropriately coupled elements of mvv_current
953
954 t_current = temperature->compute_scalar();
955 tdof = temperature->dof;
956
957 if (pstat_flag) {
958 if (pstyle == ISO) pressure->compute_scalar();
959 else pressure->compute_vector();
960 couple();
961 pressure->addstep(update->ntimestep+1);
962 }
963
964 if (pstat_flag) nh_omega_dot();
965
966 // update eta_dot
967 // update eta_press_dot
968
969 if (tstat_flag) nhc_temp_integrate();
970 if (pstat_flag && mpchain) nhc_press_integrate();
971 }
972
973 /* ---------------------------------------------------------------------- */
974
initial_integrate_respa(int,int ilevel,int)975 void FixNPTCauchy::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
976 {
977 // set timesteps by level
978
979 dtv = step_respa[ilevel];
980 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
981 dthalf = 0.5 * step_respa[ilevel];
982
983 // outermost level - update eta_dot and omega_dot, apply to v
984 // all other levels - NVE update of v
985 // x,v updates only performed for atoms in group
986
987 if (ilevel == nlevels_respa-1) {
988
989 // update eta_press_dot
990
991 if (pstat_flag && mpchain) nhc_press_integrate();
992
993 // update eta_dot
994
995 if (tstat_flag) {
996 compute_temp_target();
997 nhc_temp_integrate();
998 }
999
1000 // recompute pressure to account for change in KE
1001 // t_current is up-to-date, but compute_temperature is not
1002 // compute appropriately coupled elements of mvv_current
1003
1004 if (pstat_flag) {
1005 if (pstyle == ISO) {
1006 temperature->compute_scalar();
1007 pressure->compute_scalar();
1008 } else {
1009 temperature->compute_vector();
1010 pressure->compute_vector();
1011 }
1012 couple();
1013 pressure->addstep(update->ntimestep+1);
1014 }
1015
1016 if (pstat_flag) {
1017 compute_press_target();
1018 nh_omega_dot();
1019 nh_v_press();
1020 }
1021
1022 nve_v();
1023
1024 } else nve_v();
1025
1026 // innermost level - also update x only for atoms in group
1027 // if barostat, perform 1/2 step remap before and after
1028
1029 if (ilevel == 0) {
1030 if (pstat_flag) remap();
1031 nve_x();
1032 if (pstat_flag) remap();
1033 }
1034 }
1035
1036 /* ---------------------------------------------------------------------- */
1037
pre_force_respa(int,int ilevel,int)1038 void FixNPTCauchy::pre_force_respa(int /*vflag*/, int ilevel, int /*iloop*/)
1039 {
1040 // if barostat, redo KSpace coeffs at outermost level,
1041 // since volume has changed
1042
1043 if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag)
1044 force->kspace->setup();
1045 }
1046
1047 /* ---------------------------------------------------------------------- */
1048
final_integrate_respa(int ilevel,int)1049 void FixNPTCauchy::final_integrate_respa(int ilevel, int /*iloop*/)
1050 {
1051 // set timesteps by level
1052
1053 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
1054 dthalf = 0.5 * step_respa[ilevel];
1055
1056 // outermost level - update eta_dot and omega_dot, apply via final_integrate
1057 // all other levels - NVE update of v
1058
1059 if (ilevel == nlevels_respa-1) final_integrate();
1060 else nve_v();
1061 }
1062
1063 /* ---------------------------------------------------------------------- */
1064
couple()1065 void FixNPTCauchy::couple()
1066 {
1067 double *tensor = pressure->vector;
1068
1069 if (pstyle == ISO)
1070 p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
1071 else if (pcouple == XYZ) {
1072 double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
1073 p_current[0] = p_current[1] = p_current[2] = ave;
1074 } else if (pcouple == XY) {
1075 double ave = 0.5 * (tensor[0] + tensor[1]);
1076 p_current[0] = p_current[1] = ave;
1077 p_current[2] = tensor[2];
1078 } else if (pcouple == YZ) {
1079 double ave = 0.5 * (tensor[1] + tensor[2]);
1080 p_current[1] = p_current[2] = ave;
1081 p_current[0] = tensor[0];
1082 } else if (pcouple == XZ) {
1083 double ave = 0.5 * (tensor[0] + tensor[2]);
1084 p_current[0] = p_current[2] = ave;
1085 p_current[1] = tensor[1];
1086 } else {
1087 p_current[0] = tensor[0];
1088 p_current[1] = tensor[1];
1089 p_current[2] = tensor[2];
1090 }
1091
1092 if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2]))
1093 error->all(FLERR,"Non-numeric pressure - simulation unstable");
1094
1095 // switch order from xy-xz-yz to Voigt
1096
1097 if (pstyle == TRICLINIC) {
1098 p_current[3] = tensor[5];
1099 p_current[4] = tensor[4];
1100 p_current[5] = tensor[3];
1101
1102 if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5]))
1103 error->all(FLERR,"Non-numeric pressure - simulation unstable");
1104 }
1105 }
1106
1107 /* ----------------------------------------------------------------------
1108 change box size
1109 remap all atoms or dilate group atoms depending on allremap flag
1110 if rigid bodies exist, scale rigid body centers-of-mass
1111 ------------------------------------------------------------------------- */
1112
remap()1113 void FixNPTCauchy::remap()
1114 {
1115 int i;
1116 double oldlo,oldhi;
1117 double expfac;
1118
1119 double **x = atom->x;
1120 int *mask = atom->mask;
1121 int nlocal = atom->nlocal;
1122 double *h = domain->h;
1123
1124 // omega is not used, except for book-keeping
1125
1126 for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
1127
1128 // convert pertinent atoms and rigid bodies to lamda coords
1129
1130 if (allremap) domain->x2lamda(nlocal);
1131 else {
1132 for (i = 0; i < nlocal; i++)
1133 if (mask[i] & dilate_group_bit)
1134 domain->x2lamda(x[i],x[i]);
1135 }
1136
1137 if (nrigid)
1138 for (i = 0; i < nrigid; i++)
1139 modify->fix[rfix[i]]->deform(0);
1140
1141 // reset global and local box to new size/shape
1142
1143 // this operation corresponds to applying the
1144 // translate and scale operations
1145 // corresponding to the solution of the following ODE:
1146 //
1147 // h_dot = omega_dot * h
1148 //
1149 // where h_dot, omega_dot and h are all upper-triangular
1150 // 3x3 tensors. In Voigt notation, the elements of the
1151 // RHS product tensor are:
1152 // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
1153 //
1154 // Ordering of operations preserves time symmetry.
1155
1156 double dto2 = dto/2.0;
1157 double dto4 = dto/4.0;
1158 double dto8 = dto/8.0;
1159
1160 // off-diagonal components, first half
1161
1162 if (pstyle == TRICLINIC) {
1163
1164 if (p_flag[4]) {
1165 expfac = exp(dto8*omega_dot[0]);
1166 h[4] *= expfac;
1167 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1168 h[4] *= expfac;
1169 }
1170
1171 if (p_flag[3]) {
1172 expfac = exp(dto4*omega_dot[1]);
1173 h[3] *= expfac;
1174 h[3] += dto2*(omega_dot[3]*h[2]);
1175 h[3] *= expfac;
1176 }
1177
1178 if (p_flag[5]) {
1179 expfac = exp(dto4*omega_dot[0]);
1180 h[5] *= expfac;
1181 h[5] += dto2*(omega_dot[5]*h[1]);
1182 h[5] *= expfac;
1183 }
1184
1185 if (p_flag[4]) {
1186 expfac = exp(dto8*omega_dot[0]);
1187 h[4] *= expfac;
1188 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1189 h[4] *= expfac;
1190 }
1191 }
1192
1193 // scale diagonal components
1194 // scale tilt factors with cell, if set
1195
1196 if (p_flag[0]) {
1197 oldlo = domain->boxlo[0];
1198 oldhi = domain->boxhi[0];
1199 expfac = exp(dto*omega_dot[0]);
1200 domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
1201 domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
1202 }
1203
1204 if (p_flag[1]) {
1205 oldlo = domain->boxlo[1];
1206 oldhi = domain->boxhi[1];
1207 expfac = exp(dto*omega_dot[1]);
1208 domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
1209 domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
1210 if (scalexy) h[5] *= expfac;
1211 }
1212
1213 if (p_flag[2]) {
1214 oldlo = domain->boxlo[2];
1215 oldhi = domain->boxhi[2];
1216 expfac = exp(dto*omega_dot[2]);
1217 domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
1218 domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
1219 if (scalexz) h[4] *= expfac;
1220 if (scaleyz) h[3] *= expfac;
1221 }
1222
1223 // off-diagonal components, second half
1224
1225 if (pstyle == TRICLINIC) {
1226
1227 if (p_flag[4]) {
1228 expfac = exp(dto8*omega_dot[0]);
1229 h[4] *= expfac;
1230 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1231 h[4] *= expfac;
1232 }
1233
1234 if (p_flag[3]) {
1235 expfac = exp(dto4*omega_dot[1]);
1236 h[3] *= expfac;
1237 h[3] += dto2*(omega_dot[3]*h[2]);
1238 h[3] *= expfac;
1239 }
1240
1241 if (p_flag[5]) {
1242 expfac = exp(dto4*omega_dot[0]);
1243 h[5] *= expfac;
1244 h[5] += dto2*(omega_dot[5]*h[1]);
1245 h[5] *= expfac;
1246 }
1247
1248 if (p_flag[4]) {
1249 expfac = exp(dto8*omega_dot[0]);
1250 h[4] *= expfac;
1251 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1252 h[4] *= expfac;
1253 }
1254
1255 }
1256
1257 domain->yz = h[3];
1258 domain->xz = h[4];
1259 domain->xy = h[5];
1260
1261 // tilt factor to cell length ratio can not exceed TILTMAX in one step
1262
1263 if (domain->yz < -TILTMAX*domain->yprd ||
1264 domain->yz > TILTMAX*domain->yprd ||
1265 domain->xz < -TILTMAX*domain->xprd ||
1266 domain->xz > TILTMAX*domain->xprd ||
1267 domain->xy < -TILTMAX*domain->xprd ||
1268 domain->xy > TILTMAX*domain->xprd)
1269 error->all(FLERR,"Fix npt/cauchy has tilted box too far in one step - "
1270 "periodic cell is too far from equilibrium state");
1271
1272 domain->set_global_box();
1273 domain->set_local_box();
1274
1275 // convert pertinent atoms and rigid bodies back to box coords
1276
1277 if (allremap) domain->lamda2x(nlocal);
1278 else {
1279 for (i = 0; i < nlocal; i++)
1280 if (mask[i] & dilate_group_bit)
1281 domain->lamda2x(x[i],x[i]);
1282 }
1283
1284 if (nrigid)
1285 for (i = 0; i < nrigid; i++)
1286 modify->fix[rfix[i]]->deform(1);
1287 }
1288
1289 /* ----------------------------------------------------------------------
1290 pack entire state of Fix into one write
1291 ------------------------------------------------------------------------- */
1292
write_restart(FILE * fp)1293 void FixNPTCauchy::write_restart(FILE *fp)
1294 {
1295 int nsize = size_restart_global();
1296
1297 double *list;
1298 memory->create(list,nsize,"nh:list");
1299
1300 pack_restart_data(list);
1301
1302 if (comm->me == 0) {
1303 int size = nsize * sizeof(double);
1304 fwrite(&size,sizeof(int),1,fp);
1305 fwrite(list,sizeof(double),nsize,fp);
1306 }
1307
1308 memory->destroy(list);
1309 }
1310
1311 /* ----------------------------------------------------------------------
1312 calculate the number of data to be packed
1313 ------------------------------------------------------------------------- */
1314
size_restart_global()1315 int FixNPTCauchy::size_restart_global()
1316 {
1317 int nsize = 2;
1318 if (tstat_flag) nsize += 1 + 2*mtchain;
1319 if (pstat_flag) {
1320 nsize += 16 + 2*mpchain;
1321 if (deviatoric_flag) nsize += 6;
1322 }
1323
1324 return nsize;
1325 }
1326
1327 /* ----------------------------------------------------------------------
1328 pack restart data
1329 ------------------------------------------------------------------------- */
1330
pack_restart_data(double * list)1331 int FixNPTCauchy::pack_restart_data(double *list)
1332 {
1333 int n = 0;
1334
1335 list[n++] = tstat_flag;
1336 if (tstat_flag) {
1337 list[n++] = mtchain;
1338 for (int ich = 0; ich < mtchain; ich++)
1339 list[n++] = eta[ich];
1340 for (int ich = 0; ich < mtchain; ich++)
1341 list[n++] = eta_dot[ich];
1342 }
1343
1344 list[n++] = pstat_flag;
1345 if (pstat_flag) {
1346 list[n++] = omega[0];
1347 list[n++] = omega[1];
1348 list[n++] = omega[2];
1349 list[n++] = omega[3];
1350 list[n++] = omega[4];
1351 list[n++] = omega[5];
1352 list[n++] = omega_dot[0];
1353 list[n++] = omega_dot[1];
1354 list[n++] = omega_dot[2];
1355 list[n++] = omega_dot[3];
1356 list[n++] = omega_dot[4];
1357 list[n++] = omega_dot[5];
1358 list[n++] = vol0;
1359 list[n++] = t0;
1360 list[n++] = mpchain;
1361 if (mpchain) {
1362 for (int ich = 0; ich < mpchain; ich++)
1363 list[n++] = etap[ich];
1364 for (int ich = 0; ich < mpchain; ich++)
1365 list[n++] = etap_dot[ich];
1366 }
1367
1368 list[n++] = deviatoric_flag;
1369 if (deviatoric_flag) {
1370 list[n++] = h0_inv[0];
1371 list[n++] = h0_inv[1];
1372 list[n++] = h0_inv[2];
1373 list[n++] = h0_inv[3];
1374 list[n++] = h0_inv[4];
1375 list[n++] = h0_inv[5];
1376 }
1377 }
1378
1379 return n;
1380 }
1381
1382 /* ----------------------------------------------------------------------
1383 use state info from restart file to restart the Fix
1384 ------------------------------------------------------------------------- */
1385
restart(char * buf)1386 void FixNPTCauchy::restart(char *buf)
1387 {
1388 int n = 0;
1389 double *list = (double *) buf;
1390 int flag = static_cast<int> (list[n++]);
1391 if (flag) {
1392 int m = static_cast<int> (list[n++]);
1393 if (tstat_flag && m == mtchain) {
1394 for (int ich = 0; ich < mtchain; ich++)
1395 eta[ich] = list[n++];
1396 for (int ich = 0; ich < mtchain; ich++)
1397 eta_dot[ich] = list[n++];
1398 } else n += 2*m;
1399 }
1400 flag = static_cast<int> (list[n++]);
1401 if (flag) {
1402 omega[0] = list[n++];
1403 omega[1] = list[n++];
1404 omega[2] = list[n++];
1405 omega[3] = list[n++];
1406 omega[4] = list[n++];
1407 omega[5] = list[n++];
1408 omega_dot[0] = list[n++];
1409 omega_dot[1] = list[n++];
1410 omega_dot[2] = list[n++];
1411 omega_dot[3] = list[n++];
1412 omega_dot[4] = list[n++];
1413 omega_dot[5] = list[n++];
1414 vol0 = list[n++];
1415 t0 = list[n++];
1416 int m = static_cast<int> (list[n++]);
1417 if (pstat_flag && m == mpchain) {
1418 for (int ich = 0; ich < mpchain; ich++)
1419 etap[ich] = list[n++];
1420 for (int ich = 0; ich < mpchain; ich++)
1421 etap_dot[ich] = list[n++];
1422 } else n+=2*m;
1423 flag = static_cast<int> (list[n++]);
1424 if (flag) {
1425 h0_inv[0] = list[n++];
1426 h0_inv[1] = list[n++];
1427 h0_inv[2] = list[n++];
1428 h0_inv[3] = list[n++];
1429 h0_inv[4] = list[n++];
1430 h0_inv[5] = list[n++];
1431 }
1432 }
1433 }
1434
1435 /* ---------------------------------------------------------------------- */
1436
modify_param(int narg,char ** arg)1437 int FixNPTCauchy::modify_param(int narg, char **arg)
1438 {
1439 if (strcmp(arg[0],"temp") == 0) {
1440 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1441 if (tcomputeflag) {
1442 modify->delete_compute(id_temp);
1443 tcomputeflag = 0;
1444 }
1445 delete [] id_temp;
1446 id_temp = utils::strdup(arg[1]);
1447
1448 int icompute = modify->find_compute(arg[1]);
1449 if (icompute < 0)
1450 error->all(FLERR,"Could not find fix_modify temperature ID");
1451 temperature = modify->compute[icompute];
1452
1453 if (temperature->tempflag == 0)
1454 error->all(FLERR,
1455 "Fix_modify temperature ID does not compute temperature");
1456 if (temperature->igroup != 0 && comm->me == 0)
1457 error->warning(FLERR,"Temperature for fix modify is not for group all");
1458
1459 // reset id_temp of pressure to new temperature ID
1460
1461 if (pstat_flag) {
1462 icompute = modify->find_compute(id_press);
1463 if (icompute < 0)
1464 error->all(FLERR,"Pressure ID for fix modify does not exist");
1465 modify->compute[icompute]->reset_extra_compute_fix(id_temp);
1466 }
1467
1468 return 2;
1469
1470 } else if (strcmp(arg[0],"press") == 0) {
1471 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1472 if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
1473 if (pcomputeflag) {
1474 modify->delete_compute(id_press);
1475 pcomputeflag = 0;
1476 }
1477 delete [] id_press;
1478 id_press = utils::strdup(arg[1]);
1479
1480 int icompute = modify->find_compute(arg[1]);
1481 if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
1482 pressure = modify->compute[icompute];
1483
1484 if (pressure->pressflag == 0)
1485 error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
1486 return 2;
1487 }
1488
1489 return 0;
1490 }
1491
1492 /* ---------------------------------------------------------------------- */
1493
compute_scalar()1494 double FixNPTCauchy::compute_scalar()
1495 {
1496 int i;
1497 double volume;
1498 double energy;
1499 double kt = boltz * t_target;
1500 double lkt_press = 0.0;
1501 int ich;
1502 if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
1503 else volume = domain->xprd * domain->yprd;
1504
1505 energy = 0.0;
1506
1507 // thermostat chain energy is equivalent to Eq. (2) in
1508 // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
1509 // Sum(0.5*p_eta_k^2/Q_k,k=1,M) + L*k*T*eta_1 + Sum(k*T*eta_k,k=2,M),
1510 // where L = tdof
1511 // M = mtchain
1512 // p_eta_k = Q_k*eta_dot[k-1]
1513 // Q_1 = L*k*T/t_freq^2
1514 // Q_k = k*T/t_freq^2, k > 1
1515
1516 if (tstat_flag) {
1517 energy += ke_target * eta[0] + 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
1518 for (ich = 1; ich < mtchain; ich++)
1519 energy += kt * eta[ich] + 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
1520 }
1521
1522 // barostat energy is equivalent to Eq. (8) in
1523 // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
1524 // Sum(0.5*p_omega^2/W + P*V),
1525 // where N = natoms
1526 // p_omega = W*omega_dot
1527 // W = N*k*T/p_freq^2
1528 // sum is over barostatted dimensions
1529
1530 if (pstat_flag) {
1531 for (i = 0; i < 3; i++) {
1532 if (p_flag[i]) {
1533 energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i] +
1534 p_hydro*(volume-vol0) / (pdim*nktv2p);
1535 lkt_press += kt;
1536 }
1537 }
1538
1539 if (pstyle == TRICLINIC) {
1540 for (i = 3; i < 6; i++) {
1541 if (p_flag[i]) {
1542 energy += 0.5*omega_dot[i]*omega_dot[i]*omega_mass[i];
1543 lkt_press += kt;
1544 }
1545 }
1546 }
1547
1548 // extra contributions from thermostat chain for barostat
1549
1550 if (mpchain) {
1551 energy += lkt_press * etap[0] + 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
1552 for (ich = 1; ich < mpchain; ich++)
1553 energy += kt * etap[ich] +
1554 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
1555 }
1556
1557 // extra contribution from strain energy
1558
1559 if (deviatoric_flag) energy += compute_strain_energy();
1560 }
1561
1562 return energy;
1563 }
1564
1565 /* ----------------------------------------------------------------------
1566 return a single element of the following vectors, in this order:
1567 eta[tchain], eta_dot[tchain], omega[ndof], omega_dot[ndof]
1568 etap[pchain], etap_dot[pchain], PE_eta[tchain], KE_eta_dot[tchain]
1569 PE_omega[ndof], KE_omega_dot[ndof], PE_etap[pchain], KE_etap_dot[pchain]
1570 PE_strain[1]
1571 if no thermostat exists, related quantities are omitted from the list
1572 if no barostat exists, related quantities are omitted from the list
1573 ndof = 1,3,6 degrees of freedom for pstyle = ISO,ANISO,TRI
1574 ------------------------------------------------------------------------- */
1575
compute_vector(int n)1576 double FixNPTCauchy::compute_vector(int n)
1577 {
1578 int ilen;
1579
1580 if (tstat_flag) {
1581 ilen = mtchain;
1582 if (n < ilen) return eta[n];
1583 n -= ilen;
1584 ilen = mtchain;
1585 if (n < ilen) return eta_dot[n];
1586 n -= ilen;
1587 }
1588
1589 if (pstat_flag) {
1590 if (pstyle == ISO) {
1591 ilen = 1;
1592 if (n < ilen) return omega[n];
1593 n -= ilen;
1594 } else if (pstyle == ANISO) {
1595 ilen = 3;
1596 if (n < ilen) return omega[n];
1597 n -= ilen;
1598 } else {
1599 ilen = 6;
1600 if (n < ilen) return omega[n];
1601 n -= ilen;
1602 }
1603
1604 if (pstyle == ISO) {
1605 ilen = 1;
1606 if (n < ilen) return omega_dot[n];
1607 n -= ilen;
1608 } else if (pstyle == ANISO) {
1609 ilen = 3;
1610 if (n < ilen) return omega_dot[n];
1611 n -= ilen;
1612 } else {
1613 ilen = 6;
1614 if (n < ilen) return omega_dot[n];
1615 n -= ilen;
1616 }
1617
1618 if (mpchain) {
1619 ilen = mpchain;
1620 if (n < ilen) return etap[n];
1621 n -= ilen;
1622 ilen = mpchain;
1623 if (n < ilen) return etap_dot[n];
1624 n -= ilen;
1625 }
1626 }
1627
1628 double volume;
1629 double kt = boltz * t_target;
1630 double lkt_press = kt;
1631 int ich;
1632 if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
1633 else volume = domain->xprd * domain->yprd;
1634
1635 if (tstat_flag) {
1636 ilen = mtchain;
1637 if (n < ilen) {
1638 ich = n;
1639 if (ich == 0)
1640 return ke_target * eta[0];
1641 else
1642 return kt * eta[ich];
1643 }
1644 n -= ilen;
1645 ilen = mtchain;
1646 if (n < ilen) {
1647 ich = n;
1648 if (ich == 0)
1649 return 0.5*eta_mass[0]*eta_dot[0]*eta_dot[0];
1650 else
1651 return 0.5*eta_mass[ich]*eta_dot[ich]*eta_dot[ich];
1652 }
1653 n -= ilen;
1654 }
1655
1656 if (pstat_flag) {
1657 if (pstyle == ISO) {
1658 ilen = 1;
1659 if (n < ilen)
1660 return p_hydro*(volume-vol0) / nktv2p;
1661 n -= ilen;
1662 } else if (pstyle == ANISO) {
1663 ilen = 3;
1664 if (n < ilen) {
1665 if (p_flag[n])
1666 return p_hydro*(volume-vol0) / (pdim*nktv2p);
1667 else
1668 return 0.0;
1669 }
1670 n -= ilen;
1671 } else {
1672 ilen = 6;
1673 if (n < ilen) {
1674 if (n > 2) return 0.0;
1675 else if (p_flag[n])
1676 return p_hydro*(volume-vol0) / (pdim*nktv2p);
1677 else
1678 return 0.0;
1679 }
1680 n -= ilen;
1681 }
1682
1683 if (pstyle == ISO) {
1684 ilen = 1;
1685 if (n < ilen)
1686 return pdim*0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
1687 n -= ilen;
1688 } else if (pstyle == ANISO) {
1689 ilen = 3;
1690 if (n < ilen) {
1691 if (p_flag[n])
1692 return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
1693 else return 0.0;
1694 }
1695 n -= ilen;
1696 } else {
1697 ilen = 6;
1698 if (n < ilen) {
1699 if (p_flag[n])
1700 return 0.5*omega_dot[n]*omega_dot[n]*omega_mass[n];
1701 else return 0.0;
1702 }
1703 n -= ilen;
1704 }
1705
1706 if (mpchain) {
1707 ilen = mpchain;
1708 if (n < ilen) {
1709 ich = n;
1710 if (ich == 0) return lkt_press * etap[0];
1711 else return kt * etap[ich];
1712 }
1713 n -= ilen;
1714 ilen = mpchain;
1715 if (n < ilen) {
1716 ich = n;
1717 if (ich == 0)
1718 return 0.5*etap_mass[0]*etap_dot[0]*etap_dot[0];
1719 else
1720 return 0.5*etap_mass[ich]*etap_dot[ich]*etap_dot[ich];
1721 }
1722 n -= ilen;
1723 }
1724
1725 if (deviatoric_flag) {
1726 ilen = 1;
1727 if (n < ilen)
1728 return compute_strain_energy();
1729 n -= ilen;
1730 }
1731 }
1732
1733 return 0.0;
1734 }
1735
1736 /* ---------------------------------------------------------------------- */
1737
reset_target(double t_new)1738 void FixNPTCauchy::reset_target(double t_new)
1739 {
1740 t_target = t_start = t_stop = t_new;
1741 }
1742
1743 /* ---------------------------------------------------------------------- */
1744
reset_dt()1745 void FixNPTCauchy::reset_dt()
1746 {
1747 dtv = update->dt;
1748 dtf = 0.5 * update->dt * force->ftm2v;
1749 dthalf = 0.5 * update->dt;
1750 dt4 = 0.25 * update->dt;
1751 dt8 = 0.125 * update->dt;
1752 dto = dthalf;
1753
1754 // If using respa, then remap is performed in innermost level
1755
1756 if (utils::strmatch(update->integrate_style,"^respa"))
1757 dto = 0.5*step_respa[0];
1758
1759 if (pstat_flag)
1760 pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
1761
1762 if (tstat_flag)
1763 tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
1764 }
1765
1766 /* ----------------------------------------------------------------------
1767 extract thermostat properties
1768 ------------------------------------------------------------------------- */
1769
extract(const char * str,int & dim)1770 void *FixNPTCauchy::extract(const char *str, int &dim)
1771 {
1772 dim=0;
1773 if (tstat_flag && strcmp(str,"t_target") == 0) {
1774 return &t_target;
1775 } else if (tstat_flag && strcmp(str,"t_start") == 0) {
1776 return &t_start;
1777 } else if (tstat_flag && strcmp(str,"t_stop") == 0) {
1778 return &t_stop;
1779 } else if (tstat_flag && strcmp(str,"mtchain") == 0) {
1780 return &mtchain;
1781 } else if (pstat_flag && strcmp(str,"mpchain") == 0) {
1782 return &mtchain;
1783 }
1784 dim=1;
1785 if (tstat_flag && strcmp(str,"eta") == 0) {
1786 return η
1787 } else if (pstat_flag && strcmp(str,"etap") == 0) {
1788 return η
1789 } else if (pstat_flag && strcmp(str,"p_flag") == 0) {
1790 return &p_flag;
1791 } else if (pstat_flag && strcmp(str,"p_start") == 0) {
1792 return &p_start;
1793 } else if (pstat_flag && strcmp(str,"p_stop") == 0) {
1794 return &p_stop;
1795 } else if (pstat_flag && strcmp(str,"p_target") == 0) {
1796 return &p_target;
1797 }
1798 return nullptr;
1799 }
1800
1801 /* ----------------------------------------------------------------------
1802 perform half-step update of chain thermostat variables
1803 ------------------------------------------------------------------------- */
1804
nhc_temp_integrate()1805 void FixNPTCauchy::nhc_temp_integrate()
1806 {
1807 int ich;
1808 double expfac;
1809 double kecurrent = tdof * boltz * t_current;
1810
1811 // Update masses, to preserve initial freq, if flag set
1812
1813 if (eta_mass_flag) {
1814 eta_mass[0] = tdof * boltz * t_target / (t_freq*t_freq);
1815 for (int ich = 1; ich < mtchain; ich++)
1816 eta_mass[ich] = boltz * t_target / (t_freq*t_freq);
1817 }
1818
1819 if (eta_mass[0] > 0.0)
1820 eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
1821 else eta_dotdot[0] = 0.0;
1822
1823 double ncfac = 1.0/nc_tchain;
1824 for (int iloop = 0; iloop < nc_tchain; iloop++) {
1825
1826 for (ich = mtchain-1; ich > 0; ich--) {
1827 expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
1828 eta_dot[ich] *= expfac;
1829 eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
1830 eta_dot[ich] *= tdrag_factor;
1831 eta_dot[ich] *= expfac;
1832 }
1833
1834 expfac = exp(-ncfac*dt8*eta_dot[1]);
1835 eta_dot[0] *= expfac;
1836 eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
1837 eta_dot[0] *= tdrag_factor;
1838 eta_dot[0] *= expfac;
1839
1840 factor_eta = exp(-ncfac*dthalf*eta_dot[0]);
1841 nh_v_temp();
1842
1843 // rescale temperature due to velocity scaling
1844 // should not be necessary to explicitly recompute the temperature
1845
1846 t_current *= factor_eta*factor_eta;
1847 kecurrent = tdof * boltz * t_current;
1848
1849 if (eta_mass[0] > 0.0)
1850 eta_dotdot[0] = (kecurrent - ke_target)/eta_mass[0];
1851 else eta_dotdot[0] = 0.0;
1852
1853 for (ich = 0; ich < mtchain; ich++)
1854 eta[ich] += ncfac*dthalf*eta_dot[ich];
1855
1856 eta_dot[0] *= expfac;
1857 eta_dot[0] += eta_dotdot[0] * ncfac*dt4;
1858 eta_dot[0] *= expfac;
1859
1860 for (ich = 1; ich < mtchain; ich++) {
1861 expfac = exp(-ncfac*dt8*eta_dot[ich+1]);
1862 eta_dot[ich] *= expfac;
1863 eta_dotdot[ich] = (eta_mass[ich-1]*eta_dot[ich-1]*eta_dot[ich-1]
1864 - boltz * t_target)/eta_mass[ich];
1865 eta_dot[ich] += eta_dotdot[ich] * ncfac*dt4;
1866 eta_dot[ich] *= expfac;
1867 }
1868 }
1869 }
1870
1871 /* ----------------------------------------------------------------------
1872 perform half-step update of chain thermostat variables for barostat
1873 scale barostat velocities
1874 ------------------------------------------------------------------------- */
1875
nhc_press_integrate()1876 void FixNPTCauchy::nhc_press_integrate()
1877 {
1878 int ich,i,pdof;
1879 double expfac,factor_etap,kecurrent;
1880 double kt = boltz * t_target;
1881 double lkt_press;
1882
1883 // Update masses, to preserve initial freq, if flag set
1884
1885 if (omega_mass_flag) {
1886 double nkt = (atom->natoms + 1) * kt;
1887 for (int i = 0; i < 3; i++)
1888 if (p_flag[i])
1889 omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
1890
1891 if (pstyle == TRICLINIC) {
1892 for (int i = 3; i < 6; i++)
1893 if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
1894 }
1895 }
1896
1897 if (etap_mass_flag) {
1898 if (mpchain) {
1899 etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
1900 for (int ich = 1; ich < mpchain; ich++)
1901 etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
1902 for (int ich = 1; ich < mpchain; ich++)
1903 etap_dotdot[ich] =
1904 (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
1905 boltz * t_target) / etap_mass[ich];
1906 }
1907 }
1908
1909 kecurrent = 0.0;
1910 pdof = 0;
1911 for (i = 0; i < 3; i++)
1912 if (p_flag[i]) {
1913 kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1914 pdof++;
1915 }
1916
1917 if (pstyle == TRICLINIC) {
1918 for (i = 3; i < 6; i++)
1919 if (p_flag[i]) {
1920 kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1921 pdof++;
1922 }
1923 }
1924
1925 lkt_press = pdof * kt;
1926 etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
1927
1928 double ncfac = 1.0/nc_pchain;
1929 for (int iloop = 0; iloop < nc_pchain; iloop++) {
1930
1931 for (ich = mpchain-1; ich > 0; ich--) {
1932 expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
1933 etap_dot[ich] *= expfac;
1934 etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
1935 etap_dot[ich] *= pdrag_factor;
1936 etap_dot[ich] *= expfac;
1937 }
1938
1939 expfac = exp(-ncfac*dt8*etap_dot[1]);
1940 etap_dot[0] *= expfac;
1941 etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
1942 etap_dot[0] *= pdrag_factor;
1943 etap_dot[0] *= expfac;
1944
1945 for (ich = 0; ich < mpchain; ich++)
1946 etap[ich] += ncfac*dthalf*etap_dot[ich];
1947
1948 factor_etap = exp(-ncfac*dthalf*etap_dot[0]);
1949 for (i = 0; i < 3; i++)
1950 if (p_flag[i]) omega_dot[i] *= factor_etap;
1951
1952 if (pstyle == TRICLINIC) {
1953 for (i = 3; i < 6; i++)
1954 if (p_flag[i]) omega_dot[i] *= factor_etap;
1955 }
1956
1957 kecurrent = 0.0;
1958 for (i = 0; i < 3; i++)
1959 if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1960
1961 if (pstyle == TRICLINIC) {
1962 for (i = 3; i < 6; i++)
1963 if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1964 }
1965
1966 etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
1967
1968 etap_dot[0] *= expfac;
1969 etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
1970 etap_dot[0] *= expfac;
1971
1972 for (ich = 1; ich < mpchain; ich++) {
1973 expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
1974 etap_dot[ich] *= expfac;
1975 etap_dotdot[ich] =
1976 (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - boltz*t_target) /
1977 etap_mass[ich];
1978 etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
1979 etap_dot[ich] *= expfac;
1980 }
1981 }
1982 }
1983
1984 /* ----------------------------------------------------------------------
1985 perform half-step barostat scaling of velocities
1986 -----------------------------------------------------------------------*/
1987
nh_v_press()1988 void FixNPTCauchy::nh_v_press()
1989 {
1990 double factor[3];
1991 double **v = atom->v;
1992 int *mask = atom->mask;
1993 int nlocal = atom->nlocal;
1994 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
1995
1996 factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2));
1997 factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2));
1998 factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2));
1999
2000 if (which == NOBIAS) {
2001 for (int i = 0; i < nlocal; i++) {
2002 if (mask[i] & groupbit) {
2003 v[i][0] *= factor[0];
2004 v[i][1] *= factor[1];
2005 v[i][2] *= factor[2];
2006 if (pstyle == TRICLINIC) {
2007 v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
2008 v[i][1] += -dthalf*v[i][2]*omega_dot[3];
2009 }
2010 v[i][0] *= factor[0];
2011 v[i][1] *= factor[1];
2012 v[i][2] *= factor[2];
2013 }
2014 }
2015 } else if (which == BIAS) {
2016 for (int i = 0; i < nlocal; i++) {
2017 if (mask[i] & groupbit) {
2018 temperature->remove_bias(i,v[i]);
2019 v[i][0] *= factor[0];
2020 v[i][1] *= factor[1];
2021 v[i][2] *= factor[2];
2022 if (pstyle == TRICLINIC) {
2023 v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
2024 v[i][1] += -dthalf*v[i][2]*omega_dot[3];
2025 }
2026 v[i][0] *= factor[0];
2027 v[i][1] *= factor[1];
2028 v[i][2] *= factor[2];
2029 temperature->restore_bias(i,v[i]);
2030 }
2031 }
2032 }
2033 }
2034
2035 /* ----------------------------------------------------------------------
2036 perform half-step update of velocities
2037 -----------------------------------------------------------------------*/
2038
nve_v()2039 void FixNPTCauchy::nve_v()
2040 {
2041 double dtfm;
2042 double **v = atom->v;
2043 double **f = atom->f;
2044 double *rmass = atom->rmass;
2045 double *mass = atom->mass;
2046 int *type = atom->type;
2047 int *mask = atom->mask;
2048 int nlocal = atom->nlocal;
2049 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
2050
2051 if (rmass) {
2052 for (int i = 0; i < nlocal; i++) {
2053 if (mask[i] & groupbit) {
2054 dtfm = dtf / rmass[i];
2055 v[i][0] += dtfm*f[i][0];
2056 v[i][1] += dtfm*f[i][1];
2057 v[i][2] += dtfm*f[i][2];
2058 }
2059 }
2060 } else {
2061 for (int i = 0; i < nlocal; i++) {
2062 if (mask[i] & groupbit) {
2063 dtfm = dtf / mass[type[i]];
2064 v[i][0] += dtfm*f[i][0];
2065 v[i][1] += dtfm*f[i][1];
2066 v[i][2] += dtfm*f[i][2];
2067 }
2068 }
2069 }
2070 }
2071
2072 /* ----------------------------------------------------------------------
2073 perform full-step update of positions
2074 -----------------------------------------------------------------------*/
2075
nve_x()2076 void FixNPTCauchy::nve_x()
2077 {
2078 double **x = atom->x;
2079 double **v = atom->v;
2080 int *mask = atom->mask;
2081 int nlocal = atom->nlocal;
2082 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
2083
2084 // x update by full step only for atoms in group
2085
2086 for (int i = 0; i < nlocal; i++) {
2087 if (mask[i] & groupbit) {
2088 x[i][0] += dtv * v[i][0];
2089 x[i][1] += dtv * v[i][1];
2090 x[i][2] += dtv * v[i][2];
2091 }
2092 }
2093 }
2094
2095 /* ----------------------------------------------------------------------
2096 perform half-step thermostat scaling of velocities
2097 -----------------------------------------------------------------------*/
2098
nh_v_temp()2099 void FixNPTCauchy::nh_v_temp()
2100 {
2101 double **v = atom->v;
2102 int *mask = atom->mask;
2103 int nlocal = atom->nlocal;
2104 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
2105
2106 if (which == NOBIAS) {
2107 for (int i = 0; i < nlocal; i++) {
2108 if (mask[i] & groupbit) {
2109 v[i][0] *= factor_eta;
2110 v[i][1] *= factor_eta;
2111 v[i][2] *= factor_eta;
2112 }
2113 }
2114 } else if (which == BIAS) {
2115 for (int i = 0; i < nlocal; i++) {
2116 if (mask[i] & groupbit) {
2117 temperature->remove_bias(i,v[i]);
2118 v[i][0] *= factor_eta;
2119 v[i][1] *= factor_eta;
2120 v[i][2] *= factor_eta;
2121 temperature->restore_bias(i,v[i]);
2122 }
2123 }
2124 }
2125 }
2126
2127 /* ----------------------------------------------------------------------
2128 compute sigma tensor
2129 needed whenever p_target or h0_inv changes
2130 -----------------------------------------------------------------------*/
2131
compute_sigma()2132 void FixNPTCauchy::compute_sigma()
2133 {
2134 // if nreset_h0 > 0, reset vol0 and h0_inv
2135 // every nreset_h0 timesteps
2136
2137 if (nreset_h0 > 0) {
2138 int delta = update->ntimestep - update->beginstep;
2139 if (delta % nreset_h0 == 0) {
2140 if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
2141 else vol0 = domain->xprd * domain->yprd;
2142 h0_inv[0] = domain->h_inv[0];
2143 h0_inv[1] = domain->h_inv[1];
2144 h0_inv[2] = domain->h_inv[2];
2145 h0_inv[3] = domain->h_inv[3];
2146 h0_inv[4] = domain->h_inv[4];
2147 h0_inv[5] = domain->h_inv[5];
2148 }
2149 }
2150
2151 // generate upper-triangular half of
2152 // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t
2153 // units of sigma are are PV/L^2 e.g. atm.A
2154 //
2155 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
2156 // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
2157 // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
2158
2159 sigma[0] =
2160 vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] +
2161 p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) +
2162 h0_inv[5]*(p_target[5]*h0_inv[0] +
2163 (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) +
2164 h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] +
2165 (p_target[2]-p_hydro)*h0_inv[4]));
2166 sigma[1] =
2167 vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] +
2168 p_target[3]*h0_inv[3]) +
2169 h0_inv[3]*(p_target[3]*h0_inv[1] +
2170 (p_target[2]-p_hydro)*h0_inv[3]));
2171 sigma[2] =
2172 vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2]));
2173 sigma[3] =
2174 vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) +
2175 h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2]));
2176 sigma[4] =
2177 vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) +
2178 h0_inv[5]*(p_target[3]*h0_inv[2]) +
2179 h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2]));
2180 sigma[5] =
2181 vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) +
2182 h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) +
2183 h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3]));
2184 }
2185
2186 /* ----------------------------------------------------------------------
2187 compute strain energy
2188 -----------------------------------------------------------------------*/
2189
compute_strain_energy()2190 double FixNPTCauchy::compute_strain_energy()
2191 {
2192 // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units
2193
2194 double* h = domain->h;
2195 double d0,d1,d2;
2196
2197 d0 =
2198 sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) +
2199 sigma[5]*( h[1]*h[5]+h[3]*h[4]) +
2200 sigma[4]*( h[2]*h[4]);
2201 d1 =
2202 sigma[5]*( h[5]*h[1]+h[4]*h[3]) +
2203 sigma[1]*( h[1]*h[1]+h[3]*h[3]) +
2204 sigma[3]*( h[2]*h[3]);
2205 d2 =
2206 sigma[4]*( h[4]*h[2]) +
2207 sigma[3]*( h[3]*h[2]) +
2208 sigma[2]*( h[2]*h[2]);
2209
2210 double energy = 0.5*(d0+d1+d2)/nktv2p;
2211 return energy;
2212 }
2213
2214 /* ----------------------------------------------------------------------
2215 compute deviatoric barostat force = h*sigma*h^t
2216 -----------------------------------------------------------------------*/
2217
compute_deviatoric()2218 void FixNPTCauchy::compute_deviatoric()
2219 {
2220 // generate upper-triangular part of h*sigma*h^t
2221 // units of fdev are are PV, e.g. atm*A^3
2222 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
2223 // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
2224 // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
2225
2226 double* h = domain->h;
2227
2228 fdev[0] =
2229 h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) +
2230 h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) +
2231 h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]);
2232 fdev[1] =
2233 h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) +
2234 h[3]*( sigma[3]*h[1]+sigma[2]*h[3]);
2235 fdev[2] =
2236 h[2]*( sigma[2]*h[2]);
2237 fdev[3] =
2238 h[1]*( sigma[3]*h[2]) +
2239 h[3]*( sigma[2]*h[2]);
2240 fdev[4] =
2241 h[0]*( sigma[4]*h[2]) +
2242 h[5]*( sigma[3]*h[2]) +
2243 h[4]*( sigma[2]*h[2]);
2244 fdev[5] =
2245 h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) +
2246 h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) +
2247 h[4]*( sigma[3]*h[1]+sigma[2]*h[3]);
2248 }
2249
2250 /* ----------------------------------------------------------------------
2251 compute target temperature and kinetic energy
2252 -----------------------------------------------------------------------*/
2253
compute_temp_target()2254 void FixNPTCauchy::compute_temp_target()
2255 {
2256 double delta = update->ntimestep - update->beginstep;
2257 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2258
2259 t_target = t_start + delta * (t_stop-t_start);
2260 ke_target = tdof * boltz * t_target;
2261 }
2262
2263 /* ----------------------------------------------------------------------
2264 compute hydrostatic target pressure
2265 -----------------------------------------------------------------------*/
2266
compute_press_target()2267 void FixNPTCauchy::compute_press_target()
2268 {
2269 double delta = update->ntimestep - update->beginstep;
2270 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2271
2272 p_hydro = 0.0;
2273 for (int i = 0; i < 3; i++)
2274 if (p_flag[i]) {
2275 p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
2276 p_hydro += p_target[i];
2277 }
2278 if (pdim > 0) p_hydro /= pdim;
2279
2280 if (pstyle == TRICLINIC)
2281 for (int i = 3; i < 6; i++)
2282 p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
2283
2284 // CauchyStat: call CauchyStat to modify p_target[i] and p_hydro,
2285 // if CauchyStat enabled and pressure->vector computation has been initiated
2286
2287 if (initRUN == 1) CauchyStat();
2288 if (initRUN == 0) {
2289 for (int i=0; i < 6; i++) {
2290 h_old[i]=domain->h[i];
2291 }
2292 }
2293
2294 // when run is initialized tensor[] not available (pressure on cell wall)
2295
2296 initRUN=1;
2297
2298 // if deviatoric, recompute sigma each time p_target changes
2299
2300 if (deviatoric_flag) compute_sigma();
2301 }
2302
2303 /* ----------------------------------------------------------------------
2304 update omega_dot, omega
2305 -----------------------------------------------------------------------*/
2306
nh_omega_dot()2307 void FixNPTCauchy::nh_omega_dot()
2308 {
2309 double f_omega,volume;
2310
2311 if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
2312 else volume = domain->xprd*domain->yprd;
2313
2314 if (deviatoric_flag) compute_deviatoric();
2315
2316 mtk_term1 = 0.0;
2317 if (mtk_flag) {
2318 if (pstyle == ISO) {
2319 mtk_term1 = tdof * boltz * t_current;
2320 mtk_term1 /= pdim * atom->natoms;
2321 } else {
2322 double *mvv_current = temperature->vector;
2323 for (int i = 0; i < 3; i++)
2324 if (p_flag[i])
2325 mtk_term1 += mvv_current[i];
2326 mtk_term1 /= pdim * atom->natoms;
2327 }
2328 }
2329
2330 for (int i = 0; i < 3; i++)
2331 if (p_flag[i]) {
2332 f_omega = (p_current[i]-p_hydro)*volume /
2333 (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i];
2334 if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
2335 omega_dot[i] += f_omega*dthalf;
2336 omega_dot[i] *= pdrag_factor;
2337 }
2338
2339 mtk_term2 = 0.0;
2340 if (mtk_flag) {
2341 for (int i = 0; i < 3; i++)
2342 if (p_flag[i])
2343 mtk_term2 += omega_dot[i];
2344 if (pdim > 0) mtk_term2 /= pdim * atom->natoms;
2345 }
2346
2347 if (pstyle == TRICLINIC) {
2348 for (int i = 3; i < 6; i++) {
2349 if (p_flag[i]) {
2350 f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p);
2351 if (deviatoric_flag)
2352 f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
2353 omega_dot[i] += f_omega*dthalf;
2354 omega_dot[i] *= pdrag_factor;
2355 }
2356 }
2357 }
2358 }
2359
2360 /* ----------------------------------------------------------------------
2361 if any tilt ratios exceed limits, set flip = 1 and compute new tilt values
2362 do not flip in x or y if non-periodic (can tilt but not flip)
2363 this is b/c the box length would be changed (dramatically) by flip
2364 if yz tilt exceeded, adjust C vector by one B vector
2365 if xz tilt exceeded, adjust C vector by one A vector
2366 if xy tilt exceeded, adjust B vector by one A vector
2367 check yz first since it may change xz, then xz check comes after
2368 if any flip occurs, create new box in domain
2369 image_flip() adjusts image flags due to box shape change induced by flip
2370 remap() puts atoms outside the new box back into the new box
2371 perform irregular on atoms in lamda coords to migrate atoms to new procs
2372 important that image_flip comes before remap, since remap may change
2373 image flags to new values, making eqs in doc of Domain:image_flip incorrect
2374 ------------------------------------------------------------------------- */
2375
pre_exchange()2376 void FixNPTCauchy::pre_exchange()
2377 {
2378 double xprd = domain->xprd;
2379 double yprd = domain->yprd;
2380
2381 // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP
2382 // this avoids immediate re-flipping due to tilt oscillations
2383
2384 double xtiltmax = (0.5+DELTAFLIP)*xprd;
2385 double ytiltmax = (0.5+DELTAFLIP)*yprd;
2386
2387 int flipxy,flipxz,flipyz;
2388 flipxy = flipxz = flipyz = 0;
2389
2390 if (domain->yperiodic) {
2391 if (domain->yz < -ytiltmax) {
2392 domain->yz += yprd;
2393 domain->xz += domain->xy;
2394 flipyz = 1;
2395 } else if (domain->yz >= ytiltmax) {
2396 domain->yz -= yprd;
2397 domain->xz -= domain->xy;
2398 flipyz = -1;
2399 }
2400 }
2401
2402 if (domain->xperiodic) {
2403 if (domain->xz < -xtiltmax) {
2404 domain->xz += xprd;
2405 flipxz = 1;
2406 } else if (domain->xz >= xtiltmax) {
2407 domain->xz -= xprd;
2408 flipxz = -1;
2409 }
2410 if (domain->xy < -xtiltmax) {
2411 domain->xy += xprd;
2412 flipxy = 1;
2413 } else if (domain->xy >= xtiltmax) {
2414 domain->xy -= xprd;
2415 flipxy = -1;
2416 }
2417 }
2418
2419 int flip = 0;
2420 if (flipxy || flipxz || flipyz) flip = 1;
2421
2422 if (flip) {
2423 domain->set_global_box();
2424 domain->set_local_box();
2425
2426 domain->image_flip(flipxy,flipxz,flipyz);
2427
2428 double **x = atom->x;
2429 imageint *image = atom->image;
2430 int nlocal = atom->nlocal;
2431 for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
2432
2433 domain->x2lamda(atom->nlocal);
2434 irregular->migrate_atoms();
2435 domain->lamda2x(atom->nlocal);
2436 }
2437 }
2438
2439 /* ----------------------------------------------------------------------
2440 memory usage of Irregular
2441 ------------------------------------------------------------------------- */
2442
memory_usage()2443 double FixNPTCauchy::memory_usage()
2444 {
2445 double bytes = 0.0;
2446 if (irregular) bytes += irregular->memory_usage();
2447 return bytes;
2448 }
2449
2450 /* ----------------------------------------------------------------------
2451 initialize Cauchystat
2452 ------------------------------------------------------------------------- */
2453
CauchyStat_init()2454 void FixNPTCauchy::CauchyStat_init()
2455 {
2456 if (comm->me == 0) {
2457 std::string mesg = fmt::format("Using fix npt/cauchy with alpha={:.8f}\n",alpha);
2458 if (restartPK==1) {
2459 mesg += " (this is a continuation run)\n";
2460 } else {
2461 mesg += " (this is NOT a continuation run)\n";
2462 }
2463 utils::logmesg(lmp, mesg);
2464 }
2465
2466 if (!id_store)
2467 id_store = utils::strdup(std::string(id) + "_FIX_NH_STORE");
2468 restart_stored = modify->find_fix(id_store);
2469
2470 if (restartPK==1 && restart_stored < 0)
2471 error->all(FLERR,"Illegal npt/cauchy command. Continuation run"
2472 " must follow a previously equilibrated npt/cauchy run");
2473
2474 if (alpha<=0.0)
2475 error->all(FLERR,"Illegal fix npt/cauchy command: Alpha cannot be zero or negative.");
2476
2477 if (restart_stored < 0) {
2478 modify->add_fix(std::string(id_store) + " all STORE global 1 6");
2479 restart_stored = modify->find_fix(id_store);
2480 }
2481 init_store = (FixStore *)modify->fix[restart_stored];
2482
2483 initRUN = 0;
2484 initPK = 1;
2485
2486 #define H0(row,col) (H0[(row-1)][(col-1)])
2487 #define invH0(row,col) (invH0[(row-1)][(col-1)])
2488
2489 // initialize H0 to current box shape
2490
2491 double *h = domain->h;
2492 double *h_inv = domain->h_inv;
2493
2494 H0(1,1)=h[0]; H0(1,2)=h[5]; H0(1,3)=h[4];
2495 H0(2,1)=0.0; H0(2,2)=h[1]; H0(2,3)=h[3];
2496 H0(3,1)=0.0; H0(3,2)=0.0; H0(3,3)=h[2];
2497
2498 invH0(1,1)=h_inv[0]; invH0(1,2)=h_inv[5]; invH0(1,3)=h_inv[4];
2499 invH0(2,1)=0.0; invH0(2,2)=h_inv[1]; invH0(2,3)=h_inv[3];
2500 invH0(3,1)=0.0; invH0(3,2)=0.0; invH0(3,3)=h_inv[2];
2501
2502 CSvol0 = fabs(MathExtra::det3(H0)); //find reference volume
2503
2504 #undef H0
2505 #undef invH0
2506 }
2507
2508 /* ----------------------------------------------------------------------
2509 Cauchystat
2510 ------------------------------------------------------------------------- */
2511
CauchyStat()2512 void FixNPTCauchy::CauchyStat()
2513 {
2514 double* h = domain->h; // shape matrix in Voigt notation
2515 double* h_rate = domain->h_rate; // rate of box size/shape change in Voigt notation
2516 double H[3][3];
2517 double Hdot[3][3];
2518
2519 #define H(row,col) (H[(row-1)][(col-1)])
2520 #define Hdot(row,col) (Hdot[(row-1)][(col-1)])
2521 #define F(row,col) (F[(row-1)][(col-1)])
2522 #define Fi(row,col) (Fi[(row-1)][(col-1)])
2523 #define Fdot(row,col) (Fdot[(row-1)][(col-1)])
2524 #define invH0(row,col) (invH0[(row-1)][(col-1)])
2525 #define cauchy(row,col) (cauchy[(row-1)][(col-1)])
2526 #define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
2527 #define setPK(row,col) (setPK[(row-1)][(col-1)])
2528 #define sigmahat(row,col) (sigmahat[(row-1)][(col-1)])
2529
2530 H(1,1)=h[0]; H(1,2)=0.0; H(1,3)=0.0;
2531 H(2,1)=0.0; H(2,2)=h[1]; H(2,3)=0.0;
2532 H(3,1)=0.0; H(3,2)=0.0; H(3,3)=h[2];
2533
2534 for (int i=0;i<6;i++) {
2535 h_rate[i]=(h[i]-h_old[i])/update->dt;
2536 h_old[i]=h[i];
2537 }
2538
2539 Hdot(1,1)=h_rate[0]; Hdot(1,2)=0.0; Hdot(1,3)=0.0;
2540 Hdot(2,1)=0.0; Hdot(2,2)=h_rate[1]; Hdot(2,3)=0.0;
2541 Hdot(3,1)=0.0; Hdot(3,2)=0.0; Hdot(3,3)=h_rate[2];
2542
2543 if (domain->triclinic) {
2544 H(1,2)=h[5]; H(1,3)=h[4]; H(2,3)=h[3];
2545 Hdot(1,2)=h_rate[5]; Hdot(1,3)=h_rate[4]; Hdot(2,3)=h_rate[3];
2546 }
2547
2548 double F[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
2549 double Fi[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
2550 double Fdot[3][3]={{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
2551 double J,vol;
2552
2553 MathExtra::times3(H,invH0,F);
2554 MathExtra::times3(Hdot,invH0,Fdot);
2555 MathExtra::invert3(F,Fi);
2556 J = MathExtra::det3(F);
2557 vol=CSvol0*J;
2558
2559 double deltat;
2560 deltat = update->dt; //increment of time step
2561
2562 // Current pressure on the cell boundaries:
2563 //tensor: 0 1 2 3 4 5
2564 // x y z xy xz yz
2565
2566 double *tensor = pressure->vector;
2567
2568 double cauchy[3][3];
2569 cauchy(1,1)=-tensor[0]; cauchy(1,2)=0.0; cauchy(1,3)=0.0;
2570 cauchy(2,1)=0.0; cauchy(2,2)=-tensor[1]; cauchy(2,3)=0.0;
2571 cauchy(3,1)=0.0; cauchy(3,2)=0.0; cauchy(3,3)=-tensor[2];
2572
2573 if (domain->triclinic) {
2574 cauchy(1,2)=-tensor[3]; cauchy(1,3)=-tensor[4];
2575 cauchy(2,1)=-tensor[3]; cauchy(2,3)=-tensor[5];
2576 cauchy(3,1)=-tensor[4]; cauchy(3,2)=-tensor[5];
2577 }
2578
2579 // target pressure on the cell boundaries:
2580 // p_target: 0 1 2 3 4 5
2581 // x y z yz xz xy
2582
2583 double setcauchy[3][3];
2584 setcauchy(1,1)=-p_target[0]; setcauchy(1,2)=0.0; setcauchy(1,3)=0.0;
2585 setcauchy(2,1)=0.0; setcauchy(2,2)=-p_target[1]; setcauchy(2,3)=0.0;
2586 setcauchy(3,1)=0.0; setcauchy(3,2)=0.0; setcauchy(3,3)=-p_target[2];
2587
2588 if (domain->triclinic) {
2589 setcauchy(1,2)=-p_target[5]; setcauchy(1,3)=-p_target[4];
2590 setcauchy(2,1)=-p_target[5]; setcauchy(2,3)=-p_target[3];
2591 setcauchy(3,1)=-p_target[4]; setcauchy(3,2)=-p_target[3];
2592 }
2593
2594 //initialize:
2595 if (initPK==1) {
2596 if (restartPK==1) {
2597 double *setPKinit = init_store->astore[0];
2598 setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2];
2599 setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4];
2600 setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5];
2601 }else {
2602 setPK(1,1)=cauchy(1,1); setPK(1,2)=cauchy(1,2); setPK(1,3)=cauchy(1,3);
2603 setPK(2,1)=cauchy(2,1); setPK(2,2)=cauchy(2,2); setPK(2,3)=cauchy(2,3);
2604 setPK(3,1)=cauchy(3,1); setPK(3,2)=cauchy(3,2); setPK(3,3)=cauchy(3,3);
2605 }
2606 initPK=0;
2607 }
2608
2609 CauchyStat_Step(Fi,Fdot,cauchy,setcauchy,setPK,vol,CSvol0,deltat,alpha);
2610
2611 // use currentPK as new target:
2612 // p_target: 0 1 2 3 4 5
2613 // x y z yz xz xy
2614
2615 p_target[0]=-setPK(1,1);
2616 p_target[1]=-setPK(2,2);
2617 p_target[2]=-setPK(3,3);
2618 if (pstyle == TRICLINIC) {
2619 p_target[3]=-setPK(2,3);
2620 p_target[4]=-setPK(1,3);
2621 p_target[5]=-setPK(1,2);
2622 }
2623
2624 p_hydro = 0.0;
2625 for (int i = 0; i < 3; i++)
2626 if (p_flag[i]) {
2627 p_hydro += p_target[i];
2628 }
2629 p_hydro /= pdim;
2630
2631 // save information for Cauchystat restart
2632 double *setPKinit = init_store->astore[0];
2633 setPKinit[0] = setcauchy(1,1);
2634 setPKinit[1] = setcauchy(1,2);
2635 setPKinit[2] = setcauchy(1,3);
2636 setPKinit[3] = setcauchy(2,2);
2637 setPKinit[4] = setcauchy(2,3);
2638 setPKinit[5] = setcauchy(3,3);
2639
2640 #undef H
2641 #undef Hdot
2642 #undef F
2643 #undef Fi
2644 #undef Fdot
2645 #undef invH0
2646 #undef cauchy
2647 #undef setcauchy
2648 #undef setPK
2649 #undef sigmahat
2650
2651 }
2652
2653 /* ----------------------------------------------------------------------
2654 CauchyStat_Step
2655
2656 Inputs:
2657 Fi(3,3) : inverse of the deformation gradient
2658 Fdot(3,3) : time derivative of the deformation gradient (velocity)
2659 cauchy(3,3) : current cauchy stress tensor
2660 setcauchy(3,3) : target cauchy stress tensor
2661 alpha : parameter =0.01
2662 setPK(3,3) : current PK stress tensor, at initialization (time=0) setPK=cauchy
2663 volume : current volume of the periodic box
2664 volume0 : initial volume of the periodic box
2665 deltat : simulation step increment
2666 alpha : parameter
2667
2668 Outputs:
2669 setPK(3,3) : PK stress tensor at the next time step
2670 ------------------------------------------------------------------------- */
2671
CauchyStat_Step(double (& Fi)[3][3],double (& Fdot)[3][3],double (& cauchy)[3][3],double (& setcauchy)[3][3],double (& setPK)[3][3],double volume,double volume0,double deltat,double alpha)2672 void FixNPTCauchy::CauchyStat_Step(double (&Fi)[3][3], double (&Fdot)[3][3],
2673 double (&cauchy)[3][3], double (&setcauchy)[3][3],
2674 double (&setPK)[3][3], double volume,
2675 double volume0, double deltat, double alpha)
2676 {
2677
2678 //macros to go from c to fortran style for arrays:
2679 #define F(row,col) (F[(row-1)][(col-1)])
2680 #define Fi(row,col) (Fi[(row-1)][(col-1)])
2681 #define Fdot(row,col) (Fdot[(row-1)][(col-1)])
2682 #define cauchy(row,col) (cauchy[(row-1)][(col-1)])
2683 #define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
2684 #define setPK(row,col) (setPK[(row-1)][(col-1)])
2685 #define uv(row,col) (uv[(row-1)][(col-1)])
2686 #define deltastress(row) (deltastress[(row-1)])
2687 #define fdotvec(row) (fdotvec[(row-1)])
2688 #define dsdf(row,col) (dsdf[(row-1)][(col-1)])
2689 #define dsds(row,col) (dsds[(row-1)][(col-1)])
2690 #define deltaF(row) (deltaF[(row-1)])
2691 #define deltaPK(row) (deltaPK[(row-1)])
2692
2693 //initialize local variables:
2694 int i,j,m,n;
2695
2696 double deltastress[6],fdotvec[6];
2697 double dsdf[6][6];
2698 double dsds[6][6];
2699 double jac;
2700 double deltaF[6];
2701 double deltaPK[6];
2702
2703 // zero arrays
2704 for (i = 0; i < 6; ++i) {
2705 deltaF[i] = deltaPK[i] = deltastress[i] = fdotvec[i] = 0.0;
2706 for (j = 0; j < 6; ++j)
2707 dsdf[i][j] = dsds[i][j] = 0.0;
2708 }
2709
2710 int uv[6][2];
2711 uv(1,1)=1; uv(1,2)=1;
2712 uv(2,1)=2; uv(2,2)=2;
2713 uv(3,1)=3; uv(3,2)=3;
2714 uv(4,1)=2; uv(4,2)=3;
2715 uv(5,1)=1; uv(5,2)=3;
2716 uv(6,1)=1; uv(6,2)=2;
2717
2718 for (int ii = 1;ii <= 6;ii++) {
2719 i=uv(ii,1);
2720 j=uv(ii,2);
2721 deltastress(ii)=setcauchy(i,j)-cauchy(i,j);
2722 if (ii>3) deltastress(ii)=deltastress(ii)*2.0;
2723 fdotvec(ii)=Fdot(i,j)*deltat;
2724 }
2725
2726 for (int ii = 1;ii <= 6;ii++) {
2727 i=uv(ii,1);
2728 j=uv(ii,2);
2729 for (int jj = 1;jj <= 6;jj++) {
2730 m=uv(jj,1);
2731 n=uv(jj,2);
2732 dsds(ii,jj) = Fi(i,m)*Fi(j,n) + Fi(i,n)*Fi(j,m) + Fi(j,m)*Fi(i,n) + Fi(j,n)*Fi(i,m);
2733 for (int l = 1;l <= 3;l++) {
2734 for (int k = 1;k <= 3;k++) {
2735 dsdf(ii,jj) = dsdf(ii,jj) + cauchy(k,l)*
2736 ( Fi(i,k)*Fi(j,l)*Fi(n,m) - Fi(i,m)*Fi(j,l)*Fi(n,k) - Fi(i,k)*Fi(j,m)*Fi(n,l) );
2737 }
2738 }
2739 }
2740 }
2741
2742 jac=volume/volume0;
2743 for (int ii = 1;ii <= 6;ii++) {
2744 for (int jj = 1;jj <= 6;jj++) {
2745 dsds(ii,jj)=dsds(ii,jj)*jac/4.0;
2746 dsdf(ii,jj)=dsdf(ii,jj)*jac;
2747 }
2748 }
2749
2750 for (int ii = 1;ii <= 6;ii++) {
2751 for (int jj = 1;jj <= 6;jj++) {
2752 deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj);
2753 }
2754 }
2755
2756 for (int ii = 1;ii <= 6;ii++) {
2757 for (int jj = 1;jj <= 6;jj++) {
2758 deltaPK(ii)=deltaPK(ii)+alpha*dsds(ii,jj)*deltastress(jj);
2759 }
2760 deltaPK(ii)=deltaPK(ii)+alpha*deltaF(ii);
2761 }
2762
2763 setPK(1,1)=setPK(1,1)+deltaPK(1);
2764 setPK(2,2)=setPK(2,2)+deltaPK(2);
2765 setPK(3,3)=setPK(3,3)+deltaPK(3);
2766 setPK(2,3)=setPK(2,3)+deltaPK(4);
2767 setPK(3,2)=setPK(3,2)+deltaPK(4);
2768 setPK(1,3)=setPK(1,3)+deltaPK(5);
2769 setPK(3,1)=setPK(3,1)+deltaPK(5);
2770 setPK(1,2)=setPK(1,2)+deltaPK(6);
2771 setPK(2,1)=setPK(2,1)+deltaPK(6);
2772
2773 #undef F
2774 #undef Fi
2775 #undef Fdot
2776 #undef cauchy
2777 #undef setcauchy
2778 #undef setPK
2779 #undef uv
2780 #undef deltastress
2781 #undef fdotvec
2782 #undef dsdf
2783 #undef dsds
2784 #undef deltaF
2785 #undef deltaPK
2786
2787 }
2788