Equation of State of Elemental Crystals

Purpose: In this tutorial we describe how to find the equilibrium volume and lattice parameters of an elemental system with known crystal structure (e.g., Be with hcp structure). In the simplest case, this is achieved by generating multiple calculations with different unit cell volumes. A energy-vs-volume function E(V) is then fitted using the Birch-Murnaghan equation of state. The minimum of this function provides the equilibrium volume and bulk modulus for the elemental system with the given crystal structure.



Preparation

We are going to use the elements tool which is a python program which performs the setup and analysis of the necessary calculations. Get the instructions to download and install from the elements tool page

Example 1: Equilibrium volume of Al in fcc structure

In the following example, we want to determine the equilibrium volume and bulk modulus of fcc-Al, for different numbers of k points. In exciting the parameter ngridk defines the number of k points that we want to vary.

To get started, we create a directory at our preferred location named Al_ngridk:

$ mkdir Al_ngridk
$ cd Al_ngridk

The elements program takes one setup file as input. A sample is located in the elements directory. We copy the precast setup, named my_calcsetup.py, to our working directory and open it.
Note: In the following I will refer to the relative path to elements directory with elementspath.

$ cp elementspath/elements/my_calcsetup.py setup_Al.py
$ cd ~/Al_ngridk

Now open setup_Al.py with your preferred editor.

SETUP

The setup consists of a Python dictionary. Dictionaries are special data structures storing key value pairs in the form {'key1':value1,'key2':value2,…}. The indentation is only for better readability and has no further purpose. In our case we want to specify the values that define our structure and in addition vary ngridk.

Defining the unit cell geometry:
We want to calculate the total energy of an fcc crystal at seven different Volume steps (all cubic). Note: There should be at least 5 steps in order to fit it with Birch-Murnaghen state equation.
The crystal structure is defined by the structure value.

'structure':'fcc'

The expansion of a starting lattice parameter azero, to a series of lattice parameters with equidistant steps ('da':0.08), is performed by the scale value.

 'scale' : { 
                             'azero': 7.65, 
                             'da': 0.08, 
                             'steps': 7},

So the result would be [7.41,7.49,7.57,7.65,7.73,7.81,7.89].

Defining the species

'species':'Al'

Defining the number of k-points

'ngridk':[2,4,6,8]

In this case all our calculations will be performed using 2,4,6 and 8 k-points. Note: For cubic systems, the number of k-points is the same for every direction in space.
This makes a total of 4*7=28 calculations.

Ignoring the other predefined values for the moment, the setup file should look like the following.

 { 
        'param': { 
                  'scale' : { 
                             'azero': 7.65, 
                             'da': 0.08, 
                             'steps': 7},
                  'rmt': [2.0],
                  'rgkmax': [6],
                  'ngridk' : [2,4,6,8],
                  'swidth': [0.01],
                  'covera' :{ 
                             'coverazero': 1.6,  
                             'dcovera': 1.6/50, 
                             'steps': 7 } 
                  },
         'species': 'Al',
         'structure':'fcc',
         'mod': 'eos',
          'exectemplate':"shelcommand.xsl",
         'calculate' : 'True'}

Ready to perform the calculations!

EXECUTE & CALCULATE

In order to perform sequential calculations on a local machine, you first have to tell the program where to find the exciting binary.
For this, open the template called shelcommand.xsl, located at elementshome/templates and change the following path to fit the location of your exciting binary.

...
<xsl:text> 
      /path/to/exciting/bin/excitingsmp
      cd -
</xsl:text>
...

To automatically generate the input files and start the calculations, you have to call the elements.py script with your setup as argument.

$ python elementspath/elements/elements.py setup_Al.py

A shell script called execute is generated. Call it to perform calculations.

$ ./execute

DATA ANALYSIS

The collect_data script will extract energy values and perform the Birch-Murnaghan fit, giving equilibrium values of volume, energy, bulk modulus and its first derivative.
From within your calculation directory, execute:

