From 8ecd7e86299cdc9ac3305b2b1f3f5794f592ba1d Mon Sep 17 00:00:00 2001 From: Ofir Blumer Date: Mon, 27 Jan 2025 18:49:30 +0200 Subject: [PATCH] Fixed an error in the primitive estimator --- src/REPLICA/bosonic_exchange.cpp | 42 ++++++++++++++++-------------- src/REPLICA/bosonic_exchange.h | 6 ++--- src/REPLICA/fix_pimdb_langevin.cpp | 4 +-- src/REPLICA/fix_pimdb_nvt.cpp | 4 +-- 4 files changed, 30 insertions(+), 26 deletions(-) diff --git a/src/REPLICA/bosonic_exchange.cpp b/src/REPLICA/bosonic_exchange.cpp index c13be42915..4af76a9771 100644 --- a/src/REPLICA/bosonic_exchange.cpp +++ b/src/REPLICA/bosonic_exchange.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool beta_convention) : Pointers(lmp), - nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic), beta_convention(beta_convention){ + nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic), physical_beta_convention(beta_convention){ memory->create(temp_nbosons_array, nbosons + 1, "BosonicExchange: temp_nbosons_array"); memory->create(E_kn, (nbosons * (nbosons + 1) / 2), "BosonicExchange: E_kn"); memory->create(V, nbosons + 1, "BosonicExchange: V"); @@ -363,32 +363,36 @@ void BosonicExchange::spring_force_interior_bead(double **f) const { double BosonicExchange::prim_estimator() { - temp_nbosons_array[0] = 0.0; + double convention_correction = (physical_beta_convention ? 1 : 1.0 / np); + if (bead_num == 0) { + temp_nbosons_array[0] = 0.0; - for (int m = 1; m < nbosons + 1; ++m) { - double sig = 0.0; - temp_nbosons_array[m] = 0.0; + for (int m = 1; m < nbosons + 1; ++m) { + double sig = 0.0; + temp_nbosons_array[m] = 0.0; - double Elongest = std::numeric_limits::max(); + double Elongest = std::numeric_limits::max(); // Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309) - for (int k = m; k > 0; k--) { - Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]); - } + for (int k = m; k > 0; k--) { + Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]); + } - for (int k = m; k > 0; --k) { - double E_kn_val = get_Enk(m, k); + for (int k = m; k > 0; --k) { + double E_kn_val = get_Enk(m, k); - sig += (temp_nbosons_array[m - k] - E_kn_val) * exp(-beta * (E_kn_val + V[m - k] - Elongest)); + sig += (temp_nbosons_array[m - k] - E_kn_val) * exp(-beta * (E_kn_val + V[m - k] - Elongest)); + } + + double sig_denom_m = m * exp(-beta * (V[m] - Elongest)); + + temp_nbosons_array[m] = sig / sig_denom_m; } - - double sig_denom_m = m * exp(-beta * (V[m] - Elongest)); - - temp_nbosons_array[m] = sig / sig_denom_m; + return convention_correction * (0.5 * domain->dimension * nbosons * np / beta + temp_nbosons_array[nbosons]); + } + else { + return -convention_correction * get_bead_spring_energy(); } - - int convention_correction = (beta_convention ? 1 : np); - return 0.5 * domain->dimension * nbosons * convention_correction / beta + temp_nbosons_array[nbosons]; } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/bosonic_exchange.h b/src/REPLICA/bosonic_exchange.h index a7d5e3ff7d..cd1740cd16 100644 --- a/src/REPLICA/bosonic_exchange.h +++ b/src/REPLICA/bosonic_exchange.h @@ -27,7 +27,6 @@ namespace LAMMPS_NS { double beta, double spring_constant); double get_potential() const; - double get_interior_bead_spring_energy() const; double get_bead_spring_energy() const; void spring_force(double** f) const; @@ -37,6 +36,7 @@ namespace LAMMPS_NS { private: void evaluate_cycle_energies(); void diff_two_beads(const double* x1, int l1, const double* x2, int l2, double diff[3]) const; + double get_interior_bead_spring_energy() const; double distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2) const; double get_Enk(int m, int k) const; void set_Enk(int m, int k, double val); @@ -62,8 +62,8 @@ namespace LAMMPS_NS { // H_physical = H_reduced / P. Note however that the expressions for the various estimators are unaffected by this choice, // so as the algorithm for bosonic exchange. The code below was designed to be compatible with both conventions, // and the choice of convention only affects a single calculation within it. - // Setting the following boolian variable to false amounts to adopting the physical-beta convention. - const bool beta_convention; + // Setting the following boolian variable to false amounts to adopting the reduced-beta convention. + const bool physical_beta_convention; double spring_constant; double beta; diff --git a/src/REPLICA/fix_pimdb_langevin.cpp b/src/REPLICA/fix_pimdb_langevin.cpp index 8cb892e15c..dbd2f19bdf 100755 --- a/src/REPLICA/fix_pimdb_langevin.cpp +++ b/src/REPLICA/fix_pimdb_langevin.cpp @@ -46,7 +46,7 @@ using namespace FixConst; FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) : FixPIMDLangevin(lmp, narg, arg), nbosons(atom->nlocal), - bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, true) + bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, false) { for (int i = 3; i < narg - 1; i += 2) { if ((strcmp(arg[i], "method") == 0) && (strcmp(arg[i+1], "pimd") != 0)) { @@ -128,7 +128,7 @@ void FixPIMDBLangevin::compute_spring_energy() { void FixPIMDBLangevin::compute_t_prim() { - t_prim = (universe->iworld == 0 ? bosonic_exchange.prim_estimator() : 0); + t_prim = bosonic_exchange.prim_estimator(); } /* ---------------------------------------------------------------------- */ diff --git a/src/REPLICA/fix_pimdb_nvt.cpp b/src/REPLICA/fix_pimdb_nvt.cpp index c5920fe698..4f64c56756 100644 --- a/src/REPLICA/fix_pimdb_nvt.cpp +++ b/src/REPLICA/fix_pimdb_nvt.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : FixPIMDNVT(lmp, narg, arg), - bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, false) + bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, true) { virial = 0.; prim = 0.; @@ -56,7 +56,7 @@ void FixPIMDBNVT::pre_spring_force_estimators() { FixPIMDNVT::pre_spring_force_estimators(); spring_energy = bosonic_exchange.get_bead_spring_energy(); - prim = (universe->me == 0 ? bosonic_exchange.prim_estimator() : -bosonic_exchange.get_interior_bead_spring_energy()); + prim = bosonic_exchange.prim_estimator(); } /* ---------------------------------------------------------------------- */