Sampling dependent systematic errors in effective harmonic models
Erki Metsanurk Mattias Klintenberg



 erki.metsanurk@physics.uu.se Department of Physics and Astronomy, Uppsala University, Box 516, S-75120 Uppsala, Sweden 
 Department of Physics and Astronomy, Uppsala University, Box 516, S-75120 Uppsala, Sweden




Effective harmonic methods allow for calculating temperature dependent 
phonon frequencies by incorporating the anharmonic contributions into 
an effective harmonic Hamiltonian. The systematic errors arising from 
such an approximation are explained theoretically and quantified by 
density functional theory based numerical simulations. Two techniques 
with different approaches for sampling the finite temperature phase 
space in order to generate the force-displacement data are compared. It is 
shown that the error in free energy obtained by using either can exceed 
that obtained from 0 K harmonic lattice dynamics analysis which 
neglects the anharmonic effects.



Sampling dependent systematic errors in effective harmonic models
    M. Klintenberg
    December 30, 2023
=================================================================







§ INTRODUCTION

The harmonic approximation of the potential energy is a fundamental 
part of the analysis of vibrational properties of materials both in 
theoretical and experimental studies <cit.>. Its usefulness 
is multifold in the sense that not only does it allow for capturing the 
dominant part of the potential energy surface, but makes possible to 
derive exact analytical relations for all temperature dependent 
vibrational thermodynamic properties. In addition, the harmonic system 
can be used as a reference point for either integration- or 
perturbation-based methods in order to find the corrections due to 
anharmonicity that always exists in real systems.

The lattice dynamics approach <cit.> to calculate the normal 
vibrational modes of a crystal can break down when either very strong 
quantum mechanical effects dominate as is the case for crystalline 
helium <cit.> or more commonly due to mechanical 
instabilities, for example cubic zirconia <cit.> or 
titanium, zirconium and hafnium in BCC structure <cit.>. 
This results in some of the normal modes having imaginary vibrational 
frequencies due to decrease in the total potential energy when the 
atoms are displaced along those modes. Since the harmonic free energy 
and all the other derived thermodynamic properties are defined as 
integrals over the whole normal mode space, it would result in 
complex-valued properties and are, thus, commonly interpreted as 
undefined.

Recently, methods have been proposed to calculate free energies of such 
systems both for cases where the lattice is dynamically unstable 
<cit.> or stabilized due to temperature 
<cit.>. Both are based on partitioning of atomic 
configuration space, calculating the free energy of each region 
separately and adding those together consistent with equilibrium 
thermodynamics.

Instead of the aforementioned detailed mapping it is also possible to 
fit high temperature force-displacement data obtained from atomistic 
simulation to a truncated expansion of the potential energy surface 
after which the second order effective force constants can be used to 
calculate the temperature dependent phonon frequencies and, hence, the 
effective harmonic free energy. This is the basis for the methods of 
self-consistent ab-initio lattice dynamics (SCAILD) 
<cit.> and 
temperature-dependent effective potential (TDEP) 
<cit.>. Neither are inherently 
limited to be used for analyzing only systems with dynamical 
instabilities, which could be considered as the most severe case, but 
any kind of anharmonicity.

Currently the best way to gather the data for either of the methods is 
through density functional theory (DFT) <cit.> which in 
principle provides a parameter-free way to explore the dynamics of 
realistic, as opposed to model, systems. Ab-initio based 
thermodynamics has been shown to be a promising way to predict the 
thermal properties of materials taking all the relevant excitations, 
not limited to vibrational, into account 
<cit.>.

Without considering whether DFT with its approximations can provide 
reliable enough data, there are also other systematic and statistical 
uncertainties associated with calculating thermodynamic properties. 
Some of those, such as the limited system size, can be partially dealt 
by Fourier interpolating the information about vibrations in an 
infinite crystal <cit.>, while others such as limited 
amount of time for sampling can be dealt with by upsampled 
thermodynamic integration with Langevin dynamics (UP-TILD) or 
harmonically mapped averaging <cit.>. 
Regardless of the method used, the total uncertainty required to 
produce satisfactory results is often under 1 meV/atom. For example, a 
6 meV/atom shift on the energy scale results in 400 K (60%) 
overestimation of the FCC to BCC transition temperature for calcium 
using ab-initio calculations <cit.>.

