Calculating equation of state parameters
Objectives
- Understand that only select calculation parameters should differ when comparing different DFT calculations.
- Understand how to calculate equation of state (EOS) parameters like the equilibrium lattice parameters, the equilibrium unit cell volume, and the bulk modulus
- Use an AiiDA workflow for the EOS
Introduction
Structural properties are basic quantities one often wants to extract from DFT codes. The parameters of the Birch-Murnaghan equation of state (EOS) are examples of such quantities. The procedure to obtain them is to calculate the ground-state total energy (or better free energy) of a material for a series of test lattice constants (or volumes). The parameters are then obtained by fitting the EOS to the calculated energies.
An important aspect for this procedure is that the different DFT calculations have to be comparable on a very fine energy scale. Typically, DFT calculatons are not totally converged with respect to all parameters that are used to specify them. One thus expects to see a dependence of the energy not only with respect to the different test lattice constants, but also with respect to numerical parameters that are unrelated to the actual physics one wants to describe. Because of this, one should take care that the inputs to the DFT calculations only differ in the one parameter one actually wants to vary. For DFT codes in general one should especially keep the k-point sets and the basis set cutoffs fixed. For FLEUR, quantities like the MT radii are additional parameters to keep fixed.
The FLEUR input generator inpgen does not know anything about a user intention to compare different DFT calculations. The input it generates is always adapted to the structural and chemical specification one provides in the inpgen input. This implies that using it for different test lattice constants will result in varying FLEUR input parameters. The differing k-point sets, MT radii, basis set cutoffs will prohibit a comparison of the different calculations.
One therefore uses the input generator only once for a single test lattice constant and then modifies the generated FLEUR input such that it can be used for different scalings of the lattice constant and that it fulfills also all other demands. The FLEUR input is then modified with respect to the scaling for each of the calculations for the different test lattice constants and with these inputs one can then perform the DFT calculations.
In the following we exercise this procedure to obtain the lattice constant and the bulk modulus of a Si crystal in diamond structure.
Obtaining the data for an EOS fit for Si in diamond structure
The starting point for the calculaton is an inpgen input for some choice of test lattice constant. It is smart to choose this test lattice constant as near to the expected equilibrium lattice constant, as possible. We motivate our choice of by the experimental lattice constant of Si. We already prepared an inpgen input for this. Change to the Si subfolder and inspect it.
cd Si
cat inpSi.txt
The setup specifies a primitive fcc unit cell (latsys='cF') with a two-atom basis. This is a common way of specifying a crystal with diamond structure. The factor a0 is used to convert the input from Angstrom to Bohr radii as a length unit.
The next step is to generate a FLEUR input from this.
inpgen -f inpSi.txt
Our set of test lattice constants will range from to in steps of . The FLEUR input has to work with each of these. For this it has to be assured that in none of the calculations the MT spheres overlap. This can be achieved by scaling the radius created by inpgen by Open the inp.xml file in an editor and replace the tag
<mtSphere radius="2.17000000" gridPoints="717" logIncrement=".01600000"/>
by
<mtSphere radius="0.97*2.17000000" gridPoints="717" logIncrement=".01600000"/>
Next we create subfolders for each of the intended calculations and copy the FLEUR input files inp.xml, kpts.xml, and sym.xml to each of the folders.
for a in 0.97 0.98 0.99 1.00 1.01 1.02 1.03
do
mkdir latt-$a
cp *.xml latt-$a
done
In the final step we now start the FLEUR calculations with modified scaling of the lattice constant. The modification of the input and the starting of the calculations can be performed in a combined way with the FLEUR command line option -xmlXPath.
for a in 0.97 0.98 0.99 1.00 1.01 1.02 1.03
do
cd latt-$a
fleur_MPI -xmlXPath "/fleurInput/cell/bulkLattice/@scale=$a" #here we modify the scaling of the lattice constant
cd -
done
Note that the -xmlXPath option will not modify the inp.xml file, but only change the respective value when parsing the input. As the FLEUR input is also written out to the out.xml file, you can inspect the changes there.
Fortunately all calculations converged with only single FLEUR runs. We extract a table with the scaling and the total energy from the output of the calculations and store it in a file energies.txt.
rm energies.txt
for a in 0.97 0.98 0.99 1.00 1.01 1.02 1.03
do
echo $a `grep "total energy" latt-$a/out |tail -1|cut -d= -f 2` >> energies.txt
done
cat energies.txt
To plot the data we already prepared a small Python script employing the matplotlib library. You are free to use and adapt it to all your needs, even after the tutorial. Inspect and invoke it.
cat ../plotEnergies.py
python ../plotEnergies.py
The script will read in the energies.txt file and create a file energies.png with the plot. Take a look at the plot. If everything worked as expected you should observe a minimum of the total energy for a scaling factor of about . The place of the energy minimum defines the scaling for the equilibrium lattice constant.
Note that we did not yet fit the equation of state to our results. The fitting is typically performed by employing external scripts. The AiiDA framework provides such functionality and thus we will perform it in the following alternative EOS workflow using AiiDA.
Obtaining the EOS with an AiiDA workflow
All the steps above after the general preparation of the inp.xml file for all the calculations together with an EOS fitting are provided as the AiiDA-FLEUR eos workflow. By default it uses lattice constant scalings from to in steps of . Let us invoke it to see the outcome.
Note
The execution of the AiiDA-FLEUR EOS workflow will take a while. You can see what is going on if you sometimes invoke verdi process list in a terminal, while the workflow is running. This command provides you with a list of AiiDA processes that are still running and you see which process depends on what. With the option -a -p 1 all processes (also finished ones) of the last day will be listed.
```bash
aiida-fleur launch eos
The results of the fitting procedure are provided in the file eos.json. Inspect it.
cat eos.json
At the bottom of the file you directly see the values obtained for the bulk modulus and its derivative. In the upper part of the file you can find the equilibrium volume volume_gs together with the related scaling scaling_gs. Note that the scaling is for the volume here. To translate it to a lattice constant scaling you would have to take the 3rd root of this value.
We can also modify the default test lattice constant scalings. For this we invoke the workflow with the -wf template.json option and inspect the generated file.
aiida-fleur launch eos -wf template.json
cat wf_eos.json
Note that the here provided values are just examples that differ from the defaults we just obtained. Let us change the parameters to those we chose for the manual workflow. This means we change the number of points to and the step size to . Modify in an editor the wf_eos.json accordingly. We launch this parametrization of the AiiDA-FLEUR EOS workflow to see whether we can reproduce the calculations from the manual approach. If you want to retain the results from the preceeding AIIDA-FLEUR EOS workflow, rename the respective files or copy them to another location.
aiida-fleur launch eos -wf wf_eos.json
As the AiiDA-FLEUR energy output is provided in eV, there is a conversion factor of about between the energy values from the manual workflow stored in energies.txt to those reported in the eos.json file. Check this for a few values.