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

Corrected most of Yotam's comments, left some for discussion

Browse Source
This commit is contained in:
ofirblumer
2024-10-09 17:40:23 +03:00
parent 1fc7cca77c
commit 9d85a04573
7 changed files with 65 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files

64
src/REPLICA/bosonic_exchange.cpp
View File

@ -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<double>::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]);

21
src/REPLICA/bosonic_exchange.h
View File

@ -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;
};

11
src/REPLICA/fix_pimd_nvt.cpp
View File

@ -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++) {

6
src/REPLICA/fix_pimd_nvt.h
View File

@ -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 */

3
src/REPLICA/fix_pimdb_langevin.cpp
View File

@ -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");

76
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -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 <cmath>
#include <cstring>
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;

16
src/REPLICA/fix_pimdb_nvt.h
View File

@ -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