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 kpoint mesh and rgkmax. In addition we will benchmark the diamond as well.
Table of Contents

1. Define relevant shell variables and download scripts
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 exchangecorrelation 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:
 C_{11},
 C_{12},
 C_{44}
For systems with this symmetry, the bulk modulus B_{0} can be determined directly from the energyvsvolume curves or as a linear combination of elastic constants:
 B_{0}=(C_{11}+2*C_{12})/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 kpoints 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., workdirKP4x4x4st1) and move inside it:
$ mkdir workdirKP4x4x4st1
$ cd workdirKP4x4x4st1
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 groundstate calculation:
 kpoint mesh of 4x4x4
 rgkmax=7.0
 gmaxvr=14 (cutoff for reciprocal vector of Fourier expansion)
 stype= "MethfesselPaxton 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="UTF8"?> <?xmlstylesheet href="inputtohtml.xsl" type="text/xsl"?> <input xsi:noNamespaceSchemaLocation="EXCITINGROOT/xml/excitinginput.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchemainstance" 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="LSDAPerdewWang" stype= "MethfesselPaxton 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 workdirKP4x4x4st1 one executes SETUPelasticstrain.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 C_{11} 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 C_{44} 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 inplane strain
9 => ( eta, eta, 0, 0, 0, 0)  xy inplane 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 SETUPelasticstrain.py. Therefore, before setting up a new run please rename the workdir directory to a different name, e.g.:
$ mv workdir resultslabel0
To execute the series of calculation with input files created by SETUPelasticstrain.py you have to run the script EXECUTEelasticstrain.sh.
After the complete run, the results of the calculation for the input file workdir/input.xmli are contained in the subdirectory workdir/rundiri where i is running from 01 to the total number of strain values. The data for the energyvsstrain curves are contained in the file workdir/energyvsstrain.
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 CHECKFITenergyvsstrain.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 C_{11} 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 energyvsstrain curves.
Note: To plot the energyvsstrain curve, one can use the script GNUenergy. In this case the output file gnu.ps will contain the energyvsstrain plot.
6. Convergence test of the elastic constants vs. kpoint sampling
We repeat the whole procedure discussed in the previous Sections for the higher kpoints meshes 6x6x6, 8x8x8, and 10x10x10. First, we evaluate the behavior of the energies by using the script GNUenergy and visualizing gnu.ps or, equivalently, just by inspection of the energyvsstrain files corresponding to different kpoints meshes. In the following, we show a total energy vs. kpoint mesh plot for the deformation code 1 (strain1):
Analyzing the plot, it is obvious that the 4x4x4 kpoint 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.050.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 C_{11}.
7. Results for C_{11}, C_{12}, C_{44}, and B_{0}
The procedure described above for obtaining C_{11} should be repeated in order to get C_{44} (deformation code 7) and B_{0} (deformation code 0) as a function of the kpoint meshes 4x4x4, 6x6x6, 8x8x8, and 10x10x10. Keep in mind than the second derivative in GPa has to be divided by 3 for C_{44} and by 9 for B_{0} (see the tutorial Energy vs. strain calculations). The following table summarizes the results for the elastic constants, including also C_{12} obtained from the values of B_{0} and C_{11}:
kpoint mesh  C_{11} [GPa]  C_{12} [GPa]  C_{44} [GPa]  B_{0} [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 
Experiment^{x}  166  64  80  98 
(x) Experimental values are taken from the LandoltBornstein tables.
Conclusions: The 8x8x8 kpoint mesh for Si in the diamond structure is OK for the elastic constants, the largest deviation from experimental values is of 5% for C_{44}.
8. Convergence test of C_{11} with respect to rgkmax at fixed kpoint sampling (8x8x8)
Similarly to the kpoints 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 GNUenergy script) shows no visible change:
The value of C_{11} as function of increasing rgkmax at fixed 8x8x8 kpoint mesh:
8x8x8  rgkmax7.0  rgkmax7.5  rgkmax8.0 

C_{11} [GPa]  160  160  160 
Conclusions: The increase of rgkmax to higher values than 7 does not affect both the energyvsstrain curves and the value of C_{11}.
9. Elastic constants of diamond vs. the kpoint 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:
kpoint mesh  C_{11} [GPa]  C_{12} [GPa]  C_{44} [GPa]  B_{0} [GPa] 

4x4x4  1040  132.5  551  435 
8x8x8  1050  124.5  558  433 
10x10x10  1053  123  559  433 
Experiment^{x}  1076  125  577  452 
(x) Experimental values are taken from the LandoltBornstein tables.
Conclusions: The values of the elastic constants are converged already for an 8x8x8 kpoint mesh. The deviation from experiment is less than 4%.