# Thermal properties in the harmonic approximation¶

This section of the tutorial demonstrates how an existing FCP can be employed in conjunction with phonopy to analyze the thermal properties of a material in the harmonic approximation.

Note that this analysis by definition invokes only the second-order force constants and primarily relies on phonopy.

## Preparations¶

First a number of parameters are defined that influence the subsequent analysis. They have been moved to the top for clarity and convenience.

```"""
Calculate thermal properties within the harmonic approximation using a
hiPhive force constant potential in combination with phonopy.

Runs in approximately 30 seconds on an Intel Core i5-4670K CPU.
"""
import numpy as np
import matplotlib.pyplot as plt

from ase import Atoms
from ase.build import bulk
from hiphive import ForceConstantPotential

from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

# parameters
dim = 5  # dimension in phonopy calculation
Nq = 51  # number of q-points along each segment of the path through the BZ
mesh = [32, 32, 32]  # q-point mesh for MSD calculation
temperatures = [300, 600, 900, 1200]  # temperatures for evaluating MSD

```

Next we define the primitive structure (`prim`) along with a corresponding `PhonopyAtoms` object (`atoms_phonopy`), and initialize a `Phonopy` object (`phonopy`) that represents the chosen supercell size.

```prim = bulk('Ni')
atoms_phonopy = PhonopyAtoms(symbols=prim.get_chemical_symbols(),
scaled_positions=prim.get_scaled_positions(),
cell=prim.cell)
phonopy = Phonopy(atoms_phonopy, supercell_matrix=dim*np.eye(3),
primitive_matrix=None)

```

Now we can load the FCP generated previously from a file using the `ForceConstantPotential.read` function and retrieve the supercell (`supercell`) from the phonopy object, for which we want to set up the second-order force constant matrix. The latter is then generated using the `ForceConstantPotential.get_force_constants` method.

```fcp = ForceConstantPotential.read('fcc-nickel.fcp')
supercell = phonopy.get_supercell()
supercell = Atoms(cell=supercell.cell, numbers=supercell.numbers, pbc=True,
scaled_positions=supercell.get_scaled_positions())
fcs = fcp.get_force_constants(supercell)
print(fcs)
# access specific parts of the force constant matrix
fcs.print_force_constant((0, 1))
fcs.print_force_constant((10, 12))

```

At the end of the block, we use the `print()` function to get a compact view of the force constant matrix, which yields an output akin to the following block:

```=================== ForceConstants ===================
Orders: [2, 3, 4]
Atoms: Atoms(symbols='Ni125', pbc=True, cell=[[0.0, 8.8, 8.8], [8.8, 0.0, 8.8], [8.8, 8.8, 0.0]])

Cluster: (0, 0)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
force constant:
[[11.55459  0.      -0.     ]
[ 0.      11.55459 -0.     ]
[-0.      -0.      11.55459]]

Cluster: (0, 1)
Atom('Ni', [0.0, 0.0, 0.0], index=0)
Atom('Ni', [0.0, 1.7600000000000002, 1.7600000000000002], index=1)
force constant:
[[ 0.04545  0.      -0.     ]
[ 0.      -1.55551 -1.6022 ]
[-0.      -1.6022  -1.55551]]

...
```

Using the `ForceConstants.print_cluster` function it is furthermore possible to access the force constants for specific cluster. Here, a cluster is referred to by the indices of its sites:

```Cluster: (10, 12)
Atom('Ni', [3.5200000000000005, 0.0, 3.5200000000000005], index=10)
Atom('Ni', [3.5200000000000005, 3.5200000000000005, 7.040000000000001], index=12)
force constant:
[[ 0.00008 -0.       0.     ]
[-0.       0.002    0.00463]
[ 0.       0.00463  0.002  ]]
```

This concludes the part that involves hiPhive.

## Harmonic analysis¶

### Mean-square displacements¶

Now we can conduct an analysis of the mean-square displacement (and potentially other thermal properties) at several temperatures.

```phonopy.set_force_constants(fcs.get_fc_array(order=2))
phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=False)
phonopy.set_thermal_displacements(temperatures=temperatures)
_, msds = phonopy.get_thermal_displacements()
msds = np.sum(msds, axis=1)  # sum up the MSD over x,y,z
for temperature, msd in zip(temperatures, msds):
print('T = {:4d} K    MSD = {:.5f} A**2'.format(temperature, msd))

```

This yields:

```T =  300 K    MSD = 0.01326 A**2
T =  600 K    MSD = 0.02550 A**2
T =  900 K    MSD = 0.03797 A**2
T = 1200 K    MSD = 0.05049 A**2
```

In the molecular dynamics (MD) tutorial section these data are compared with MD results.

### Phonon dispersion¶

It is equally straightforward to extract the phonon dispersion. After one has defined a path through the Brillouin zone (BZ) this only requires two lines of code thanks to phonopy.