Given the requirement of high accuracy for any method used to calculate 
thermodynamic properties it is necessary to understand how the models 
behave under certain conditions. In this study we analyze the accuracy 
of effective harmonic models by comparing SCAILD and TDEP both 
theoretically and numerically at different temperatures. It has been 
shown previously that the difference between the free energies 
predicted by the two methods increases with temperature 
<cit.>, however, no theoretical explanation was provided 
what in particular could cause the discrepancy between these two very 
similar models.



§ METHODOLOGY



 §.§ Harmonic approximation at 0 K and finite temperatures




In order to simplify the comparison of the two methods a slight 
reformulation of the theory is needed compared to that described in the 
original works. The starting point is the force-displacement relation 
through a 3N × 3N force constant matrix Φ, so for any 
atomic configuration c the forces are given by


    f⃗_⃗c⃗ = - Φ(θ⃗) u⃗_⃗c⃗


In a similar fashion to Hellman  <cit.>, the 
elements of Φ are expressed as linear combinations of parameters 
θ⃗ as the force constants of crystalline systems are in general 
not independent, but constrained by rotational, 
translational and inversion symmetry of the lattice and by the 
requirement that force appears on any atom under rigid translation 
of the whole crystal. For the supercell used in this study, the number 
of force constants are reduced from 82944 to 52.

This allows to rewrite Equation <ref> as


    f⃗_⃗c⃗ = C(u⃗_⃗c⃗) θ⃗


where the elements of matrix C_c are linear combinations of 
displacements the coefficients of which depend on the choice of
the supercell and the aforementioned symmetry constraints.
The parameters θ⃗ can now be found with linear least
squares method using Moore-Penrose inverse of C


    θ⃗ = C^+(u⃗) f⃗


Where we have omitted the indices c since both f⃗ and
u⃗ can contain the forces and displacements of any
number of atomic configurations. The more anharmonic the system is,
the more configurations are needed in order to keep the uncertainty
of the fit sufficiently small. By expressing θ⃗ as a function of 
temperature we obtain a temperature dependent effective harmonic force
constant matrix from which the phonon frequencies and eigenvectors
can be calculated.











When working only with the vibrational modes commensurate with the 
supercell it is more convenient to consider the supercell itself to be 
the unit cell and calculate the exact normal mode frequencies which 
correspond to q⃗ = (0,0,0),  only at the Γ-point of
the reciprocal lattice of the supercell. In 
that case the dynamical matrix for the whole system is


    D = M^-1/2Φ(θ⃗) M^-1/2 = Q(θ⃗) Ω^2 (θ⃗) Q(θ⃗)^T


where M is a diagonal matrix of the masses of the atoms, Ω^2 a 
diagonal matrix of eigenvalues and Q a matrix the columns of which 
are the eigenvectors.

By substituting Equation <ref> into Equation <ref>
we get the following general equation that applies both to SCAILD and TDEP

    Q(θ⃗)^T  M^-1/2f⃗_⃗c⃗ = - Ω(θ⃗)^2 Q(θ⃗)^T M^1/2u⃗_⃗c⃗


One of the approximations of SCAILD is that the eigenvectors are 
temperature independent and can, thus, be obtained from a quick 0 K 
phonon calculation. Whether this holds, depends on the symmetry 
properties of the system and the choice of the supercell, which 
determines the grid of sampled q-points <cit.>. For the 
system considered in this work, this approximation does not have 
significant effect on the results.

Equation <ref> then becomes


    ϕ⃗_⃗c⃗ = - Ω^2 υ⃗_⃗c⃗


where ϕ⃗_⃗c⃗ = Q^T M^-1/2f⃗_⃗c⃗ is the normal 
mode force and υ⃗_⃗c⃗ = Q^T M^1/2u⃗_⃗c⃗ the 
normal mode displacement. The θ⃗ dependence is removed 
since we are interested in Ω, which can be obtained from a
simple linear least squares fit, and how its elements are related 
to the parameters becomes irrelevant.

