Energy vs. Strain Calculations

by Pasquale Pavone for exciting carbon

Purpose: In this tutorial you will learn how to set up and execute a series of calculations for strained structures. Additionally, it will be explained how to obtain the derivatives of the energy-vs-strain curves at zero strain and how these quantities are related to elastic constants.

### 0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here is a list of the scripts which are relevant for this tutorial with a short description.

• SETUP-elastic-strain.py: Python script for generating strained structures.
• EXECUTE-elastic-strain.sh: (Bash) shell script for running a series of exciting calculations.
• CHECKFIT-energy-vs-strain.py: Python script for extracting derivatives at zero strain of energy-vs-strain curves.
• PLOT-energy.py: Python visualization tool for energy-vs-strain curves.
• PLOT-status.py: Python visualization tool for RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
• PLOT-maxforce.py: Python visualization tool for the maximum amplitude of the force on the atoms during relaxation.
• PLOT-checkderiv.py: Python visualization tool for the calculation of derivatives at zero strain using the fit of energy-vs-strain curves.
• PLOT-optimized-geometry.py: Python visualization tool for relaxed coordinates of atoms in the unit cell.

From now on the symbol will indicate the shell prompt. Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries. ### 1. Theoretical background The energy of a crystal depends on the state of distorsion in which the crystal is found. Usually, the "measure" of the state of distortion is given either in terms of the physical-strain matrix: (1) \begin{align} \pmb{\varepsilon} = \left( \begin{array}{ccc} \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz} \\ \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz} \\ \varepsilon_{zx} & \varepsilon_{zy} & \varepsilon_{zz} \end{array} \right)\,, \end{align} or in terms of the Lagrangian-strain matrix: (2) \begin{align} \pmb{\eta} = \left( \begin{array}{ccc} \eta_{xx} & \eta_{xy} & \eta_{xz} \\ \eta_{yx} & \eta_{yy} & \eta_{yz} \\ \eta_{zx} & \eta_{zy} & \eta_{zz} \end{array} \right)\,, \end{align} where\pmb{\varepsilon}$and$\pmb{\eta}are symmetric matrices and are related to each other by the relationship (3) \begin{align} \pmb{\eta} = \pmb{\varepsilon}+ \frac{1}{2}\pmb{\varepsilon}^2\:. \end{align} The actual deformation of the crystal is given by (4) \begin{align} \pmb{r}'=(\pmb{1} +\pmb{\varepsilon})\cdot\pmb{r}\:. \end{align} Using the previous definitions, the total energy of a crystal can be written as a Taylor expansion in terms of powers of the Lagrangian strain (5) \begin{align} E[\pmb{\eta}]=E_0 + \sum_{ij} \sigma^{(0)}_{ij} \,\eta_{ij} + \frac{1}{2}\sum_{ij,kl} C^{(2)}_{ij,kl}\,\eta_{ij}\, \eta_{kl} + \frac{1}{6}\sum_{ij,kl,mn} C^{(3)}_{ij,kl,mn}\,\eta_{ij}\, \eta_{kl} \,\eta_{mn} + \cdots \end{align} where\sigma^0_{ij}=0is the reference (unstrained) configuration is the equilibrium one. The other coefficients in this Taylor expansion are defined as elastic constants of different order, e.g., the elastic constants of the 2nd order (called also linear elastic constants) are (6) \begin{align} C^{(2)}_{ij,kl}=\left. \frac{\partial^2 E[\pmb{\eta}]}{\partial\eta_{ij}\,\partial\eta_{kl}} \right|_{\pmb{\eta}=0}\,. \end{align} The results above are often expressed in terms of the Voigt notation, which is a way to represent a symmetric tensor by reducing its order. Thus, following this notation the strain matrices above can be represented as vectors in a six-dimentional space: (7) \begin{align} \pmb{\eta}_{\rm Voigt}\equiv\pmb{\widetilde\eta}= (\eta_1,\eta_2,\eta_3,\eta_4,\eta_5,\eta_6) \equiv (\eta_{xx},\eta_{yy},\eta_{zz},2\eta_{yz},2\eta_{xz},2\eta_{xy})\,. \end{align} The Taylor expansion of the elastic energy can be rewritten using Voigt notation as (8) \begin{align} E[\pmb{\widetilde\eta}]=E_0 + \pmb{\sigma}^{(0)}\cdot\pmb{\widetilde\eta} + \frac{1}{2}\pmb{\widetilde\eta}\cdot\pmb{C}^{(2)}\cdot\pmb{\widetilde\eta} + \cdots \end{align} where\pmb{C}^{(2)}$is the 6$\times6 matrix (9) \begin{align} \pmb{C}^{(2)} \equiv \pmb{C}= \left( \begin{array}{cccccc} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16}\\ C_{21} & C_{22} & C_{23} & C_{24} & C_{25} & C_{26}\\ C_{31} & C_{32} & C_{33} & C_{34} & C_{35} & C_{36}\\ C_{41} & C_{42} & C_{43} & C_{44} & C_{45} & C_{46}\\ C_{51} & C_{52} & C_{53} & C_{54} & C_{55} & C_{56}\\ C_{61} & C_{62} & C_{63} & C_{64} & C_{65} & C_{66} \end{array} \right)\,. \end{align} For a given crystal structure, the number of independent elastic-constants components can be reduced using the crystal symmetry. For instance, for cubic systems the matrix of the elastic constants reduces to (10) \begin{align} \pmb{C}_{\rm cubic}= \left( \begin{array}{cccccc} C_{11} & C_{12} & C_{12} & 0 & 0 & 0\\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0\\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0\\ 0 & 0 & 0 & C_{44} & 0 & 0\\ 0 & 0 & 0 & 0 & C_{44} & 0\\ 0 & 0 & 0 & 0 & 0 & C_{44} \end{array} \right)\,. \end{align} ### 2. Set up the calculations ##### i) Preparation of the input file The first step is to create a directory for each system that you want to investigate. Here, we consider the calculation of the energy-vs-strain curves for carbon in the diamond structure. However, the procedure we show you is valid for any system. Thus, we will create a directory diamond-elastic-strain and we move inside it.  mkdir diamond-elastic-strain
$cd diamond-elastic-strain  Inside this directory, we create (or copy from a previous calculation) the file input.xml corresponding to a calculation for the equilibrium structure of diamond. This file could look like the following. <input> <title>Diamond: Equilibrium structure</title> <structure speciespath="$EXCITINGROOT/species">
<crystal scale="6.714">
<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="C.xml" rmt="1.25">
<atom coord="0.00 0.00 0.00" />
<atom coord="0.25 0.25 0.25" />
</species>
</structure>

