By fitting the energy-volume curve, we can determine several equilibrium properties, including the equilibrium energy $E_0$, equilibrium volume $V_0$, equilibrium bulk modulus $B_0$, and its derivative $B_0'$. These properties are then utilized in the Einstein model to provide an initial estimation for the thermodynamic properties, namely the heat capacity $C_v$ and the free energy $F$.
We start by importing matplotlib, numpy and the pyiron project class.
import matplotlib.pyplot as plt
import numpy as np
from pyiron import Project
%matplotlib inline
In the next step we create a project, by specifying its relative path/name.
pr = Project('thermo')
To analyse the energy volume dependence a single super cell is sufficient, so we create an iron super cell as an example.
basis = pr.create.structure.bulk(name='Fe', cubic=True, a=2.75)
basis.plot3d()
NGLWidget()
Let's first see how we can get the energy as a function of the volume mainly based on the knowledge you acquired in the first tutorial notebook. For this, we are going to vary the volume of our initial structure by up to 5 %.
Energy volume curves are commonly calculated with ab-initio simulation codes. In this example, we use GPAW, but the functionality below would work just as well with any DFT code, just by changing the job name prefix Gpaw
to the DFT code of your choice (e.g. Vasp
, Sphinx
...). In this tutorial, we select an energy cutoff of 320 eV, but the simulation would run without changing any parameter. For the rest of the DFT parameters, you can look up functions starting with job.set_*
.
# Exponent of 1/3 to translate volume strain of 5 % to line strain
line_strain_list = np.linspace(0.95, 1.05, 7)**(1 / 3) - 1
for strain in line_strain_list:
dft = pr.create.job.Gpaw(job_name=('sample', strain))
dft.set_encut(320.0) # Optional, among other plane wave dft parameters
dft.structure = basis.apply_strain(strain, return_box=True)
dft.run()
The job sample_m0d01695243 was saved and received the ID: 19694133 The job sample_m0d01123691 was saved and received the ID: 19694134 The job sample_m0d00558671 was saved and received the ID: 19694135 The job sample_0d0 was saved and received the ID: 19694136 The job sample_0d00552497 was saved and received the ID: 19694137 The job sample_0d01098989 was saved and received the ID: 19694138 The job sample_0d01639636 was saved and received the ID: 19694139
If your environment allows you to launch calculations to a queueing system, it might be useful to set
dft.server.queue = "name_of_the_queue"
dft.server.cores = number_of_cores
dft.server.run_time = run_time_in_seconds
Queueing system setup can be found in the pyiron docs, the homepage for which is here
results = {'energy': [], 'volume': []}
for strain in line_strain_list:
dft = pr.load(('sample', strain))
results['volume'].append(dft.structure.get_volume())
results['energy'].append(dft.output.energy_pot[-1])
plt.ylabel('Potential energy eV')
plt.xlabel(r'Volume $\mathrm{\AA}^3$')
plt.grid()
plt.scatter(results['volume'], results['energy'], label='DFT results');
coeff = np.polyfit(results['volume'], results['energy'], 3)
volume_mesh = np.linspace(np.min(results['volume']), np.max(results['volume']), 1000)
energy_mesh = np.polyval(coeff, volume_mesh)
plt.plot(volume_mesh, energy_mesh, 'r--', label='3rd order polynomial fit');
plt.legend();
So the results look good enough for us to start extract physical parameters:
equi_volume = np.roots(np.polyder(coeff)).min()
print('Equilibrium volume:', np.round(equi_volume, decimals=3), 'A^3')
equi_bulk_mod = np.polyval(np.polyder(coeff, 2), equi_volume) * equi_volume
print('Equilibrium bulk modulus:', np.round(equi_bulk_mod, decimals=3), 'eV/A^3')
Equilibrium volume: 21.19 A^3 Equilibrium bulk modulus: 1.753 eV/A^3
Translation of eV/A^3 to GPa takes a factor of around 160, meaning the bulk modulus measured here is around 280 GPa, which is comparable to experimental results.
So far, we examined the procedure for conducting an energy-volume curve analysis through the manual creation of jobs. It is great that it requires only around 10 to 20 lines of code, but since it's such a common routine, we also offer the master class Murnaghan
for this purpose, which offers an array of enhanced functionalities designed to streamline the process.
In order to launch a set of DFT calculations with varying volumes and do the analysis afterwards, the Murnaghan job can be used.
dft = pr.create.job.Gpaw(job_name='gpaw')
dft.set_encut(320) # Again, optional
dft.structure = basis.copy()
murn = dft.create_job(job_type=pr.job_type.Murnaghan, job_name='murn')
murn.input
Parameter | Value | Comment | |
---|---|---|---|
0 | num_points | 11 | number of sample points |
1 | fit_type | polynomial | ['polynomial', 'birch', 'birchmurnaghan', 'mur... |
2 | fit_order | 3 | order of the fit polynom |
3 | vol_range | 0.1 | relative volume variation around volume define... |
4 | axes | [x, y, z] | Axes along which the strain will be applied |
5 | strains | None | List of strains that should be calculated. If... |
In principle you can run the simulation without touching the input, but since we chose 7 points and volume strain of 5%, let's take these parameters here as well:
murn.input['num_points'] = 7
murn.input['vol_range'] = 0.05
murn.run()
The job murn was saved and received the ID: 19694140 The job murn_0_95 was saved and received the ID: 19694141 The job murn_0_9666667 was saved and received the ID: 19694142 The job murn_0_9833333 was saved and received the ID: 19694143 The job murn_1_0 was saved and received the ID: 19694144 The job murn_1_0166667 was saved and received the ID: 19694145 The job murn_1_0333333 was saved and received the ID: 19694146 The job murn_1_05 was saved and received the ID: 19694147
murn.plot();
In the first section of this notebook, we took a look at a polynomial fit (because it's simple!). In pyiron, we offer multiple ways to fit the points: Vinet, Murnaghan, Birch-Murnaghan and polynomial. We show the comparison below:
fit_dict = {
"Vinet": murn.fit_vinet(),
"Murnaghan": murn.fit_murnaghan(),
"Birch-Murnaghan": murn.fit_birch_murnaghan(),
"Polynomial": murn.fit_polynomial(),
}
names = list(fit_dict.keys())
volumes = [d["volume_eq"] for d in fit_dict.values()]
plt.ylim([min(volumes) - np.std(volumes), max(volumes) + np.std(volumes)])
plt.ylabel(r"Equilibrium volume $\mathrm{\AA}^3$")
plt.bar(names, volumes);
We can see that the fitting method gives a slight difference in the results. It is important to keep in mind that in DFT calculations the equilibrium volume can easily vary by more than one percent, depending on the type of pseudo-potential used. In this regard, we can consider all of these results to be acceptably correct.
names = list(fit_dict.keys())
volumes = [d["bulkmodul_eq"] for d in fit_dict.values()]
plt.ylim([min(volumes) - np.std(volumes), max(volumes) + np.std(volumes)])
plt.ylabel("Equilibrium bulk modulus GPa")
plt.bar(names, volumes);
We launched the DFT calculations above using the energy cutoff of 320 eV. However, in order to obtain a physically meaningful result, the energy cutoff (among other plane-wave DFT parameters) should be varied to see the convergence of the output values we are interested in. As the energy cutoff refers to the maximum kinetic energy that is allowed for the plane wave basis functions used to describe the electronic wavefunctions, choosing a higher energy cutoff means including more plane waves in the basis set, which can improve the accuracy of the calculations but also increases the computational cost. Conversely, choosing a lower energy cutoff reduces the computational cost but can lead to less accurate results. Therefore, the energy cutoff is an important parameter in plane wave DFT calculations, and it must be carefully chosen to balance accuracy and computational efficiency for a particular system or calculation. Here, we take a look at the partial convergence of the equilibrium volume and bulk modulus. For the sake of quick DFT convergence the energy cutoff is chosen between 270 eV and 320 eV (and not higher), but in the scientific context you should see how far you have to go in order to see a true convergence, which also applies for other plane-wave DFT parameters (e.g. k-point mesh)
encut_lst = np.linspace(270, 320, 6)
volume_lst = []
bulk_modulus_lst = []
for encut in encut_lst:
dft = pr.create.job.Gpaw(job_name=('dft', encut))
dft.set_encut(encut)
dft.structure = basis.copy()
murn = dft.create_job(job_type=pr.job_type.Murnaghan, job_name=dft.job_name.replace("dft", "murn"))
murn.input['num_points'] = 7
murn.input['vol_range'] = 0.05
murn.run()
volume_lst.append(murn.fit_vinet()["volume_eq"])
bulk_modulus_lst.append(murn.fit_vinet()["bulkmodul_eq"])
The job murn_270d0 was saved and received the ID: 19694148 The job murn_270d0_0_95 was saved and received the ID: 19694149 The job murn_270d0_0_9666667 was saved and received the ID: 19694150 The job murn_270d0_0_9833333 was saved and received the ID: 19694151 The job murn_270d0_1_0 was saved and received the ID: 19694152 The job murn_270d0_1_0166667 was saved and received the ID: 19694153 The job murn_270d0_1_0333333 was saved and received the ID: 19694154 The job murn_270d0_1_05 was saved and received the ID: 19694155 The job murn_280d0 was saved and received the ID: 19694156 The job murn_280d0_0_95 was saved and received the ID: 19694157 The job murn_280d0_0_9666667 was saved and received the ID: 19694158 The job murn_280d0_0_9833333 was saved and received the ID: 19694159 The job murn_280d0_1_0 was saved and received the ID: 19694160 The job murn_280d0_1_0166667 was saved and received the ID: 19694161 The job murn_280d0_1_0333333 was saved and received the ID: 19694162 The job murn_280d0_1_05 was saved and received the ID: 19694163 The job murn_290d0 was saved and received the ID: 19694164 The job murn_290d0_0_95 was saved and received the ID: 19694165 The job murn_290d0_0_9666667 was saved and received the ID: 19694166 The job murn_290d0_0_9833333 was saved and received the ID: 19694167 The job murn_290d0_1_0 was saved and received the ID: 19694168 The job murn_290d0_1_0166667 was saved and received the ID: 19694169 The job murn_290d0_1_0333333 was saved and received the ID: 19694170 The job murn_290d0_1_05 was saved and received the ID: 19694171 The job murn_300d0 was saved and received the ID: 19694172 The job murn_300d0_0_95 was saved and received the ID: 19694173 The job murn_300d0_0_9666667 was saved and received the ID: 19694174 The job murn_300d0_0_9833333 was saved and received the ID: 19694175 The job murn_300d0_1_0 was saved and received the ID: 19694176 The job murn_300d0_1_0166667 was saved and received the ID: 19694177 The job murn_300d0_1_0333333 was saved and received the ID: 19694178 The job murn_300d0_1_05 was saved and received the ID: 19694179 The job murn_310d0 was saved and received the ID: 19694180 The job murn_310d0_0_95 was saved and received the ID: 19694181 The job murn_310d0_0_9666667 was saved and received the ID: 19694182 The job murn_310d0_0_9833333 was saved and received the ID: 19694183 The job murn_310d0_1_0 was saved and received the ID: 19694184 The job murn_310d0_1_0166667 was saved and received the ID: 19694185 The job murn_310d0_1_0333333 was saved and received the ID: 19694186 The job murn_310d0_1_05 was saved and received the ID: 19694187 The job murn_320d0 was saved and received the ID: 19694188 The job murn_320d0_0_95 was saved and received the ID: 19694189 The job murn_320d0_0_9666667 was saved and received the ID: 19694190 The job murn_320d0_0_9833333 was saved and received the ID: 19694191 The job murn_320d0_1_0 was saved and received the ID: 19694192 The job murn_320d0_1_0166667 was saved and received the ID: 19694193 The job murn_320d0_1_0333333 was saved and received the ID: 19694194 The job murn_320d0_1_05 was saved and received the ID: 19694195
plt.subplot(2, 1, 1)
plt.grid()
plt.ylabel(r"Equilibrium volume $\mathrm{\AA}^3$")
plt.plot(encut_lst, volume_lst, "-o");
plt.subplot(2, 1, 2)
plt.grid()
plt.ylabel("Equilibrium bulk modulus GPa")
plt.xlabel("Energy cut-off eV")
plt.plot(encut_lst, bulk_modulus_lst, "-o");