For a perfectly harmonic crystal there exists one and only one 
Ω^2. For anharmonic crystals however the uniqueness is lost due 
to the higher order terms in the Hamiltonian becoming relevant. A 
typical distribution of the normal mode force-displacement relation for 
a single phonon mode is depicted in Figures <ref> 
and <ref>. The shape is that of a Gaussian due to both 
the nature of constant temperature sampling of the configurations and 
the anharmonicity of the system. The modes that are symmetrically 
equivalent can be placed on the same plot thereby increasing the amount 
of data, thereby reducing the uncertainty of the linear fit and 
retaining the symmetry properties of the vibrations.

As we have shown, TDEP is mathematically equivalent to SCAILD for 
systems with temperature independent phonon eigenvectors. The linear 
least squares fit is done in real space for the former and in normal 
mode space for the latter. 

Our approach differs slightly from the 
original formulation of SCAILD where sampling of the normal mode 
displacements, in order to obtain the forces, was limited to certain 
discrete values that give the same mean squared displacement as a quantum 
harmonic oscillator at a constant temperature, 


    u⃗ = M^-1/2 Q d⃗


where


    d_i = ±√(ħ/Ω_i,i[1/2 + ( expħΩ_i,i/k_B T -1 )^-1])


That makes it possible to take the average of Ω^2 over all the sampled 
configurations which means the average slope of all the lines from the 
origin to every point in the force-displacement plot. This is not 
mathematically equivalent to a linear least squares fit through all the 
points. In addition, when the displacements for SCAILD are sampled from 
a Gaussian distribution, as done in this work, the most likely 
displacement for a given mode is zero while the force due to the 
anharmonicity is nonzero which leads to many large positive and 
negative Ω^2 values which are not likely to cancel out leading 
to a large error in the estimation of the average value.



 §.§ Sampling of the configuration space


The major difference between TDEP and SCAILD is the way the 
displacements are obtained. For TDEP it is a constant temperature 
molecular dynamics simulation consistent with the Hamiltonian of the 
anharmonic system whereas for SCAILD the displacements are sampled from 
the distribution self-consistent with the current best fitted harmonic 
Hamiltonian. The major advantage of the latter is the reduction in the 
number of calculations due to more efficient generation of uncorrelated 
samples as there exists an analytical expression that gives correct 
constant temperature statistics.

Using the samples generated by a canonical molecular dynamics or Monte 
Carlo simulation can be interpreted as finding the harmonic potential 
that reproduces the forces given by the anharmonic Hamiltonian for 
constant temperature with the smallest error in least square sense. It 
is not unlike using the method of force matching quite commonly 
employed in order to fit empirical interatomic potentials to DFT 
reference data <cit.>. Lack of transferability is often 
an issue with empirical potentials and, as we explain below, it also 
plays an important role in the effective harmonic models.

It is slightly less obvious what the sampling used in SCAILD method 
yields. In addition to the resulting model being harmonic, the 
displacements themselves are generated assuming non-interacting phonon 
modes. It has been suggested that this results in probing the wrong 
parts of the phase space <cit.> which is certainly true 
when the actual anharmonic system is concerned, however the goal here 
is not to perfectly reproduce the whole anharmonic Hamiltonian, but to 
pack the anharmonic vibrational information into a truncated Taylor 
expansion of the potential. Nevertheless by careful reasoning one can 
identify the possible effects of the harmonic sampling when compared to 
the anharmonic one. Since in the former case all the correlations in 
the motion of the normal modes are ignored, the probability of sampling 
higher energy configurations increases which results in larger magnitudes
for the sampled forces and in turn higher effective vibrational
frequencies which increase with temperature.

While it is clear that the two sampling methods are different, it does 
not automatically follow that one provides a better free energy 
estimation than the other. In either case we obtain a harmonic 
Hamiltonian both of which are just approximations and so are the 
predicted free energies. While there might exist an effective harmonic 
model that gives the exact free energy for the system, it is not 
immediately clear whether either of the methods should do that.



 §.§ Thermodynamic integration and free energy perturbation


Since we expect SCAILD and TDEP to give different harmonic force 
constant matrices and therefore also different free energies, the full 
vibrational free energy needs to be calculated for a proper comparison. 
Both methods provide an excellent reference system for free energy 
perturbation (FEP) and thermodynamic integration (TI) when compared to 
the Einstein crystal <cit.> which has smaller 
overlap with the phase space of the anharmonic system, or a 0 K 
harmonic model <cit.> which has undefined free 
energy in case of phonon modes with imaginary frequencies.

