Raw outputs for potential and force components

examples/features/path_integral_molecular_dynamics/pimd+mts_liquid_water/input.xml

@@ -1,8 +1,11 @@
 <simulation mode='md' verbosity='high'>
   <output prefix='simulation'>
     <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved, temperature{kelvin}, kinetic_cv, potential, pressure_cv{megapascal}, pot_component_raw(0), pot_component_raw(1) ] </properties>
+    <properties stride='10' filename='pot_lr'>  [ pot_component(1;0), pot_component(1;1), pot_component(1;2), pot_component(1;3) ] </properties> 
+    <properties stride='10' filename='pot_lr-raw'>  [ pot_component_raw(1;0), pot_component_raw(1;1) ] </properties>
     <trajectory filename='pos' stride='10'> positions </trajectory>
-    <trajectory filename='f-0c' stride='10'> forces_component_raw(0) </trajectory>
+    <trajectory filename='f_sr_c' stride='10'> forces_component_raw(0) </trajectory>
+    <trajectory filename='f_lr-raw' stride='10'> forces_component_raw(0,0) </trajectory>
     <checkpoint stride='200'/>
   </output>
   <total_steps>100</total_steps>
@@ -22,14 +25,16 @@
     </initialize>
     <forces>
        <!--
-         MTS setup - apply the fast (short-range) force in the inner loop and the slow (long-range) force in the outer loop. 
-         Note that if the outer loop contains a *correction* to the inner loop the weights should be 
-         [-1,1] (fast force) and [1,0] (slow force)
+           MTS setup - apply the fast (short-range) force in the inner loop and the slow (long-range) force in the 
+           outer loop. Also does ring-polymer contraction, by computing the long-range force on just two beads. 
+
+           Note that if the outer loop contains a *correction* to the inner loop the weights should be 
+           [-1,1] (fast force) and [1,0] (slow force)         
        -->
       <force forcefield='short_range'>
         <mts_weights>[0,1]</mts_weights>
       </force>
-      <force forcefield='long_range'>
+      <force forcefield='long_range' nbeads='2'>
         <mts_weights>[1,0]</mts_weights>
       </force>
     </forces>

examples/temp/pes/pswater/combine-extras/input.xml

@@ -7,8 +7,8 @@
     <trajectory filename='kin' format='xyz' stride='2'> kinetic_cv </trajectory>
     <trajectory filename='dip' stride='2' extra_type='dipole'> extras </trajectory>
     <trajectory filename='raw' stride='2' extra_type='raw'> extras </trajectory>    
-    <trajectory filename='dip0' stride='2' extra_type='dipole'> extras_component(0) </trajectory>
-    <trajectory filename='dip1' stride='2' extra_type='dipole'> extras_component(1) </trajectory>
+    <trajectory filename='dip0' stride='2' extra_type='dipole'> extras_component_raw(0) </trajectory>
+    <trajectory filename='dip1' stride='2' extra_type='dipole'> extras_component_raw(1) </trajectory>
     <checkpoint stride='100' filename='chk' overwrite='true'/>
     <checkpoint stride='2000' filename='restart' overwrite='false'/>
   </output>

examples/temp/pes/pswater/nvt/input.xml

@@ -7,7 +7,7 @@
     <trajectory filename='kin' format='xyz' stride='2'> kinetic_cv </trajectory>
     <trajectory filename='dip' stride='2' extra_type='dipole'> extras </trajectory>
     <trajectory filename='raw' stride='2' extra_type='raw'> extras </trajectory>    
-    <trajectory filename='dip0' stride='2' extra_type='dipole'> extras_component(0) </trajectory>
+    <trajectory filename='dip0' stride='2' extra_type='dipole'> extras_component_raw(0) </trajectory>
     <checkpoint stride='100' filename='chk' overwrite='true'/>
     <checkpoint stride='2000' filename='restart' overwrite='false'/>
   </output>

ipi/engine/forces.py

@@ -1161,41 +1161,44 @@ def get_vir(self):
             vir += v
         return vir
 
-    def pots_component(self, index, weighted=True):
+    def pots_component(self, index, weighted=True, interpolate=True):
         """Fetches the index^th component of the total potential."""
-        if weighted:
-            if self.mforces[index].weight != 0:
-                return self.mforces[index].weight * self.mrpc[index].b2tob1(
-                    self.mforces[index].pots
-                )
+
+        if weighted and self.mforces[index].weight == 0:
+            if interpolate:
+                return np.zeros(self.mrpc[index].nbeads1)
             else:
-                return 0
+                return np.zeros(self.mrpc[index].nbeads2)
         else:
-            return self.mrpc[index].b2tob1(self.mforces[index].pots)
-
-    def forces_component(self, index, weighted=True):
+            pots = dstrip(self.mforces[index].pots).copy()
+            if weighted:
+                pots *= self.mforces[index].weight
+            if interpolate:
+                pots = self.mrpc[index].b2tob1(pots)
+            return pots
+
+    def forces_component(self, index, weighted=True, interpolate=True):
         """Fetches the index^th component of the total force."""
-        if weighted:
-            if self.mforces[index].weight != 0:
-                return self.mforces[index].weight * self.mrpc[index].b2tob1(
-                    dstrip(self.mforces[index].f)
-                )
+
+        if weighted and self.mforces[index].weight == 0:
+            if interpolate:
+                return np.zeros((self.mrpc[index].nbeads1, 3 * self.natoms))
             else:
