How to Calculate the Stress Tensor

by Rostam Golesorkhtabar & Pasquale Pavone for exciting oxygen

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. The application of STRESS-exciting to the determination of the stress tensor for a non-equilibrium configuration of boron nitride in the wurtzite phase (w-BN) 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.

• STRESS-exciting-setup.py: Python script for generating structures at different strains.
• STRESS-exciting-submit.sh: (Bash) shell script for running a series of exciting calculations.
• STRESS-exciting-analyze.py: Python script for fitting the energy-vs-strain and CVe-vs-strain curves.
• STRESS-exciting-clean.sh: (Bash) shell script for cleaning unnecessary files.
• STRESS-exciting-result.py: Python script for calculating the stress components.
• exciting2sgroup.py: Python 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.

All strains considered in this tutorial are physical strains (ε).

Extra requirement: Tool for space-group determination

The scripts in this tutorial use the sgroup tool. If you have not done before, this tool should be downloaded and installed.

---

☛ How to check if sgroup is already installed

In order to check if sgroup is already installed in the computer you are working, write on the command line:

$ sgroup

Two possible messages may appear now on the screen.

• First case: sgroup is already installed:

No input file specified.

• Second case: sgroup must be installed:

bash: sgroup: command not found

---

☛ How to install sgroup

The code sgroup is a utility which allows to determine the space group and symmetry operations of a crystal structure.

• Tool name: sgroup
• Authors: B.Z. Yanchitsky, A.N. Timoshevskii
• Download sgroup here

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

$ tar xfvp adon_v1_0.tar.gz
$ cd SpaceGroups
$ make
$ cp sgroup $EXCITINGTOOLS

Now, you have everything you need for starting and performing general lattice optimization.

---

1. Choise of a reference configuration for w-BN

In the next sections, we will calculate the stress tensor of a given reference crystal configuration. Here, we show how this configuration can be chosen. If the reference configuration is at the theoretical equilibrium, the calculated stress tensor should vanish as an obvious consequence of the definition of equilibrium.

In order to make the calculation "more interesting", we decide to consider a non-equilibrium reference configuration with, e.g., a non vanishing pressure. We notice here that the pressure is related to the trace of the stress tensor as:

(1)

\begin{align} \nonumber P = -\frac{1}{3}\sum_{\alpha = 1}^{3} \sigma^{\phantom{I}}_{\alpha \alpha}\,. \end{align}

We obtain the non-equilibrium configuration by applying to an initial configuration, chosen at the experimental equilibrium, a uniform volume strain. In order to do this, we procede as alredy seen in General lattice optimization. Following the procedure described therein, we will be able to generate new configurations and to determine the corresponding pressures.

The example which is shown here is boron nitride in the wurtzite structure (w-BN), which is the most stable phase of boron nitride. A sketch of this phase of BN is shown in the following:

In order to execute the preliminary calculations, we create a working directory BN-reference and move inside it.

$ mkdir wBN-reference
$ cd wBN-reference

Here, copy the following lines inside a file input.xml corresponding to the experimetal equilibrium structure of w-BN.

<input>
 
   <title> w-BN experimental equilibrium structure</title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="4.8188">
         <basevect>  0.5  0.866025403  0.0000 </basevect>
         <basevect> -0.5  0.866025403  0.0000 </basevect>
         <basevect>  0.0  0.000000000  1.6353 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0000"/>
         <atom coord="0.333333333  0.333333333  0.5000"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.3749"/>
         <atom coord="0.333333333  0.333333333  0.8749"/>
      </species>
 
   </structure>
 
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE"
      ngridk="6 6 4"
      rgkmax="6.5">
   </groundstate>
 
   <relax
      method="harmonic"
      taunewton="2.0"/>
 
</input>

In this file, the attribute method of the relax element has been set to "harmonic". This is a convenient choice since there is only a single internal parameter which can be optimized for this structure (i.e., the relative distance along the z-direction between the two Nitrogen atoms)-

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

---

☛ Complete procedure for generating structures and estimating pressures

First, we generate input files for a series of volumes using the script OPTIMIZE-lattice.py.

$ OPTIMIZE-lattice.py 

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

     Number and name of space group: 194 (P 63/m m c)       
     Hexagonal I structure in the Laue classification.

     Which parameter would you like to optimize?
     1 ... volume                               
     2 ... c/a ratio with constant volume       
>>>> Please choose '1' or '2': 1

