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
19 #include "fix_rattle.h"
20
21 #include "atom.h"
22 #include "comm.h"
23 #include "domain.h"
24 #include "error.h"
25 #include "force.h"
26 #include "math_extra.h"
27 #include "memory.h"
28 #include "modify.h"
29 #include "update.h"
30
31 #include <cmath>
32 #include <cstring>
33
34 using namespace LAMMPS_NS;
35 using namespace FixConst;
36 using namespace MathExtra;
37
38 // set RATTLE_DEBUG 1 to check constraints at end of timestep
39
40 #define RATTLE_DEBUG 0
41
42 // set RATTLE_RAISE_ERROR 1 if you want this fix to raise
43 // an error if the constraints cannot be satisfied
44
45 #define RATTLE_RAISE_ERROR 0
46
47 // You can enable velocity and coordinate checks separately
48 // by setting RATTLE_TEST_VEL/POS true
49
50 #define RATTLE_TEST_VEL false
51 #define RATTLE_TEST_POS false
52
53 enum{V,VP,XSHAKE};
54
55 /* ---------------------------------------------------------------------- */
56
FixRattle(LAMMPS * lmp,int narg,char ** arg)57 FixRattle::FixRattle(LAMMPS *lmp, int narg, char **arg) :
58 FixShake(lmp, narg, arg)
59 {
60 rattle = 1;
61
62 // allocate memory for unconstrained velocity update
63
64 vp = nullptr;
65 FixRattle::grow_arrays(atom->nmax);
66
67 // default communication mode
68 // necessary for compatibility with SHAKE
69 // see pack_forward and unpack_forward
70
71 comm_mode = XSHAKE;
72 vflag_post_force = 0;
73
74 verr_max = 0;
75 derr_max = 0;
76 }
77
78 /* ---------------------------------------------------------------------- */
79
~FixRattle()80 FixRattle::~FixRattle()
81 {
82 memory->destroy(vp);
83
84 if (RATTLE_DEBUG) {
85
86 // communicate maximum distance error
87
88 double global_derr_max, global_verr_max;
89 int npid;
90
91 MPI_Reduce(&derr_max, &global_derr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
92 MPI_Reduce(&verr_max, &global_verr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
93
94 MPI_Comm_rank (world, &npid); // Find out process rank
95
96 if (npid == 0 && screen) {
97 fprintf(screen, "RATTLE: Maximum overall relative position error ( (r_ij-d_ij)/d_ij ): %.10g\n", global_derr_max);
98 fprintf(screen, "RATTLE: Maximum overall absolute velocity error (r_ij * v_ij): %.10g\n", global_verr_max);
99 }
100 }
101 }
102
103 /* ---------------------------------------------------------------------- */
104
setmask()105 int FixRattle::setmask()
106 {
107 int mask = 0;
108 mask |= PRE_NEIGHBOR;
109 mask |= POST_FORCE;
110 mask |= POST_FORCE_RESPA;
111 mask |= FINAL_INTEGRATE;
112 mask |= FINAL_INTEGRATE_RESPA;
113 if (RATTLE_DEBUG) mask |= END_OF_STEP;
114 return mask;
115 }
116
117 /* ----------------------------------------------------------------------
118 initialize RATTLE and check that this is the last final_integrate fix
119 ------------------------------------------------------------------------- */
120
init()121 void FixRattle::init() {
122
123 // initialize SHAKE first
124
125 FixShake::init();
126
127 // show a warning if any final-integrate fix comes after this one
128
129 int after = 0;
130 int flag = 0;
131 for (int i = 0; i < modify->nfix; i++) {
132 if (strcmp(id,modify->fix[i]->id) == 0) after = 1;
133 else if ((modify->fmask[i] & FINAL_INTEGRATE) && after) flag = 1;
134 }
135
136 if (flag && comm->me == 0)
137 error->warning(FLERR,
138 "Fix rattle should come after all other integration fixes ");
139 }
140
141 /* ----------------------------------------------------------------------
142 This method carries out an unconstrained velocity update first and
143 then applies the velocity corrections directly (v and vp are modified).
144 ------------------------------------------------------------------------- */
145
post_force(int vflag)146 void FixRattle::post_force(int vflag)
147 {
148 // remember vflag for the coordinate correction in this->final_integrate
149
150 vflag_post_force = vflag;
151
152 // unconstrained velocity update by half a timestep
153 // similar to FixShake::unconstrained_update()
154
155 update_v_half_nocons();
156
157 // communicate the unconstrained velocities
158
159 if (nprocs > 1) {
160 comm_mode = VP;
161 comm->forward_comm_fix(this);
162 }
163
164 // correct the velocity for each molecule accordingly
165
166 int m;
167 for (int i = 0; i < nlist; i++) {
168 m = list[i];
169 if (shake_flag[m] == 2) vrattle2(m);
170 else if (shake_flag[m] == 3) vrattle3(m);
171 else if (shake_flag[m] == 4) vrattle4(m);
172 else vrattle3angle(m);
173 }
174 }
175
176 /* ---------------------------------------------------------------------- */
177
post_force_respa(int vflag,int ilevel,int)178 void FixRattle::post_force_respa(int vflag, int ilevel, int /*iloop*/)
179 {
180 // remember vflag for the coordinate correction in this->final_integrate
181
182 vflag_post_force = vflag;
183
184 // unconstrained velocity update by half a timestep
185 // similar to FixShake::unconstrained_update()
186
187 update_v_half_nocons_respa(ilevel);
188
189 // communicate the unconstrained velocities
190
191 if (nprocs > 1) {
192 comm_mode = VP;
193 comm->forward_comm_fix(this);
194 }
195
196 // correct the velocity for each molecule accordingly
197
198 int m;
199 for (int i = 0; i < nlist; i++) {
200 m = list[i];
201 if (shake_flag[m] == 2) vrattle2(m);
202 else if (shake_flag[m] == 3) vrattle3(m);
203 else if (shake_flag[m] == 4) vrattle4(m);
204 else vrattle3angle(m);
205 }
206 }
207
208 /* ----------------------------------------------------------------------
209 let SHAKE calculate the constraining forces for the coordinates
210 ------------------------------------------------------------------------- */
211
final_integrate()212 void FixRattle::final_integrate()
213 {
214 comm_mode = XSHAKE;
215 FixShake::post_force(vflag_post_force);
216 }
217
218 /* ---------------------------------------------------------------------- */
219
final_integrate_respa(int ilevel,int iloop)220 void FixRattle::final_integrate_respa(int ilevel, int iloop)
221 {
222 comm_mode = XSHAKE;
223 FixShake::post_force_respa(vflag_post_force, ilevel, iloop);
224 }
225
226 /* ----------------------------------------------------------------------
227 correct velocities of molecule m with 2 constraints bonds and 1 angle
228 ------------------------------------------------------------------------- */
229
vrattle3angle(int m)230 void FixRattle::vrattle3angle(int m)
231 {
232 tagint i0,i1,i2;
233 double c[3], l[3], a[3][3], r01[3], imass[3],
234 r02[3], r12[3], vp01[3], vp02[3], vp12[3];
235
236 // local atom IDs and constraint distances
237
238 i0 = atom->map(shake_atom[m][0]);
239 i1 = atom->map(shake_atom[m][1]);
240 i2 = atom->map(shake_atom[m][2]);
241
242 // r01,r02,r12 = distance vec between atoms
243
244 MathExtra::sub3(x[i1],x[i0],r01);
245 MathExtra::sub3(x[i2],x[i0],r02);
246 MathExtra::sub3(x[i2],x[i1],r12);
247
248 // take into account periodicity
249
250 domain->minimum_image(r01);
251 domain->minimum_image(r02);
252 domain->minimum_image(r12);
253
254 // v01,v02,v12 = velocity differences
255
256 MathExtra::sub3(vp[i1],vp[i0],vp01);
257 MathExtra::sub3(vp[i2],vp[i0],vp02);
258 MathExtra::sub3(vp[i2],vp[i1],vp12);
259
260 // matrix coeffs and rhs for lamda equations
261
262 if (rmass) {
263 imass[0] = 1.0/rmass[i0];
264 imass[1] = 1.0/rmass[i1];
265 imass[2] = 1.0/rmass[i2];
266 } else {
267 imass[0] = 1.0/mass[type[i0]];
268 imass[1] = 1.0/mass[type[i1]];
269 imass[2] = 1.0/mass[type[i2]];
270 }
271
272 // setup matrix
273
274 a[0][0] = (imass[1] + imass[0]) * MathExtra::dot3(r01,r01);
275 a[0][1] = (imass[0] ) * MathExtra::dot3(r01,r02);
276 a[0][2] = (-imass[1] ) * MathExtra::dot3(r01,r12);
277 a[1][0] = a[0][1];
278 a[1][1] = (imass[0] + imass[2]) * MathExtra::dot3(r02,r02);
279 a[1][2] = (imass[2] ) * MathExtra::dot3(r02,r12);
280 a[2][0] = a[0][2];
281 a[2][1] = a[1][2];
282 a[2][2] = (imass[2] + imass[1]) * MathExtra::dot3(r12,r12);
283
284 // sestup RHS
285
286 c[0] = -MathExtra::dot3(vp01,r01);
287 c[1] = -MathExtra::dot3(vp02,r02);
288 c[2] = -MathExtra::dot3(vp12,r12);
289
290 // calculate the inverse matrix exactly
291
292 solve3x3exactly(a,c,l);
293
294 // add corrections to the velocities if processor owns atom
295
296 if (i0 < nlocal) {
297 for (int k=0; k<3; k++)
298 v[i0][k] -= imass[0]* ( l[0] * r01[k] + l[1] * r02[k] );
299 }
300 if (i1 < nlocal) {
301 for (int k=0; k<3; k++)
302 v[i1][k] -= imass[1] * ( -l[0] * r01[k] + l[2] * r12[k] );
303 }
304 if (i2 < nlocal) {
305 for (int k=0; k<3; k++)
306 v[i2][k] -= imass[2] * ( -l[1] * r02[k] - l[2] * r12[k] );
307 }
308 }
309
310 /* ---------------------------------------------------------------------- */
311
vrattle2(int m)312 void FixRattle::vrattle2(int m)
313 {
314 tagint i0, i1;
315 double imass[2], r01[3], vp01[3];
316
317 // local atom IDs and constraint distances
318
319 i0 = atom->map(shake_atom[m][0]);
320 i1 = atom->map(shake_atom[m][1]);
321
322 // r01 = distance vec between atoms, with PBC
323
324 MathExtra::sub3(x[i1],x[i0],r01);
325 domain->minimum_image(r01);
326
327 // v01 = distance vectors for velocities
328
329 MathExtra::sub3(vp[i1],vp[i0],vp01);
330
331 // matrix coeffs and rhs for lamda equations
332
333 if (rmass) {
334 imass[0] = 1.0/rmass[i0];
335 imass[1] = 1.0/rmass[i1];
336 } else {
337 imass[0] = 1.0/mass[type[i0]];
338 imass[1] = 1.0/mass[type[i1]];
339 }
340
341 // Lagrange multiplier: exact solution
342
343 double l01 = - MathExtra::dot3(r01,vp01) /
344 (MathExtra::dot3(r01,r01) * (imass[0] + imass[1]));
345
346 // add corrections to the velocities if the process owns this atom
347
348 if (i0 < nlocal) {
349 for (int k=0; k<3; k++)
350 v[i0][k] -= imass[0] * l01 * r01[k];
351 }
352 if (i1 < nlocal) {
353 for (int k=0; k<3; k++)
354 v[i1][k] -= imass[1] * (-l01) * r01[k];
355 }
356 }
357
358 /* ---------------------------------------------------------------------- */
359
vrattle3(int m)360 void FixRattle::vrattle3(int m)
361 {
362 tagint i0,i1,i2;
363 double imass[3], r01[3], r02[3], vp01[3], vp02[3],
364 a[2][2],c[2],l[2];
365
366 // local atom IDs and constraint distances
367
368 i0 = atom->map(shake_atom[m][0]);
369 i1 = atom->map(shake_atom[m][1]);
370 i2 = atom->map(shake_atom[m][2]);
371
372 // r01,r02 = distance vec between atoms, with PBC
373
374 MathExtra::sub3(x[i1],x[i0],r01);
375 MathExtra::sub3(x[i2],x[i0],r02);
376
377 domain->minimum_image(r01);
378 domain->minimum_image(r02);
379
380 // vp01,vp02 = distance vectors between velocities
381
382 MathExtra::sub3(vp[i1],vp[i0],vp01);
383 MathExtra::sub3(vp[i2],vp[i0],vp02);
384
385 if (rmass) {
386 imass[0] = 1.0/rmass[i0];
387 imass[1] = 1.0/rmass[i1];
388 imass[2] = 1.0/rmass[i2];
389 } else {
390 imass[0] = 1.0/mass[type[i0]];
391 imass[1] = 1.0/mass[type[i1]];
392 imass[2] = 1.0/mass[type[i2]];
393 }
394
395 // setup matrix
396
397 a[0][0] = (imass[1] + imass[0]) * MathExtra::dot3(r01,r01);
398 a[0][1] = (imass[0] ) * MathExtra::dot3(r01,r02);
399 a[1][0] = a[0][1];
400 a[1][1] = (imass[0] + imass[2]) * MathExtra::dot3(r02,r02);
401
402 // setup RHS
403
404 c[0] = - MathExtra::dot3(vp01,r01);
405 c[1] = - MathExtra::dot3(vp02,r02);
406
407 // calculate the inverse 2x2 matrix exactly
408
409 solve2x2exactly(a,c,l);
410
411 // add corrections to the velocities if the process owns this atom
412
413 if (i0 < nlocal) {
414 for (int k=0; k<3; k++)
415 v[i0][k] -= imass[0] * ( l[0] * r01[k] + l[1] * r02[k] );
416 }
417 if (i1 < nlocal)
418 for (int k=0; k<3; k++) {
419 v[i1][k] -= imass[1] * ( -l[0] * r01[k] );
420 }
421 if (i2 < nlocal) {
422 for (int k=0; k<3; k++)
423 v[i2][k] -= imass[2] * ( -l[1] * r02[k] );
424 }
425 }
426
427 /* ---------------------------------------------------------------------- */
428
vrattle4(int m)429 void FixRattle::vrattle4(int m)
430 {
431 tagint i0,i1,i2,i3;
432 double imass[4], c[3], l[3], a[3][3],
433 r01[3], r02[3], r03[3], vp01[3], vp02[3], vp03[3];
434
435 // local atom IDs and constraint distances
436
437 i0 = atom->map(shake_atom[m][0]);
438 i1 = atom->map(shake_atom[m][1]);
439 i2 = atom->map(shake_atom[m][2]);
440 i3 = atom->map(shake_atom[m][3]);
441
442 // r01,r02,r12 = distance vec between atoms, with PBC
443
444 MathExtra::sub3(x[i1],x[i0],r01);
445 MathExtra::sub3(x[i2],x[i0],r02);
446 MathExtra::sub3(x[i3],x[i0],r03);
447
448 domain->minimum_image(r01);
449 domain->minimum_image(r02);
450 domain->minimum_image(r03);
451
452 // vp01,vp02,vp03 = distance vectors between velocities
453
454 MathExtra::sub3(vp[i1],vp[i0],vp01);
455 MathExtra::sub3(vp[i2],vp[i0],vp02);
456 MathExtra::sub3(vp[i3],vp[i0],vp03);
457
458 // matrix coeffs and rhs for lamda equations
459
460 if (rmass) {
461 imass[0] = 1.0/rmass[i0];
462 imass[1] = 1.0/rmass[i1];
463 imass[2] = 1.0/rmass[i2];
464 imass[3] = 1.0/rmass[i3];
465 } else {
466 imass[0] = 1.0/mass[type[i0]];
467 imass[1] = 1.0/mass[type[i1]];
468 imass[2] = 1.0/mass[type[i2]];
469 imass[3] = 1.0/mass[type[i3]];
470 }
471
472 // setup matrix
473
474 a[0][0] = (imass[0] + imass[1]) * MathExtra::dot3(r01,r01);
475 a[0][1] = (imass[0] ) * MathExtra::dot3(r01,r02);
476 a[0][2] = (imass[0] ) * MathExtra::dot3(r01,r03);
477 a[1][0] = a[0][1];
478 a[1][1] = (imass[0] + imass[2]) * MathExtra::dot3(r02,r02);
479 a[1][2] = (imass[0] ) * MathExtra::dot3(r02,r03);
480 a[2][0] = a[0][2];
481 a[2][1] = a[1][2];
482 a[2][2] = (imass[0] + imass[3]) * MathExtra::dot3(r03,r03);
483
484 // setup RHS
485
486 c[0] = - MathExtra::dot3(vp01,r01);
487 c[1] = - MathExtra::dot3(vp02,r02);
488 c[2] = - MathExtra::dot3(vp03,r03);
489
490 // calculate the inverse 3x3 matrix exactly
491
492 solve3x3exactly(a,c,l);
493
494 // add corrections to the velocities if the process owns this atom
495
496 if (i0 < nlocal) {
497 for (int k=0; k<3; k++)
498 v[i0][k] -= imass[0] * ( l[0] * r01[k] + l[1] * r02[k] + l[2] * r03[k]);
499 }
500 if (i1 < nlocal) {
501 for (int k=0; k<3; k++)
502 v[i1][k] -= imass[1] * (-l[0] * r01[k]);
503 }
504 if (i2 < nlocal) {
505 for (int k=0; k<3; k++)
506 v[i2][k] -= imass[2] * ( -l[1] * r02[k]);
507 }
508 if (i3 < nlocal) {
509 for (int k=0; k<3; k++)
510 v[i3][k] -= imass[3] * ( -l[2] * r03[k]);
511 }
512 }
513
514 /* ---------------------------------------------------------------------- */
515
solve2x2exactly(const double a[][2],const double c[],double l[])516 void FixRattle::solve2x2exactly(const double a[][2],
517 const double c[], double l[])
518 {
519 double determ, determinv;
520
521 // calculate the determinant of the matrix
522
523 determ = a[0][0] * a[1][1] - a[0][1] * a[1][0];
524
525 // check if matrix is actually invertible
526
527 if (determ == 0.0) error->one(FLERR,"Rattle determinant = 0.0");
528 determinv = 1.0/determ;
529
530 // Calculate the solution: (l01, l02)^T = A^(-1) * c
531
532 l[0] = determinv * ( a[1][1] * c[0] - a[0][1] * c[1]);
533 l[1] = determinv * (-a[1][0] * c[0] + a[0][0] * c[1]);
534 }
535
536 /* ---------------------------------------------------------------------- */
537
solve3x3exactly(const double a[][3],const double c[],double l[])538 void FixRattle::solve3x3exactly(const double a[][3],
539 const double c[], double l[])
540 {
541 double ai[3][3];
542 double determ, determinv;
543
544 // calculate the determinant of the matrix
545
546 determ = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] +
547 a[0][2]*a[1][0]*a[2][1] - a[0][0]*a[1][2]*a[2][1] -
548 a[0][1]*a[1][0]*a[2][2] - a[0][2]*a[1][1]*a[2][0];
549
550 // check if matrix is actually invertible
551
552 if (determ == 0.0) error->one(FLERR,"Rattle determinant = 0.0");
553
554 // calculate the inverse 3x3 matrix: A^(-1) = (ai_jk)
555
556 determinv = 1.0/determ;
557 ai[0][0] = determinv * (a[1][1]*a[2][2] - a[1][2]*a[2][1]);
558 ai[0][1] = -determinv * (a[0][1]*a[2][2] - a[0][2]*a[2][1]);
559 ai[0][2] = determinv * (a[0][1]*a[1][2] - a[0][2]*a[1][1]);
560 ai[1][0] = -determinv * (a[1][0]*a[2][2] - a[1][2]*a[2][0]);
561 ai[1][1] = determinv * (a[0][0]*a[2][2] - a[0][2]*a[2][0]);
562 ai[1][2] = -determinv * (a[0][0]*a[1][2] - a[0][2]*a[1][0]);
563 ai[2][0] = determinv * (a[1][0]*a[2][1] - a[1][1]*a[2][0]);
564 ai[2][1] = -determinv * (a[0][0]*a[2][1] - a[0][1]*a[2][0]);
565 ai[2][2] = determinv * (a[0][0]*a[1][1] - a[0][1]*a[1][0]);
566
567 // calculate the solution: (l01, l02, l12)^T = A^(-1) * c
568
569 for (int i=0; i<3; i++) {
570 l[i] = 0;
571 for (int j=0; j<3; j++)
572 l[i] += ai[i][j] * c[j];
573 }
574 }
575
576 /* ---------------------------------------------------------------------- */
577
reset_dt()578 void FixRattle::reset_dt()
579 {
580 FixShake::reset_dt();
581 }
582
583 /* ----------------------------------------------------------------------
584 carry out an unconstrained velocity update (vp is modified)
585 ------------------------------------------------------------------------- */
586
update_v_half_nocons()587 void FixRattle::update_v_half_nocons()
588 {
589 const double dtfv = 0.5 * update->dt * force->ftm2v;
590 double dtfvinvm;
591
592 if (rmass) {
593 for (int i = 0; i < nlocal; i++) {
594 if (shake_flag[i]) {
595 dtfvinvm = dtfv / rmass[i];
596 for (int k=0; k<3; k++)
597 vp[i][k] = v[i][k] + dtfvinvm * f[i][k];
598 }
599 else
600 vp[i][0] = vp[i][1] = vp[i][2] = 0;
601 }
602 }
603 else {
604 for (int i = 0; i < nlocal; i++) {
605 dtfvinvm = dtfv/mass[type[i]];
606 if (shake_flag[i]) {
607 for (int k=0; k<3; k++)
608 vp[i][k] = v[i][k] + dtfvinvm * f[i][k];
609 }
610 else
611 vp[i][0] = vp[i][1] = vp[i][2] = 0;
612 }
613 }
614 }
615
616 /* ---------------------------------------------------------------------- */
617
update_v_half_nocons_respa(int)618 void FixRattle::update_v_half_nocons_respa(int /*ilevel*/)
619 {
620 // carry out unconstrained velocity update
621
622 update_v_half_nocons();
623 }
624
625 /* ----------------------------------------------------------------------
626 memory usage of local atom-based arrays
627 ------------------------------------------------------------------------- */
628
memory_usage()629 double FixRattle::memory_usage()
630 {
631 int nmax = atom->nmax;
632 double bytes = FixShake::memory_usage();
633 bytes += (double)nmax*3 * sizeof(double);
634 return bytes;
635 }
636
637 /* ----------------------------------------------------------------------
638 allocate local atom-based arrays
639 ------------------------------------------------------------------------- */
640
grow_arrays(int nmax)641 void FixRattle::grow_arrays(int nmax)
642 {
643 FixShake::grow_arrays(nmax);
644 memory->destroy(vp);
645 memory->create(vp,nmax,3,"rattle:vp");
646 }
647
648 /* ---------------------------------------------------------------------- */
649
pack_forward_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)650 int FixRattle::pack_forward_comm(int n, int *list, double *buf,
651 int pbc_flag, int *pbc)
652 {
653 int i,j,m;
654 m = 0;
655
656 switch (comm_mode) {
657 case XSHAKE:
658 m = FixShake::pack_forward_comm(n, list, buf, pbc_flag, pbc);
659 break;
660 case VP:
661 for (i = 0; i < n; i++) {
662 j = list[i];
663 buf[m++] = vp[j][0];
664 buf[m++] = vp[j][1];
665 buf[m++] = vp[j][2];
666 }
667 break;
668
669 case V:
670 for (i = 0; i < n; i++) {
671 j = list[i];
672 buf[m++] = v[j][0];
673 buf[m++] = v[j][1];
674 buf[m++] = v[j][2];
675 }
676 break;
677 }
678 return m;
679 }
680
681 /* ---------------------------------------------------------------------- */
682
unpack_forward_comm(int n,int first,double * buf)683 void FixRattle::unpack_forward_comm(int n, int first, double *buf)
684 {
685 int i, m, last;
686 m = 0;
687 last = first + n;
688
689 switch (comm_mode) {
690 case XSHAKE:
691 FixShake::unpack_forward_comm(n, first,buf);
692 break;
693
694 case VP:
695 for (i = first; i < last; i++) {
696 vp[i][0] = buf[m++];
697 vp[i][1] = buf[m++];
698 vp[i][2] = buf[m++];
699 }
700 break;
701
702 case V:
703 for (i = first; i < last; i++) {
704 v[i][0] = buf[m++];
705 v[i][1] = buf[m++];
706 v[i][2] = buf[m++];
707 }
708 break;
709 }
710 }
711
712
713 /* ----------------------------------------------------------------------
714 Let shake calculate new constraining forces for the coordinates;
715 As opposed to the regular shake call, this method is usually called from
716 end_of_step fixes after the second velocity integration has happened.
717 ------------------------------------------------------------------------- */
718
shake_end_of_step(int vflag)719 void FixRattle::shake_end_of_step(int vflag) {
720
721 if (nprocs > 1) {
722 comm_mode = V;
723 comm->forward_comm_fix(this);
724 }
725
726 comm_mode = XSHAKE;
727 FixShake::shake_end_of_step(vflag);
728 }
729
730
731 /* ----------------------------------------------------------------------
732 Let shake calculate new constraining forces and correct the
733 coordinates. Nothing to do for rattle here.
734 ------------------------------------------------------------------------- */
735
correct_coordinates(int vflag)736 void FixRattle::correct_coordinates(int vflag) {
737 comm_mode = XSHAKE;
738 FixShake::correct_coordinates(vflag);
739 }
740
741
742 /* ----------------------------------------------------------------------
743 Remove the velocity component along any bond.
744 ------------------------------------------------------------------------- */
745
correct_velocities()746 void FixRattle::correct_velocities() {
747
748 // Copy current velocities instead of unconstrained_update, because the correction
749 // should happen instantaneously and not after the next half step.
750
751 for (int i = 0; i < atom->nlocal; i++) {
752 if (shake_flag[i]) {
753 for (int k=0; k<3; k++)
754 vp[i][k] = v[i][k];
755 }
756 else
757 vp[i][0] = vp[i][1] = vp[i][2] = 0;
758 }
759
760 // communicate the unconstrained velocities
761
762 if (nprocs > 1) {
763 comm_mode = VP;
764 comm->forward_comm_fix(this);
765 }
766
767 // correct the velocity for each molecule accordingly
768
769 int m;
770 for (int i = 0; i < nlist; i++) {
771 m = list[i];
772 if (shake_flag[m] == 2) vrattle2(m);
773 else if (shake_flag[m] == 3) vrattle3(m);
774 else if (shake_flag[m] == 4) vrattle4(m);
775 else vrattle3angle(m);
776 }
777 }
778
779
780 /* ----------------------------------------------------------------------
781 DEBUGGING methods
782 The functions below allow you to check whether the
783 coordinate and velocity constraints are satisfied at the
784 end of the timestep
785 only enabled if RATTLE_DEBUG is set to 1 at top of file
786 checkX tests if shakeX and vrattleX worked as expected
787 ------------------------------------------------------------------------- */
788
end_of_step()789 void FixRattle::end_of_step()
790 {
791 if (nprocs > 1) {
792 comm_mode = V;
793 comm->forward_comm_fix(this);
794 }
795 if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL) && RATTLE_RAISE_ERROR) {
796 error->one(FLERR, "Rattle failed ");
797 }
798 }
799
800
801 /* ---------------------------------------------------------------------- */
802
check_constraints(double ** v,bool checkr,bool checkv)803 bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
804 {
805 int m;
806 bool ret = true;
807 int i=0;
808 while (i < nlist && ret) {
809 m = list[i];
810 if (shake_flag[m] == 2) ret = check2(v,m,checkr,checkv);
811 else if (shake_flag[m] == 3) ret = check3(v,m,checkr,checkv);
812 else if (shake_flag[m] == 4) ret = check4(v,m,checkr,checkv);
813 else ret = check3angle(v,m,checkr,checkv);
814 i++;
815 if (!RATTLE_RAISE_ERROR) ret = true;
816 }
817 return ret;
818 }
819
820 /* ---------------------------------------------------------------------- */
821
check2(double ** v,int m,bool checkr,bool checkv)822 bool FixRattle::check2(double **v, int m, bool checkr, bool checkv)
823 {
824 bool stat;
825 double r01[3],v01[3];
826 const double tol = tolerance;
827 double bond1 = bond_distance[shake_type[m][0]];
828
829 tagint i0 = atom->map(shake_atom[m][0]);
830 tagint i1 = atom->map(shake_atom[m][1]);
831
832 MathExtra::sub3(x[i1],x[i0],r01);
833 domain->minimum_image(r01);
834 MathExtra::sub3(v[i1],v[i0],v01);
835
836 stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol));
837 if (!stat)
838 error->one(FLERR,"Coordinate constraints are not satisfied "
839 "up to desired tolerance ");
840
841 stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol));
842 if (!stat)
843 error->one(FLERR,"Velocity constraints are not satisfied "
844 "up to desired tolerance ");
845 return stat;
846 }
847
848 /* ---------------------------------------------------------------------- */
849
check3(double ** v,int m,bool checkr,bool checkv)850 bool FixRattle::check3(double **v, int m, bool checkr, bool checkv)
851 {
852 bool stat;
853 tagint i0,i1,i2;
854 double r01[3], r02[3], v01[3], v02[3];
855 const double tol = tolerance;
856 double bond1 = bond_distance[shake_type[m][0]];
857 double bond2 = bond_distance[shake_type[m][1]];
858
859 i0 = atom->map(shake_atom[m][0]);
860 i1 = atom->map(shake_atom[m][1]);
861 i2 = atom->map(shake_atom[m][2]);
862
863 MathExtra::sub3(x[i1],x[i0],r01);
864 MathExtra::sub3(x[i2],x[i0],r02);
865
866 domain->minimum_image(r01);
867 domain->minimum_image(r02);
868
869 MathExtra::sub3(v[i1],v[i0],v01);
870 MathExtra::sub3(v[i2],v[i0],v02);
871
872 stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
873 fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol));
874 if (!stat)
875 error->one(FLERR,"Coordinate constraints are not satisfied "
876 "up to desired tolerance ");
877
878 stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
879 fabs(MathExtra::dot3(r02,v02)) > tol));
880 if (!stat)
881 error->one(FLERR,"Velocity constraints are not satisfied "
882 "up to desired tolerance ");
883 return stat;
884 }
885
886 /* ---------------------------------------------------------------------- */
887
check4(double ** v,int m,bool checkr,bool checkv)888 bool FixRattle::check4(double **v, int m, bool checkr, bool checkv)
889 {
890 bool stat = true;
891 const double tol = tolerance;
892 double r01[3], r02[3], r03[3], v01[3], v02[3], v03[3];
893
894 int i0 = atom->map(shake_atom[m][0]);
895 int i1 = atom->map(shake_atom[m][1]);
896 int i2 = atom->map(shake_atom[m][2]);
897 int i3 = atom->map(shake_atom[m][3]);
898 double bond1 = bond_distance[shake_type[m][0]];
899 double bond2 = bond_distance[shake_type[m][1]];
900 double bond3 = bond_distance[shake_type[m][2]];
901
902 MathExtra::sub3(x[i1],x[i0],r01);
903 MathExtra::sub3(x[i2],x[i0],r02);
904 MathExtra::sub3(x[i3],x[i0],r03);
905
906 domain->minimum_image(r01);
907 domain->minimum_image(r02);
908 domain->minimum_image(r03);
909
910 MathExtra::sub3(v[i1],v[i0],v01);
911 MathExtra::sub3(v[i2],v[i0],v02);
912 MathExtra::sub3(v[i3],v[i0],v03);
913
914 stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
915 fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol ||
916 fabs(sqrt(MathExtra::dot3(r03,r03))-bond3) > tol));
917 if (!stat)
918 error->one(FLERR,"Coordinate constraints are not satisfied "
919 "up to desired tolerance ");
920
921 stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
922 fabs(MathExtra::dot3(r02,v02)) > tol ||
923 fabs(MathExtra::dot3(r03,v03)) > tol));
924 if (!stat)
925 error->one(FLERR,"Velocity constraints are not satisfied "
926 "up to desired tolerance ");
927 return stat;
928 }
929
930 /* ---------------------------------------------------------------------- */
931
check3angle(double ** v,int m,bool checkr,bool checkv)932 bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
933 {
934 bool stat = true;
935 const double tol = tolerance;
936 double r01[3], r02[3], r12[3], v01[3], v02[3], v12[3];
937
938 int i0 = atom->map(shake_atom[m][0]);
939 int i1 = atom->map(shake_atom[m][1]);
940 int i2 = atom->map(shake_atom[m][2]);
941 double bond1 = bond_distance[shake_type[m][0]];
942 double bond2 = bond_distance[shake_type[m][1]];
943 double bond12 = angle_distance[shake_type[m][2]];
944
945 MathExtra::sub3(x[i1],x[i0],r01);
946 MathExtra::sub3(x[i2],x[i0],r02);
947 MathExtra::sub3(x[i2],x[i1],r12);
948
949 domain->minimum_image(r01);
950 domain->minimum_image(r02);
951 domain->minimum_image(r12);
952
953 MathExtra::sub3(v[i1],v[i0],v01);
954 MathExtra::sub3(v[i2],v[i0],v02);
955 MathExtra::sub3(v[i2],v[i1],v12);
956
957
958
959 double db1 = fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1);
960 double db2 = fabs(sqrt(MathExtra::dot3(r02,r02))-bond2);
961 double db12 = fabs(sqrt(MathExtra::dot3(r12,r12))-bond12);
962
963
964 stat = !(checkr && (db1 > tol ||
965 db2 > tol ||
966 db12 > tol));
967
968 if (derr_max < db1/bond1) derr_max = db1/bond1;
969 if (derr_max < db2/bond2) derr_max = db2/bond2;
970 if (derr_max < db12/bond12) derr_max = db12/bond12;
971
972
973 if (!stat && RATTLE_RAISE_ERROR)
974 error->one(FLERR,"Coordinate constraints are not satisfied "
975 "up to desired tolerance ");
976
977 double dv1 = fabs(MathExtra::dot3(r01,v01));
978 double dv2 = fabs(MathExtra::dot3(r02,v02));
979 double dv12 = fabs(MathExtra::dot3(r12,v12));
980
981 if (verr_max < dv1) verr_max = dv1;
982 if (verr_max < dv2) verr_max = dv2;
983 if (verr_max < dv12) verr_max = dv12;
984
985
986 stat = !(checkv && (dv1 > tol ||
987 dv2 > tol ||
988 dv12> tol));
989
990
991 if (!stat && RATTLE_RAISE_ERROR)
992 error->one(FLERR,"Velocity constraints are not satisfied "
993 "up to desired tolerance!");
994
995
996 return stat;
997 }
998