-                return np.zeros((self.nbeads, self.natoms * 3), float)
+                return np.zeros((self.mrpc[index].nbeads2, 3 * self.natoms))
         else:
-            return self.mrpc[index].b2tob1(dstrip(self.mforces[index].f))
+            forces = dstrip(self.mforces[index].f).copy()
+            if weighted:
+                forces *= self.mforces[index].weight
+            if interpolate:
+                forces = self.mrpc[index].b2tob1(forces)
+            return forces
 
     def extras_component(self, index):
-        """Fetches extras that are computed for one specific force component."""
+        """Fetches extras that are computed for one specific force component.
+
+        Does not attempt to apply weights or interpolate, always returns raw stuff.
+        """
 
-        if self.nbeads != self.mforces[index].nbeads:
-            raise ValueError(
-                "Cannot fetch extras for a component when using ring polymer contraction"
-            )
-        if self.mforces[index].weight == 0:
-            raise ValueError(
-                "Cannot fetch extras for a component that has not been computed because of zero weight"
-            )
         return self.mforces[index].extras
 
     def forcesvirs_4th_order(self, index):

ipi/engine/properties.py

@@ -361,7 +361,8 @@ def __init__(self):
                 "longhelp": """The contribution to the system potential from one of the force components.
                        Takes one mandatory argument index (zero-based) that indicates which component of the
                        potential must be returned. The optional argument 'bead' will print the potential associated
-                       with the specified bead. If the potential is weighed, the weight will be applied. """,
+                       with the specified bead (interpolated to the full ring polymer). 
+                       If the potential is weighed, the weight will be applied. """,
                 "func": (
                     lambda index, bead="-1": (
                         self.forces.pots_component(int(index)).sum() / self.beads.nbeads
@@ -376,14 +377,16 @@ def __init__(self):
                 "longhelp": """The contribution to the system potential from one of the
                        force components. Takes one mandatory argument index (zero-based) that indicates
                        which component of the potential must be returned. The optional argument 'bead'
-                       will print the potential associated with the specified bead. Potential weights
-                       will not be applied. """,
+                       will print the potential associated with the specified bead, at the level of 
+                       discretization of the given component. Potential weights will not be applied. """,
                 "func": (
                     lambda index, bead="-1": (
-                        self.forces.pots_component(int(index), False).sum()
+                        self.forces.pots_component(int(index), False, False).sum()
                         / self.beads.nbeads
                         if int(bead) < 0
-                        else self.forces.pots_component(int(index), False)[int(bead)]
+                        else self.forces.pots_component(int(index), False, False)[
+                            int(bead)
+                        ]
                     )
                 ),
             },
@@ -2838,7 +2841,7 @@ def __init__(self):
                 "help": """The contribution to the system forces from one of the force components.
                        Takes one mandatory argument index (zero-based) that indicates which component of the
                        potential must be returned. The optional argument 'bead' will print the potential associated
-                       with the specified bead, otherwise the centoid force is computed. 
+                       with the specified bead (interpolated to the full ring polymer), otherwise the centoid force is computed. 
                        If the potential is weighed, the weight will be applied. """,
                 "func": lambda index, bead="-1": (
                     self.system.forces.forces_component(int(index)).sum(axis=0)
@@ -2852,20 +2855,24 @@ def __init__(self):
                 "help": """The contribution to the system forces from one of the force components.
                        Takes one mandatory argument index (zero-based) that indicates which component of the
                        potential must be returned. The optional argument 'bead' will print the potential associated
-                       with the specified bead, otherwise the centoid force is computed. 
-                       If the potential is weighed, the weight will be applied. """,
+                       with the specified bead (with the level of discretization of the component), otherwise the 
+                       centoid force is computed. The weight of the potential is not applied. """,
                 "func": lambda index, bead="-1": (
-                    self.system.forces.forces_component(int(index), False).sum(axis=0)
+                    self.system.forces.forces_component(int(index), False, False).sum(
+                        axis=0
+                    )
                     / self.system.beads.nbeads
                     if int(bead) < 0
-                    else self.system.forces.forces_component(int(index), False)[
+                    else self.system.forces.forces_component(int(index), False, False)[
                         int(bead)
                     ]
                 ),
             },
-            "extras_component": {
-                "help": """The additional data returned by the client code, printed verbatim or expanded as a dictionary. See "extras". 
-                           Fetches the extras from a specific force component, indicated in parentheses [extras_component(idx)]. """,
+            "extras_component_raw": {
+                "help": """The additional data returned by the client code, printed verbatim or expanded 
+                           as a dictionary. See "extras". 
+                           Fetches the extras from a specific force component, indicated in parentheses 
+                           [extras_component(idx)]. Never applies weighting or contraction""",
                 "func": (lambda idx: self.system.forces.extras_component(int(idx))),
             },
             "extras_bias": {

ipi/utils/nmtransform.py

@@ -247,6 +247,9 @@ def __init__(self, nbeads1, nbeads2, open_paths=None):
            nbeads2: The rescaled number of beads.
         """
 
+        self.nbeads1 = nbeads1
+        self.nbeads2 = nbeads2
+
         if open_paths is None:
             open_paths = []
         self._open = open_paths