Energy volume curve

Theory

Fitting the energy volume curve allows to calculate the equilibrium energy $E_0$, the equilirbium volume $V_0$, the equilibrium bulk modulus $B_0$ and its derivative $B^{'}_0$. These quantities can then be used as part of the Einstein model to get an initial prediction for the thermodynamik properties, the heat capacity $C_v$ and the free energy $F$.

Initialisation

We start by importing matplotlib, numpy and the pyiron project class.

In [1]:
%matplotlib inline 
import matplotlib.pyplot as plt
import numpy as np
from pyiron import Project

In the next step we create a project, by specifying the name of the project.

In [2]:
pr = Project(path='thermo')

Atomistic structure

To analyse the energy volume dependence a single super cell is sufficient, so we create an iron super cell as an example.

In [3]:
basis = pr.create_structure(element='Fe', bravais_basis='bcc', lattice_constant=2.75)
basis.plot3d()

Calculation

Energy volume curves are commonly calculated with ab initio codes, so we use VASP in this example. But we focus on the generic commands so the same example works with any DFT code. We choose 'vasp' as job name prefix, select an energy cut off of $320 eV$ and assign the basis to the job. Afterwards we apply the corresponding strain.

In [4]:
for strain in np.linspace(0.95, 1.05, 7):
    strain_str = str(strain).replace('.', '_')
    job_vasp_strain = pr.create_job(job_type=pr.job_type.Gpaw, job_name='gpaw_' + strain_str)
    job_vasp_strain.set_encut(320.0)
    job_vasp_strain.structure = basis.copy()
    job_vasp_strain.structure.set_cell(cell=basis.cell * strain ** (1/3), scale_atoms=True)
    job_vasp_strain.run()
The job gpaw_0_95 was saved and received the ID: 1
The job gpaw_0_9666666666666667 was saved and received the ID: 2
The job gpaw_0_9833333333333333 was saved and received the ID: 3
The job gpaw_1_0 was saved and received the ID: 4
The job gpaw_1_0166666666666666 was saved and received the ID: 5
The job gpaw_1_0333333333333334 was saved and received the ID: 6
The job gpaw_1_05 was saved and received the ID: 7

As these are simple calculation, there is no need to submit them to the queuing sytem. We can confirm the status of the calculation with the job_table. If the status of each job is marked as finished, then we can continue with the next step.

In [5]:
pr.job_table()
Out[5]:
id status chemicalformula job subjob projectpath project timestart timestop totalcputime computer hamilton hamversion parentid masterid
0 1 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/ 2020-10-02 17:24:24.200176 None None [email protected]#1 GpawJob None None None
1 2 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/ 2020-10-02 17:26:45.417210 None None [email protected]#1 GpawJob None None None
2 3 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/ 2020-10-02 17:28:37.112334 None None [email protected]#1 GpawJob None None None
3 4 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/ 2020-10-02 17:30:26.714705 None None [email protected]#1 GpawJob None None None
4 5 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/ 2020-10-02 17:31:58.800251 None None [email protected]#1 GpawJob None None None
5 6 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/ 2020-10-02 17:34:47.304029 None None [email protected]#1 GpawJob None None None
6 7 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/ 2020-10-02 17:36:23.322563 None None [email protected]#1 GpawJob None None None

Analysis

We aggregate the data for further processing in two separated lists, one for the volumes and one for the energies. To do so we iterate over the jobs within the project, filter the job names which contain the string 'vasp' and from those extract the final volume and the final energy.

In [6]:
volume_lst, energy_lst = zip(*[[job['output/generic/volume'][-1], job['output/generic/energy_pot'][-1]] 
                               for job in pr.iter_jobs(convert_to_object=False) if 'gpaw' in job.job_name])   

We plot the aggregated data using matplotlib.

In [7]:
plt.plot(volume_lst, energy_lst, 'x-')
plt.xlabel('Volume ($\AA ^ 3$)')
plt.ylabel('Energy (eV)')
Out[7]:
Text(0, 0.5, 'Energy (eV)')

Encut Dependence

To extend the complexity of our simulation protocol we can not only iterate over different strains but also different energy cutoffs. For this we use multiple sub projects to structure the data. And we summarize the previous code in multiple functions to maintain a high level of readability. The first function calculates a specific strained configuration for an specifc energy cut off, while the second function analyses the different strained calculations for a specific energy cutoff and returns the list of energy volume pairs.

Functions

In [8]:
def vasp_calculation_for_strain(pr, basis, strain, encut):
    strain_str = str(strain).replace('.', '_')
    job_vasp_strain = pr.create_job(job_type=pr.job_type.Gpaw, job_name='gpaw_' + strain_str)
    job_vasp_strain.set_encut(encut)
    job_vasp_strain.structure = basis.copy()
    job_vasp_strain.structure.set_cell(cell=basis.cell * strain ** (1/3), scale_atoms=True)
    job_vasp_strain.run()
In [9]:
def energy_volume_pairs(pr):
    volume_lst, energy_lst = zip(*[[job['output/generic/volume'][-1], job['output/generic/energy_pot'][-1]] 
                               for job in pr.iter_jobs(convert_to_object=False) if 'gpaw' in job.job_name])
    return volume_lst, energy_lst

Calculation

With these functions we can structure our code and implement the additional for loop to include multiple energy cutoffs.

In [10]:
for encut in np.linspace(270, 320, 6):
    encut_str = 'encut_' + str(int(encut))
    pr_encut = pr.open(encut_str)
    for strain in np.linspace(0.95, 1.05, 7):
        vasp_calculation_for_strain(pr=pr_encut, 
                                    basis=basis, 
                                    strain=strain, 
                                    encut=encut)
The job gpaw_0_95 was saved and received the ID: 8
The job gpaw_0_9666666666666667 was saved and received the ID: 9
The job gpaw_0_9833333333333333 was saved and received the ID: 10
The job gpaw_1_0 was saved and received the ID: 11
The job gpaw_1_0166666666666666 was saved and received the ID: 12
The job gpaw_1_0333333333333334 was saved and received the ID: 13
The job gpaw_1_05 was saved and received the ID: 14
The job gpaw_0_95 was saved and received the ID: 15
The job gpaw_0_9666666666666667 was saved and received the ID: 16
The job gpaw_0_9833333333333333 was saved and received the ID: 17
The job gpaw_1_0 was saved and received the ID: 18
The job gpaw_1_0166666666666666 was saved and received the ID: 19
The job gpaw_1_0333333333333334 was saved and received the ID: 20
The job gpaw_1_05 was saved and received the ID: 21
The job gpaw_0_95 was saved and received the ID: 22
The job gpaw_0_9666666666666667 was saved and received the ID: 23
The job gpaw_0_9833333333333333 was saved and received the ID: 24
The job gpaw_1_0 was saved and received the ID: 25
The job gpaw_1_0166666666666666 was saved and received the ID: 26
The job gpaw_1_0333333333333334 was saved and received the ID: 27
The job gpaw_1_05 was saved and received the ID: 28
The job gpaw_0_95 was saved and received the ID: 29
The job gpaw_0_9666666666666667 was saved and received the ID: 30
The job gpaw_0_9833333333333333 was saved and received the ID: 31
The job gpaw_1_0 was saved and received the ID: 32
The job gpaw_1_0166666666666666 was saved and received the ID: 33
The job gpaw_1_0333333333333334 was saved and received the ID: 34
The job gpaw_1_05 was saved and received the ID: 35
The job gpaw_0_95 was saved and received the ID: 36
The job gpaw_0_9666666666666667 was saved and received the ID: 37
The job gpaw_0_9833333333333333 was saved and received the ID: 38
The job gpaw_1_0 was saved and received the ID: 39
The job gpaw_1_0166666666666666 was saved and received the ID: 40
The job gpaw_1_0333333333333334 was saved and received the ID: 41
The job gpaw_1_05 was saved and received the ID: 42
The job gpaw_0_95 was saved and received the ID: 43
The job gpaw_0_9666666666666667 was saved and received the ID: 44
The job gpaw_0_9833333333333333 was saved and received the ID: 45
The job gpaw_1_0 was saved and received the ID: 46
The job gpaw_1_0166666666666666 was saved and received the ID: 47
The job gpaw_1_0333333333333334 was saved and received the ID: 48
The job gpaw_1_05 was saved and received the ID: 49

Analysis

The analysis is structured in a similar way. Here we use iter_groups() to iterate over the existing subprojects within our project and plot the individual energy volume curves using the functions defined above.

In [11]:
for pr_encut in pr.iter_groups():
    volume_lst, energy_lst = energy_volume_pairs(pr_encut)
    plt.plot(volume_lst, energy_lst, 'x-', label=pr_encut.base_name)
plt.xlabel('Volume ($\AA ^ 3$)')
plt.ylabel('Energy (eV)')
plt.legend()
Out[11]:
<matplotlib.legend.Legend at 0x7f638828d6d0>

Fitting

After we created multiple datasets we can now start to fit the converged results. While it is possible to fit the results using a simple polynomial fit we prefer to use the phyiscally motivated birch murnaghan equation or the vinet equation. For this we create the Murnaghan object and use it is fitting functionality:

In [12]:
murn = pr.create_job(job_type=pr.job_type.Murnaghan, job_name='murn')

Birch Marnaghan

In [13]:
[e0, b0, bP, v0], [e0_error, b0_error, bP_error, v0_error] = murn._fit_leastsq(volume_lst=volume_lst,
                                                                               energy_lst=energy_lst,
                                                                               fittype='birchmurnaghan')
[e0, b0, bP, v0] 
Out[13]:
[-16.623387215037408, 280.9875784436634, 4.060730693834813, 21.19199044120968]

Vinet

In [14]:
[e0, b0, bP, v0], [e0_error, b0_error, bP_error, v0_error] = murn._fit_leastsq(volume_lst=volume_lst,
                                                                               energy_lst=energy_lst,
                                                                               fittype='vinet')
[e0, b0, bP, v0] 
Out[14]:
[-16.623384845899427, 280.93805771100557, 4.105492272090299, 21.19185363600345]

We see that both equation of states give slightly different results, with overall good agreement. To validate the agreement we plot the with with the original data.

In [15]:
vol_lst = np.linspace(np.min(volume_lst), np.max(volume_lst), 1000)
plt.plot(volume_lst, energy_lst, label='dft')
plt.plot(vol_lst, murn.fit_module.vinet_energy(vol_lst, e0, b0/ 160.21766208, bP, v0), label='vinet')
plt.xlabel('Volume ($\AA ^ 3$)')
plt.ylabel('Energy (eV)')
plt.legend()
Out[15]:
<matplotlib.legend.Legend at 0x7f638078d150>

Murnaghan Module

Besides the fitting capabilities the Murnaghan module can also be used to run a set of calculations. For this we define a reference job, which can be either a Vasp calculation or any other pyiron job type and then specify the input parameters for the Murnaghan job.

In [16]:
job_vasp_strain = pr.create_job(job_type=pr.job_type.Gpaw, job_name='gpaw')
job_vasp_strain.set_encut(320)
job_vasp_strain.structure = basis.copy()
In [17]:
murn = pr.create_job(job_type=pr.job_type.Murnaghan, job_name='murn')
murn.ref_job = job_vasp_strain
murn.input
Out[17]:
Parameter Value Comment
0 num_points 11 number of sample points
1 fit_type polynomial ['polynomial', 'birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']
2 fit_order 3 order of the fit polynom
3 vol_range 0.1 relative volume variation around volume defined by ref_ham

We modify the input parameters to agree with the settings used in the examples above and execute the simulation by calling the run command on the murnaghan job object.

In [18]:
murn.input['num_points'] = 7
murn.input['vol_range'] = 0.05
In [19]:
type(murn.structure)
Out[19]:
ase.atoms.Atoms
In [20]:
pr.job_table()
Out[20]:
id status chemicalformula job subjob projectpath project timestart timestop totalcputime computer hamilton hamversion parentid masterid
0 1 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/ 2020-10-02 17:24:24.200176 None None [email protected]#1 GpawJob None None None
1 2 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/ 2020-10-02 17:26:45.417210 None None [email protected]#1 GpawJob None None None
2 3 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/ 2020-10-02 17:28:37.112334 None None [email protected]#1 GpawJob None None None
3 4 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/ 2020-10-02 17:30:26.714705 None None [email protected]#1 GpawJob None None None
4 5 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/ 2020-10-02 17:31:58.800251 None None [email protected]#1 GpawJob None None None
5 6 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/ 2020-10-02 17:34:47.304029 None None [email protected]#1 GpawJob None None None
6 7 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/ 2020-10-02 17:36:23.322563 None None [email protected]#1 GpawJob None None None
7 8 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:38:00.805999 None None [email protected]#1 GpawJob None None None
8 9 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:40:27.023982 None None [email protected]#1 GpawJob None None None
9 10 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:42:55.820191 None None [email protected]#1 GpawJob None None None
10 11 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:45:10.442772 None None [email protected]#1 GpawJob None None None
11 12 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:47:59.450726 None None [email protected]#1 GpawJob None None None
12 13 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:51:07.518608 None None [email protected]#1 GpawJob None None None
13 14 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_270/ 2020-10-02 17:54:45.224784 None None [email protected]#1 GpawJob None None None
14 15 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_280/ 2020-10-02 17:58:11.528057 None None [email protected]#1 GpawJob None None None
15 16 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:00:11.919363 None None [email protected]#1 GpawJob None None None
16 17 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:02:25.229474 None None [email protected]#1 GpawJob None None None
17 18 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:05:13.598633 None None [email protected]#1 GpawJob None None None
18 19 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:08:06.130672 None None [email protected]#1 GpawJob None None None
19 20 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:11:21.717226 None None [email protected]#1 GpawJob None None None
20 21 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_280/ 2020-10-02 18:14:19.003564 None None [email protected]#1 GpawJob None None None
21 22 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:16:48.228097 None None [email protected]#1 GpawJob None None None
22 23 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:19:42.602848 None None [email protected]#1 GpawJob None None None
23 24 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:22:33.879253 None None [email protected]#1 GpawJob None None None
24 25 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:25:12.937586 None None [email protected]#1 GpawJob None None None
25 26 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:27:18.445423 None None [email protected]#1 GpawJob None None None
26 27 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:29:51.108935 None None [email protected]#1 GpawJob None None None
27 28 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_290/ 2020-10-02 18:33:09.633165 None None [email protected]#1 GpawJob None None None
28 29 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:35:24.100573 None None [email protected]#1 GpawJob None None None
29 30 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:37:51.134146 None None [email protected]#1 GpawJob None None None
30 31 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:40:15.407176 None None [email protected]#1 GpawJob None None None
31 32 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:42:27.007123 None None [email protected]#1 GpawJob None None None
32 33 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:45:20.422390 None None [email protected]#1 GpawJob None None None
33 34 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:47:26.819490 None None [email protected]#1 GpawJob None None None
34 35 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_300/ 2020-10-02 18:49:24.101232 None None [email protected]#1 GpawJob None None None
35 36 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_310/ 2020-10-02 18:51:11.902579 None None [email protected]#1 GpawJob None None None
36 37 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_310/ 2020-10-02 18:53:53.696423 None None [email protected]#1 GpawJob None None None
37 38 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_310/ 2020-10-02 18:55:31.120613 None None [email protected]#1 GpawJob None None None
38 39 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_310/ 2020-10-02 18:57:12.122217 None None [email protected]#1 GpawJob None None None
39 40 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_310/ 2020-10-02 18:58:42.202686 None None [email protected]#1 GpawJob None None None
40 41 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_310/ 2020-10-02 19:00:25.512077 None None [email protected]#1 GpawJob None None None
41 42 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_310/ 2020-10-02 19:02:08.222775 None None [email protected]#1 GpawJob None None None
42 43 finished None gpaw_0_95 /gpaw_0_95 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:03:42.105391 None None [email protected]#1 GpawJob None None None
43 44 finished None gpaw_0_9666666666666667 /gpaw_0_9666666666666667 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:05:35.407201 None None [email protected]#1 GpawJob None None None
44 45 finished None gpaw_0_9833333333333333 /gpaw_0_9833333333333333 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:07:21.099215 None None [email protected]#1 GpawJob None None None
45 46 finished None gpaw_1_0 /gpaw_1_0 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:08:51.322535 None None [email protected]#1 GpawJob None None None
46 47 finished None gpaw_1_0166666666666666 /gpaw_1_0166666666666666 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:10:14.501550 None None [email protected]#1 GpawJob None None None
47 48 finished None gpaw_1_0333333333333334 /gpaw_1_0333333333333334 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:11:57.005028 None None [email protected]#1 GpawJob None None None
48 49 finished None gpaw_1_05 /gpaw_1_05 None /home/jovyan/thermo/encut_320/ 2020-10-02 19:13:34.919295 None None [email protected]#1 GpawJob None None None
In [21]:
murn.run()
The job murn was saved and received the ID: 50
/srv/conda/envs/notebook/lib/python3.7/site-packages/ase/cell.py:17: FutureWarning: Cell object will no longer have pbc
  warnings.warn(deprecation_msg, FutureWarning)
The job strain_0_95 was saved and received the ID: 51
The job strain_0_9666667 was saved and received the ID: 52
The job strain_0_9833333 was saved and received the ID: 53
The job strain_1_0 was saved and received the ID: 54
The job strain_1_0166667 was saved and received the ID: 55
The job strain_1_0333333 was saved and received the ID: 56
The job strain_1_05 was saved and received the ID: 57
job_id:  51 finished
job_id:  52 finished
job_id:  53 finished
job_id:  54 finished
job_id:  55 finished
job_id:  56 finished
job_id:  57 finished

Afterwards we can use the build in capabilites to plot the resulting energy volume curve and fit different equations of state to the calculated energy volume pairs.

In [22]:
murn.output_to_pandas()
Out[22]:
volume energy error id equilibrium_b_prime equilibrium_bulk_modulus equilibrium_energy equilibrium_volume
0 19.757031 -16.527632 0.0 51 4.704621 280.897787 -16.623369 21.189923
1 20.103646 -16.569446 0.0 52 4.704621 280.897787 -16.623369 21.189923
2 20.450260 -16.599599 0.0 53 4.704621 280.897787 -16.623369 21.189923
3 20.796875 -16.616336 0.0 54 4.704621 280.897787 -16.623369 21.189923
4 21.143490 -16.623997 0.0 55 4.704621 280.897787 -16.623369 21.189923
5 21.490104 -16.619111 0.0 56 4.704621 280.897787 -16.623369 21.189923
6 21.836719 -16.607260 0.0 57 4.704621 280.897787 -16.623369 21.189923
In [23]:
murn.plot()
In [24]:
murn.fit_vinet()
Out[24]:
{'fit_type': 'vinet',
 'volume_eq': 21.19185363600345,
 'energy_eq': -16.623384845899427,
 'bulkmodul_eq': 280.93805771100557,
 'b_prime_eq': 4.105492272090299,
 'least_square_error': array([4.97985872e-04, 1.01905536e+01, 1.63735940e+00, 8.58104678e-03])}

Common mistakes

Not copying the basis

It is important to copy the basis before applying the strain, as the strain has to be applied on the initial structure, not the previous structure:

In [25]:
volume_lst_with_copy = []
for strain in np.linspace(0.95, 1.05, 7):
    basis_copy = basis.copy()
    basis_copy.set_cell(cell=basis.cell * strain ** (1/3), scale_atoms=True)
    volume_lst_with_copy.append(basis_copy.get_volume())
In [26]:
basis_copy = basis.copy()
volume_lst_without_copy = []
for strain in np.linspace(0.95, 1.05, 7):
    basis_copy.set_cell(cell=basis_copy.cell * strain ** (1/3), scale_atoms=True)
    volume_lst_without_copy.append(basis_copy.get_volume())
In [27]:
volume_lst_with_copy, volume_lst_without_copy
Out[27]:
([19.757031250000004,
  20.10364583333333,
  20.450260416666666,
  20.796874999999996,
  21.143489583333338,
  21.490104166666654,
  21.83671875000001],
 [19.757031250000004,
  19.098463541666664,
  18.780155815972222,
  18.780155815972222,
  19.09315841290509,
  19.729597026668593,
  20.716076878002024])

Rescaling the cell

Another common issue is the rescaling of the supercell, there are multiple options to choose from. We used the option to scale the atoms with the supercell.

In [28]:
basis_copy = basis.copy()
strain = 0.5
basis_copy.set_cell(cell=basis_copy.cell * strain ** (1/3), scale_atoms=True)
basis_copy.plot3d()

A nother typical case is rescaling the cell to increase the distance between the atoms or add vacuum. But that is not what we want to fit an energy volume curve.

In [29]:
basis_copy = basis.copy()
strain = 0.5
basis_copy.set_cell(cell=basis_copy.cell * strain ** (1/3), scale_atoms=False)
basis_copy.plot3d()

The same can be achieved by setting the basis to relative coordinates.

In [30]:
basis_copy = basis.copy()
strain = 0.5
basis_copy.set_relative()
basis_copy.cell *= strain ** (1/3)
basis_copy.plot3d()
In [31]:
basis_copy = basis.copy()
strain = 0.5
basis_copy.cell *= strain ** (1/3)
basis_copy.plot3d()
In [ ]: