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 $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.
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.
 OPTIMIZElattice.py: Python script. A manager program which calls the setup and analyze scripts.
 OPTIMIZEsetup.py: Python script for generating structures at different volume/strains. This script is used within the lattice script.
 OPTIMIZEanalyze.py: Python script for fitting the energyvsvolume and energyvsstrain curves. This script is called by the lattice script.
 OPTIMIZEsubmit.sh: (Bash) shell script for running a series of exciting calculations.
 LATTICEparameters.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 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. 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 OPTIMIZElattice.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 $\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 the system that you want to investigate. In this tutorial, we consider as an example the calculation of energyvsvolume and energyvsstrain 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
$ SETUPexcitingroot.sh Be_opt.xml
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
$ LATTICEparameters.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 $c/a$ ratio.
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 OPTIMIZElattice.py.
$ OPTIMIZElattice.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=V_{initial}·(1+ϵ)^{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.
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 volParameters volxml vol_01 vol_02 vol_03 vol_04 vol_05
$
To execute the calculations, you have to run the script OPTIMIZEsubmit.sh. If you do so, the screen output will be similar to the following.
$ OPTIMIZEsubmit.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 OPTIMIZElattice.py python script in the current directory.
$ OPTIMIZElattice.py
>>>> Murnaghan or BirchMurnaghan EOS: [M/B]
At this point, the script is asking whether you desire to use a Murnaghan (M) or BirchMurnaghan (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
$ OPTIMIZElattice.py
>>>> Murnaghan or BirchMurnaghan 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: "BMoptimized.xml".
=====================================================================
$
Here, if you list the contents of directory, you will see the script has generated some new files
$ ls
BMoptimized.xml BM_eos.png volParameters vol_02 vol_05
BM_eos.eps INFO_VOL volxml vol_03
BM_eos.out energyvsvolume 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 BMoptimized.xml. The values lattice parameters $a, b, c, \alpha, \beta,$ and $\gamma$ at this step can be visualized using the command
$ LATTICEparameters.sh BMoptimized.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 $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 1VOL (first step, optimizing only the volume). Then, you should copy the BMoptimized.xml file to the current directory with the new name 1optimized.xml. This file will be used as the input file in the next step.
$ cd ..
$ mv VOL 1VOL
$ cp 1VOL/BMoptimized.xml 1optimized.xml
STEP2: Optimizing the $c/a$ ratio
In order to perform the next optimization step, you have to run again OPTIMIZElattice.py in the Be_OPT directory, typing entries like in the following.
$ OPTIMIZElattice.py
>>>> Please enter the exciting input file name: 1optimized.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 1optimized.xml. Now, move to the COA directory and run OPTIMIZEsubmit.sh.
$ cd COA
$ OPTIMIZEsubmit.sh
When calculations are finished, run again OPTIMIZElattice.py.
$ OPTIMIZElattice.py
=====================================================================
Optimized lattice parameter saved into the file: "coaoptimized.xml".
=====================================================================
$
In this case, the optimization is performed using a fourthorder 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 $c/a$ ratio). The optimized structure is saved in the coaoptimized.xml file. The current values of the lattice parameter are obtained once more by running the script LATTICEparameters.sh
$ LATTICEparameters.sh coaoptimized.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 $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 2COA (second optimization step, varying only $c/a$). Then, copy the coaoptimized.xml file to the current directory with the name 2optimized.xml.
$ cd ..
$ mv COA 2COA
$ cp 2COA/coaoptimized.xml 2optimized.xml
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 OPTIMIZElattice.py and using as entries values 2optimized.xml, 1, 0.03, and 11 in the given order. After having performed the calculation (running the script OPTIMIZEsubmit.sh inside the directory VOL), you run OPTIMIZElattice.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 BMoptimized.xml to the parent directory under the name 3optimized.xml.
$ cd ..
$ mv VOL 3VOL
$ cp 3VOL/BMoptimized.xml 3optimized.xml
STEP4: Optimizing again the $c/a$ ratio
Proceed in a similar way to to STEP2. Run the script OPTIMIZElattice.py using as entries values 3optimized.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 $c/a$ ratio has been completed and the new optimized structure has been saved in coaoptimized.xml.
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 4th iteration you reached an overall convergence of the lattice parameters $a$, $b$, and $c$ of about 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).
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.