This tutorial was prepared for the 2023 SSCHA School by Đorđe Dangić. You can see here the video os the hands-on session:

The material needed for this tutorial can be downloaded here.

In previous lessons we saw how to calculate vibrational properties of material using SSCHA. Now we will use this acquired knowledge to calculate lattice thermal conductivity of materials. We will need dynamical matrices (auxiliary ones, not hessians) and the third order force constants (we already calculated them when we checked the dynamical stability of the system). With these we can calculate materials’ harmonic (phonon frequencies and phonon group velocities) and anharmonic properties (phonon lifetimes and spectral functions) which is all we need to calculate lattice thermal conductivity.

Lattice thermal conductivity of silicon

As a first exercise let’s calculate lattice thermal conductivity of silicon. Silicon is very harmonic material which means it’s lattice thermal conductivity is very high. This also makes it a good test case to check the equivalence of Green-Kubo and Boltzmann transport equation approaches in the limit of vanishing anharmonicity. To speed up the calculation we will use Tersoff potential to obtain the second and third order force constants. We will do this with this simple script:

import numpy as np
from quippy.potential import Potential
from ase import Atoms
import ase.io
from ase.eos import calculate_eos
from ase.units import kJ
from ase.phonons import Phonons as AsePhonons
from ase.constraints import ExpCellFilter
from ase.optimize import BFGS, QuasiNewton
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.Structure
import sscha, sscha.Ensemble, sscha.SchaMinimizer, sscha.Relax


# This function will use ASE to give us a starting
# guess for dynamical matrices
def get_starting_dynamical_matrices(structure_filename,
        potential, supercell):
    atoms = ase.io.read(structure_filename)
    atoms.set_calculator(potential)
    ecf = ExpCellFilter(atoms)
    qn = QuasiNewton(ecf)
    qn.run(fmax=0.0005)

    structure = CC.Structure.Structure()
    structure.generate_from_ase_atoms(atoms, get_masses = True)
    dyn = CC.Phonons.compute_phonons_finite_displacements(structure,
    potential, supercell = supercell)
    dyn.Symmetrize()
    dyn.ForcePositiveDefinite()

    eos = calculate_eos(atoms)
    v0, e0, B = eos.fit()
    bulk = B / kJ * 1.0e24

    return dyn, bulk


# Our input variables
temperature = 100.0
nconf = 1000
max_pop = 1000

# Load in Tersoff potential
pot = Potential('IP Tersoff', param_filename=
'../06_the_SSCHA_with_machine_learning_potentials/ip.parms.Tersoff.xml')
supercell = tuple((4*np.ones(3, dtype=int)).tolist())
dyn, bulk = get_starting_dynamical_matrices(
'../06_the_SSCHA_with_machine_learning_potentials/POSCAR', pot, supercell)


# Generate the ensemble and the minimizer objects
ensemble = sscha.Ensemble.Ensemble(dyn, T0=temperature,
supercell = dyn.GetSupercell())
ensemble.generate(N = nconf)
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)
minimizer.min_step_dyn = 0.1
minimizer.kong_liu_ratio = 0.5
minimizer.meaningful_factor = 0.001
minimizer.max_ka = 1000

# Relax structure
relax = sscha.Relax.SSCHA(minimizer, ase_calculator = pot,
N_configs = nconf, max_pop = max_pop, save_ensemble = True)
relax.vc_relax(static_bulk_modulus = bulk,
ensemble_loc = "directory_of_the_ensemble")

# Generate ensemble for third-order FC with the relaxed dynamical matrices
new_ensemble = sscha.Ensemble.Ensemble(relax.minim.dyn, T0=temperature,
supercell = relax.minim.dyn.GetSupercell())
new_ensemble.generate(N = nconf*5)
new_ensemble.compute_ensemble(pot, compute_stress = True,
stress_numerical = False, cluster = None, verbose = True)
# We minimize the free energy with this new ensemble
new_minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(new_ensemble)
new_minimizer.minim_struct = False
new_minimizer.set_minimization_step(0.1)
new_minimizer.meaningful_factor = 0.001
new_minimizer.max_ka = 10000
new_minimizer.init()
new_minimizer.run()
new_minimizer.dyn.save_qe('final_dyn')
# Update weights with a new dynamical matrice
new_ensemble.update_weights(new_minimizer.dyn, temperature)

