Structural optimization using forces and stress¶

This tutorial provides an overview of the methods available in Siesta for structural optimization, which are quite varied in their nature, scope, and sophistication.

Skip the background information - Take me to the exercises!

background Background Information¶

Geometry optimization¶

Note

Be aware that in this context, the words geometry/structure, and optimization/minimization are often used interchangeably. The same thing happens for minimization algorithm/minimizer/optimizer.

One of the most common tasks in computational material science is to optimize a structure. By optimization, we mean to find the structure that has the lowest possible energy; i.e., the one that is the most stable. Since this should be an energy minimum, and the forces are opposite to the energy gradients, it also means that the forces on the atoms should be zero.

There are several ways to do this, but most of them involve using the forces over the atoms to move them in the direction of said forces, or in a direction that is related to the forces:

\[\mathbf{r}_{a}(k+1) = \mathbf{r}_{a}(k) - \mathbf{d}_{a}( \mathbf{F}_{a}(k) )\]

Where \(\mathbf{r}_{a}\) is the position of atom \(a\) and \(\mathbf{F}_{a}\) is the force over that atom. \(\mathbf{d}\) indicates a step in a direction which depends on the force. Finally, \(k\) indicates the current iteration step; essentially, this minimization is performed iteratively until the forces fall below a certain threshold (i.e., until the forces are so small that the structure is at a minimum).

Since SIESTA already provides the forces over the atoms (calculated via the Hellman-Feynman theorem), it is just a matter of choosing an appropriate algorithm and any additional parameters that algorithm might need. Note that the algorithms presented here will not avoid local minima; for that, you would require methods that are based on molecular dynamics (MD).

Note

Since we will be using the forces over the atoms to get an optimal structure, this means that all parameters related to the quality of your calculation need to be set properly beforehand. This includes the basis set, the mesh cut-off, the k-point sampling and the level of theory.

For the purposes of this tutorial we will be using very low quality settings to speed up the calculations; but you should not use this in any serious calculation.

Variable Cell Optimizations¶

Since we use SIESTA to calculate structures for periodic systems, one might also need to optimize not only the atomic positions, but also the lattice vectors themselves. This is a bit more complex than the simple geometry case, as it also requires the optimization of the stress tensor. This tensor is defined as:

\[\sigma_{ij} = \frac{1}{V} \frac{\partial E}{\partial e_{ij}}\]

Where \(e_{ij}\), the strain tensor, indicates the change in the lattice vectors. If we take the lattice vectors as the columns of a matrix \(h\), then the strain tensor appears as:

\[h' = h + e.h = (I + e).h\]

Where \(I\) is the identity matrix. Now, in order to couple the movement of the atoms with the changes of the unit cell, we re-define the atomic coordinates as re-scaled by the lattice vectors:

\[ \begin{align}\begin{aligned}\mathbf{s}_{a} = \sum_{j} r_{a} \mathbf{a^{*}}_{j} | \mathbf{a}_{j} |\\\mathbf{t}_{a} = \sum_{j} f_{a} \frac{ \mathbf{a}_{j} }{ | \mathbf{a}_{j} | }\end{aligned}\end{align} \]

Where \(\mathbf{s}_{a}\) and \(\mathbf{t}_{a}\) are the scaled atomic coordinates and forces, respectively. The vectors \(\mathbf{a^{*}}\) are defined as the reciprocal lattice vectors without a \(2\pi\) factor, so that:

\[\mathbf{a}_{i} \cdot \mathbf{a^{*}}_{j} = \delta_{ij}\]

The last thing we need to define is the variables that will represent our cell. For the cell itself, we will use as “coordinates” the strain tensor; while as gradient we will use the stress tensor mutiplied by the volume:

\[ \begin{align}\begin{aligned}\mathbf{s}_{cell} = cell - cell_{initial}\\\mathbf{t}_{cell} = - \sigma . V\end{aligned}\end{align} \]

We then use all the \(\mathbf{s(k)}\) combined as the variables for our minimization algorithm, and \(\mathbf{t(k)}\) as the gradients. Once we have our new \(\mathbf{s(k+1)}\) coordinates, we recover the cell vectors, and re-scale again the atomic coordinates to get \(\mathbf{r}_{a}(k+1)\).

If one wants to minimize the lattice under a given external pressure or stress, the equations above are modified so that the gradient now includes the target stress:

\[\mathbf{t}_{cell} = - (\sigma - \sigma_{target}) . V\]

Since the new coordinates depend on the lattice vectors themselves, variable cell optimizations usually take a lot more to converge than regular geometry optimizations.

basic Conjugate gradients with SiH¶

We will use an example of a cubic box of 64 Silicon atoms, with an extra H atom placed in the middle of a Si-Si bond. The current structure is not quite optimal, and the H atom will try to ‘push out’ its neighbors to lower the total energy and minimize the forces.

As usual, we have our fdf input file and our pseudopotentials:

Pseudos: Si.psml, H.psml
Input file: sih.fdf

A peek at the fdf file

SystemName          H in 64-atom silicon
SystemLabel         sih

NumberOfAtoms       65
NumberOfSpecies     2
%block ChemicalSpeciesLabel
 1  14  Si
 2   1  H
%endblock ChemicalSpeciesLabel

# Note that this basis is truly small and only set as an example.
PAO.BasisSize       SZ
PAO.EnergyShift     300 meV

# Warning: if you change the basis set size, you also need to change the
# number of eigenstates!
NumberOfEigenstates  150
MaxSCFIterations     50
DM.MixingWeight      0.3
DM.NumberPulay       4
DM.Tolerance         1.d-3
DM.UseSaveDM         T

MeshCutoff             40.0 Ry
SolutionMethod         diagon
ElectronicTemperature  25 meV

MD.TypeOfRun     CG
MD.Steps         50
MD.MaxForceTol   0.05 eV/Ang

WriteCoorXmol T

# System coordinates and cell vectors.
LatticeConstant     5.430 Ang
%block LatticeVectors
2.0 0.0 0.0
0.0 2.0 0.0
0.0 0.0 2.0
%endblock LatticeVectors

AtomicCoordinatesFormat  ScaledCartesian
%block AtomicCoordinatesAndAtomicSpecies
    0.000    0.000    0.000   1 #  Si  1
    0.250    0.250    0.250   1 #  Si  2
    0.000    0.500    0.500   1 #  Si  3
    0.250    0.750    0.750   1 #  Si  4
    0.500    0.000    0.500   1 #  Si  5
    0.750    0.250    0.750   1 #  Si  6
    0.500    0.500    0.000   1 #  Si  7
    0.750    0.750    0.250   1 #  Si  8
    1.000    0.000    0.000   1 #  Si  9
    1.250    0.250    0.250   1 #  Si 10
    1.000    0.500    0.500   1 #  Si 11
    1.250    0.750    0.750   1 #  Si 12
    1.500    0.000    0.500   1 #  Si 13
    1.750    0.250    0.750   1 #  Si 14
    1.500    0.500    0.000   1 #  Si 15
    1.750    0.750    0.250   1 #  Si 16
    0.000    1.000    0.000   1 #  Si 17
    0.250    1.250    0.250   1 #  Si 18
    0.000    1.500    0.500   1 #  Si 19
    0.250    1.750    0.750   1 #  Si 20
    0.500    1.000    0.500   1 #  Si 21
    0.750    1.250    0.750   1 #  Si 22
    0.500    1.500    0.000   1 #  Si 23
    0.750    1.750    0.250   1 #  Si 24
    0.000    0.000    1.000   1 #  Si 25
    0.250    0.250    1.250   1 #  Si 26
    0.000    0.500    1.500   1 #  Si 27
    0.250    0.750    1.750   1 #  Si 28
    0.500    0.000    1.500   1 #  Si 29
    0.750    0.250    1.750   1 #  Si 30
    0.500    0.500    1.000   1 #  Si 31
    0.750    0.750    1.250   1 #  Si 32
    1.000    1.000    0.000   1 #  Si 33
    1.250    1.250    0.250   1 #  Si 34
    1.000    1.500    0.500   1 #  Si 35
    1.250    1.750    0.750   1 #  Si 36
    1.500    1.000    0.500   1 #  Si 37
    1.750    1.250    0.750   1 #  Si 38
    1.500    1.500    0.000   1 #  Si 39
    1.750    1.750    0.250   1 #  Si 40
    1.000    0.000    1.000   1 #  Si 41
    1.250    0.250    1.250   1 #  Si 42
    1.000    0.500    1.500   1 #  Si 43
    1.250    0.750    1.750   1 #  Si 44
    1.500    0.000    1.500   1 #  Si 45
    1.750    0.250    1.750   1 #  Si 46
    1.500    0.500    1.000   1 #  Si 47
    1.750    0.750    1.250   1 #  Si 48
    0.000    1.000    1.000   1 #  Si 49
    0.250    1.250    1.250   1 #  Si 50
    0.000    1.500    1.500   1 #  Si 51
    0.250    1.750    1.750   1 #  Si 52
    0.500    1.000    1.500   1 #  Si 53
    0.750    1.250    1.750   1 #  Si 54
    0.500    1.500    1.000   1 #  Si 55
    0.750    1.750    1.250   1 #  Si 56
    1.000    1.000    1.000   1 #  Si 57
    1.250    1.250    1.250   1 #  Si 58
    1.000    1.500    1.500   1 #  Si 59
    1.250    1.750    1.750   1 #  Si 60
    1.500    1.000    1.500   1 #  Si 61
    1.750    1.250    1.750   1 #  Si 62
    1.500    1.500    1.000   1 #  Si 63
    1.750    1.750    1.250   1 #  Si 64
    1.125    1.125    1.125   2 #  H  65
%endblock AtomicCoordinatesAndAtomicSpecies

For the purposes of running this example quickly, we are using a SZ basis set and a very low cutoff. You might want to move beyond this later, but there should not be qualitative changes in the results.

First, let’s have a look at the options for geometry relaxation:

MD.TypeOfRun     CG
MD.Steps         50
MD.MaxForceTol   0.05 eV/Ang

CG stands for Conjugate Gradients, and it’s the default algorithm for geometry optimization within SIESTA. You can run the input file as usual:

$ siesta < sih.fdf > sih.out

The program will run for 13 geometry optimization steps; that is, steps in which the structure is modified in response to the forces until the cycle converges. The tolerance for convergence is given by the MD.MaxForceTol option showed above (note the explicit units), and it is slighty above the default in SIESTA (0.04 eV/Ang). MD.Steps indicates the maximum number of steps for convergence; SIESTA will stop and output a non-converged structure if this maximum is reached.

You might want to see the evolution of the total energy of the system during the run. You can scroll through the file or, as a shortcut, use the grep command (assuming you have redirected the output to file sih.out):

$ grep enth sih.out

There is a gain of more than 5 eV in the relaxation. If you want to plot these energies with gnuplot, you can use this plot script. You can use it by providing the output file as an argument:

$ ./plot_ene.sh sih.out

The image will be stored in a file called energy_iteration.png.

Now that we now that our structure has a much lower energy, thus being more stable, we can have a look at how things changed from the starting point. One way to look at this is via the .BONDS and .BONDS_FINAL files. Both contain the distances between each atom and its closest neighbours; the first file refers to the starting point and the second to the last point in the optimization (i.e. the converged structure). In particular, we will have a look at the neighbours of the lone H atom; since it is the last atom in our coordinates file, we can use the tail command to get the last few lines.

For the initial structure:

$ tail *.BONDS

And for the converged structure:

tail *.BONDS_FINAL

Compare the coordinates for the H atoms and the distances to its closest neighbours, the Silicon atoms 57 and 58. Without visualizing the structure, Can you deduce what’s happening to the structure during the optimization?

CG is not the only possible minimization algorithm within SIESTA. Other options for MD.TypeOfRun include, for example, the Fast Inertial Relaxation Engine (FIRE) and Broyden (MD.TypeOfRun Broyden, or MD.TypeOfRun FIRE). Edit your fdf file to use the Broyden method. Then try the FIRE method. How do they perform with respect to CG?

Performance for the FIRE method is a bit system-dependent; in this case, it is not really useful out of the box. If you want to play around with the method, you could try modifying the MD.Fire.TimeStep option.

Now that you have tried three different algorithms, pick the one which is fastest and lower the force tolerance to 0.01 eV/Ang (MD.MaxForceTol 0.01 eV/Ang). How many steps does it take to converge now? How much lower is the final energy? How different is the geometry?

basic Variable Cell Optimization using Si8¶

In the example above the atoms are moved in response to the forces, but the lattice vectors stay fixed. If you have scrolled through the output file you might have seen lines mentioning “stress”; for example:

Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      -78.32      -78.32      -78.32       -2.81       -2.81       -2.81

A non-zero stress means that the lattice vectors are not optimal, and the energy can be lowered by optimizing them. The relaxation algorithm will use now two sets of variables: atomic coordinates and lattice vectors. Both are moved in response to the forces and the stress, until these fall below a set tolerance.

The example we will use is an 8-atom Si cell (the ‘conventional cubic cell’). You can find the files here:

Pseudos: Si.psml
Input file: sih.fdf

Note the following options:

MD.TypeOfRun          CG
MD.Steps              100
MD.VariableCell       T
MD.MaxForceTol        0.1 eV/Ang
MD.MaxStressTol       0.1 GPa
MD.TargetPressure     0.0 GPa

First, we have enabled a variable cell optimization via the MD.VariableCell option. See that the tolerances for force and stress have different physical dimensions. The last line tells the program to aim for a zero (diagonal) stress (in practice, “atmospheric pressure”). See MD.MaxStressTol and MD.TargetPressure for more details.

Note

As with conventional geometry optimization, you can use the Broyden and FIRE methods with variable cell. In this simple example we will only use CG, but you can explore other methods if you are having trouble with your own systems.

Take note of the initial values for the simulation box:

LatticeConstant     5.535 Ang
%block LatticeVectors
  1.150  0.200  0.000
  0.000  1.050  0.000
 -0.100  0.000  0.900
%endblock LatticeVectors

We see that the first two cell vectors are too large, so the system is under tensile stress along the x and y directions. Conversely, the third cell vector is too small (compressive stress).

Run the example as usual:

$ siesta < si8.fdf > si8.out

You can now have a look at the stress components during the geometry optimization. You can do this by grepping for the word ‘oigt’ in the output file:

$ grep oigt si8.out

Output

Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       50.56       27.57      -41.29       51.64      -94.90      132.57
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       84.95       42.39      -23.99       21.24      -63.84      108.45
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      122.76       61.05       -4.00      -21.67      -23.96       71.85
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      117.94       41.47       20.47      -54.09       24.94       -9.88
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      127.13       62.67        5.25      -31.11       -9.55       51.67
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       99.12       32.94        0.69      -21.01       -4.87        9.26
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       56.55      -12.32       -8.11       -6.35        5.46      -48.15
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       76.56        9.16       -3.80      -13.46       -0.00      -21.27
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       36.04      -11.57      -30.23       -2.38       -6.85      -10.56
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -6.27      -38.36      -65.62       13.20      -22.08       16.62
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      -44.87      -50.19     -124.24       36.28      -44.69       51.10
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -8.00      -39.13      -67.10       13.98      -22.83       17.61
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      -10.19        1.18      -75.23       -8.52       10.95       64.68
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -5.99       -5.69      -77.98       -1.87        1.11       53.96
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):        3.70       10.90      -55.74       -2.97        3.45       34.84
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       10.93       25.08      -33.07       -3.49        5.96       16.90
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       17.22       37.56      -10.43       -4.34        8.20        0.34
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       22.45       48.59       11.86       -4.68        9.99      -14.66
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       17.67       38.46       -8.60       -4.47        8.41       -0.96
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -4.02       -0.14      -17.34       40.72       -9.49      -26.04
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):        8.95       23.53      -11.09       13.69        0.17      -11.65
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -4.54        6.90      -12.83        6.63        3.18       -2.19
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):      -30.80      -22.07      -16.88       -6.77        9.17       15.28
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -7.83        3.05      -13.36        4.93        3.77        0.01
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -7.97      -13.62       -7.67      -15.19       -1.74      -16.98
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -7.49       -3.28      -11.08       -2.90        1.83       -6.58
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -0.98        0.77       -2.00       -0.91       -2.63       -0.64
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):        9.68        7.26       12.89        2.34       -9.45        8.91
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -0.28        1.23       -1.07       -0.80       -3.10       -0.01
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):        5.41        2.95        8.25        0.33        6.47       -1.42
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):        1.09        1.62        1.09       -0.68       -0.94       -0.39
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -1.49       -0.88       -1.06       -1.75        1.17       -0.03
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -0.13        0.48        0.07       -1.13        0.02       -0.24
Stress tensor Voigt[x,y,z,yz,xz,xy] (kbar):       -0.59       -0.38       -0.65        0.08       -0.47       -0.28