<groundstate
ngridk="8 8 8"
swidth="0.0001"
gmaxvr="14"
xctype="GGA_PBE_SOL">
</groundstate>

<relax/>

</input>


Please, remember that the input file for an exciting calculation must always be called input.xml.

Be sure to set the correct path for the exciting root directory (indicated in this example by $EXCITINGROOT) to the one pointing to the place where the exciting directory is placed. In order to do this, use the command $ SETUP-excitingroot.sh


Be sure to have in your file the appropriate command for performing the structure optimization: Deforming your system may change the relative positions of the atoms in the unit cell.

   <relax/>

##### ii) Generation of input files for distorted structures

All strains considered in this tutorial are Lagrangian strains.

In order to generate input files for a series of distorted structure, you have to run the script SETUP-elastic-strain.py. Notice that the script SETUP-elastic-strain.py always generates a working directory containing input files for different strains. Results of the current calculations will be also stored in the working directory. The directory name can be specified by adding the name in the command line.

$SETUP-elastic-strain.py DIRECTORYNAME  If no name is given, the script use the default name workdir. Very important: The working directory is overwritten each time you execute the script SETUP-elastic-strain.py. Therefore, choose different names for different calculations. The script SETUP-elastic-strain.py produces the following output on the screen (using deformation-0 as working directory). $ SETUP-elastic-strain.py deformation-0

Enter maximum Lagrangian strain [smax] >>>> 0.10

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

