diff --git a/examples/profiling/remd_scpimd_gas/init.xyz b/examples/profiling/remd_scpimd_gas/init.xyz
index f2f3870f8..f009e9fc7 100644
--- a/examples/profiling/remd_scpimd_gas/init.xyz
+++ b/examples/profiling/remd_scpimd_gas/init.xyz
@@ -1,6 +1,5 @@
 8
-# CELL(abcABC):   50.0 50.0 50.0 90.00000    90.00000    90.00000  
-H 0 0 0
+# CELL(abcABC):   10.0 10.0 10.0 90.00000    90.00000    90.00000  
 H 0 0 1
 H 0 0 2
 H 0 0 3
@@ -8,3 +7,4 @@ H 0 0 4
 H 0 0 5
 H 0 0 6
 H 0 0 7
+H 0 0 8
diff --git a/examples/profiling/remd_scpimd_gas/input.xml b/examples/profiling/remd_scpimd_gas/input.xml
index 2fb7345c7..fe7ee2697 100644
--- a/examples/profiling/remd_scpimd_gas/input.xml
+++ b/examples/profiling/remd_scpimd_gas/input.xml
@@ -17,12 +17,12 @@
    </prng>
   <system_template>
      <labels>   [ PREFIX, FILENAME, TEMPERATURE, PRESSURE ] </labels> 
-     <instance> [ REF00, init.xyz, 300, 1 ] </instance>
-     <instance> [ REF01, init.xyz, 310, 1 ] </instance>
-     <instance> [ REF02, init.xyz, 320, 1 ] </instance>
-     <instance> [ REF10, init.xyz, 300, 10 ] </instance>
-     <instance> [ REF11, init.xyz, 310, 10 ] </instance>
-     <instance> [ REF12, init.xyz, 320, 10 ] </instance>
+     <instance> [ REF00, init.xyz, 1, 1 ] </instance>
+     <instance> [ REF01, init.xyz, 2, 1 ] </instance>
+     <instance> [ REF02, init.xyz, 5, 1 ] </instance>
+     <instance> [ REF10, init.xyz, 1, 10 ] </instance>
+     <instance> [ REF11, init.xyz, 2, 10 ] </instance>
+     <instance> [ REF12, init.xyz, 5, 10 ] </instance>
      <template>          
           <system prefix="PREFIX">
             <initialize nbeads="4">
@@ -38,21 +38,21 @@
             <motion mode="dynamics">
               <dynamics mode="scnpt">
                 <barostat mode='sc-isotropic'>
-                  <tau units='femtosecond'> 200</tau>
+                  <tau> 10</tau>
                   <thermostat mode='langevin'>
-                    <tau units='femtosecond'> 100</tau>
+                    <tau> 1</tau>
                   </thermostat>
-                  <h0> [ 32.84024, 0, 0, 0, 32.84024, 0, 0, 0, 32.84024 ]</h0>
+                  <h0> [ 10, 0, 0, 0, 10, 0, 0, 0, 10 ]</h0>
                 </barostat>
                 <thermostat mode='langevin'>
-                  <tau units='femtosecond'> 100 </tau>
+                  <tau> 1 </tau>
                 </thermostat>
-                <timestep units="femtosecond"> 1.0 </timestep>
+                <timestep> 0.01 </timestep>
               </dynamics>
             </motion>
             <ensemble>
-              <temperature units="kelvin"> TEMPERATURE </temperature>
-              <pressure units="megapascal"> PRESSURE </pressure>
+              <temperature> TEMPERATURE </temperature>
+              <pressure> PRESSURE </pressure>
               <bias> <force forcefield="dummy.2"/> </bias>
               <bias_weights> [0.0] </bias_weights>
             </ensemble>
diff --git a/examples/profiling/remd_scpimd_noff/input.xml b/examples/profiling/remd_scpimd_noff/input.xml
index c824bd96e..3ac128633 100644
--- a/examples/profiling/remd_scpimd_noff/input.xml
+++ b/examples/profiling/remd_scpimd_noff/input.xml
@@ -15,12 +15,12 @@
    </prng>
   <system_template>
      <labels>   [ PREFIX, FILENAME, TEMPERATURE, PRESSURE ] </labels> 
