for helium version of Exciting
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 the elastic constants.
|
Table of Contents
|
1. Define relevant shell variables and download scripts
Very important: Before starting, the following shell variables must be set by the user:
EXCITINGROOT = Directory where exciting has been downloaded, e.g.: /home/user/exciting .
EXCITINGRUNDIR = Directory where the user runs all the calculations (must be existing!!), e.g.: /home/user/tempdir .
EXCITINGSCRIPTS = Directory where the scripts are downloaded (not necessary if the downloaded files are copied in /usr/bin), e.g.: /home/user/scriptsdir .
The setting of these variables can be done in a bash shell by typing (from now on the symbol $ will indicate the shell prompt):
$ export EXCITINGROOT=/local_path_to_exciting_root
$ export EXCITINGRUNDIR=/local_run_directory
$ export EXCITINGSCRIPTS=/local_tutorial_scripts_bin_directory
$ export PATH=$PATH:$EXCITINGSCRIPTS
The directory names /local_path_to_exciting_root, /local_run_directory, and /local_tutorial_scripts_bin_directory are dummy names which must be explicitly changed by the user to the appropriate ones.
As a second step, it is necessary to download the scripts which are used in this tutorial. Here is a list of these scripts with a short description, click on the script name to download it:
- 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 the derivatives at zero strain of the energy-vs-strain curves.
- GNU-energy: Gnuplot visualization tool for the energy-vs-strain curves.
- GNU-status: Gnuplot visualization tool for the RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
- GNU-forces: Gnuplot visualization tool for the maximum amplitude of the force on the atoms during relaxation.
- GNU-check: Gnuplot visualization tool for the fit of the derivatives of the energy-vs-strain curves at zero strain.
- GNU-geometry: Gnuplot visualization tool for the relaxed coordinates of the atoms in the unit cell (implemented only for cells containing 2 atoms)
All the downloaded scrips must be executable, if this is not the case use, e.g.:
$ chmod u+x SETUP-elastic-strain.py
Requirements: Bash shell. Python numpy, lxml, and sys libraries. Gnuplot (visualization tools).
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. In this tutorial, we consider as an example 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 and we move inside it:
$ mkdir diamond
$ cd diamond
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:
<?xml version="1.0" encoding="UTF-8"?> <?xml-stylesheet href="inputtohtml.xsl" type="text/xsl"?> <input xsi:noNamespaceSchemaLocation="/local_path_to_exciting_root/xml/excitinginput.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" > <title>Diamond: equilibrium structure</title> <structure speciespath="/local_path_to_exciting_root/species"> <crystal scale="6.7468"> <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="10 10 10" xctype="GGAPerdew-Burke-Ernzerhof" stype= "Methfessel-Paxton 1" swidth="0.001" gmaxvr="18" lmaxvr="8" rgkmax="7.0"> </groundstate> <structureoptimization></structureoptimization> </input>
Please, remind that the input file for an exciting calculation must be always called input.xml.
Be sure to:
- set the correct path for the exciting root directory (/local_path_to_exciting_root) to the one pointing to the place where the exciting directory is placed:
... <input xsi:noNamespaceSchemaLocation="/local_path_to_exciting_root/xml/excitinginput.xsd" ... > ... <structure speciespath="/local_path_to_exciting_root/species"> ...
- have in your file the appropriate commands for performing the structure optimization: Deforming your system may change the relative positions of the atoms in the unit cell:
<structureoptimization></structureoptimization>
ii) Generation of input files for distorted structures
All the 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, which will produce the following output on the screen:
$ SETUP-elastic-strain.py
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
------------------------------------------------------------------------
Enter deformation code >>>> 0
$
In this example, the (on screen) input entries are preceded by the symbol »». The 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 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 SETUP-elastic-strain.py. Therefore, before setting up a new run please rename the workdir directory to a different name, e.g.:
$ mv workdir results-label-0
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:
$ EXECUTE-elastic-strain.sh
Running exciting for file input-01.xml ----------------------------------
...
Run completed for file input-11.xml -------------------------------------
$
After the complete run, the results of the calculation for the input file workdir/input.xml-i are contained in the subdirectory workdir/rundir-i where i is running from 01 to the total number of strain values. The data for the energy-vs-strain curves are contained in the file workdir/energy-vs-strain.
4. Post-processing: Extract energy derivatives
In order to analyse the results of the calculations, you first have to move to the working directory.
$ cd workdir
At this point, you can use the python script CHECKFIT-energy-vs-strain.py for extracting the derivatives of the energy-vs-strain curves at zero strain:
$ CHECKFIT-energy-vs-strain.py
Enter maximum strain for the interpolation >>>> 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 ==> 4291.32 [GPa]
Polynomial of order 3 ==> 4291.32 [GPa]
Polynomial of order 4 ==> 3890.14 [GPa]
Polynomial of order 5 ==> 3890.14 [GPa]
Polynomial of order 6 ==> 3900.85 [GPa]
Polynomial of order 7 ==> 3900.85 [GPa]
$
In this example, the input entries are preceded by the symbol »». The 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.
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 PostScript output file is always named gnu.ps.
i) GNU-energy
This script allows for the visualization of the energy-vs-strain curve. It is executed as follows:
$ GNU-energy
The script generates the PostScript file gnu.ps. In the following, we display the result for the example mentioned above (carbon in the diamond structure, volume-strain deformations up to 10%):
ii) GNU-check
This is a very important script 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,
- the maximum degree of the polynomial used in the fitting procedure.
The script is executed as follows:
$ GNU-check
An example of the script output is the following:
The most appropriate value for the derivative in this example is represented by the horizontal plateau (at about 3900 GPa) of the curves corresponding to the highest order fits.
iii) GNU-status
Gnuplot 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:
$ GNU-status
Enter label [r or 01,02,...] >>>> 01
$
In this example, the input entry, preceded by the symbol »» (in our example 01) represents the label of the calculation for which you would like to plot the RMS deviations. In particular, the choice r refers to the currently running calculation and, e.g., 02 to the calculation already saved in the directory workdir/rundir-02. An example of the script output is the following:
iv) GNU-forces
Gnuplot visualization tool for the maximum amplitude of the force on the atoms during relaxation (only shown if relaxation is present). It is executed as follows:
$ GNU-forces
Enter label [r or 01,02,...] >>>> 01
$
The input entry definition is the same as for the script GNU-status. An example of the script output (for the deformation label 7) is the following:
The red points show the calculated value at each iteration, whereas the blue line indicates the target value of the maximum amplitude of the force for stopping the relaxation.
v) GNU-geometry
Gnuplot visualization tool for the relaxed coordinates of the atoms in the unit cell (implemented only for cells containing 2 atoms). It is executed as follows:
$ GNU-geometry
An example of the script output (for the deformation label 7) is the following:
The vector (a1, a2, a3) represents (in terms of crystal coordinates) the relaxed vector joining the two atoms in the unit cell of the crystal . Similarly, (ref1, ref2, ref3) indicates the corresponding unrelaxed vector.
Remark on the visualization scripts
All the visualization scripts accept two on-the-command-line entries specifying the y-axis limiting values. So that, e.g., the command
$ GNU-check 3850 3950
will visualize the results for the derivative in the y-axis range between 3850 and 3950 GPa .
6. Post-processing: How to derive elastic constants
The second derivatives of the energy-vs-strain curves calculated at zero strain 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
