# Tutorial:Geometry optimization

In this tutorial we will show how to perform geometry optimizations using Octopus.

## Contents

## Methane molecule

### Manual optimization

We will start with the methane molecule. Before using any sophisticated algorithm for the geometry optimization, it is quite instructive to do it first "by hand". Of course, this is only practical for systems with very few degrees of freedom, but this is the case for methane, as we only need to optimize the CH bond-length. We will use the following files.

#### inp

`CalculationMode`

= gs`UnitsOutput`

= eV_Angstrom`Radius`

= 3.5*angstrom`Spacing`

= 0.18*angstrom CH = 1.2*angstrom %`Coordinates`

"C" | 0 | 0 | 0 "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) %

This input file is similar to ones used in other tutorials. We have seen in the Total energy convergence tutorial that for these values of the spacing and radius the total energy was converged to within 0.1 eV.

#### ch.sh

You can run the code by hand, changing the bond-length for each calculation, but you can also use the following script:

#!/bin/bash echo "#CH distance Total Energy" > ch.log list="0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30" for CH in $list do sed -ie "/CH = /s|[0-9.]\+|$CH|" inp octopus >& out-$CH energy=`grep Total static/info | head -1 | cut -d "=" -f 2` echo $CH $energy >> ch.log rm -rf restart done

#### Results

Once you run the calculations, you should obtain values very similar to the following ones:

#CH distance Total Energy 0.90 -215.22314794 0.95 -217.00685915 1.00 -218.09391896 1.05 -218.64710365 1.10 -218.78366713 1.15 -218.61270501 1.20 -218.20244112 1.25 -217.61312454 1.30 -216.8944146

You can see on the right how this curve looks like when plotted. By simple inspection we realize that the equilibrium CH distance that we obtain is around 1.1 Å. This should be compared to the experimental value of 1.094 Å. From this data can you calculate the frequency of the vibrational breathing mode of methane?

### FIRE

Since changing the bond-lengths by hand is not exactly the simplest method to perform the geometry optimization, we will now use automated algorithms for doing precisely that. Several of these algorithms are available in Octopus. The default and recommend one is the the FIRE method^{[1]}. Octopus can also use the GSL optimization routines, so all the methods provided by this library are available.

#### Input

First we will run a simple geometry optimization using the following input file:

`FromScratch`

= yes`CalculationMode`

= go`UnitsOutput`

= eV_Angstrom`BoxShape`

= sphere`Radius`

= 5.0*angstrom`Spacing`

= 0.18*angstrom CH = 1.2*angstrom %`Coordinates`

"C" | 0 | 0 | 0 "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) %

In this case we are using default values for all the input variables related to the geometry optimization. Compared to the input file above, we have changed the `CalculationMode`

to `go`. Because the default minimum box of spheres centered on the atoms has some caveats when performing a geometry optimization (see bellow for more details), we also changed the box shape to a single sphere with a 5 Å radius. Like always, convergence of the final results should be checked with respect to the mesh parameters. In this case you can trust us that a radius of 5 Å is sufficient.

The geometry in the input file will be used as a starting point. Note that although a single variable CH was used to specify the structure, the coordinates will be optimized independently.

#### Output

When running Octopus in geometry optimization mode the output will be a bit different from the normal ground-state mode. Since each minimization step requires several SCF cycles, only one line of information for each SCF iteration will printed:

... Info: Starting SCF iteration. ETA: .......1......2.......3......4......5.......6......7.......8......9......0 iter 1 : etot -2.27004788E+02 : abs_dens 1.96E+00 : etime 0.2s ETA: .......1......2.......3......4......5.......6......7.......8......9......0 ... iter 14 : etot -2.18209977E+02 : abs_dens 5.03E-06 : etime 0.2s ETA: .......1......2.......3......4......5.......6......7.......8......9......0 iter 15 : etot -2.18209908E+02 : abs_dens 1.70E-06 : etime 0.3s Info: Writing states. 2021/08/27 at 20:52:58 Info: Finished writing states. 2021/08/27 at 20:52:58 Info: SCF converged in 15 iterations ** Warning: ** Some of the states are not fully converged! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 1 ++++++++++++++++++++++ Energy = -218.2099082587 eV Max force = 2.5239574375 eV/A Max dr = 0.0002435249 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ...

At the end of each minimization step the code will print the total energy, the maximum force acting on the ions and the maximum displacement of the ions from the previous geometry to the next one. Octopus will also print the current geometry to a file named ** go.XXXX.xyz** (XXXX is the minimization iteration number) in a directory named

**. In the working directory there is also a file named**

`geom`**that accumulates the geometries of each SCF calculation. The results are also summarized in**

`work-geom.xyz`**. The most recent individual geometry is in**

`geom/optimization.log`**.**

`last.xyz`After some minimization steps the code will exit.

... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 32 ++++++++++++++++++++++ Energy = -218.7910310374 eV Max force = 0.0513885602 eV/A Max dr = 0.0052755143 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Writing final coordinates to min.xyz Info: Finished writing information to 'restart/gs'. ...