There are off-diagonal components of the strain, leading to ‘shear’components of the stress. Note the signs and sizes of these initial stress components, and watch how they move towards zero. You can see the stress If you want to plot these values with gnuplot, you can use this plot script. It works as the energy plot script, and will store the results in a file called stress_iteration.png.

In addition, you can check the evolution of the energy as before:

$ grep enth si8.out

How do the stress and the energy evolve during the relaxation?

Have a look at the final shape for the cell vectors and the atomic positions. Look for “outcell” and “outcoor” in the last few lines of the output file, before the final energy decomposition. Have a look at the cell vector modules and angles. Is the shape what you expected? What happened to the Si atom coordinates?

basic Using Constraints¶

We will now have a look at how to use constraints in SIESTA. Constraints are used to fix the position of some atoms, or to keep the cell vectors fixed. You can find a full description for Geometry.Constraints in the manual.

You can find the input files here:

Pseudos: Si.psml
Input file: sih.fdf

The options are the same as before, with a new addition at the bottom:

%block Geometry.Constraints
  position from 1 to 4
  cellvector B
%endblock Geometry.Constraints

In this case, we are constraining the position of the first four atoms, and the shape of the second cell vector (hence the “B”). Vector constraints can only be applied if the cell vector has a single cartesian component; thus, we cannot use them for the first or third cell vectors (“A” and “C”).

