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

Fixed an error in the primitive estimator

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

42
src/REPLICA/bosonic_exchange.cpp
View File

@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
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), beta_convention(beta_convention){
nbosons(nbosons), np(np), bead_num(bead_num), apply_minimum_image(mic), physical_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");
@ -363,32 +363,36 @@ void BosonicExchange::spring_force_interior_bead(double **f) const {
double BosonicExchange::prim_estimator()
{
temp_nbosons_array[0] = 0.0;
double convention_correction = (physical_beta_convention ? 1 : 1.0 / np);
if (bead_num == 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;
for (int m = 1; m < nbosons + 1; ++m) {
double sig = 0.0;
temp_nbosons_array[m] = 0.0;
double Elongest = std::numeric_limits<double>::max();
double Elongest = std::numeric_limits<double>::max();
// Numerical stability (Xiong & Xiong, doi:10.1103/PhysRevE.106.025309)
for (int k = m; k > 0; k--) {
Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]);
}
for (int k = m; k > 0; k--) {
Elongest = std::min(Elongest, get_Enk(m, k) + V[m - k]);
}
for (int k = m; k > 0; --k) {
double E_kn_val = get_Enk(m, k);
for (int k = m; k > 0; --k) {
double E_kn_val = get_Enk(m, k);
sig += (temp_nbosons_array[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));
temp_nbosons_array[m] = sig / sig_denom_m;
}
double sig_denom_m = m * exp(-beta * (V[m] - Elongest));
temp_nbosons_array[m] = sig / sig_denom_m;
return convention_correction * (0.5 * domain->dimension * nbosons * np / beta + temp_nbosons_array[nbosons]);
}
else {
return -convention_correction * get_bead_spring_energy();
}
int convention_correction = (beta_convention ? 1 : np);
return 0.5 * domain->dimension * nbosons * convention_correction / beta + temp_nbosons_array[nbosons];
}
/* ---------------------------------------------------------------------- */

6
src/REPLICA/bosonic_exchange.h
View File

@ -27,7 +27,6 @@ namespace LAMMPS_NS {
double beta, double spring_constant);
double get_potential() const;
double get_interior_bead_spring_energy() const;
double get_bead_spring_energy() const;
void spring_force(double** f) const;
@ -37,6 +36,7 @@ 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 get_interior_bead_spring_energy() 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);
@ -62,8 +62,8 @@ namespace LAMMPS_NS {
// 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;
// Setting the following boolian variable to false amounts to adopting the reduced-beta convention.
const bool physical_beta_convention;
double spring_constant;
double beta;

4
src/REPLICA/fix_pimdb_langevin.cpp
View File

@ -46,7 +46,7 @@ using namespace FixConst;
FixPIMDBLangevin::FixPIMDBLangevin(LAMMPS *lmp, int narg, char **arg) :
FixPIMDLangevin(lmp, narg, arg), nbosons(atom->nlocal),
bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, true)
bosonic_exchange(lmp, atom->nlocal, np, universe->me, false, false)
{
for (int i = 3; i < narg - 1; i += 2) {
if ((strcmp(arg[i], "method") == 0) && (strcmp(arg[i+1], "pimd") != 0)) {
@ -128,7 +128,7 @@ void FixPIMDBLangevin::compute_spring_energy() {
void FixPIMDBLangevin::compute_t_prim()
{
t_prim = (universe->iworld == 0 ? bosonic_exchange.prim_estimator() : 0);
t_prim = bosonic_exchange.prim_estimator();
}
/* ---------------------------------------------------------------------- */

4
src/REPLICA/fix_pimdb_nvt.cpp
View File

@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
FixPIMDBNVT::FixPIMDBNVT(LAMMPS *lmp, int narg, char **arg) :
FixPIMDNVT(lmp, narg, arg),
bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, false)
bosonic_exchange(lmp, atom->nlocal, np, universe->me, true, true)
{
virial = 0.;
prim = 0.;
@ -56,7 +56,7 @@ void FixPIMDBNVT::pre_spring_force_estimators()
{
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());
prim = bosonic_exchange.prim_estimator();
}
/* ---------------------------------------------------------------------- */