Purpose: In this tutorial you will learn how to set up and execute the optimization of the relative position of the atoms in the unit cell of a crystal. Examples will be presented for the diamond crystal. Additionally, it will be explained how to visualize the relative atomic coordinates during the relaxation process.
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.
- EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
- PLOT-relaxdistance.py: Python visualization tool for following relative atomic coordinates of atoms during the relaxation process.
- PLOT-status.py: Python visualization tool for RMS deviations of the SCF potential as a function of the iteration number during the SCF loop.
- PLOT-maxforce.py: Python visualization tool for the maximum amplitude of the force on the atoms during relaxation.
From now on the symbol $ will indicate the shell prompt.
Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.
1. Diamond: Set up the calculations
As a first example, we consider the relaxation of the internal coordinates of the carbon atoms in the diamond structure (fcc with two atoms in the unit cell). The same formalism can be easily applied to other materials which crystallize in the diamond structure, e.g., silicon, germanium, and α-tin. In fact, we already know that the optimized atomic coordinates (in crystal coordinates) of the two atoms in the diamond structure are (apart from a common translation)
- (0, 0, 0) for the first atom,
- (0.25, 0,25, 0.25) for the second atom.
In this example, we will verify this fact. The first step is to create a directory for each system that you want to investigate. Thus, we will create a directory diamond-relax and we move inside it.
$ mkdir diamond-relax $ cd diamond-relax
Inside this directory, we create the file input.xml. This file could look like the following.
<input> <title>Diamond: Internal-structure optimization</title> <structure speciespath="$EXCITINGROOT/species"> <crystal scale="6.719"> <basevect> 0.5 0.5 0.0 </basevect> <basevect> 0.5 0.0 0.5 </basevect> <basevect> 0.0 0.5 0.5 </basevect> </crystal> <species speciesfile="C.xml" rmt="1.20"> <atom coord="0.00 0.00 0.00" /> <atom coord="0.30 0.30 0.30" /> </species> </structure> <groundstate ngridk="6 6 6" swidth="0.0001" xctype="GGA_PBE_SOL"> </groundstate> <relax></relax> </input>
Please, remember that the input file for an exciting calculation must always be called input.xml.
N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command
By examining the input file, we notice, in particular, that the chosen lattice constant is
... <crystal scale="6.719"> ...
which corresponds to the equilibrium one (calculated with ngridk = "10 10 10"). Furthermore, the initial atomic positions are set according to
... <species speciesfile="C.xml" rmt="1.20"> <atom coord="0.00 0.00 0.00" /> <atom coord="0.30 0.30 0.30" /> </species> ...
Finally, in order to perform the structure optimization, we must insert in the input file, after the SCF calculation performed in groundstate, the following line.
... <relax></relax> ...
2. Diamond: Perform the calculation
To execute the calculation with the input file we created in the previous section, you have to run the script EXECUTE-single.sh. Notice that the script EXECUTE-single.sh always generates a working directory containing results of the current calculations. The directory name can be specified by adding the name in the command line.
$ EXECUTE-single.sh DIRECTORYNAME
If no name is given, the script use the default name rundir-00. Very important: The working directory is overwritten each time you execute the script EXECUTE-single.sh. Therefore, choose different names for different calculations.
The script EXECUTE-single.sh produces the following output on the screen (using relax-1 as working directory).
$ EXECUTE-single.sh relax-1 ===> Output directory is "relax-1" <=== Running exciting for file input.xml ------------------------------------- ... Run completed for file input.xml ---------------------------------------- $
3. Analysis of the results: Visualization tools
All the scripts mentioned here produce as output both a PostScript (PLOT.ps) and a PNG (PLOT.png) file.
Python visualization tool to follow the relative atomic coordinates of the two atoms in the diamond unit cell during the relaxation process. It is executed as follows.
$ PLOT-relaxdistance.py DIRECTORYNAME [ATOM1 ATOM2 YMIN YMAX]
Here, the first entry on the command line, (DIRECTORYNAME), must be specified and represents the label (i.e., the name of the working directory) for which you would like to follow the relaxation of the relative atomic coordinates.
The script PLOT-relaxdistance.py accepts four more optional online entries indicated above between the square brackets (notice that, in order to execute the script, the square brackets must be omitted). The optional entries, ATOM1, ATOM2, YMIN, and YMAX, which can be specified on the calling line have the following meaning:
- The entries ATOM1 and ATOM2 may assume integer values from 1 to the number of atoms in the unit cell. The script will plot the relative position of the atom with index ATOM2 with respect to the position of the atom with index ATOM1. Default values are 1 and 2 for ATOM1 and ATOM2, respectively.
- Assigning numerical values to the entries YMIN and YMAX, you can set the minimum (YMIN) and the maximum (YMAX) value on the vertical axis, respectively.
Here, (Δ1, Δ2, Δ3) represents the 3D vector joining the 2 atoms in the unit cell, expressed in crystal coordinates.
Python visualization tool for the root mean square (RMS) deviations of the effective SCF potential as a function of the iteration number during the SCF loop. It is executed as follows.
$ PLOT-status.py DIRECTORYNAME
The different line segments correspond to SCF calculations for different geometries during the relaxation procedure.
Python visualization tool for the maximum amplitude of the force on the atoms during relaxation. It is executed as follows.
$ PLOT-maxforce.py DIRECTORYNAME
The red points show the calculated value at each iteration, whereas the blue line indicates the target value (set by the attribute epsforce inside relax) for the maximum amplitude of the force for stopping the relaxation.
4. Diamond: Additional exercises
i) Change the relaxation method
exciting can perform the optimization of the position of the atoms in the unit cell using three different methods. The application of one of the three available methods can be triggered by setting the attribute method (see Input Reference) of the element relax to one of the following choices.
Change the input file by explicitly adding the attributes method and taunewton inside relax and choose for method a value different from the default, e.g.
... <relax method="harmonic" taunewton="0.2"> </relax> ...
- How many optimization runs are needed for reaching the convergence for the different methods?
- How much CPU time is required by the different methods?
- What is the influence on the results (number of optimization runs and CPU time) of the value of the attribute taunewton?
- What is the most efficient method for this example for diamond?
ii) Change the initial position
Change the initial position of the second atom in the unit cell according to the following (non symmetrical choice).
... <species speciesfile="C.xml" rmt="1.20"> <atom coord="0.00 0.00 0.00" /> <atom coord="0.30 0.20 0.37" /> </species> ...
Perform a new calculation.
- Why do you expect a longer execution time in this case?
iii) Change the lattice parameter
Set the cubic lattice constant to a smaller value, as in the following example.
... <crystal scale="5.0"> ...
Due to the smaller volume (and distances), set the muffin-tin radius of the carbon atoms to a smaller value.
... <species speciesfile="C.xml" rmt="1.05"> ...
Set the initial coordinates of the second atom according to
... <atom coord="0.00 0.00 0.00" /> <atom coord="0.45 0.45 0.45" /> ...
Include the new parameter gmaxvr in the groundstate element (this addition is needed due to the strongly reduced value of the lattice constants which is used in this calculations, see Input Reference).
... <groundstate ngridk="6 6 6" gmaxvr="15" swidth="0.0001"> </groundstate> ...
Perform a new calculation.
- What are the equilibrium relative atomic coordinates in this case?
- What happens if the initial coordinates of the second atom are set to (0.26, 0.26, 0.26)? (Hint: set rmt = "0.95")
- What is the physical meaning of this result?