------------------------------------------------------------------------
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
11 =>  ( eta,    0,    0,  eta,    0,    0)  | mixed strain
12 =>  ( eta,    0,    0,    0,  eta,    0)  | mixed strain
13 =>  ( eta,    0,    0,    0,    0,  eta)  | mixed strain
14 =>  ( eta,  eta,    0,  eta,    0,    0)  | mixed strain
------------------------------------------------------------------------

Enter deformation code >>>> 0

$  In this example, (on screen) input entries are preceded by the symbol ">>>>". Entry values must be typed on the screen when requested. The first entry (in our example 0.10) represents the absolute value of the maximum strain for which we want to perform the calculation. The second entry (11) is the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one. The third (last) entry (0) is a self-explained label indicating the type of deformation. The latter is always referred to 2-dimensional strain tensors in the Voigt notation (so that, e.g., a strain value of 0.10 corresponds, for the choice 1 of the deformation code, to a linear deformation of 10% along the x direction). After running the script, a directory called deformation-0 is created, which contains input files for different strain values. ### 3. Execute the calculations 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. If a name for the working directory has been specified, then you must give it here, too. $ EXECUTE-elastic-strain.sh deformation-0

===> Output directory is "deformation-0" <===

Running exciting for file input-01.xml ----------------------------------

...

Run completed for file input-11.xml -------------------------------------

$  After the complete run, move to the working directory deformation-0. $ cd deformation-0


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

### 4. Post-processing: Extract energy derivatives

At this point, inside the directory deformation-0, you can use the python script CHECKFIT-energy-vs-strain.py for extracting derivatives at zero strain of energy-vs-strain curves.

$CHECKFIT-energy-vs-strain.py Enter maximum strain for the fit >>>> 0.10 Enter the order of derivative >>>> 2 ########################################### Fit data----------------------------------- Deformation code ==> 0 Deformation label ==> EEE000 Maximum value of the strain ==> 0.10000000 Number of strain values used ==> 11 Fit results for the derivative of order 2 Polynomial of order 2 ==> 4467.25 [GPa] Polynomial of order 3 ==> 4467.25 [GPa] Polynomial of order 4 ==> 4053.37 [GPa] Polynomial of order 5 ==> 4053.37 [GPa] Polynomial of order 6 ==> 4060.24 [GPa] Polynomial of order 7 ==> 4060.24 [GPa] ###########################################$


In this example, input entries are preceded by the symbol ">>>>". Entry values must be typed on the screen when requested. The first entry (in our example 0.10) represents the absolute value of the maximum strain for which we want to perform the calculation. The second entry (2) is the order of the derivative that we want to obtain.

The script generates the output files check-energy-derivatives and order-of-derivative, which can be used in the post-processing analysis. Results of this script can be analyzed using the visualization tool PLOT-checkderiv.py.

### 5. Post-processing: Visualization tools

All the scripts mentioned here must be executed in the directory where the energy-vs-strain, check-energy-derivatives, and order-of-derivative files are located. The scripts produce as output a PostScript file named PLOT.ps as well as a png file (PLOT.png).

i) PLOT-energy.py

This script allows for the visualization of the energy-vs-strain curve. It is executed as follows.

$PLOT-energy.py  In the following, we display the result for the example mentioned above (carbon in the diamond structure, volume-strain deformations up to 10%). ii) PLOT-checkderiv.py This is a very important tool that allows to represent the dependence of the calculated derivatives of the energy-vs-strain curve on • the range of points included in the fitting procedure ("maximum lagrangian strain"), • the maximum degree of the polynomial used in the fitting procedure ("n"). The script PLOT-checkderiv.py requires as input the check-energy-derivatives and order-of-derivative files generated by CHECKFIT-energy-vs-strain.py and is executed as follows. $ PLOT-checkderiv.py YMIN YMAX

One or two optional entries can be specified on the calling line. Assigning numerical values to these two entries, you can set the minimum (YMIN) and the maximum (YMAX) value on the vertical axis, respectively. An example of the script output is the following.

