ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Corrected PIMD-B after CR

Browse Source
This commit is contained in:
Ofir Blumer
2025-01-27 15:30:24 +02:00
parent b2def45011
commit d91a75a9af
8 changed files with 49 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files

35
src/REPLICA/bosonic_exchange.cpp
View File

@ -32,33 +32,18 @@
using namespace LAMMPS_NS;
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool ipy_convention) :
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) {
nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic), 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");
memory->create(V_backwards, nbosons + 1, "BosonicExchange: V_backwards");
memory->create(connection_probabilities, nbosons * nbosons, "BosonicExchange: connection probabilities");
// In the "i-Pi convention" [J. Chem. Phys. 133, 124104 (2010); also J. Chem. Phys. 74, 4078-4095 (1981)],
// the Boltzmann exponents have the form exp[-(beta/P)H], where H is the classical Hamiltonian of the
// ring polymers. This results in a canonical distribution at P times the physical temperature.
// In contrast, "Tuckerman's convention" [J. Chem. Phys. 99, 2796-2808 (1993)] uses weights of the form exp(-beta*H),
// such that the temperature of the canonical ensemble coincides with the physical temperature.
// Notably, the classical Hamiltonians of the two conventions differ, with the spring constant
// in the i-Pi convention being P times larger than that in Tuckerman's convention. Additionally, the i-Pi convention
// lacks a 1/P prefactor in front of the external potential. The Hamiltonians of the two conventions are related through
// H_tuckerman = H_ipi / 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 Tuckerman's convention.
ipy_convention = ipy_convention;
// CR: I'm not sure the credit in the naming is the best, let's discuss
}
void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
double beta, double kT, double spring_constant) {
double beta, double spring_constant) {
this->x = x;
this->x_prev = x_prev;
this->x_next = x_next;
@ -104,7 +89,7 @@ void BosonicExchange::diff_two_beads(const double* x1, int l1, const double* x2,
/* ---------------------------------------------------------------------- */
double BosonicExchange::distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2) {
double BosonicExchange::distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2) const {
double diff[3];
diff_two_beads(x1, l1, x2, l2, diff);
return diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2];
@ -150,7 +135,7 @@ void BosonicExchange::evaluate_cycle_energies()
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_Enk(int m, int k) {
double BosonicExchange::get_Enk(int m, int k) const {
int end_of_m = m * (m + 1) / 2;
return E_kn[end_of_m - k];
}
@ -234,8 +219,7 @@ double BosonicExchange::get_potential() const {
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_total_spring_energy_for_bead() {
// CR: this function doesn't make sense for the last bead because of exchange?
double BosonicExchange::get_interior_bead_spring_energy() const {
double spring_energy_for_bead = 0.;
for (int i = 0; i < nbosons; i++) {
spring_energy_for_bead += 0.5 * spring_constant * distance_squared_two_beads(x, i, x_prev, i);
@ -245,8 +229,9 @@ double BosonicExchange::get_total_spring_energy_for_bead() {
/* ---------------------------------------------------------------------- */
double BosonicExchange::get_Vn(int n) const {
return V[n];
double BosonicExchange::get_bead_spring_energy() const {
double spring_energy_for_bead = (bead_num == 0 ? get_potential() : get_interior_bead_spring_energy());
return spring_energy_for_bead;
}
/* ---------------------------------------------------------------------- */
@ -402,7 +387,7 @@ double BosonicExchange::prim_estimator()
temp_nbosons_array[m] = sig / sig_denom_m;
}
int convention_correction = (ipy_convention ? 1 : np);
int convention_correction = (beta_convention ? 1 : np);
return 0.5 * domain->dimension * nbosons * convention_correction / beta + temp_nbosons_array[nbosons];
}

32
src/REPLICA/bosonic_exchange.h
View File

@ -20,17 +20,15 @@ namespace LAMMPS_NS {
class BosonicExchange : protected Pointers {
public:
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool ipy_convention);
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool beta_convention);
~BosonicExchange();
// CR: remove parameter kT
void prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
double beta, double kT, double spring_constant);
double beta, double spring_constant);
double get_potential() const;
double get_total_spring_energy_for_bead(); // CR: const
// CR: no usages of this function? Should at most return the potential on all particles, the rest should be encapsulated in the class (method should be private)
double get_Vn(int n) const;
double get_interior_bead_spring_energy() const;
double get_bead_spring_energy() const;
void spring_force(double** f) const;
@ -39,8 +37,8 @@ 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 distance_squared_two_beads(const double* x1, int l1, const double* x2, int l2); // CR: const
double get_Enk(int m, int k); // CR: 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);
void evaluate_connection_probabilities();
void spring_force_last_bead(double** f) const;
@ -53,6 +51,19 @@ namespace LAMMPS_NS {
const int np;
const int bead_num;
const bool apply_minimum_image;
// In the "reduced-beta convention" [e.g. in J. Chem. Phys. 133, 124104 (2010); also J. Chem. Phys. 74, 4078-4095 (1981)],
// the Boltzmann exponents have the form exp[-(beta/P)H], where H is the classical Hamiltonian of the
// ring polymers. This results in a canonical distribution at P times the physical temperature.
// In contrast, the "physical-beta convention" [e.g. in J. Chem. Phys. 99, 2796-2808 (1993)] uses weights of the form exp(-beta*H),
// such that the temperature of the canonical ensemble coincides with the physical temperature.
// Notably, the classical Hamiltonians of the two conventions differ, with the spring constant
// in the reduced-beta convention being P times larger than that in the physical-beta convention. Additionally, the reduced-beta convention
// lacks a 1/P prefactor in front of the external potential. The Hamiltonians of the two conventions are related through
// 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;
double spring_constant;
double beta;
@ -66,12 +77,7 @@ namespace LAMMPS_NS {
double* connection_probabilities;
double* temp_nbosons_array;
bool ipy_convention; // CR: const. Also put earlier with all the other fields
// CR: spelled iPi (but perhaps rename anyways)
};
}
#endif

5
src/REPLICA/fix_pimd_nvt.cpp
View File

@ -282,8 +282,7 @@ void FixPIMDNVT::post_force(int /*flag*/)
for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
prepare_coordinates();
// CR: rename to pre_spring_force_estimators() (better conveys the reason for this function)
estimate_energies();
pre_spring_force_estimators();
spring_force();
if (method == CMD || method == NMPIMD) {
@ -540,7 +539,7 @@ void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector)
/* ---------------------------------------------------------------------- */
void FixPIMDNVT::estimate_energies(){
void FixPIMDNVT::pre_spring_force_estimators(){
vir_estimator();
}

2
src/REPLICA/fix_pimd_nvt.h
View File

@ -65,7 +65,7 @@ class FixPIMDNVT : public Fix {
double omega_np, fbond, spring_energy, sp, virial;
int x_last, x_next;
virtual void prepare_coordinates();
virtual void estimate_energies();
virtual void pre_spring_force_estimators();
virtual void spring_force();
void vir_estimator();

17
src/REPLICA/fix_pimdb_langevin.cpp
View File

@ -45,7 +45,7 @@ using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) :
FixPIMDLangevin(lmp, narg, arg),
FixPIMDLangevin(lmp, narg, arg), nbosons(atom->nlocal),
bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, true)
{
for (int i = 3; i < narg - 1; i += 2) {
@ -69,11 +69,8 @@ FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) :
method = PIMD;
// CR: let's add a test that it fails if the user asks for a property larger than this (also in nvt)
size_vector = 6;
nbosons = atom->nlocal;
memory->create(f_tag_order, nbosons, 3, "FixPIMDBLangevin:f_tag_order");
}
@ -96,7 +93,7 @@ void FixPIMDBLangevin::prepare_coordinates()
bosonic_exchange.prepare_with_coordinates(me_bead_positions,
last_bead_positions, next_bead_positions,
beta_np, 1 / beta, ff);
beta_np, ff);
}
/* ---------------------------------------------------------------------- */
@ -122,8 +119,7 @@ void FixPIMDBLangevin::spring_force() {
/* ---------------------------------------------------------------------- */
void FixPIMDBLangevin::compute_spring_energy() {
// CR: in both cases you call bosonic_exchange.something, so put this logic inside bosonic_exchange
se_bead = (universe->iworld == 0 ? bosonic_exchange.get_potential() : bosonic_exchange.get_total_spring_energy_for_bead());
se_bead = bosonic_exchange.get_bead_spring_energy();
MPI_Allreduce(&se_bead, &total_spring_energy, 1, MPI_DOUBLE, MPI_SUM, universe->uworld);
total_spring_energy /= universe->procs_per_world[universe->iworld];
}
@ -132,9 +128,7 @@ void FixPIMDBLangevin::compute_spring_energy() {
void FixPIMDBLangevin::compute_t_prim()
{
if (universe->iworld == 0)
t_prim = bosonic_exchange.prim_estimator();
// CR: else assign 0 (to not be dependent on the context of the call, and this is a choice you're making here)
t_prim = (universe->iworld == 0 ? bosonic_exchange.prim_estimator() : 0);
}
/* ---------------------------------------------------------------------- */
@ -144,5 +138,8 @@ double FixPIMDBLangevin::compute_vector(int n)
if (0 <= n && n < 6) {
return FixPIMDLangevin::compute_vector(n);
}
else {
error->universe_all(FLERR, "Fix only has 6 outputs!");
}
return 0.0;
}

2
src/REPLICA/fix_pimdb_langevin.h
View File

@ -39,7 +39,7 @@ protected:
void spring_force() override;
private:
int nbosons; // CR: const
const int nbosons;
BosonicExchange bosonic_exchange;
double** f_tag_order;
};

22
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -52,20 +52,11 @@ FixPIMDBNVT::~FixPIMDBNVT() {
/* ---------------------------------------------------------------------- */
void FixPIMDBNVT::estimate_energies()
void FixPIMDBNVT::pre_spring_force_estimators()
{
// CR: call super function instead of calling again to vir_estimator
vir_estimator();
if (universe->me == 0)
{
prim = bosonic_exchange.prim_estimator();
spring_energy = bosonic_exchange.get_potential();
}
else {
double total_spring_energy_for_bead = bosonic_exchange.get_total_spring_energy_for_bead();
prim = -total_spring_energy_for_bead;
spring_energy = total_spring_energy_for_bead;
}
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());
}
/* ---------------------------------------------------------------------- */
@ -77,7 +68,7 @@ void FixPIMDBNVT::prepare_coordinates()
double *xlast = buf_beads[x_last];
double *xnext = buf_beads[x_next];
double ff = fbond * atom->mass[atom->type[0]];
bosonic_exchange.prepare_with_coordinates(*x, xlast, xnext, beta, 1 / beta, -ff);
bosonic_exchange.prepare_with_coordinates(*x, xlast, xnext, beta, -ff);
}
void FixPIMDBNVT::spring_force()
@ -98,6 +89,9 @@ double FixPIMDBNVT::compute_vector(int n)
if (n == 3) {
return prim;
}
else {
error->universe_all(FLERR, "Fix only has 4 outputs!");
}
return 0.0;
}

2
src/REPLICA/fix_pimdb_nvt.h
View File

@ -34,7 +34,7 @@ class FixPIMDBNVT : public FixPIMDNVT {
protected:
void prepare_coordinates() override;
void spring_force() override;
void estimate_energies() override;
void pre_spring_force_estimators() override;
private:
BosonicExchange bosonic_exchange;