:sequential_nav: next .. _tutorial-basic-structure-optimization: Structural optimization using forces and stress =============================================== .. dropdown:: **Do you want to get all tutorial materials on your local machine?** :color: primary You can get all files and set up your local environment. :ref:`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. .. button-ref:: tutorial-basic-structure-optimization-start :color: success :outline: :expand: :shadow: Take me to the exercises! :bdg-info:`background` Background Information --------------------------------------------- Geometry optimization ^^^^^^^^^^^^^^^^^^^^^ .. dropdown:: Additional information: Slides :color: info 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: .. math:: \mathbf{r}_{a}(k+1) = \mathbf{r}_{a}(k) - \mathbf{d}_{a}( \mathbf{F}_{a}(k) ) Where :math:`\mathbf{r}_{a}` is the position of atom :math:`a` and :math:`\mathbf{F}_{a}` is the force over that atom. :math:`\mathbf{d}` indicates a step in a direction which depends on the force. Finally, :math:`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: .. math:: \sigma_{ij} = \frac{1}{V} \frac{\partial E}{\partial e_{ij}} Where :math:`e_{ij}`, the strain tensor, indicates the change in the lattice vectors. If we take the lattice vectors as the columns of a matrix :math:`h`, then the strain tensor appears as: .. math:: h' = h + e.h = (I + e).h Where :math:`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: .. math:: \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} | } Where :math:`\mathbf{s}_{a}` and :math:`\mathbf{t}_{a}` are the scaled atomic coordinates and forces, respectively. The vectors :math:`\mathbf{a^{*}}` are defined as the reciprocal lattice vectors without a :math:`2\pi` factor, so that: .. math:: \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: .. math:: \mathbf{s}_{cell} = cell - cell_{initial} \mathbf{t}_{cell} = - \sigma . V We then use all the :math:`\mathbf{s(k)}` combined as the variables for our minimization algorithm, and :math:`\mathbf{t(k)}` as the gradients. Once we have our new :math:`\mathbf{s(k+1)}` coordinates, we recover the cell vectors, and re-scale again the atomic coordinates to get :math:`\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: .. math:: \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. .. _tutorial-basic-structure-optimization-start: :bdg-success:`basic` Conjugate gradients with SiH ------------------------------------------------- .. image:: images/sys-sih.png :width: 600px :align: center .. dropdown:: **Are you under a local environment?** :color: primary 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: * Pseudos: :download:`Si.psml`, :download:`H.psml` * Input file: :download:`sih.fdf` .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/01-SiH/sih.fdf :linenos: 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: .. literalinclude:: work-files/01-SiH/sih.fdf :linenos: :lineno-match: :start-at: MD.TypeOfRun :end-at: MD.MaxForceTol CG stands for `Conjugate Gradients `_, and it's the default algorithm for geometry optimization within SIESTA. You can run the input file as usual: .. code-block:: console $ 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 :ref:`MD.MaxForceTol` option showed above (note the explicit units), and it is slighty above the default in SIESTA (0.04 eV/Ang). :ref:`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*): .. code-block:: console $ grep enth sih.out .. dropdown:: Output .. code-block:: 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 :download:`plot script`. You can use it by providing the output file as an argument: .. code-block:: console $ ./plot_ene.sh sih.out .. dropdown:: Output .. image:: images/cg1.png :width: 600px :align: center 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: .. code-block:: console $ tail *.BONDS .. dropdown:: Output .. code-block:: 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 .. dropdown:: Output .. code-block:: 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?* .. dropdown:: Answer: What's happening to the structure? :color: success 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 :ref:`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?* .. dropdown:: Answer: How do Broyden and FIRES perform? :color: success 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. .. image:: images/broyden1.png :width: 600px :align: center 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. .. image:: images/fire1.png :width: 600px :align: center 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 :ref:`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?* .. dropdown:: Answer: Evaluating tighter convergence. :color: success 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 .. dropdown:: Additional information: Visualizing your geometry optimization. :color: info Enabling the :ref:`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 `_ . .. dropdown:: Additional information: Using LUA to get more minimization engines. :color: info These optimization engines are implemented within SIESTA itself, and are generally sufficient. Alternatively, it is possible to use the :ref:`Lua scripting engine` embedded in SIESTA to access other algorithms, for example L-BFGS. :bdg-success:`basic` Variable Cell Optimization using Si8 --------------------------------------------------------- .. image:: images/sys-si8.png :width: 400px :align: center 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. .. dropdown:: **Are you under a local environment?** :color: primary 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: * Pseudos: :download:`Si.psml` * Input file: :download:`sih.fdf` .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/02-varcell_cg/si8.fdf :linenos: Note the following options: .. literalinclude:: work-files/02-varcell_cg/si8.fdf :linenos: :lineno-match: :start-at: MD.TypeOfRun :end-at: MD.TargetPressure First, we have enabled a variable cell optimization via the :ref:`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 :ref:`MD.MaxStressTol` and :ref:`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: .. literalinclude:: work-files/02-varcell_cg/si8.fdf :linenos: :lineno-match: :start-at: LatticeConstant :end-at: %endblock 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: .. code-block:: console $ 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: .. code-block:: console $ grep oigt si8.out .. dropdown:: Output .. image:: images/varcell_stress.png :width: 600px :align: center .. code-block:: 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 :download:`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: .. code-block:: console $ grep enth si8.out .. dropdown:: Output .. image:: images/varcell_ene.png :width: 600px :align: center .. code-block:: 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?* .. dropdown:: Answer: How do the stress and the energy evolve during the relaxation? :color: success 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?* .. dropdown:: Answer: Final cell vectors and atomic positions. :color: success 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. .. dropdown:: Additional information: Variable cell and vacuum regions. :color: info 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. .. dropdown:: Additional information: Using Quenched MD to relax. :color: info 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 :ref:`molecular dynamics tutorial`. :bdg-success:`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 :ref:`Geometry.Constraints` in the manual. .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *03-constraints* You can find the input files here: * Pseudos: :download:`Si.psml` * Input file: :download:`sih.fdf` .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/03-constraints/si8.fdf :linenos: The options are the same as before, with a new addition at the bottom: .. literalinclude:: work-files/03-constraints/si8.fdf :linenos: :lineno-match: :start-at: %block Geometry.Constraints :end-at: %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: .. code-block:: console $ 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?* .. dropdown:: Answer: Final cell vectors and atomic positions. :color: success 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: .. code-block:: console $ grep " constrained" si8.out .. dropdown:: Output .. code-block:: 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 :download:`this script` to plot the maximum constrained force per iteration. It works like the other scripts before: .. code-block:: console $ ./plot_constr.sh si8.out And you will get a file called *max_force_iteration.png*: .. image:: images/max_force_iteration.png :width: 600px :align: center :bdg-warning:`intermediate` Using Z-matrix constraints in a water molecule -------------------------------------------------------------------------- .. image:: images/sys-wat.png :width: 400px :align: center .. dropdown:: **Are you under a local environment?** :color: primary 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: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`h2oZ.fdf` .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/04-h2o_zmatrix/h2oZ.fdf :linenos: Consider this section in file *h2oZ.fdf*: .. literalinclude:: work-files/04-h2o_zmatrix/h2oZ.fdf :linenos: :lineno-match: :start-at: ZM.UnitsLength :end-at: %endblock Zmatrix The first two options, :ref:`ZM.UnitsLength` and :ref:`ZM.UnitsAngle`, indicate the units for the bond distances and angles. The next block, :ref:`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: .. code-block:: console $ 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 :download:`h2oZ_full.fdf` file provided here. Then run the example as before: .. code-block:: console $ siesta < h2oZ_full.fdf > h2oZ_full.out *Which are the final values for the bond lengths and the bond angle?* .. dropdown:: Answer: Final bond lengths and angle. :color: success 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°?* .. dropdown:: Answer: Optimized bond length with fixed angle. :color: success 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: .. literalinclude:: work-files/04-h2o_zmatrix/h2oZ.fdf :linenos: :lineno-match: :start-at: ZM.ForceTolLength :end-at: ZM.MaxDisplAngle The first two variables, :ref:`ZM.ForceTolLength` and :ref:`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, :ref:`ZM.MaxDisplLength` and :ref:`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.