Both FEP and TI allow to calculate free energy differences from 
ensemble averages of energy differences. There are several variations 
of each. In this work we used a linear path through a parameter 
λ from the effective harmonic system to the anharmonic one, in 
our case DFT system as follows:


    U(λ) = λ U_DFT + (1-λ) U_EH


To get the free energy difference from thermodynamic integration
we calculate


    Δ F = ∫_0^1 ⟨∂ U(λ)/∂λ⟩_λ dλ =
     ∫_0^1 ⟨ U_DFT - U_EH⟩_λ dλ


In practice there are possibilities either to sample λ on a 
discrete, but not necessarily equidistant grid or perform adiabatic 
switching where λ changes continuously throughout the 
simulation.

Since SCAILD and TDEP sample a lot of configurations and provide both 
the harmonic and DFT energies, it is also in principle possible to 
calculate the anharmonic free energy without any extra steps through 
free energy perturbation.


    Δ F = -1/β·ln⟨exp [ - β (U_B - U_A) ] ⟩ _A


For TDEP, U_A are the potential energies from a DFT-MD calculation 
and U_B the potential energies calculated using the fitted harmonic 
model. For SCAILD, U_A are the harmonic potential energies calculated 
using the converged fit and U_B the potential energies for the same 
structures obtained from DFT. FEP can also be done using stratified 
calculations. Similarly to thermodynamic integration several 
intermediate simulations using mixed Hamiltonians of the two endpoints 
can be done and the free energy differences between each step 
calculated using Equation <ref> and subsequently summed together 
<cit.>.




 §.§ Calculations





Cubic zirconia was chosen for the system to be studied for several 
reasons. It has been studied both using 0 K <cit.>, self 
consistent lattice dynamics <cit.> and other methods 
<cit.>. At ambient pressure it is not dynamically stable as 
indicated by large amount of imaginary phonon modes which disappear 
when a SCAILD analysis is performed. Another option for achieving 
dynamical stability for cubic ZrO2 is to put it under high pressure 
which we opted to do in order to retain the ability to compute the free 
energy from 0 K phonon calculation for comparison.

As explained above, the theory simplifies a lot in cases where the 
eigenvectors can be assumed not to depend on temperature. ZrO2 contains 
two types of atoms with relatively large difference in atomic weights 
which should theoretically allow for the possibility of temperature 
dependent eigenvectors. The size of the supercell was chosen to be 2 
× 2 × 2 of the conventional unit cell with lattice parameter 
of 9.7 Å containing a total of 96 atoms. It is not guaranteed to give 
a converged free energy with respect to the size, but is a good balance 
between capturing the overall vibrational, both harmonic and 
anharmonic, behavior while allowing for sufficiently long molecular 
dynamics simulations.

All phonon calculations were done using phonopy <cit.> without 
applying non-analytical term correction. For SCAILD and TDEP our own 
implementations were used. For SCAILD we acquired the eigenvectors from 
phonopy after which the displacements were generated for every normal mode 
according to a Gaussian distribution with the same mean square 
displacement as given by theory for a classical harmonic oscillator at 
a constant temperature,  in Equation <ref>, d⃗
is obtained from


    σ_i^2 = ⟨ d_i^2  ⟩ = k_B T/Ω_i,i^2


The sampling was from a classical distribution in order to be directly 
comparable to classical DFT molecular dynamics, although at the 
investigated temperatures we expect the effect on the results to be 
relatively small.

As shown in Figure <ref> 
the changes in free energy between consecutive iterations can be orders 
of magnitude smaller than the changes over many iterations so instead 
of basing the convergence on the former we ran the simulations for a 
fixed number of (500) steps.

All the molecular dynamics simulations for TDEP were performed using 
Nosé-Hoover thermostatting with the Nosé mass corresponding to 
approximately a time period of 80 fs and velocity Verlet integration 
with a timestep of 1 fs. The symmetry requirements for the force 
constant matrix and crystal were taken into account when performing the 
linear fit of the displacement-force data. No additional cutoff 
distance in addition to that determined by the supercell was used for 
the force constants.

