Tutorial:Geometry optimization
From OctopusWiki
Note: Octopus is not particularly optimized with regard to geometry optimization. You can do it, nevertheless, at your own risk ;)
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. One gets:
CH distance Total Energy 0.90 -215.24450503 0.95 -217.02493406 1.00 -218.12498285 1.05 -218.60328522 1.10 -218.71401944 1.15 -218.57738023 1.20 -218.21486251 1.25 -217.62467183 1.30 -216.91046026
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= goUnits= eV_Angstromradius= 3.5spacing= 0.175 CH = 1.2 %Coordinates"C" | 0 | 0 | 0 | no "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) | no "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) | no "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) | no "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) | no %
The geometry in the input file will be used as a starting point.
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 iter 1 : etot -2.19100040E+02 : abs_dens 2.57E-01 : etime 1.9s iter 2 : etot -2.18572949E+02 : abs_dens 1.55E-01 : etime 2.7s ... iter 7 : etot -2.18214852E+02 : abs_dens 1.93E-05 : etime 1.9s iter 8 : etot -2.18214852E+02 : abs_dens 3.48E-06 : etime 0.9s Info: SCF converged in 8 iterations Info: Starting SCF iteration iter 1 : etot -2.21092641E+02 : abs_dens 2.95E-01 : etime 2.0s iter 2 : etot -2.18895408E+02 : abs_dens 1.06E-01 : etime 2.2s ... iter 11 : etot -2.18676099E+02 : abs_dens 1.53E-05 : etime 0.8s iter 12 : etot -2.18676099E+02 : abs_dens 5.17E-06 : etime 0.9s Info: SCF converged in 12 iterations ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++ MINIMIZATION ITER #: 1 ++++++++++++++++++++++ Energy = -218.6762601569 Max force = 0.5711774044 Max dr = 0.0763801209 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ...
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 contains the geometry of the current SCF cycle.
After some minimization steps the code will exit. 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
C 0.000000 0.000000 0.000000
H 0.631716 0.631716 0.631716
H -0.631716 -0.631716 0.631716
H 0.631716 -0.631716 -0.631716
H -0.631716 0.631716 -0.631716
and so we find an equilibrium CH distance of 1.0942 Å.
Na3
Let's now try something different: a small sodium cluster. This time we will use a more sophisticated method:
CalculationMode= goUnits= eV_Angstromradius= 9.0spacing= 0.3XYZCoordinates= "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 = -16.8992796305
Max force = 0.0255296037
Max dr = 0.1337945103
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Calculation ended on 2008/01/19 at 16:11:32
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, 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.
