Skip to content

Harmonic Sampling#

Warnings

  • Never ever rely on estimates obtained by the following method without checking for the relevance of anharmonic effects (see next tutorial).
  • The tutorial assumes you are familiar with performing phonon calculations.

Thermodynamic sampling via MD simulations#

As outlined in the MD tutorial, MD simulations provide a way to sample phase space in order to obtain thermodynamic expectation values of observables,

\begin{align} \langle O \rangle = \lim_{t_0 \rightarrow \infty} \frac{1}{t_0} \int_0^{t_0} {\rm d} t ~ O \left( \mathbf{R} (t), {\bf P} (t)\right) \label{eq:O2}~, \end{align}

where O \left( \mathbf{R} (t), {\bf P} (t) \right) denotes the instantaneous value of the observable O obtained for the phase-space configuration \{{\bf R} (t), {\bf P} (t) \} evaluated with respect to the the full many-body Hamiltonian \mathcal H ({\bf P}, {\bf R}). For example, O could be the instantaneous pressure as used in the MD tutorial.

Circumventing the dynamical propagation by generating samples from a harmonic distribution#

In MD simulations, the many-body potential \mathcal V needs to be evaluated both for propagating the dynamical evolution of the system, and for evaluating the observable along the trajectory. On a conceptual level, only the evaluation of the observable is necessary to provide an estimate of the thermodynamic expectation value if representative samples can be generated by other means.

Following [West2006], we now present a way to create these samples based on the harmonic approximation. In this approximation, the potential is given by

\begin{align} \mathcal V^{(2)} ({\bf R}) = \frac{1}{2} \sum_{I, J} \Delta {\bf R}_I \cdot \Phi_{IJ} \Delta {\bf R}_J~, \label{eq:V2} \end{align}

with the 3 \times 3 force constant matrices \Phi_{IJ} and atomic displacements \Delta {\bf R}_I. The Newton equations of motion therefore read

\begin{align} \left( \begin{array}{c} M_{1} \Delta \ddot{\mathbf{R}}_{1} \\ M_{2} \Delta \mathbf{R}_{2} \\ \vdots \\ M_{N} \Delta \ddot{\mathbf{R}}_{N} \end{array}\right)=\left(\begin{array}{cccc} \Phi_{11} & \Phi_{12} & \cdots & \Phi_{1 N} \\ \Phi_{21} & \Phi_{22} & \cdots & \Phi_{2 N} \\ \vdots & \vdots & \ddots & \vdots \\ \Phi_{N 1} & \Phi_{N 2} & \cdots & \Phi_{N N} \end{array}\right)\left(\begin{array}{c} \Delta \mathbf{R}_{1} \\ \Delta \mathbf{R}_{2} \\ \vdots \\ \Delta \mathbf{R}_{N} \end{array}\right)~, \end{align}

and can be solved analytically, yielding

\begin{align} \Delta {\bf R}_I (t) = \frac{1}{\sqrt{M_I}} \sum_s {\bf e}_{sI} A_s \sin (\omega_s t + \phi_s)~, \label{eq:R2(t)} \end{align}

where {\bf e}_{sI} denotes the eigenvector of the dynamical matrix D_{IJ} = \Phi_{IJ} / \sqrt{M_I M_J} corresponding to atom I and phonon mode s and \omega_s the respective frequency, cf. the phonon tutorial. The amplitudes A_s and phases \phi_s are fixed by the initial condition \{{\bf R} (t_0), {\bf P} (t_0) \}, where in thermal equilibrium [Dove1993]

\begin{align} \left\langle A_{s}\right\rangle = \sqrt{ \frac{\hbar}{\omega_{s}}\left(n_{\mathrm{B}}(\omega_s, T)+\frac{1}{2}\right) } ~\stackrel{k_{\rm B} T \,\gg\, \hbar \omega_s}{\longrightarrow}~ \frac{\sqrt{2 k_{\rm B} T}}{\omega_s}~. \label{eq:As} \end{align}

Noting that \sin (\omega_s t + \phi_s) is essentially a random number in [-1, 1], and mimicking thermal fluctuations, samples can be generated by

\begin{align} \Delta {\bf R}_{I} (n) =\frac{1}{\sqrt{M_{I}}} \sum_{s} \zeta_{s} (n) \left\langle A_{s}\right\rangle {\bf e}_{s I}~, \label{eq:samples} \end{align}

where \langle A_s \rangle is given by Eq. \eqref{eq:As} and each \zeta_s (n) is a normally distributed random number, where n labels the sample.

Example: LDA-Si at 300K#

Obtain force constants#

