Full depend object implementation of ExchangePotential (i-pi#311)

Refactor the implementation of the ExchangePotential class using depend objects. I think we can now consider that this Fixes i-pi#249 

---------

Co-authored-by: yotamfe <yotam.feldman@gmail.com>

README.md

@@ -22,12 +22,14 @@ Quick Setup
 
 To use i-PI with an existing driver, install and update using `pip`:
 
-Last version::
+Last version:
+
 ```bash
 python -m pip install git+https://github.com/i-pi/i-pi.git
 ```
 
-Last Release::
+Last Release:
+
 ```bash
 pip install -U ipi
 ```

ipi/engine/normalmodes.py

@@ -333,13 +333,6 @@ def bind(self, ensemble, motion, beads=None, forces=None):
             dependencies=[self._pnm, self.beads._sm3, self._nm_factor],
         )
 
-        self._exchange = depend_value(
-            name="exchange",
-            value=None,
-            func=self.get_exchange,
-            dependencies=[self._bosons, self.beads._m3, self._omegan2],
-        )
-
         self._vspring_and_fspring = depend_value(
             name="v_and_fs",
             value=[None, None],
@@ -349,11 +342,19 @@ def bind(self, ensemble, motion, beads=None, forces=None):
                 self._omegak,
                 self._o_omegak,
                 self.beads._m3,
-                self._exchange,
                 self.beads._q,
             ],
         )
 
+        if len(self.bosons) > 0:
+            self.exchange_potential = ExchangePotential(self.nbeads, len(self.bosons))
+            self.exchange_potential.bind(self.beads, self.ensemble, self)
+            self._vspring_and_fspring.add_dependency(
+                self.exchange_potential._vspring_and_fspring
+            )
+        else:
+            self.exchange_potential = None
+
         # just return split potential and force for ease of access
         self._vspring = depend_value(
             name="vspring",
@@ -399,24 +400,6 @@ def resolve_bosons(self):
 
         return bosons_array
 
-    def get_exchange(self):
-        """Sets up an ExchangePotential object to compute bosonic springs"""
-
-        if len(self.bosons) == 0:
-            return None
-
-        masses = dstrip(self.beads.m)[self.bosons]
-        if len(set(masses)) > 1:
-            raise ValueError(
-                "Bosons must have the same mass, found %s for bosons %s"
-                % (str(masses), str(self.bosons))
-            )
-        boson_mass = masses[0]
-        betaP = 1.0 / (self.nbeads * units.Constants.kb * self.ensemble.temp)
-        return ExchangePotential(
-            len(self.bosons), self.nbeads, boson_mass, dstrip(self.omegan2), betaP
-        )
-
     def get_omegan(self):
         """Returns the effective vibrational frequency for the interaction
         between replicas.
@@ -744,12 +727,7 @@ def get_vspring_and_fspring(self):
             fspring += self.transform.nm2b(fspringnm)
 
         if len(self.bosons) > 0:
-            # positions of only the boson atoms
-            qbosons = dstrip(self.beads.q).reshape((self.nbeads, self.natoms, 3))[
-                :, self.bosons, :
-            ]
-            self.exchange.set_coordinates(qbosons)
-            vspring_b, fspring_b = self.exchange.get_vspring_and_fspring()
+            vspring_b, fspring_b = self.exchange_potential.vspring_and_fspring
 
             vspring += vspring_b
             # overwrites the spring forces for the bosonic particles

ipi/engine/properties.py

@@ -1334,7 +1334,7 @@ def get_kintd(self, atom=""):
             atom, iatom, latom, skip_atom_indices=set(self.nm.bosons)
         )
         if bosons_included:
-            res += self.nm.exchange.get_kinetic_td()
+            res += self.nm.exchange_potential.kinetic_td
             ncount += len(bosons_included)
 
         if ncount == 0:
@@ -2705,19 +2705,19 @@ def get_ti_term(self, atom=""):
         return ti
 
     def get_exchange_distinct_prob(self):
-        if self.nm.exchange is None:
+        if self.nm.exchange_potential is None:
             raise Exception("No bosons found for exchange_distinct_prob")
-        return self.nm.exchange.get_distinct_probability()
+        return self.nm.exchange_potential.distinct_probability
 
     def get_exchange_longest_prob(self):
-        if self.nm.exchange is None:
+        if self.nm.exchange_potential is None:
             raise Exception("No bosons found for exchange_all_prob")
-        return self.nm.exchange.get_longest_probability()
+        return self.nm.exchange_potential.longest_probability
 
     def get_fermionic_sign(self):
-        if self.nm.exchange is None:
+        if self.nm.exchange_potential is None:
             raise Exception("No bosons found for fermionic_sign")
-        return self.nm.exchange.get_fermionic_sign()
+        return self.nm.exchange_potential.fermionic_sign
 
 
 class Trajectories: