From 9d85a045734dae7deeb3215f4eb3f7203c56e7b2 Mon Sep 17 00:00:00 2001 From: ofirblumer Date: Wed, 9 Oct 2024 17:40:23 +0300 Subject: [PATCH] Corrected most of Yotam's comments, left some for discussion --- src/REPLICA/bosonic_exchange.cpp | 64 ++++++++----------------- src/REPLICA/bosonic_exchange.h | 21 ++++----- src/REPLICA/fix_pimd_nvt.cpp | 11 +++-- src/REPLICA/fix_pimd_nvt.h | 6 ++- src/REPLICA/fix_pimdb_langevin.cpp | 3 +- src/REPLICA/fix_pimdb_nvt.cpp | 76 +++++++----------------------- src/REPLICA/fix_pimdb_nvt.h | 16 ++----- 7 files changed, 65 insertions(+), 132 deletions(-) diff --git a/src/REPLICA/bosonic_exchange.cpp b/src/REPLICA/bosonic_exchange.cpp index 0903add217..cf859c709e 100644 --- a/src/REPLICA/bosonic_exchange.cpp +++ b/src/REPLICA/bosonic_exchange.cpp @@ -7,19 +7,15 @@ using namespace LAMMPS_NS; -BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic) : +BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool iPyConvention) : Pointers(lmp), nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic) { - // CR: temp_nbosons_array needs to be of size nbosons + 1 (my bad, I thought we fixed it here too) - memory->create(temp_nbosons_array, nbosons, "BosonicExchange: temp_nbosons_array"); - memory->create(separate_atom_spring, nbosons, "BosonicExchange: separate_atom_spring"); + 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"); memory->create(V_backwards, nbosons + 1, "BosonicExchange: V_backwards"); memory->create(connection_probabilities, nbosons * nbosons, "BosonicExchange: connection probabilities"); - memory->create(prim_est, nbosons + 1, "BosonicExchange: prim_est"); - // CR: why is this stored in an array, rather than computed on the fly upon get_spring_energy()? - memory->create(spring_energy, nbosons, "BosonicExchange: spring_energy"); + iPyConvention = iPyConvention; } void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next, @@ -28,14 +24,7 @@ void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_ this->x_prev = x_prev; this->x_next = x_next; this->beta = beta; - // CR: confusing that it may or may not be 1/beta. It's used only for the primitive estimator, right? - // CR: If we can't think of a better way, having two versions of the primitive estimator is better imo. - // CR: Or just one version suffices, and the code outside can decide if it needs to divide by P? - // CR: Let's talk about the math here - this->kT = kT; this->spring_constant = spring_constant; - // CR: remove comments - // evaluate_cycle_energies(); if (bead_num == 0 || bead_num == np - 1) { // exterior beads @@ -45,23 +34,19 @@ void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_ evaluate_connection_probabilities(); } - // CR: If spring_energy is not stored in an array, this can be removed. if (0 != bead_num) { - Evaluate_spring_energy(); + calc_total_spring_energy_for_bead(); } } /* ---------------------------------------------------------------------- */ BosonicExchange::~BosonicExchange() { - memory->destroy(prim_est); memory->destroy(connection_probabilities); memory->destroy(V_backwards); memory->destroy(V); memory->destroy(E_kn); - memory->destroy(separate_atom_spring); memory->destroy(temp_nbosons_array); - memory->destroy(spring_energy); } /* ---------------------------------------------------------------------- */ @@ -133,6 +118,7 @@ void BosonicExchange::evaluate_cycle_energies() double BosonicExchange::get_Enk(int m, int k) { // CR: Frankly, there is no reason for this layout. If we can organize it in a more reasonable way, that // CR: would be nice. Not a requirement. + // OB: It indeed seemed like black magic to me at the begining. int end_of_m = m * (m + 1) / 2; return E_kn[end_of_m - k]; } @@ -208,7 +194,6 @@ void BosonicExchange::Evaluate_V_backwards() { V_backwards[0] = V[nbosons]; } - /* ---------------------------------------------------------------------- */ double BosonicExchange::get_potential() const { @@ -217,20 +202,15 @@ double BosonicExchange::get_potential() const { /* ---------------------------------------------------------------------- */ -// YF: this is the spring energy of what? total energy over the beads of the same index? better to have an indicative name -double BosonicExchange::get_spring_energy() const { - double summed_energy = 0; - for (int i = 0; i < nbosons; i++) { - summed_energy += spring_energy[i]; - } - return summed_energy; +double BosonicExchange::get_total_spring_energy_for_bead() { + return spring_energy_for_bead; } - /* ---------------------------------------------------------------------- */ -void BosonicExchange::Evaluate_spring_energy() { +void BosonicExchange::calc_total_spring_energy_for_bead() { + spring_energy_for_bead = 0.; for (int i = 0; i < nbosons; i++) { - spring_energy[i] = 0.5 * spring_constant * distance_squared_two_beads(x, i, x_prev, i); + spring_energy_for_bead += 0.5 * spring_constant * distance_squared_two_beads(x, i, x_prev, i); } } @@ -242,12 +222,6 @@ double BosonicExchange::get_Vn(int n) const { /* ---------------------------------------------------------------------- */ -double BosonicExchange::get_E_kn_serial_order(int i) const { - return E_kn[i]; -} - -/* ---------------------------------------------------------------------- */ - void BosonicExchange::spring_force(double** f) { if (bead_num == np - 1) { spring_force_last_bead(f); @@ -378,12 +352,11 @@ void BosonicExchange::spring_force_interior_bead(double **f) double BosonicExchange::prim_estimator() { - // CR: can use temp_nbosons_array instead of allocating one just for this calculation, - // CR: the important thing is the value you return, intermediate calculations are discarded. - prim_est[0] = 0.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; // Numerical stability (Xiong & Xiong method) double Elongest = std::numeric_limits::max(); @@ -395,23 +368,26 @@ double BosonicExchange::prim_estimator() for (int k = m; k > 0; --k) { double E_kn_val = get_Enk(m, k); - sig += (prim_est[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)); - prim_est[m] = sig / sig_denom_m; + temp_nbosons_array[m] = sig / sig_denom_m; } - - return 0.5 * domain->dimension * nbosons * np * kT + prim_est[nbosons]; + + // OB: is this a better solution? + int convention_correction = (iPyConvention ? 1 : np); + return 0.5 * domain->dimension * nbosons * convention_correction / beta + temp_nbosons_array[nbosons]; } /* ---------------------------------------------------------------------- */ -double BosonicExchange::vir_estimator(double **x, double **f) +double BosonicExchange::vir_estimator(double **x, double **f) const { // CR: I like the decision that the bosonic exchange is responsible for the logic of this estimator. // CR: I don't know how I feel about the code duplication. + // OB: You mean the code duplication? Do you want to better arrange the original NVT and take out the estimator from spring_force? double virial = 0; for (int i = 0; i < nbosons; i++) { virial += -0.5 * (x[i][0] * f[i][0] + x[i][1] * f[i][1] + x[i][2] * f[i][2]); diff --git a/src/REPLICA/bosonic_exchange.h b/src/REPLICA/bosonic_exchange.h index 74a7eb1b32..5a64f9c987 100644 --- a/src/REPLICA/bosonic_exchange.h +++ b/src/REPLICA/bosonic_exchange.h @@ -20,24 +20,23 @@ namespace LAMMPS_NS { class BosonicExchange : protected Pointers { public: - BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic); + BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool iPyConvention); ~BosonicExchange(); void prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next, double beta, double kT, double spring_constant); double get_potential() const; - double get_spring_energy() const; + double get_total_spring_energy_for_bead(); double get_Vn(int n) const; - // CR: I think this method is not used by anyone and can be removed? (probably my bad) - double get_E_kn_serial_order(int i) const; // CR: actually this function can also be const? (probably my bad) + // OB: I don't think it can (see below) void spring_force(double** f); double prim_estimator(); - // CR: add const to arguments and to function - double vir_estimator(double **x, double **f); + // OB: I left **x and **f not const such that they'll be compatible with atom->x and atom->f + double vir_estimator(double **x, double **f) const; private: void evaluate_cycle_energies(); @@ -47,12 +46,13 @@ namespace LAMMPS_NS { void set_Enk(int m, int k, double val); void evaluate_connection_probabilities(); // CR: actually these force functions can also be const? (probably my bad) + // OB: as written now they cannot because they use "diff_two_beads" - it gave me an error... void spring_force_last_bead(double** f); void spring_force_first_bead(double** f); void spring_force_interior_bead(double** f); void Evaluate_VBn(); void Evaluate_V_backwards(); - void Evaluate_spring_energy(); + void calc_total_spring_energy_for_bead(); const int nbosons; const int np; @@ -61,7 +61,6 @@ namespace LAMMPS_NS { double spring_constant; double beta; - double kT; const double* x; const double* x_prev; const double* x_next; @@ -72,11 +71,9 @@ namespace LAMMPS_NS { double* connection_probabilities; double* temp_nbosons_array; - // CR: no longer used, purge! how wonderful - double* separate_atom_spring; + double spring_energy_for_bead; - double *prim_est; - double *spring_energy; + bool iPyConvention; }; diff --git a/src/REPLICA/fix_pimd_nvt.cpp b/src/REPLICA/fix_pimd_nvt.cpp index cbd8ec02f1..107fcb6537 100644 --- a/src/REPLICA/fix_pimd_nvt.cpp +++ b/src/REPLICA/fix_pimd_nvt.cpp @@ -43,8 +43,6 @@ using namespace MathConst; using MathSpecial::powint; -enum { PIMD, NMPIMD, CMD }; - /* ---------------------------------------------------------------------- */ FixPIMDNVT::FixPIMDNVT(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) @@ -220,7 +218,7 @@ void FixPIMDNVT::init() const double Plank = force->hplanck; double hbar = Plank / (2.0 * MY_PI) * sp; - double beta = 1.0 / (Boltzmann * nhc_temp); + beta = 1.0 / (Boltzmann * nhc_temp); double _fbond = 1.0 * np / (beta * beta * hbar * hbar); omega_np = sqrt((double) np) / (hbar * beta) * sqrt(force->mvv2e); @@ -277,6 +275,7 @@ void FixPIMDNVT::post_force(int /*flag*/) for (int j = 0; j < 3; j++) atom->f[i][j] /= np; comm_exec(atom->x); + kinetic_estimators(); spring_force(); if (method == CMD || method == NMPIMD) { @@ -533,6 +532,10 @@ void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector) /* ---------------------------------------------------------------------- */ +void FixPIMDNVT::kinetic_estimators(){} + +/* ---------------------------------------------------------------------- */ + void FixPIMDNVT::spring_force() { spring_energy = 0.0; @@ -545,7 +548,7 @@ void FixPIMDNVT::spring_force() double *xlast = buf_beads[x_last]; double *xnext = buf_beads[x_next]; - // CR: add space back (so that there'd be no diff from the main branch) + virial = 0.0; for (int i = 0; i < nlocal; i++) { diff --git a/src/REPLICA/fix_pimd_nvt.h b/src/REPLICA/fix_pimd_nvt.h index 2cd1397562..bf3dab88e2 100644 --- a/src/REPLICA/fix_pimd_nvt.h +++ b/src/REPLICA/fix_pimd_nvt.h @@ -25,6 +25,8 @@ FixStyle(pimd/nvt,FixPIMDNVT); namespace LAMMPS_NS { +enum { PIMD, NMPIMD, CMD }; + class FixPIMDNVT : public Fix { public: FixPIMDNVT(class LAMMPS *, int, char **); @@ -63,7 +65,8 @@ class FixPIMDNVT : public Fix { double omega_np, fbond, spring_energy, sp, virial; int x_last, x_next; - void spring_force(); + virtual void kinetic_estimators(); + virtual void spring_force(); /* fictitious mass */ @@ -101,6 +104,7 @@ class FixPIMDNVT : public Fix { int nhc_nchain; bool nhc_ready; double nhc_temp, dtv, dtf, t_sys; + double beta; double **nhc_eta; /* coordinates of NH chains for ring-polymer beads */ double **nhc_eta_dot; /* velocities of NH chains */ diff --git a/src/REPLICA/fix_pimdb_langevin.cpp b/src/REPLICA/fix_pimdb_langevin.cpp index e1f18bffe1..e224902233 100755 --- a/src/REPLICA/fix_pimdb_langevin.cpp +++ b/src/REPLICA/fix_pimdb_langevin.cpp @@ -47,8 +47,7 @@ enum{PIMD,NMPIMD,CMD}; FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) : FixPIMDLangevin(lmp, narg, arg), - bosonic_exchange(lmp, atom->nlocal, np, universe->me, - false) + bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, true) { if (method != PIMD) { error->universe_all(FLERR, "Method not supported in fix pimdb/langevin; only method PIMD"); diff --git a/src/REPLICA/fix_pimdb_nvt.cpp b/src/REPLICA/fix_pimdb_nvt.cpp index d1c472a670..7a70ac721d 100644 --- a/src/REPLICA/fix_pimdb_nvt.cpp +++ b/src/REPLICA/fix_pimdb_nvt.cpp @@ -22,40 +22,19 @@ ------------------------------------------------------------------------- */ #include "fix_pimdb_nvt.h" - #include "atom.h" -#include "comm.h" -#include "domain.h" #include "error.h" #include "force.h" -#include "math_const.h" -#include "math_special.h" -#include "memory.h" #include "universe.h" -#include "update.h" - -#include -#include using namespace LAMMPS_NS; -// CR: where are you using these? -using namespace FixConst; -using namespace MathConst; - -// CR: where are you using this? -using MathSpecial::powint; - -// CR: should not define the same enum in two places -enum{PIMD,NMPIMD,CMD}; /* ---------------------------------------------------------------------- */ FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : FixPIMDNVT(lmp, narg, arg), - bosonic_exchange(lmp, atom->nlocal, np, universe->me,false) + bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, false) { - // CR: better to define it as a field in PIMDNVT and use the same (other duplicates the code for that) - beta = 1.0 / force->boltz / nhc_temp; virial = 0.; prim = 0.; spring_energy = 0.; @@ -69,19 +48,13 @@ FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) : FixPIMDBNVT::~FixPIMDBNVT() { } -// CR: add the convention for separating methods -void FixPIMDBNVT::post_force(int /*flag*/) -{ - double **x = atom->x; - double **f = atom->f; - for (int i = 0; i < atom->nlocal; i++) // YF: always indicate the block with curly braces - for (int j = 0; j < 3; j++) atom->f[i][j] /= np; // YF: here too, and newline - // CR: all of this is present also in FixPIMDNVT. Why not override spring_force() instead? - // CR: For the estimators, you can create virtual methods in fix_pimd_nvt and override them here - comm_exec(atom->x); - virial = bosonic_exchange.vir_estimator(x, f); - if (0 == universe->me) // CR: universe->me == 0 +/* ---------------------------------------------------------------------- */ + +void FixPIMDBNVT::kinetic_estimators() +{ + virial = bosonic_exchange.vir_estimator(atom->x, atom->f); + if (universe->me == 0) { prim = bosonic_exchange.prim_estimator(); spring_energy = bosonic_exchange.get_potential(); @@ -90,32 +63,19 @@ void FixPIMDBNVT::post_force(int /*flag*/) // CR: ha, tricky, nice. But where does the constant come into play? // CR: Meaning, what is the "meaning" of the prim value? Because it's not that the sum gives you the kinetic energy, // CR: right? - prim = -bosonic_exchange.get_spring_energy(); - spring_energy = bosonic_exchange.get_spring_energy(); - } - spring_force(x, f); - - // CR: all this too is redundant here, the logic should appear only in fix_pimd_nvt - if (method == NMPIMD) - { - /* forward comm for the force on ghost atoms */ - - nmpimd_fill(atom->f); - - /* inter-partition comm */ - - comm_exec(atom->f); - - /* normal-mode transform */ - - nmpimd_transform(buf_beads, atom->f, M_f2fp[universe->iworld]); + // OB: The sum over the output files from all beads do give the kinetic energy with this implementation... + prim = -bosonic_exchange.get_total_spring_energy_for_bead(); + spring_energy = bosonic_exchange.get_total_spring_energy_for_bead(); } } /* ---------------------------------------------------------------------- */ -void FixPIMDBNVT::spring_force(double **x, double **f) +void FixPIMDBNVT::spring_force() { + double **x = atom->x; + double **f = atom->f; + double *xlast = buf_beads[x_last]; double *xnext = buf_beads[x_next]; double ff = fbond * atom->mass[atom->type[0]]; @@ -125,12 +85,12 @@ void FixPIMDBNVT::spring_force(double **x, double **f) } /* ---------------------------------------------------------------------- */ + double FixPIMDBNVT::compute_vector(int n) { - // CR: call base function for the values that you didn't add - if (n == 0) return spring_energy; - if (n == 1) return t_sys; - if (n == 2) return virial; + if (0 <= n && n < 3) { + return FixPIMDNVT::compute_vector(n); + } // CR: needs to be added also to the documentation. // CR: Reminds that we need to add documentation about the entire bosonic fix if (n == 3) return prim; diff --git a/src/REPLICA/fix_pimdb_nvt.h b/src/REPLICA/fix_pimdb_nvt.h index e39e501660..dffc7b67f2 100644 --- a/src/REPLICA/fix_pimdb_nvt.h +++ b/src/REPLICA/fix_pimdb_nvt.h @@ -29,25 +29,19 @@ class FixPIMDBNVT : public FixPIMDNVT { public: FixPIMDBNVT(class LAMMPS *, int, char **); ~FixPIMDBNVT(); - void post_force(int) override; + // void post_force(int) override; double compute_vector(int) override; protected: - void spring_force(double **x, double **f); // CR: override? - // CR: if it's not used in the case or derived classes, or expected to be, should be private, not protected + void spring_force() override; + void kinetic_estimators() override; + + private: void vir_estimator(double **x, double **f); - // CR: this too BosonicExchange bosonic_exchange; - double beta; - // CR: shouldn't redefine fields that exist in the base class - double virial; double prim; - // CR: here too - double spring_energy; }; - - } // namespace LAMMPS_NS #endif