We will re-use the calculation from the previous tutorial on phonon calculations in the conventional cell with 8 atoms. If you didn't run these calculations, you find the respective calculations in our reference repository.

As start, perform the postprocess in the phonopy folder:

cd phonopy
vibes output phonopy -bs

Check the bandstructure in output/bandstructure.pdf for plausibility.

Remap the force constants to the supercell#

The force constants are written to output/FORCE_CONSTANTS which is a condensed representation in (N_{\rm prim}, N, 3, 3) shape, where N_{\rm prim} is the number of atoms in the primitive cell and N is the number of atoms in the supercell. For creating samples, they need to be mapped to a full 3 N \times 3N shape. This can be done with the CLI tool vibes utils fc remap,

cd output
vibes utils fc remap

which will create FORCE_CONSTANTS_remapped as a plain, numpy readable text file.

The current folder …/phonopy/output should now contain the following files:

>>> ls
bandstructure.pdf  band.yaml  FORCE_CONSTANTS  FORCE_CONSTANTS_remapped  geometry.in.primitive  geometry.in.supercell

We are now ready to create samples according to Eq. \eqref{eq:samples}.

Create samples#

Samples can be created with the CLI tool utils create-samples. We now create 10 samples according to Eq. \eqref{eq:samples} at a temperature of 300\,{\rm K} by running

vibes utils create-samples geometry.in.supercell -fc FORCE_CONSTANTS_remapped -n 10 -T 300

The geometry files are written as geometry.in.supercell.0300K.??? to the current directory. Let's move them to a folder called samples_20K:

mkdir samples_300K
mv geometry.in.supercell.* samples_300K

Create a new working directory where you

  • copy the files geometry.in.primitive and geometry.in.supercell,
  • copy your samples-folder samples_300K.

Move to that directory. Remark: copying the geometry.in.primitive is optional. However, you should always know from where you started, it might save some hours of unnecesary work at some point :-).

Compute the samples#

In order to compute pressure in all samples, prepare an aims.in like this:

[files]
geometries:                    samples_300K/geometry.in.*
primitive:                     geometry.in.primitive
supercell:                     geometry.in.supercell

[calculator]
name:                          aims
socketio:                      true
workdir:                       aims_300K

[calculator.parameters]
xc:                            pw-lda
compute_analytical_stress:     true

[calculator.kpoints]
density:                       2

[calculator.basissets]
default:                       light

This will run a calculation for all structures found in samples_300K/geometry.in.* as discussed earlier in the tutorial on singlepoint calculations, and write the results to a trajectory.son file in the directory aims_300K.

Run the calculation with

vibes run singlepoint aims.in | tee log.aims

This will take a few minutes depending on your machine.

Compare sampling vs. MD#

Similar to MD you can create a trajectory dataset from trajectory.son by running vibes output md:

cd  aims_300K
vibes output md

Since this is not an MD, the time axis in the dataset will simply label the samples instead of corresponding to an actual simulation time.

The dataset in trajectory.nc can be inspected similar to the MD case. For example, run

import xarray as xr
from ase.units import GPa


ds = xr.load_dataset("trajectory.nc")

p = ds.pressure_potential.to_series() / GPa

ax = p.plot(marker="x", lw=0)

p.expanding().mean().plot(ax=ax, color="k")

ax.set_xlabel('Sample number')
Pressure plot

image

The mean pressure and standard error1 are

\begin{align*} \langle p_{\rm Pot} (300\,{\rm K}) \rangle^{(2)} = (-0.11 \pm 0.05)\,{\rm GPa}~, \end{align*}

which is converged with a precision of \sim 20\,\%. However, the superscript \langle \cdot \rangle^{(2)} reminds us that we used the harmonic approximation to create the samples. Indeed, the pressure found by harmonic sampling is about 40\,\% larger than the MD reference computed in the previous tutorial. Yet, the two values coincide within their error margins. Further sampling on both sides would be needed to compute the actual difference between the two sampling techniques.

Take Home Messages#

The sampling trick can be used to rapidly estimate the expectation value of a static observable and typically converges 1-2 orders of magnitude faster than a MD simulation2.

The sampling trick only works if the harmonic approximation is a good reference for the system of interest. This might or might not be the case and is difficult to tell a priori. We present a scheme to estimate the importance of anharmonic effects in the next chapter of the tutorial.


  1. The samples are uncorrelated by construction, so \tau = 1 and \tilde N_t = N_t in this case. 

  2. Remember that the correlation length was about 200 steps in the MD case while the samples created from the harmonic distribution are uncorrelated per construction.