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: Axel Kohlmeyer (Temple U)
17 ------------------------------------------------------------------------- */
18
19 #include "omp_compat.h"
20 #include "fix_nh_omp.h"
21 #include <cmath>
22 #include "atom.h"
23 #include "compute.h"
24 #include "domain.h"
25 #include "error.h"
26 #include "modify.h"
27
28 using namespace LAMMPS_NS;
29 using namespace FixConst;
30
31 enum{NOBIAS,BIAS};
32 enum{ISO,ANISO,TRICLINIC};
33
34 #define TILTMAX 1.5
35
36 typedef struct { double x,y,z; } dbl3_t;
37
38 /* ----------------------------------------------------------------------
39 change box size
40 remap all atoms or dilate group atoms depending on allremap flag
41 if rigid bodies exist, scale rigid body centers-of-mass
42 ------------------------------------------------------------------------- */
43
remap()44 void FixNHOMP::remap()
45 {
46 double oldlo,oldhi,expfac;
47
48 double * const * _noalias const x = atom->x;
49 const int * _noalias const mask = atom->mask;
50 const int nlocal = atom->nlocal;
51 double * _noalias const h = domain->h;
52
53 // omega is not used, except for book-keeping
54
55 for (int i = 0; i < 6; i++) omega[i] += dto*omega_dot[i];
56
57 // convert pertinent atoms and rigid bodies to lamda coords
58
59 if (allremap) domain->x2lamda(nlocal);
60 else {
61 #if defined (_OPENMP)
62 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
63 #endif
64 for (int i = 0; i < nlocal; i++)
65 if (mask[i] & dilate_group_bit)
66 domain->x2lamda(x[i],x[i]);
67 }
68
69 if (nrigid)
70 for (int i = 0; i < nrigid; i++)
71 modify->fix[rfix[i]]->deform(0);
72
73 // reset global and local box to new size/shape
74
75 // this operation corresponds to applying the
76 // translate and scale operations
77 // corresponding to the solution of the following ODE:
78 //
79 // h_dot = omega_dot * h
80 //
81 // where h_dot, omega_dot and h are all upper-triangular
82 // 3x3 tensors. In Voigt notation, the elements of the
83 // RHS product tensor are:
84 // h_dot = [0*0, 1*1, 2*2, 1*3+3*2, 0*4+5*3+4*2, 0*5+5*1]
85 //
86 // Ordering of operations preserves time symmetry.
87
88 const double dto2 = dto/2.0;
89 const double dto4 = dto/4.0;
90 const double dto8 = dto/8.0;
91
92 // off-diagonal components, first half
93
94 if (pstyle == TRICLINIC) {
95
96 if (p_flag[4]) {
97 expfac = exp(dto8*omega_dot[0]);
98 h[4] *= expfac;
99 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
100 h[4] *= expfac;
101 }
102
103 if (p_flag[3]) {
104 expfac = exp(dto4*omega_dot[1]);
105 h[3] *= expfac;
106 h[3] += dto2*(omega_dot[3]*h[2]);
107 h[3] *= expfac;
108 }
109
110 if (p_flag[5]) {
111 expfac = exp(dto4*omega_dot[0]);
112 h[5] *= expfac;
113 h[5] += dto2*(omega_dot[5]*h[1]);
114 h[5] *= expfac;
115 }
116
117 if (p_flag[4]) {
118 expfac = exp(dto8*omega_dot[0]);
119 h[4] *= expfac;
120 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
121 h[4] *= expfac;
122 }
123 }
124
125 // scale diagonal components
126 // scale tilt factors with cell, if set
127
128 if (p_flag[0]) {
129 oldlo = domain->boxlo[0];
130 oldhi = domain->boxhi[0];
131 expfac = exp(dto*omega_dot[0]);
132 domain->boxlo[0] = (oldlo-fixedpoint[0])*expfac + fixedpoint[0];
133 domain->boxhi[0] = (oldhi-fixedpoint[0])*expfac + fixedpoint[0];
134 }
135
136 if (p_flag[1]) {
137 oldlo = domain->boxlo[1];
138 oldhi = domain->boxhi[1];
139 expfac = exp(dto*omega_dot[1]);
140 domain->boxlo[1] = (oldlo-fixedpoint[1])*expfac + fixedpoint[1];
141 domain->boxhi[1] = (oldhi-fixedpoint[1])*expfac + fixedpoint[1];
142 if (scalexy) h[5] *= expfac;
143 }
144
145 if (p_flag[2]) {
146 oldlo = domain->boxlo[2];
147 oldhi = domain->boxhi[2];
148 expfac = exp(dto*omega_dot[2]);
149 domain->boxlo[2] = (oldlo-fixedpoint[2])*expfac + fixedpoint[2];
150 domain->boxhi[2] = (oldhi-fixedpoint[2])*expfac + fixedpoint[2];
151 if (scalexz) h[4] *= expfac;
152 if (scaleyz) h[3] *= expfac;
153 }
154
155 // off-diagonal components, second half
156
157 if (pstyle == TRICLINIC) {
158
159 if (p_flag[4]) {
160 expfac = exp(dto8*omega_dot[0]);
161 h[4] *= expfac;
162 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
163 h[4] *= expfac;
164 }
165
166 if (p_flag[3]) {
167 expfac = exp(dto4*omega_dot[1]);
168 h[3] *= expfac;
169 h[3] += dto2*(omega_dot[3]*h[2]);
170 h[3] *= expfac;
171 }
172
173 if (p_flag[5]) {
174 expfac = exp(dto4*omega_dot[0]);
175 h[5] *= expfac;
176 h[5] += dto2*(omega_dot[5]*h[1]);
177 h[5] *= expfac;
178 }
179
180 if (p_flag[4]) {
181 expfac = exp(dto8*omega_dot[0]);
182 h[4] *= expfac;
183 h[4] += dto4*(omega_dot[5]*h[3]+omega_dot[4]*h[2]);
184 h[4] *= expfac;
185 }
186
187 }
188
189 domain->yz = h[3];
190 domain->xz = h[4];
191 domain->xy = h[5];
192
193 // tilt factor to cell length ratio can not exceed TILTMAX in one step
194
195 if (domain->yz < -TILTMAX*domain->yprd ||
196 domain->yz > TILTMAX*domain->yprd ||
197 domain->xz < -TILTMAX*domain->xprd ||
198 domain->xz > TILTMAX*domain->xprd ||
199 domain->xy < -TILTMAX*domain->xprd ||
200 domain->xy > TILTMAX*domain->xprd)
201 error->all(FLERR,"Fix npt/nph has tilted box too far in one step - "
202 "periodic cell is too far from equilibrium state");
203
204 domain->set_global_box();
205 domain->set_local_box();
206
207 // convert pertinent atoms and rigid bodies back to box coords
208
209 if (allremap) domain->lamda2x(nlocal);
210 else {
211 #if defined (_OPENMP)
212 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
213 #endif
214 for (int i = 0; i < nlocal; i++)
215 if (mask[i] & dilate_group_bit)
216 domain->lamda2x(x[i],x[i]);
217 }
218
219 if (nrigid)
220 for (int i = 0; i < nrigid; i++)
221 modify->fix[rfix[i]]->deform(1);
222 }
223
224
225 /* ----------------------------------------------------------------------
226 perform half-step barostat scaling of velocities
227 -----------------------------------------------------------------------*/
228
nh_v_press()229 void FixNHOMP::nh_v_press()
230 {
231 const double factor0 = exp(-dt4*(omega_dot[0]+mtk_term2));
232 const double factor1 = exp(-dt4*(omega_dot[1]+mtk_term2));
233 const double factor2 = exp(-dt4*(omega_dot[2]+mtk_term2));
234 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
235 const int * _noalias const mask = atom->mask;
236 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
237
238 if (which == NOBIAS) {
239 #if defined(_OPENMP)
240 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
241 #endif
242 for (int i = 0; i < nlocal; i++) {
243 if (mask[i] & groupbit) {
244 v[i].x *= factor0;
245 v[i].y *= factor1;
246 v[i].z *= factor2;
247 if (pstyle == TRICLINIC) {
248 v[i].x += -dthalf*(v[i].y*omega_dot[5] + v[i].z*omega_dot[4]);
249 v[i].y += -dthalf*v[i].z*omega_dot[3];
250 }
251 v[i].x *= factor0;
252 v[i].y *= factor1;
253 v[i].z *= factor2;
254 }
255 }
256 } else if (which == BIAS) {
257 #if defined(_OPENMP)
258 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
259 #endif
260 for (int i = 0; i < nlocal; i++) {
261 double buf[3];
262 if (mask[i] & groupbit) {
263 temperature->remove_bias_thr(i,&v[i].x,buf);
264 v[i].x *= factor0;
265 v[i].y *= factor1;
266 v[i].z *= factor2;
267 if (pstyle == TRICLINIC) {
268 v[i].x += -dthalf*(v[i].y*omega_dot[5] + v[i].z*omega_dot[4]);
269 v[i].y += -dthalf*v[i].z*omega_dot[3];
270 }
271 v[i].x *= factor0;
272 v[i].y *= factor1;
273 v[i].z *= factor2;
274 temperature->restore_bias_thr(i,&v[i].x,buf);
275 }
276 }
277 }
278 }
279
280 /* ----------------------------------------------------------------------
281 perform half-step update of velocities
282 -----------------------------------------------------------------------*/
283
nve_v()284 void FixNHOMP::nve_v()
285 {
286 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
287 const dbl3_t * _noalias const f = (dbl3_t *) atom->f[0];
288 const int * _noalias const mask = atom->mask;
289 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
290
291 if (atom->rmass) {
292 const double * _noalias const rmass = atom->rmass;
293 #if defined(_OPENMP)
294 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
295 #endif
296 for (int i = 0; i < nlocal; i++) {
297 if (mask[i] & groupbit) {
298 const double dtfm = dtf / rmass[i];
299 v[i].x += dtfm*f[i].x;
300 v[i].y += dtfm*f[i].y;
301 v[i].z += dtfm*f[i].z;
302 }
303 }
304 } else {
305 const double *_noalias const mass = atom->mass;
306 const int * _noalias const type = atom->type;
307 #if defined(_OPENMP)
308 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
309 #endif
310 for (int i = 0; i < nlocal; i++) {
311 if (mask[i] & groupbit) {
312 const double dtfm = dtf / mass[type[i]];
313 v[i].x += dtfm*f[i].x;
314 v[i].y += dtfm*f[i].y;
315 v[i].z += dtfm*f[i].z;
316 }
317 }
318 }
319 }
320
321 /* ----------------------------------------------------------------------
322 perform full-step update of positions
323 -----------------------------------------------------------------------*/
324
nve_x()325 void FixNHOMP::nve_x()
326 {
327 dbl3_t * _noalias const x = (dbl3_t *) atom->x[0];
328 const dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
329 const int * _noalias const mask = atom->mask;
330 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
331
332 // x update by full step only for atoms in group
333
334 #if defined(_OPENMP)
335 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
336 #endif
337 for (int i = 0; i < nlocal; i++) {
338 if (mask[i] & groupbit) {
339 x[i].x += dtv * v[i].x;
340 x[i].y += dtv * v[i].y;
341 x[i].z += dtv * v[i].z;
342 }
343 }
344 }
345 /* ----------------------------------------------------------------------
346 perform half-step thermostat scaling of velocities
347 -----------------------------------------------------------------------*/
348
nh_v_temp()349 void FixNHOMP::nh_v_temp()
350 {
351 dbl3_t * _noalias const v = (dbl3_t *) atom->v[0];
352 const int * _noalias const mask = atom->mask;
353 const int nlocal = (igroup == atom->firstgroup) ? atom->nfirst : atom->nlocal;
354
355 if (which == NOBIAS) {
356 #if defined(_OPENMP)
357 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
358 #endif
359 for (int i = 0; i < nlocal; i++) {
360 if (mask[i] & groupbit) {
361 v[i].x *= factor_eta;
362 v[i].y *= factor_eta;
363 v[i].z *= factor_eta;
364 }
365 }
366 } else if (which == BIAS) {
367 #if defined(_OPENMP)
368 #pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
369 #endif
370 for (int i = 0; i < nlocal; i++) {
371 double buf[3];
372 if (mask[i] & groupbit) {
373 temperature->remove_bias_thr(i,&v[i].x,buf);
374 v[i].x *= factor_eta;
375 v[i].y *= factor_eta;
376 v[i].z *= factor_eta;
377 temperature->restore_bias_thr(i,&v[i].x,buf);
378 }
379 }
380 }
381 }
382
383