How to Calculate the Stress Tensor

by Rostam Golesorkhtabar for exciting boron

Purpose: In this tutorial you will learn how to set up and execute exciting calculations using the STRESS-exciting interface, which allows to obtain full stress tensors of crystal systems for any crystal structure. In addition, the application of STRESS-exciting to the determination of stress tensor for the hexagonal Be is explicitly presented.

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.

  • Python script for generating structures at different strains.
  • (Bash) shell script for running a series of exciting calculations.
  • Python script for fitting the energy-vs-strain and CVe-vs-strain curves.
  • (Bash) shell script for cleaning unnecessary files.
  • Python script for calculating the stress components.
  • exciting2sgroup.xsl: xsl script for converting the exciting input to an sgroup input file.
  • Grace.par: File containing xmgrace parameters needed by visualization tools.

Requirements: Bash shell. Python numpy, lxml, math, sys and os libraries. xmgrace (visualization tool).

From now on the symbol $ will indicate the shell prompt.

Extra requirement: Tool for space-group determination

STRESS-exciting uses the sgroup tool. If you have not done before, this tool should be downloaded and installed. The code sgroup is a utility which allows to determine the space group and symmetry operations of a crystal structure.

After the download, you will get a tar.gz file, go to the directory where you saved this file and execute the following commands.

$ tar xfvp DownloadedFile.tar.gz
$ cd SpaceGroups
$ make

Now, you have everything you need for starting and performing the stress calculations.

1. 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. In this tutorial, we consider as an example Be in hexagonal structure (the procedure we use is, nevertheless, valid for any system). Thus, we will create a directory Be_stress and we move inside it.

$ mkdir Be_stress
$ cd Be_stress

Inside this directory, we create or copy an structure file which can be named Be_stress.xml. In order to create an input file corresponding to a non equilibrium structure, we consider the input file used as a starting point in the General Lattice Optimization tutorial and modify the structure described in this file using a deformation which is an hydrostatic deformation with strain value of -0.03 %. Notice once more that for this structure the stress tensor is not vanishing. The input file file could look like the following.

   <title>Be: Stress Calculation</title>
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="4.300">
         <basevect>  0.970000000  0.000000000  0.000000000 </basevect>
         <basevect> -0.485000000  0.840044638  0.000000000 </basevect>
         <basevect>  0.000000000  0.000000000  1.455000000 </basevect>
      <species speciesfile="Be.xml" rmt="1.95">
         <atom coord="0.66666667 0.33333333 0.75000000"/>
         <atom coord="0.33333333 0.66666667 0.25000000"/>
   <groundstate ngridk="6 6 4">

This file can be saved with any name. In this tutorial is not necessary to rename the exciting input file as input.xml, because this file is the input of STRESS-exciting and not of exciting itself. Please, notice that the input file for a direct exciting calculation must be always 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

$ Be_stress.xml
ii) Hydrostatic pressure

The hydrostatic pressure of this structure can be calculated as in the example shown in General Lattice Optimization using the Birch-Murnaghan equation of state after the first optimization step. If you have already performed the calculation presented in General Lattice Optimization, you should have a file which called BM_eos.out in the 1-VOL directory. It contains all information which is related to the Birch-Murnaghan fit including the hydrostatic pressure. The content of this file should be similar to the following:

 === Birch-Murnaghan eos =========================
 Fit accuracy: 
     Log(Final residue in [Ha]): -6.04 

 Final parameters:
     E_min = -29.36096 [Ha]
     V_min = 106.5871 [Bohr^3]
     B_0   = 121.951 [GPa]
     B'    = 3.553

 Volume     E_dft-E_eos     Pressure [GPa]
94.2633     +0.00000009     +18.615 
98.7043     -0.00000041     +10.735
103.2857    +0.00000065     +4.057
108.0004    -0.00000046     -1.569
112.8597    +0.00000012     -6.298

As can be seen, the hydrostatic pressure of the initial structure (the first of the 5 structures which are calculated, corresponding to a volume of 94.2633 Bohr3) is +18.615 GPa. This pressure is related to the components of the stress tensor, $\sigma_{\alpha \beta}$, as following:

\begin{align} \nonumber P = -\frac{1}{3}\sum_{\alpha = 1}^{3} \sigma_{\alpha \alpha}\,. \end{align}

As a useful check, at the end of this tutorial, when all stress components are calculated, the hydrostatic pressure derived from the stress tensor should also have a value around 18.615 GPa or 186.15 kbar.

iii) Generation of input files for distorted structures

All strains considered in this tutorial are physical strains.

In order to generate input files for a series of distorted structures, you have to run the script, which will produce the following output on the screen.


>>>> Please enter the exciting input file name: Be_stress.xml

     Number and name of space group: 194 (P 63/m m c)
     Hexagonal I structure in the Laue classification.       
     This structure has 2 independent stress components.