>>>> Please enter the maximum physical strain value  
     The suggested value is between 0.001 and 0.050: 0.03
     The maximum physical strain is 0.03

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

$

Entry values must be typed on the screen when requested. In this case, entries are the following.

1. input.xml, the name of the input file you have created.
2. 1, in order to change only the volume.
3. 0.03, the absolute value of the maximum strain for which we want to perform the calculation. Notice that in this case the physical strain ε is defined by the relationship V=V_initial·(1+ε)³, where V_initial is the volume of the unit cell defined in input.xml.
4. 7, the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one.

Now, you must move to the directory VOL and execute the calculation there using the script OPTIMIZE-submit.sh.

$ cd VOL
$ OPTIMIZE-submit.sh 

SCF calculation of "vol_01" starts --------------------------------
   Elapsed time = 1m51s

SCF calculation of "vol_02" starts --------------------------------
   Elapsed time = 1m25s

...

SCF calculation of "vol_07" starts --------------------------------
   Elapsed time = 0m52s

$

When calculations finished, run again the OPTIMIZE-lattice.py script in the current directory choosing the Birch-Murnaghan equation-of-state option.

$ OPTIMIZE-lattice.py

>>>> Murnaghan or Birch-Murnaghan EOS: [M/B] B

 ===================================================================== 
  Fit accuracy:
     Log(Final residue in [Ha]): -4.47

 Final parameters:
     E_min = -159.3893675 [Ha]
     V_min = 161.2181 [Bohr^3]
     B_0   = 369.848 [GPa]
     B'    = 4.079

 Optimized lattice parameter saved into the file: "BM-optimized.xml".
 =====================================================================

$

---

After performing the calculations described above, type the command

$ tail -n8 BM_eos.out

The result will list the volumes and pressures corresponding to the 7 configurations that have been calculated.

 Volume      E_dft-E_eos     Pressure [GPa]
144.6305     +0.00000725     +50.128
149.1499     -0.00002318     +33.731
153.7625     +0.00002142     +19.289
158.4739     -0.00000211     +6.576
163.2709     +0.00000104     -4.56
168.1687     -0.00000902     -14.324
173.1635     +0.00000460     -22.854

In the next we will focus only on the first of these configurations, the one with the smallest volume, which will be quoted in the next as reference configuration.

 Volume      E_dft-E_eos     Pressure [GPa]
144.6305     +0.00000725     +50.128

For this configuration we obtain a positive hydrostatic pressure of about 50.1 GPa.

In the next, we will calculate directly the stress tensor of this reference configuration.

---

2. Set up and perform a calculation of the stress for the reference configuration

To start the preparation of the calculation of the stress, we move to the parent directory, create a new directory wBN-stress and move inside it.

$ cd ../../
$ mkdir wBN-stress
$ cd wBN-stress

Inside this directory, we copy the input files corresponding to the reference configuration:

$ cp ../wBN-reference/VOL/vol_01/input.xml ./
$ cp ../wBN-reference/VOL/vol_01/geometry_opt.xml ./

Notice that the atomic positions (element atom) given inside the file input.xml are not the optimized ones, which are given instead in the file geometry_opt.xml. This can be seen using the following command:

$ grep "atom coord" *.xml

As a result of this command you will get on the screen :

geometry_opt.xml:      <atom coord="0.0000000000      0.0000000000      0.0000000000"/>
geometry_opt.xml:      <atom coord="0.3333333330      0.3333333330      0.5000000000"/>
geometry_opt.xml:      <atom coord="0.0000000000      0.0000000000      0.3745472900"/>
geometry_opt.xml:      <atom coord="0.3333333330      0.3333333330      0.8745472900"/>
input.xml:         <atom coord="0.000000000  0.000000000  0.0000"/>
input.xml:         <atom coord="0.333333333  0.333333333  0.5000"/>
input.xml:         <atom coord="0.000000000  0.000000000  0.3749"/>
input.xml:         <atom coord="0.333333333  0.333333333  0.8749"/>

Here, we can notice to the differences in the last components of the 3rd and 4th atom. By cutting the atomic positions contained inside geometry_opt.xml and pasting them into input.xml, you should get the following input file.

