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

cr

Browse Source
This commit is contained in:
yotamfe
2024-10-09 15:15:36 +03:00
parent c645a6b841
commit 1fc7cca77c
5 changed files with 43 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/REPLICA/bosonic_exchange.cpp
View File

@ -10,6 +10,7 @@ using namespace LAMMPS_NS;
BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num, bool mic) :
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(E_kn, (nbosons * (nbosons + 1) / 2), "BosonicExchange: E_kn");
@ -17,6 +18,7 @@ BosonicExchange::BosonicExchange(LAMMPS *lmp, int nbosons, int np, int bead_num,
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");
}
@ -26,8 +28,13 @@ 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) {
@ -38,6 +45,7 @@ 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();
}
@ -123,6 +131,8 @@ 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.
int end_of_m = m * (m + 1) / 2;
return E_kn[end_of_m - k];
}
@ -207,6 +217,7 @@ 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++) {
@ -367,6 +378,8 @@ 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;
for (int m = 1; m < nbosons + 1; ++m) {
@ -397,6 +410,8 @@ double BosonicExchange::prim_estimator()
double BosonicExchange::vir_estimator(double **x, double **f)
{
// 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.
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]);

5
src/REPLICA/bosonic_exchange.h
View File

@ -29,11 +29,14 @@ namespace LAMMPS_NS {
double get_potential() const;
double get_spring_energy() const;
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)
void spring_force(double** f);
double prim_estimator();
// CR: add const to arguments and to function
double vir_estimator(double **x, double **f);
private:
@ -43,6 +46,7 @@ namespace LAMMPS_NS {
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)
void spring_force_last_bead(double** f);
void spring_force_first_bead(double** f);
void spring_force_interior_bead(double** f);
@ -68,6 +72,7 @@ namespace LAMMPS_NS {
double* connection_probabilities;
double* temp_nbosons_array;
// CR: no longer used, purge! how wonderful
double* separate_atom_spring;
double *prim_est;

1
src/REPLICA/fix_pimd_nvt.cpp
View File

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

20
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -38,11 +38,14 @@
#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};
/* ---------------------------------------------------------------------- */
@ -51,6 +54,7 @@ FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) :
FixPIMDNVT(lmp, narg, arg),
bosonic_exchange(lmp, atom->nlocal, np, universe->me,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.;
@ -65,26 +69,33 @@ 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++)
for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
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)
if (0 == universe->me) // CR: universe->me == 0
{
prim = bosonic_exchange.prim_estimator();
spring_energy = bosonic_exchange.get_potential();
}
else {
// 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 */
@ -116,9 +127,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;
// 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;
return 0.0;
}

6
src/REPLICA/fix_pimdb_nvt.h
View File

@ -33,12 +33,16 @@ class FixPIMDBNVT : public FixPIMDNVT {
double compute_vector(int) override;
protected:
void spring_force(double **x, double **f);
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 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;
};