For both methods the simulations were done at 600, 1100 and 1600 K. At 
higher temperatures the oxygen atoms started migrating and occasionally 
creating vacancy-interstitial pairs which complicates the analysis not 
only for SCAILD and TDEP, but also for the following thermodynamic 
integration which set the upper limit used in this work. The number of 
MD steps was 20000 at 600 and 1100 K, and 30000 at 1600 K.

Thermodynamic integration was done for λ values of 0, 0.1, 
0.25, 0.5, 0.75, 0.9 and 1. For TDEP λ=1 data is already 
acquired before the fitting and no extra calculation is needed. The 
same data was also used for SCAILD at λ=1. For SCAILD at 
λ=0 a separate fast MD run was done using the harmonic 
Hamiltonian after which a subset of uncorrelated samples was chosen for 
which the energy was calculated using DFT. Whereas Nosé-Hoover 
thermostat was used for TDEP mostly because it automatically retains 
the position of the center of mass of the system and therefore 
simplifies the fitting, Langevin dynamics was used for the 
thermodynamic integration steps in order to avoid the possibility of 
incorrect sampling of the phase space when the system is close to 
harmonic <cit.>. The drift of the center of mass was not 
problematic for TI since only the energy differences were obtained 
instead of displacements.

All DFT calculations were done using VASP 
<cit.> with 
projector-augmented wave method <cit.> with a plane wave 
cutoff energy of 400 eV and GGA-PBE <cit.> XC-functional, 
<cit.>. The number of electrons treated explicitly was 12 
and 6 for Zr and O respectively. A  2 × 2 × 2 
Monkhorst-Pack Γ-centered k-point grid <cit.> and 
Gaussian smearing with σ=0.05 were used.

The choice for the k-point grid was based on the balance between 
sufficient accuracy and computational time. Fits for SCAILD and TDEP 
based on shorter runs with 3 × 3 × 3 Γ-centered grid 
at 1600 K resulted in roughly 1 meV/atom difference. For thermodynamic 
integration a check was done by recalculating the λ=0.5 point 
for TDEP at 1600 K using a 4 × 4 × 4 gamma-centered grid. 
The difference seen was within the statistical uncertainty of the 
original calculation so no extra up-sampling as done using the UP-TILD 
<cit.> method was performed. 

As shown by Hellman <cit.>, the relative errors in forces 
decrease by orders of magnitude as the magnitude of the forces 
increases to that of the thermally excited states. Based on that the 0 
K harmonic calculations including obtaining the eigenvectors for SCAILD 
were done using a denser 7 × 7 × 7 k-point grid.




§ RESULTS










The vibrational free energies calculated using all of the methods are 
shown in Table <ref>. The full anharmonic vibrational 
free energies are in good agreement. Even at the highest considered 
temperature, 1600 K, all stratified free energy perturbation and 
thermodynamic integration results are within 1 meV/atom. Both SCAILD 
and TDEP are however off by about 20 meV/atom and, perhaps 
surprisingly, even the free energy calculated from harmonic lattice 
dynamics is a lot closer to the value obtained from thermodynamic 
integration than those of the either effective harmonic method.

The error in free energies predicted by SCAILD and TDEP is opposite in 
sign. This is also evident from the thermodynamic integration results 
shown in Figure <ref>. Since TDEP force constant matrix 
is fitted to the DFT-MD results also the average potential energy 
difference between DFT and TDEP is minimal at λ=1. It is also 
seen that setting the average potential energies to be equal as done in 
Ref  (Equation 21) would do little to reduce the 
error. In this case the whole TI curve would be shifted down by roughly 
3 meV/atom for the temperature of 1600 K which accounts for only 15 % 
of the total error. Similar conclusions can be drawn for SCAILD for 
which the agreement with DFT is best at λ=0 with the error due 
to mismatched average potential energies being slightly larger and if 
taken into account provides a better estimation for the full 
vibrational free energy than TDEP.

Free energy perturbation was not able to correct all of the discrepancy 
possibly due to the limited number of samples, however the stratified 
version gives free energies that agree very well with TI results. This 
is expected since the same input data was used for both methods.

