More explanation

examples/heat-capacity/data/input.xml

@@ -2,7 +2,7 @@
   <output prefix='simulation'>
       <properties stride='2' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, scaledcoords(fd_delta=5e-3) ] </properties>
   </output>
-  <total_steps> 5000 </total_steps>
+  <total_steps> 2000 </total_steps>
   <prng>
     <seed> 32342 </seed>
   </prng>

examples/heat-capacity/heat-capacity.py

@@ -36,13 +36,17 @@
 # ----------------------------
 #
 # This follows the same steps as the ``path-integrals`` example. One important
-# difference is that we will request the ``scaledcoords`` output, which
+# difference is that we will request the ``scaledcoords`` output to the relevant
+# section of the ``i-PI`` input XML file, which
 # contains estimators that can be used to calculate the total energy and
 # heat capacity as defined in this `paper <https://arxiv.org/abs/physics/0505109>`_.
-# In order to do this, the ``scaledcoords`` output is added to the relevant section
-# of the ``i-PI`` input XML file.
+#
+# The input file is shown below. It should be noted that the ``scaledcoords``
+# is given a finite differences displacement as a parameter. This is necessary
+# as the estimators require higher order derivatives of the potential energy,
+# which are calculated using finite differences.
 
-# Open and read the XML file
+# Open and show the XML file
 with open("data/input.xml", "r") as file:
     xml_content = file.read()
 print(xml_content)
@@ -90,21 +94,31 @@
 # <documentation https://ipi-code.org/i-pi/output-tags.html>,
 # the two quantities returned by the ``scaledcoords`` output are ``eps_v``
 # and ``eps_v'``, defined in the aforementioned
-# `paper <https://arxiv.org/abs/physics/0505109>`_. The same paper contains
-# the formulas to calculate the total energy and heat capacity from these
-# estimators.
+# `paper <https://arxiv.org/abs/physics/0505109>`_.
+# 
+# These estimators (:math:`\eps_v` and :math:`\eps_v'`) are derived in the
+# "scaled coordinates" formalism, which is a useful trick to avoid the
+# growth of the error in the instantaneous values of the estimators with
+# the number of beads used in the path integral simulation.
+# 
+# The same paper contains the formulas to calculate the total energy and
+# heat capacity from these estimators:
+# 
+# .. math::
+#   E = \eps_v
+#   C_V = k_B \beta^2 \left( \langle \eps_v^2 \rangle - \langle
+#       \eps_v \rangle^2 - \langle \eps_v' \rangle \right)
 #
 # First the energy:
 
-eps_v = np.loadtxt("simulation.out")[100:, 6]
-eps_v_prime = np.loadtxt("simulation.out")[100:, 7]
-# discarding the first 100 steps which are highly non-equilibrated
+eps_v = np.loadtxt("simulation.out")[:, 6]
+eps_v_prime = np.loadtxt("simulation.out")[:, 7]
 
 energy_estimator = eps_v  # from the paper
 
 fix, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True)
 ax.plot(
-    output_data["time"][100:],
+    output_data["time"],
     energy_estimator,
     "b",
     label="Total energy$",