Convergence Test for Elastic Constants

for helium version of Exciting

Purpose: In this tutorial you will find some benchmark calculations for testing the convergence of the elastic constants of Si in diamond structure as a function of parameters such as the k-point mesh and rgkmax. In addition we will benchmark the diamond as well.

Very important: We suppose you already set up all the variables and executables as described in Section 1. of the tutorial Energy vs. strain calculations.

### 2. General considerations

We will consider here as test case Silicon in the cubic diamond structure. The calculations will use the LSDA functional for the exchange-correlation term. We will use as cubic lattice constant the value a=10.210 bohr (as can be verified using the tutorial Volume Optimization for Cubic Systems).

The cubic diamond structure has three independent elastic constants:

• C11,
• C12,
• C44

For systems with this symmetry, the bulk modulus B0 can be determined directly from the energy-vs-volume curves or as a linear combination of elastic constants:

• B0=(C11+2*C12)/3

In order to calculate all three independent elastic constants, we need to consider three deformation types. In the following, we will alwsays refer to the deformation codes defined in Section 2.ii of the tutorial Energy vs. strain calculations. In particular we will use the deformation types corresponding to the following codes:

• deformation code 0 ==> volume dilatation and compression
• deformation code 1 ==> uniaxial strain
• deformation code 7 ==> shear strain along the (111) direction

The procedure we will follow is the same is for all three deformations. As Si is a semiconductor, we do not need large k-points meshes, therefore, we will start with the investigation of the following meshes:

• 4x4x4
• 8x8x8
• 10x10x10
• … (if needed).

For the moment, we keep constant the value of the rgkmax cutoff (see Input Reference for more details) to rgkmax=7.0.

### 3. Setting the input file

The first step is to create a working directory (e.g., workdir-KP4x4x4-st1) and move inside it:

``````\$ mkdir workdir-KP4x4x4-st1
\$ cd workdir-KP4x4x4-st1```
```

Inside this directory, we create a file input.xml corresponding to a calculation for the diamond structure of Si with equlibrium lattice constant of 10.210 Bohr, with the following choice of parameters for the ground-state calculation:
• k-point mesh of 4x4x4
• rgkmax=7.0
• gmaxvr=14 (cutoff for reciprocal vector of Fourier expansion)
• stype= "Methfessel-Paxton 1" (smearing method)
• swidth="0.001" (actual value of the smearing)

Furthermore, we include structure optimization letting the internal positions of the atoms relax upon deformations.

The initial input.xml file should now look like this:

```<?xml version="1.0" encoding="UTF-8"?>
<?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?>

<input xsi:noNamespaceSchemaLocation="EXCITINGROOT/xml/excitinginput.xsd"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsltpath="EXCITINGROOT/xml/">
<title>Silicon</title>
<structure speciespath="EXCITINGROOT/species">
<crystal scale="10.210">
<basevect> 0.5    0.5     0.0 </basevect>
<basevect> 0.5    0.0     0.5 </basevect>
<basevect> 0.0    0.5     0.5 </basevect>
</crystal>
<species speciesfile="Si.xml" rmt="2.0">
<atom coord="0.00 0.00 0.00" />
<atom coord="0.25 0.25 0.25" />
</species>
</structure>
<groundstate ngridk="4 4 4"
xctype="LSDAPerdew-Wang"
stype= "Methfessel-Paxton 1"
swidth="0.001"
gmaxvr="14"
rgkmax="7.0">
</groundstate>
<structureoptimization></structureoptimization>
</input>
```

where the variable EXCITINGROOT must be substituted by its actual value (see, e.g., the tutorial Energy vs. strain calculations).

### 4. Setting up and performing calculations for the deformation code 1

In this directory workdir-KP4x4x4-st1 one executes SETUP-elastic-strain.py, choosing 0.1 for the maximum strain value and 21 for number of strains. This means that our structure will be deformed according to Lagrangian strain up to a maximum (minimum) value of 0.1 (-0.1) and within these limits 21 equally spaced deformations will be generated. The last entry to be given is the deformation code. We choose the code 1, (equally we could have chosen the equivalent codes 2 or 3 for cubic systems) that is directly related to the C11 elastic constant. The computations for this deformation code are more demanding than for the code 0 (volume optimization) and less demanding than the code 4 (for the C44 elastic constant).

The screen output corresponding to our choices will be:

``````Enter maximum Lagrangian strain [smax] >>>> 0.1

Enter the number of strain values in [-smax,smax] >>>> 21

------------------------------------------------------------------------
List of deformation codes for strains in Voigt notation
------------------------------------------------------------------------
0 =>  ( eta,  eta,  eta,    0,    0,    0)  | volume strain
1 =>  ( eta,    0,    0,    0,    0,    0)  | linear strain along x
2 =>  (   0,  eta,    0,    0,    0,    0)  | linear strain along y
3 =>  (   0,    0,  eta,    0,    0,    0)  | linear strain along z
4 =>  (   0,    0,    0,  eta,    0,    0)  | yz shear strain
5 =>  (   0,    0,    0,    0,  eta,    0)  | xz shear strain
6 =>  (   0,    0,    0,    0,    0,  eta)  | xy shear strain
7 =>  (   0,    0,    0,  eta,  eta,  eta)  | shear strain along (111)
8 =>  ( eta,  eta,    0,    0,    0,    0)  | xy in-plane strain
9 =>  ( eta, -eta,    0,    0,    0,    0)  | xy in-plane shear strain
10 =>  ( eta,  eta,  eta,  eta,  eta,  eta)  | global strain
------------------------------------------------------------------------

Enter deformation code >>>> 1```
```

After running the script, a directory called workdir is created, which contains the input files for the different strain values.

