# Tutorial:Geometry optimization

Note 1: FIRE method is implemented and recommended to perform geometry optimization. Note 2: For versions prior to 5.0.0 we would always recommend if possible to do your relaxations in a different DFT code, and then use the structure for Octopus.

## Methane molecule

To start we will go back to the Methane molecule from the ground-state tutorial.

Now that we have fixed the parameters of the mesh we can proceed and optimize the bond length. It is instructive to do it first "by hand" (well, if you wish you can modify the script we used for before to do this series of runs automatically). One gets:

# CH distance Total Energy 0.90 -215.22445018 0.95 -217.00783029 1.00 -218.09447045 1.05 -218.64691808 1.10 -218.78606597 1.15 -218.61467189 1.20 -218.20446189 1.25 -217.61518791 1.30 -216.89632549

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?

Changing the bond-lengths by hand is not exactly the simplest method to perform geometry optimization. Fortunately, there are a number of automated algorithms for doing precisely that. To use them, you should change `CalculationMode`

to `go`. Octopus relies on the GSL optimization routines, so all the methods provided by this library are available through the variable `GOMethod`

.

First we will run a simple geometry optimization using the default values by just changing the `CalculationMode`

to `go`. The input file should look like this:

`CalculationMode`

= go`UnitsOutput`

= eV_Angstrom`GOMethod`

= fire`ExperimentalFeatures`

= yes`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) %

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.

In the output, you will see this warning:

** Warning: ** GOMethod = fire is under development. ** It might not work or produce wrong results.

This is related to the fact we specified `ExperimentalFeatures`

= yes. Without that, the calculation would have stopped instead of writing a warning. Yes, the FIRE method^{[1]}
is experimental, but it is pretty safe here in the sense that the forces will be correct regardless, so if anything goes wrong you will know.
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:

... iter 21 : etot -2.18202487E+02 : abs_dens 2.00E-07 : force 4.98E-07 : etime 0.7s Info: Writing states. 2016/09/12 at 15:31:08 Info: Finished writing states. 2016/09/12 at 15:31:08 Info: SCF converged in 21 iterations ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 1 ++++++++++++++++++++++ Energy = -218.2024874817 eV Max force = 2.5304116861 eV/A Max dr = 0.0708697876 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 #: 34 ++++++++++++++++++++++ Energy = -218.7861551998 eV Max force = 0.0396564762 eV/A Max dr = 0.0005971382 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.633363 0.633363 0.633363 H -0.633363 -0.633363 0.633363 H 0.633363 -0.633363 -0.633363 H -0.633363 0.633363 -0.633363

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`## Na

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

`CalculationMode`

= go`UnitsOutput`

= eV_Angstrom`radius`

= 9.0*angstrom`spacing`

= 0.3*angstrom`XYZCoordinates`

= "inp.xyz"`GOMethod`

= cg_bfgs

and start from a random geometry (file ** inp.xyz**):

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

After some time the code should stop successfully:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 3 ++++++++++++++++++++++ Energy = -17.0614891400 eV Max force = 0.0351654788 eV/A Max dr = 0.1529689428 A ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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`

), so let's set it to 0.01 and rerun the code (don't forget to copy the contents of file ** min.xyz** to file

**!). This time, after some minimization steps Octopus should return the following error message:**

`inp.xyz`**************************** FATAL ERROR ***************************** *** Fatal Error (description follows) *-------------------------------------------------------------------- * Error occurred during the GSL minimization procedure: * iteration is not making progress towards solution **********************************************************************

This means the minimization algorithm 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.

If you have ` BoxShape = minimum`, it is a good idea to restart at least once before accepting the final structure, because force calculation will be compromised if the atoms move enough during the relaxation process that the minimum mesh for the new positions differs significantly from the mesh for the original positions. But when you restart, the mesh is constructed again consistently with the atoms.

You can also try `GOObjective`

= minimize_forces with `GOMethod`

= steep.

It is recommended to use `FilterPotentials`

= filter_TS for geometry optimization to reduce eggbox effects (which is already set on by default)

## References

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