The previous plots can be used to determine the best range of deformations and order of polynomial fit for each distortion. By analyzing the plot, we note that curves corresponding to the higher order of the polynomial used in the fit show a horizontal plateau at about 4060 GPa. This can be assumed to be the converged value for the second derivative, from the point of view of the fit (further information on this topic can be found here). For this distortion type, this value equals 9 times the bulk modulus. Thus, the extracted value of the bulk modulus is about 451 GPa.

iii) PLOT-status.py

Python visualization tool for the RMS deviations of the effective SCF potential as a function of the iteration number during the SCF loop. It is executed as follows.

$PLOT-status.py LABEL  Here, the entry on the command line, LABEL, must be specified and represents the name of the directory where an exciting calculation has been or is still running. The PLOT-status.py script is particularly useful in the latter case because one can follow "live" the convergence of the calculation. An example of the PostScript output of the script is the following. Different line segments correspond to SCF calculations for different geometries during the relaxation. iv) PLOT-maxforce.py Python visualization tool for the maximum amplitude of the force on atoms during relaxation. It is useful for deformations which allow for internal relaxation of atomic positions, e.g., for the deformation with the code 7. It is executed as follows. $ PLOT-maxforce.py LABEL


The input entry definition is the same as for the script PLOT-status.py. If the symmetry of the deformation applied to the crystal is such that no extra force is applied to the atoms (e.g., as it happens for deformation 1 and 2) the output of the script PLOT-maxforce.py will be

Either data file not (yet) ready for visualization
or maximum force target reached already at the initial configuration.

Otherwise a plot will be produced as output. An example of the script output (for the deformation type 7) is the following.

The red points show the calculated value at each optimization step, whereas the blue line indicates the target value of the maximum amplitude of the force for stopping the relaxation.

v) PLOT-optimized-geometry.py

Python visualization tool for showing the optimized geometry compared to the reference (unrelaxed) geometry for the relative atomic coordinates of two atoms in the unit cell as a function of Lagrangian strain. It is useful for deformations which allow for internal relaxation of atomic positions, e.g., for the deformation with the code 7. It is executed as follows.

$PLOT-optimized-geometry.py ATOM1 ATOM2 YMIN YMAX  Up to four optional entries can be specified on the calling line. The first two entries are numeric labels of the atoms referring to the listing of the atoms in the input file input.xml. For instance, setting ATOM1 = 2 and ATOM2 = 6 (for a cell containing at least 6 atoms!) means that one considers the second and the sixth atoms, as listed in input.xml (including atoms of all species). If ATOM1 and ATOM2 are not specified, their default values are 1 and 2, respectively. In case ATOM1 and ATOM2 are explicitly given, you can set the minimum (YMIN) and the maximum (YMAX) value on the vertical axis, respectively. An example of the PostScript output of the script for the deformation code 7 is the following. Here, (Δ1, Δ2, Δ3) and (Δ1ref, Δ2ref, Δ3ref) represent the position difference vector,$\bf r$ATOM2 -$\bf r\$ATOM1, expressed in lattice coordinates, for the optimized geometry and the unrelaxed (reference) case, respectively.

### 6. Post-processing: How to derive elastic constants

Second derivatives calculated at zero strain of energy-vs-strain curves are combinations of the elastic constants Cij where the indexes i,j=1,2,…,6 are given in the Voigt notation. In the example that we are considering here, carbon in the cubic diamond structure, only 3 different elastic constants are non vanishing

• C11
• C12
• C44

In order to extract these three elastic constants, three different deformation types must be used. For cubic systems the best choice is represented by the following deformation types

• Volume strain (in our script corresponding to the label 0)
• Uniaxial strain in the 100 direction (label 1)
• Shear strain along the 111 direction (label 7)

Which in turns correspond to the following combination of elastic constants:

• label 0: 3 C11+ 6 C12 = 9 B0
• label 1: C11
• label 7: 3 C44

where B0 is the bulk modulus.

Experimental reference values for diamond:

• C11 = 1076 GPa
• C12 = 125 GPa
• C44 = 577 GPa
• B0 = 452 GPa