Tutorial:Geometry optimization

From OctopusWiki
Jump to: navigation, search

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

Methane Energy vs CH.png

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
  "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 geom. In the working directory there is also a file named work-geom.xyz that accumulates the geometries of each SCF calculation. The results are also summarized in geom/optimization.log. The most recent individual geometry is in 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:

 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 work-geom.xyz as an XYZ file into XCrySDen.


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):


   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 inp.xyz!). This time, after some minimization steps Octopus should return the following error message:

**************************** 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)


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

Previous Tutorial:Periodic systems - Next Tutorial:Basic QOCT

Back to Tutorial