Simple Examples of Structure Optimization¶

by Pasquale Pavone for exciting sodium

(Jupyter notebook by Tejus Rohatgi, Mara Voiculescu & Martin Kuban)

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.

Table of Contents

0. Before Starting

1. Diamond: Set up the calculations

2. Diamond: Perform the calculations

3. Analysis of the results: Visualization tools

4. Diamond: Additional exercises

0. Before Starting¶

Read the following paragraphs before starting with the rest of this tutorial!

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source $EXCITINGROOT/tools/excitingjupyter/venv/excitingvenv/bin/activate

Here is a list of the scripts which are relevant for this tutorial with a short description.

execute.single.py: Python 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.

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

(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.

%%bash
mkdir -p 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
      taubfgs="0.3" />
 
</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

%%bash
cd diamond-relax
python3 -m excitingscripts.setup.excitingroot
cd ..

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
      taubfgs="0.3" />
...

The meaning of the attribute taubfgs will be given in Section 4.

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.py as shown in the following.

%%bash
cd diamond-relax
python3 -m excitingscripts.execute.single 
cd ..

3. Analysis of the results: Visualization tools¶

All the scripts mentioned here produce as output a PNG (PLOT.png) file.

i) plot.relaxdistance.py

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.

%%bash
cd diamond-relax
python3 -m excitingscripts.plot.relaxdistance
cd ../

the PNG output of the script should look like the following.

No description has been provided for this image

Here, (Δ1, Δ2, Δ3) represents the 3D vector joining the 2 atoms in the unit cell, expressed in crystal coordinates. Notice that in this example, due of the symmetry of the system and of the initial configuration, the three components (Δ1, Δ2, Δ3) are identical, resulting in a single visible curve.

The following optional entries, which can be specified on the command are also supported by this script.

python3 -m excitingscripts.plot.relaxdistance DIRECTORYNAME ATOM1 ATOM2 YMIN YMAX

Here, the first entry on the command line, (DIRECTORYNAME),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 default value is the current directory where the script is executed.

The next entries, ATOM1, ATOM2, YMIN, and YMAX, 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. As an example, using the command

ii) plot.status.py

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.

%%bash
cd diamond-relax
python3 -m excitingscripts.plot.status 
cd ..

Here, the entry on the command line, the DIRECTORYNAME, can be specified and represents the name of the directory where an exciting calculation has been completed or it is still running. The default value is the current directory where the script is executed. The plot.status.py script is particularly useful in the latter case because one can follow "live" the convergence of the calculation. As an example, using the command

the PNG output of the script should look like the following.

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.

%%bash
cd diamond-relax
python3 -m excitingscripts.plot.maxforce 
cd ..

The input entry definition is the same as for the script plot.status.py. As an example, using the command

the PNG output of the script should look like the following.

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.

Method | Description
:--- | :--- "bfgs" | (default) Limited memory algorithm for bound constrained optimization, see Byrd, et al. SIAM J. Scientific Computing 16, 1190 (1995). This method requires high accuracy for the determination of the total energy in dependence of the maximum allowed for the residual atomic force. For this reason, the default value of the attribute epsengy (defined in groundstate) is decreased in order to be at least equal to the value of the attribute epsforce of the relax element multiplied by a factor 0.02. The range of variation of each degree of freedom is set by the attribute taubfgs. | "newton" |Simple (Newton-like) method. At each step of a structural optimization run, each atom is displaced proportionally to the force acting on it. The proportionality constant is given by the value of the attribute taunewton. The optimal choice for this parameter is system dependent and should be adjusted in order to have the most efficient optimization. | "harmonic" |Method based on the combination of the "newton" method and the harmonic approximation. Contrary to "newton", all cartesian components of the position vector of each atom are treated independently. At the each optimization step each cartesian component of the position vector of each atom is updated using the same algorithm as in "newton" ("newton step") unless the "harmonic condition" is fulfilled. In this case ("harmonic step"), atomic positions are updated according to a linear interpolation of the forces at the last two optimization steps. The "harmonic" method is of general validity and it is particularly efficient when the atomic configuration is close to the optimized one and the internal degrees of freedom are weakly coupled. |

Change the input file by explicitly changing the attribute method and adding taunewton inside relax and choose for method a value different from the default, e.g.

...
   <relax
      method="harmonic"
      taunewton="0.2" />
...

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 most efficient method for this example for diamond?
What does it change if you execute an input file with the following parameters?

...
   <relax
      method="harmonic"
      taunewton="0.4" />
...

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>
...
   <relax
      taubfgs="0.3" />
...

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?