From 77597cf1001585872f6099484d97dd5a0040f6ef Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Sat, 21 Sep 2024 16:09:41 +0200 Subject: [PATCH] Updated heat capacity example to my satisfaction --- examples/heat-capacity/data/in.lmp | 2 +- examples/heat-capacity/data/input.xml | 7 +- examples/heat-capacity/heat-capacity.py | 181 ++++++++++++------- examples/pi-metad/pi-metad.py | 2 +- src/__pycache__/get_examples.cpython-312.pyc | Bin 1506 -> 0 bytes 5 files changed, 125 insertions(+), 67 deletions(-) delete mode 100644 src/__pycache__/get_examples.cpython-312.pyc diff --git a/examples/heat-capacity/data/in.lmp b/examples/heat-capacity/data/in.lmp index 3e3667f2..61e74927 100644 --- a/examples/heat-capacity/data/in.lmp +++ b/examples/heat-capacity/data/in.lmp @@ -3,7 +3,7 @@ atom_style full pair_style lj/cut/tip4p/long 1 2 1 1 0.278072379 17.001 # high-pppm precision and shift to get meaningful fd estimates -kspace_style pppm/tip4p 1e-6 +kspace_style pppm/tip4p 1e-5 pair_modify shift yes bond_style class2 angle_style harmonic diff --git a/examples/heat-capacity/data/input.xml b/examples/heat-capacity/data/input.xml index 8a9c2067..ae6b8209 100644 --- a/examples/heat-capacity/data/input.xml +++ b/examples/heat-capacity/data/input.xml @@ -1,6 +1,9 @@ - - [ step, time{picosecond}, conserved, potential, kinetic_cv, scaledcoords(fd_delta=5e-3) ] + + + [ step, time{picosecond}, conserved, potential, kinetic_cv, + scaledcoords(fd_delta=5e-3) ] + 2000 diff --git a/examples/heat-capacity/heat-capacity.py b/examples/heat-capacity/heat-capacity.py index 0c36ca04..e0174ea2 100644 --- a/examples/heat-capacity/heat-capacity.py +++ b/examples/heat-capacity/heat-capacity.py @@ -2,7 +2,9 @@ Quantum heat capacity of water ============================== -:Authors: Michele Ceriotti `@ceriottm `_ +:Authors: +Filippo Bigi `@frostedoyster `_; +Michele Ceriotti `@ceriottm `_ This example shows how to estimate the heat capacity of liquid water from a path integral molecular dynamics simulation. The dynamics are @@ -12,8 +14,10 @@ model `_. """ +import os import subprocess import time +import xml.etree.ElementTree as ET import ipi import matplotlib.pyplot as plt @@ -21,55 +25,79 @@ # %% -# A non-trivial estimator -# ----------------------- +# A non-trivial energy estimator +# ------------------------------ # -# As introduced in the ``path-integrals`` example, path-integral estimators +# As introduced in the +# `path-integrals example +# `_, +# path-integral estimators # for observables that depend on momenta are generally not trivial to compute. # -# In this example, we will focus on the heat capacity, which is one such +# In this example, we will focus on the constant-volume heat capacity, +# :math:`c_V`, which is one such # observable, and we will calculate it for liquid water at room temperature. +# Because of the presence of high-frequency vibrations, many of the nuclear +# degrees of freedom are trapped in the vibrational ground state, which reduces +# substantially the heat capacity from the value that would be obtained +# in the classical limit. See `this review +# `_ +# for an overview of the impact of quantum nuclei on the properties of water. +# From a computational perspective, this means that it is necessary to use +# specialized simulations and estimators to evaluate the correct value. # %% # Running the PIMD calculation -# ---------------------------- +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # This follows the same steps as the ``path-integrals`` example. One important # 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 `_. +# heat capacity as following +# `Yamamoto, J. Chem. Phys. (2005) `_. # -# The input file is shown below. It should be noted that the ``scaledcoords`` +# The input file is shown below. It should be noted that ``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 show the XML file -with open("data/input.xml", "r") as file: - xml_content = file.read() -print(xml_content) +# which are calculated using finite differences. This also means that +# evaluating the estimator adds substantial overhead (so it is wise to only +# compute it every few simulation steps, to eliminate correlations between +# snapshots) and that one should be careful to use well-converged simulation +# parameters to avoid discontinuities and noise (for instance, we increase +# the accuracy of the particle-mesh electrostatic calculation, and use a +# shifted Lennard-Jones potential to avoid a discontinuity at the cutoff). + +# Open and show the relevant part of the input +xmlroot = ET.parse("data/input.xml").getroot() +print(" " + ET.tostring(xmlroot.find(".//properties"), encoding="unicode")) # %% # We launch the i-PI and LAMMPS processes, exactly as in the # ``path-integrals`` example. -ipi_process = subprocess.Popen(["i-pi", "data/input.xml"]) -time.sleep(2) # wait for i-PI to start -lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(2)] - -ipi_process.wait() -lmp_process[0].wait() -lmp_process[1].wait() - -output_data, output_desc = ipi.read_output("simulation.out") +# don't rerun if the outputs already exist +ipi_process = None +if not os.path.exists("water-cv.out"): + ipi_process = subprocess.Popen(["i-pi", "data/input.xml"]) + time.sleep(2) # wait for i-PI to start + lmp_process = [subprocess.Popen(["lmp", "-in", "data/in.lmp"]) for i in range(2)] +# %% +# Skip this cell if you want to run in the background +if ipi_process is not None: + ipi_process.wait() + lmp_process[0].wait() + lmp_process[1].wait() # %% +# Analyzing the results +# ~~~~~~~~~~~~~~~~~~~~~ # Let's plot the potential and conserved energy as a function of time, # just to check that the simulation ran sensibly. +output_data, output_desc = ipi.read_output("water-cv.out") fix, ax = plt.subplots(1, 1, figsize=(4, 3), constrained_layout=True) ax.plot( output_data["time"], @@ -84,19 +112,18 @@ label="Conserved, $H$", ) ax.set_xlabel(r"$t$ / ps") -ax.set_ylabel(r"energy / eV") +ax.set_ylabel(r"energy / a.u.") ax.legend() plt.show() -plt.close() # %% -# As described in the ``i-PI`` -# , +# As described in the `i-PI documentation +# `_, # the two quantities returned by the ``scaledcoords`` output are ``eps_v`` # and ``eps_v'``, defined in the aforementioned # `paper `_. # -# These estimators (:math:`\eps_v` and :math:`\eps_v'`) are derived in the +# These estimators (:math:`\epsilon_v` and :math:`\epsilon_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. @@ -105,9 +132,9 @@ # heat capacity from these estimators: # # .. math:: -# E = \langle \eps_v \rangle -# C_V = k_B \beta^2 \left( \langle \eps_v^2 \rangle - \langle -# \eps_v \rangle^2 - \langle \eps_v' \rangle \right) +# E = \langle \epsilon_v \rangle \quad +# C_V = k_B \beta^2 \left( \langle \epsilon_v^2 \rangle - \langle +# \epsilon_v \rangle^2 - \langle \epsilon_v' \rangle \right) # # First the energy, whose estimator will be compared to the total energy # calculated as the sum of the potential and kinetic energy estimators. @@ -115,8 +142,8 @@ # estimator (the "centroid virial" estimator), the two total energies are # the same. -eps_v = np.loadtxt("simulation.out")[:, 5] -eps_v_prime = np.loadtxt("simulation.out")[:, 6] +eps_v = output_data["scaledcoords(fd_delta=5e-3)"][:, 0] +eps_v_prime = output_data["scaledcoords(fd_delta=5e-3)"][:, 1] energy_estimator = eps_v # first formula @@ -124,7 +151,8 @@ ax.plot( output_data["time"], energy_estimator - energy_estimator[0], - label="Total energy (scaled coordinates estimator)", + "b-", + label="scaled coordinates estimator", ) ax.plot( output_data["time"][:], @@ -134,35 +162,51 @@ + output_data["kinetic_cv"] - output_data["kinetic_cv"][0] ), - label="Total energy (potential + kinetic)", + "r.", + label="potential + virial kinetic", ) ax.set_xlabel(r"$t$ / ps") -ax.set_ylabel(r"$E / a.u.$") +ax.set_ylabel(r"total energy / a.u.") ax.legend() plt.show() # %% -# And, finally, the heat capacity: +# And, finally, the heat capacity. Note that normally the simulation +# requires a few ps for equilibration. Here we discard a few dozen steps +# to eliminate the initial jump, which is due to the relaxation of the +# ring polymers starting from a single atomic configuration. # i-PI scaledcoords outputs are in atomic units (see docs) kB = 3.16681e-6 # Boltzmann constant in atomic units T = 298.0 # temperature in K, as defined in the input file beta = 1.0 / (kB * T) +skip = 20 heat_capacity = ( # second formula - kB * (beta**2) * (np.mean(eps_v**2) - np.mean(eps_v) ** 2 - np.mean(eps_v_prime)) + kB + * (beta**2) + * ( + np.mean(eps_v[skip:] ** 2) + - np.mean(eps_v[skip:]) ** 2 + - np.mean(eps_v_prime[skip:]) + ) ) heat_capacity_per_molecule = heat_capacity / 32 # 32 molecules in the simulation print(f"Heat capacity (per water molecule): {(heat_capacity_per_molecule/kB):.2f} kB") +# %% +# You may recognize that the first part of the estimator is reminiscent +# of the classical estimator for the heat capacity as the fluctuations of the +# (quantum) total energy, which in this case however requires a correction given +# by the mean of the second part of the scaled-coordinates estimator. # %% -# Finding an error estimate -# ------------------------- +# Estimating the statistical error +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # # Especially with such an underconverged simulation, it is important to -# estimate the error in the heat capacity. +# estimate the statistical error in the heat capacity. # # Generally, errors on measurements are computed # as "standard errors", i.e. the standard deviation of a series of data points @@ -187,32 +231,38 @@ def autocorrelation_time(x): # %% -# We can now calculate the error on the heat capacity estimate. +# Using these helper functions, we can now calculate the error on the various +# parts of the heat capacity estimator. Note also the autocorrelation times, that +# are just a little larger than one, indicating that the stride used to print out +# the estimators is appropriate (as there is little correlation between the samples). # Autocorrelation times (i.e. number of steps needed to have independent samples) -autocorr_time_error_eps_v = autocorrelation_time(eps_v) -autocorr_time_error_eps_v_squared = autocorrelation_time(eps_v**2) -autocorr_time_error_eps_v_prime = autocorrelation_time(eps_v_prime) +autocorr_time_error_delta_eps_v = autocorrelation_time( + (eps_v[skip:] - eps_v[skip:].mean()) ** 2 +) +autocorr_time_error_eps_v_prime = autocorrelation_time(eps_v_prime[skip:]) + +print( + f""" +Autocorrelation times (in number of samples): +(eps-)^2: {autocorr_time_error_delta_eps_v:.2f} +eps': {autocorr_time_error_eps_v_prime:.2f} +""" +) # Effective number of samples -effective_samples_eps_v = len(eps_v) / autocorr_time_error_eps_v -effective_samples_eps_v_squared = len(eps_v) / autocorr_time_error_eps_v_squared -effective_samples_eps_v_prime = len(eps_v) / autocorr_time_error_eps_v_prime +effective_samples_delta_eps_v = len(eps_v[skip:]) / autocorr_time_error_delta_eps_v +effective_samples_eps_v_prime = len(eps_v[skip:]) / autocorr_time_error_eps_v_prime # Standard errors using the effective number of samples -error_eps_v = np.std(eps_v) / np.sqrt(effective_samples_eps_v) -error_eps_v_squared = np.std(eps_v**2) / np.sqrt(effective_samples_eps_v_squared) -error_eps_v_prime = np.std(eps_v_prime) / np.sqrt(effective_samples_eps_v_prime) +error_delta_eps_v = np.std((eps_v[skip:] - eps_v[skip:].mean()) ** 2) / np.sqrt( + effective_samples_delta_eps_v +) +error_eps_v_prime = np.std(eps_v_prime[skip:]) / np.sqrt(effective_samples_eps_v_prime) -# Error on the heat capacity +# Error on the heat capacity (assuming quadrature sum) error_heat_capacity = ( - kB - * (beta**2) - * np.sqrt( # error propagation in the sum of terms - error_eps_v_squared**2 - + (2.0 * np.mean(eps_v) * error_eps_v) ** 2 # error of **2 - + error_eps_v_prime**2 - ) + kB * (beta**2) * np.sqrt(error_delta_eps_v**2 + error_eps_v_prime**2) ) error_heat_capacity_per_molecule = ( @@ -227,7 +277,10 @@ def autocorrelation_time(x): # %% # The obtained heat capacity is consistent with the values from the literature -# (see ...). +# (see e.g. `Ceriotti et al., J. Chem. Phys. (2011) +# `_ +# where the convergence of the heat capacity with number of beads is shown for the same +# water model used in this example). # However, the error is quite large, which is expected given the short simulation time. # To reduce the error, one would need to run a longer simulation. Other important error # sources, which are not accounted for in the error estimate, are the finite size of the @@ -235,7 +288,9 @@ def autocorrelation_time(x): # reliable results. # # In a realistic simulation, up to a few 100s of picoseconds might be needed to reduce -# the sampling error to a small value (1-10% of the heat capacity). For water at room +# the sampling error to a small value (1-2% of the heat capacity). For water at room # temperature, you will need 32 beads at the very least (8 were used in this example). -# It is more difficult to give a general rule for the system size, but a few hundred -# water molecules would be a reasonable guess (32 were used in this example). +# It is more difficult to give a general rule for the system size: (quantum) energy +# fluctuations are usually localized, but to guarantee accurate sampling of the +# liquid structure, a few hundred water molecules would be a reasonable guess +# (32 were used in this example). diff --git a/examples/pi-metad/pi-metad.py b/examples/pi-metad/pi-metad.py index 1dc87f51..52a17382 100644 --- a/examples/pi-metad/pi-metad.py +++ b/examples/pi-metad/pi-metad.py @@ -467,7 +467,7 @@ def moving_average(arr, window_size): # ----------------------------------- # # You can see `this recipe -# `_ +# `_ # for a brief introduction to path integral simulations with `i-PI`. # From a practical perspective, very little needs to change with respect # to the classical case. diff --git a/src/__pycache__/get_examples.cpython-312.pyc b/src/__pycache__/get_examples.cpython-312.pyc deleted file mode 100644 index bb90b39b9fb0e31bd0e6f08ce79e61ed0d411f03..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1506 zcmb7EUuYaf7@wV;+uh6EHCVx13c@^C>8i#IvNXPogjSbR;6!2XSU|n={ZLe(d*Wzxigr zZ+_o5`)w{a0tEei?f2?C5P)CgOG5g<*nM09;0AC(8MulI8%kM;H7vtet7Vm_btBm} z)b-5JPOL*BbIJNp9}+n_j7ThEY={aax$1^Eg`~t&P9O4}2W3pVDA`{gaZ_)Da>mt! z%#xg|FP;E^9f5#I&Tj0b+0UbZWw;DYqg^aH%xiMYDRvmKIl-@%=FZ@dGcTy&u;p>P z(jd4>8MblF;|r|{wjJ)#0Q@3a4H7zobg`09tAclb5aI^l zgTn(*kCzV$J^%!Ok3{=EfeA1#fi7@hP4JsD09*kaCHv}HVhwdYJ+u?+KXyEz`IqoMYR$~+Qcfu`$)i)~nV!qg&EYtTw$%w8zIe)hH5QUqxjp{C6jBBNr5B(QxF83AR! z^qES8s-D1H!iw*oEG$r;6okiXHfh13E-BcY`l2MynQ$mwtOzcI%qi3ew~`Vko6Ax5 z-wq}RE-N@ZE5r);LjyA}Dqml`eet^!-yGX4y>aJS?}hn4RhZT8J$v#x?CKqT+sJPj z`HzpTkFAYukDFWL<`<9k#tV0ir*~vrAHbCM*A5V7D0{YEREx*afl)~bcRVVmg>hh1 z_Qygj49HcujE^aqZ%;{DRg!Iarue_dxt$$gG9hs`)d|Y@ zTb4^5%VP34O#U&L{FyO%aUw(9%0}FLRHt4L>0)m7Ccv^}-z_Ub=jP_lGYM&2UtESv zAj8CRswwU=XAjABnL950zA^1nx78ps>`~#9!xOfJ=&OpN{0v_C31s3Y-PeFI^6t`G zOWXR`mOi!$*Y!1hJ-e3u^xUoaoAbA{?qzb_sm@e)x--3s sR_nctxnk^SNYQqVK;`-V2te8Q+TW|=KcaN!$@eGjqI_Qk=vZ9aKa2iGf&c&j