1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 /* ----------------------------------------------------------------------
15 Package FixPIMD
16 Purpose Quantum Path Integral Algorithm for Quantum Chemistry
17 Copyright Voth Group @ University of Chicago
18 Authors Chris Knight & Yuxing Peng (yuxing at uchicago.edu)
19
20 Updated Oct-01-2011
21 Version 1.0
22 ------------------------------------------------------------------------- */
23
24 #include "fix_pimd.h"
25
26 #include "atom.h"
27 #include "comm.h"
28 #include "domain.h"
29 #include "error.h"
30 #include "force.h"
31 #include "math_const.h"
32 #include "memory.h"
33 #include "universe.h"
34 #include "update.h"
35
36 #include <cmath>
37 #include <cstring>
38
39 using namespace LAMMPS_NS;
40 using namespace FixConst;
41 using namespace MathConst;
42
43 enum { PIMD, NMPIMD, CMD };
44
45 /* ---------------------------------------------------------------------- */
46
FixPIMD(LAMMPS * lmp,int narg,char ** arg)47 FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
48 {
49 max_nsend = 0;
50 tag_send = nullptr;
51 buf_send = nullptr;
52
53 max_nlocal = 0;
54 buf_recv = nullptr;
55 buf_beads = nullptr;
56
57 size_plan = 0;
58 plan_send = plan_recv = nullptr;
59
60 M_x2xp = M_xp2x = M_f2fp = M_fp2f = nullptr;
61 lam = nullptr;
62 mode_index = nullptr;
63
64 mass = nullptr;
65
66 array_atom = nullptr;
67 nhc_eta = nullptr;
68 nhc_eta_dot = nullptr;
69 nhc_eta_dotdot = nullptr;
70 nhc_eta_mass = nullptr;
71
72 method = PIMD;
73 fmass = 1.0;
74 nhc_temp = 298.15;
75 nhc_nchain = 2;
76 sp = 1.0;
77
78 for (int i = 3; i < narg - 1; i += 2) {
79 if (strcmp(arg[i], "method") == 0) {
80 if (strcmp(arg[i + 1], "pimd") == 0)
81 method = PIMD;
82 else if (strcmp(arg[i + 1], "nmpimd") == 0)
83 method = NMPIMD;
84 else if (strcmp(arg[i + 1], "cmd") == 0)
85 method = CMD;
86 else
87 error->universe_all(FLERR, "Unknown method parameter for fix pimd");
88 } else if (strcmp(arg[i], "fmass") == 0) {
89 fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
90 if (fmass < 0.0 || fmass > 1.0)
91 error->universe_all(FLERR, "Invalid fmass value for fix pimd");
92 } else if (strcmp(arg[i], "sp") == 0) {
93 sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
94 if (fmass < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd");
95 } else if (strcmp(arg[i], "temp") == 0) {
96 nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
97 if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd");
98 } else if (strcmp(arg[i], "nhc") == 0) {
99 nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp);
100 if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd");
101 } else
102 error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i]));
103 }
104
105 /* Initiation */
106
107 size_peratom_cols = 12 * nhc_nchain + 3;
108
109 nhc_offset_one_1 = 3 * nhc_nchain;
110 nhc_offset_one_2 = 3 * nhc_nchain + 3;
111 nhc_size_one_1 = sizeof(double) * nhc_offset_one_1;
112 nhc_size_one_2 = sizeof(double) * nhc_offset_one_2;
113
114 restart_peratom = 1;
115 peratom_flag = 1;
116 peratom_freq = 1;
117
118 global_freq = 1;
119 vector_flag = 1;
120 size_vector = 2;
121 extvector = 1;
122 comm_forward = 3;
123
124 atom->add_callback(Atom::GROW); // Call LAMMPS to allocate memory for per-atom array
125 atom->add_callback(Atom::RESTART); // Call LAMMPS to re-assign restart-data for per-atom array
126
127 grow_arrays(atom->nmax);
128
129 // some initilizations
130
131 nhc_ready = false;
132 }
133
134 /* ---------------------------------------------------------------------- */
~FixPIMD()135 FixPIMD::~FixPIMD()
136 {
137 delete[] mass;
138 atom->delete_callback(id, Atom::GROW);
139 atom->delete_callback(id, Atom::RESTART);
140
141 memory->destroy(M_x2xp);
142 memory->destroy(M_xp2x);
143 memory->destroy(M_f2fp);
144 memory->destroy(M_fp2f);
145 memory->sfree(lam);
146
147 if (buf_beads)
148 for (int i = 0; i < np; i++) memory->sfree(buf_beads[i]);
149 delete[] buf_beads;
150 delete[] plan_send;
151 delete[] plan_recv;
152 delete[] mode_index;
153
154 memory->sfree(tag_send);
155 memory->sfree(buf_send);
156 memory->sfree(buf_recv);
157
158 memory->destroy(array_atom);
159 memory->destroy(nhc_eta);
160 memory->destroy(nhc_eta_dot);
161 memory->destroy(nhc_eta_dotdot);
162 memory->destroy(nhc_eta_mass);
163 }
164
165 /* ---------------------------------------------------------------------- */
setmask()166 int FixPIMD::setmask()
167 {
168 int mask = 0;
169 mask |= POST_FORCE;
170 mask |= INITIAL_INTEGRATE;
171 mask |= FINAL_INTEGRATE;
172 return mask;
173 }
174
175 /* ---------------------------------------------------------------------- */
176
init()177 void FixPIMD::init()
178 {
179 if (atom->map_style == Atom::MAP_NONE)
180 error->all(FLERR, "Fix pimd requires an atom map, see atom_modify");
181
182 if (universe->me == 0 && universe->uscreen)
183 fprintf(universe->uscreen, "Fix pimd initializing Path-Integral ...\n");
184
185 // prepare the constants
186
187 np = universe->nworlds;
188 inverse_np = 1.0 / np;
189
190 /* The first solution for the force constant, using SI units
191
192 const double Boltzmann = 1.3806488E-23; // SI unit: J/K
193 const double Plank = 6.6260755E-34; // SI unit: m^2 kg / s
194
195 double hbar = Plank / ( 2.0 * MY_PI ) * sp;
196 double beta = 1.0 / ( Boltzmann * input.nh_temp);
197
198 // - P / ( beta^2 * hbar^2) SI unit: s^-2
199 double _fbond = -1.0 / (beta*beta*hbar*hbar) * input.nbeads;
200
201 // convert the units: s^-2 -> (kcal/mol) / (g/mol) / (A^2)
202 fbond = _fbond * 4.184E+26;
203
204 */
205
206 /* The current solution, using LAMMPS internal real units */
207
208 const double Boltzmann = force->boltz;
209 const double Plank = force->hplanck;
210
211 double hbar = Plank / (2.0 * MY_PI);
212 double beta = 1.0 / (Boltzmann * nhc_temp);
213 double _fbond = 1.0 * np / (beta * beta * hbar * hbar);
214
215 omega_np = sqrt(np) / (hbar * beta) * sqrt(force->mvv2e);
216 fbond = -_fbond * force->mvv2e;
217
218 if (universe->me == 0)
219 printf("Fix pimd -P/(beta^2 * hbar^2) = %20.7lE (kcal/mol/A^2)\n\n", fbond);
220
221 dtv = update->dt;
222 dtf = 0.5 * update->dt * force->ftm2v;
223
224 comm_init();
225
226 mass = new double[atom->ntypes + 1];
227
228 if (method == CMD || method == NMPIMD)
229 nmpimd_init();
230 else
231 for (int i = 1; i <= atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass;
232
233 if (!nhc_ready) nhc_init();
234 }
235
236 /* ---------------------------------------------------------------------- */
237
setup(int vflag)238 void FixPIMD::setup(int vflag)
239 {
240 if (universe->me == 0 && universe->uscreen)
241 fprintf(universe->uscreen, "Setting up Path-Integral ...\n");
242
243 post_force(vflag);
244 }
245
246 /* ---------------------------------------------------------------------- */
247
initial_integrate(int)248 void FixPIMD::initial_integrate(int /*vflag*/)
249 {
250 nhc_update_v();
251 nhc_update_x();
252 }
253
254 /* ---------------------------------------------------------------------- */
255
final_integrate()256 void FixPIMD::final_integrate()
257 {
258 nhc_update_v();
259 }
260
261 /* ---------------------------------------------------------------------- */
262
post_force(int)263 void FixPIMD::post_force(int /*flag*/)
264 {
265 for (int i = 0; i < atom->nlocal; i++)
266 for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
267
268 comm_exec(atom->x);
269 spring_force();
270
271 if (method == CMD || method == NMPIMD) {
272 /* forward comm for the force on ghost atoms */
273
274 nmpimd_fill(atom->f);
275
276 /* inter-partition comm */
277
278 comm_exec(atom->f);
279
280 /* normal-mode transform */
281
282 nmpimd_transform(buf_beads, atom->f, M_f2fp[universe->iworld]);
283 }
284 }
285
286 /* ----------------------------------------------------------------------
287 Nose-Hoover Chains
288 ------------------------------------------------------------------------- */
289
nhc_init()290 void FixPIMD::nhc_init()
291 {
292 double tau = 1.0 / omega_np;
293 double KT = force->boltz * nhc_temp;
294
295 double mass0 = KT * tau * tau;
296 int max = 3 * atom->nlocal;
297
298 for (int i = 0; i < max; i++) {
299 for (int ichain = 0; ichain < nhc_nchain; ichain++) {
300 nhc_eta[i][ichain] = 0.0;
301 nhc_eta_dot[i][ichain] = 0.0;
302 nhc_eta_dot[i][ichain] = 0.0;
303 nhc_eta_dotdot[i][ichain] = 0.0;
304 nhc_eta_mass[i][ichain] = mass0;
305 if ((method == CMD || method == NMPIMD) && universe->iworld == 0)
306 ;
307 else
308 nhc_eta_mass[i][ichain] *= fmass;
309 }
310
311 nhc_eta_dot[i][nhc_nchain] = 0.0;
312
313 for (int ichain = 1; ichain < nhc_nchain; ichain++)
314 nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain - 1] * nhc_eta_dot[i][ichain - 1] *
315 nhc_eta_dot[i][ichain - 1] * force->mvv2e -
316 KT) /
317 nhc_eta_mass[i][ichain];
318 }
319
320 // Zero NH acceleration for CMD
321
322 if (method == CMD && universe->iworld == 0)
323 for (int i = 0; i < max; i++)
324 for (int ichain = 0; ichain < nhc_nchain; ichain++) nhc_eta_dotdot[i][ichain] = 0.0;
325
326 nhc_ready = true;
327 }
328
329 /* ---------------------------------------------------------------------- */
330
nhc_update_x()331 void FixPIMD::nhc_update_x()
332 {
333 int n = atom->nlocal;
334 double **x = atom->x;
335 double **v = atom->v;
336
337 if (method == CMD || method == NMPIMD) {
338 nmpimd_fill(atom->v);
339 comm_exec(atom->v);
340
341 /* borrow the space of atom->f to store v in cartisian */
342
343 v = atom->f;
344 nmpimd_transform(buf_beads, v, M_xp2x[universe->iworld]);
345 }
346
347 for (int i = 0; i < n; i++) {
348 x[i][0] += dtv * v[i][0];
349 x[i][1] += dtv * v[i][1];
350 x[i][2] += dtv * v[i][2];
351 }
352 }
353
354 /* ---------------------------------------------------------------------- */
355
nhc_update_v()356 void FixPIMD::nhc_update_v()
357 {
358 int n = atom->nlocal;
359 int *type = atom->type;
360 double **v = atom->v;
361 double **f = atom->f;
362
363 for (int i = 0; i < n; i++) {
364 double dtfm = dtf / mass[type[i]];
365 v[i][0] += dtfm * f[i][0];
366 v[i][1] += dtfm * f[i][1];
367 v[i][2] += dtfm * f[i][2];
368 }
369
370 t_sys = 0.0;
371 if (method == CMD && universe->iworld == 0) return;
372
373 double expfac;
374 int nmax = 3 * atom->nlocal;
375 double KT = force->boltz * nhc_temp;
376 double kecurrent, t_current;
377
378 double dthalf = 0.5 * update->dt;
379 double dt4 = 0.25 * update->dt;
380 double dt8 = 0.125 * update->dt;
381
382 for (int i = 0; i < nmax; i++) {
383 int iatm = i / 3;
384 int idim = i % 3;
385
386 double *vv = v[iatm];
387
388 kecurrent = mass[type[iatm]] * vv[idim] * vv[idim] * force->mvv2e;
389 t_current = kecurrent / force->boltz;
390
391 double *eta = nhc_eta[i];
392 double *eta_dot = nhc_eta_dot[i];
393 double *eta_dotdot = nhc_eta_dotdot[i];
394
395 eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
396
397 for (int ichain = nhc_nchain - 1; ichain > 0; ichain--) {
398 expfac = exp(-dt8 * eta_dot[ichain + 1]);
399 eta_dot[ichain] *= expfac;
400 eta_dot[ichain] += eta_dotdot[ichain] * dt4;
401 eta_dot[ichain] *= expfac;
402 }
403
404 expfac = exp(-dt8 * eta_dot[1]);
405 eta_dot[0] *= expfac;
406 eta_dot[0] += eta_dotdot[0] * dt4;
407 eta_dot[0] *= expfac;
408
409 // Update particle velocities half-step
410
411 double factor_eta = exp(-dthalf * eta_dot[0]);
412 vv[idim] *= factor_eta;
413
414 t_current *= (factor_eta * factor_eta);
415 kecurrent = force->boltz * t_current;
416 eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
417
418 for (int ichain = 0; ichain < nhc_nchain; ichain++) eta[ichain] += dthalf * eta_dot[ichain];
419
420 eta_dot[0] *= expfac;
421 eta_dot[0] += eta_dotdot[0] * dt4;
422 eta_dot[0] *= expfac;
423
424 for (int ichain = 1; ichain < nhc_nchain; ichain++) {
425 expfac = exp(-dt8 * eta_dot[ichain + 1]);
426 eta_dot[ichain] *= expfac;
427 eta_dotdot[ichain] =
428 (nhc_eta_mass[i][ichain - 1] * eta_dot[ichain - 1] * eta_dot[ichain - 1] - KT) /
429 nhc_eta_mass[i][ichain];
430 eta_dot[ichain] += eta_dotdot[ichain] * dt4;
431 eta_dot[ichain] *= expfac;
432 }
433
434 t_sys += t_current;
435 }
436
437 t_sys /= nmax;
438 }
439
440 /* ----------------------------------------------------------------------
441 Normal Mode PIMD
442 ------------------------------------------------------------------------- */
443
nmpimd_init()444 void FixPIMD::nmpimd_init()
445 {
446 memory->create(M_x2xp, np, np, "fix_feynman:M_x2xp");
447 memory->create(M_xp2x, np, np, "fix_feynman:M_xp2x");
448 memory->create(M_f2fp, np, np, "fix_feynman:M_f2fp");
449 memory->create(M_fp2f, np, np, "fix_feynman:M_fp2f");
450
451 lam = (double *) memory->smalloc(sizeof(double) * np, "FixPIMD::lam");
452
453 // Set up eigenvalues
454
455 lam[0] = 0.0;
456 if (np % 2 == 0) lam[np - 1] = 4.0 * np;
457
458 for (int i = 2; i <= np / 2; i++) {
459 lam[2 * i - 3] = lam[2 * i - 2] = 2.0 * np * (1.0 - 1.0 * cos(2.0 * MY_PI * (i - 1) / np));
460 }
461
462 // Set up eigenvectors for non-degenerated modes
463
464 for (int i = 0; i < np; i++) {
465 M_x2xp[0][i] = 1.0 / np;
466 if (np % 2 == 0) M_x2xp[np - 1][i] = 1.0 / np * pow(-1.0, i);
467 }
468
469 // Set up eigenvectors for degenerated modes
470
471 for (int i = 0; i < (np - 1) / 2; i++)
472 for (int j = 0; j < np; j++) {
473 M_x2xp[2 * i + 1][j] = sqrt(2.0) * cos(2.0 * MY_PI * (i + 1) * j / np) / np;
474 M_x2xp[2 * i + 2][j] = -sqrt(2.0) * sin(2.0 * MY_PI * (i + 1) * j / np) / np;
475 }
476
477 // Set up Ut
478
479 for (int i = 0; i < np; i++)
480 for (int j = 0; j < np; j++) {
481 M_xp2x[i][j] = M_x2xp[j][i] * np;
482 M_f2fp[i][j] = M_x2xp[i][j] * np;
483 M_fp2f[i][j] = M_xp2x[i][j];
484 }
485
486 // Set up masses
487
488 int iworld = universe->iworld;
489
490 for (int i = 1; i <= atom->ntypes; i++) {
491 mass[i] = atom->mass[i];
492
493 if (iworld) {
494 mass[i] *= lam[iworld];
495 mass[i] *= fmass;
496 }
497 }
498 }
499
500 /* ---------------------------------------------------------------------- */
501
nmpimd_fill(double ** ptr)502 void FixPIMD::nmpimd_fill(double **ptr)
503 {
504 comm_ptr = ptr;
505 comm->forward_comm_fix(this);
506 }
507
508 /* ---------------------------------------------------------------------- */
509
nmpimd_transform(double ** src,double ** des,double * vector)510 void FixPIMD::nmpimd_transform(double **src, double **des, double *vector)
511 {
512 int n = atom->nlocal;
513 int m = 0;
514
515 for (int i = 0; i < n; i++)
516 for (int d = 0; d < 3; d++) {
517 des[i][d] = 0.0;
518 for (int j = 0; j < np; j++) { des[i][d] += (src[j][m] * vector[j]); }
519 m++;
520 }
521 }
522
523 /* ---------------------------------------------------------------------- */
524
spring_force()525 void FixPIMD::spring_force()
526 {
527 spring_energy = 0.0;
528
529 double **x = atom->x;
530 double **f = atom->f;
531 double *_mass = atom->mass;
532 int *type = atom->type;
533 int nlocal = atom->nlocal;
534
535 double *xlast = buf_beads[x_last];
536 double *xnext = buf_beads[x_next];
537
538 for (int i = 0; i < nlocal; i++) {
539 double delx1 = xlast[0] - x[i][0];
540 double dely1 = xlast[1] - x[i][1];
541 double delz1 = xlast[2] - x[i][2];
542 xlast += 3;
543 domain->minimum_image(delx1, dely1, delz1);
544
545 double delx2 = xnext[0] - x[i][0];
546 double dely2 = xnext[1] - x[i][1];
547 double delz2 = xnext[2] - x[i][2];
548 xnext += 3;
549 domain->minimum_image(delx2, dely2, delz2);
550
551 double ff = fbond * _mass[type[i]];
552
553 double dx = delx1 + delx2;
554 double dy = dely1 + dely2;
555 double dz = delz1 + delz2;
556
557 f[i][0] -= (dx) *ff;
558 f[i][1] -= (dy) *ff;
559 f[i][2] -= (dz) *ff;
560
561 spring_energy += (dx * dx + dy * dy + dz * dz);
562 }
563 }
564
565 /* ----------------------------------------------------------------------
566 Comm operations
567 ------------------------------------------------------------------------- */
568
comm_init()569 void FixPIMD::comm_init()
570 {
571 if (size_plan) {
572 delete[] plan_send;
573 delete[] plan_recv;
574 }
575
576 if (method == PIMD) {
577 size_plan = 2;
578 plan_send = new int[2];
579 plan_recv = new int[2];
580 mode_index = new int[2];
581
582 int rank_last = universe->me - comm->nprocs;
583 int rank_next = universe->me + comm->nprocs;
584 if (rank_last < 0) rank_last += universe->nprocs;
585 if (rank_next >= universe->nprocs) rank_next -= universe->nprocs;
586
587 plan_send[0] = rank_next;
588 plan_send[1] = rank_last;
589 plan_recv[0] = rank_last;
590 plan_recv[1] = rank_next;
591
592 mode_index[0] = 0;
593 mode_index[1] = 1;
594 x_last = 1;
595 x_next = 0;
596 } else {
597 size_plan = np - 1;
598 plan_send = new int[size_plan];
599 plan_recv = new int[size_plan];
600 mode_index = new int[size_plan];
601
602 for (int i = 0; i < size_plan; i++) {
603 plan_send[i] = universe->me + comm->nprocs * (i + 1);
604 if (plan_send[i] >= universe->nprocs) plan_send[i] -= universe->nprocs;
605
606 plan_recv[i] = universe->me - comm->nprocs * (i + 1);
607 if (plan_recv[i] < 0) plan_recv[i] += universe->nprocs;
608
609 mode_index[i] = (universe->iworld + i + 1) % (universe->nworlds);
610 }
611
612 x_next = (universe->iworld + 1 + universe->nworlds) % (universe->nworlds);
613 x_last = (universe->iworld - 1 + universe->nworlds) % (universe->nworlds);
614 }
615
616 if (buf_beads) {
617 for (int i = 0; i < np; i++)
618 if (buf_beads[i]) delete[] buf_beads[i];
619 delete[] buf_beads;
620 }
621
622 buf_beads = new double *[np];
623 for (int i = 0; i < np; i++) buf_beads[i] = nullptr;
624 }
625
626 /* ---------------------------------------------------------------------- */
627
comm_exec(double ** ptr)628 void FixPIMD::comm_exec(double **ptr)
629 {
630 int nlocal = atom->nlocal;
631
632 if (nlocal > max_nlocal) {
633 max_nlocal = nlocal + 200;
634 int size = sizeof(double) * max_nlocal * 3;
635 buf_recv = (double *) memory->srealloc(buf_recv, size, "FixPIMD:x_recv");
636
637 for (int i = 0; i < np; i++)
638 buf_beads[i] = (double *) memory->srealloc(buf_beads[i], size, "FixPIMD:x_beads[i]");
639 }
640
641 // copy local positions
642
643 memcpy(buf_beads[universe->iworld], &(ptr[0][0]), sizeof(double) * nlocal * 3);
644
645 // go over comm plans
646
647 for (int iplan = 0; iplan < size_plan; iplan++) {
648 // sendrecv nlocal
649
650 int nsend;
651
652 MPI_Sendrecv(&(nlocal), 1, MPI_INT, plan_send[iplan], 0, &(nsend), 1, MPI_INT, plan_recv[iplan],
653 0, universe->uworld, MPI_STATUS_IGNORE);
654
655 // allocate arrays
656
657 if (nsend > max_nsend) {
658 max_nsend = nsend + 200;
659 tag_send =
660 (tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMD:tag_send");
661 buf_send =
662 (double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3, "FixPIMD:x_send");
663 }
664
665 // send tags
666
667 MPI_Sendrecv(atom->tag, nlocal, MPI_LMP_TAGINT, plan_send[iplan], 0, tag_send, nsend,
668 MPI_LMP_TAGINT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
669
670 // wrap positions
671
672 double *wrap_ptr = buf_send;
673 int ncpy = sizeof(double) * 3;
674
675 for (int i = 0; i < nsend; i++) {
676 int index = atom->map(tag_send[i]);
677
678 if (index < 0) {
679 char error_line[256];
680
681 sprintf(error_line,
682 "Atom " TAGINT_FORMAT " is missing at world [%d] "
683 "rank [%d] required by rank [%d] (" TAGINT_FORMAT ", " TAGINT_FORMAT
684 ", " TAGINT_FORMAT ").\n",
685 tag_send[i], universe->iworld, comm->me, plan_recv[iplan], atom->tag[0],
686 atom->tag[1], atom->tag[2]);
687
688 error->universe_one(FLERR, error_line);
689 }
690
691 memcpy(wrap_ptr, ptr[index], ncpy);
692 wrap_ptr += 3;
693 }
694
695 // sendrecv x
696
697 MPI_Sendrecv(buf_send, nsend * 3, MPI_DOUBLE, plan_recv[iplan], 0, buf_recv, nlocal * 3,
698 MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
699
700 // copy x
701
702 memcpy(buf_beads[mode_index[iplan]], buf_recv, sizeof(double) * nlocal * 3);
703 }
704 }
705
706 /* ---------------------------------------------------------------------- */
707
pack_forward_comm(int n,int * list,double * buf,int,int *)708 int FixPIMD::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
709 {
710 int i, j, m;
711
712 m = 0;
713
714 for (i = 0; i < n; i++) {
715 j = list[i];
716 buf[m++] = comm_ptr[j][0];
717 buf[m++] = comm_ptr[j][1];
718 buf[m++] = comm_ptr[j][2];
719 }
720
721 return m;
722 }
723
724 /* ---------------------------------------------------------------------- */
725
unpack_forward_comm(int n,int first,double * buf)726 void FixPIMD::unpack_forward_comm(int n, int first, double *buf)
727 {
728 int i, m, last;
729
730 m = 0;
731 last = first + n;
732 for (i = first; i < last; i++) {
733 comm_ptr[i][0] = buf[m++];
734 comm_ptr[i][1] = buf[m++];
735 comm_ptr[i][2] = buf[m++];
736 }
737 }
738
739 /* ----------------------------------------------------------------------
740 Memory operations
741 ------------------------------------------------------------------------- */
742
memory_usage()743 double FixPIMD::memory_usage()
744 {
745 return (double) atom->nmax * size_peratom_cols * sizeof(double);
746 }
747
748 /* ---------------------------------------------------------------------- */
749
grow_arrays(int nmax)750 void FixPIMD::grow_arrays(int nmax)
751 {
752 if (nmax == 0) return;
753 int count = nmax * 3;
754
755 memory->grow(array_atom, nmax, size_peratom_cols, "FixPIMD::array_atom");
756 memory->grow(nhc_eta, count, nhc_nchain, "FixPIMD::nh_eta");
757 memory->grow(nhc_eta_dot, count, nhc_nchain + 1, "FixPIMD::nh_eta_dot");
758 memory->grow(nhc_eta_dotdot, count, nhc_nchain, "FixPIMD::nh_eta_dotdot");
759 memory->grow(nhc_eta_mass, count, nhc_nchain, "FixPIMD::nh_eta_mass");
760 }
761
762 /* ---------------------------------------------------------------------- */
763
copy_arrays(int i,int j,int)764 void FixPIMD::copy_arrays(int i, int j, int /*delflag*/)
765 {
766 int i_pos = i * 3;
767 int j_pos = j * 3;
768
769 memcpy(nhc_eta[j_pos], nhc_eta[i_pos], nhc_size_one_1);
770 memcpy(nhc_eta_dot[j_pos], nhc_eta_dot[i_pos], nhc_size_one_2);
771 memcpy(nhc_eta_dotdot[j_pos], nhc_eta_dotdot[i_pos], nhc_size_one_1);
772 memcpy(nhc_eta_mass[j_pos], nhc_eta_mass[i_pos], nhc_size_one_1);
773 }
774
775 /* ---------------------------------------------------------------------- */
776
pack_exchange(int i,double * buf)777 int FixPIMD::pack_exchange(int i, double *buf)
778 {
779 int offset = 0;
780 int pos = i * 3;
781
782 memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1);
783 offset += nhc_offset_one_1;
784 memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2);
785 offset += nhc_offset_one_2;
786 memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1);
787 offset += nhc_offset_one_1;
788 memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1);
789 offset += nhc_offset_one_1;
790
791 return size_peratom_cols;
792 }
793
794 /* ---------------------------------------------------------------------- */
795
unpack_exchange(int nlocal,double * buf)796 int FixPIMD::unpack_exchange(int nlocal, double *buf)
797 {
798 int offset = 0;
799 int pos = nlocal * 3;
800
801 memcpy(nhc_eta[pos], buf + offset, nhc_size_one_1);
802 offset += nhc_offset_one_1;
803 memcpy(nhc_eta_dot[pos], buf + offset, nhc_size_one_2);
804 offset += nhc_offset_one_2;
805 memcpy(nhc_eta_dotdot[pos], buf + offset, nhc_size_one_1);
806 offset += nhc_offset_one_1;
807 memcpy(nhc_eta_mass[pos], buf + offset, nhc_size_one_1);
808 offset += nhc_offset_one_1;
809
810 return size_peratom_cols;
811 }
812
813 /* ---------------------------------------------------------------------- */
814
pack_restart(int i,double * buf)815 int FixPIMD::pack_restart(int i, double *buf)
816 {
817 int offset = 0;
818 int pos = i * 3;
819 // pack buf[0] this way because other fixes unpack it
820 buf[offset++] = size_peratom_cols + 1;
821
822 memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1);
823 offset += nhc_offset_one_1;
824 memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2);
825 offset += nhc_offset_one_2;
826 memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1);
827 offset += nhc_offset_one_1;
828 memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1);
829 offset += nhc_offset_one_1;
830
831 return size_peratom_cols + 1;
832 }
833
834 /* ---------------------------------------------------------------------- */
835
unpack_restart(int nlocal,int nth)836 void FixPIMD::unpack_restart(int nlocal, int nth)
837 {
838 double **extra = atom->extra;
839
840 // skip to Nth set of extra values
841 // unpack the Nth first values this way because other fixes pack them
842
843 int m = 0;
844 for (int i = 0; i < nth; i++) m += static_cast<int>(extra[nlocal][m]);
845 m++;
846
847 int pos = nlocal * 3;
848
849 memcpy(nhc_eta[pos], extra[nlocal] + m, nhc_size_one_1);
850 m += nhc_offset_one_1;
851 memcpy(nhc_eta_dot[pos], extra[nlocal] + m, nhc_size_one_2);
852 m += nhc_offset_one_2;
853 memcpy(nhc_eta_dotdot[pos], extra[nlocal] + m, nhc_size_one_1);
854 m += nhc_offset_one_1;
855 memcpy(nhc_eta_mass[pos], extra[nlocal] + m, nhc_size_one_1);
856 m += nhc_offset_one_1;
857
858 nhc_ready = true;
859 }
860
861 /* ---------------------------------------------------------------------- */
862
maxsize_restart()863 int FixPIMD::maxsize_restart()
864 {
865 return size_peratom_cols + 1;
866 }
867
868 /* ---------------------------------------------------------------------- */
869
size_restart(int)870 int FixPIMD::size_restart(int /*nlocal*/)
871 {
872 return size_peratom_cols + 1;
873 }
874
875 /* ---------------------------------------------------------------------- */
876
compute_vector(int n)877 double FixPIMD::compute_vector(int n)
878 {
879 if (n == 0) { return spring_energy; }
880 if (n == 1) { return t_sys; }
881 return 0.0;
882 }
883