>>>> Please enter the maximum amount of strain 
     The suggested value is between 0.0010 and 0.0100 : 0.002
     The maximum amount of strain is 0.002

>>>> Please enter the number of the distorted structures [odd number > 4]: 21
     The number of the distorted structures is 21


Entry values must be typed on the screen when requested. In this case, the entries are the following.
  1. Be_stress.xml, the name of the input file.
  2. 0.002, the absolute value of the maximum strain for which we want to perform the calculation.
  3. 21, the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one.

The script generates some new directories and files, which you can list using the following command (in the Be_stress directory).

$ ls
.  ..  Be_stress.xml  Distorted_Parameters  Dst01  Dst02  INFO_Stress  Structures_exciting  sgroup.out

You may also list the single directories, e.g., if you list Dst01 you will see the following (sub)directories.

$ ls Dst01
.   Dst01_01  Dst01_03  Dst01_05  Dst01_07  Dst01_09  Dst01_11  Dst01_13  Dst01_15  Dst01_17  Dst01_19  Dst01_21
..  Dst01_02  Dst01_04  Dst01_06  Dst01_08  Dst01_10  Dst01_12  Dst01_14  Dst01_16  Dst01_18  Dst01_20 

3. Execute the calculations

To execute the series of calculations, run (in the Be_stress directory) the script, which will produce an output on the screen similar to the following.


        | SCF calculation of "Dst01_01" starts |
   Elapsed time = 0m7s
Thu Jul 24 13:52:44 CEST 2014

        | SCF calculation of "Dst01_02" starts |
   Elapsed time = 0m7s
Thu Jul 24 13:52:51 CEST 2014

        | SCF calculation of "Dst01_03" starts |


        | SCF calculation of "Dst02_21" starts |
   Elapsed time = 0m7s
Thu Jul 24 13:57:37 CEST 2014

After the complete run, stress components can be calculated.

4. Analyzing the calculations

Stress components are obtained from energy calculations using a similar method to the one illustrated in Energy vs. strain calculations. Within this approach, first derivatives of the energy curves are evaluated with the help of a curve fitting. Calculations performed above were producing 21 points for each energy curve. The quality of the fitting procedure can be improved by increasing the number of data points per energy curve.

Now, we analyze our calculations. The script allows the analysis of the dependence of the calculated derivatives of the energy-vs-strain curve on

  1. the range of distortions included in the fitting procedure (the x axis in xmgrace plots),
  2. the degree of the polynomial fit used in the fitting procedure (different color curves in xmgrace plots).

The script is executed as follows.


At this point, four xmgrace plots will appear on your screen automatically (for more information on how to deal with xmgrace plots, see Xmgrace: A Quickstart).

  • Results for the distortion Dst01
Dst01_d1E.png Dst01_CVe.png
  • Results for the distortion Dst02
Dst02_d1E.png Dst02_CVe.png

The previous plots can be used to determine the best range of deformations and order of polynomial fit for each distortion.

As an example, we analyze the first plot, corresponding to the distortion Dst01. Distortion types are listed in the file Distorted_Parameter. By examining this file, we can see that the Dst01 distortion corresponds to an applied physical strain in the Voigt notation with the form (η,η,η,0,0,0), where η is a strain parameter. This deformation type is directly connected with the hydrostatic pressure. For each distorsion type, two plots appear. The first plot contains the first derivative of the energy with respect to the strain parameter η as a function of the maximum strain and of the order of polynomial fit. In a similar way, the second plot shows the corresponding cross-validation error. By analyzing the first plot for the Dst01 distortion, we notice that curves corresponding to the lower order of the polynomial used in the fit show a horizontal plateau at about -555 kbar. This can be assumed to be the converged value for the first 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 the negative of 3 times the pressure. Thus, the extracted value of the pressure is 185 kbar.

The script generates the file, which will be discussed and used in the next section.

5. Numerical values of the stress components

In order to obtain the numerical values of the stress tensor, you should use the following procedure. The first step is to edit the file, which should have the form

Dst01    eta_max    Fit_order
Dst02    eta_max    Fit_order

In each row of this file, you should insert a value for eta_max and Fit_order. According to the result of the visual analysis of the previous figures, these two values are extracted by considering the plateau regions of the corresponding plot: eta_max is a value of maximum distortion located in the plateau region of the curve representing the polynomial fit of order Fit_order. In our case, we can choose for each distortion 0.002 and 1 as the best distortion range and polynomial fit, respectively. Thus, the file should contain

Dst01     0.002     1
Dst02     0.002     1

Now, you run the script.


Finally, numerical values of the stress components are written in the output file STRESS.OUT.

In order to save storage space on disk, you may use the command


to delete unnecessary output files.

  • Repeat all calculations using a denser k-point grid and a larger number of strain points.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License