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 shell variables and download scripts
Read the following paragraphs before starting with the rest of this tutorial!
Before starting, be sure that relevant shell variables are already defined and that the excitingscripts directory has already been downloaded, as specified in Tutorial scripts and environment variables. 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.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.
From now on the symbol $ will indicate the shell prompt.
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. 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
- Website: http://cpc.cs.qub.ac.uk/summaries/adon
- Download sgroup here (open the link, push the Download botton, and save the file containg the sgroup package)
After the download, you will get a tar.gz file, go to the directory where you saved this file and execute the following commands.
$ tar xfvp DownloadedFile.tar.gz $ cd SpaceGroups $ make $ cp sgroup $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.
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.
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.3100"> <basevect> 1.00000000 0.00000000 0.00000000 </basevect> <basevect> -0.50000000 0.86602540 0.00000000 </basevect> <basevect> 0.00000000 0.00000000 1.51000000 </basevect> </crystal> <species speciesfile="Be.xml" rmt="1.45"> <atom coord="0.66666666 0.33333334 0.75000000"/> <atom coord="0.33333334 0.66666666 0.25000000"/> </species> </structure> <groundstate ngridk="6 6 4" gmaxvr="14" maxscl="20"> </groundstate> <structureoptimization/> </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
3. Lattice optimization of hexagonal crystals
In the next, we illustrate the iterative procedure for performing the optimization of the crystal stucture 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.
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.010, 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
To execute the calculations, you have to run the script OPTIMIZE-submit.sh. If you do so, you see the following messages on the screen.
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.
At this point, the script is asking whether you desire to use a Murnaghan (M) or Birch-Murnaghan (B) equation of state for extracting the equilibrium parameters such as the equilibrium volume and the bulk modulus. For more details about the equation of state see here.
Therefore, if you type B or b, you will get as a result
Moreover, the script generates a plot (PostScript file vol.ps) 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. They do not coincide with the "exact" bulk modulus and bulk modulus pressure derivative of the crystal. Indeed, these "exact" values should be obtained by fitting, with an equation of state (Birch or Murnaghan), the function E=E(V), where for each volume V, E(V) is the energy obtained by optimizing, at that given V, all other lattice and internal parameters.
A file corresponding to an exciting input file for the optimized geometry is created with the name vol-optimized.xml. If you are interested to check how accurate are the calculated equilibrium parameters at this step, you can find more information here.
Now, you can delete unnecessary files by executing the command
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 vol-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.
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.
Now, move to the COA directory and run OPTIMIZE-submit.sh.
When calculations are finished, run again OPTIMIZE-lattice.py.
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.ps) should look like the following.
Delete unnecessary files.
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.
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 using the following entries.
After having performed the calculation, 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 vol-optimized.xml to the parent directory under the name 3-VOL.xml.
STEP4: Optimizing again the c/a ratio
Repeat the procedure in STEP2 with the following input.
In a similar way to STEP2, 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 10-1 [Bohr3]
- The c/a ratio is converged within 2$\times$10-4.
- The energy at the minimum is converged within 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 behaviour of the equilibrium parameters.