Optimize the bond length of a hydrogen molecule

This section shows how we can use KSSOLV2.0 to optimize the geometry of the simplest molecule, a hydrogen molecule with two hydrogen atoms. This example is provided in the examples directory as the m-file h2_scf_opt.m

In the first part of this script, we set up the atomic species, unit supercell containing the molecule, and geometry of the molecular. In particular, we define the (x,y,z) coordinates of the hydrogen atoms as

xyzlist = [
  1.5     0.0      0.0
  0.0     0.0      0.0
];

The bond length is 1.5 Bohr, which is not the optimal length (1.4 Bohr).

We set the planewave cutoff energy to 100.0 Hartree for this calculation.

mol = Molecule('supercell',C, 'atomlist',atomlist, 'xyzlist' ,xyzlist, ...
               'ecutwfc',100.0, 'name','h2' );

We perform a baseline KSDFT calculation at the non-optimal bond length to obtain the ground state energy associated with this geometry

[mol,H,X,info] = scf(mol);

After 10 SCF iteration, convergence is reached (up to \(10^{-7}\) default error tolerance in potential).

...
SCF iter  10:
Rel Vtot Err    =            5.453e-08
Total Energy    = -1.1673374036056e+00

Note that the total energy of the molecule in this geometry is -1.1673374036056 Hartree.

We then call the relaxatom function to optimize the geometry of the atoms. Before calling this function, we set some of the parameters used in both the ground state DFT calculation and geometry optimization. In particular, we choose to use the MATLAB unconstrained minimization solver fminunc (which is part of MATLAB’s optimization toolbox) to perform the geometry optimization, and set the convergence tolerance to be \(10^{-4}\) for the gradient norm.

ksopts = setksopt;
ksopts = setksopt(ksopts,'scftol',1e-7,...
                 'cgtol',1e-8,...
                 'maxscfiter',100, ...
                 'maxcgiter',100, ...
                 'relaxmethod','fminunc', ...
                 'relaxtol', 1e-4);

relaxatom performs a sequence of SCF calculations to evaluate the total energy (objective function) and forces (gradient of the objective function with respect to the atomic coordinates) associated with each geometry configuration as it is optimized (relaxed). We can check the total energy and gradient norm for each of these runs, for example

...
Starting SCF4M calculation for h2...
...
||HX-XD||_F     =            7.239e-09
fval = -1.168262761165295e+00
norm(g) =   2.042e-03
...

When relaxatom is terminated, we see

||HX-XD||_F     =            7.418e-09
fval = -1.168265790996790e+00
norm(g) =   7.318e-07

Local minimum found.

Optimization completed because the size of the gradient is less than
the selected value of the optimality tolerance.

We then print out the bond length for the optimized geometry

optimized bond length =   1.423e+00 (Bohr)

This is much closer to the true optimal bond length of 1.4 Bohr. It is not completely optimized because the planewave energy cutoff we used is not sufficiently large.

The total energy associated with this geometry is

Etot            = -1.1682657909968e+00

which is less than the total energy -1.1673374036056 that we started with.