by Rostam Golesorkhtabar & Pasquale Pavone for exciting oxygen
Purpose: In this tutorial, you will learn how to optimize a general crystal structure. An explicit example is given for hexagonal structures. Here, you will set up and execute a series of calculations for different volumes (at constant
|
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.
- OPTIMIZE-lattice.py: Python script. A manager program which calls the setup and analyze scripts.
- OPTIMIZE-setup.py: Python script for generating structures at different volume/strains. This script is used within the lattice script.
- OPTIMIZE-analyze.py: Python script for fitting the energy-vs-volume and energy-vs-strain curves. This script is called by the lattice script.
- OPTIMIZE-submit.sh: (Bash) shell script for running a series of exciting calculations.
- LATTICE-parameters.sh: (Bash) shell script for extracting the values of the lattice parameters out of any exciting input file.
- exciting2sgroup.py: Python script for converting an exciting input file to an input file for the program sgroup.
From now on the symbol $ will indicate the shell prompt.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.
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.
1. Basics of lattice optimization
The lattice of a general crystal structure is determined by giving six lattice parameters,
In order to perform the optimization of the energy of a crystal with respect to all lattice parameters, we can use an iterative cyclic procedure where in turns one parameter is varied and the remaining 5 are kept fixed. The procedure can be repeated until the obtained equilibrium parameters do not vary anymore within a desired accuracy. In particular, the minimization procedure as performed by the script OPTIMIZE-lattice.py will contain the following cycles.
- Minimization with respect to the volume
V (by applying an isotropic strain). - Minimization with respect to the
b/a ratio (all other parameters are fixed). - Minimization with respect to the
c/a ratio (all other parameters are fixed). - Minimization with respect to the angle
α between theb andc axes (all other parameters are fixed). - Minimization with respect to the angle
β between thea andc axes (all other parameters are fixed). - Minimization with respect to the angle
γ between thea andb axes (all other parameters are fixed).
In the example reported in this tutorial, we consider a crystal with
hexagonal crystal structure. In this case there are only two free
parameters, the volume
2. Preparation of the input file
The first step is to create a directory for the system that you want to investigate. In this tutorial, we consider as an example the calculation of energy-vs-volume and energy-vs-strain curves for Be in the hexagonal structure. Therefore, we create a directory Be_OPT and we move inside it.
$ mkdir Be_OPT
$ cd Be_OPT
Inside this directory, we create or copy an input file for hexagonal Be with the name Be_opt.xml. In this tutorial is not necessary to rename the input file as input.xml, because this file is the input of a specific script used in this tutorial and not of exciting itself. Please, notice that the input file for a direct exciting calculation must be always called input.xml.
This file should look like the following.
<input> <title>Be: Lattice optimization</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
$ SETUP-excitingroot.sh Be_opt.xml
The values of the initial lattice parameters,
$ LATTICE-parameters.sh Be_opt.xml
As an output on the screen you should get
Bravais lattice: Hexagonal
a b c alpha beta gamma
4.30000 4.30000 6.45000 90.00000 90.00000 120.00000
Please note:
- The values of the parameter ngridk in the input file have been chosen in order to allow for a fast calculation. Of course, the complete optimization procedure presented in this tutorial should be repeated for more accurate values of ngridk and the other computational parameters up to the desidered accuracy in the resulting lattice constants.
3. Lattice optimization of hexagonal crystals
In the next, we illustrate the iterative procedure for performing the
optimization of the crystal structure of hexagonal Beryllium. The only
relevant parameters in this case are the volume of the unit cell and the
STEP1: Optimizing the volume
At the first step, we optimize the energy with respect to the volume. In order to generate input files for a series of volumes you have to use the script OPTIMIZE-lattice.py.
$ OPTIMIZE-lattice.py
>>>> Please enter the exciting input file name: Be_opt.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]: 5
The number of the distorted structures is 5
$
Entry values must be typed on the screen when requested. In this case, entries are the following.
- Be_opt.xml, the name of the input file you have created.
- 1, the type of optimization you choose.
- 0.03, the absolute value of the maximum strain, ϵmax, for which we want to perform the calculation. Notice that in this case the physical strain ϵ is defined by the relationship V=Vinitial·(1+ϵ)3, where Vinitial is the volume of the unit cell defined in Be_opt.xml.
- 5, 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. If you list the contents of directory, you should see something like
$ cd VOL
$ ls
INFO_VOL vol-Parameters vol-xml vol_01 vol_02 vol_03 vol_04 vol_05
$
To execute the calculations, you have to run the script OPTIMIZE-submit.sh. If you do so, the screen output will be similar to the following.
$ OPTIMIZE-submit.sh
SCF calculation of "vol_01" starts --------------------------------
Elapsed time = 0m5s
SCF calculation of "vol_02" starts --------------------------------
Elapsed time = 0m5s
...
SCF calculation of "vol_05" starts --------------------------------
Elapsed time = 0m5s
$
When calculations finished to run, results can be analyzed. In order to do this, you have to run again the OPTIMIZE-lattice.py python script in the current directory.
$ OPTIMIZE-lattice.py
>>>> Murnaghan or Birch-Murnaghan EOS: [M/B]
At this point, the script is asking whether you desire to use a Murnaghan (M) or Birch-Murnaghan (B) equation of state for extracting equilibrium parameters such as the equilibrium volume and bulk modulus.
Therefore, if you type B or b, you will get as a result
$ OPTIMIZE-lattice.py
>>>> Murnaghan or Birch-Murnaghan EOS: [M/B] b
=====================================================================
Fit accuracy:
Log(Final residue in [Ha]): -6.03
Final parameters:
E_min = -29.3609602 [Ha]
V_min = 106.587 [Bohr^3]
B_0 = 121.952 [GPa]
B' = 3.553
Optimized lattice parameter saved into the file: "BM-optimized.xml".
=====================================================================
$
Here, if you list the contents of directory, you will see the script has generated some new files
$ ls
BM-optimized.xml BM_eos.png vol-Parameters vol_02 vol_05
BM_eos.eps INFO_VOL vol-xml vol_03
BM_eos.out energy-vs-volume vol_01 vol_04
$
The BM_eos.out file, for instance, contains all final fit parameters and pressures which are calculated from equation of state fitting procedure. Also, the script generates a plot, PostScript (BM_eos.eps) or PNG (BM_eos.png) file, which looks like the following.
On this plot, you can also find the optimized values of the parameters appearing in the equation of state (minimum energy, equilibrium volume, bulk modulus, and bulk modulus pressure derivative).
Please note:
- The bulk modulus and bulk modulus pressure derivative which are derived here have to be interpreted only as fitting parameters. This is due to the fact that in this tutorial we are changing the volume, however, we are keeping all other lattice parameter constant. In order to obtain the real physical bulk modulus and bulk modulus pressure derivative, one has to fully optimize at each given volume with respect to all other lattice and internal parameters (e.g., in the case of hexagonal beryllium, with respect to the
c/a ratio, too).- The visual analysis of the plots is very important. The user should always check them at each step. In particular, if the minimum lies ouside the displayed region, the calculation should be restarted with more appropriate values of the initial parameters (e.g., volume or
c/a ratio). The optimal situation is when the minimum of the energy curves is located in the middle of the investigated region.- If the difference in energy between the calculated points and the fit is larger than the final required accuracy in the energy, the calculation should be restarted with more appropriate computational parameters. In particular, one should consider the number of k points (ngridk), the value of rgkmax, and the accuracy in the calculated total energy (epsengy).
A file corresponding to an exciting input file for the optimized geometry is created with the name BM-optimized.xml. The values lattice parameters
$ LATTICE-parameters.sh BM-optimized.xml
As an output on the screen you should get
Bravais lattice: Hexagonal
a b c alpha beta gamma
4.34538 4.34538 6.51806 90.00000 90.00000 120.00000
Due to the fact that in this step the volume only has been optimized,
the relative variation with respect to the initial values of the three
lattice parameters
At this point, you have performed the first optimization step by varying only the volume. In order to be prepared for the next step, you should move now to the parent directory and rename the VOL directory to 1-VOL (first step, optimizing only the volume). Then, you should copy the BM-optimized.xml file to the current directory with the new name 1-optimized.xml. This file will be used as the input file in the next step.
$ cd ..
$ mv VOL 1-VOL
$ cp 1-VOL/BM-optimized.xml 1-optimized.xml
STEP2: Optimizing the c/a ratio
In order to perform the next optimization step, you have to run again OPTIMIZE-lattice.py in the Be_OPT directory, typing entries like in the following.
$ OPTIMIZE-lattice.py
>>>> Please enter the exciting input file name: 1-optimized.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': 2
>>>> 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]: 5
The number of the distorted structures is 5
$
Notice that in this case the physical strain ϵ is defined by the relationship (c/a)=(c/a)initial·(1+ϵ), where (c/a)initial corresponds to the unit cell defined in 1-optimized.xml. Now, move to the COA directory and run OPTIMIZE-submit.sh.
$ cd COA
$ OPTIMIZE-submit.sh
When calculations are finished, run again OPTIMIZE-lattice.py.
$ OPTIMIZE-lattice.py
=====================================================================
Optimized lattice parameter saved into the file: "coa-optimized.xml".
=====================================================================
$
In this case, the optimization is performed using a fourth-order polynomial fit for calculating the minimum energy and the corresponding strain. The resulting output coa.png (also available as PostScript file as coa.eps) should look like the following.
At this point, you have completed the second optimization step (by varying only the
$ LATTICE-parameters.sh coa-optimized.xml
Bravais lattice: Hexagonal
a b c alpha beta gamma
4.32416 4.32416 6.58219 90.00000 90.00000 120.00000
$
Here, you can notice that by changing the
$ cd ..
$ mv COA 2-COA
$ cp 2-COA/coa-optimized.xml 2-optimized.xml
STEP3: Optimizing again the volume
Notice: Due to the fact that the volume and
Repeat now the procedure already explained in STEP1, running the script OPTIMIZE-lattice.py and using as entries values 2-optimized.xml, 1, 0.03, and 11 in the given order. After having performed the calculation (running the script OPTIMIZE-submit.sh inside the directory VOL), you run OPTIMIZE-lattice.py and get the following plot.
At this point, you have optimized the volume for the second time. Follow the last part of STEP1 and copy the file BM-optimized.xml to the parent directory under the name 3-optimized.xml.
$ cd ..
$ mv VOL 3-VOL
$ cp 3-VOL/BM-optimized.xml 3-optimized.xml
STEP4: Optimizing again the c/a ratio
Proceed in a similar way to to STEP2. Run the script OPTIMIZE-lattice.py using as entries values 3-optimized.xml, 2, 0.03, and 11 in the given order. Using the same procedure as in the previous steps, you will end up with the following plot.
At this point, the second optimization of the
Reaching convergence
In the following table you find a summary of the results for the lattice parameters of the first 4 optimization steps.
|
As you can see from the previous table, at the 4-th iteration you reached an overall convergence of the lattice parameters
In some case, it may happen that after some new optimization steps the values of the lattice parameters oscillate without reaching convergence. Then, you should improve the accuracy in the energy of the single calculations by correspondingly changing computational parameters such as ngridk, rgkmax, and epsengy.
Note on the accuracy of the calculations
The quality of the results you can obtain for the structural
optimization of a crystal using the procedure illustrated in this
tutorial strongly depends on the accuracy in the determination of the
energy of each single calculations. For instance, if one performs a
minimization step of the
<groundstate ... ngridk="3 3 2" rgkmax="4.0" ...> </groundstate>
a result is obtained which is similar to the one illustrated in the next plot, where no global minimum can be clearly identified.
In order to obtain reasonable results, a general rule is that the accuracy of the energy for each calculations should be much smaller than the energy difference between two different calculations (two different points in the above plot).
Final step: Internal optimization
The procedure shown in this tutorial optimizes the cell parameters of a crystal. If the crystal symmetry allows for internal relaxation of the atomic positions (i.e., without changing the cell parameters) a single further calculation is necessary to optimize all degrees of freedom of the crystal. One has to
- copy the last optimized input xml file to input.xml;
- be sure that the element relax is present in input.xml;
- perform a single calculation of exciting.
http://exciting.wikidot.com/oxygen-general-lattice-optimization
Hello, my name is Dr Abderrahmane Reggad. I'm a 53 year old and I am a researcher in materials' sciences. 
0 Comments