Small tweaks

examples/heat-capacity/README.rst

@@ -1,7 +1,7 @@
-Running, analizing and visualizing a path integral MD simulation
-================================================================
+Heat capacity of water from path integral molecular dynamics
+============================================================
 
-This example shows how to run a path integral molecular dynamics 
-simulation using ``i-PI``, analyze the output and visualize the 
-trajectory in ``chemiscope``. It uses `LAMMPS <http://lammps.org>`_
-as the driver to simulate the q-TIP4P/f water model. 
+This example shows how to calculate the constant-volume heat
+capacity from a path integral molecular dynamics simulation 
+using ``i-PI``. `LAMMPS <http://lammps.org>`_ is used
+as the driver for the evaluation of the q-TIP4P/f water model.

examples/heat-capacity/data/input.xml

@@ -1,9 +1,6 @@
 <simulation verbosity='medium' safe_stride='100'>
   <output prefix='simulation'>
     <properties stride='1' filename='out'>  [ step, time{picosecond}, conserved{electronvolt}, temperature{kelvin}, kinetic_cv{electronvolt}, potential{electronvolt}, scaledcoords ] </properties>
-    <trajectory filename='pos' stride='20'> positions </trajectory>
-    <trajectory filename='kin' stride='20'> kinetic_cv </trajectory>
-    <trajectory filename='kod' stride='20'> kinetic_od </trajectory>
   </output>
   <total_steps> 200 </total_steps>
   <prng>

examples/heat-capacity/heat-capacity.py

@@ -1,6 +1,6 @@
 """
-Path integral molecular dynamics
-================================
+Quantum heat capacity of water
+==============================
 
 :Authors: Michele Ceriotti `@ceriottm <https://github.com/ceriottm/>`_
 
@@ -21,8 +21,8 @@
 
 
 # %%
-# Quantum heat capacity of water
-# ------------------------------
+# A non-trivial estimator
+# -----------------------
 #
 # As introduced in the ``path-integrals`` example, path-integral estimators
 # for observables that depend on momenta are generally not trivial to compute.
@@ -89,6 +89,7 @@
 ax.set_xlabel(r"$t$ / ps")
 ax.set_ylabel(r"energy / eV")
 ax.legend()
+plt.show()
 
 # %%
 # We will now plot the values of the energy and heat capacity estimators
@@ -102,21 +103,23 @@
 #
 # First the energy:
 
-eps_v = np.loadtxt("simulation.out")[:, 6]
-eps_v_prime = np.loadtxt("simulation.out")[:, 7]
+eps_v = np.loadtxt("simulation.out")[20:, 6]
+eps_v_prime = np.loadtxt("simulation.out")[20:, 7]
+# discarding the first 20 steps which are highly non-equilibrated
 
 energy_estimator = eps_v  # from the paper
 
 fix, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True)
 ax.plot(
-    output_data["time"],
+    output_data["time"][20:],
     energy_estimator,
     "b",
     label="Total energy$",
 )
 ax.set_xlabel(r"$t$ / ps")
 ax.set_ylabel(r"$E / a.u.$")
 ax.legend()
+plt.show()
 
 # %%
 # And, finally, the heat capacity: