From e6d31485335cb3f239d3d722d0bddd151bbd4cc3 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Thu, 9 Nov 2023 23:34:01 -0500 Subject: [PATCH 01/11] add method=pimd support for Langevin thermostat --- src/REPLICA/fix_pimd_langevin.cpp | 120 +++++++++++++++++++----------- 1 file changed, 76 insertions(+), 44 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index cffaf327e4..d4e2b6e526 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -50,7 +50,7 @@ using namespace FixConst; using MathConst::MY_PI; using MathConst::THIRD; -enum { NMPIMD }; +enum { PIMD, NMPIMD }; enum { PHYSICAL, NORMAL }; enum { BAOAB, OBABO }; enum { ISO, ANISO, TRICLINIC }; @@ -121,6 +121,8 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[i], "method") == 0) { if (strcmp(arg[i + 1], "nmpimd") == 0) method = NMPIMD; + elif (strcmp(arg[i + 1], "pimd") == 0) + method = PIMD; else error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin"); } else if (strcmp(arg[i], "integrator") == 0) { @@ -159,7 +161,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : error->universe_all(FLERR, "Invalid fmass value for fix pimd/langevin"); } else if (strcmp(arg[i], "sp") == 0) { sp = utils::numeric(FLERR, arg[i + 1], false, lmp); - if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/nvt"); + if (sp < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd/langevin"); } else if (strcmp(arg[i], "fmmode") == 0) { if (strcmp(arg[i + 1], "physical") == 0) fmmode = PHYSICAL; @@ -170,9 +172,11 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : "Unknown fictitious mass mode for fix pimd/langevin. Only physical " "mass and normal mode mass are supported!"); } else if (strcmp(arg[i], "scale") == 0) { + if (method == PIMD) + error->universe_all(FLERR, "The scale parameter of the PILE_L thermostat is not supported for method pimd. Delete scale parameter if you do want to use method pimd."); pilescale = utils::numeric(FLERR, arg[i + 1], false, lmp); if (pilescale < 0.0) - error->universe_all(FLERR, "Invalid pile scale value for fix pimd/langevin"); + error->universe_all(FLERR, "Invalid PILE_L scale value for fix pimd/langevin"); } else if (strcmp(arg[i], "temp") == 0) { temp = utils::numeric(FLERR, arg[i + 1], false, lmp); if (temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd/langevin"); @@ -245,6 +249,12 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : error->universe_all( FLERR, fmt::format("Must not use pressure coupling with {} ensemble", Ensembles[ensemble])); + if (method == PIMD && pstat_flag) + error->universe_all(FLERR, "Pressure control has not been supported for method pimd yet. Please set method to nmpimd."); + + if (method == PIMD && fmmode == NORMAL) + error->universe_all(FLERR, "Normal mode mass is not supported for method pimd. Please set method to nmpimd."); + /* Initiation */ global_freq = 1; @@ -890,20 +900,22 @@ void FixPIMDLangevin::langevin_init() _omega_k = new double[np]; Lan_c = new double[np]; Lan_s = new double[np]; - if (fmmode == PHYSICAL) { - for (int i = 0; i < np; i++) { - _omega_k[i] = _omega_np * sqrt(lam[i]) / sqrt(fmass); - Lan_c[i] = cos(sqrt(lam[i]) * _omega_np_dt_half); - Lan_s[i] = sin(sqrt(lam[i]) * _omega_np_dt_half); + if (method == NMPIMD) { + if (fmmode == PHYSICAL) { + for (int i = 0; i < np; i++) { + _omega_k[i] = _omega_np * sqrt(lam[i]) / sqrt(fmass); + Lan_c[i] = cos(sqrt(lam[i]) * _omega_np_dt_half); + Lan_s[i] = sin(sqrt(lam[i]) * _omega_np_dt_half); + } + } else if (fmmode == NORMAL) { + for (int i = 0; i < np; i++) { + _omega_k[i] = _omega_np / sqrt(fmass); + Lan_c[i] = cos(_omega_np_dt_half); + Lan_s[i] = sin(_omega_np_dt_half); + } + } else { + error->universe_all(FLERR, "Unknown fmmode setting; only physical and normal are supported!"); } - } else if (fmmode == NORMAL) { - for (int i = 0; i < np; i++) { - _omega_k[i] = _omega_np / sqrt(fmass); - Lan_c[i] = cos(_omega_np_dt_half); - Lan_s[i] = sin(_omega_np_dt_half); - } - } else { - error->universe_all(FLERR, "Unknown fmmode setting; only physical and normal are supported!"); } if (tau > 0) @@ -925,27 +937,35 @@ void FixPIMDLangevin::langevin_init() if (thermostat == PILE_L) { std::string out = "\nInitializing PI Langevin equation thermostat...\n"; out += "Bead ID | omega | tau | c1 | c2\n"; - tau_k = new double[np]; - c1_k = new double[np]; - c2_k = new double[np]; - tau_k[0] = tau; - c1_k[0] = c1; - c2_k[0] = c2; - for (int i = 1; i < np; i++) { - tau_k[i] = 0.5 / pilescale / _omega_k[i]; - if (integrator == OBABO) - c1_k[i] = exp(-0.5 * update->dt / tau_k[i]); - else if (integrator == BAOAB) - c1_k[i] = exp(-1.0 * update->dt / tau_k[i]); - else - error->universe_all(FLERR, - "Unknown integrator parameter for fix pimd/langevin. Only obabo and " - "baoab integrators are supported!"); - c2_k[i] = sqrt(1.0 - c1_k[i] * c1_k[i]); - } - for (int i = 0; i < np; i++) { - out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_k[i], tau_k[i], - c1_k[i], c2_k[i]); + if (method == NMPIMD) { + tau_k = new double[np]; + c1_k = new double[np]; + c2_k = new double[np]; + tau_k[0] = tau; + c1_k[0] = c1; + c2_k[0] = c2; + for (int i = 1; i < np; i++) { + tau_k[i] = 0.5 / pilescale / _omega_k[i]; + if (integrator == OBABO) + c1_k[i] = exp(-0.5 * update->dt / tau_k[i]); + else if (integrator == BAOAB) + c1_k[i] = exp(-1.0 * update->dt / tau_k[i]); + else + error->universe_all(FLERR, + "Unknown integrator parameter for fix pimd/langevin. Only obabo and " + "baoab integrators are supported!"); + c2_k[i] = sqrt(1.0 - c1_k[i] * c1_k[i]); + } + for (int i = 0; i < np; i++) { + out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_k[i], tau_k[i], + c1_k[i], c2_k[i]); + } + } + else if (method == PIMD) { + for (int i = 0; i < np; i++) { + out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_np / sqrt(fmass), tau, + c1, c2); + } } if (thermostat == PILE_L) out += "PILE_L thermostat successfully initialized!\n"; out += "\n"; @@ -961,13 +981,25 @@ void FixPIMDLangevin::o_step() int *type = atom->type; double beta_np = 1.0 / force->boltz / Lan_temp * inverse_np * force->mvv2e; if (thermostat == PILE_L) { - for (int i = 0; i < nlocal; i++) { - atom->v[i][0] = c1_k[universe->iworld] * atom->v[i][0] + - c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); - atom->v[i][1] = c1_k[universe->iworld] * atom->v[i][1] + - c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); - atom->v[i][2] = c1_k[universe->iworld] * atom->v[i][2] + - c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + if (method == NMPIMD) { + for (int i = 0; i < nlocal; i++) { + atom->v[i][0] = c1_k[universe->iworld] * atom->v[i][0] + + c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][1] = c1_k[universe->iworld] * atom->v[i][1] + + c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][2] = c1_k[universe->iworld] * atom->v[i][2] + + c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + } + } + else if (method == PIMD) { + for (int i = 0; i < nlocal; i++) { + atom->v[i][0] = c1 * atom->v[i][0] + + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][1] = c1 * atom->v[i][1] + + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][2] = c1 * atom->v[i][2] + + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + } } } } From f413d395a52c54b9e0e47e0adcd0c3a297aaba75 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Fri, 10 Nov 2023 00:20:51 -0500 Subject: [PATCH 02/11] add support for method=PIMD; add q_step function --- src/REPLICA/fix_pimd_langevin.cpp | 74 ++++++++++++++++++++++++++----- src/REPLICA/fix_pimd_langevin.h | 4 ++ 2 files changed, 68 insertions(+), 10 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index d4e2b6e526..0068546f44 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -538,11 +538,18 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]); else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); + qc_step(); + a_step(); + qc_step(); + a_step(); + } + else if (method == PIMD) { + q_step(); + q_step(); + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } - qc_step(); - a_step(); - qc_step(); - a_step(); } else if (integrator == BAOAB) { if (pstat_flag) { compute_totke(); @@ -556,18 +563,32 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]); else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); + qc_step(); + a_step(); + } + else if (method == PIMD) { + q_step(); + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } - qc_step(); - a_step(); if (tstat_flag) { o_step(); if (removecomflag) remove_com_motion(); if (pstat_flag) press_o_step(); } - qc_step(); - a_step(); + if (method == NMPIMD) { + qc_step(); + a_step(); + } + else if (method == PIMD) { + q_step(); + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } } else { - error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin"); + error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab integrators are supported!"); } collect_xc(); compute_spring_energy(); @@ -614,6 +635,7 @@ void FixPIMDLangevin::final_integrate() void FixPIMDLangevin::post_force(int /*flag*/) { + if (method == NMPIMD) { if (atom->nmax > maxunwrap) reallocate_x_unwrap(); if (atom->nmax > maxxc) reallocate_xc(); int nlocal = atom->nlocal; @@ -638,6 +660,7 @@ void FixPIMDLangevin::post_force(int /*flag*/) compute_vir(); compute_cvir(); compute_t_vir(); + } compute_pote(); if (method == NMPIMD) { inter_replica_comm(f); @@ -646,6 +669,12 @@ void FixPIMDLangevin::post_force(int /*flag*/) else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, f, M_x2xp[universe->iworld]); } + else if (method == PIMD) { + spring_force(); + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } c_pe->addstep(update->ntimestep + 1); c_press->addstep(update->ntimestep + 1); } @@ -660,6 +689,8 @@ void FixPIMDLangevin::end_of_step() if (pstat_flag) compute_totenthalpy(); } +/* ---------------------------------------------------------------------- */ + void FixPIMDLangevin::collect_xc() { int nlocal = atom->nlocal; @@ -693,7 +724,9 @@ void FixPIMDLangevin::collect_xc() void FixPIMDLangevin::b_step() { - + // used for both NMPIMD and PIMD + // For NMPIMD, force only includes the contribution of external potential. + // For PIMD, force includes the contributions of external potential and spring force. int n = atom->nlocal; int *type = atom->type; double **v = atom->v; @@ -711,6 +744,8 @@ void FixPIMDLangevin::b_step() void FixPIMDLangevin::qc_step() { + // used for NMPIMD + // evolve the centroid mode int nlocal = atom->nlocal; double **x = atom->x; double **v = atom->v; @@ -775,6 +810,8 @@ void FixPIMDLangevin::qc_step() void FixPIMDLangevin::a_step() { + // used for NMPIMD + // use analytical solution of harmonic oscillator to evolve the non-centroid modes int n = atom->nlocal; double **x = atom->x; double **v = atom->v; @@ -806,6 +843,14 @@ void FixPIMDLangevin::a_step() /* ---------------------------------------------------------------------- */ +void FixPIMDLangevin::q_step() +{ + // used for PIMD + // evolve all beads +} + +/* ---------------------------------------------------------------------- */ + void FixPIMDLangevin::baro_init() { vw[0] = vw[1] = vw[2] = vw[3] = vw[4] = vw[5] = 0.0; @@ -1086,6 +1131,12 @@ void FixPIMDLangevin::nmpimd_transform(double **src, double **des, double *vecto } } +/* ---------------------------------------------------------------------- */ + +void FixPIMDLangevin::spring_force() +{ +} + /* ---------------------------------------------------------------------- Comm operations ------------------------------------------------------------------------- */ @@ -1111,6 +1162,9 @@ void FixPIMDLangevin::comm_init() planrecv[i] = universe->root_proc[irecv]; modeindex[i] = irecv; } + + x_next = (universe->iworld + 1 + universe->nworlds) % (universe->nworlds); + x_last = (universe->iworld - 1 + universe->nworlds) % (universe->nworlds); } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/fix_pimd_langevin.h b/src/REPLICA/fix_pimd_langevin.h index 9730f65376..0f21b908b0 100644 --- a/src/REPLICA/fix_pimd_langevin.h +++ b/src/REPLICA/fix_pimd_langevin.h @@ -77,6 +77,8 @@ class FixPIMDLangevin : public Fix { int me, nprocs, ireplica, nreplica, nprocs_universe; int ntotal, maxlocal; + int x_last, x_next; + int cmode; int sizeplan; int *plansend, *planrecv; @@ -93,6 +95,7 @@ class FixPIMDLangevin : public Fix { void comm_init(); void inter_replica_comm(double **ptr); + void spring_force(); /* normal-mode operations */ @@ -121,6 +124,7 @@ class FixPIMDLangevin : public Fix { a_step(); // integrate for dt/2 according to A part (non-centroid mode, harmonic force between replicas) void qc_step(); // integrate for dt/2 for the centroid mode (x <- x + v * dt/2) void o_step(); // integrate for dt according to O part (O-U process, for thermostating) + void q_step(); // integrate for dt/2 for all the beads (x <- x + v * dt/2) /* Bussi-Zykova-Parrinello barostat */ From 63935194191f8b5d06005647d8c28914cc475b73 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Fri, 10 Nov 2023 16:07:00 -0500 Subject: [PATCH 03/11] q_step function --- src/REPLICA/fix_pimd_langevin.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 0068546f44..5a2df1bc75 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -847,6 +847,17 @@ void FixPIMDLangevin::q_step() { // used for PIMD // evolve all beads + int nlocal = atom->nlocal; + double **x = atom->x; + double **v = atom->v; + + if (!pstat_flag) { + for (int i = 0; i < nlocal; i++) { + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } } /* ---------------------------------------------------------------------- */ From 98a0f43c9bbeac7f1688b073710d6b8d4e3aea09 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Sat, 11 Nov 2023 01:04:13 -0500 Subject: [PATCH 04/11] add spring_force function --- src/REPLICA/fix_pimd_langevin.cpp | 82 ++++++++++++++++++++++++++----- 1 file changed, 69 insertions(+), 13 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 5a2df1bc75..6ef74c2a1a 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -121,7 +121,7 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : if (strcmp(arg[i], "method") == 0) { if (strcmp(arg[i + 1], "nmpimd") == 0) method = NMPIMD; - elif (strcmp(arg[i + 1], "pimd") == 0) + else if (strcmp(arg[i + 1], "pimd") == 0) method = PIMD; else error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin"); @@ -488,6 +488,13 @@ void FixPIMDLangevin::setup(int vflag) else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); } + else if (method == PIMD) { + inter_replica_comm(x); + spring_force(); + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } collect_xc(); compute_spring_energy(); compute_t_prim(); @@ -591,9 +598,16 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab integrators are supported!"); } collect_xc(); - compute_spring_energy(); - compute_t_prim(); - compute_p_prim(); + if (method == PIMD) { + inter_replica_comm(x); + spring_force(); + compute_spring_energy(); + } + if (method == NMPIMD) { + compute_spring_energy(); + compute_t_prim(); + compute_p_prim(); + } if (method == NMPIMD) { inter_replica_comm(x); @@ -635,14 +649,15 @@ void FixPIMDLangevin::final_integrate() void FixPIMDLangevin::post_force(int /*flag*/) { - if (method == NMPIMD) { - if (atom->nmax > maxunwrap) reallocate_x_unwrap(); - if (atom->nmax > maxxc) reallocate_xc(); int nlocal = atom->nlocal; double **x = atom->x; double **f = atom->f; imageint *image = atom->image; tagint *tag = atom->tag; + + if (method == NMPIMD) { + if (atom->nmax > maxunwrap) reallocate_x_unwrap(); + if (atom->nmax > maxxc) reallocate_xc(); for (int i = 0; i < nlocal; i++) { x_unwrap[i][0] = x[i][0]; x_unwrap[i][1] = x[i][1]; @@ -669,12 +684,7 @@ void FixPIMDLangevin::post_force(int /*flag*/) else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, f, M_x2xp[universe->iworld]); } - else if (method == PIMD) { - spring_force(); - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); - } + c_pe->addstep(update->ntimestep + 1); c_press->addstep(update->ntimestep + 1); } @@ -1146,6 +1156,41 @@ void FixPIMDLangevin::nmpimd_transform(double **src, double **des, double *vecto void FixPIMDLangevin::spring_force() { + spring_energy = 0.0; + + double **x = atom->x; + double **f = atom->f; + double* _mass = atom->mass; + int* type = atom->type; + int nlocal = atom->nlocal; + tagint* tagtmp = atom->tag; + + int *mask = atom->mask; + + for (int i=0; iuworld); total_spring_energy /= universe->procs_per_world[universe->iworld]; + } + else if (method == PIMD) { + total_spring_energy = se_bead = 0.0; + MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld); + total_spring_energy /= universe->procs_per_world[universe->iworld]; + } + else { + error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } } /* ---------------------------------------------------------------------- */ From ba32afc06e293384653f4856069575268ae47e9c Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Sat, 11 Nov 2023 01:45:36 -0500 Subject: [PATCH 05/11] fix spring_energy --- src/REPLICA/fix_pimd_langevin.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 6ef74c2a1a..96cf4bef35 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -1188,7 +1188,7 @@ void FixPIMDLangevin::spring_force() f[i][1] -= (dy) * ff; f[i][2] -= (dz) * ff; - spring_energy += ff * (delx2*delx2+dely2*dely2+delz2*delz2); + spring_energy += 0.5 * ff * (delx2*delx2+dely2*dely2+delz2*delz2); } } } From 4ef27552c4492b92c3ea889309ccd2ec3e2c8245 Mon Sep 17 00:00:00 2001 From: Yifan Li Date: Sat, 11 Nov 2023 03:10:00 -0500 Subject: [PATCH 06/11] fix spring_force()'s position --- src/REPLICA/fix_pimd_langevin.cpp | 28 ++++++++++++++++++++-------- 1 file changed, 20 insertions(+), 8 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 96cf4bef35..4e94b77a86 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -598,11 +598,7 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab integrators are supported!"); } collect_xc(); - if (method == PIMD) { - inter_replica_comm(x); - spring_force(); - compute_spring_energy(); - } + if (method == NMPIMD) { compute_spring_energy(); compute_t_prim(); @@ -676,6 +672,18 @@ void FixPIMDLangevin::post_force(int /*flag*/) compute_cvir(); compute_t_vir(); } + + if (method == PIMD) { + if (mapflag) { + for (int i = 0; i < nlocal; i++) { domain->unmap(x[i], image[i]); } + } + inter_replica_comm(x); + spring_force(); + compute_spring_energy(); + if (mapflag) { + for (int i = 0; i < nlocal; i++) { domain->unmap_inv(x[i], image[i]); } + } + } compute_pote(); if (method == NMPIMD) { inter_replica_comm(f); @@ -1165,7 +1173,10 @@ void FixPIMDLangevin::spring_force() int nlocal = atom->nlocal; tagint* tagtmp = atom->tag; + // printf("iworld = %d, x_last = %d, x_next = %d\n", universe->iworld, x_last, x_next); int *mask = atom->mask; + + // int idx_tmp = atom->map(1); for (int i=0; i Date: Sat, 11 Nov 2023 04:10:31 -0500 Subject: [PATCH 07/11] update document for method=pimd --- doc/src/fix_pimd.rst | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/doc/src/fix_pimd.rst b/doc/src/fix_pimd.rst index 5b51b97c52..6abac408ca 100644 --- a/doc/src/fix_pimd.rst +++ b/doc/src/fix_pimd.rst @@ -31,7 +31,7 @@ Syntax .. parsed-literal:: *keywords* = *method* or *integrator* or *ensemble* or *fmmode* or *fmass* or *scale* or *temp* or *thermostat* or *tau* or *iso* or *aniso* or *barostat* or *taup* or *fixcom* or *lj* - *method* value = *nmpimd* + *method* value = *nmpimd* (default) or *pimd* *integrator* value = *obabo* or *baoab* *fmmode* value = *physical* or *normal* *fmass* value = scaling factor on mass @@ -137,9 +137,6 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics the real particle. .. note:: - Fix pimd/langevin only supports *method* value *nmpimd*. This should be enough - for most PIMD applications for quantum thermodynamics purpose. - Motion of the centroid can be effectively uncoupled from the other normal modes by scaling the fictitious masses to achieve a partial adiabatic separation. This is called a Centroid Molecular Dynamics @@ -151,6 +148,10 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics only the k > 0 modes are thermostatted, not the centroid degrees of freedom. +Fix pimd/langevin supports *method* value *nmpimd* and *pimd*. The default value is *nmpimd*. +If *method* is *nmpimd*, the normal mode representation is used to integrate the equations of motion. The exact solution of harmonic oscillator is used to propagate the free ring polymer part of the Hamiltonian. +If *method* is *pimd*, the Cartesian representation is used to integrate the equations of motion. The harmonic force is added to the total force of the system, and the numerical integrator is used to propagate the Hamiltonian. + The keyword *integrator* specifies the Trotter splitting method used by *fix pimd/langevin*. See :ref:`(Liu) ` for a discussion on the OBABO and BAOAB splitting schemes. Typically either of the two should work fine. From 3d14e2e0e23b1b7ee4efb74b71613c129910260e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Nov 2023 06:05:07 -0500 Subject: [PATCH 08/11] whitespace --- src/REPLICA/fix_pimd_langevin.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index 4e94b77a86..ba02b7a184 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -1177,25 +1177,25 @@ void FixPIMDLangevin::spring_force() int *mask = atom->mask; // int idx_tmp = atom->map(1); - + for (int i=0; i Date: Sat, 11 Nov 2023 06:09:15 -0500 Subject: [PATCH 09/11] doc tweaks --- doc/src/fix_pimd.rst | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_pimd.rst b/doc/src/fix_pimd.rst index 6abac408ca..91c5e58add 100644 --- a/doc/src/fix_pimd.rst +++ b/doc/src/fix_pimd.rst @@ -137,6 +137,7 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics the real particle. .. note:: + Motion of the centroid can be effectively uncoupled from the other normal modes by scaling the fictitious masses to achieve a partial adiabatic separation. This is called a Centroid Molecular Dynamics @@ -148,9 +149,15 @@ normal-mode PIMD. A value of *cmd* is for centroid molecular dynamics only the k > 0 modes are thermostatted, not the centroid degrees of freedom. -Fix pimd/langevin supports *method* value *nmpimd* and *pimd*. The default value is *nmpimd*. -If *method* is *nmpimd*, the normal mode representation is used to integrate the equations of motion. The exact solution of harmonic oscillator is used to propagate the free ring polymer part of the Hamiltonian. -If *method* is *pimd*, the Cartesian representation is used to integrate the equations of motion. The harmonic force is added to the total force of the system, and the numerical integrator is used to propagate the Hamiltonian. +.. versionadded:: TBD + + Mode *pimd* added to fix pimd/langevin. + +Fix pimd/langevin supports the *method* values *nmpimd* and *pimd*. The default value is *nmpimd*. +If *method* is *nmpimd*, the normal mode representation is used to integrate the equations of motion. +The exact solution of harmonic oscillator is used to propagate the free ring polymer part of the Hamiltonian. +If *method* is *pimd*, the Cartesian representation is used to integrate the equations of motion. +The harmonic force is added to the total force of the system, and the numerical integrator is used to propagate the Hamiltonian. The keyword *integrator* specifies the Trotter splitting method used by *fix pimd/langevin*. See :ref:`(Liu) ` for a discussion on the OBABO and BAOAB splitting schemes. Typically @@ -208,6 +215,7 @@ The keyword *thermostat* reads *style* and *seed* of thermostat for fix style *p be *PILE_L* (path integral Langevin equation local thermostat, as described in :ref:`Ceriotti `), and *seed* should a positive integer number, which serves as the seed of the pseudo random number generator. .. note:: + The fix style *pimd/langevin* uses the stochastic PILE_L thermostat to control temperature. This thermostat works on the normal modes of the ring polymer. The *tau* parameter controls the centroid mode, and the *scale* parameter controls the non-centroid modes. @@ -270,6 +278,7 @@ related tasks for each of the partitions, e.g. read_restart system_${ibead}.restart2 .. note:: + Fix *pimd/langevin* dumps the Cartesian coordinates, but dumps the velocities and forces in the normal mode representation. If the Cartesian velocities and forces are needed, it is easy to perform the transformation when doing post-processing. From 9ef1b2d64d7a97af80cb20462026533430781dfc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 11 Nov 2023 06:20:55 -0500 Subject: [PATCH 10/11] apply clang-format --- src/REPLICA/fix_pimd_langevin.cpp | 195 +++++++++++++++--------------- 1 file changed, 100 insertions(+), 95 deletions(-) diff --git a/src/REPLICA/fix_pimd_langevin.cpp b/src/REPLICA/fix_pimd_langevin.cpp index ba02b7a184..d328420ce9 100644 --- a/src/REPLICA/fix_pimd_langevin.cpp +++ b/src/REPLICA/fix_pimd_langevin.cpp @@ -173,7 +173,10 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : "mass and normal mode mass are supported!"); } else if (strcmp(arg[i], "scale") == 0) { if (method == PIMD) - error->universe_all(FLERR, "The scale parameter of the PILE_L thermostat is not supported for method pimd. Delete scale parameter if you do want to use method pimd."); + error->universe_all( + FLERR, + "The scale parameter of the PILE_L thermostat is not supported for method pimd. Delete " + "scale parameter if you do want to use method pimd."); pilescale = utils::numeric(FLERR, arg[i + 1], false, lmp); if (pilescale < 0.0) error->universe_all(FLERR, "Invalid PILE_L scale value for fix pimd/langevin"); @@ -250,10 +253,13 @@ FixPIMDLangevin::FixPIMDLangevin(LAMMPS *lmp, int narg, char **arg) : FLERR, fmt::format("Must not use pressure coupling with {} ensemble", Ensembles[ensemble])); if (method == PIMD && pstat_flag) - error->universe_all(FLERR, "Pressure control has not been supported for method pimd yet. Please set method to nmpimd."); + error->universe_all(FLERR, + "Pressure control has not been supported for method pimd yet. Please set " + "method to nmpimd."); if (method == PIMD && fmmode == NORMAL) - error->universe_all(FLERR, "Normal mode mass is not supported for method pimd. Please set method to nmpimd."); + error->universe_all( + FLERR, "Normal mode mass is not supported for method pimd. Please set method to nmpimd."); /* Initiation */ @@ -487,13 +493,13 @@ void FixPIMDLangevin::setup(int vflag) nmpimd_transform(bufsortedall, x, M_x2xp[universe->iworld]); else if (cmode == MULTI_PROC) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); - } - else if (method == PIMD) { + } else if (method == PIMD) { inter_replica_comm(x); spring_force(); - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } else { + error->universe_all( + FLERR, + "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } collect_xc(); compute_spring_energy(); @@ -549,13 +555,13 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) a_step(); qc_step(); a_step(); - } - else if (method == PIMD) { + } else if (method == PIMD) { q_step(); q_step(); - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } else { + error->universe_all( + FLERR, + "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } } else if (integrator == BAOAB) { if (pstat_flag) { @@ -572,12 +578,12 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]); qc_step(); a_step(); - } - else if (method == PIMD) { + } else if (method == PIMD) { q_step(); - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } else { + error->universe_all( + FLERR, + "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } if (tstat_flag) { o_step(); @@ -587,15 +593,17 @@ void FixPIMDLangevin::initial_integrate(int /*vflag*/) if (method == NMPIMD) { qc_step(); a_step(); - } - else if (method == PIMD) { + } else if (method == PIMD) { q_step(); - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + } else { + error->universe_all( + FLERR, + "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } } else { - error->universe_all(FLERR, "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab integrators are supported!"); + error->universe_all(FLERR, + "Unknown integrator parameter for fix pimd/langevin. Only obabo and baoab " + "integrators are supported!"); } collect_xc(); @@ -652,25 +660,25 @@ void FixPIMDLangevin::post_force(int /*flag*/) tagint *tag = atom->tag; if (method == NMPIMD) { - if (atom->nmax > maxunwrap) reallocate_x_unwrap(); - if (atom->nmax > maxxc) reallocate_xc(); - for (int i = 0; i < nlocal; i++) { - x_unwrap[i][0] = x[i][0]; - x_unwrap[i][1] = x[i][1]; - x_unwrap[i][2] = x[i][2]; - } - if (mapflag) { - for (int i = 0; i < nlocal; i++) { domain->unmap(x_unwrap[i], image[i]); } - } - for (int i = 0; i < nlocal; i++) { - xc[i][0] = xcall[3 * (tag[i] - 1) + 0]; - xc[i][1] = xcall[3 * (tag[i] - 1) + 1]; - xc[i][2] = xcall[3 * (tag[i] - 1) + 2]; - } + if (atom->nmax > maxunwrap) reallocate_x_unwrap(); + if (atom->nmax > maxxc) reallocate_xc(); + for (int i = 0; i < nlocal; i++) { + x_unwrap[i][0] = x[i][0]; + x_unwrap[i][1] = x[i][1]; + x_unwrap[i][2] = x[i][2]; + } + if (mapflag) { + for (int i = 0; i < nlocal; i++) { domain->unmap(x_unwrap[i], image[i]); } + } + for (int i = 0; i < nlocal; i++) { + xc[i][0] = xcall[3 * (tag[i] - 1) + 0]; + xc[i][1] = xcall[3 * (tag[i] - 1) + 1]; + xc[i][2] = xcall[3 * (tag[i] - 1) + 2]; + } - compute_vir(); - compute_cvir(); - compute_t_vir(); + compute_vir(); + compute_cvir(); + compute_t_vir(); } if (method == PIMD) { @@ -1032,13 +1040,12 @@ void FixPIMDLangevin::langevin_init() } for (int i = 0; i < np; i++) { out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_k[i], tau_k[i], - c1_k[i], c2_k[i]); + c1_k[i], c2_k[i]); } - } - else if (method == PIMD) { + } else if (method == PIMD) { for (int i = 0; i < np; i++) { - out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_np / sqrt(fmass), tau, - c1, c2); + out += fmt::format(" {:d} {:.8e} {:.8e} {:.8e} {:.8e}\n", i, _omega_np / sqrt(fmass), + tau, c1, c2); } } if (thermostat == PILE_L) out += "PILE_L thermostat successfully initialized!\n"; @@ -1064,15 +1071,14 @@ void FixPIMDLangevin::o_step() atom->v[i][2] = c1_k[universe->iworld] * atom->v[i][2] + c2_k[universe->iworld] * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); } - } - else if (method == PIMD) { + } else if (method == PIMD) { for (int i = 0; i < nlocal; i++) { - atom->v[i][0] = c1 * atom->v[i][0] + - c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); - atom->v[i][1] = c1 * atom->v[i][1] + - c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); - atom->v[i][2] = c1 * atom->v[i][2] + - c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][0] = + c1 * atom->v[i][0] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][1] = + c1 * atom->v[i][1] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); + atom->v[i][2] = + c1 * atom->v[i][2] + c2 * sqrt(1.0 / mass[type[i]] / beta_np) * random->gaussian(); } } } @@ -1168,39 +1174,38 @@ void FixPIMDLangevin::spring_force() double **x = atom->x; double **f = atom->f; - double* _mass = atom->mass; - int* type = atom->type; + double *_mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; - tagint* tagtmp = atom->tag; + tagint *tagtmp = atom->tag; // printf("iworld = %d, x_last = %d, x_next = %d\n", universe->iworld, x_last, x_next); int *mask = atom->mask; // int idx_tmp = atom->map(1); - for (int i=0; ix; - double *_mass = atom->mass; - int *type = atom->type; - int nlocal = atom->nlocal; + double **x = atom->x; + double *_mass = atom->mass; + int *type = atom->type; + int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) { - spring_energy += 0.5 * _mass[type[i]] * fbond * lam[universe->iworld] * - (x[i][0] * x[i][0] + x[i][1] * x[i][1] + x[i][2] * x[i][2]); - } - MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world); - MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld); - total_spring_energy /= universe->procs_per_world[universe->iworld]; - } - else if (method == PIMD) { + for (int i = 0; i < nlocal; i++) { + spring_energy += 0.5 * _mass[type[i]] * fbond * lam[universe->iworld] * + (x[i][0] * x[i][0] + x[i][1] * x[i][1] + x[i][2] * x[i][2]); + } + MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld); + total_spring_energy /= universe->procs_per_world[universe->iworld]; + } else if (method == PIMD) { total_spring_energy = se_bead = 0.0; MPI_Allreduce(&spring_energy, &se_bead, 1, MPI_DOUBLE, MPI_SUM, world); MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld); - total_spring_energy /= universe->procs_per_world[universe->iworld]; - } - else { - error->universe_all(FLERR, "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); + total_spring_energy /= universe->procs_per_world[universe->iworld]; + } else { + error->universe_all( + FLERR, + "Unknown method parameter for fix pimd/langevin. Only nmpimd and pimd are supported!"); } } From 44b126a87d9800b1ff97a9d4b2a0ec59a93847aa Mon Sep 17 00:00:00 2001 From: jbcouli Date: Thu, 16 Nov 2023 10:35:02 -0700 Subject: [PATCH 11/11] correct typo and link in body particles doc --- doc/src/Howto_body.rst | 6 +++--- doc/src/pair_body_rounded_polyhedron.rst | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/Howto_body.rst b/doc/src/Howto_body.rst index 88fa2d9c97..115b7797c8 100644 --- a/doc/src/Howto_body.rst +++ b/doc/src/Howto_body.rst @@ -170,9 +170,9 @@ with this body style to compute body/body and body/non-body interactions. The *rounded/polygon* body style represents body particles as a 2d polygon with a variable number of N vertices. This style can only be used for 2d models; see the :doc:`boundary ` command. See the -"pair_style body/rounded/polygon" page for a diagram of two -squares with rounded circles at the vertices. Special cases for N = 1 -(circle) and N = 2 (rod with rounded ends) can also be specified. +:doc:`pair_style body/rounded/polygon ` page for +a diagram of two squares with rounded circles at the vertices. Special cases +for N = 1 (circle) and N = 2 (rod with rounded ends) can also be specified. One use of this body style is for 2d discrete element models, as described in :ref:`Fraige `. diff --git a/doc/src/pair_body_rounded_polyhedron.rst b/doc/src/pair_body_rounded_polyhedron.rst index f2f7c1676a..b3eaf72321 100644 --- a/doc/src/pair_body_rounded_polyhedron.rst +++ b/doc/src/pair_body_rounded_polyhedron.rst @@ -40,7 +40,7 @@ rounded/polyhedron particles. This pairwise interaction between the rounded polyhedra is described in :ref:`Wang `, where a polyhedron does not have sharp corners and edges, but is rounded at its vertices and edges by spheres -centered on each vertex with a specified diameter. The edges if the +centered on each vertex with a specified diameter. The edges of the polyhedron are defined between pairs of adjacent vertices. Its faces are defined by a loop of edges. The sphere diameter for each polygon is specified in the data file read by the :doc:`read data `