Structural optimization using forces and stress¶
Are you using this tutorial for a school?
Ask tutors if you need to set up your local environment. Do that now before proceeding.
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¶
Additional information: Slides
These slides from previous presentations might help improve your understanding of the concepts presented in this tutorial:
Presentation by Emilio Artacho
Presentation by Marivi Fernadez-Serra
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:
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:
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:
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:
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:
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:
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:
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¶

Are you under a school environment?
Enter directory 01-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:
A peek at the fdf file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | 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:
28 29 30 | 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
Output
Target enthalpy (eV/cell) -7361.1235
Target enthalpy (eV/cell) -7363.6491
Target enthalpy (eV/cell) -7366.1559
Target enthalpy (eV/cell) -7367.6264
Target enthalpy (eV/cell) -7367.5108
Target enthalpy (eV/cell) -7367.7042
Target enthalpy (eV/cell) -7367.3450
Target enthalpy (eV/cell) -7368.0054
Target enthalpy (eV/cell) -7367.9459
Target enthalpy (eV/cell) -7368.0551
Target enthalpy (eV/cell) -7368.0636
Target enthalpy (eV/cell) -7368.0763
Target enthalpy (eV/cell) -7368.0768
Target enthalpy (eV/cell) -7368.0804
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
Output
63 Si 2.3513 Ang. Really at: 15.3918 15.3918 10.2612
Neighbors of: 65 H at: 11.5439 11.5439 11.5439
58 Si 1.1756 Ang. Really at: 12.8265 12.8265 12.8265
57 Si 1.1756 Ang. Really at: 10.2612 10.2612 10.2612
32 Si 2.9586 Ang. Really at: 7.6959 7.6959 12.8265
61 Si 2.9586 Ang. Really at: 15.3918 10.2612 15.3918
63 Si 2.9586 Ang. Really at: 15.3918 15.3918 10.2612
22 Si 2.9586 Ang. Really at: 7.6959 12.8265 7.6959
59 Si 2.9586 Ang. Really at: 10.2612 15.3918 15.3918
12 Si 2.9586 Ang. Really at: 12.8265 7.6959 7.6959
And for the converged structure:
tail *.BONDS_FINAL
Output
45 Si 2.3490 Ang. Really at: 15.3905 20.5085 15.3905
Neighbors of: 65 H at: 11.5439 11.5439 11.5439
58 Si 1.6624 Ang. Really at: 13.3576 13.3576 13.3576
57 Si 1.6624 Ang. Really at: 9.7301 9.7301 9.7301
63 Si 3.0411 Ang. Really at: 15.4946 15.4946 10.1984
32 Si 3.0411 Ang. Really at: 7.5932 7.5932 12.8893
59 Si 3.0411 Ang. Really at: 10.1984 15.4946 15.4946
12 Si 3.0411 Ang. Really at: 12.8893 7.5932 7.5932
22 Si 3.0411 Ang. Really at: 7.5932 12.8893 7.5932
61 Si 3.0411 Ang. Really at: 15.4946 10.1984 15.4946
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?
Answer: What’s happening to the structure?
The coordinates for hydrogen do not change, but both silicon atoms increase their distance. If you watch closely, you will see that both silicons are in exact opposite directions; thus, they are both getting “pushed” symmetrically by the H atom. Since the H atom has equal forces in both directions, it does not move.
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?
Answer: How do Broyden and FIRES perform?
Broyden usually takes fewer steps to converge; in this case, 6 steps. The final
energy is -7368.0762
, which is only a difference of 4 meV with respect to
the one we saw for CG. This slight difference comes from the fact that the
tolerance we are using is not stringent enough.

If we have a look at the sih.BONDS_FINAL file, you can see that the difference in geometry is minimal:
Neighbors of: 65 H at: 11.5439 11.5439 11.5439
58 Si 1.6610 Ang. Really at: 13.3560 13.3560 13.3560
57 Si 1.6610 Ang. Really at: 9.7317 9.7317 9.7317
63 Si 3.0366 Ang. Really at: 15.4874 15.4874 10.1929
32 Si 3.0366 Ang. Really at: 7.6003 7.6003 12.8948
61 Si 3.0366 Ang. Really at: 15.4874 10.1929 15.4874
22 Si 3.0366 Ang. Really at: 7.6003 12.8948 7.6003
59 Si 3.0366 Ang. Really at: 10.1929 15.4874 15.4874
12 Si 3.0366 Ang. Really at: 12.8948 7.6003 7.6003
FIRE in this case takes longer than CG, about 23 steps.