# Calculate Hessian and the third order tensor (return_d3 = True)
dyn_hessian, d3_tensor = new_ensemble.get_free_energy_hessian(include_v4 = False,
    get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)
dyn_hessian.save_qe('hessian_dyn')

Here we used 4x4x4 supercell. You will need to converge results with respect to the size of the supercell. A good check for the convergence could be the decay of the second and third order force constants with the distance. Now that we have second and third order force constants, we can calculate lattice thermal conductivity. For this we provide following script:

from __future__ import print_function
from __future__ import division

import numpy as np
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.ForceTensor
import cellconstructor.ThermalConductivity
import time

dyn_prefix = 'final_dyn'
nqirr = 8

SSCHA_TO_MS = cellconstructor.ThermalConductivity.SSCHA_TO_MS
RY_TO_THZ = cellconstructor.ThermalConductivity.SSCHA_TO_THZ
dyn = CC.Phonons.Phonons(dyn_prefix, nqirr)

supercell = dyn.GetSupercell()

fc3 = CC.ForceTensor.Tensor3(dyn.structure,
dyn.structure.generate_supercell(supercell), supercell)

d3 = np.load("d3.npy")
fc3.SetupFromTensor(d3)
fc3 = CC.ThermalConductivity.centering_fc3(fc3)

mesh = [10,10,10]
smear = 0.03/RY_TO_THZ

tc = CC.ThermalConductivity.ThermalConductivity(dyn, fc3,
kpoint_grid = mesh, scattering_grid = mesh, smearing_scale = None,
smearing_type = 'constant', cp_mode = 'quantum', off_diag = False)

temperatures = np.linspace(200,1200,10,dtype=float)
start_time = time.time()
tc.setup_harmonic_properties(smear)
tc.write_harmonic_properties_to_file()

tc.calculate_kappa(mode = 'SRTA', temperatures = temperatures,
write_lifetimes = True, gauss_smearing = True, offdiag_mode = 'wigner',
kappa_filename = 'Thermal_conductivity_SRTA', lf_method = 'fortran-LA')

print('Calculated SSCHA kappa in: ', time.time() - start_time)

tc.calculate_kappa(mode = 'GK', write_lineshapes=False,
ne = 1000, temperatures = temperatures,
kappa_filename = 'Thermal_conductivity_GK')

print('Calculated SSCHA kappa in: ', time.time() - start_time)
# Save ThermalConductivity object for postprocessing.
tc.save_pickle()

Important parts of the script are:

  • We define mesh on which we calculate phonon properties to be the same as the mesh we are calculating scattering processes (variable mesh). This does not have to be true. In most cases scattering_grid can be much smaller than kpoint_grid. Converge your results with respect to both grids.

  • We use smearing approach to satisfy energy conservation laws. There are two ways: constant and adaptive. In the case of smearing_type = 'constant' we have to provide smearing value in Ry as the argument to setup_harmonic_properties function. In case we choose adaptive smearing, the smearing constant will be different for different phonon modes. We still can define global variable smearing_scale with which we multiply precomputed smearing constants. smearing_scale = 1.0 works pretty well in most cases. Converge your results with respect to smearing variables.

  • off_diag variable defines whether we are doing calculation with what was termed as coherent transport. This will be important for highly anharmonic materials with large bunching of phonon modes.

  • Function calculate_kappa does most of the work. Here we will describe main options:

    • mode defines which method to use to calculate lattice thermal conductivity. Options are SRTA which is Boltzmann transport equation solution in single relaxation time approximation and GK (Dangic et al.) which is Green-Kubo method that uses phonon spectral functions instead of phonon lifetimes. These two modes should give similar results in low anharmonicity materials, but different in strongly anharmonic ones.

    • gauss_smearing defines how we treat energy conservation in the calculation of self energy. If True it will use Gaussian functions, if False it will use Lorentzian functions. In case of Gaussian smearing real part of the self energy is calculated using Kramers-Kronig transformation.

    • offdiag_mode defines how we calculate coherent transport if mode = 'SRTA'. Two options: wigner (Simoncelli et al.) and gk (Isaeva et al.). If mode is GK, coherent transport is included naturally.

    • lf_method defines how lifetimes are calculated in case mode = 'SRTA'. In short you want to keep fortan-, and then add LA or P. These should give more or less same results. Additional option is SC where we solve phonon lifetimes self-consistently, meaning we account for the phonon lineshifts.

    • ne defines the number of frequency steps if we are calculating phonon lineshapes. Also important in case of lf_method = 'SC' because we solve self-consistent equation on a grid of frequency values linearly interpolating real and imaginary part. Larger is better. Converge your results with respect to ne.