Run the example as usual:

$ siesta < si8.fdf > si8.out

Again, have a look at the final shape for the cell vectors and the atomic positions (look for “outcell” and “outcoor”). Have a look at the cell vectors now. What happened to them? What happened to the Si atom coordinates now?

An important factor for constrained simulations is that the forces over constrained atoms (and the stress over the constrained cell vectors) must be zero. For the stress, we can verify this in the last part of the output file:

siesta: Stress tensor (static) (eV/Ang**3):
siesta:     0.000308    0.067632   -0.000418
siesta:     0.067631    0.030605   -0.005606
siesta:    -0.000415   -0.005607   -0.000535

siesta: Constrained stress tensor (static) (eV/Ang**3):
siesta:     0.000308    0.000000   -0.000416
siesta:     0.000000    0.000000    0.000000
siesta:    -0.000416    0.000000   -0.000535

Note that the constrained stress tensor is much more symmetric, since, again, it ignores the components related to the second lattice vector. The constrained forces cannot be obtained from the output file, but rather you can have a look at the .FAC and .FA files. The si8.FA file contains the forces on the atoms before applying constraints (in eV/Ang), and looks like this:

8
-0.667091721E+00   0.870915839E+00   0.132590789E+01
 0.658235184E+00  -0.866434730E+00  -0.132492768E+01
