by Rostam Golesorkhtabar, Pasquale Pavone, & Konstantin Lion for exciting nitrogen
Purpose: In this tutorial you will learn how to set up and execute exciting calculations using the STRESSexciting interface, which allows to obtain full stress tensors of crystal systems for any crystal structure. In addition, the application of STRESSexciting to the determination of stress tensor for the hexagonal Be is explicitly presented.
Table of Contents

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.
 STRESSexcitingsetup.py: Python script for generating structures at different strains.
 STRESSexcitingsubmit.sh: (Bash) shell script for running a series of exciting calculations.
 STRESSexcitinganalyze.py: Python script for fitting the energyvsstrain and CVevsstrain curves.
 STRESSexcitingclean.sh: (Bash) shell script for cleaning unnecessary files.
 STRESSexcitingresult.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.
Extra requirement: Tool for spacegroup determination
The scripts in this tutorial use the sgroup tool. If you have not done before, this tool should be downloaded and installed.
1. Set up and perform a reference calculation
In order to test the calculation of the stress tensor presented in the next sections, we set up and perform a calculation for the system we have chosen in a reference configuration. The example which is shown here is the hexagona structure of Be under an applied pressure. First, create a working directory Bereference and move inside it.
$ mkdir Bereference
$ cd Bereference
Here, copy the following lines inside a file input.xml.
<input> <title>Be: Stress calculation</title> <structure speciespath="$EXCITINGROOT/species"> <crystal scale="4.300"> <basevect> 1.00000000 0.00000000 0.00000000 </basevect> <basevect> 0.50000000 0.86602540 0.00000000 </basevect> <basevect> 0.00000000 0.00000000 1.50000000 </basevect> </crystal> <species speciesfile="Be.xml" rmt="1.95"> <atom coord="0.66666667 0.33333333 0.75000000"/> <atom coord="0.33333333 0.66666667 0.25000000"/> </species> </structure> <groundstate ngridk="6 6 4" xctype="GGA_PBE_SOL"> </groundstate> <relax/> </input>
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
$ SETUPexcitingroot.sh
Now, we procede as alredy seen in General lattice optimization in order to generate new structures and determine the corressponding pressure.
Now, you should have a file which called BM_eos.out inside the current directory. It contains all information which is related to the BirchMurnaghan fit including the hydrostatic pressure corresponding to the calculated structures. The content of this file should be similar to the following:
=== BirchMurnaghan 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.948 [GPa]
B' = 3.554
=================================================
Volume E_dftE_eos Pressure [GPa]
94.2633 +0.00000009 +18.616
98.7043 0.00000041 +10.735
103.2857 +0.00000065 +4.057
108.0004 0.00000047 1.569
112.8597 +0.00000012 6.298
As can be seen, the last 5 lines of this file show the volumes of the calculated structures (first column) and the corresponding hydrostatic pressure (third column). The first structure, (with a volume of 94.2633 Bohr^{3}) corresponds to a hydrostatic pressure of +18.616 GPa. In the following sections, this structure will be denoted as reference configuration.
We notice here that the pressure is related to the components of the stress tensor, $\sigma_{\alpha \beta}$, as:
(1)Next, we will calculate directly the stress tensor of this reference configuration and will use the equation above to (re)calculate the corresponding hydrostatic pressure.
2. Set up and perform a calculation of the stress for the reference configuration
Before continuing, we remind here that all strains considered in this tutorial are physical strains. To start the preparation of the calculation of the stress we move to the parent directory, create a new directory Bestress and move inside it.
$ cd ../
$ mkdir Bestress
$ cd Bestress
Inside this directory, we copy the input file corresponding to the reference configuration:
$ cp ../Bereference/VOL/vol_01/input.xml ./
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 STRESSexcitingsetup.py, which will produce the following output on the screen.
$ STRESSexcitingsetup.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.
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.
 input.xml, the name of the input file.
 0.002, the absolute value of the maximum strain for which we want to perform the calculation.
 21, 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 STRESSexcitingsubmit.sh, which will produce an output on the screen similar to the following.
$ STRESSexcitingsubmit.sh
++
 SCF calculation of "Dst01_01" starts 
++
Elapsed time = 0m5s
Mon 16 Jul 16:53:43 CEST 2018
++
 SCF calculation of "Dst01_02" starts 
++
Elapsed time = 0m5s
Mon 16 Jul 16:53:48 CEST 2018
++
 SCF calculation of "Dst01_03" starts 
++
...
++
 SCF calculation of "Dst02_21" starts 
++
Elapsed time = 0m5s
Mon 16 Jul 16:57:29 CEST 2018
$
After the complete run, stress components 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 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 STRESSexcitinganalyze.py allows the analysis of the dependence of the calculated derivatives of the energyvsstrain curve on
 the range of distortions included in the fitting procedure (the x axis in xmgrace plots),
 the degree of the polynomial fit used in the fitting procedure (different color curves in xmgrace plots).
The STRESSexcitinganalyze.py script is executed as follows.
$ STRESSexcitinganalyze.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 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). In a similar way, the results for the Dst02 distortion can be analyzed, too.
The script STRESSexcitinganalyze.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 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 STRESS.IN file will now contain
Dst01 0.002 1
Dst02 0.002 1
Now, you run the STRESSexcitingresult.py script.
$ STRESSexcitingresult.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
Tue Jul 17 12:38:53 2018
The hydrostatic pressure and stress tensor in kilobar (kbar)
P = 184.841019 [kbar]
150.863125 0.000000 0.000000
0.000000 150.863125 0.000000
0.000000 0.000000 252.796807
The hydrostatic pressure and stress tensor in gigapascal (GPa)
P = 18.484102 [GPa]
15.086313 0.000000 0.000000
0.000000 15.086313 0.000000
0.000000 0.000000 25.279681
From this output, we can extract a value of +18.484 GPa for the hydrostatic pressure of the reference configuration, which is in very good agreement with the value obtained in Section 1.
Exercise
 Repeat all calculations using a denser kpoint grid and a larger number of strain points.