Geometry optimization

We use DFTK and the GeometryOptimization package to find the minimal-energy bond length of the $H_2$ molecule. First we set up an appropriate DFTKCalculator (see AtomsCalculators integration), for which we use the LDA model just like in the Tutorial for silicon in combination with a pseudodojo pseudopotential (see Pseudopotentials).

using DFTK
using PseudoPotentialData

pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
calc = DFTKCalculator(;
    model_kwargs = (; functionals=LDA(), pseudopotentials),  # model_DFT keyword arguments
    basis_kwargs = (; kgrid=[1, 1, 1], Ecut=10)  # PlaneWaveBasis keyword arguments
)
DFTKCalculator(functionals=Xc(lda_x, lda_c_pw), pseudopotentials=PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf"), Ecut=10, kgrid=[1, 1, 1])

Next we set up an initial hydrogen molecule within a box of vacuum. We use the parameters of the equivalent tutorial from ABINIT and DFTK's AtomsBase integration to setup the hydrogen molecule. We employ a simulation box of 10 bohr times 10 bohr times 10 bohr and a pseudodojo pseudopotential.

using LinearAlgebra
using Unitful
using UnitfulAtomic

r0 = 1.4   # Initial bond length in Bohr
a  = 10.0  # Box size in Bohr

cell_vectors = [[a, 0, 0]u"bohr", [0, a, 0]u"bohr", [0, 0, a]u"bohr"]
h2_crude = periodic_system([:H => [0, 0, 0.]u"bohr",
                            :H => [r0, 0, 0]u"bohr"],
                           cell_vectors)
FlexibleSystem(H₂, periodicity = TTT):
    cell_vectors      : [      10        0        0;
                                0       10        0;
                                0        0       10]u"a₀"

    Atom(H,  [       0,        0,        0]u"a₀")
    Atom(H,  [     1.4,        0,        0]u"a₀")

Finally we call minimize_energy! to start the geometry optimisation. We use verbosity=2 to get some insight into the minimisation. With verbosity=1 only a summarising table would be printed and with verbosity=0 (default) the minimisation would be quiet.

using GeometryOptimization
results = minimize_energy!(h2_crude, calc; tol_forces=2e-6, verbosity=2)
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.111214155265                   -0.83   0.80    9.0   26.5ms
  2   -1.117209536450       -2.22       -2.01   0.80    1.0   46.7ms
  3   -1.117250517502       -4.39       -2.71   0.80    1.0   17.9ms
  4   -1.117251163626       -6.19       -3.56   0.80    1.0   18.2ms
  5   -1.117251189451       -7.59       -4.54   0.80    1.0   18.4ms
  6   -1.117251190088       -9.20       -4.87   0.80    2.0   21.0ms
  7   -1.117251190142      -10.27       -5.74   0.80    1.0   19.4ms
  8   -1.117251190143      -12.23       -6.13   0.80    1.0   35.6ms
  9   -1.117251190143      -12.73       -7.05   0.80    1.0   20.4ms
 10   -1.117251190143      -14.54       -8.01   0.80    1.0   20.0ms
 11   -1.117251190143      -14.88       -8.81   0.80    2.0   22.1ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.117251190143                   -9.32   0.80    1.0    116ms
  2   -1.117251190143      -15.65      -10.41   0.80    1.0   19.5ms

Geometry optimisation convergence (in atomic units)
┌─────┬─────────────────┬───────────┬─────────────┬────────┐
│   n           Energy  log10(ΔE)   max(Force)   Δtime │
├─────┼─────────────────┼───────────┼─────────────┼────────┤
│   0 │ -1.117251190143 │           │   0.0269318 │  1.33s │
└─────┴─────────────────┴───────────┴─────────────┴────────┘

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.117517998060                   -2.46   0.80    1.0   11.2ms
  2   -1.117520586051       -5.59       -3.49   0.80    1.0   24.3ms
  3   -1.117520612985       -7.57       -4.23   0.80    1.0   18.6ms
  4   -1.117520613449       -9.33       -5.14   0.80    1.0   19.1ms
  5   -1.117520613464      -10.81       -5.92   0.80    1.0   19.2ms
  6   -1.117520613465      -12.13       -6.42   0.80    1.0   23.4ms
  7   -1.117520613465      -13.32       -7.39   0.80    1.0   19.9ms
  8   -1.117520613465      -14.75       -8.18   0.80    1.0   20.3ms
  9   -1.117520613465   +  -14.65       -8.97   0.80    1.0   25.8ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118247264592                   -1.69   0.80    1.0   11.5ms
  2   -1.118339043934       -4.04       -2.70   0.80    1.0   18.1ms
  3   -1.118340011113       -6.01       -3.44   0.80    1.0   24.5ms
  4   -1.118340028554       -7.76       -4.36   0.80    1.0   18.8ms
  5   -1.118340029138       -9.23       -5.12   0.80    1.0   19.4ms
  6   -1.118340029165      -10.57       -5.62   0.80    1.0    113ms
  7   -1.118340029167      -11.73       -6.68   0.80    1.0   20.8ms
  8   -1.118340029167      -13.57       -7.58   0.80    1.0   20.8ms
  9   -1.118340029167      -15.35       -8.12   0.80    1.0   20.5ms
 10   -1.118340029167   +  -14.88       -8.55   0.80    1.0    927ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118340029167                   -9.05   0.80    1.0   11.3ms
  2   -1.118340029167   +    -Inf      -10.36   0.80    1.0   18.1ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118340029167                  -11.04   0.80    1.0   11.2ms
  2   -1.118340029167      -14.70      -11.41   0.80    1.0   18.2ms

┌─────┬─────────────────┬───────────┬─────────────┬────────┐
│   n           Energy  log10(ΔE)   max(Force)   Δtime │
├─────┼─────────────────┼───────────┼─────────────┼────────┤
│   1 │ -1.118340029167 │           │  0.00297055 │  1.72s │
└─────┴─────────────────┴───────────┴─────────────┴────────┘

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118346721507                   -3.09   0.80    1.0   11.3ms
  2   -1.118346864759       -6.84       -4.11   0.80    1.0   17.9ms
  3   -1.118346866236       -8.83       -4.86   0.80    1.0   18.5ms
  4   -1.118346866262      -10.59       -5.78   0.80    1.0   18.9ms
  5   -1.118346866262      -12.08       -6.55   0.80    1.0   18.9ms
  6   -1.118346866262      -13.44       -7.04   0.80    1.0   19.5ms
  7   -1.118346866262      -15.65       -8.13   0.80    1.0   20.0ms
  8   -1.118346866262      -14.75       -8.95   0.80    1.0   19.9ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118354724004                   -2.61   0.80    1.0   11.0ms
  2   -1.118356087056       -5.87       -3.62   0.80    1.0   17.9ms
  3   -1.118356101131       -7.85       -4.37   0.80    1.0   18.1ms
  4   -1.118356101377       -9.61       -5.29   0.80    1.0   18.5ms
  5   -1.118356101385      -11.10       -6.06   0.80    1.0   18.9ms
  6   -1.118356101386      -12.45       -6.55   0.80    1.0   19.0ms
  7   -1.118356101386      -13.57       -7.65   0.80    1.0   19.8ms
  8   -1.118356101386      -15.05       -8.45   0.80    1.0   19.7ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356101386                   -9.24   0.80    1.0   11.1ms
  2   -1.118356101386      -15.35       -9.82   0.80    1.0   24.3ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356101386                  -10.62   0.80    1.0   24.0ms
  2   -1.118356101386      -14.54      -11.46   0.80    1.0   18.4ms

┌─────┬─────────────────┬───────────┬─────────────┬────────┐
│   n           Energy  log10(ΔE)   max(Force)   Δtime │
├─────┼─────────────────┼───────────┼─────────────┼────────┤
│   2 │ -1.118356101386 │     -4.79 │  4.47268e-5 │  609ms │
└─────┴─────────────────┴───────────┴─────────────┴────────┘

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356102848                   -4.94   0.80    1.0   11.4ms
  2   -1.118356102877      -10.54       -5.96   0.80    1.0   18.0ms
  3   -1.118356102877      -12.52       -6.70   0.80    1.0   18.1ms
  4   -1.118356102877      -14.33       -7.62   0.80    1.0   18.6ms
  5   -1.118356102877      -15.35       -8.40   0.80    1.0   19.0ms
  6   -1.118356102877   +  -15.05       -8.89   0.80    1.0   19.3ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356104771                   -4.40   0.80    1.0   11.3ms
  2   -1.118356105116       -9.46       -5.42   0.80    1.0   18.3ms
  3   -1.118356105120      -11.45       -6.17   0.80    1.0   18.7ms
  4   -1.118356105120      -13.18       -7.09   0.80    1.0   18.8ms
  5   -1.118356105120      -14.48       -7.86   0.80    1.0   20.0ms
  6   -1.118356105120   +  -15.05       -8.35   0.80    1.0   19.3ms
  7   -1.118356105120   +  -15.35       -9.45   0.80    1.0   27.5ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356105120                  -10.08   0.80    1.0   11.3ms
  2   -1.118356105120   +  -14.24      -11.02   0.80    1.0   18.2ms
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -1.118356105120                  -11.67   0.80    1.0   11.1ms
  2   -1.118356105120   +  -14.33      -12.44   0.80    1.0   18.1ms

┌─────┬─────────────────┬───────────┬─────────────┬────────┐
│   n           Energy  log10(ΔE)   max(Force)   Δtime │
├─────┼─────────────────┼───────────┼─────────────┼────────┤
│   3 │ -1.118356105120 │     -8.43 │  1.11492e-8 │  533ms │
└─────┴─────────────────┴───────────┴─────────────┴────────┘

Structure after optimisation (note that the atom has wrapped around)

results.system
FlexibleSystem(H₂, periodicity = TTT):
    cell_vectors      : [      10        0        0;
                                0       10        0;
                                0        0       10]u"a₀"

    Atom(H,  [-0.0431829, 4.45861e-12, 2.32107e-12]u"a₀")
    Atom(H,  [ 1.44318, 4.30856e-12, 2.25565e-12]u"a₀")

Compute final bond length:

rmin = norm(position(results.system[1]) - position(results.system[2]))
println("Optimal bond length: ", rmin)
Optimal bond length: 1.486365864270879 a₀

Our results (1.486 Bohr) agrees with the equivalent tutorial from ABINIT.