In this tutorial we will show how to perform geometry optimizations using Octopus.
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.
= gs = eV_Angstrom = 3.5*angstrom = 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) %
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.
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
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?
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. Octopus can also use the GSL optimization routines, so all the methods provided by this library are available.
First we will run a simple geometry optimization using the following input file:
= yes = go = eV_Angstrom = sphere = 5.0*angstrom = 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) %
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 theto 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.
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 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 #: 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 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.
Here is the input file to use
= yes = go = eV_Angstrom = 9.0*angstrom = 0.3*angstrom = "inp.xyz" = 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.
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 (= 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 min.xyz file inp.xyz and running Octopus again without changing input file.
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 thevariable (and ). Its default value is 0.051 eV/Å (0.001 H/b), so let's set it to 0.01 eV/Å:
= yes = go = eV_Angstrom = 9.0*angstrom = 0.3*angstrom = "inp.xyz" = cg_bfgs = 0.01*eV/angstrom
Copy the contents of file min.xyz to file inp.xyz and rerun the code. This time, after some minimization steps Octopus should return the following error message:
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ 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.
- Erik Bitzek, Pekka Koskinen, Franz Gähler, Michael Moseler, and Peter Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97 170201 (2006)