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: Peter Wirnsberger (University of Cambridge)
17
18 This source file implements the asymmetric version of the enhanced heat
19 exchange (eHEX/a) algorithm. The paper is available for download on
20 arXiv: https://arxiv.org/pdf/1507.07081.pdf.
21
22 This file is based on fix_heat.cpp written by Paul Crozier (SNL)
23 which implements the heat exchange (HEX) algorithm.
24 ------------------------------------------------------------------------- */
25
26 #include "fix_ehex.h"
27
28 #include "atom.h"
29 #include "domain.h"
30 #include "error.h"
31 #include "force.h"
32 #include "group.h"
33 #include "memory.h"
34 #include "modify.h"
35 #include "region.h"
36 #include "update.h"
37
38 #include <cmath>
39 #include <cstring>
40
41 #include "fix_shake.h"
42
43 using namespace LAMMPS_NS;
44 using namespace FixConst;
45
46 enum{CONSTANT,EQUAL,ATOM};
47
48 /* ---------------------------------------------------------------------- */
49
FixEHEX(LAMMPS * lmp,int narg,char ** arg)50 FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
51 idregion(nullptr), x(nullptr), f(nullptr), v(nullptr),
52 mass(nullptr), rmass(nullptr), type(nullptr), scalingmask(nullptr)
53 {
54 MPI_Comm_rank(world, &me);
55
56 // check
57 if (narg < 4) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
58
59 scalar_flag = 1;
60 global_freq = 1;
61 extscalar = 0;
62
63 // apply fix every nevery timesteps
64
65 nevery = utils::inumeric(FLERR,arg[3],false,lmp);
66
67 if (nevery <= 0) error->all(FLERR,"Illegal fix ehex command");
68
69 // heat flux into the reservoir
70
71 heat_input = utils::numeric(FLERR,arg[4],false,lmp);
72
73 // optional args
74
75 iregion = -1;
76
77 // NOTE: constraints are deactivated by default
78
79 constraints = 0;
80
81 // NOTE: cluster rescaling is deactivated by default
82
83 cluster = 0;
84
85 // NOTE: hex = 1 means that no coordinate correction is applied in which case eHEX reduces to HEX
86
87 hex = 0;
88
89 int iarg = 5;
90 while (iarg < narg) {
91
92 if (strcmp(arg[iarg],"region") == 0) {
93 if (iarg+2 > narg) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters ");
94 iregion = domain->find_region(arg[iarg+1]);
95 if (iregion == -1)
96 error->all(FLERR,"Region ID for fix ehex does not exist");
97 idregion = utils::strdup(arg[iarg+1]);
98 iarg += 2;
99 }
100
101 // apply constraints (shake/rattle) at the end of the timestep
102
103 else if (strcmp(arg[iarg], "constrain") == 0) {
104 constraints = 1;
105 iarg += 1;
106 }
107
108 // rescale only if the entire molecule is contained within the region
109
110 else if (strcmp(arg[iarg], "com") == 0) {
111 cluster = 1;
112 iarg += 1;
113 }
114
115 // don't apply a coordinate correction if this keyword is specified
116
117 else if (strcmp(arg[iarg], "hex") == 0) {
118 hex = 1;
119 iarg+= 1;
120 }
121 else
122 error->all(FLERR, "Illegal fix ehex keyword ");
123 }
124
125 // check options
126
127 if (cluster && !constraints)
128 error->all(FLERR, "You can only use the keyword 'com' together with the keyword 'constrain' ");
129
130 scale = 1.0;
131 scalingmask = nullptr;
132 FixEHEX::grow_arrays(atom->nmax);
133 atom->add_callback(Atom::GROW);
134
135 }
136
137
138 /* ---------------------------------------------------------------------- */
139
grow_arrays(int nmax)140 void FixEHEX::grow_arrays(int nmax) {
141 memory->grow(scalingmask, nmax,"ehex:scalingmask");
142 }
143
144
145 /* ---------------------------------------------------------------------- */
146
~FixEHEX()147 FixEHEX::~FixEHEX()
148 {
149 atom->delete_callback(id,Atom::GROW);
150 delete [] idregion;
151 memory->destroy(scalingmask);
152
153 }
154
155 /* ---------------------------------------------------------------------- */
156
setmask()157 int FixEHEX::setmask()
158 {
159 int mask = 0;
160 mask |= END_OF_STEP;
161 return mask;
162 }
163
164 /* ---------------------------------------------------------------------- */
165
init()166 void FixEHEX::init()
167 {
168 // set index and check validity of region
169
170 if (iregion >= 0) {
171 iregion = domain->find_region(idregion);
172 if (iregion == -1)
173 error->all(FLERR,"Region ID for fix ehex does not exist");
174 }
175
176 // cannot have 0 atoms in group
177
178 if (group->count(igroup) == 0)
179 error->all(FLERR,"Fix ehex group has no atoms");
180
181 fshake = nullptr;
182 if (constraints) {
183
184 // check if constraining algorithm is used (FixRattle inherits from FixShake)
185
186 int cnt_shake = 0;
187 int id_shake;
188 for (int i = 0; i < modify->nfix; i++) {
189 if (strcmp("rattle", modify->fix[i]->style) == 0 ||
190 strcmp("shake", modify->fix[i]->style) == 0) {
191 cnt_shake++;
192 id_shake = i;
193 }
194 }
195
196 if (cnt_shake > 1)
197 error->all(FLERR,"Multiple instances of fix shake/rattle detected (not supported yet)");
198 else if (cnt_shake == 1) {
199 fshake = ((FixShake*) modify->fix[id_shake]);
200 }
201 else if (cnt_shake == 0)
202 error->all(FLERR, "Fix ehex was configured with keyword constrain, but shake/rattle was not defined");
203 }
204 }
205
206
207
208 /* ---------------------------------------------------------------------- */
209
210
end_of_step()211 void FixEHEX::end_of_step() {
212 // store local pointers
213
214 x = atom->x;
215 f = atom->f;
216 v = atom->v;
217 mass = atom->mass;
218 rmass = atom->rmass;
219 type = atom->type;
220 nlocal = atom->nlocal;
221
222 // determine which sites are to be rescaled
223
224 update_scalingmask();
225
226 // rescale velocities
227
228 rescale();
229
230 // if required use shake/rattle to correct coordinates and velocities
231
232 if (constraints && fshake)
233 fshake->shake_end_of_step(0);
234 }
235
236
237
238 /* ----------------------------------------------------------------------
239 Iterate over all atoms, rescale the velocities and apply coordinate
240 corrections.
241 ------------------------------------------------------------------------- */
242
rescale()243 void FixEHEX::rescale() {
244 double Kr, Ke, escale;
245 double vsub[3],vcm[3], sfr[3];
246 double mi;
247 double dt;
248 double F, mr, epsr_ik, sfvr, eta_ik;
249
250 dt = update->dt;
251
252 // calculate centre of mass properties
253
254 com_properties(vcm, sfr, &sfvr, &Ke, &Kr, &masstotal);
255
256 // heat flux into the reservoir
257
258 F = heat_input * force->ftm2v * nevery;
259
260 // total mass
261
262 mr = masstotal;
263
264 // energy scaling factor
265
266 escale = 1. + (F*dt)/Kr;
267
268 // safety check for kinetic energy
269
270 if (escale < 0.0) error->all(FLERR,"Fix ehex kinetic energy went negative");
271
272 scale = sqrt(escale);
273 vsub[0] = (scale-1.0) * vcm[0];
274 vsub[1] = (scale-1.0) * vcm[1];
275 vsub[2] = (scale-1.0) * vcm[2];
276
277 for (int i = 0; i < nlocal; i++) {
278
279 if (scalingmask[i]) {
280
281 mi = (rmass) ? rmass[i] : mass[type[i]];
282
283 for (int k=0; k<3; k++) {
284
285 // apply coordinate correction unless running in hex mode
286
287 if (!hex) {
288
289 // epsr_ik implements Eq. (20) in the paper
290
291 eta_ik = mi * F/(2.*Kr) * (v[i][k] - vcm[k]);
292 epsr_ik = eta_ik / (mi*Kr) * (F/48. + sfvr/6.*force->ftm2v) - F/(12.*Kr) * (f[i][k]/mi - sfr[k]/mr)*force->ftm2v;
293
294 x[i][k] -= dt*dt*dt * epsr_ik;
295 }
296
297 // rescale the velocity
298
299 v[i][k] = scale*v[i][k] - vsub[k];
300 }
301 }
302 }
303 }
304
305
306 /* ---------------------------------------------------------------------- */
307
compute_scalar()308 double FixEHEX::compute_scalar()
309 {
310 return scale;
311 }
312
313 /* ----------------------------------------------------------------------
314 memory usage of local atom-based arrays
315 ------------------------------------------------------------------------- */
316
memory_usage()317 double FixEHEX::memory_usage()
318 {
319 double bytes = 0.0;
320 bytes += (double)atom->nmax * sizeof(double);
321 return bytes;
322 }
323
324
325 /* ----------------------------------------------------------------------
326 Update the array scalingmask depending on which individual atoms
327 will be rescaled or not.
328 ------------------------------------------------------------------------- */
329
update_scalingmask()330 void FixEHEX::update_scalingmask() {
331 int m;
332 int lid;
333 bool stat;
334 int nsites;
335
336 // prematch region
337
338 Region *region = nullptr;
339 if (iregion >= 0) {
340 region = domain->regions[iregion];
341 region->prematch();
342 }
343
344 // only rescale molecules whose center of mass if fully contained in the region
345
346 if (cluster) {
347
348 // loop over all clusters
349
350 for (int i=0; i < fshake->nlist; i++) {
351
352 // cluster id
353
354 m = fshake->list[i];
355
356 // check if the centre of mass of the cluster is inside the region
357 // if region == nullptr, just check the group information of all sites
358
359 if (fshake->shake_flag[m] == 1) nsites = 3;
360 else if (fshake->shake_flag[m] == 2) nsites = 2;
361 else if (fshake->shake_flag[m] == 3) nsites = 3;
362 else if (fshake->shake_flag[m] == 4) nsites = 4;
363 else nsites = 0;
364
365 if (nsites == 0) {
366 error->all(FLERR,"Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist");
367 }
368
369 stat = check_cluster(fshake->shake_atom[m], nsites, region);
370
371 for (int l=0; l < nsites; l++) {
372 lid = atom->map(fshake->shake_atom[m][l]);
373 scalingmask[lid] = stat;
374 }
375 }
376
377 // check atoms that do not belong to any cluster
378
379 for (int i=0; i<atom->nlocal; i++) {
380 if (fshake->shake_flag[i] == 0)
381 scalingmask[i] = rescale_atom(i,region);
382 }
383
384 }
385
386 // no clusters, just individual sites (e.g. monatomic system or flexible molecules)
387
388 else {
389 for (int i=0; i<atom->nlocal; i++)
390 scalingmask[i] = rescale_atom(i,region);
391 }
392
393 }
394
395
396 /* ----------------------------------------------------------------------
397 Check if the centre of mass of the cluster to be constrained is
398 inside the region.
399 ------------------------------------------------------------------------- */
400
check_cluster(tagint * shake_atom,int n,Region * region)401 bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) {
402
403 // IMPORTANT NOTE: If any site of the cluster belongs to a group
404 // which should not be rescaled than all of the sites
405 // will be ignored!
406
407 double **x = atom->x;
408 double * rmass = atom->rmass;
409 double * mass = atom->mass;
410 int * type = atom->type;
411 int * mask = atom->mask;
412 double xcom[3], xtemp[3];
413 double mcluster, mi;
414 bool stat;
415 int lid[4];
416
417 // accumulate mass and centre of mass position
418
419 stat = true;
420 xcom[0] = 0.;
421 xcom[1] = 0.;
422 xcom[2] = 0.;
423 mcluster = 0;
424
425 for (int i = 0; i < n; i++) {
426
427 // get local id
428
429 lid[i] = atom->map(shake_atom[i]);
430
431 // check if all sites of the cluster belong to the correct group
432
433 stat = stat && (mask[lid[i]] & groupbit);
434
435 if (region && stat) {
436
437 // check if reduced mass is used
438
439 mi = (rmass) ? rmass[lid[i]] : mass[type[lid[i]]];
440 mcluster += mi;
441
442 // accumulate centre of mass
443 // NOTE: you can either use unwrapped coordinates or take site x[lid[0]] as reference,
444 // i.e. reconstruct the molecule around this site and calculate the com.
445
446 for (int k=0; k<3; k++)
447 xtemp[k] = x[lid[i]][k] - x[lid[0]][k];
448
449 // take into account pbc
450
451 domain->minimum_image(xtemp);
452
453 for (int k=0; k<3; k++)
454 xcom[k] += mi * (x[lid[0]][k] + xtemp[k]) ;
455 }
456 }
457
458 // check if centre of mass is inside the region (if specified)
459
460 if (region && stat) {
461
462 // check mass
463
464 if (mcluster < 1.e-14) {
465 error->all(FLERR, "Fix ehex shake cluster has almost zero mass.");
466 }
467
468 // divide by total mass
469
470 for (int k=0; k<3; k++)
471 xcom[k] = xcom[k]/mcluster;
472
473 // apply periodic boundary conditions (centre of mass could be outside the box)
474 // and check if molecule is inside the region
475
476 domain->remap(xcom);
477 stat = stat && region->match(xcom[0], xcom[1], xcom[2]);
478 }
479
480 return stat;
481 }
482
483
484 /* ----------------------------------------------------------------------
485 Check if atom i has the correct group and is inside the region.
486 ------------------------------------------------------------------------- */
487
rescale_atom(int i,Region * region)488 bool FixEHEX::rescale_atom(int i, Region*region) {
489 bool stat;
490 double x_r[3];
491
492 // check mask and group
493
494 stat = (atom->mask[i] & groupbit);
495
496 if (region) {
497
498 x_r[0] = atom->x[i][0];
499 x_r[1] = atom->x[i][1];
500 x_r[2] = atom->x[i][2];
501
502 // apply periodic boundary conditions
503
504 domain->remap(x_r);
505
506 // check if the atom is in the group/region
507
508 stat = stat && region->match(x_r[0],x_r[1],x_r[2]);
509 }
510
511 return stat;
512 }
513
514 /* ----------------------------------------------------------------------
515 Calculate global properties of the atoms inside the reservoir.
516 (e.g. com velocity, kinetic energy, total mass,...)
517 ------------------------------------------------------------------------- */
518
com_properties(double * vr,double * sfr,double * sfvr,double * K,double * Kr,double * mr)519 void FixEHEX::com_properties(double * vr, double * sfr, double *sfvr, double *K, double *Kr, double *mr) {
520 double ** f = atom->f;
521 double ** v = atom->v;
522 int nlocal = atom->nlocal;
523 double *rmass= atom->rmass;
524 double *mass = atom->mass;
525 int *type = atom->type;
526 double l_vr[3];
527 double l_mr;
528 double l_sfr[3];
529 double l_sfvr;
530 double l_K;
531 double mi;
532 double l_buf[9];
533 double buf[9];
534
535 // calculate partial sums
536
537 l_vr[0] = l_vr[1] = l_vr[2] = 0;
538 l_sfr[0] = l_sfr[1] = l_sfr[2] = 0;
539 l_sfvr = 0;
540 l_mr = 0;
541 l_K = 0;
542
543 for (int i = 0; i < nlocal; i++) {
544 if (scalingmask[i]) {
545
546 // check if reduced mass is used
547
548 mi = (rmass) ? rmass[i] : mass[type[i]];
549
550 // accumulate total mass
551
552 l_mr += mi;
553
554 // accumulate kinetic energy
555
556 l_K += mi/2. * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
557
558 // sum_j f_j * v_j
559
560 l_sfvr += f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2];
561
562 // accumulate com velocity and sum of forces
563
564 for (int k=0; k<3; k++) {
565 l_vr[k] += mi * v[i][k];
566 l_sfr[k]+= f[i][k];
567 }
568 }
569 }
570
571 // reduce sums
572
573 l_buf[0] = l_vr[0];
574 l_buf[1] = l_vr[1];
575 l_buf[2] = l_vr[2];
576 l_buf[3] = l_K;
577 l_buf[4] = l_mr;
578 l_buf[5] = l_sfr[0];
579 l_buf[6] = l_sfr[1];
580 l_buf[7] = l_sfr[2];
581 l_buf[8] = l_sfvr;
582
583 MPI_Allreduce(l_buf, buf, 9, MPI_DOUBLE, MPI_SUM, world);
584
585 // total mass of region
586
587 *mr = buf[4];
588
589 if (*mr < 1.e-14) {
590 error->all(FLERR, "Fix ehex error mass of region is close to zero");
591 }
592
593 // total kinetic energy of region
594
595 *K = buf[3];
596
597 // centre of mass velocity of region
598
599 vr[0] = buf[0]/(*mr);
600 vr[1] = buf[1]/(*mr);
601 vr[2] = buf[2]/(*mr);
602
603 // sum of forces
604
605 sfr[0] = buf[5];
606 sfr[1] = buf[6];
607 sfr[2] = buf[7];
608
609 // calculate non-translational kinetic energy
610
611 *Kr = *K - 0.5* (*mr) * (vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]);
612
613 // calculate sum_j f_j * (v_j - v_r) = sum_j f_j * v_j - v_r * sum_j f_j
614
615 *sfvr = buf[8] - (vr[0]*sfr[0] + vr[1]*sfr[1] + vr[2]*sfr[2]);
616 }
617
618