If the minimization was successful a file named ** min.xyz** will be written in the working directory containing the minimized geometry. In this case it should look like this:

5 units: A C -0.000000 -0.000000 -0.000000 H 0.633663 0.633663 0.633663 H -0.633663 -0.633663 0.633663 H 0.633663 -0.633663 -0.633663 H -0.633663 0.633663 -0.633663

Has the symmetry been retained? What is the equilibrium bond length determined? Check the forces in ** static/info**. You can also visualize the optimization process as an animation by loading

**as an XYZ file into XCrySDen.**

`work-geom.xyz`## Sodium trimer

Let's now try something different: a small sodium cluster. This time we will use an alternate method related to conjugate gradients.

#### Input

Here is the input file to use

`FromScratch`

= yes`CalculationMode`

= go`UnitsOutput`

= eV_Angstrom`Radius`

= 9.0*angstrom`Spacing`

= 0.3*angstrom`XYZCoordinates`

= "inp.xyz"`GOMethod`

= cg_bfgs

and the corresponding ** inp.xyz** file

3 Na -1.65104552837 -1.51912734276 -1.77061890357 Na -0.17818042320 -0.81880768425 -2.09945089196 Na 1.36469690888 0.914070162416 0.40877139068

In this case the starting geometry is a random geometry.

#### Output

After some time the code should stop successfully. The information about the last minimization iteration should look like this:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 3 ++++++++++++++++++++++ Energy = -16.8812926632 eV Max force = 0.0281998327 eV/A Max dr = 0.1322052321 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

You might have noticed that in this case the code might perform several SCF cycles for each minimization iteration. This is normal, as the minimization algorithm used here requires the evaluation of the energy and forces at several geometries per minimization iteration. The minimized geometry can be found in the ** min.xyz** file:

3 units: A Na -2.269092 -1.807901 -1.617872 Na 0.457329 -0.515180 -2.237397 Na 1.347055 0.899173 0.394161

#### Caveats of using a minimum box

For this example we are using the minimum box (` BoxShape = minimum`). This box shape is the one used by default, as it is usually the most efficient one, but it requires some care when using it for a geometry optimization. The reason is that the box is generated at the beginning of the calculation and is not updated when the atoms move. This means that at the end of the optimization, the centers of the spheres composing the box will not coincide anymore with the position of the atoms. This is a problem if the atoms move enough such that their initial and final positions differ significantly. If this is the case, the best is to restart the calculation using the previously minimized geometry as the starting geometry, so that the box is constructed again consistently with the atoms. In this example, this can be done by copying the contents of the

**file**

`min.xyz`**and running Octopus again without changing input file.**

`inp.xyz`This time the calculation should converge after only one iteration:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 1 ++++++++++++++++++++++ Energy = -16.8896899780 eV Max force = 0.0447462812 eV/A Max dr = 0.0450282495 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

You can now compare the old ** min.xyz** file with the new one:

3 units: A Na -2.252939 -1.788719 -1.588720 Na 0.484412 -0.498432 -2.232198 Na 1.302295 0.863198 0.362397

Although the differences are not too large, there are also not negligible.

#### Further improving the geometry

Imagine that you now want to improve the geometry obtained and reduce the value of the forces acting on the ions. The stopping criterion is controlled by the `GOTolerance`

variable (and `GOMinimumMove`

). Its default value is 0.051 eV/Å (0.001 H/b), so let's set it to 0.01 eV/Å:

`FromScratch`

= yes`CalculationMode`

= go`UnitsOutput`

= eV_Angstrom`Radius`

= 9.0*angstrom`Spacing`

= 0.3*angstrom`XYZCoordinates`

= "inp.xyz"`GOMethod`

= cg_bfgs`GOTolerance`

= 0.01*eV/angstrom

Copy the contents of file ** min.xyz** to file

**and rerun the code. This time, after some minimization steps Octopus should return the following error message:**

`inp.xyz`++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 21 ++++++++++++++++++++++ Energy = -16.9105971375 eV Max force = 0.0137311212 eV/A Max dr = 0.0000000000 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ **************************** FATAL ERROR ***************************** *** Fatal Error (description follows) *-------------------------------------------------------------------- * Error occurred during the GSL minimization procedure: * iteration is not making progress towards solution **********************************************************************

This means that, although the minimization algorithm managed to get very close (~ 0.014 eV/Å), it was unable to improve the solutions up to the desired accuracy. Unfortunately there are not many things that can be done in this case. The most simple one is just to restart the calculation starting from the last geometry (use the last ** go.XXXX.xyz** file) and see if the code is able to improve further the geometry.

## References

- ↑
Erik Bitzek, Pekka Koskinen, Franz Gähler, Michael Moseler, and Peter Gumbsch,
*Structural relaxation made simple*, Phys. Rev. Lett.**97**170201 (2006)