```def get_band(q_start, q_stop, N):
""" Return path between q_start and q_stop """
return np.array([q_start + (q_stop-q_start)*i/(N-1) for i in range(N)])

G2X = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0]), Nq)
X2K2G = get_band(np.array([0.5, 0.5, 1.0]), np.array([0, 0, 0]), Nq)
G2L = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0.5]), Nq)
bands = [G2X, X2K2G, G2L]

# get phonon dispersion
phonopy.set_band_structure(bands)
qvecs, qnorms, freqs, _ = phonopy.get_band_structure()

```

Finally, we can plot the phonon dispersion using matplotlib.

```fig = plt.figure()

kpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
kpts_labels = ['\$\\Gamma\$', 'X', '\$\\Gamma\$', 'L']

plt.axvline(x=kpts[1], color='k', linewidth=0.9)
plt.axvline(x=kpts[2], color='k', linewidth=0.9)

for q, freq, in zip(qnorms, freqs):
plt.plot(q, freq, color='b', linewidth=2.0)

plt.xlabel('Wave vector \$\\vec{q}\$', fontsize=14.0)
plt.ylabel('Frequency \$\\omega\$ (THz)', fontsize=14.0)
plt.xticks(kpts, kpts_labels, fontsize=14.0)
plt.xlim([0.0, qnorms[-1][-1]])
plt.ylim([0.0, 12.0])

plt.tight_layout()
plt.savefig('phonon_dispersion.pdf')
```

This yields the following dispersion, which is virtually indistinguishable from the dispersion one obtains when using the EMT calculator directly.

## Source code¶

The complete source code is available in `examples/tutorial/3_compute_harmonic_thermal_properties.py`

```"""
Calculate thermal properties within the harmonic approximation using a
hiPhive force constant potential in combination with phonopy.

Runs in approximately 30 seconds on an Intel Core i5-4670K CPU.
"""
import numpy as np
import matplotlib.pyplot as plt

from ase import Atoms
from ase.build import bulk
from hiphive import ForceConstantPotential

from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms

# parameters
dim = 5  # dimension in phonopy calculation
Nq = 51  # number of q-points along each segment of the path through the BZ
mesh = [32, 32, 32]  # q-point mesh for MSD calculation
temperatures = [300, 600, 900, 1200]  # temperatures for evaluating MSD

# set up phonopy
prim = bulk('Ni')
atoms_phonopy = PhonopyAtoms(symbols=prim.get_chemical_symbols(),
scaled_positions=prim.get_scaled_positions(),
cell=prim.cell)
phonopy = Phonopy(atoms_phonopy, supercell_matrix=dim*np.eye(3),
primitive_matrix=None)

# get force constants
supercell = phonopy.get_supercell()
supercell = Atoms(cell=supercell.cell, numbers=supercell.numbers, pbc=True,
scaled_positions=supercell.get_scaled_positions())
fcs = fcp.get_force_constants(supercell)
print(fcs)
# access specific parts of the force constant matrix
fcs.print_force_constant((0, 1))
fcs.print_force_constant((10, 12))

# get mean square displacements
phonopy.set_force_constants(fcs.get_fc_array(order=2))
phonopy.set_mesh(mesh, is_eigenvectors=True, is_mesh_symmetry=False)
phonopy.set_thermal_displacements(temperatures=temperatures)
_, msds = phonopy.get_thermal_displacements()
msds = np.sum(msds, axis=1)  # sum up the MSD over x,y,z
for temperature, msd in zip(temperatures, msds):
print('T = {:4d} K    MSD = {:.5f} A**2'.format(temperature, msd))

# FCC q-point paths
def get_band(q_start, q_stop, N):
""" Return path between q_start and q_stop """
return np.array([q_start + (q_stop-q_start)*i/(N-1) for i in range(N)])

G2X = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0]), Nq)
X2K2G = get_band(np.array([0.5, 0.5, 1.0]), np.array([0, 0, 0]), Nq)
G2L = get_band(np.array([0, 0, 0]), np.array([0.5, 0.5, 0.5]), Nq)
bands = [G2X, X2K2G, G2L]

# get phonon dispersion
phonopy.set_band_structure(bands)
qvecs, qnorms, freqs, _ = phonopy.get_band_structure()

# plot dispersion
fig = plt.figure()

kpts = [0.0, qnorms[0][-1], qnorms[1][-1], qnorms[2][-1]]
kpts_labels = ['\$\\Gamma\$', 'X', '\$\\Gamma\$', 'L']

plt.axvline(x=kpts[1], color='k', linewidth=0.9)
plt.axvline(x=kpts[2], color='k', linewidth=0.9)

for q, freq, in zip(qnorms, freqs):
plt.plot(q, freq, color='b', linewidth=2.0)

plt.xlabel('Wave vector \$\\vec{q}\$', fontsize=14.0)
plt.ylabel('Frequency \$\\omega\$ (THz)', fontsize=14.0)
plt.xticks(kpts, kpts_labels, fontsize=14.0)
plt.xlim([0.0, qnorms[-1][-1]])
plt.ylim([0.0, 12.0])

plt.tight_layout()
plt.savefig('phonon_dispersion.pdf')
```