However, FIRE really shines when structures much harder to converge, so it is
a bit system dependent. As with Broyden, the differences in final energy
(-7368.0790
) and coordinates are minimal:
Neighbors of: 65 H at: 11.5439 11.5439 11.5439
57 Si 1.6637 Ang. Really at: 9.7287 9.7287 9.7287
58 Si 1.6637 Ang. Really at: 13.3590 13.3590 13.3590
12 Si 3.0384 Ang. Really at: 12.8908 7.5972 7.5972
59 Si 3.0384 Ang. Really at: 10.1969 15.4906 15.4906
22 Si 3.0384 Ang. Really at: 7.5972 12.8908 7.5972
61 Si 3.0384 Ang. Really at: 15.4906 10.1969 15.4906
32 Si 3.0384 Ang. Really at: 7.5972 7.5972 12.8908
63 Si 3.0384 Ang. Really at: 15.4906 15.4906 10.1969
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?
Answer: Evaluating tighter convergence.
According to the previous results, we have picked the Broyden algorithm, as
it provides faster convergence than both CG and FIRE. With the tighter, we now
take 10 minimization steps instead of 6. The final energy is now -7368.0843
,
which is 8.1 meV lower than the previous result. The difference in geometry is
also minimal, with differences well under 0.01 Ang:
Neighbors of: 65 H at: 11.5439 11.5439 11.5439
58 Si 1.6636 Ang. Really at: 13.3589 13.3589 13.3589
57 Si 1.6636 Ang. Really at: 9.7288 9.7288 9.7288
61 Si 3.0464 Ang. Really at: 15.5030 10.2054 15.5030
59 Si 3.0464 Ang. Really at: 10.2054 15.5030 15.5030
22 Si 3.0464 Ang. Really at: 7.5848 12.8823 7.5848
63 Si 3.0464 Ang. Really at: 15.5030 15.5030 10.2054
12 Si 3.0464 Ang. Really at: 12.8823 7.5848 7.5848
32 Si 3.0464 Ang. Really at: 7.5848 7.5848 12.8823
Additional information: Visualizing your geometry optimization.
Enabling the WriteCoorXmol option (WriteCoorXmol T
)
will produce an XYZ-formatted file with the *.ANI extension.
For geometry optimizations, the file *.ANI contains all the structures generated during the relaxation. It can be processed by various graphical tools, such as XCrysden , VMD or Ovito .
Additional information: Using LUA to get more minimization engines.
These optimization engines are implemented within SIESTA itself, and are generally sufficient. Alternatively, it is possible to use the Lua scripting engine embedded in SIESTA to access other algorithms, for example L-BFGS.
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.
Are you under a school environment?
Enter directory 02-varcell_cg
The example we will use is an 8-atom Si cell (the ‘conventional cubic cell’). You can find the files here:
A peek at the fdf file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 | SystemName Si - 8 atoms supercell SystemLabel si8 NumberOfAtoms 8 NumberOfSpecies 1 %block ChemicalSpeciesLabel 1 14 Si %endblock ChemicalSpeciesLabel PAO.BasisSize SZ PAO.EnergyShift 300 meV 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 MeshCutoff 30.0 Ry %block kgrid_Monkhorst_Pack 2 0 0 0.0 0 2 0 0.0 0 0 2 0.0 %endblock kgrid_Monkhorst_Pack MaxSCFIterations 50 DM.MixingWeight 0.3 DM.NumberPulay 3 DM.Tolerance 1.d-3 DM.UseSaveDM SolutionMethod diagon ElectronicTemperature 25 meV 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 WriteMDHistory T AtomicCoordinatesFormat ScaledbyLatticeVectors %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 %endblock AtomicCoordinatesAndAtomicSpecies |
Note the following options:
42 43 44 45 46 47 | 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:
16 17 18 19 20 21 | 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
Output

