Got rid of the calc function in bosonic_exchange, the calculation is performed in the get function.

Added prepare_coordinates in pimd/nvt and pimd/langevin for that purpose

src/REPLICA/bosonic_exchange.cpp

@@ -70,16 +70,6 @@ void BosonicExchange::prepare_with_coordinates(const double* x, const double* x_
         Evaluate_V_backwards();
         evaluate_connection_probabilities();
     }
-
-    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();
-    }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -243,15 +233,11 @@ double BosonicExchange::get_potential() const {
 /* ---------------------------------------------------------------------- */
 
 double BosonicExchange::get_total_spring_energy_for_bead() {
-    return spring_energy_for_bead;
-}
-/* ---------------------------------------------------------------------- */
-
-void BosonicExchange::calc_total_spring_energy_for_bead() {
-    spring_energy_for_bead = 0.;
+    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);
     }
+    return spring_energy_for_bead;
 }
 
 /* ---------------------------------------------------------------------- */

src/REPLICA/bosonic_exchange.h

@@ -46,7 +46,6 @@ namespace LAMMPS_NS {
         void spring_force_interior_bead(double** f) const;
         void Evaluate_VBn();
         void Evaluate_V_backwards();
-        void calc_total_spring_energy_for_bead();
 
         const int nbosons;
         const int np;
@@ -65,7 +64,6 @@ namespace LAMMPS_NS {
         double* connection_probabilities;
 
         double* temp_nbosons_array;
-        double spring_energy_for_bead;
 
         bool ipy_convention;
     };

src/REPLICA/fix_pimd_langevin.cpp

@@ -513,7 +513,7 @@ void FixPIMDLangevin::setup(int vflag)
     else if (cmode == MULTI_PROC)
       nmpimd_transform(bufbeads, x, M_x2xp[universe->iworld]);
   } else if (method == PIMD) {
-    inter_replica_comm(x);
+    prepare_coordinates();
     spring_force();
   } else {
     error->universe_all(
@@ -670,6 +670,13 @@ void FixPIMDLangevin::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
+void FixPIMDLangevin::prepare_coordinates()
+{
+  inter_replica_comm(atom->x);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixPIMDLangevin::post_force(int /*flag*/)
 {
   int nlocal = atom->nlocal;
@@ -703,7 +710,7 @@ void FixPIMDLangevin::post_force(int /*flag*/)
     if (mapflag) {
       for (int i = 0; i < nlocal; i++) { domain->unmap(x[i], image[i]); }
     }
-    inter_replica_comm(x);
+    prepare_coordinates();   
     spring_force();
     compute_spring_energy();
     compute_t_prim();

src/REPLICA/fix_pimd_langevin.h

@@ -103,6 +103,7 @@ class FixPIMDLangevin : public Fix {
   int *counts, *displacements;
 
   void comm_init();
+  virtual void prepare_coordinates();
   void inter_replica_comm(double **ptr);
   void virtual spring_force();
 

src/REPLICA/fix_pimd_nvt.cpp

@@ -269,12 +269,19 @@ void FixPIMDNVT::final_integrate()
 
 /* ---------------------------------------------------------------------- */
 
+void FixPIMDNVT::prepare_coordinates()
+{
+  comm_exec(atom->x);
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixPIMDNVT::post_force(int /*flag*/)
 {
   for (int i = 0; i < atom->nlocal; i++)
     for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
 
-  comm_exec(atom->x);
+  prepare_coordinates();
   estimate_energies();
   spring_force();
 

src/REPLICA/fix_pimd_nvt.h

@@ -64,7 +64,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 spring_force();
   void vir_estimator();

src/REPLICA/fix_pimdb_langevin.cpp

@@ -84,16 +84,23 @@ FixPIMDBLangevin::~FixPIMDBLangevin() {
 
 /* ---------------------------------------------------------------------- */
 
-void FixPIMDBLangevin::spring_force() {
-    double ff = fbond * atom->mass[atom->type[0]]; // TODO: ensure that all masses are the same
-    int nlocal = atom->nlocal;
-    double* me_bead_positions = *(atom->x);
-    double* last_bead_positions = &bufsortedall[x_last * nlocal][0];
-    double* next_bead_positions = &bufsortedall[x_next * nlocal][0];
+void FixPIMDBLangevin::prepare_coordinates()
+{
+  inter_replica_comm(atom->x);
+  double ff = fbond * atom->mass[atom->type[0]];
+  int nlocal = atom->nlocal;
+  double* me_bead_positions = *(atom->x);
+  double* last_bead_positions = &bufsortedall[x_last * nlocal][0];
+  double* next_bead_positions = &bufsortedall[x_next * nlocal][0];
 
-    bosonic_exchange.prepare_with_coordinates(me_bead_positions,
+  bosonic_exchange.prepare_with_coordinates(me_bead_positions,
                                               last_bead_positions, next_bead_positions,
                                               beta_np, 1 / beta, ff);
+}
+ 
+/* ---------------------------------------------------------------------- */
+
+void FixPIMDBLangevin::spring_force() {
 
     for (int i = 0; i < nbosons; i++) {
         f_tag_order[i][0] = 0.0;

src/REPLICA/fix_pimdb_langevin.h

@@ -26,7 +26,7 @@ FixStyle(pimdb/langevin,FixPIMDBLangevin);
 namespace LAMMPS_NS {
 
 class FixPIMDBLangevin : public FixPIMDLangevin {
- public:
+public:
     FixPIMDBLangevin(class LAMMPS *, int, char **);
     ~FixPIMDBLangevin();
 
@@ -34,7 +34,8 @@ class FixPIMDBLangevin : public FixPIMDLangevin {
     void compute_spring_energy() override;
     void compute_t_prim() override;
 
-        protected:
+protected:
+    void prepare_coordinates() override;
     void spring_force() override;
 
 private:

src/REPLICA/fix_pimdb_nvt.cpp

@@ -69,16 +69,20 @@ void FixPIMDBNVT::estimate_energies()
 
 /* ---------------------------------------------------------------------- */
 
-void FixPIMDBNVT::spring_force()
+void FixPIMDBNVT::prepare_coordinates()
 {
+  comm_exec(atom->x);
   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]]; 
-  
+  double ff = fbond * atom->mass[atom->type[0]];
   bosonic_exchange.prepare_with_coordinates(*x, xlast, xnext, beta, 1 / beta, -ff);
+}
+
+void FixPIMDBNVT::spring_force()
+{
+  double **f = atom->f;
+
   bosonic_exchange.spring_force(f);
 }
 

src/REPLICA/fix_pimdb_nvt.h

@@ -32,6 +32,7 @@ class FixPIMDBNVT : public FixPIMDNVT {
     double compute_vector(int) override;
 
  protected:
+    void prepare_coordinates() override;
     void spring_force() override;
     void estimate_energies() override;