This calculation should take a few minutes. The results are save in the kappa_filename.

Question:

If we check results we see that SRTA and GK results are different. Why? How can we improve this calculation?

Lattice thermal conductivity of GeTe

As a second example we will calculate lattice thermal conductivity of GeTe. GeTe is a highly anharmonic material with a phase transition from rhombohedral to cubic phase at around 700 K. This means its lattice thermal conductivity is very low. Additionally, it should show difference between SRTA and GK methods.

For SSCHA minimization we can calculate atomic properties using Gaussian Approximation Potential developed for this material. However, in the interest of time we provided the dynamical matrices calculated at 0 K and the third order force constantsin the folder 09_Thermal_conductivity_calculations_with_the_SSCHA.

Exercise:

Calculate lattice thermal conductivity of GeTe up to 1200 K (sample temperature from 300 K every 200 K). Is there a difference between GK and SRTA methods?

Exercise:

Check if coherent transport has an influence on thermal conductivity in this material system.

Finally, in case we want to do some postprocessing we can load in the previously saved ThermalConductivity object and access all previously calculated data. For example, we can calculate phonon density of states calculated with auxiliary force constants and the one calculated with phonon lineshapes:

from __future__ import print_function
from __future__ import division

# Import the modules to read the dynamical matrix
import numpy as np
import cellconstructor as CC
import cellconstructor.Phonons
import cellconstructor.ForceTensor
import cellconstructor.ThermalConductivity
import time
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

tc = CC.ThermalConductivity.load_thermal_conductivity()

# See at which temperatures we calculated stuff
tc.what_temperatures()

key = list(tc.lineshapes.keys()) # Get Ts for lineshapes

# DOS calculated from auxiliary force constants
harm_dos = tc.get_dos()
# Temperature dependent DOS calculated from lineshapes
# first two arrays are raw data
# second two is gaussian smoothed results \
#for the distance between energy points de
anharm_dos = tc.get_dos_from_lineshapes(float(key[-1]), de = 0.1)


# Plot results
fig = plt.figure(figsize=(6.4, 4.8))
gs1 = gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs1[0, 0])
ax.plot(harm_dos[0], harm_dos[1], 'k-',
zorder=0, label = 'Harmonic')
ax.plot(anharm_dos[0], anharm_dos[1], 'r-',
zorder=0, label = 'Anharmonic raw @ ' + key[-1] + ' K')
ax.plot(anharm_dos[2], anharm_dos[3], 'b-',
zorder=0, label = 'Anharmonic smooth @ ' + key[-1] + ' K')
ax.set_xlabel('Frequency (THz)')
ax.set_ylabel('Density of states')
ax.legend(loc = 'upper right')
ax.set_ylim(bottom = 0.0)
fig.savefig('test.pdf')
plt.show()

Additionally, if we want to check a specific phonon lineshape (for example at \(\Gamma\) point) we can do it with a bit of hacking:

for iqpt in range(tc.nkpt):
    if(np.linalg.norm(tc.k_points[iqpt]) == 0.0):
        break

energies = np.arange(len(tc.lineshapes[key[-1]][0,0]),
dtype=float)*tc.delta_omega + tc.delta_omega
tc.write_lineshape('Lineshape_at_Gamma',
tc.lineshapes[key[-1]][iqpt], iqpt, energies, 'no')