` 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
```

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*.

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
```

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*.

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 |