$ python elementspath/collect_data.py

DATA VISUALIZATION

$ python elementspath/plot.py

A shell prompt will appear. Typing eos, will start a interactive plotting window.

What do you want to plot?
For energy vs. volume type: eos
For convergence type: conv
>>>eos

Next you will have to specify which parameters you want to plot. Simply type all to display all your calculations.

Specify which values of ngridk to plot. e.g. 2,4
>>>all
eos_Al.png

This is our resulting equation of state plot.

The calculated values of minimal volume, energy and bulk modulus are located in eos_data.xml. Open eos_data.py with an editor.
Look for the line at the end of each graph element. The following would be the result for ngridk = 8.

<graph bulk_mod="85.5558115532" d_bulk_mod="4.88761268985" equi_a="7.53077128166" equi_volume="106.772246883" min_energy="-241.919361513" norm_res_vect="5.9241543164e-05">

Now let us look at the convergence with respect to ngridk.
Execute plot.py again and type conv.

conv_Al.png

When comparing with experimental data, we find a good agreement of the lattice parameters. But there is a rather big difference in the bulk modulus.
We could improve our results by adding more volume points to our calculation and take a smaller interval in volume (now that we know the location of the minimum).

experimental calculated
B [GPa] 76 86
a [Bohr] 7.65 7.53

Example 2: Bulk-modulus of Be in hcp structure

The procedure we worked through in the previous example, is similar for every elemental cubic system.
Things change a little if we want to investigate hexagonal systems (or in general, systems with lower symmetry).
In that case we have more than one degree of freedom, and for that to vary more then one geometry parameter.

For Be in hcp structure, we have to vary two parameters (a=b, c). These are the lattice parameter a and the ratio of the two different lattice parameters c/a.

SETUP

Defining the unit cell geometry:
The crystal structure is defined by the structure value.

'structure':'hcp'
'scale' : { 
                      'azero': 4.319, 
                      'da': 0.08, 
                      'steps': 7}

In addition to the scale value, described in exercise 1, we have to generate a series of c/a values now. The expansion of a starting c/a ratio coverazero, to a series with equidistant steps ('dcovera':0.08), is performed by the covera value.

'covera' :{ 
                      'coverazero': 1.6,  
                      'dcovera': 1.6/30, 
                      'steps': 7 }

This generates a series where c/a is varied for seven constant volume steps.

Defining the species

'species':'Be'

Defining the number of k-points

We do not perform any convergence test here, because it would be too time consuming. If we would want to calculate the equation of state with 4 different values for ngridk, like above, we would have to perform 7*7*4=196 calculations.
So we change the value for the number of k-points to:

'ngridk':[8]

EXECUTE & CALCULATE

Same procedure as exercise 1

DATA ANALYSIS

Same procedure as exercise 1

DATA VISUALIZATION

$ python elementspath/plot.py

A shell prompt will appear. Typing covera will start a interactive plotting window.

What do you want to plot?
For energy vs. volume type: eos
For energy vs. c/a type: covera
For convergence type: conv
>>>covera
Specify which values of ngridk to plot. e.g. 2,4
>>>all
eos_Be.png

From these curves, we extracted the minima and fitted the energy-volume data with the Birch-Murnaghan equation.

Execute plot.py again and type eos.

coa_Be.png

This is our resulting equation of state plot.

The calculated values of minimal volume, energy and bulk modulus are located in eos_data.xml. Open eos_data.xml with an editor.
Look for the line at the end of the file.

<graph equi_a="4.19018292334" equi_coa="1.59797950525">
<eos bulk_mod="131.020361438" d_bulk_mod="3.25990672278" equi_volume="101.81242565" min_energy="-29.2092970004" param="" />

If we compare our results with experimental data, we find a rather good agreement:

experimental calculated
B [GPa] 130 131
a [Bohr] 4.32 4.19
c [Bohr] 6.77 6.69
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License