General Lattice Optimization¶
by Rostam Golesorkhtabar & Pasquale Pavone for exciting sodium
(Jupyter notebook by Tejus Rohatgi, Mara Voiculescu & Martin Kuban)
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 $c/a$ ratio) and for different $c/a$ ratios (at constant volume) for Be in the hexagonal structure. The tools which are used in this tutorial are applicable for any crystal type.
0. Before Starting¶
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:
source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate
Here is a list of the scripts which are relevant for this tutorial with a short description.
optimize.setup.py: Python script for generating structures at different volume/strains.optimize.analyze.py: Python script for fitting the energy-vs-volume and energy-vs-strain curves.optimize.submit.py: Python script for running a series ofexcitingcalculations.lattice.parameters.py: Python script for extracting the values of the lattice parameters out of anyexcitinginput file.exciting2sgroup.py: Python script for converting anexcitinginput file to an input file for the program sgroup.
$\Rightarrow$ 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
$\Rightarrow$ 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
-
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 $EXCITINGROOT/tools/excitingscripts
Now, you have everything you need for starting and performing general lattice optimization.
1. Basics of lattice optimization¶
The lattice of a general crystal structure is determined by giving six lattice parameters, $a$, $b$, $c$, $\alpha$, $\beta$, and $\gamma$. The first 3 parameters are connected with the length of the 3 primitive vectors of the crystal; the last 3 are the angles between the primitive vectors. The definition of the lattice parameters can be seen in the following figure taken from the wikipedia page on lattice constants.
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 excitingscripts.optimize.lattice 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 $\alpha$ between the $b$ and $c$ axes (all other parameters are fixed).
- Minimization with respect to the angle $\beta$ between the $a$ and $c$ axes (all other parameters are fixed).
- Minimization with respect to the angle $\gamma$ between the $a$ and $b$ 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 $V$ and the $c/a$ ratio. In the next, we show how to perform the lattice optimization with respect to these two parameters.
2. 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 energy-vs-volume and energy-vs-strain curves for Be in the hexagonal structure. Therefore, we create a directory Be_OPT.
%%bash
mkdir -p 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
%%bash
cd Be_OPT
python3 -m excitingscripts.setup.excitingroot -i Be_opt.xml -o Be_opt.xml
cd ..
The values of the initial lattice parameters, $a$, $b$, $c$, $\alpha$, $\beta$, and $\gamma$, corresponding to the initial input file Be_opt.xml can be extracted using the following command
%%bash
cd Be_OPT
python3 -m excitingscripts.lattice.parameters -i Be_opt.xml
cd ../
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
- The values of the parameter
ngridkin 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 ofngridkand 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 $c/a$ ratio.
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 excitingscripts.optimize.setup.
%%bash
cd Be_OPT
python3 -m excitingscripts.optimize.setup -i Be_opt.xml --opt 1 --s-max 0.03 -n 5
cd ../
The entry values must be inlcuded when running the script. In this case, the 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, $\epsilon_{max}$, for which we want to perform the calculation. Notice that in this case the physical strain $\epsilon$ is defined by the relationship $V=V_{initial} \cdot (1+\epsilon)^3$, where $V_{initial}$ 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.
If you do so, the screen output will be similar to the following:
Number and name of space group: 194 (P 63/m m c)
Hexagonal I structure in the Laue classification.
Parameters to optimize:
1 ... volume
2 ... c/a ratio with constant volume
1 -> Parameter to optimize is volume
Maximum physical strain is 0.03
Number of distorted structures is 5
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:
INFO_VOL
source.xml
vol-Parameters
vol_1
vol_2
vol_3
vol_4
vol_5
%%bash
cd Be_OPT/VOL
ls
cd ../../
To execute the calculations, you have to run the script excitingscripts.optimize.submit. If you do so, the screen output will be similar to the following.
%%bash
cd Be_OPT/VOL
python3 -m excitingscripts.optimize.submit
cd ../../
When calculations finished to run, results can be analyzed. In order to do this, you have to run again the optimize.analyze.py python script in the current directory.
%%bash
cd Be_OPT/VOL
python3 -m excitingscripts.optimize.analyze --eos b
cd ../../
At this point, you need to specify 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 as an entry on the command line, you will get as a result
=====================================================================
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
%%bash
cd Be_OPT/VOL
ls
cd ../../
You should see something like:
BM_eos.out
BM_eos.png
INFO_VOL
energy-vs-volume
source.xml
vol-Parameters
vol-optimized.xml
vol_1
vol_2
vol_3
vol_4
vol_5
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, a 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).
- 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 outside 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 ofrgkmax, 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 $a$, $b$, $c$, $\alpha$, $\beta$, and $\gamma$ at this step can be visualized using the command:
%%bash
cd Be_OPT/VOL
python3 -m excitingscripts.lattice.parameters -i BM-optimized.xml
cd ../../
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 a, b, and c is the same. If you are interested to check how accurate the calculated equilibrium parameters at this step are, you can find more information here.
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.
%%bash
cd Be_OPT
mv VOL 1-VOL
cp 1-VOL/BM-optimized.xml 1-optimized.xml
cd ../
%%bash
cd Be_OPT
python3 -m excitingscripts.optimize.setup -i 1-optimized.xml --opt 2 --s-max 0.03 -n 5
cd ../
Notice that in this case the physical strain $\epsilon$ is defined by the relationship $(c/a)=(c/a)_{initial}\cdot (1+\epsilon)$, where $(c/a)_{initial}$ corresponds to the unit cell defined in 1-optimized.xml. Now, move to the COA directory and run optimize.submit.py.
%%bash
cd Be_OPT/COA
python3 -m excitingscripts.optimize.submit
cd ../../
When calculations are finished, run again optimize.analyze.py.
%%bash
cd Be_OPT/COA
python3 -m excitingscripts.optimize.analyze
cd ../../
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 should look like the following.
At this point, you have completed the second optimization step (by varying only the $c/a$ ratio). The optimized structure is saved in the coa-optimized.xml file. The current values of the lattice parameter are obtained once more by running the script lattice.parameters.py
%%bash
cd Be_OPT/COA
python3 -m excitingscripts.lattice.parameters -i coa-optimized.xml
cd ../../
As an output on the screen you should get
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 $c/a$ ratio at fixed $V$, the values of $a$,$b$ and $c$ differ from the ones at the previous step. Now, you should move out to the parent directory, rename COA to 2-COA (second optimization step, varying only $c/a$). Then, copy the coa-optimized.xml file to the current directory with the name 2-optimized.xml.
%%bash
cd Be_OPT
mv COA 2-COA
cp 2-COA/coa-optimized.xml 2-optimized.xml
cd ../
STEP3: Optimizing again the volume
Notice: Due to the fact that the volume and $c/a$ ratio have been already optimized once, we can increase the number of calculations at each step in order to achieve a better fitting accuracy.
Repeat now the procedure already explained in STEP1, running the script optimize.setup.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.py inside the directory VOL), you run optimize.analyze.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 vol-optimized.xml to the parent directory under the name 3-optimized.xml.
At this point, the second optimization of the c/a ratio has been completed and the new optimized structure has been saved in coa-optimized.xml.
Reaching convergence
In the following table you find a summary of the results for the lattice parameters of the first 4 optimization steps.
| Step | OPT | amin [Bohr] | bmin [Bohr] | cmin [Bohr] | αmin [degrees] | βmin [degrees] | γmin [degrees] |
|---|---|---|---|---|---|---|---|
| 0 | 4.30000 | 4.30000 | 6.45000 | 90.00000 | 90.00000 | 120.00000 | |
| 1 | $V$ | 4.34538 | 4.34538 | 6.51806 | 90.00000 | 90.00000 | 120.00000 |
| 2 | $c/a$ | 4.32416 | 4.32416 | 6.58219 | 90.00000 | 90.00000 | 120.00000 |
| 3 | $V$ | 4.32127 | 4.32127 | 6.57779 | 90.00000 | 90.00000 | 120.00000 |
| 4 | $c/a$ | 4.32088 | 4.32088 | 6.57897 | 90.00000 | 90.00000 | 120.00000 |
As you can see from the previous table, at the 4-th iteration you reached an overall convergence of the lattice parameters $a$, $b$, and $c$ of about $\mathbf{1 \times 10^{-3}}$ Bohr. If this result corresponds to the desired accuracy, you can stop the optimization procedure. Otherwise, you proceed with the next step and, using the new results, you check again the convergence behaviour of the equilibrium 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 $c/a$ ratio starting from the initial input file Be_opt.xml, and the parameters ngridk and rgkmax are set as follows
<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).
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
relaxis present in input.xml; - perform a single calculation of
exciting.