tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Accepted some of Yotam's code review, there are still open discussions (search OB)

Browse Source
This commit is contained in:
Ofir Blumer
2024-12-23 15:50:27 +02:00
parent 6905bc736b
commit df44ee9504
6 changed files with 39 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

37
src/REPLICA/bosonic_exchange.cpp
View File

@ -7,7 +7,7 @@
using namespace LAMMPS_NS;
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool iPyConvention) :
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic, bool ipy_convention) :
Pointers(lmp),
nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic) {
memory->create(temp_nbosons_array, nbosons + 1, "BosonicExchange: temp_nbosons_array");
@ -15,9 +15,17 @@ BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num,
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");
// CR: use underscore naming convention, not camelCase
// CR: need to document the meaning, because it's tricky and weird
iPyConvention = iPyConvention;
// 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.
// Setting the following boolian variable to false amounts to adopting Tuckerman's convention.
ipy_convention = ipy_convention;
}
void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
@ -36,11 +44,13 @@ void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_
evaluate_connection_probabilities();
}
// CR: bead_num != 0
if (0 != bead_num) {
if (bead_num != 0) {
// CR: why not calculate the total spring force inside get_total_spring_energy_for_bead()?
// CR: Is this a quantity that is queries more than once per step?
// CR: see also comment in fix_pimdb_nvt.cpp
// OB: It should be calculate here, if it is calculated later the energies turn up weird.
// I am not 100% remember why, but I remember there was a problem with the order of this calculation
// and the current implementation works and is also consistent with the distinguishable nvt code.
calc_total_spring_energy_for_bead();
}
}
@ -58,7 +68,7 @@ BosonicExchange::~BosonicExchange() {
/* ---------------------------------------------------------------------- */
void BosonicExchange::diff_two_beads(const double* x1, int l1, const double* x2, int l2,
double diff[3]) {
double diff[3]) const {
l1 = l1 % nbosons;
l2 = l2 % nbosons;
double delx2 = x2[3 * l2 + 0] - x1[3 * l1 + 0];
@ -259,8 +269,7 @@ void BosonicExchange::evaluate_connection_probabilities() {
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_last_bead(double** f)
{
void BosonicExchange::spring_force_last_bead(double** f) const {
const double* x_first_bead = x_next;
const double* x_last_bead = x;
@ -294,8 +303,7 @@ void BosonicExchange::spring_force_last_bead(double** f)
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_first_bead(double** f)
{
void BosonicExchange::spring_force_first_bead(double** f) const {
const double* x_first_bead = x;
const double* x_last_bead = x_prev;
@ -329,8 +337,7 @@ void BosonicExchange::spring_force_first_bead(double** f)
/* ---------------------------------------------------------------------- */
void BosonicExchange::spring_force_interior_bead(double **f)
{
void BosonicExchange::spring_force_interior_bead(double **f) const {
for (int l = 0; l < nbosons; l++) {
double sum_x = 0.0;
double sum_y = 0.0;
@ -382,8 +389,8 @@ double BosonicExchange::prim_estimator()
temp_nbosons_array[m] = sig / sig_denom_m;
}
int convention_correction = (iPyConvention ? 1 : np);
int convention_correction = (ipy_convention ? 1 : np);
return 0.5 * domain->dimension * nbosons * convention_correction / beta + temp_nbosons_array[nbosons];
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */

19
src/REPLICA/bosonic_exchange.h
View File

@ -20,7 +20,7 @@ namespace LAMMPS_NS {
class BosonicExchange : protected Pointers {
public:
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool iPyConvention);
BosonicExchange(class LAMMPS *, int nbosons, int np, int bead_num, bool mic, bool ipy_convention);
~BosonicExchange();
void prepare_with_coordinates(const double* x, const double* x_prev, const double* x_next,
@ -30,25 +30,20 @@ namespace LAMMPS_NS {
double get_total_spring_energy_for_bead();
double get_Vn(int n) 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();
private:
void evaluate_cycle_energies();
void diff_two_beads(const double* x1, int l1, const double* x2, int l2, double diff[3]);
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);
double get_Enk(int m, int k);
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...
// CR: diff_two_beads can be const too?
void spring_force_last_bead(double** f);
void spring_force_first_bead(double** f);
void spring_force_interior_bead(double** f);
void spring_force_last_bead(double** f) const;
void spring_force_first_bead(double** f) const;
void spring_force_interior_bead(double** f) const;
void Evaluate_VBn();
void Evaluate_V_backwards();
void calc_total_spring_energy_for_bead();
@ -72,10 +67,10 @@ namespace LAMMPS_NS {
double* temp_nbosons_array;
double spring_energy_for_bead;
bool iPyConvention;
bool ipy_convention;
};
}
#endif
#endif

4
src/REPLICA/fix_pimd_nvt.cpp
View File

@ -275,7 +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();
estimate_energies();
spring_force();
if (method == CMD || method == NMPIMD) {
@ -532,7 +532,7 @@ void FixPIMDNVT::nmpimd_transform(double **src, double **des, double *vector)
/* ---------------------------------------------------------------------- */
void FixPIMDNVT::kinetic_estimators(){
void FixPIMDNVT::estimate_energies(){
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 kinetic_estimators();
virtual void estimate_energies();
virtual void spring_force();
void vir_estimator();

12
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -34,6 +34,7 @@ using namespace LAMMPS_NS;
FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) :
FixPIMDNVT(lmp, narg, arg),
// CR: apply mic (compatible with previous implementation, and with pimd_nvt)
// OB: Did you mean apply minimum image convention? Because "bosonic_exchange" takes care of it...
bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, false)
{
virial = 0.;
@ -52,19 +53,18 @@ FixPIMDBNVT::~FixPIMDBNVT() {
/* ---------------------------------------------------------------------- */
void FixPIMDBNVT::kinetic_estimators()
void FixPIMDBNVT::estimate_energies()
{
vir_estimator();
if (universe->me == 0)
{
prim = bosonic_exchange.prim_estimator();
// CR: spring_energy is not a kinetic estimator?... rename? separate?
spring_energy = bosonic_exchange.get_potential();
}
else {
// CR: Keep the quantity in local variable to avoid recomputing between prim and spring_energy
prim = -bosonic_exchange.get_total_spring_energy_for_bead();
spring_energy = bosonic_exchange.get_total_spring_energy_for_bead();
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;
}
}
@ -94,4 +94,4 @@ double FixPIMDBNVT::compute_vector(int n)
// CR: Reminds that we need to add documentation about the entire bosonic fix
if (n == 3) return prim;
return 0.0;
}
}

4
src/REPLICA/fix_pimdb_nvt.h
View File

@ -29,13 +29,11 @@ class FixPIMDBNVT : public FixPIMDNVT {
public:
FixPIMDBNVT(class LAMMPS *, int, char **);
~FixPIMDBNVT();
// CR: remove comment
// void post_force(int) override;
double compute_vector(int) override;
protected:
void spring_force() override;
void kinetic_estimators() override;
void estimate_energies() override;
private:
BosonicExchange bosonic_exchange;