Notice that the directory workdir will be overwritten each time you execute the script SETUP-elastic-strain.py. Therefore, before setting up a new run please rename the workdir directory to a different name, e.g.:

````\$ mv workdir results-label-0`
```

To execute the series of calculation with input files created by SETUP-elastic-strain.py you have to run the script EXECUTE-elastic-strain.sh.

After the complete run, the results of the calculation for the input file workdir/input.xml-i are contained in the subdirectory workdir/rundir-i where i is running from 01 to the total number of strain values. The data for the energy-vs-strain curves are contained in the file workdir/energy-vs-strain.

In a very similar way, you can proceed for calculations for the deformation codes 0 and 7.

### 5. Evaluating the elastic constants

We move now to the directory workdir. Here, we execute the script CHECKFIT-energy-vs-strain.py and enter 0.1 for the maximum strain and 2 for the second derivative of the energy (this choice correspond to the calculation the linear elastic constants). As a output, one gets the values of the second derivative of the energy for the maximum strain of 0.1, obtained using polynomials of different order for the fit of the calculated data.

The screen output corresponding to our choices will be:

``````Enter maximum strain for the interpolation >>>> 0.1

Enter the order of derivative >>>> 2

###########################################

Fit data-----------------------------------

Deformation code             ==> 1
Deformation label            ==> E00000
Maximum value of the strain  ==> 0.10000000
Number of strain values used ==> 21

Fit results for the derivative of order   2

Polynomial of order  2 ==>     154.66 [GPa]
Polynomial of order  3 ==>     154.66 [GPa]
Polynomial of order  4 ==>     156.68 [GPa]
Polynomial of order  5 ==>     156.68 [GPa]
Polynomial of order  6 ==>     164.00 [GPa]
Polynomial of order  7 ==>     164.00 [GPa]```
```

The script produces also a PostScript output in the file gnu.ps.

In this example, we get 156.7 and 164 GPa from the 4th and 6th order of the polynomial used in the fit, respectively. We can repeat the execution of this script and change the range of the strains to be used (e.g., setting the first entry to 0.05, we take into account only values of -0.05, -0.04, …, 0.04, 0.05). The results of the fitting procedure can be visualized using gnuplot in order to estimate the value of the second derivative of the energy (here, it is directly C11 in GPa). See the tutorial Energy vs. strain calculations for obtaining more information on how to extract the converged (from the point of view of the fit) value for the derivatives of the energy-vs-strain curves.

Note: To plot the energy-vs-strain curve, one can use the script GNU-energy. In this case the output file gnu.ps will contain the energy-vs-strain plot.

### 6. Convergence test of the elastic constants vs. k-point sampling

We repeat the whole procedure discussed in the previous Sections for the higher k-points meshes 6x6x6, 8x8x8, and 10x10x10. First, we evaluate the behavior of the energies by using the script GNU-energy and visualizing gnu.ps or, equivalently, just by inspection of the energy-vs-strain files corresponding to different k-points meshes. In the following, we show a total energy vs. k-point mesh plot for the deformation code 1 (strain-1): Analyzing the plot, it is obvious that the 4x4x4 k-point mesh is not satisfactory. However, we are mainly interested in convergence of the elastic constants, therefore, we will evaluate it graphically for the case of the 10x10x10 mesh: The 2nd, 4th, and 6th order polynomial fit within the strain range (0.05-0.08) exhibit a plateau values around 160.7 GPa, which can be considered the converged (from the point of view of the fit) value of C11.

### 7. Results for C11, C12, C44, and B0

The procedure described above for obtaining C11 should be repeated in order to get C44 (deformation code 7) and B0 (deformation code 0) as a function of the k-point meshes 4x4x4, 6x6x6, 8x8x8, and 10x10x10. Keep in mind than the second derivative in GPa has to be divided by 3 for C44 and by 9 for B0 (see the tutorial Energy vs. strain calculations). The following table summarizes the results for the elastic constants, including also C12 obtained from the values of B0 and C11:

k-point mesh C11 [GPa] C12 [GPa] C44 [GPa] B0 [GPa]
4x4x4 164 66.5 72 99
6x6x6 160 65.5 76 97
8x8x8 160 65 76 96.5
10x10x10 160.5 64 77 96
16x16x16 161.5 63 76 96
Experimentx 166 64 80 98

(x) Experimental values are taken from the Landolt-Bornstein tables.

Conclusions: The 8x8x8 k-point mesh for Si in the diamond structure is OK for the elastic constants, the largest deviation from experimental values is of 5% for C44.

### 8. Convergence test of C11 with respect to rgkmax at fixed k-point sampling (8x8x8)

Similarly to the k-points convergence tests, we set up and perform calculations with increased rgkmax, i.e., rgkmax=7.5 and 8.0. The calculated total energy plot (using the GNU-energy script) shows no visible change: The value of C11 as function of increasing rgkmax at fixed 8x8x8 k-point mesh:

8x8x8 rgkmax-7.0 rgkmax-7.5 rgkmax-8.0
C11 [GPa] 160 160 160

Conclusions: The increase of rgkmax to higher values than 7 does not affect both the energy-vs-strain curves and the value of C11.

### 9. Elastic constants of diamond vs. the k-point sampling.

We use the procedure described in Sections 3. and 4. with constant rgkmax=7.0 for the calculation of the elastic constants of Carbon in the diamond structure. The results are summarized below:

k-point mesh C11 [GPa] C12 [GPa] C44 [GPa] B0 [GPa]
4x4x4 1040 132.5 551 435
8x8x8 1050 124.5 558 433
10x10x10 1053 123 559 433
Experimentx 1076 125 577 452

(x) Experimental values are taken from the Landolt-Bornstein tables.

Conclusions: The values of the elastic constants are converged already for an 8x8x8 k-point mesh. The deviation from experiment is less than 4%.