<input>
 
   <title> w-BN reference crystal configuration</title>
 
   <structure speciespath="$EXCITINGROOT/species">
 
      <crystal scale="4.8188">
         <basevect>    0.4850000000000000    0.8400446409100000    0.0000000000000000 </basevect>
         <basevect>   -0.4850000000000000    0.8400446409100000    0.0000000000000000 </basevect>
         <basevect>    0.0000000000000000    0.0000000000000000    1.5862410000000000 </basevect>
      </crystal>
      <species speciesfile="B.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.0000"/>
         <atom coord="0.333333333  0.333333333  0.5000"/>
      </species>
      <species speciesfile="N.xml" rmt="1.25">
         <atom coord="0.000000000  0.000000000  0.3745472900"/> 
         <atom coord="0.333333333  0.333333333  0.8745472900"/>
      </species>
 
   </structure>
 
   <groundstate do="fromscratch" xctype="GGA_PBE" ngridk="6 6 4" rgkmax="6.5">
   </groundstate>
 
   <relax method="harmonic" taunewton="2.0"/>
 
</input>

Notice that we have changed in this file the title element. Furthermore, you should update the name of species directory by using, e.g., the command

$ SETUP-excitingroot.sh

i) Generation of input files for distorted structures

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

$ STRESS-exciting-setup.py 

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

     Number and name of space group: 186 (P 63 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.01 and 0.05 : 0.05
     The maximum amount of strain is 0.05

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

$

Entry values must be typed on the screen when requested. In this case, the entries are the following.

1. input.xml, the name of the input file.
2. 0.05, the absolute value of the maximum strain for which we want to perform the calculation.
3. 11, the number of deformed structures equally spaced in strain, which are generated between the maximum negative strain and the maximum positive one.

ii) Execute the calculations

To execute the series of calculations, run the script STRESS-exciting-submit.sh, which will produce an output on the screen similar to the following.

$ STRESS-exciting-submit.sh 

        +--------------------------------------+
        | SCF calculation of "Dst01_01" starts |
        +--------------------------------------+
   Elapsed time = 1m33s

        +--------------------------------------+
        | SCF calculation of "Dst01_02" starts |
        +--------------------------------------+
   Elapsed time = 0m59s

        +--------------------------------------+
        | SCF calculation of "Dst01_03" starts |
        +--------------------------------------+

...

        +--------------------------------------+
        | SCF calculation of "Dst02_11" starts |
        +--------------------------------------+
   Elapsed time = 1m22s
$

After the complete run, the components of the stress tensor can be calculated.

iii) 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 11 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 STRESS-exciting-analyze.py 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 STRESS-exciting-analyze.py script is executed as follows.

$ STRESS-exciting-analyze.py

At this point, two 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

• Results for the distortion Dst02

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_Parameters. 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, a plot appear. This 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.

Analyzing the plot for the Dst01 distortion, we notice that curves corresponding to the higher-order polynomial fit used show a clear trend to become flat at a value between -1490 and -1480 kbar. This value 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). In a similar way, the results for the Dst02 distortion can be analyzed, too.

The script STRESS-exciting-analyze.py generates the file STRESS.IN, which will be discussed and used in the next section.

iv)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 STRESS.IN, which should have the form

Dst01    eps_max    Fit_order
Dst02    eps_max    Fit_order

In each row of this file, you should insert a value for eps_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: eps_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.05 and 5 as the best distortion range and polynomial fit, respectively. Thus, the STRESS.IN file will now contain

Dst01     0.05     5
Dst02     0.05     5

Now, you run the STRESS-exciting-result.py script.

$ STRESS-exciting-result.py

Finally, numerical values of the stress components are written in the output file STRESS.OUT. In this case the content of STRESS.OUT should be similar to

 _______________________________________________

 Hydrostatic pressure:  P = 495.9 [kbar]
 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

 Stress tensor [kbar]:

   [        -449.1           0.0           0.0 ]
   [           0.0        -449.1           0.0 ]
   [           0.0           0.0        -589.5 ]
 _______________________________________________

 Hydrostatic pressure:  P = 49.59 [Gpa]
 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

 Stress tensor [GPa]:

   [        -44.91          0.00          0.00 ]
   [          0.00        -44.91          0.00 ]
   [          0.00          0.00        -58.95 ]
 _______________________________________________

From this output, we can extract a positive value of about 49.6 GPa for the hydrostatic pressure of the reference configuration, which is in good agreement with the value estimated in Section 1.

Exercise

• Check the convergence of the results for the stress tensor with respect to the number of k-points, the parameter rgkmax, and the number of strain points used for the fit.