-0.667091725E+00   0.870915847E+00   0.132590789E+01
 0.658235190E+00  -0.866434723E+00  -0.132492767E+01
-0.461767750E-02  -0.555508348E-02  -0.310101580E-02
 0.599223267E-02   0.367068527E-02   0.218699862E-02
-0.461767155E-02  -0.555506740E-02  -0.310102455E-02
 0.599224312E-02   0.367068567E-02   0.218699625E-02

Note how the forces over the first four atoms (the ones we are constrained) are actually huge when compared to the rest. Now, if we have a look at the si8.FAC file, we can see that those forces are set to zero after applying the constraints:

8
 0.000000000E+00   0.000000000E+00   0.000000000E+00
 0.000000000E+00   0.000000000E+00   0.000000000E+00
 0.000000000E+00   0.000000000E+00   0.000000000E+00
 0.000000000E+00   0.000000000E+00   0.000000000E+00
-0.461767750E-02  -0.555508348E-02  -0.310101580E-02
 0.599223267E-02   0.367068527E-02   0.218699862E-02
-0.461767155E-02  -0.555506740E-02  -0.310102455E-02
 0.599224312E-02   0.367068567E-02   0.218699625E-02

You can also get the maximum value for the forces over the constrained atoms by grepping the word ” constrained” in the output file:

$ grep " constrained" si8.out

intermediate Using Z-matrix constraints in a water molecule¶

Sometimes the important structural degrees of freedom that we want to optimize are not easily represented in cartesian coordinates. In molecules, for example, bond lengths, angles, and torsion angles are much more relevant than particular values of the cartesian coordinates. For many years, chemists have used more appropriate ways to represent molecular structure, and the Zmatrix is one of them.

You can find the input files here:

Pseudos: O.psml, H.psml
Input files: h2oZ.fdf

Consider this section in file h2oZ.fdf:

ZM.UnitsLength Ang
ZM.UnitsAngle deg
%block Zmatrix
molecule_cartesian
  1   0 0 0   0.0   0.0  0.0        0 0 0
  2   1 0 0   1.0  90.0 37.743919   0 0 0
  2   1 2 0   1.0 106.0 90.0        0 0 0
%endblock Zmatrix

The first two options, ZM.UnitsLength and ZM.UnitsAngle, indicate the units for the bond distances and angles. The next block, Zmatrix, contains the Z-matrix for the system.

In the first section, we have first four integer numbers in which we indicate the atomic species index and the reference for the bond length, the angle, and the dihedral or torsional angle. Since the first atom indicates the origin, it is not referenced against any other atom. Thus the values of zero for everything. Things become clearer once we have a look at the third atom: it is a hydrogen atom (hence the species index “2”), which has a bond length of 1.0 with respect to atom 1, has an angle of 106° with respect to atoms 1 and 2, and has a dihedral of 90° with respect to atoms 1, 2 and the origin. Ignore the last three zeros for now.

If we run the example as it is, you will notice that the system will not relax, and will run only a single geometry optimization step. This is because the entire block is considered as a set of constants, therefore, there is nothing to optimize since bonds and angles are fixed. In order to allow them to change we need to introduce a new section in the ZMatrix block called “variables”:

%block Zmatrix
  molecule_cartesian
    1   0 0 0   0.0  0.0  0.0       0 0 0
    2   1 0 0   1.0 90.0 37.743919  0 0 0
    2   1 2 0   1.0  HOH 90.0       0 1 0
  variables
    HOH 106.0
%endblock Zmatrix

The variable “HOH” is introduced in the block, and it is used in the angle for the second H atom. Note the new value “1” in the last line of the block, which indicates that the angle is now a variable that can be changed.

Modify the input file as indicated above, and run the example:

$ siesta < h2oZ.fdf > h2oZ.out

We can check the final value for the Z-matrix variables in the output file called ZMATRIX_VARS:

HOH            103.719265373056

We can see that the angle has now changed. By doing this same thing for the bonds, we can now enable a full optimization of the molecule:

%block Zmatrix
molecule_cartesian
  1 0 0 0   0.0  0.0  0.0      0 0 0
  2 1 0 0   HO1 90.0 37.743919 1 0 0
  2 1 2 0   HO2  HOH 90.0      1 1 0
variables
    HO1 1.0
    HO2 1.0
    HOH 106.0
%endblock Zmatrix

Modify the fdf file, or use this h2oZ_full.fdf file provided here. Then run the example as before:

$ siesta < h2oZ_full.fdf > h2oZ_full.out

Which are the final values for the bond lengths and the bond angle?

Now that you are familiar with the format, how would you change the file to keep the angle constant? Do the final bond lengths change if this angle is kept constant at 106°?

To finish our look at the Z-Matrix optimization, lets have a look at the last few variables:

ZM.ForceTolLength  0.04 eV/Ang
ZM.ForceTolAngle   0.0001 eV/deg
ZM.MaxDisplLength  0.1 Ang
ZM.MaxDisplAngle   20.0 deg

The first two variables, ZM.ForceTolLength and ZM.ForceTolAngle, are used to set the tolerance for the forces over the bond lengths and the angles, in the same vein as the tolerance used for CG optimization with cartesian coordinates.

The last two variables, ZM.MaxDisplLength and ZM.MaxDisplAngle, indicate the maximum displacement allowed for bonds and angles when performing the optimization. For example, if at any step of the optimization, any of the bonds would change by more than 0.1 Ang, the change is kept at 0.1 Ang instead.