In either case shifting the potential energy does not result in change 
in the force constants and therefore the phonon dispersions which are 
shown in Figures <ref> for TDEP and 
<ref> for SCAILD. The general trend is for TDEP to 
predict a decrease and for SCAILD an increase in frequencies over all 
of the reciprocal space. A notable exception to that is the lowest 
frequency mode at X-point which would be unstable without the external 
pressure. Both methods predict an increase in the frequency albeit 
SCAILD from 2.9 THz at 0 K to 6.1 THz, and TDEP to 4.4 THz at 1600 K.
We must note that since the calculations were done at fixed volume,
the typical lowering of the frequencies due to thermal expansion
is not taken into account.

The reason why TDEP, at least when up to second order force constants 
are considered, cannot provide a better estimate than SCAILD is that an 
effective harmonic Hamiltonian fitted at either λ=0 or at 
λ=1 exhibits similar non-transferability. Whereas the average 
energy difference between DFT and EH model can be forced to be zero at one of 
the endpoints of the TI curve, its absolute value can only 
monotonically increase as TI is performed and the sampled atomic 
configurations will start to differ from those used for the fit. This 
can be proven using Gibbs-Bogoliubov inequality <cit.>, which 
states that


    (∂^2F/∂λ^2)_NVT≤ 0


As the anharmonicity increases, the transferability and, thus, the 
estimate for the full vibrational free energy becomes worse. This can 
also be explained as follows: let us consider a full infinite expansion 
of the DFT Hamiltonian. Starting from a harmonic model we carry out the 
TDEP process, that is perform molecular dynamics runs at λ=1, 
and for each subsequent run increase the number of terms in the 
expansion. For up to quadratic term the effective harmonic model can be 
fitted exactly and Δ U = 0 for any λ. As we add more 
terms Δ U at λ = 0 starts to deviate from 0 more and 
more while being kept 0 at λ = 1. Whereas LD ignores higher 
than second order terms, both SCAILD and TDEP incorporate those in an 
effective manner into the harmonic Hamiltonian and as shown, this can 
lead to a much larger error. In addition, it is not obvious that adding 
for example the third order force constants to the effective model 
Hamiltionan necessarily reduces the error, since the second order force 
constants may remain relatively unaffected.




§ CONCLUSIONS

Given temperature independent eigenvectors, the effective harmonic 
models obtained either through SCAILD or TDEP produce two types of 
systematic errors. One is due to the different choice of the ensemble 
of displacements and forces used to fit the data. Based on our results 
neither could be considered better than the other, the reason being 
that in either case it is not the remaining explicitly anharmonic 
energy obtained through thermodynamic integration that is minimized, 
but the potential energy difference at one of the endpoints of the 
λ-integration curve. The other systematic error, common to 
both, is introduced by trying to reproduce the anharmonic interactions 
in an effective harmonic manner.

The free energies predicted by either method diverge as the temperature 
increases with the general trend being SCAILD predicting an increase in 
the phonon frequencies and higher free energy than the exact value as 
opposed to TDEP which predicts decreasing frequencies and lower free 
energy. For the system considered in this work both give worse estimate 
of the free energy than 0 K harmonic model. Therefore relying on one or 
another can introduce large errors when investigating the vibrational 
thermodynamic properties of materials, especially when the results of 
either are compared with those obtained through other methods, for 
example, when using the quasiharmonic approximation to calculate the 
free energy of one phase, while one of the effective harmonic methods 
for another, in order to compute the phase stability. In addition, due 
to the non-linearity of the systematic error with respect to 
temperature, it is possible that the error in other thermodynamic 
quantities derived from free energy, such as heat capacity, will have 
even larger errors. Free energy perturbation can be used without extra 
computational cost in order to estimate the error, but since it may not 
account for all of the difference between the approximate and true free 
energy, thermodynamic integration or any other method that takes all of 
the anharmonicity explicitly into account has to be used.




§ ACKNOWLEDGMENTS

The research leading to these results has received funding from the 
Swedish Centre for Nuclear Technology (SKC). Computational resources 
were provided by Swedish National Infrastructure for Computing (SNIC).
Comments by A. Samanta and A. Tamm from LLNL are greatly appreciated.




