1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing author: W. Michael Brown (Intel)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_nh_intel.h"
20
21 #include "atom.h"
22 #include "domain.h"
23 #include "error.h"
24 #include "force.h"
25 #include "memory.h"
26 #include "modify.h"
27 #include "neighbor.h"
28 #include "update.h"
29
30 #include <cmath>
31
32 using namespace LAMMPS_NS;
33 using namespace FixConst;
34
35 #define TILTMAX 1.5
36
37 enum{NOBIAS,BIAS};
38 enum{ISO,ANISO,TRICLINIC};
39
40 typedef struct { double x,y,z; } dbl3_t;
41
42 /* ----------------------------------------------------------------------
43 NVT,NPH,NPT integrators for improved Nose-Hoover equations of motion
44 ---------------------------------------------------------------------- */
45
FixNHIntel(LAMMPS * lmp,int narg,char ** arg)46 FixNHIntel::FixNHIntel(LAMMPS *lmp, int narg, char **arg) :
47 FixNH(lmp, narg, arg)
48 {
49 _dtfm = 0;
50 _nlocal3 = 0;
51 _nlocal_max = 0;
52 }
53
54 /* ---------------------------------------------------------------------- */
55
~FixNHIntel()56 FixNHIntel::~FixNHIntel()
57 {
58 }
59
60 /* ---------------------------------------------------------------------- */
61
setup(int vflag)62 void FixNHIntel::setup(int vflag)
63 {
64 FixNH::setup(vflag);
65 reset_dt();
66 }
67
68 /* ----------------------------------------------------------------------
69 change box size
70 remap all atoms or dilate group atoms depending on allremap flag
71 if rigid bodies exist, scale rigid body centers-of-mass
72 ------------------------------------------------------------------------- */
73
remap()74 void FixNHIntel::remap()
75 {
76 double oldlo,oldhi;
77 double expfac;
78
79 dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
80 int *mask = atom->mask;
81 int nlocal = atom->nlocal;
82 double *h = domain->h;
83
84 // omega is not used, except for book-keeping
85
86 for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
87
88 // convert pertinent atoms and rigid bodies to lamda coords
89 const double hi0 = domain->h_inv[0];
90 const double hi1 = domain->h_inv[1];
91 const double hi2 = domain->h_inv[2];
92 const double hi3 = domain->h_inv[3];
93 const double hi4 = domain->h_inv[4];
94 const double hi5 = domain->h_inv[5];
95 const double b0 = domain->boxlo[0];
96 const double b1 = domain->boxlo[1];
97 const double b2 = domain->boxlo[2];
98
99 if (allremap) {
100 #if defined(LMP_SIMD_COMPILER)
101 #if defined(USE_OMP_SIMD)
102 #pragma omp simd
103 #else
104 #pragma simd
105 #endif
106 #pragma vector aligned
107 #endif
108 for (int i = 0; i < nlocal; i++) {
109 const double d0 = x[i].x - b0;
110 const double d1 = x[i].y - b1;
111 const double d2 = x[i].z - b2;
112 x[i].x = hi0*d0 + hi5*d1 + hi4*d2;
113 x[i].y = hi1*d1 + hi3*d2;
114 x[i].z = hi2*d2;
115 }
116 } else {
117 #if defined(LMP_SIMD_COMPILER)
118 #if defined(USE_OMP_SIMD)
119 #pragma omp simd
120 #else
121 #pragma simd
122 #endif
123 #pragma vector aligned
124 #endif
125 for (int i = 0; i < nlocal; i++) {
126 if (mask[i] & dilate_group_bit) {
127 const double d0 = x[i].x - b0;
128 const double d1 = x[i].y - b1;
129 const double d2 = x[i].z - b2;
130 x[i].x = hi0*d0 + hi5*d1 + hi4*d2;
131 x[i].y = hi1*d1 + hi3*d2;
132 x[i].z = hi2*d2;
133 }
134 }
135 }
136
137 if (nrigid)
138 for (int i = 0; i < nrigid; i++)
139 modify->fix[rfix[i]]->deform(0);
140
141 // reset global and local box to new size/shape
142
143 // this operation corresponds to applying the
144 // translate and scale operations
145 // corresponding to the solution of the following ODE:
146 //
147 // h_dot = omega_dot * h
148 //
149 // where h_dot, omega_dot and h are all upper-triangular
150 // 3x3 tensors. In Voigt notation, the elements of the
151 // RHS product tensor are:
152 // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
153 //
154 // Ordering of operations preserves time symmetry.
155
156 double dto2 = dto/2.0;
157 double dto4 = dto/4.0;
158 double dto8 = dto/8.0;
159
160 // off-diagonal components, first half
161
162 if (pstyle == TRICLINIC) {
163
164 if (p_flag[4]) {
165 expfac = exp(dto8*omega_dot[0]);
166 h[4] *= expfac;
167 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
168 h[4] *= expfac;
169 }
170
171 if (p_flag[3]) {
172 expfac = exp(dto4*omega_dot[1]);
173 h[3] *= expfac;
174 h[3] += dto2*(omega_dot[3]*h[2]);
175 h[3] *= expfac;
176 }
177
178 if (p_flag[5]) {
179 expfac = exp(dto4*omega_dot[0]);
180 h[5] *= expfac;
181 h[5] += dto2*(omega_dot[5]*h[1]);
182 h[5] *= expfac;
183 }
184
185 if (p_flag[4]) {
186 expfac = exp(dto8*omega_dot[0]);
187 h[4] *= expfac;
188 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
189 h[4] *= expfac;
190 }
191 }
192
193 // scale diagonal components
194 // scale tilt factors with cell, if set
195
196 if (p_flag[0]) {
197 oldlo = domain->boxlo[0];
198 oldhi = domain->boxhi[0];
199 expfac = exp(dto*omega_dot[0]);
200 domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
201 domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
202 }
203
204 if (p_flag[1]) {
205 oldlo = domain->boxlo[1];
206 oldhi = domain->boxhi[1];
207 expfac = exp(dto*omega_dot[1]);
208 domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
209 domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
210 if (scalexy) h[5] *= expfac;
211 }
212
213 if (p_flag[2]) {
214 oldlo = domain->boxlo[2];
215 oldhi = domain->boxhi[2];
216 expfac = exp(dto*omega_dot[2]);
217 domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
218 domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
219 if (scalexz) h[4] *= expfac;
220 if (scaleyz) h[3] *= expfac;
221 }
222
223 // off-diagonal components, second half
224
225 if (pstyle == TRICLINIC) {
226
227 if (p_flag[4]) {
228 expfac = exp(dto8*omega_dot[0]);
229 h[4] *= expfac;
230 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
231 h[4] *= expfac;
232 }
233
234 if (p_flag[3]) {
235 expfac = exp(dto4*omega_dot[1]);
236 h[3] *= expfac;
237 h[3] += dto2*(omega_dot[3]*h[2]);
238 h[3] *= expfac;
239 }
240
241 if (p_flag[5]) {
242 expfac = exp(dto4*omega_dot[0]);
243 h[5] *= expfac;
244 h[5] += dto2*(omega_dot[5]*h[1]);
245 h[5] *= expfac;
246 }
247
248 if (p_flag[4]) {
249 expfac = exp(dto8*omega_dot[0]);
250 h[4] *= expfac;
251 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
252 h[4] *= expfac;
253 }
254
255 }
256
257 domain->yz = h[3];
258 domain->xz = h[4];
259 domain->xy = h[5];
260
261 // tilt factor to cell length ratio can not exceed TILTMAX in one step
262
263 if (domain->yz < -TILTMAX*domain->yprd ||
264 domain->yz > TILTMAX*domain->yprd ||
265 domain->xz < -TILTMAX*domain->xprd ||
266 domain->xz > TILTMAX*domain->xprd ||
267 domain->xy < -TILTMAX*domain->xprd ||
268 domain->xy > TILTMAX*domain->xprd)
269 error->all(FLERR,"Fix npt/nph has tilted box too far in one step - "
270 "periodic cell is too far from equilibrium state");
271
272 domain->set_global_box();
273 domain->set_local_box();
274
275 // convert pertinent atoms and rigid bodies back to box coords
276 const double h0 = domain->h[0];
277 const double h1 = domain->h[1];
278 const double h2 = domain->h[2];
279 const double h3 = domain->h[3];
280 const double h4 = domain->h[4];
281 const double h5 = domain->h[5];
282 const double nb0 = domain->boxlo[0];
283 const double nb1 = domain->boxlo[1];
284 const double nb2 = domain->boxlo[2];
285
286 if (allremap) {
287 #if defined(LMP_SIMD_COMPILER)
288 #if defined(USE_OMP_SIMD)
289 #pragma omp simd
290 #else
291 #pragma simd
292 #endif
293 #pragma vector aligned
294 #endif
295 for (int i = 0; i < nlocal; i++) {
296 x[i].x = h0*x[i].x + h5*x[i].y + h4*x[i].z + nb0;
297 x[i].y = h1*x[i].y + h3*x[i].z + nb1;
298 x[i].z = h2*x[i].z + nb2;
299 }
300 } else {
301 #if defined(LMP_SIMD_COMPILER)
302 #if defined(USE_OMP_SIMD)
303 #pragma omp simd
304 #else
305 #pragma simd
306 #endif
307 #pragma vector aligned
308 #endif
309 for (int i = 0; i < nlocal; i++) {
310 if (mask[i] & dilate_group_bit) {
311 x[i].x = h0*x[i].x + h5*x[i].y + h4*x[i].z + nb0;
312 x[i].y = h1*x[i].y + h3*x[i].z + nb1;
313 x[i].z = h2*x[i].z + nb2;
314 }
315 }
316 }
317
318 if (nrigid)
319 for (int i = 0; i < nrigid; i++)
320 modify->fix[rfix[i]]->deform(1);
321 }
322
323 /* ---------------------------------------------------------------------- */
324
reset_dt()325 void FixNHIntel::reset_dt()
326 {
327 dtv = update->dt;
328 dtf = 0.5 * update->dt * force->ftm2v;
329 dthalf = 0.5 * update->dt;
330 dt4 = 0.25 * update->dt;
331 dt8 = 0.125 * update->dt;
332 dto = dthalf;
333
334 // If using respa, then remap is performed in innermost level
335
336 if (utils::strmatch(update->integrate_style,"^respa"))
337 dto = 0.5*step_respa[0];
338
339 if (pstat_flag)
340 pdrag_factor = 1.0 - (update->dt * p_freq_max * drag / nc_pchain);
341
342 if (tstat_flag)
343 tdrag_factor = 1.0 - (update->dt * t_freq * drag / nc_tchain);
344
345 const int * const mask = atom->mask;
346 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst :
347 atom->nlocal;
348
349 if (nlocal > _nlocal_max) {
350 if (_nlocal_max) memory->destroy(_dtfm);
351 _nlocal_max = static_cast<int>(1.20 * nlocal);
352 memory->create(_dtfm, _nlocal_max * 3, "fix_nh_intel:dtfm");
353 }
354
355 _nlocal3 = nlocal * 3;
356
357 if (igroup == 0) {
358 if (atom->rmass) {
359 const double * const rmass = atom->rmass;
360 int n = 0;
361 for (int i = 0; i < nlocal; i++) {
362 _dtfm[n++] = dtf / rmass[i];
363 _dtfm[n++] = dtf / rmass[i];
364 _dtfm[n++] = dtf / rmass[i];
365 }
366 } else {
367 const double * const mass = atom->mass;
368 const int * const type = atom->type;
369 int n = 0;
370 for (int i = 0; i < nlocal; i++) {
371 _dtfm[n++] = dtf / mass[type[i]];
372 _dtfm[n++] = dtf / mass[type[i]];
373 _dtfm[n++] = dtf / mass[type[i]];
374 }
375 }
376 } else {
377 if (atom->rmass) {
378 const double * const rmass = atom->rmass;
379 int n = 0;
380 for (int i = 0; i < nlocal; i++)
381 if (mask[i] & groupbit) {
382 _dtfm[n++] = dtf / rmass[i];
383 _dtfm[n++] = dtf / rmass[i];
384 _dtfm[n++] = dtf / rmass[i];
385 } else {
386 _dtfm[n++] = 0.0;
387 _dtfm[n++] = 0.0;
388 _dtfm[n++] = 0.0;
389 }
390 } else {
391 const double * const mass = atom->mass;
392 const int * const type = atom->type;
393 int n = 0;
394 for (int i = 0; i < nlocal; i++)
395 if (mask[i] & groupbit) {
396 _dtfm[n++] = dtf / mass[type[i]];
397 _dtfm[n++] = dtf / mass[type[i]];
398 _dtfm[n++] = dtf / mass[type[i]];
399 } else {
400 _dtfm[n++] = 0.0;
401 _dtfm[n++] = 0.0;
402 _dtfm[n++] = 0.0;
403 }
404 }
405 }
406 }
407
408 /* ----------------------------------------------------------------------
409 perform half-step barostat scaling of velocities
410 -----------------------------------------------------------------------*/
411
nh_v_press()412 void FixNHIntel::nh_v_press()
413 {
414 if (pstyle == TRICLINIC || which == BIAS) {
415 FixNH::nh_v_press();
416 return;
417 }
418
419 dbl3_t * _noalias const v = (dbl3_t *)atom->v[0];
420 int *mask = atom->mask;
421 int nlocal = atom->nlocal;
422 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
423
424 double f0 = exp(-dt4*(omega_dot[0]+mtk_term2));
425 double f1 = exp(-dt4*(omega_dot[1]+mtk_term2));
426 double f2 = exp(-dt4*(omega_dot[2]+mtk_term2));
427 f0 *= f0;
428 f1 *= f1;
429 f2 *= f2;
430
431 if (igroup == 0) {
432 #if defined(LMP_SIMD_COMPILER)
433 #if defined(USE_OMP_SIMD)
434 #pragma omp simd
435 #else
436 #pragma simd
437 #endif
438 #pragma vector aligned
439 #endif
440 for (int i = 0; i < nlocal; i++) {
441 v[i].x *= f0;
442 v[i].y *= f1;
443 v[i].z *= f2;
444 }
445 } else {
446 #if defined(LMP_SIMD_COMPILER)
447 #if defined(USE_OMP_SIMD)
448 #pragma omp simd
449 #else
450 #pragma simd
451 #endif
452 #pragma vector aligned
453 #endif
454 for (int i = 0; i < nlocal; i++) {
455 if (mask[i] & groupbit) {
456 v[i].x *= f0;
457 v[i].y *= f1;
458 v[i].z *= f2;
459 }
460 }
461 }
462 }
463
464 /* ----------------------------------------------------------------------
465 perform half-step update of velocities
466 -----------------------------------------------------------------------*/
467
nve_v()468 void FixNHIntel::nve_v()
469 {
470 if (neighbor->ago == 0) reset_dt();
471 double * _noalias const v = atom->v[0];
472 const double * _noalias const f = atom->f[0];
473 #if defined(LMP_SIMD_COMPILER)
474 #if defined(USE_OMP_SIMD)
475 #pragma omp simd
476 #else
477 #pragma simd
478 #endif
479 #pragma vector aligned
480 #endif
481 for (int i = 0; i < _nlocal3; i++)
482 v[i] += _dtfm[i] * f[i];
483 }
484
485 /* ----------------------------------------------------------------------
486 perform full-step update of positions
487 -----------------------------------------------------------------------*/
488
nve_x()489 void FixNHIntel::nve_x()
490 {
491 double * _noalias const x = atom->x[0];
492 double * _noalias const v = atom->v[0];
493
494 // x update by full step only for atoms in group
495
496 if (igroup == 0) {
497 #if defined(LMP_SIMD_COMPILER)
498 #if defined(USE_OMP_SIMD)
499 #pragma omp simd
500 #else
501 #pragma simd
502 #endif
503 #pragma vector aligned
504 #endif
505 for (int i = 0; i < _nlocal3; i++)
506 x[i] += dtv * v[i];
507 } else {
508 #if defined(LMP_SIMD_COMPILER)
509 #if defined(USE_OMP_SIMD)
510 #pragma omp simd
511 #else
512 #pragma simd
513 #endif
514 #pragma vector aligned
515 #endif
516 for (int i = 0; i < _nlocal3; i++) {
517 if (_dtfm[i] != 0.0)
518 x[i] += dtv * v[i];
519 }
520 }
521 }
522
523 /* ----------------------------------------------------------------------
524 perform half-step thermostat scaling of velocities
525 -----------------------------------------------------------------------*/
526
nh_v_temp()527 void FixNHIntel::nh_v_temp()
528 {
529 if (which == BIAS) {
530 FixNH::nh_v_temp();
531 return;
532 }
533
534 double * _noalias const v = atom->v[0];
535
536 if (igroup == 0) {
537 #if defined(LMP_SIMD_COMPILER)
538 #if defined(USE_OMP_SIMD)
539 #pragma omp simd
540 #else
541 #pragma simd
542 #endif
543 #pragma vector aligned
544 #endif
545 for (int i = 0; i < _nlocal3; i++)
546 v[i] *= factor_eta;
547 } else {
548 #if defined(LMP_SIMD_COMPILER)
549 #if defined(USE_OMP_SIMD)
550 #pragma omp simd
551 #else
552 #pragma simd
553 #endif
554 #pragma vector aligned
555 #endif
556 for (int i = 0; i < _nlocal3; i++) {
557 if (_dtfm[i] != 0.0)
558 v[i] *= factor_eta;
559 }
560 }
561 }
562
memory_usage()563 double FixNHIntel::memory_usage()
564 {
565 return FixNH::memory_usage() + _nlocal_max * 3 * sizeof(double);
566 }
567