Structure Optimization for Cubic Systems

for lithium version of Exciting


Purpose: In this tutorial you will learn how to set up and execute a structural minimization for a cubic system. Additionally, it will be explained how to visualize the relative atomic coordinates during the relaxation process.



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.

  • 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. Set up the calculations

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 relaxation of the internal coordinates of carbon (nevertheless, the same formalism can be easily applied to silicon or germanium) in the diamond structure. In fact, we already know that the optimized atomic coordinates (crystal coordinates) of this 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 tutorial, we will verify this fact. Thus, we will create a directory diamond-structure-optimization and we move inside it.

$ mkdir diamond-structure-optimization
$ cd diamond-structure-optimization

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"
               maxscl="20">
  </groundstate>
 
  <structureoptimization></structureoptimization>
 
</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

$ SETUP-excitingroot.sh

By examining the input file, we notice, in particular, that the chosen lattice constant is

    ...
    <crystal scale="6.719">
    ...

which correspond 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.

  ...
  <structureoptimization></structureoptimization>
  ...

2. 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 optimization-1 as workking directory).

$ EXECUTE-single.sh optimization-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 a PostScript file named PLOT.ps .

i) PLOT-relaxdistance.py

Python visualization tool for following 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 LABEL [YMIN YMAX]

Here, the first entry on the command line, (LABEL), 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. In particular, the choice r refers to the currently running calculation and, more generally, LABEL to the calculation already saved in the directory LABEL.

The script PLOT-relaxdistance.py accepts two more optional online entries (YMIN and YMAX) that can be specified on the calling line. Assigning numerical values to these two entries, you can set the minimum (YMIN) and the maximum (YMAX) value on the vertical axis, respectively. An example of the PostScript output of the script is the following.
c-relaxdistance.png

Here, (Δ1, Δ2, Δ3) represent the vector joining the 2 atoms in the unit cell, expressed in crystal coordinates.

ii) PLOT-status.py

Python visualization tool for the 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 LABEL
The input entry definition is the same as for the first entry for the script PLOT-relaxdistance.py. An example of the script output is the following.
c-status.png

The different line segments correspond to SCF calculations for different geometries during the relaxation procedure.

iii) PLOT-maxforce.py

Python visualization tool for the maximum amplitude of the force on the atoms during relaxation. It is executed as follows.

$ PLOT-maxforce.py LABEL
The input entry definition is the same as for the script PLOT-status.py. An example of the script output is the following.
c-maxforce.png

The red points show the calculated value at each iteration, whereas the blue line indicates the target value of the maximum amplitude of the force for stopping the relaxation.


4. Additional exercises

i) 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?

ii) 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="14"
               swidth="0.0001"
               maxscl="20">
  </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="1.00")
  • What is the physical meaning of this result?

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License