-     <instance> [ REF00, init.xyz, 300, 1 ] </instance>
-     <instance> [ REF01, init.xyz, 310, 1 ] </instance>
-     <instance> [ REF02, init.xyz, 320, 1 ] </instance>
-     <instance> [ REF10, init.xyz, 300, 10 ] </instance>
-     <instance> [ REF11, init.xyz, 310, 10 ] </instance>
-     <instance> [ REF12, init.xyz, 320, 10 ] </instance>
+     <instance> [ REF00, init.xyz, 1, 1 ] </instance>
+     <instance> [ REF01, init.xyz, 2, 1 ] </instance>
+     <instance> [ REF02, init.xyz, 5, 1 ] </instance>
+     <instance> [ REF10, init.xyz, 1, 10 ] </instance>
+     <instance> [ REF11, init.xyz, 2, 10 ] </instance>
+     <instance> [ REF12, init.xyz, 5, 10 ] </instance>
      <template>          
           <system prefix="PREFIX">
             <initialize nbeads="4">
@@ -36,21 +36,21 @@
             <motion mode="dynamics">
               <dynamics mode="scnpt">
                 <barostat mode='sc-isotropic'>
-                  <tau units='femtosecond'> 200</tau>
+                  <tau> 10</tau>
                   <thermostat mode='langevin'>
-                    <tau units='femtosecond'> 100</tau>
+                    <tau> 1</tau>
                   </thermostat>
-                  <h0> [ 32.84024, 0, 0, 0, 32.84024, 0, 0, 0, 32.84024 ]</h0>
+                  <h0> [ 10, 0, 0, 0, 10, 0, 0, 0, 10 ]</h0>
                 </barostat>
                 <thermostat mode='langevin'>
-                  <tau units='femtosecond'> 100 </tau>
+                  <tau> 1 </tau>
                 </thermostat>
-                <timestep units="femtosecond"> 1.0 </timestep>
+                <timestep> 0.01 </timestep>
               </dynamics>
             </motion>
             <ensemble>
-              <temperature units="kelvin"> TEMPERATURE </temperature>
-              <pressure units="megapascal"> PRESSURE </pressure>
+              <temperature> TEMPERATURE </temperature>
+              <pressure> PRESSURE </pressure>
               <bias> <force forcefield="dummy.2"/> </bias>
               <bias_weights> [0.0] </bias_weights>
             </ensemble>
diff --git a/ipi/engine/forces.py b/ipi/engine/forces.py
index 0b123123d..7a67bd2e4 100644
--- a/ipi/engine/forces.py
+++ b/ipi/engine/forces.py
@@ -1120,10 +1120,14 @@ def forcesvirs_4th_order(self, index):
         # calculates the finite displacement.
         fbase = dstrip(self.f)
         eps = self.mforces[index].epsilon
-        delta = np.abs(eps) / np.sqrt(
-            (fbase / self.beads.m3 * fbase / self.beads.m3).sum()
-            / (self.nbeads * self.natoms)
-        )
+        foverm = np.sqrt((fbase / self.beads.m3 * fbase / self.beads.m3).sum() / (self.nbeads * self.natoms))
+        if np.abs(foverm)>1e-20:
+            delta = np.abs(eps) / np.sqrt(
+                (fbase / self.beads.m3 * fbase / self.beads.m3).sum()
+                / (self.nbeads * self.natoms)
+            )
+        else: # defaults to eps if otherwise we'd get an 1/0
+            delta = np.abs(eps) 
         dq = delta * fbase / self.beads.m3
 
         # stores the force component.
@@ -1141,11 +1145,11 @@ def forcesvirs_4th_order(self, index):
             # This makes the FD incorrect. A fix could be to ALWAYS compute the odd and the even
             # beads on two different ring polymers of P / 2 beads, but centered difference + RPC
             # seems to be a neater solution. Anyway, the cost of a CNT-DIFF in comparison to a FWD-DIFF
-            # is marginal in when RPC + MTS is used.
+            # is marginal when RPC + MTS is used.
 
             if self.mforces[index].nbeads != self.nbeads:
                 warning(
-                    "ERROR: high order PIMD + RPC works with a centered finite difference only! (Uness you find an elegant solution :))"
+                    "ERROR: high order PIMD + RPC works with a centered finite difference only! (Unless you find an elegant solution :))"
                 )
                 exit()