Structure relaxation¶
This tutorial demonstrates how to relax a structure and calculate an energy-volume curve with only a few commands.
First we need to import a few support functions for creating the structure (bulk
), the CPUNEP
calculator and the relax_structure
function from calorine as well as functionality from matplotlib
and pandas
that we will use for plotting.
[1]:
import numpy as np
from ase.build import bulk
from ase.units import GPa
from calorine.calculators import CPUNEP
from calorine.tools import relax_structure
from matplotlib import pyplot as plt
from pandas import DataFrame
Next we create a structure, set up a CPUNEP
calculator instance, and attach the latter to the structure.
[2]:
structure = bulk('PbTe', crystalstructure='rocksalt', a=6.6)
calc = CPUNEP('PbTe_NEP3.txt')
structure.set_calculator(calc)
Now we can relax the structure by calling the relax_structure
function. To demonstrate the effect of the relaxation, we print the pressure before and after the relaxation.
[3]:
pressure = -np.sum(structure.get_stress()[:3]) / 3 / GPa
print(f'pressure before: {pressure:.2f} GPa')
relax_structure(structure)
pressure = -np.sum(structure.get_stress()[:3]) / 3 / GPa
print(f'pressure after: {pressure:.2f} GPa')
pressure before: 1.25 GPa
pressure after: 0.00 GPa
We can readily calculate the energy-volume curve from 80% to 105% of the equilibrium volume. To this end, we use again the relax_structure
function but set the constant_volume
argument to True
. This implies that the volume remains constant during the relaxation while the positions and the (volume constrained) cell metric are allowed to change. The results are first compiled into a list, which is then converted into a pandas DataFrame
object. The latter is not necessary as we
could simply use a list but it is often convenient when dealing with more complex data sets.
[4]:
data = []
for volsc in np.arange(0.8, 1.05, 0.01):
s = structure.copy()
s.cell *= volsc ** (1 / 3) # the cubic root converts from volume strain to linear strain
s.calc = calc
relax_structure(s, constant_volume=True)
data.append(dict(volume=s.get_volume() / len(structure),
energy=s.get_potential_energy() / len(structure),
pressure=-np.sum(s.get_stress()[:3]) / 3 / GPa))
df = DataFrame.from_dict(data)
Finally, we plot the energy-volume curve along with the pressure-volume curve.
[5]:
fig, ax = plt.subplots(figsize=(4.5, 3.5), dpi=90)
ax.plot(df.volume, df.energy, 'o-', label='energy')
ax2 = ax.twinx()
ax2.plot(df.volume, df.pressure, 'x-', color='C1', label='pressure')
ax.set_xlabel('Volume (Å$^3$/atom)')
ax.set_ylabel('Energy (eV/atom)')
ax2.set_ylabel('Pressure (GPa)')
fig.legend(loc='upper right', bbox_to_anchor=(0.85, 0.85));
