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.
- 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.
- OPTIMIZE-clean.sh: (Bash) shell script for cleaning unnecessary files.
- exciting2sgroup.xsl: xsl 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, $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 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 parameter are fixed).
- Minimization with respect to the $c/a$ ratio (all other parameter are fixed).
- Minimization with respect to the angle $\alpha$ between the $b$ and $c$ axes (all other parameter are fixed).
- Minimization with respect to the angle $\beta$ between the $a$ and $c$ axes (all other parameter are fixed).
- Minimization with respect to the angle $\gamma$ between the $a$ and $b$ axes (all other parameter 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 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 exciting input file for hexagonal Be with the name Be_opt.xml name. This file could 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>
This file can be saved with any name. In this tutorial is not necessary to rename the exciting input file as input.xml, because this file is the input of the script OPTIMIZE-lattice.py and not of exciting itself. Please, notice that the input file for a direct exciting calculation must be always called input.xml.
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 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 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 for which we want to perform the calculation.
- 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 = 0m3s SCF calculation of "vol_02" starts -------------------------------- Elapsed time = 0m3s ... SCF calculation of "vol_05" starts -------------------------------- Elapsed time = 0m3s $
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.04 Final parameters: E_min = -29.36096 [Ha] V_min = 106.5871 [Bohr^3] B_0 = 121.95 [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).
- 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. 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-VOL.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-VOL.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-VOL.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 $
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 plot (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 coa-optimized.xml file. Similar to the previous step, 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-COA.xml.
$ cd .. $ mv COA 2-COA $ cp 2-COA/coa-optimized.xml 2-COA.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 choose smaller range of distortion for the next steps.
Repeat now the procedure already explained in STEP1, running the script OPTIMIZE-lattice.py and using as entries values 2-COA.xml, 1, 0.005, 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-VOL.xml.
$ cd .. $ mv VOL 3-VOL $ cp 3-VOL/BM-optimized.xml 3-VOL.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-VOL.xml, 2, 0.005, 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 coa-optimized.xml.
In the following table you find a summary of the result of the first 4 optimization steps.
As you can see from the previous table, at the 4-th iteration you reached the following convergence.
- The equilibrium volume is converged within 2$\times$10-1 Bohr3.
- The c/a ratio is converged within 4$\times$10-4.
- The energy at the minimum seems to be converged within 3$\times$10-4 mHa. Indeed, such a small value should be considered an artifact of the optimization procedure, which assumes that the calculated total energies are exact. However, the accuracy in the determination of the minimum energy cannot be smaller than the accuracy of the total energy in a single SCF calculation. For the calculations performed in this tutorial, total energies are calculated with the default value of the accuracy, i.e., 10-3 mHa.
If these results correspond 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 behavior of the equilibrium parameters.
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 input file Be_opt.xml to input.xml;
- inside input.xml, set the cell parameters to the optimized ones;
- be sure that the element relax is present in input.xml;
- perform a single calculation of exciting.