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: Mark Stevens (SNL), Aidan Thompson (SNL), Zheng Gong (ENS Lyon)
17 ---------------------------------------------------------------------------------------- */
18
19 #include "fix_tgnh_drude.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_drude.h"
28 #include "force.h"
29 #include "irregular.h"
30 #include "kspace.h"
31 #include "memory.h"
32 #include "modify.h"
33 #include "neighbor.h"
34 #include "respa.h"
35 #include "update.h"
36
37 #include <cstring>
38 #include <cmath>
39
40 using namespace LAMMPS_NS;
41 using namespace FixConst;
42
43 #define DELTAFLIP 0.1
44 #define TILTMAX 1.5
45
46 enum{NOBIAS,BIAS};
47 enum{NONE,XYZ,XY,YZ,XZ};
48 enum{ISO,ANISO,TRICLINIC};
49
50 /* ----------------------------------------------------------------------
51 NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
52 ---------------------------------------------------------------------- */
53
FixTGNHDrude(LAMMPS * lmp,int narg,char ** arg)54 FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) :
55 Fix(lmp, narg, arg),
56 rfix(nullptr), irregular(nullptr), id_temp(nullptr), id_press(nullptr),
57 etamol(nullptr), etamol_dot(nullptr), etamol_dotdot(nullptr), etamol_mass(nullptr),
58 etaint(nullptr), etaint_dot(nullptr), etaint_dotdot(nullptr), etaint_mass(nullptr),
59 etadrude(nullptr), etadrude_dot(nullptr), etadrude_dotdot(nullptr), etadrude_mass(nullptr),
60 etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr)
61 {
62 if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");
63
64 restart_global = 1;
65 dynamic_group_allow = 0;
66 time_integrate = 1;
67 scalar_flag = 1;
68 vector_flag = 1;
69 global_freq = 1;
70 extscalar = 1;
71 extvector = 0;
72 ecouple_flag = 1;
73
74 // default values
75
76 pcouple = NONE;
77 mtchain = mpchain = 3;
78 nc_tchain = nc_pchain = 1;
79 mtk_flag = 1;
80 deviatoric_flag = 0;
81 nreset_h0 = 0;
82 flipflag = 1;
83
84 tcomputeflag = 0;
85 pcomputeflag = 0;
86 id_temp = nullptr;
87 id_press = nullptr;
88
89 // turn on tilt factor scaling, whenever applicable
90
91 dimension = domain->dimension;
92
93 scaleyz = scalexz = scalexy = 0;
94 if (domain->yperiodic && domain->xy != 0.0) scalexy = 1;
95 if (domain->zperiodic && dimension == 3) {
96 if (domain->yz != 0.0) scaleyz = 1;
97 if (domain->xz != 0.0) scalexz = 1;
98 }
99
100 // set fixed-point to default = center of cell
101
102 fixedpoint[0] = 0.5*(domain->boxlo[0]+domain->boxhi[0]);
103 fixedpoint[1] = 0.5*(domain->boxlo[1]+domain->boxhi[1]);
104 fixedpoint[2] = 0.5*(domain->boxlo[2]+domain->boxhi[2]);
105
106 tstat_flag = 0;
107 double t_period = 0.0, tdrude_period = 0.0;
108
109 double p_period[6];
110 for (int i = 0; i < 6; i++) {
111 p_start[i] = p_stop[i] = p_period[i] = p_target[i] = 0.0;
112 p_flag[i] = 0;
113 }
114
115 // process keywords
116
117 int iarg = 3;
118
119 while (iarg < narg) {
120 if (strcmp(arg[iarg],"temp") == 0) {
121 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
122 tstat_flag = 1;
123 t_start = utils::numeric(FLERR,arg[iarg+1],false,lmp);
124 t_target = t_start;
125 t_stop = utils::numeric(FLERR,arg[iarg+2],false,lmp);
126 t_period = utils::numeric(FLERR,arg[iarg+3],false,lmp);
127 if (t_start <= 0.0 || t_stop <= 0.0)
128 error->all(FLERR,
129 "Target temperature for fix nvt/npt/nph cannot be 0.0");
130 tdrude_target = utils::numeric(FLERR,arg[iarg+4],false,lmp);
131 tdrude_period = utils::numeric(FLERR,arg[iarg+5],false,lmp);
132 iarg += 6;
133
134 } else if (strcmp(arg[iarg],"iso") == 0) {
135 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
136 pcouple = XYZ;
137 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
138 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
139 p_period[0] = p_period[1] = p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
140 p_flag[0] = p_flag[1] = p_flag[2] = 1;
141 if (dimension == 2) {
142 p_start[2] = p_stop[2] = p_period[2] = 0.0;
143 p_flag[2] = 0;
144 }
145 iarg += 4;
146 } else if (strcmp(arg[iarg],"aniso") == 0) {
147 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
148 pcouple = NONE;
149 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
150 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
151 p_period[0] = p_period[1] = p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
152 p_flag[0] = p_flag[1] = p_flag[2] = 1;
153 if (dimension == 2) {
154 p_start[2] = p_stop[2] = p_period[2] = 0.0;
155 p_flag[2] = 0;
156 }
157 iarg += 4;
158 } else if (strcmp(arg[iarg],"tri") == 0) {
159 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
160 pcouple = NONE;
161 scalexy = scalexz = scaleyz = 0;
162 p_start[0] = p_start[1] = p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
163 p_stop[0] = p_stop[1] = p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
164 p_period[0] = p_period[1] = p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
165 p_flag[0] = p_flag[1] = p_flag[2] = 1;
166 p_start[3] = p_start[4] = p_start[5] = 0.0;
167 p_stop[3] = p_stop[4] = p_stop[5] = 0.0;
168 p_period[3] = p_period[4] = p_period[5] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
169 p_flag[3] = p_flag[4] = p_flag[5] = 1;
170 if (dimension == 2) {
171 p_start[2] = p_stop[2] = p_period[2] = 0.0;
172 p_flag[2] = 0;
173 p_start[3] = p_stop[3] = p_period[3] = 0.0;
174 p_flag[3] = 0;
175 p_start[4] = p_stop[4] = p_period[4] = 0.0;
176 p_flag[4] = 0;
177 }
178 iarg += 4;
179 } else if (strcmp(arg[iarg],"x") == 0) {
180 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
181 p_start[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
182 p_stop[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
183 p_period[0] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
184 p_flag[0] = 1;
185 deviatoric_flag = 1;
186 iarg += 4;
187 } else if (strcmp(arg[iarg],"y") == 0) {
188 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
189 p_start[1] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
190 p_stop[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
191 p_period[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
192 p_flag[1] = 1;
193 deviatoric_flag = 1;
194 iarg += 4;
195 } else if (strcmp(arg[iarg],"z") == 0) {
196 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
197 p_start[2] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
198 p_stop[2] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
199 p_period[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
200 p_flag[2] = 1;
201 deviatoric_flag = 1;
202 iarg += 4;
203 if (dimension == 2)
204 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
205
206 } else if (strcmp(arg[iarg],"yz") == 0) {
207 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
208 p_start[3] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
209 p_stop[3] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
210 p_period[3] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
211 p_flag[3] = 1;
212 deviatoric_flag = 1;
213 scaleyz = 0;
214 iarg += 4;
215 if (dimension == 2)
216 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
217 } else if (strcmp(arg[iarg],"xz") == 0) {
218 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
219 p_start[4] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
220 p_stop[4] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
221 p_period[4] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
222 p_flag[4] = 1;
223 deviatoric_flag = 1;
224 scalexz = 0;
225 iarg += 4;
226 if (dimension == 2)
227 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
228 } else if (strcmp(arg[iarg],"xy") == 0) {
229 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
230 p_start[5] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
231 p_stop[5] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
232 p_period[5] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
233 p_flag[5] = 1;
234 deviatoric_flag = 1;
235 scalexy = 0;
236 iarg += 4;
237
238 } else if (strcmp(arg[iarg],"couple") == 0) {
239 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
240 if (strcmp(arg[iarg+1],"xyz") == 0) pcouple = XYZ;
241 else if (strcmp(arg[iarg+1],"xy") == 0) pcouple = XY;
242 else if (strcmp(arg[iarg+1],"yz") == 0) pcouple = YZ;
243 else if (strcmp(arg[iarg+1],"xz") == 0) pcouple = XZ;
244 else if (strcmp(arg[iarg+1],"none") == 0) pcouple = NONE;
245 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
246 iarg += 2;
247 } else if (strcmp(arg[iarg],"tchain") == 0) {
248 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
249 mtchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
250 if (mtchain < 1) error->all(FLERR,"Illegal fix nvt/npt/nph command");
251 iarg += 2;
252 } else if (strcmp(arg[iarg],"pchain") == 0) {
253 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
254 mpchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
255 if (mpchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
256 iarg += 2;
257 } else if (strcmp(arg[iarg],"mtk") == 0) {
258 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
259 if (strcmp(arg[iarg+1],"yes") == 0) mtk_flag = 1;
260 else if (strcmp(arg[iarg+1],"no") == 0) mtk_flag = 0;
261 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
262 iarg += 2;
263 } else if (strcmp(arg[iarg],"tloop") == 0) {
264 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
265 nc_tchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
266 if (nc_tchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
267 iarg += 2;
268 } else if (strcmp(arg[iarg],"ploop") == 0) {
269 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
270 nc_pchain = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
271 if (nc_pchain < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
272 iarg += 2;
273 } else if (strcmp(arg[iarg],"nreset") == 0) {
274 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
275 nreset_h0 = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
276 if (nreset_h0 < 0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
277 iarg += 2;
278 } else if (strcmp(arg[iarg],"scalexy") == 0) {
279 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
280 if (strcmp(arg[iarg+1],"yes") == 0) scalexy = 1;
281 else if (strcmp(arg[iarg+1],"no") == 0) scalexy = 0;
282 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
283 iarg += 2;
284 } else if (strcmp(arg[iarg],"scalexz") == 0) {
285 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
286 if (strcmp(arg[iarg+1],"yes") == 0) scalexz = 1;
287 else if (strcmp(arg[iarg+1],"no") == 0) scalexz = 0;
288 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
289 iarg += 2;
290 } else if (strcmp(arg[iarg],"scaleyz") == 0) {
291 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
292 if (strcmp(arg[iarg+1],"yes") == 0) scaleyz = 1;
293 else if (strcmp(arg[iarg+1],"no") == 0) scaleyz = 0;
294 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
295 iarg += 2;
296 } else if (strcmp(arg[iarg],"flip") == 0) {
297 if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
298 if (strcmp(arg[iarg+1],"yes") == 0) flipflag = 1;
299 else if (strcmp(arg[iarg+1],"no") == 0) flipflag = 0;
300 else error->all(FLERR,"Illegal fix nvt/npt/nph command");
301 iarg += 2;
302 } else if (strcmp(arg[iarg],"fixedpoint") == 0) {
303 if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
304 fixedpoint[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
305 fixedpoint[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
306 fixedpoint[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
307 iarg += 4;
308 } else error->all(FLERR,"Illegal fix nvt/npt/nph command");
309 }
310
311 // error checks
312
313 if (dimension == 2 && (p_flag[2] || p_flag[3] || p_flag[4]))
314 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
315 if (dimension == 2 && (pcouple == YZ || pcouple == XZ))
316 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
317 if (dimension == 2 && (scalexz == 1 || scaleyz == 1 ))
318 error->all(FLERR,"Invalid fix nvt/npt/nph command for a 2d simulation");
319
320 if (pcouple == XYZ && (p_flag[0] == 0 || p_flag[1] == 0))
321 error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
322 if (pcouple == XYZ && dimension == 3 && p_flag[2] == 0)
323 error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
324 if (pcouple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
325 error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
326 if (pcouple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
327 error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
328 if (pcouple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
329 error->all(FLERR,"Invalid fix nvt/npt/nph command pressure settings");
330
331 // require periodicity in tensile dimension
332
333 if (p_flag[0] && domain->xperiodic == 0)
334 error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
335 if (p_flag[1] && domain->yperiodic == 0)
336 error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
337 if (p_flag[2] && domain->zperiodic == 0)
338 error->all(FLERR,"Cannot use fix nvt/npt/nph on a non-periodic dimension");
339
340 // require periodicity in 2nd dim of off-diagonal tilt component
341
342 if (p_flag[3] && domain->zperiodic == 0)
343 error->all(FLERR,
344 "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
345 if (p_flag[4] && domain->zperiodic == 0)
346 error->all(FLERR,
347 "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
348 if (p_flag[5] && domain->yperiodic == 0)
349 error->all(FLERR,
350 "Cannot use fix nvt/npt/nph on a 2nd non-periodic dimension");
351
352 if (scaleyz == 1 && domain->zperiodic == 0)
353 error->all(FLERR,"Cannot use fix nvt/npt/nph "
354 "with yz scaling when z is non-periodic dimension");
355 if (scalexz == 1 && domain->zperiodic == 0)
356 error->all(FLERR,"Cannot use fix nvt/npt/nph "
357 "with xz scaling when z is non-periodic dimension");
358 if (scalexy == 1 && domain->yperiodic == 0)
359 error->all(FLERR,"Cannot use fix nvt/npt/nph "
360 "with xy scaling when y is non-periodic dimension");
361
362 if (p_flag[3] && scaleyz == 1)
363 error->all(FLERR,"Cannot use fix nvt/npt/nph with "
364 "both yz dynamics and yz scaling");
365 if (p_flag[4] && scalexz == 1)
366 error->all(FLERR,"Cannot use fix nvt/npt/nph with "
367 "both xz dynamics and xz scaling");
368 if (p_flag[5] && scalexy == 1)
369 error->all(FLERR,"Cannot use fix nvt/npt/nph with "
370 "both xy dynamics and xy scaling");
371
372 if (!domain->triclinic && (p_flag[3] || p_flag[4] || p_flag[5]))
373 error->all(FLERR,"Can not specify Pxy/Pxz/Pyz in "
374 "fix nvt/npt/nph with non-triclinic box");
375
376 if (pcouple == XYZ && dimension == 3 &&
377 (p_start[0] != p_start[1] || p_start[0] != p_start[2] ||
378 p_stop[0] != p_stop[1] || p_stop[0] != p_stop[2] ||
379 p_period[0] != p_period[1] || p_period[0] != p_period[2]))
380 error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
381 if (pcouple == XYZ && dimension == 2 &&
382 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
383 p_period[0] != p_period[1]))
384 error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
385 if (pcouple == XY &&
386 (p_start[0] != p_start[1] || p_stop[0] != p_stop[1] ||
387 p_period[0] != p_period[1]))
388 error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
389 if (pcouple == YZ &&
390 (p_start[1] != p_start[2] || p_stop[1] != p_stop[2] ||
391 p_period[1] != p_period[2]))
392 error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
393 if (pcouple == XZ &&
394 (p_start[0] != p_start[2] || p_stop[0] != p_stop[2] ||
395 p_period[0] != p_period[2]))
396 error->all(FLERR,"Invalid fix nvt/npt/nph pressure settings");
397
398 if ((tstat_flag && t_period <= 0.0) ||
399 (p_flag[0] && p_period[0] <= 0.0) ||
400 (p_flag[1] && p_period[1] <= 0.0) ||
401 (p_flag[2] && p_period[2] <= 0.0) ||
402 (p_flag[3] && p_period[3] <= 0.0) ||
403 (p_flag[4] && p_period[4] <= 0.0) ||
404 (p_flag[5] && p_period[5] <= 0.0))
405 error->all(FLERR,"Fix nvt/npt/nph damping parameters must be > 0.0");
406
407 // set pstat_flag and box change and restart_pbc variables
408
409 pre_exchange_flag = 0;
410 pstat_flag = 0;
411 pstyle = ISO;
412
413 for (int i = 0; i < 6; i++)
414 if (p_flag[i]) pstat_flag = 1;
415
416 if (pstat_flag) {
417 if (p_flag[0]) box_change |= BOX_CHANGE_X;
418 if (p_flag[1]) box_change |= BOX_CHANGE_Y;
419 if (p_flag[2]) box_change |= BOX_CHANGE_Z;
420 if (p_flag[3]) box_change |= BOX_CHANGE_YZ;
421 if (p_flag[4]) box_change |= BOX_CHANGE_XZ;
422 if (p_flag[5]) box_change |= BOX_CHANGE_XY;
423 no_change_box = 1;
424
425 // pstyle = TRICLINIC if any off-diagonal term is controlled -> 6 dof
426 // else pstyle = ISO if XYZ coupling or XY coupling in 2d -> 1 dof
427 // else pstyle = ANISO -> 3 dof
428
429 if (p_flag[3] || p_flag[4] || p_flag[5]) pstyle = TRICLINIC;
430 else if (pcouple == XYZ || (dimension == 2 && pcouple == XY)) pstyle = ISO;
431 else pstyle = ANISO;
432
433 // pre_exchange only required if flips can occur due to shape changes
434
435 if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
436 pre_exchange_flag = pre_exchange_migrate = 1;
437 if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0))
438 pre_exchange_flag = pre_exchange_migrate = 1;
439 }
440
441 // convert input periods to frequencies
442
443 t_freq = tdrude_freq = 0.0;
444 p_freq[0] = p_freq[1] = p_freq[2] = p_freq[3] = p_freq[4] = p_freq[5] = 0.0;
445
446 if (tstat_flag) {
447 t_freq = 1.0 / t_period;
448 tdrude_freq = 1.0 / tdrude_period;
449 }
450 if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
451 if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
452 if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
453 if (p_flag[3]) p_freq[3] = 1.0 / p_period[3];
454 if (p_flag[4]) p_freq[4] = 1.0 / p_period[4];
455 if (p_flag[5]) p_freq[5] = 1.0 / p_period[5];
456
457 // Nose/Hoover temp and pressure init
458
459 size_vector = 3;
460
461 if (tstat_flag) {
462 int ich;
463
464 etaint = new double[mtchain];
465 // add one extra dummy thermostat for eta_dot, set to zero
466 etaint_dot = new double[mtchain+1];
467 etaint_dot[mtchain] = 0.0;
468 etaint_dotdot = new double[mtchain];
469 for (ich = 0; ich < mtchain; ich++) {
470 etaint[ich] = etaint_dot[ich] = etaint_dotdot[ich] = 0.0;
471 }
472 etaint_mass = new double[mtchain];
473
474 etamol = new double[mtchain];
475 // add one extra dummy thermostat for eta_dot, set to zero
476 etamol_dot = new double[mtchain+1];
477 etamol_dot[mtchain] = 0.0;
478 etamol_dotdot = new double[mtchain];
479 for (ich = 0; ich < mtchain; ich++) {
480 etamol[ich] = etamol_dot[ich] = etamol_dotdot[ich] = 0.0;
481 }
482 etamol_mass = new double[mtchain];
483
484 etadrude = new double[mtchain];
485 // add one extra dummy thermostat for eta_dot, set to zero
486 etadrude_dot = new double[mtchain+1];
487 etadrude_dot[mtchain] = 0.0;
488 etadrude_dotdot = new double[mtchain];
489 for (ich = 0; ich < mtchain; ich++) {
490 etadrude[ich] = etadrude_dot[ich] = etadrude_dotdot[ich] = 0.0;
491 }
492 etadrude_mass = new double[mtchain];
493 }
494
495 if (pstat_flag) {
496 omega[0] = omega[1] = omega[2] = 0.0;
497 omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
498 omega_mass[0] = omega_mass[1] = omega_mass[2] = 0.0;
499 omega[3] = omega[4] = omega[5] = 0.0;
500 omega_dot[3] = omega_dot[4] = omega_dot[5] = 0.0;
501 omega_mass[3] = omega_mass[4] = omega_mass[5] = 0.0;
502
503 if (mpchain) {
504 int ich;
505 etap = new double[mpchain];
506
507 // add one extra dummy thermostat, set to zero
508
509 etap_dot = new double[mpchain+1];
510 etap_dot[mpchain] = 0.0;
511 etap_dotdot = new double[mpchain];
512 for (ich = 0; ich < mpchain; ich++) {
513 etap[ich] = etap_dot[ich] =
514 etap_dotdot[ich] = 0.0;
515 }
516 etap_mass = new double[mpchain];
517 }
518 }
519
520 nrigid = 0;
521 rfix = nullptr;
522
523 if (pre_exchange_flag) irregular = new Irregular(lmp);
524 else irregular = nullptr;
525
526 // initialize vol0,t0 to zero to signal uninitialized
527 // values then assigned in init(), if necessary
528
529 vol0 = t0 = 0.0;
530
531 // find fix drude
532 int ifix;
533 for (ifix = 0; ifix < modify->nfix; ifix++)
534 if (strcmp(modify->fix[ifix]->style,"drude") == 0) break;
535 if (ifix == modify->nfix) error->all(FLERR, "fix tgnh/drude requires fix drude");
536 fix_drude = (FixDrude *) modify->fix[ifix];
537
538 // make sure ghost atoms have velocity
539 if (!comm->ghost_velocity)
540 error->all(FLERR,"fix tgnh/drude requires ghost velocities. Use comm_modify vel yes");
541 }
542
543 /* ---------------------------------------------------------------------- */
544
~FixTGNHDrude()545 FixTGNHDrude::~FixTGNHDrude()
546 {
547 if (copymode) return;
548
549 delete [] rfix;
550
551 delete irregular;
552
553 // delete temperature and pressure if fix created them
554
555 if (tcomputeflag) modify->delete_compute(id_temp);
556 delete [] id_temp;
557
558 if (tstat_flag) {
559 delete [] etaint;
560 delete [] etaint_dot;
561 delete [] etaint_dotdot;
562 delete [] etaint_mass;
563 delete [] etamol;
564 delete [] etamol_dot;
565 delete [] etamol_dotdot;
566 delete [] etamol_mass;
567 delete [] etadrude;
568 delete [] etadrude_dot;
569 delete [] etadrude_dotdot;
570 delete [] etadrude_mass;
571 }
572
573 if (pstat_flag) {
574 if (pcomputeflag) modify->delete_compute(id_press);
575 delete [] id_press;
576 if (mpchain) {
577 delete [] etap;
578 delete [] etap_dot;
579 delete [] etap_dotdot;
580 delete [] etap_mass;
581 }
582 }
583 }
584
585 /* ---------------------------------------------------------------------- */
586
setmask()587 int FixTGNHDrude::setmask()
588 {
589 int mask = 0;
590 mask |= INITIAL_INTEGRATE;
591 mask |= FINAL_INTEGRATE;
592 mask |= INITIAL_INTEGRATE_RESPA;
593 mask |= PRE_FORCE_RESPA;
594 mask |= FINAL_INTEGRATE_RESPA;
595 if (pre_exchange_flag) mask |= PRE_EXCHANGE;
596 return mask;
597 }
598
599 /* ---------------------------------------------------------------------- */
600
init()601 void FixTGNHDrude::init()
602 {
603 // ensure no conflict with fix deform
604
605 if (pstat_flag)
606 for (int i = 0; i < modify->nfix; i++)
607 if (strcmp(modify->fix[i]->style,"deform") == 0) {
608 int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
609 if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
610 (p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) ||
611 (p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5]))
612 error->all(FLERR,"Cannot use fix npt and fix deform on "
613 "same component of stress tensor");
614 }
615
616 // set temperature and pressure ptrs
617
618 int icompute = modify->find_compute(id_temp);
619 if (icompute < 0)
620 error->all(FLERR,"Temperature ID for fix nvt/npt does not exist");
621 temperature = modify->compute[icompute];
622
623 if (temperature->tempbias) which = BIAS;
624 else which = NOBIAS;
625
626 if (pstat_flag) {
627 icompute = modify->find_compute(id_press);
628 if (icompute < 0)
629 error->all(FLERR,"Pressure ID for fix npt/nph does not exist");
630 pressure = modify->compute[icompute];
631 }
632
633 // set timesteps and frequencies
634
635 dtv = update->dt;
636 dtf = 0.5 * update->dt * force->ftm2v;
637 dthalf = 0.5 * update->dt;
638 dt4 = 0.25 * update->dt;
639 dt8 = 0.125 * update->dt;
640 dto = dthalf;
641
642 p_freq_max = 0.0;
643 if (pstat_flag) {
644 p_freq_max = MAX(p_freq[0],p_freq[1]);
645 p_freq_max = MAX(p_freq_max,p_freq[2]);
646 if (pstyle == TRICLINIC) {
647 p_freq_max = MAX(p_freq_max,p_freq[3]);
648 p_freq_max = MAX(p_freq_max,p_freq[4]);
649 p_freq_max = MAX(p_freq_max,p_freq[5]);
650 }
651 }
652
653 // tally the number of dimensions that are barostatted
654 // set initial volume and reference cell, if not already done
655
656 if (pstat_flag) {
657 pdim = p_flag[0] + p_flag[1] + p_flag[2];
658 if (vol0 == 0.0) {
659 if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
660 else vol0 = domain->xprd * domain->yprd;
661 h0_inv[0] = domain->h_inv[0];
662 h0_inv[1] = domain->h_inv[1];
663 h0_inv[2] = domain->h_inv[2];
664 h0_inv[3] = domain->h_inv[3];
665 h0_inv[4] = domain->h_inv[4];
666 h0_inv[5] = domain->h_inv[5];
667 }
668 }
669
670 boltz = force->boltz;
671 nktv2p = force->nktv2p;
672
673 if (force->kspace) kspace_flag = 1;
674 else kspace_flag = 0;
675
676 if (utils::strmatch(update->integrate_style,"^respa")) {
677 nlevels_respa = ((Respa *) update->integrate)->nlevels;
678 step_respa = ((Respa *) update->integrate)->step;
679 dto = 0.5*step_respa[0];
680 }
681
682 // detect if any rigid fixes exist so rigid bodies move when box is remapped
683 // rfix[] = indices to each fix rigid
684
685 delete [] rfix;
686 nrigid = 0;
687 rfix = nullptr;
688
689 for (int i = 0; i < modify->nfix; i++)
690 if (modify->fix[i]->rigid_flag) nrigid++;
691 if (nrigid) {
692 rfix = new int[nrigid];
693 nrigid = 0;
694 for (int i = 0; i < modify->nfix; i++)
695 if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
696 }
697 }
698
699 /* ----------------------------------------------------------------------
700 compute T,P before integrator starts
701 ------------------------------------------------------------------------- */
702
setup_mol_mass_dof()703 void FixTGNHDrude::setup_mol_mass_dof() {
704 double *mass = atom->mass;
705 int *mask = atom->mask;
706 tagint *molecule = atom->molecule;
707 int *type = atom->type;
708 int *drudetype = fix_drude->drudetype;
709 int n_drude, n_drude_tmp = 0;
710 tagint id_mol = 0, n_mol_in_group = 0;
711
712 for (int i = 0; i < atom->nlocal; i++) {
713 // molecule id starts from 1. max(id_mol) equals to the number of molecules in the system
714 id_mol = std::max(id_mol, molecule[i]);
715 if (mask[i] & groupbit) {
716 if (drudetype[type[i]] == DRUDE_TYPE)
717 n_drude_tmp++;
718 }
719 }
720 MPI_Allreduce(&n_drude_tmp, &n_drude, 1, MPI_LMP_TAGINT, MPI_SUM, world);
721 MPI_Allreduce(&id_mol, &n_mol, 1, MPI_LMP_TAGINT, MPI_MAX, world);
722
723 // use flag_mol to determine the number of molecules in the fix group
724 int *flag_mol = new int[n_mol + 1];
725 int *flag_mol_tmp = new int[n_mol + 1];
726 memset(flag_mol_tmp, 0, sizeof(int) * (n_mol + 1));
727 for (int i = 0; i < atom->nlocal; i++) {
728 if (mask[i] & groupbit) {
729 flag_mol_tmp[molecule[i]] = 1;
730 }
731 }
732 MPI_Allreduce(flag_mol_tmp, flag_mol, n_mol + 1, MPI_INT, MPI_SUM, world);
733 for (int i = 1; i < n_mol + 1; i++) {
734 if (flag_mol[i])
735 n_mol_in_group++;
736 }
737 delete[] flag_mol;
738 delete[] flag_mol_tmp;
739
740 // length of v_mol set to n_mol+1, so that the subscript start from 1, we can call v_mol[n_mol]
741 memory->create(v_mol, n_mol + 1, 3, "fix_tgnh_drude::v_mol");
742 memory->create(v_mol_tmp, n_mol + 1, 3, "fix_tgnh_drude::v_mol_tmp");
743 memory->create(mass_mol, n_mol + 1, "fix_tgnh_drude::mass_mol");
744
745 double *mass_tmp = new double[n_mol + 1];
746 memset(mass_tmp, 0, sizeof(double) * (n_mol + 1));
747 for (int i = 0; i < atom->nlocal; i++) {
748 id_mol = molecule[i];
749 mass_tmp[id_mol] += mass[type[i]];
750 }
751 MPI_Allreduce(mass_tmp, mass_mol, n_mol + 1, MPI_DOUBLE, MPI_SUM, world);
752 delete[] mass_tmp;
753
754 // DOFs
755 t_current = temperature->compute_scalar();
756 tdof = temperature->dof;
757 // remove DOFs of COM translational motion based on the number of molecules in the group
758 dof_mol = 3.0 * n_mol_in_group - 3.0 * n_mol_in_group / n_mol;
759 dof_drude = 3.0 * n_drude;
760 dof_int = tdof - dof_mol - dof_drude;
761
762 if (comm->me == 0) {
763 if (screen) {
764 fprintf(screen, "TGNHC thermostat for Drude model\n");
765 fprintf(screen, " DOFs of molecules, atoms and dipoles: %.1f %.1f %.1f\n",
766 dof_mol, dof_int, dof_drude);
767 }
768 if (logfile) {
769 fprintf(logfile, "TGNHC thermostat for Drude model\n");
770 fprintf(logfile, " DOFs of molecules, atoms and dipoles: %.1f %.1f %.1f\n",
771 dof_mol, dof_int, dof_drude);
772 }
773 }
774 if (dof_mol <=0 or dof_int <=0 or dof_drude <=0)
775 error->all(FLERR, "TGNHC thermostat requires DOFs of molecules, atoms and dipoles larger than 0");
776 }
777
setup(int)778 void FixTGNHDrude::setup(int /*vflag*/)
779 {
780 setup_mol_mass_dof();
781 // t_target is needed by NVT and NPT in compute_scalar()
782 // If no thermostat or using fix nphug,
783 // t_target must be defined by other means.
784
785 if (tstat_flag && strstr(style,"nphug") == nullptr) {
786 compute_temp_target();
787 } else if (pstat_flag) {
788
789 // t0 = reference temperature for masses
790 // cannot be done in init() b/c temperature cannot be called there
791 // is b/c Modify::init() inits computes after fixes due to dof dependence
792 // guesstimate a unit-dependent t0 if actual T = 0.0
793 // if it was read in from a restart file, leave it be
794
795 if (t0 == 0.0) {
796 t0 = temperature->compute_scalar();
797 if (t0 == 0.0) {
798 if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
799 else t0 = 300.0;
800 }
801 }
802 t_target = t0;
803 }
804
805 if (pstat_flag) compute_press_target();
806
807 if (pstat_flag) {
808 if (pstyle == ISO) pressure->compute_scalar();
809 else pressure->compute_vector();
810 couple();
811 pressure->addstep(update->ntimestep+1);
812 }
813
814 // masses and initial forces on thermostat variables
815
816 if (tstat_flag) {
817 etaint_mass[0] = ke2int_target / (t_freq * t_freq);
818 etamol_mass[0] = ke2mol_target / (t_freq * t_freq);
819 etadrude_mass[0] = ke2drude_target / (tdrude_freq * tdrude_freq);
820 for (int ich = 1; ich < mtchain; ich++) {
821 etaint_mass[ich] = boltz * t_target / (t_freq * t_freq);
822 etamol_mass[ich] = boltz * t_target / (t_freq * t_freq);
823 etadrude_mass[ich] = boltz * tdrude_target / (tdrude_freq * tdrude_freq);
824
825 etaint_dotdot[ich] = (etaint_mass[ich - 1] * etaint_dot[ich - 1] * etaint_dot[ich - 1] -
826 boltz * t_target) / etaint_mass[ich];
827 etamol_dotdot[ich] = (etamol_mass[ich - 1] * etamol_dot[ich - 1] * etamol_dot[ich - 1] -
828 boltz * t_target) / etamol_mass[ich];
829 etadrude_dotdot[ich] = (etadrude_mass[ich - 1] * etadrude_dot[ich - 1] * etadrude_dot[ich - 1] -
830 boltz * tdrude_target) / etadrude_mass[ich];
831 }
832 }
833
834 // masses and initial forces on barostat variables
835
836 if (pstat_flag) {
837 double kt = boltz * t_target;
838 double nkt = (atom->natoms + 1) * kt;
839
840 for (int i = 0; i < 3; i++)
841 if (p_flag[i])
842 omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
843
844 if (pstyle == TRICLINIC) {
845 for (int i = 3; i < 6; i++)
846 if (p_flag[i]) omega_mass[i] = nkt/(p_freq[i]*p_freq[i]);
847 }
848
849 // masses and initial forces on barostat thermostat variables
850
851 if (mpchain) {
852 etap_mass[0] = boltz * t_target / (p_freq_max*p_freq_max);
853 for (int ich = 1; ich < mpchain; ich++)
854 etap_mass[ich] = boltz * t_target / (p_freq_max*p_freq_max);
855 for (int ich = 1; ich < mpchain; ich++)
856 etap_dotdot[ich] =
857 (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] -
858 boltz * t_target) / etap_mass[ich];
859 }
860 }
861 }
862
863 /* ----------------------------------------------------------------------
864 1st half of Verlet update
865 ------------------------------------------------------------------------- */
866
initial_integrate(int)867 void FixTGNHDrude::initial_integrate(int /*vflag*/)
868 {
869 // update eta_press_dot
870
871 if (pstat_flag && mpchain) nhc_press_integrate();
872
873 // update eta_dot
874
875 if (tstat_flag) {
876 compute_temp_target();
877 nhc_temp_integrate();
878 }
879
880 // need to recompute pressure to account for change in KE
881 // t_current is up-to-date, but compute_temperature is not
882 // compute appropriately coupled elements of mvv_current
883
884 if (pstat_flag) {
885 if (pstyle == ISO) {
886 temperature->compute_scalar();
887 pressure->compute_scalar();
888 } else {
889 temperature->compute_vector();
890 pressure->compute_vector();
891 }
892 couple();
893 pressure->addstep(update->ntimestep+1);
894 }
895
896 if (pstat_flag) {
897 compute_press_target();
898 nh_omega_dot();
899 nh_v_press();
900 }
901
902 nve_v();
903
904 // remap simulation box by 1/2 step
905
906 if (pstat_flag) remap();
907
908 nve_x();
909
910 // remap simulation box by 1/2 step
911 // redo KSpace coeffs since volume has changed
912
913 if (pstat_flag) {
914 remap();
915 if (kspace_flag) force->kspace->setup();
916 }
917 }
918
919 /* ----------------------------------------------------------------------
920 2nd half of Verlet update
921 ------------------------------------------------------------------------- */
922
final_integrate()923 void FixTGNHDrude::final_integrate()
924 {
925 nve_v();
926
927 // re-compute temp before nh_v_press()
928 // only needed for temperature computes with BIAS on reneighboring steps:
929 // b/c some biases store per-atom values (e.g. temp/profile)
930 // per-atom values are invalid if reneigh/comm occurred
931 // since temp->compute() in initial_integrate()
932
933 if (which == BIAS && neighbor->ago == 0)
934 t_current = temperature->compute_scalar();
935
936 if (pstat_flag) nh_v_press();
937
938 // compute new T,P after velocities rescaled by nh_v_press()
939 // compute appropriately coupled elements of mvv_current
940
941 t_current = temperature->compute_scalar();
942 tdof = temperature->dof;
943
944 // need to recompute pressure to account for change in KE
945 // t_current is up-to-date, but compute_temperature is not
946 // compute appropriately coupled elements of mvv_current
947
948 if (pstat_flag) {
949 if (pstyle == ISO) pressure->compute_scalar();
950 else {
951 temperature->compute_vector();
952 pressure->compute_vector();
953 }
954 couple();
955 pressure->addstep(update->ntimestep+1);
956 }
957
958 if (pstat_flag) nh_omega_dot();
959
960 // update eta_dot
961 // update eta_press_dot
962
963 if (tstat_flag) nhc_temp_integrate();
964 if (pstat_flag && mpchain) nhc_press_integrate();
965 }
966
967 /* ---------------------------------------------------------------------- */
968
initial_integrate_respa(int,int ilevel,int)969 void FixTGNHDrude::initial_integrate_respa(int /*vflag*/, int ilevel, int /*iloop*/)
970 {
971 // set timesteps by level
972
973 dtv = step_respa[ilevel];
974 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
975 dthalf = 0.5 * step_respa[ilevel];
976
977 // outermost level - update eta_dot and omega_dot, apply to v
978 // all other levels - NVE update of v
979 // x,v updates only performed for atoms in group
980
981 if (ilevel == nlevels_respa-1) {
982
983 // update eta_press_dot
984
985 if (pstat_flag && mpchain) nhc_press_integrate();
986
987 // update eta_dot
988
989 if (tstat_flag) {
990 compute_temp_target();
991 nhc_temp_integrate();
992 }
993
994 // recompute pressure to account for change in KE
995 // t_current is up-to-date, but compute_temperature is not
996 // compute appropriately coupled elements of mvv_current
997
998 if (pstat_flag) {
999 if (pstyle == ISO) {
1000 temperature->compute_scalar();
1001 pressure->compute_scalar();
1002 } else {
1003 temperature->compute_vector();
1004 pressure->compute_vector();
1005 }
1006 couple();
1007 pressure->addstep(update->ntimestep+1);
1008 }
1009
1010 if (pstat_flag) {
1011 compute_press_target();
1012 nh_omega_dot();
1013 nh_v_press();
1014 }
1015
1016 nve_v();
1017
1018 } else nve_v();
1019
1020 // innermost level - also update x only for atoms in group
1021 // if barostat, perform 1/2 step remap before and after
1022
1023 if (ilevel == 0) {
1024 if (pstat_flag) remap();
1025 nve_x();
1026 if (pstat_flag) remap();
1027 }
1028 }
1029
1030 /* ---------------------------------------------------------------------- */
1031
pre_force_respa(int,int ilevel,int)1032 void FixTGNHDrude::pre_force_respa(int /*vflag*/, int ilevel, int /*iloop*/)
1033 {
1034 // if barostat, redo KSpace coeffs at outermost level,
1035 // since volume has changed
1036
1037 if (ilevel == nlevels_respa-1 && kspace_flag && pstat_flag)
1038 force->kspace->setup();
1039 }
1040
1041 /* ---------------------------------------------------------------------- */
1042
final_integrate_respa(int ilevel,int)1043 void FixTGNHDrude::final_integrate_respa(int ilevel, int /*iloop*/)
1044 {
1045 // set timesteps by level
1046
1047 dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
1048 dthalf = 0.5 * step_respa[ilevel];
1049
1050 // outermost level - update eta_dot and omega_dot, apply via final_integrate
1051 // all other levels - NVE update of v
1052
1053 if (ilevel == nlevels_respa-1) final_integrate();
1054 else nve_v();
1055 }
1056
1057 /* ---------------------------------------------------------------------- */
1058
couple()1059 void FixTGNHDrude::couple()
1060 {
1061 double *tensor = pressure->vector;
1062
1063 if (pstyle == ISO)
1064 p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
1065 else if (pcouple == XYZ) {
1066 double ave = 1.0/3.0 * (tensor[0] + tensor[1] + tensor[2]);
1067 p_current[0] = p_current[1] = p_current[2] = ave;
1068 } else if (pcouple == XY) {
1069 double ave = 0.5 * (tensor[0] + tensor[1]);
1070 p_current[0] = p_current[1] = ave;
1071 p_current[2] = tensor[2];
1072 } else if (pcouple == YZ) {
1073 double ave = 0.5 * (tensor[1] + tensor[2]);
1074 p_current[1] = p_current[2] = ave;
1075 p_current[0] = tensor[0];
1076 } else if (pcouple == XZ) {
1077 double ave = 0.5 * (tensor[0] + tensor[2]);
1078 p_current[0] = p_current[2] = ave;
1079 p_current[1] = tensor[1];
1080 } else {
1081 p_current[0] = tensor[0];
1082 p_current[1] = tensor[1];
1083 p_current[2] = tensor[2];
1084 }
1085
1086 if (!std::isfinite(p_current[0]) || !std::isfinite(p_current[1]) || !std::isfinite(p_current[2]))
1087 error->all(FLERR,"Non-numeric pressure - simulation unstable");
1088
1089 // switch order from xy-xz-yz to Voigt
1090
1091 if (pstyle == TRICLINIC) {
1092 p_current[3] = tensor[5];
1093 p_current[4] = tensor[4];
1094 p_current[5] = tensor[3];
1095
1096 if (!std::isfinite(p_current[3]) || !std::isfinite(p_current[4]) || !std::isfinite(p_current[5]))
1097 error->all(FLERR,"Non-numeric pressure - simulation unstable");
1098 }
1099 }
1100
1101 /* ----------------------------------------------------------------------
1102 change box size
1103 remap all atoms or dilate group atoms depending on allremap flag
1104 if rigid bodies exist, scale rigid body centers-of-mass
1105 ------------------------------------------------------------------------- */
1106
remap()1107 void FixTGNHDrude::remap()
1108 {
1109 int i;
1110 double oldlo,oldhi;
1111 double expfac;
1112
1113 int nlocal = atom->nlocal;
1114 double *h = domain->h;
1115
1116 // omega is not used, except for book-keeping
1117
1118 for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
1119
1120 // convert pertinent atoms and rigid bodies to lamda coords
1121
1122 domain->x2lamda(nlocal);
1123
1124 if (nrigid)
1125 for (i = 0; i < nrigid; i++)
1126 modify->fix[rfix[i]]->deform(0);
1127
1128 // reset global and local box to new size/shape
1129
1130 // this operation corresponds to applying the
1131 // translate and scale operations
1132 // corresponding to the solution of the following ODE:
1133 //
1134 // h_dot = omega_dot * h
1135 //
1136 // where h_dot, omega_dot and h are all upper-triangular
1137 // 3x3 tensors. In Voigt notation, the elements of the
1138 // RHS product tensor are:
1139 // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
1140 //
1141 // Ordering of operations preserves time symmetry.
1142
1143 double dto2 = dto/2.0;
1144 double dto4 = dto/4.0;
1145 double dto8 = dto/8.0;
1146
1147 // off-diagonal components, first half
1148
1149 if (pstyle == TRICLINIC) {
1150
1151 if (p_flag[4]) {
1152 expfac = exp(dto8*omega_dot[0]);
1153 h[4] *= expfac;
1154 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1155 h[4] *= expfac;
1156 }
1157
1158 if (p_flag[3]) {
1159 expfac = exp(dto4*omega_dot[1]);
1160 h[3] *= expfac;
1161 h[3] += dto2*(omega_dot[3]*h[2]);
1162 h[3] *= expfac;
1163 }
1164
1165 if (p_flag[5]) {
1166 expfac = exp(dto4*omega_dot[0]);
1167 h[5] *= expfac;
1168 h[5] += dto2*(omega_dot[5]*h[1]);
1169 h[5] *= expfac;
1170 }
1171
1172 if (p_flag[4]) {
1173 expfac = exp(dto8*omega_dot[0]);
1174 h[4] *= expfac;
1175 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1176 h[4] *= expfac;
1177 }
1178 }
1179
1180 // scale diagonal components
1181 // scale tilt factors with cell, if set
1182
1183 if (p_flag[0]) {
1184 oldlo = domain->boxlo[0];
1185 oldhi = domain->boxhi[0];
1186 expfac = exp(dto*omega_dot[0]);
1187 domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
1188 domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
1189 }
1190
1191 if (p_flag[1]) {
1192 oldlo = domain->boxlo[1];
1193 oldhi = domain->boxhi[1];
1194 expfac = exp(dto*omega_dot[1]);
1195 domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
1196 domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
1197 if (scalexy) h[5] *= expfac;
1198 }
1199
1200 if (p_flag[2]) {
1201 oldlo = domain->boxlo[2];
1202 oldhi = domain->boxhi[2];
1203 expfac = exp(dto*omega_dot[2]);
1204 domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
1205 domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
1206 if (scalexz) h[4] *= expfac;
1207 if (scaleyz) h[3] *= expfac;
1208 }
1209
1210 // off-diagonal components, second half
1211
1212 if (pstyle == TRICLINIC) {
1213
1214 if (p_flag[4]) {
1215 expfac = exp(dto8*omega_dot[0]);
1216 h[4] *= expfac;
1217 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1218 h[4] *= expfac;
1219 }
1220
1221 if (p_flag[3]) {
1222 expfac = exp(dto4*omega_dot[1]);
1223 h[3] *= expfac;
1224 h[3] += dto2*(omega_dot[3]*h[2]);
1225 h[3] *= expfac;
1226 }
1227
1228 if (p_flag[5]) {
1229 expfac = exp(dto4*omega_dot[0]);
1230 h[5] *= expfac;
1231 h[5] += dto2*(omega_dot[5]*h[1]);
1232 h[5] *= expfac;
1233 }
1234
1235 if (p_flag[4]) {
1236 expfac = exp(dto8*omega_dot[0]);
1237 h[4] *= expfac;
1238 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
1239 h[4] *= expfac;
1240 }
1241
1242 }
1243
1244 domain->yz = h[3];
1245 domain->xz = h[4];
1246 domain->xy = h[5];
1247
1248 // tilt factor to cell length ratio can not exceed TILTMAX in one step
1249
1250 if (domain->yz < -TILTMAX*domain->yprd ||
1251 domain->yz > TILTMAX*domain->yprd ||
1252 domain->xz < -TILTMAX*domain->xprd ||
1253 domain->xz > TILTMAX*domain->xprd ||
1254 domain->xy < -TILTMAX*domain->xprd ||
1255 domain->xy > TILTMAX*domain->xprd)
1256 error->all(FLERR,"Fix npt/nph has tilted box too far in one step - "
1257 "periodic cell is too far from equilibrium state");
1258
1259 domain->set_global_box();
1260 domain->set_local_box();
1261
1262 // convert pertinent atoms and rigid bodies back to box coords
1263
1264 domain->lamda2x(nlocal);
1265
1266 if (nrigid)
1267 for (i = 0; i < nrigid; i++)
1268 modify->fix[rfix[i]]->deform(1);
1269 }
1270
1271 /* ----------------------------------------------------------------------
1272 pack entire state of Fix into one write
1273 ------------------------------------------------------------------------- */
1274
write_restart(FILE * fp)1275 void FixTGNHDrude::write_restart(FILE *fp)
1276 {
1277 int nsize = size_restart_global();
1278
1279 double *list;
1280 memory->create(list,nsize,"nh:list");
1281
1282 pack_restart_data(list);
1283
1284 if (comm->me == 0) {
1285 int size = nsize * sizeof(double);
1286 fwrite(&size,sizeof(int),1,fp);
1287 fwrite(list,sizeof(double),nsize,fp);
1288 }
1289
1290 memory->destroy(list);
1291 }
1292
1293 /* ----------------------------------------------------------------------
1294 calculate the number of data to be packed
1295 ------------------------------------------------------------------------- */
1296
size_restart_global()1297 int FixTGNHDrude::size_restart_global()
1298 {
1299 int nsize = 2;
1300 if (tstat_flag) nsize += 1 + 6*mtchain;
1301 if (pstat_flag) {
1302 nsize += 16 + 2*mpchain;
1303 if (deviatoric_flag) nsize += 6;
1304 }
1305
1306 return nsize;
1307 }
1308
1309 /* ----------------------------------------------------------------------
1310 pack restart data
1311 ------------------------------------------------------------------------- */
1312
pack_restart_data(double * list)1313 int FixTGNHDrude::pack_restart_data(double *list)
1314 {
1315 int n = 0;
1316
1317 list[n++] = tstat_flag;
1318 if (tstat_flag) {
1319 list[n++] = mtchain;
1320 for (int ich = 0; ich < mtchain; ich++) {
1321 list[n++] = etamol[ich];
1322 list[n++] = etaint[ich];
1323 list[n++] = etadrude[ich];
1324 }
1325 for (int ich = 0; ich < mtchain; ich++) {
1326 list[n++] = etamol_dot[ich];
1327 list[n++] = etaint_dot[ich];
1328 list[n++] = etadrude_dot[ich];
1329 }
1330 }
1331
1332 list[n++] = pstat_flag;
1333 if (pstat_flag) {
1334 list[n++] = omega[0];
1335 list[n++] = omega[1];
1336 list[n++] = omega[2];
1337 list[n++] = omega[3];
1338 list[n++] = omega[4];
1339 list[n++] = omega[5];
1340 list[n++] = omega_dot[0];
1341 list[n++] = omega_dot[1];
1342 list[n++] = omega_dot[2];
1343 list[n++] = omega_dot[3];
1344 list[n++] = omega_dot[4];
1345 list[n++] = omega_dot[5];
1346 list[n++] = vol0;
1347 list[n++] = t0;
1348 list[n++] = mpchain;
1349 if (mpchain) {
1350 for (int ich = 0; ich < mpchain; ich++)
1351 list[n++] = etap[ich];
1352 for (int ich = 0; ich < mpchain; ich++)
1353 list[n++] = etap_dot[ich];
1354 }
1355
1356 list[n++] = deviatoric_flag;
1357 if (deviatoric_flag) {
1358 list[n++] = h0_inv[0];
1359 list[n++] = h0_inv[1];
1360 list[n++] = h0_inv[2];
1361 list[n++] = h0_inv[3];
1362 list[n++] = h0_inv[4];
1363 list[n++] = h0_inv[5];
1364 }
1365 }
1366
1367 return n;
1368 }
1369
1370 /* ----------------------------------------------------------------------
1371 use state info from restart file to restart the Fix
1372 ------------------------------------------------------------------------- */
1373
restart(char * buf)1374 void FixTGNHDrude::restart(char *buf)
1375 {
1376 int n = 0;
1377 double *list = (double *) buf;
1378 int flag = static_cast<int> (list[n++]);
1379 if (flag) {
1380 int m = static_cast<int> (list[n++]);
1381 if (tstat_flag && m == mtchain) {
1382 for (int ich = 0; ich < mtchain; ich++) {
1383 etamol[ich] = list[n++];
1384 etaint[ich] = list[n++];
1385 etadrude[ich] = list[n++];
1386 }
1387 for (int ich = 0; ich < mtchain; ich++) {
1388 etamol_dot[ich] = list[n++];
1389 etaint_dot[ich] = list[n++];
1390 etadrude_dot[ich] = list[n++];
1391 }
1392 } else n += 2*m;
1393 }
1394 flag = static_cast<int> (list[n++]);
1395 if (flag) {
1396 omega[0] = list[n++];
1397 omega[1] = list[n++];
1398 omega[2] = list[n++];
1399 omega[3] = list[n++];
1400 omega[4] = list[n++];
1401 omega[5] = list[n++];
1402 omega_dot[0] = list[n++];
1403 omega_dot[1] = list[n++];
1404 omega_dot[2] = list[n++];
1405 omega_dot[3] = list[n++];
1406 omega_dot[4] = list[n++];
1407 omega_dot[5] = list[n++];
1408 vol0 = list[n++];
1409 t0 = list[n++];
1410 int m = static_cast<int> (list[n++]);
1411 if (pstat_flag && m == mpchain) {
1412 for (int ich = 0; ich < mpchain; ich++)
1413 etap[ich] = list[n++];
1414 for (int ich = 0; ich < mpchain; ich++)
1415 etap_dot[ich] = list[n++];
1416 } else n+=2*m;
1417 flag = static_cast<int> (list[n++]);
1418 if (flag) {
1419 h0_inv[0] = list[n++];
1420 h0_inv[1] = list[n++];
1421 h0_inv[2] = list[n++];
1422 h0_inv[3] = list[n++];
1423 h0_inv[4] = list[n++];
1424 h0_inv[5] = list[n++];
1425 }
1426 }
1427 }
1428
1429 /* ---------------------------------------------------------------------- */
1430
modify_param(int narg,char ** arg)1431 int FixTGNHDrude::modify_param(int narg, char **arg)
1432 {
1433 if (strcmp(arg[0],"temp") == 0) {
1434 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1435 if (tcomputeflag) {
1436 modify->delete_compute(id_temp);
1437 tcomputeflag = 0;
1438 }
1439 delete [] id_temp;
1440 id_temp = utils::strdup(arg[1]);
1441
1442 int icompute = modify->find_compute(arg[1]);
1443 if (icompute < 0)
1444 error->all(FLERR,"Could not find fix_modify temperature ID");
1445 temperature = modify->compute[icompute];
1446
1447 if (temperature->tempflag == 0)
1448 error->all(FLERR,
1449 "Fix_modify temperature ID does not compute temperature");
1450 if (temperature->igroup != 0 && comm->me == 0)
1451 error->warning(FLERR,"Temperature for fix modify is not for group all");
1452
1453 // reset id_temp of pressure to new temperature ID
1454
1455 if (pstat_flag) {
1456 icompute = modify->find_compute(id_press);
1457 if (icompute < 0)
1458 error->all(FLERR,"Pressure ID for fix modify does not exist");
1459 modify->compute[icompute]->reset_extra_compute_fix(id_temp);
1460 }
1461
1462 return 2;
1463
1464 } else if (strcmp(arg[0],"press") == 0) {
1465 if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
1466 if (!pstat_flag) error->all(FLERR,"Illegal fix_modify command");
1467 if (pcomputeflag) {
1468 modify->delete_compute(id_press);
1469 pcomputeflag = 0;
1470 }
1471 delete [] id_press;
1472 id_press = utils::strdup(arg[1]);
1473
1474 int icompute = modify->find_compute(arg[1]);
1475 if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");
1476 pressure = modify->compute[icompute];
1477
1478 if (pressure->pressflag == 0)
1479 error->all(FLERR,"Fix_modify pressure ID does not compute pressure");
1480 return 2;
1481 }
1482
1483 return 0;
1484 }
1485
1486 /* ---------------------------------------------------------------------- */
1487
compute_scalar()1488 double FixTGNHDrude::compute_scalar()
1489 {
1490 int i;
1491 double volume;
1492 double energy;
1493 double kt = boltz * t_target;
1494 double kt_drude = boltz * tdrude_target;
1495 double lkt_press = 0.0;
1496 int ich;
1497 if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
1498 else volume = domain->xprd * domain->yprd;
1499
1500 energy = 0.0;
1501
1502 // thermostat chain energy is equivalent to Eq. (2) in
1503 // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
1504 // 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),
1505 // where L = tdof
1506 // M = mtchain
1507 // p_eta_k = Q_k*eta_dot[k-1]
1508 // Q_1 = L*k*T/t_freq^2
1509 // Q_k = k*T/t_freq^2, k > 1
1510
1511 if (tstat_flag) {
1512 energy += ke2mol_target * etamol[0] + 0.5 * etamol_mass[0] * etamol_dot[0] * etamol_dot[0];
1513 energy += ke2int_target * etaint[0] + 0.5 * etaint_mass[0] * etaint_dot[0] * etaint_dot[0];
1514 energy += ke2drude_target * etadrude[0] + 0.5 * etadrude_mass[0] * etadrude_dot[0] * etadrude_dot[0];
1515 for (ich = 1; ich < mtchain; ich++) {
1516 energy += kt * etamol[ich] + 0.5*etamol_mass[ich]*etamol_dot[ich]*etamol_dot[ich];
1517 energy += kt * etaint[ich] + 0.5*etaint_mass[ich]*etaint_dot[ich]*etaint_dot[ich];
1518 energy += kt_drude * etadrude[ich] + 0.5*etadrude_mass[ich]*etadrude_dot[ich]*etadrude_dot[ich];
1519 }
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
compute_vector(int n)1567 double FixTGNHDrude::compute_vector(int n)
1568 {
1569 if (!temp_computed_end_of_step)
1570 compute_temp_mol_int_drude(true);
1571 switch (n) {
1572 case 0:
1573 return t_mol;
1574 case 1:
1575 return t_int;
1576 case 2:
1577 return t_drude;
1578 default:
1579 return 0.0;
1580 }
1581 }
1582
1583 /* ---------------------------------------------------------------------- */
1584
reset_target(double t_new)1585 void FixTGNHDrude::reset_target(double t_new)
1586 {
1587 t_target = t_start = t_stop = t_new;
1588 }
1589
1590 /* ---------------------------------------------------------------------- */
1591
reset_dt()1592 void FixTGNHDrude::reset_dt()
1593 {
1594 dtv = update->dt;
1595 dtf = 0.5 * update->dt * force->ftm2v;
1596 dthalf = 0.5 * update->dt;
1597 dt4 = 0.25 * update->dt;
1598 dt8 = 0.125 * update->dt;
1599 dto = dthalf;
1600
1601 // If using respa, then remap is performed in innermost level
1602
1603 if (utils::strmatch(update->integrate_style,"^respa"))
1604 dto = 0.5*step_respa[0];
1605 }
1606
compute_temp_mol_int_drude(bool end_of_step)1607 void FixTGNHDrude::compute_temp_mol_int_drude(bool end_of_step) {
1608 double **v = atom->v;
1609 double *mass = atom->mass;
1610 tagint *molecule = atom->molecule;
1611 int *type = atom->type;
1612 int *mask = atom->mask;
1613 int *drudetype = fix_drude->drudetype;
1614 tagint *drudeid = fix_drude->drudeid;
1615 int imol, ci, di;
1616 double mass_com, mass_reduced, mass_core, mass_drude;
1617 double vint, vcom, vrel;
1618 // use array instead of two numbers to save MPI_Allreduce()
1619 double ke2_int_drude_tmp[2] = {0.0, 0.0};
1620 double ke2_int_drude[2];
1621
1622 memset(*v_mol_tmp, 0, sizeof(double) * (n_mol + 1) * 3); // the length of v_mol is n_mol+1
1623
1624 /**
1625 * If there are velocity bias, need to remove them before calculate kinetic energies
1626 */
1627 for (int i = 0; i < atom->nlocal; i++) {
1628 if (mask[i] & groupbit) {
1629 if (which == BIAS)
1630 temperature->remove_bias(i, v[i]);
1631
1632 imol = molecule[i];
1633 for (int k = 0; k < 3; k++)
1634 v_mol_tmp[imol][k] += v[i][k] * mass[type[i]];
1635
1636 if (which == BIAS)
1637 temperature->restore_bias(i, v[i]);
1638 }
1639 }
1640 MPI_Allreduce(*v_mol_tmp, *v_mol, (n_mol + 1) * 3, MPI_DOUBLE, MPI_SUM, world);
1641
1642 ke2mol = 0;
1643 for (int i = 1; i < n_mol + 1; i++) {
1644 for (int k = 0; k < 3; k++) {
1645 v_mol[i][k] /= mass_mol[i];
1646 ke2mol += mass_mol[i] * (v_mol[i][k] * v_mol[i][k]);
1647 }
1648 }
1649 ke2mol *= force->mvv2e;
1650 t_mol = ke2mol / dof_mol / boltz;
1651
1652 /**
1653 * Have to call remove_bias at the innermost loop, because drude atom may be a ghost
1654 */
1655 for (int i = 0; i < atom->nlocal; i++) {
1656 if (mask[i] & groupbit) {
1657 imol = molecule[i];
1658 if (drudetype[type[i]] == NOPOL_TYPE) {
1659 if (which == BIAS)
1660 temperature->remove_bias(i, v[i]);
1661 for (int k = 0; k < 3; k++) {
1662 vint = v[i][k] - v_mol[imol][k];
1663 ke2_int_drude_tmp[0] += mass[type[i]] * vint * vint;
1664 }
1665 if (which == BIAS)
1666 temperature->restore_bias(i, v[i]);
1667 } else if (drudetype[type[i]] == CORE_TYPE) {
1668 /**
1669 * have to use closet_image()
1670 * even though all images have the same velocity and it's sort of read-only
1671 * but the bias velocity may depends on it's position like in compute vis/pp
1672 */
1673 ci = i;
1674 di = domain->closest_image(i, atom->map(drudeid[i]));
1675 if (which == BIAS) {
1676 temperature->remove_bias(ci, v[ci]);
1677 temperature->remove_bias(di, v[di]);
1678 }
1679 mass_core = mass[type[ci]];
1680 mass_drude = mass[type[di]];
1681 mass_com = mass_core + mass_drude;
1682 mass_reduced = mass_core * mass_drude / mass_com;
1683 for (int k = 0; k < 3; k++) {
1684 vcom = (mass_core * v[ci][k] + mass_drude * v[di][k]) / mass_com;
1685 vint = vcom - v_mol[imol][k];
1686 ke2_int_drude_tmp[0] += mass_com * vint * vint;
1687 vrel = v[di][k] - v[ci][k];
1688 ke2_int_drude_tmp[1] += mass_reduced * vrel * vrel;
1689 }
1690 if (which == BIAS) {
1691 temperature->restore_bias(ci, v[ci]);
1692 temperature->restore_bias(di, v[di]);
1693 }
1694 }
1695 }
1696 }
1697 MPI_Allreduce(ke2_int_drude_tmp, ke2_int_drude, 2, MPI_DOUBLE, MPI_SUM, world);
1698 ke2int = ke2_int_drude[0] * force->mvv2e;
1699 ke2drude = ke2_int_drude[1] * force->mvv2e;
1700 t_int = ke2int / dof_int / boltz;
1701 t_drude = ke2drude / dof_drude / boltz;
1702
1703 temp_computed_end_of_step = end_of_step;
1704 }
1705
1706 /* ----------------------------------------------------------------------
1707 perform half-step update of chain thermostat variables
1708 ------------------------------------------------------------------------- */
1709
nhc_temp_integrate()1710 void FixTGNHDrude::nhc_temp_integrate()
1711 {
1712 compute_temp_mol_int_drude(false);
1713
1714 // update masses of thermostat in case target temperature changes
1715 etamol_mass[0] = ke2mol_target / (t_freq*t_freq);
1716 etaint_mass[0] = ke2int_target / (t_freq*t_freq);
1717 for (int ich = 1; ich < mtchain; ich++) {
1718 etamol_mass[ich] = boltz * t_target / (t_freq*t_freq);
1719 etaint_mass[ich] = boltz * t_target / (t_freq*t_freq);
1720 }
1721
1722 // thermostat for molecular COM
1723 factor_eta_mol = propagate(etamol, etamol_dot, etamol_dotdot, etamol_mass,
1724 ke2mol, ke2mol_target, t_target);
1725 factor_eta_int = propagate(etaint, etaint_dot, etaint_dotdot, etaint_mass,
1726 ke2int, ke2int_target, t_target);
1727 factor_eta_drude = propagate(etadrude, etadrude_dot, etadrude_dotdot, etadrude_mass,
1728 ke2drude, ke2drude_target, tdrude_target);
1729
1730 nh_v_temp();
1731 }
1732
propagate(double * eta,double * eta_dot,double * eta_dotdot,const double * eta_mass,const double & ke2,const double & ke2_target,const double & tt) const1733 double FixTGNHDrude::propagate(double *eta, double *eta_dot, double *eta_dotdot, const double *eta_mass,
1734 const double &ke2, const double &ke2_target, const double &tt) const {
1735 int ich;
1736 double expfac;
1737 double ncfac = 1.0 / nc_tchain;
1738 double factor_eta = 1.0;
1739
1740 eta_dotdot[0] = (ke2 - ke2_target) / eta_mass[0];
1741 for (int iloop = 0; iloop < nc_tchain; iloop++) {
1742 for (ich = mtchain - 1; ich > 0; ich--) {
1743 expfac = exp(-ncfac * dt8 * eta_dot[ich + 1]);
1744 eta_dot[ich] *= expfac;
1745 eta_dot[ich] += eta_dotdot[ich] * ncfac * dt4;
1746 eta_dot[ich] *= expfac;
1747 }
1748 expfac = exp(-ncfac * dt8 * eta_dot[1]);
1749 eta_dot[0] *= expfac;
1750 eta_dot[0] += eta_dotdot[0] * ncfac * dt4;
1751 eta_dot[0] *= expfac;
1752 factor_eta *= exp(-ncfac * dthalf * eta_dot[0]);
1753
1754 for (ich = 0; ich < mtchain; ich++)
1755 eta[ich] += ncfac * dthalf * eta_dot[ich];
1756
1757 eta_dotdot[0] = (ke2 * factor_eta * factor_eta - ke2_target) / eta_mass[0];
1758 eta_dot[0] *= expfac;
1759 eta_dot[0] += eta_dotdot[0] * ncfac * dt4;
1760 eta_dot[0] *= expfac;
1761 for (ich = 1; ich < mtchain; ich++) {
1762 expfac = exp(-ncfac * dt8 * eta_dot[ich + 1]);
1763 eta_dot[ich] *= expfac;
1764 eta_dotdot[ich] = (eta_mass[ich - 1] * eta_dot[ich - 1] * eta_dot[ich - 1]
1765 - boltz * tt) / eta_mass[ich];
1766 eta_dot[ich] += eta_dotdot[ich] * ncfac * dt4;
1767 eta_dot[ich] *= expfac;
1768 }
1769 }
1770 return factor_eta;
1771 }
1772
1773 /* ----------------------------------------------------------------------
1774 perform half-step update of chain thermostat variables for barostat
1775 scale barostat velocities
1776 ------------------------------------------------------------------------- */
1777
nhc_press_integrate()1778 void FixTGNHDrude::nhc_press_integrate()
1779 {
1780 int ich,i,pdof;
1781 double expfac,factor_etap,kecurrent;
1782 double kt = boltz * t_target;
1783 double lkt_press;
1784
1785 // Update masses, to preserve initial freq, if t_target changed
1786 double nkt = (atom->natoms + 1) * kt;
1787 for (int i = 0; i < 3; i++)
1788 if (p_flag[i])
1789 omega_mass[i] = nkt / (p_freq[i] * p_freq[i]);
1790 if (pstyle == TRICLINIC) {
1791 for (int i = 3; i < 6; i++)
1792 if (p_flag[i]) omega_mass[i] = nkt / (p_freq[i] * p_freq[i]);
1793 }
1794 if (mpchain) {
1795 etap_mass[0] = kt / (p_freq_max * p_freq_max);
1796 for (int ich = 1; ich < mpchain; ich++)
1797 etap_mass[ich] = kt / (p_freq_max * p_freq_max);
1798 for (int ich = 1; ich < mpchain; ich++)
1799 etap_dotdot[ich] = (etap_mass[ich - 1] * etap_dot[ich - 1] * etap_dot[ich - 1]
1800 - kt) / etap_mass[ich];
1801 }
1802
1803 kecurrent = 0.0;
1804 pdof = 0;
1805 for (i = 0; i < 3; i++)
1806 if (p_flag[i]) {
1807 kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1808 pdof++;
1809 }
1810
1811 if (pstyle == TRICLINIC) {
1812 for (i = 3; i < 6; i++)
1813 if (p_flag[i]) {
1814 kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1815 pdof++;
1816 }
1817 }
1818
1819 if (pstyle == ISO) lkt_press = kt;
1820 else lkt_press = pdof * kt;
1821 etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
1822
1823 double ncfac = 1.0/nc_pchain;
1824 for (int iloop = 0; iloop < nc_pchain; iloop++) {
1825
1826 for (ich = mpchain-1; ich > 0; ich--) {
1827 expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
1828 etap_dot[ich] *= expfac;
1829 etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
1830 etap_dot[ich] *= expfac;
1831 }
1832
1833 expfac = exp(-ncfac*dt8*etap_dot[1]);
1834 etap_dot[0] *= expfac;
1835 etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
1836 etap_dot[0] *= expfac;
1837
1838 for (ich = 0; ich < mpchain; ich++)
1839 etap[ich] += ncfac*dthalf*etap_dot[ich];
1840
1841 factor_etap = exp(-ncfac*dthalf*etap_dot[0]);
1842 for (i = 0; i < 3; i++)
1843 if (p_flag[i]) omega_dot[i] *= factor_etap;
1844
1845 if (pstyle == TRICLINIC) {
1846 for (i = 3; i < 6; i++)
1847 if (p_flag[i]) omega_dot[i] *= factor_etap;
1848 }
1849
1850 kecurrent = 0.0;
1851 for (i = 0; i < 3; i++)
1852 if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1853
1854 if (pstyle == TRICLINIC) {
1855 for (i = 3; i < 6; i++)
1856 if (p_flag[i]) kecurrent += omega_mass[i]*omega_dot[i]*omega_dot[i];
1857 }
1858
1859 etap_dotdot[0] = (kecurrent - lkt_press)/etap_mass[0];
1860
1861 etap_dot[0] *= expfac;
1862 etap_dot[0] += etap_dotdot[0] * ncfac*dt4;
1863 etap_dot[0] *= expfac;
1864
1865 for (ich = 1; ich < mpchain; ich++) {
1866 expfac = exp(-ncfac*dt8*etap_dot[ich+1]);
1867 etap_dot[ich] *= expfac;
1868 etap_dotdot[ich] =
1869 (etap_mass[ich-1]*etap_dot[ich-1]*etap_dot[ich-1] - kt) / etap_mass[ich];
1870 etap_dot[ich] += etap_dotdot[ich] * ncfac*dt4;
1871 etap_dot[ich] *= expfac;
1872 }
1873 }
1874 }
1875
1876 /* ----------------------------------------------------------------------
1877 perform half-step barostat scaling of velocities
1878 -----------------------------------------------------------------------*/
1879
nh_v_press()1880 void FixTGNHDrude::nh_v_press()
1881 {
1882 double factor[3];
1883 double **v = atom->v;
1884 int *mask = atom->mask;
1885 int nlocal = atom->nlocal;
1886 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
1887
1888 factor[0] = exp(-dt4*(omega_dot[0]+mtk_term2));
1889 factor[1] = exp(-dt4*(omega_dot[1]+mtk_term2));
1890 factor[2] = exp(-dt4*(omega_dot[2]+mtk_term2));
1891
1892 if (which == NOBIAS) {
1893 for (int i = 0; i < nlocal; i++) {
1894 if (mask[i] & groupbit) {
1895 v[i][0] *= factor[0];
1896 v[i][1] *= factor[1];
1897 v[i][2] *= factor[2];
1898 if (pstyle == TRICLINIC) {
1899 v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
1900 v[i][1] += -dthalf*v[i][2]*omega_dot[3];
1901 }
1902 v[i][0] *= factor[0];
1903 v[i][1] *= factor[1];
1904 v[i][2] *= factor[2];
1905 }
1906 }
1907 } else if (which == BIAS) {
1908 for (int i = 0; i < nlocal; i++) {
1909 if (mask[i] & groupbit) {
1910 temperature->remove_bias(i,v[i]);
1911 v[i][0] *= factor[0];
1912 v[i][1] *= factor[1];
1913 v[i][2] *= factor[2];
1914 if (pstyle == TRICLINIC) {
1915 v[i][0] += -dthalf*(v[i][1]*omega_dot[5] + v[i][2]*omega_dot[4]);
1916 v[i][1] += -dthalf*v[i][2]*omega_dot[3];
1917 }
1918 v[i][0] *= factor[0];
1919 v[i][1] *= factor[1];
1920 v[i][2] *= factor[2];
1921 temperature->restore_bias(i,v[i]);
1922 }
1923 }
1924 }
1925 }
1926
1927 /* ----------------------------------------------------------------------
1928 perform half-step update of velocities
1929 -----------------------------------------------------------------------*/
1930
nve_v()1931 void FixTGNHDrude::nve_v()
1932 {
1933 double dtfm;
1934 double **v = atom->v;
1935 double **f = atom->f;
1936 double *rmass = atom->rmass;
1937 double *mass = atom->mass;
1938 int *type = atom->type;
1939 int *mask = atom->mask;
1940 int nlocal = atom->nlocal;
1941 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
1942
1943 if (rmass) {
1944 for (int i = 0; i < nlocal; i++) {
1945 if (mask[i] & groupbit) {
1946 dtfm = dtf / rmass[i];
1947 v[i][0] += dtfm*f[i][0];
1948 v[i][1] += dtfm*f[i][1];
1949 v[i][2] += dtfm*f[i][2];
1950 }
1951 }
1952 } else {
1953 for (int i = 0; i < nlocal; i++) {
1954 if (mask[i] & groupbit) {
1955 dtfm = dtf / mass[type[i]];
1956 v[i][0] += dtfm*f[i][0];
1957 v[i][1] += dtfm*f[i][1];
1958 v[i][2] += dtfm*f[i][2];
1959 }
1960 }
1961 }
1962 }
1963
1964 /* ----------------------------------------------------------------------
1965 perform full-step update of positions
1966 -----------------------------------------------------------------------*/
1967
nve_x()1968 void FixTGNHDrude::nve_x()
1969 {
1970 double **x = atom->x;
1971 double **v = atom->v;
1972 int *mask = atom->mask;
1973 int nlocal = atom->nlocal;
1974 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
1975
1976 // x update by full step only for atoms in group
1977
1978 for (int i = 0; i < nlocal; i++) {
1979 if (mask[i] & groupbit) {
1980 x[i][0] += dtv * v[i][0];
1981 x[i][1] += dtv * v[i][1];
1982 x[i][2] += dtv * v[i][2];
1983 }
1984 }
1985 }
1986
1987 /* ----------------------------------------------------------------------
1988 perform half-step thermostat scaling of velocities
1989 -----------------------------------------------------------------------*/
1990
nh_v_temp()1991 void FixTGNHDrude::nh_v_temp()
1992 {
1993 double **v = atom->v;
1994 double *mass = atom->mass;
1995 int *mask = atom->mask;
1996 int *type = atom->type;
1997 tagint *molecule = atom->molecule;
1998 int *drudetype = fix_drude->drudetype;
1999 tagint *drudeid = fix_drude->drudeid;
2000
2001 int imol, i, j, ci, di, itype;
2002 double mass_com, mass_core, mass_drude;
2003 double vint, vcom, vrel;
2004
2005 /**
2006 * If there are velocity bias, need to remove them before scale velocity
2007 * Have to call remove_bias at the innermost loop, because drude atom may be a ghost
2008 */
2009 for (i = 0; i < atom->nlocal; i++) {
2010 if (mask[i] & groupbit) {
2011 imol = molecule[i];
2012 itype = drudetype[type[i]];
2013 if (itype == NOPOL_TYPE) {
2014 if (which == BIAS)
2015 temperature->remove_bias(i, v[i]);
2016 for (int k = 0; k < 3; k++) {
2017 vint = v[i][k] - v_mol[imol][k];
2018 vint *= factor_eta_int;
2019 v[i][k] = v_mol[imol][k] * factor_eta_mol + vint;
2020 }
2021 if (which == BIAS)
2022 temperature->restore_bias(i, v[i]);
2023 } else {
2024 // have to use closest_image() because we are manipulating the velocity
2025 j = domain->closest_image(i, atom->map(drudeid[i]));
2026 if (itype == DRUDE_TYPE && j < atom->nlocal) continue;
2027 if (itype == CORE_TYPE) {
2028 ci = i;
2029 di = j;
2030 } else {
2031 ci = j;
2032 di = i;
2033 }
2034 if (which == BIAS) {
2035 temperature->remove_bias(ci, v[ci]);
2036 temperature->remove_bias(di, v[di]);
2037 }
2038 mass_core = mass[type[ci]];
2039 mass_drude = mass[type[di]];
2040 mass_com = mass_core + mass_drude;
2041 for (int k = 0; k < 3; k++) {
2042 vcom = (mass_core * v[ci][k] + mass_drude * v[di][k]) / mass_com;
2043 vint = vcom - v_mol[imol][k];
2044 vrel = v[di][k] - v[ci][k];
2045 vint *= factor_eta_int;
2046 vrel *= factor_eta_drude;
2047 v[ci][k] = v_mol[imol][k] * factor_eta_mol + vint - vrel * mass_drude / mass_com;
2048 v[di][k] = v_mol[imol][k] * factor_eta_mol + vint + vrel * mass_core / mass_com;
2049 }
2050 if (which == BIAS) {
2051 temperature->restore_bias(ci, v[ci]);
2052 temperature->restore_bias(di, v[di]);
2053 }
2054 }
2055 }
2056 }
2057 }
2058
2059 /* ----------------------------------------------------------------------
2060 compute sigma tensor
2061 needed whenever p_target or h0_inv changes
2062 -----------------------------------------------------------------------*/
2063
compute_sigma()2064 void FixTGNHDrude::compute_sigma()
2065 {
2066 // if nreset_h0 > 0, reset vol0 and h0_inv
2067 // every nreset_h0 timesteps
2068
2069 if (nreset_h0 > 0) {
2070 int delta = update->ntimestep - update->beginstep;
2071 if (delta % nreset_h0 == 0) {
2072 if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
2073 else vol0 = domain->xprd * domain->yprd;
2074 h0_inv[0] = domain->h_inv[0];
2075 h0_inv[1] = domain->h_inv[1];
2076 h0_inv[2] = domain->h_inv[2];
2077 h0_inv[3] = domain->h_inv[3];
2078 h0_inv[4] = domain->h_inv[4];
2079 h0_inv[5] = domain->h_inv[5];
2080 }
2081 }
2082
2083 // generate upper-triangular half of
2084 // sigma = vol0*h0inv*(p_target-p_hydro)*h0inv^t
2085 // units of sigma are are PV/L^2 e.g. atm.A
2086 //
2087 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
2088 // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
2089 // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
2090
2091 sigma[0] =
2092 vol0*(h0_inv[0]*((p_target[0]-p_hydro)*h0_inv[0] +
2093 p_target[5]*h0_inv[5]+p_target[4]*h0_inv[4]) +
2094 h0_inv[5]*(p_target[5]*h0_inv[0] +
2095 (p_target[1]-p_hydro)*h0_inv[5]+p_target[3]*h0_inv[4]) +
2096 h0_inv[4]*(p_target[4]*h0_inv[0]+p_target[3]*h0_inv[5] +
2097 (p_target[2]-p_hydro)*h0_inv[4]));
2098 sigma[1] =
2099 vol0*(h0_inv[1]*((p_target[1]-p_hydro)*h0_inv[1] +
2100 p_target[3]*h0_inv[3]) +
2101 h0_inv[3]*(p_target[3]*h0_inv[1] +
2102 (p_target[2]-p_hydro)*h0_inv[3]));
2103 sigma[2] =
2104 vol0*(h0_inv[2]*((p_target[2]-p_hydro)*h0_inv[2]));
2105 sigma[3] =
2106 vol0*(h0_inv[1]*(p_target[3]*h0_inv[2]) +
2107 h0_inv[3]*((p_target[2]-p_hydro)*h0_inv[2]));
2108 sigma[4] =
2109 vol0*(h0_inv[0]*(p_target[4]*h0_inv[2]) +
2110 h0_inv[5]*(p_target[3]*h0_inv[2]) +
2111 h0_inv[4]*((p_target[2]-p_hydro)*h0_inv[2]));
2112 sigma[5] =
2113 vol0*(h0_inv[0]*(p_target[5]*h0_inv[1]+p_target[4]*h0_inv[3]) +
2114 h0_inv[5]*((p_target[1]-p_hydro)*h0_inv[1]+p_target[3]*h0_inv[3]) +
2115 h0_inv[4]*(p_target[3]*h0_inv[1]+(p_target[2]-p_hydro)*h0_inv[3]));
2116 }
2117
2118 /* ----------------------------------------------------------------------
2119 compute strain energy
2120 -----------------------------------------------------------------------*/
2121
compute_strain_energy()2122 double FixTGNHDrude::compute_strain_energy()
2123 {
2124 // compute strain energy = 0.5*Tr(sigma*h*h^t) in energy units
2125
2126 double* h = domain->h;
2127 double d0,d1,d2;
2128
2129 d0 =
2130 sigma[0]*(h[0]*h[0]+h[5]*h[5]+h[4]*h[4]) +
2131 sigma[5]*( h[1]*h[5]+h[3]*h[4]) +
2132 sigma[4]*( h[2]*h[4]);
2133 d1 =
2134 sigma[5]*( h[5]*h[1]+h[4]*h[3]) +
2135 sigma[1]*( h[1]*h[1]+h[3]*h[3]) +
2136 sigma[3]*( h[2]*h[3]);
2137 d2 =
2138 sigma[4]*( h[4]*h[2]) +
2139 sigma[3]*( h[3]*h[2]) +
2140 sigma[2]*( h[2]*h[2]);
2141
2142 double energy = 0.5*(d0+d1+d2)/nktv2p;
2143 return energy;
2144 }
2145
2146 /* ----------------------------------------------------------------------
2147 compute deviatoric barostat force = h*sigma*h^t
2148 -----------------------------------------------------------------------*/
2149
compute_deviatoric()2150 void FixTGNHDrude::compute_deviatoric()
2151 {
2152 // generate upper-triangular part of h*sigma*h^t
2153 // units of fdev are are PV, e.g. atm*A^3
2154 // [ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
2155 // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
2156 // [ 4 3 2 ] [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
2157
2158 double* h = domain->h;
2159
2160 fdev[0] =
2161 h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) +
2162 h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) +
2163 h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]);
2164 fdev[1] =
2165 h[1]*( sigma[1]*h[1]+sigma[3]*h[3]) +
2166 h[3]*( sigma[3]*h[1]+sigma[2]*h[3]);
2167 fdev[2] =
2168 h[2]*( sigma[2]*h[2]);
2169 fdev[3] =
2170 h[1]*( sigma[3]*h[2]) +
2171 h[3]*( sigma[2]*h[2]);
2172 fdev[4] =
2173 h[0]*( sigma[4]*h[2]) +
2174 h[5]*( sigma[3]*h[2]) +
2175 h[4]*( sigma[2]*h[2]);
2176 fdev[5] =
2177 h[0]*( sigma[5]*h[1]+sigma[4]*h[3]) +
2178 h[5]*( sigma[1]*h[1]+sigma[3]*h[3]) +
2179 h[4]*( sigma[3]*h[1]+sigma[2]*h[3]);
2180 }
2181
2182 /* ----------------------------------------------------------------------
2183 compute target temperature and kinetic energy
2184 -----------------------------------------------------------------------*/
2185
compute_temp_target()2186 void FixTGNHDrude::compute_temp_target()
2187 {
2188 double delta = update->ntimestep - update->beginstep;
2189 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2190
2191 t_target = t_start + delta * (t_stop-t_start);
2192 ke2mol_target = dof_mol * boltz * t_target;
2193 ke2int_target = dof_int * boltz * t_target;
2194 ke2drude_target = dof_drude * boltz * tdrude_target;
2195 }
2196
2197 /* ----------------------------------------------------------------------
2198 compute hydrostatic target pressure
2199 -----------------------------------------------------------------------*/
2200
compute_press_target()2201 void FixTGNHDrude::compute_press_target()
2202 {
2203 double delta = update->ntimestep - update->beginstep;
2204 if (delta != 0.0) delta /= update->endstep - update->beginstep;
2205
2206 p_hydro = 0.0;
2207 for (int i = 0; i < 3; i++)
2208 if (p_flag[i]) {
2209 p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
2210 p_hydro += p_target[i];
2211 }
2212 if (pdim > 0) p_hydro /= pdim;
2213
2214 if (pstyle == TRICLINIC)
2215 for (int i = 3; i < 6; i++)
2216 p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
2217
2218 // if deviatoric, recompute sigma each time p_target changes
2219
2220 if (deviatoric_flag) compute_sigma();
2221 }
2222
2223 /* ----------------------------------------------------------------------
2224 update omega_dot, omega
2225 -----------------------------------------------------------------------*/
2226
nh_omega_dot()2227 void FixTGNHDrude::nh_omega_dot()
2228 {
2229 double f_omega,volume;
2230
2231 if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
2232 else volume = domain->xprd*domain->yprd;
2233
2234 if (deviatoric_flag) compute_deviatoric();
2235
2236 mtk_term1 = 0.0;
2237 if (mtk_flag) {
2238 if (pstyle == ISO) {
2239 mtk_term1 = tdof * boltz * t_current;
2240 mtk_term1 /= pdim * atom->natoms;
2241 } else {
2242 double *mvv_current = temperature->vector;
2243 for (int i = 0; i < 3; i++)
2244 if (p_flag[i])
2245 mtk_term1 += mvv_current[i];
2246 mtk_term1 /= pdim * atom->natoms;
2247 }
2248 }
2249
2250 for (int i = 0; i < 3; i++)
2251 if (p_flag[i]) {
2252 f_omega = (p_current[i]-p_hydro)*volume /
2253 (omega_mass[i] * nktv2p) + mtk_term1 / omega_mass[i];
2254 if (deviatoric_flag) f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
2255 omega_dot[i] += f_omega*dthalf;
2256 }
2257
2258 mtk_term2 = 0.0;
2259 if (mtk_flag) {
2260 for (int i = 0; i < 3; i++)
2261 if (p_flag[i])
2262 mtk_term2 += omega_dot[i];
2263 if (pdim > 0) mtk_term2 /= pdim * atom->natoms;
2264 }
2265
2266 if (pstyle == TRICLINIC) {
2267 for (int i = 3; i < 6; i++) {
2268 if (p_flag[i]) {
2269 f_omega = p_current[i]*volume/(omega_mass[i] * nktv2p);
2270 if (deviatoric_flag)
2271 f_omega -= fdev[i]/(omega_mass[i] * nktv2p);
2272 omega_dot[i] += f_omega*dthalf;
2273 }
2274 }
2275 }
2276 }
2277
2278 /* ----------------------------------------------------------------------
2279 if any tilt ratios exceed limits, set flip = 1 and compute new tilt values
2280 do not flip in x or y if non-periodic (can tilt but not flip)
2281 this is b/c the box length would be changed (dramatically) by flip
2282 if yz tilt exceeded, adjust C vector by one B vector
2283 if xz tilt exceeded, adjust C vector by one A vector
2284 if xy tilt exceeded, adjust B vector by one A vector
2285 check yz first since it may change xz, then xz check comes after
2286 if any flip occurs, create new box in domain
2287 image_flip() adjusts image flags due to box shape change induced by flip
2288 remap() puts atoms outside the new box back into the new box
2289 perform irregular on atoms in lamda coords to migrate atoms to new procs
2290 important that image_flip comes before remap, since remap may change
2291 image flags to new values, making eqs in doc of Domain:image_flip incorrect
2292 ------------------------------------------------------------------------- */
2293
pre_exchange()2294 void FixTGNHDrude::pre_exchange()
2295 {
2296 double xprd = domain->xprd;
2297 double yprd = domain->yprd;
2298
2299 // flip is only triggered when tilt exceeds 0.5 by DELTAFLIP
2300 // this avoids immediate re-flipping due to tilt oscillations
2301
2302 double xtiltmax = (0.5+DELTAFLIP)*xprd;
2303 double ytiltmax = (0.5+DELTAFLIP)*yprd;
2304
2305 int flipxy,flipxz,flipyz;
2306 flipxy = flipxz = flipyz = 0;
2307
2308 if (domain->yperiodic) {
2309 if (domain->yz < -ytiltmax) {
2310 domain->yz += yprd;
2311 domain->xz += domain->xy;
2312 flipyz = 1;
2313 } else if (domain->yz >= ytiltmax) {
2314 domain->yz -= yprd;
2315 domain->xz -= domain->xy;
2316 flipyz = -1;
2317 }
2318 }
2319
2320 if (domain->xperiodic) {
2321 if (domain->xz < -xtiltmax) {
2322 domain->xz += xprd;
2323 flipxz = 1;
2324 } else if (domain->xz >= xtiltmax) {
2325 domain->xz -= xprd;
2326 flipxz = -1;
2327 }
2328 if (domain->xy < -xtiltmax) {
2329 domain->xy += xprd;
2330 flipxy = 1;
2331 } else if (domain->xy >= xtiltmax) {
2332 domain->xy -= xprd;
2333 flipxy = -1;
2334 }
2335 }
2336
2337 int flip = 0;
2338 if (flipxy || flipxz || flipyz) flip = 1;
2339
2340 if (flip) {
2341 domain->set_global_box();
2342 domain->set_local_box();
2343
2344 domain->image_flip(flipxy,flipxz,flipyz);
2345
2346 double **x = atom->x;
2347 imageint *image = atom->image;
2348 int nlocal = atom->nlocal;
2349 for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
2350
2351 domain->x2lamda(atom->nlocal);
2352 irregular->migrate_atoms();
2353 domain->lamda2x(atom->nlocal);
2354 }
2355 }
2356
2357 /* ----------------------------------------------------------------------
2358 memory usage of Irregular
2359 ------------------------------------------------------------------------- */
2360
memory_usage()2361 double FixTGNHDrude::memory_usage()
2362 {
2363 double bytes = 0.0;
2364 if (irregular) bytes += irregular->memory_usage();
2365 return bytes;
2366 }
2367