Target enthalpy (eV/cell) -916.1832
Target enthalpy (eV/cell) -917.0587
Target enthalpy (eV/cell) -917.7618
Target enthalpy (eV/cell) -917.3919
Target enthalpy (eV/cell) -917.8380
Target enthalpy (eV/cell) -918.3276
Target enthalpy (eV/cell) -918.3588
Target enthalpy (eV/cell) -918.4119
Target enthalpy (eV/cell) -918.7417
Target enthalpy (eV/cell) -918.8595
Target enthalpy (eV/cell) -918.7842
Target enthalpy (eV/cell) -918.8600
Target enthalpy (eV/cell) -919.0721
Target enthalpy (eV/cell) -919.0755
Target enthalpy (eV/cell) -919.3231
Target enthalpy (eV/cell) -919.4762
Target enthalpy (eV/cell) -919.5342
Target enthalpy (eV/cell) -919.5006
Target enthalpy (eV/cell) -919.5347
Target enthalpy (eV/cell) -919.4793
Target enthalpy (eV/cell) -919.5881
Target enthalpy (eV/cell) -919.6368
Target enthalpy (eV/cell) -919.5803
Target enthalpy (eV/cell) -919.6385
Target enthalpy (eV/cell) -919.6228
Target enthalpy (eV/cell) -919.6494
Target enthalpy (eV/cell) -919.6572
Target enthalpy (eV/cell) -919.6394
Target enthalpy (eV/cell) -919.6572
Target enthalpy (eV/cell) -919.6491
Target enthalpy (eV/cell) -919.6580
Target enthalpy (eV/cell) -919.6580
Target enthalpy (eV/cell) -919.6582
Target enthalpy (eV/cell) -919.6583
How do the stress and the energy evolve during the relaxation?
Answer: How do the stress and the energy evolve during the relaxation?
While both the stress and the energy decrease across the simulation, they do not do so monotonically. This is the typical case for a variable cell run in general. In most cases, the energy can decrease monotonically if the quality of the simulation is good in terms of basis set and mesh cutoff, but it is not guaranteed.
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?
Answer: Final cell vectors and atomic positions.
Let’s first have a look at the atomic positions:
outcoor: Relaxed atomic coordinates (fractional):
-0.00023685 0.00006112 -0.00083770 1 1 Si
0.24986519 0.24997244 0.24915731 1 2 Si
-0.00023685 0.50006112 0.49916230 1 3 Si
0.24986519 0.74997244 0.74915731 1 4 Si
0.49976315 0.00006112 0.49916230 1 5 Si
0.74986519 0.24997244 0.74915731 1 6 Si
0.49976315 0.50006112 -0.00083770 1 7 Si
0.74986519 0.74997244 0.24915731 1 8 Si
You can see that they did not change much from the starting position. This is somewhat expected, since the initial positions correspond to a high-symmetry configuration; in particular, these are the standard positions for Si in the diamond structure.
Regarding the final lattice vectors:
outcell: Unit cell vectors (Ang):
5.531028 0.503687 0.269460
-0.503545 5.539740 -0.039584
-0.272517 0.014028 5.554600
outcell: Cell vector modules (Ang) : 5.560448 5.562719 5.561298
outcell: Cell angles (23,13,12) (deg): 90.0091 90.0065 90.0102
outcell: Cell volume (Ang**3) : 172.0177
You can see here that even if the cartesian components of the cell vectors look strange, they actually correspond to a cubic cell. Notice the vector modules and the angles between them. Essentially, the cell has rotated during the optimization, but it is still a cubic cell.
Additional information: Variable cell and vacuum regions.
When doing a variable cell optimization in 2D/slab systems, take into account that the size of the vacuum might change during the optimization. You might want to freeze or constraint the lattice vector in the non-periodic direction; you can use the constraints explained below for this purpose.
Additional information: Using Quenched MD to relax.
Variable cell optimization can be quite slow, especially for large or complex systems. In some cases, an alternative approach is to use a Quenched Molecular Dynamics. You can see more about this under the molecular dynamics tutorial.
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.
Are you under a school environment?
Enter directory 03-constraints
You can find the input files here:
A peek at the fdf file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 | SystemName Si - 8 atoms supercell SystemLabel si8 NumberOfAtoms 8 NumberOfSpecies 1 %block ChemicalSpeciesLabel 1 14 Si %endblock ChemicalSpeciesLabel PAO.BasisSize SZ PAO.EnergyShift 300 meV 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 MeshCutoff 30.0 Ry %block kgrid_Monkhorst_Pack 2 0 0 0.0 0 2 0 0.0 0 0 2 0.0 %endblock kgrid_Monkhorst_Pack MaxSCFIterations 50 DM.MixingWeight 0.3 DM.NumberPulay 3 DM.Tolerance 1.d-3 DM.UseSaveDM SolutionMethod diagon ElectronicTemperature 25 meV 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 WriteMDHistory T AtomicCoordinatesFormat ScaledbyLatticeVectors %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 %endblock AtomicCoordinatesAndAtomicSpecies %block Geometry.Constraints position from 1 to 4 cellvector B %endblock Geometry.Constraints |
The options are the same as before, with a new addition at the bottom:
63 64 65 66 | %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?
Answer: Final cell vectors and atomic positions.
Let’s first have a look at the atomic positions:
outcoor: Relaxed atomic coordinates (fractional):
0.00000000 0.00000000 0.00000000 1 1 Si
0.25000000 0.25000000 0.25000000 1 2 Si
-0.00000000 0.50000000 0.50000000 1 3 Si
0.25000000 0.75000000 0.75000000 1 4 Si
0.50562329 -0.01537356 0.53639122 1 5 Si
0.74441312 0.26554844 0.71384696 1 6 Si
0.50562329 0.48462644 0.03639122 1 7 Si
0.74441312 0.76554844 0.21384696 1 8 Si
You can see now that the coordinates for the first four atoms have not changed at all. We can also say the same thing for the second lattice vector:
outcell: Unit cell vectors (Ang):
5.692548 1.107000 0.139511
0.000000 5.811750 0.000000
-0.385822 0.000000 5.202928
outcell: Cell vector modules (Ang) : 5.800863 5.811750 5.217214
outcell: Cell angles (23,13,12) (deg): 90.0000 92.7849 78.9986
outcell: Cell volume (Ang**3) : 172.4448
Note that the cell is not cubic anymore!
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
1 -0.667091721E+00 0.870915839E+00 0.132590789E+01
2 0.658235184E+00 -0.866434730E+00 -0.132492768E+01
3 -0.667091725E+00 0.870915847E+00 0.132590789E+01
4 0.658235190E+00 -0.866434723E+00 -0.132492767E+01
5 -0.461767750E-02 -0.555508348E-02 -0.310101580E-02
6 0.599223267E-02 0.367068527E-02 0.218699862E-02
7 -0.461767155E-02 -0.555506740E-02 -0.310102455E-02
8 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
1 0.000000000E+00 0.000000000E+00 0.000000000E+00
2 0.000000000E+00 0.000000000E+00 0.000000000E+00
3 0.000000000E+00 0.000000000E+00 0.000000000E+00
4 0.000000000E+00 0.000000000E+00 0.000000000E+00
5 -0.461767750E-02 -0.555508348E-02 -0.310101580E-02
6 0.599223267E-02 0.367068527E-02 0.218699862E-02
7 -0.461767155E-02 -0.555506740E-02 -0.310102455E-02
8 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
Output
Max 2.146535 constrained
Max 1.568286 constrained
Max 0.828636 constrained
Max 0.917629 constrained
Max 0.836277 constrained
Max 0.394626 constrained
Max 0.954425 constrained
Max 0.749100 constrained
Max 0.190859 constrained
Max 0.468002 constrained
Max 0.235445 constrained
Max 0.425601 constrained
Max 0.191185 constrained
Max 0.171746 constrained
Max 0.316498 constrained
Max 0.162930 constrained
Max 0.467104 constrained
Max 0.153301 constrained
Max 0.082862 constrained
Max 0.081149 constrained
Max 0.214416 constrained
Max 0.005992 constrained
You can also use this script
to plot the maximum constrained force per iteration. It works like the other
scripts before:
$ ./plot_constr.sh si8.out
And you will get a file called max_force_iteration.png:

intermediate Using Z-matrix constraints in a water molecule¶

Are you under a school environment?
Enter the directory 04-h2o_zmatrix
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:
A peek at the fdf file
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 | # # Z-matrix specification for a H20 molecule # SystemName h2oZ SystemLabel h2oZ NumberOfSpecies 2 NumberOfAtoms 3 %block ChemicalSpeciesLabel 1 8 O 2 1 H %endblock ChemicalSpeciesLabel AtomicCoordinatesFormat NotScaledCartesianAng 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 MD.TypeOfRun CG MD.Steps 20 ZM.ForceTolLength 0.04 eV/Ang ZM.ForceTolAngle 0.0001 eV/deg ZM.MaxDisplLength 0.1 Ang ZM.MaxDisplAngle 20.0 deg MeshCutoff 120 Ry PAO.EnergyShift 0.01 Ry DM.NumberPulay 5 DM.MixingWeight 0.3 |
Consider this section in file h2oZ.fdf:
16 17 18 19 20 21 22 23 | 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?
Answer: Final bond lengths and angle.
We can get this information, again from the ZMATRIX_VARS file:
HO1 0.976622194379286
HO2 0.976646141366239
HOH 104.955795530240
As you can see the bonds have shortened a bit, and the angle has decreased to 104.96°. Does this make sense? Consider that the experimental values are 0.96 Ang for the bond lengths and 104.5° for 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°?
Answer: Optimized bond length with fixed angle.
The new ZMatrix block should look like this:
%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 106.0 90.0 1 0 0
variables
HO1 1.0
HO2 1.0
%endblock Zmatrix
The final bond lengths, again from the ZMATRIX_VARS file, are:
HO1 0.976415483377268
HO2 0.976194280048438
Note that the bonds are very similar to the ones obtained before.
To finish our look at the Z-Matrix optimization, lets have a look at the last few variables:
27 28 29 30 | 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.