Simple Examples of Structure Optimization

by Pasquale Pavone for exciting fluorine

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

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

$ SETUP-excitingroot.sh

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

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.

$ 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. As an example, using the command
$ PLOT-relaxdistance.py relax-1
the PostScript od PNG output of the script should look like the following.
c-relaxdistance.png

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.

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.

$ PLOT-status.py DIRECTORYNAME

Here, the entry on the command line, DIRECTORYNAME, must be specified and represents the name of the directory where an exciting calculation has been completed or it is still running. 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

$ PLOT-status.py relax-1
the PostScript od PNG output of the script should look like 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 DIRECTORYNAME

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

$ PLOT-maxforce.py relax-1

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

c-maxforce.png

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 (see Input Reference).
"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 (see Input Reference). The optimal choice for this parameter is system dependent and should be adjusted in order to have the most efficient optimization. For further details see Input Reference.
"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" (see Input Reference) 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 (further details can be found in Input Reference). 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?
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License