We use DFTK and the GeometryOptimization
package to find the minimal-energy bond length of the H2 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)
nothing # hide
n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.110282381972 -0.82 0.80 8.0 4.08s 2 -1.117191702987 -2.16 -1.82 0.80 1.0 2.73s 3 -1.117250264002 -4.23 -2.66 0.80 1.0 19.3ms 4 -1.117251216249 -6.02 -3.30 0.80 1.0 21.9ms 5 -1.117251300360 -7.08 -3.77 0.80 1.0 33.4ms 6 -1.117251310001 -8.02 -4.86 0.80 1.0 19.3ms 7 -1.117251310375 -9.43 -5.14 0.80 2.0 23.8ms 8 -1.117251310379 -11.34 -5.59 0.80 1.0 20.4ms 9 -1.117251310381 -11.80 -6.54 0.80 1.0 23.3ms 10 -1.117251310381 -13.30 -7.26 0.80 2.0 25.6ms 11 -1.117251310381 -15.18 -7.62 0.80 1.0 21.8ms 12 -1.117251310381 + -15.05 -7.90 0.80 1.0 22.2ms 13 -1.117251310381 -14.35 -8.59 0.80 1.0 21.9ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.117251310381 -8.94 0.80 1.0 2.13s 2 -1.117251310381 + -14.61 -9.91 0.80 1.0 18.5ms Geometry optimisation convergence (in atomic units) ┌─────┬─────────────────┬───────────┬────────────┬────────┐ │ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │ ├─────┼─────────────────┼───────────┼────────────┼────────┤ │ 0 │ -1.117251310381 │ │ 0.0269317 │ 22.4s │ └─────┴─────────────────┴───────────┴────────────┴────────┘ n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.117518116179 -2.46 0.80 1.0 38.1ms 2 -1.117520704471 -5.59 -3.49 0.80 1.0 16.6ms 3 -1.117520731407 -7.57 -4.23 0.80 1.0 17.6ms 4 -1.117520731867 -9.34 -5.14 0.80 1.0 31.7ms 5 -1.117520731883 -10.81 -5.93 0.80 1.0 17.5ms 6 -1.117520731883 -12.16 -6.43 0.80 1.0 17.9ms 7 -1.117520731883 -13.37 -7.42 0.80 1.0 20.5ms 8 -1.117520731883 + -Inf -8.25 0.80 1.0 20.9ms 9 -1.117520731883 + -14.57 -8.88 0.80 1.0 18.9ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118247378325 -1.69 0.80 1.0 13.5ms 2 -1.118339157117 -4.04 -2.70 0.80 1.0 16.4ms 3 -1.118340124290 -6.01 -3.44 0.80 1.0 18.5ms 4 -1.118340141731 -7.76 -4.36 0.80 1.0 19.0ms 5 -1.118340142315 -9.23 -5.12 0.80 1.0 19.2ms 6 -1.118340142342 -10.57 -5.62 0.80 1.0 19.6ms 7 -1.118340142344 -11.73 -6.68 0.80 1.0 18.2ms 8 -1.118340142344 -13.58 -7.58 0.80 1.0 20.7ms 9 -1.118340142344 + -15.05 -8.12 0.80 1.0 20.9ms 10 -1.118340142344 -14.95 -8.55 0.80 1.0 21.1ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118340142344 -9.05 0.80 1.0 10.5ms 2 -1.118340142344 + -14.88 -10.36 0.80 1.0 22.5ms ┌─────┬─────────────────┬───────────┬────────────┬────────┐ │ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │ ├─────┼─────────────────┼───────────┼────────────┼────────┤ │ 1 │ -1.118340142344 │ -2.96 │ 0.00297053 │ 1.72s │ └─────┴─────────────────┴───────────┴────────────┴────────┘ n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118346834631 -3.09 0.80 1.0 10.5ms 2 -1.118346977882 -6.84 -4.11 0.80 1.0 22.2ms 3 -1.118346979359 -8.83 -4.86 0.80 1.0 16.8ms 4 -1.118346979384 -10.59 -5.78 0.80 1.0 19.4ms 5 -1.118346979385 -12.08 -6.55 0.80 1.0 17.4ms 6 -1.118346979385 -13.43 -7.04 0.80 1.0 19.8ms 7 -1.118346979385 -14.75 -8.13 0.80 1.0 20.0ms 8 -1.118346979385 -14.57 -8.95 0.80 1.0 20.3ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118354837055 -2.61 0.80 1.0 10.3ms 2 -1.118356200096 -5.87 -3.62 0.80 1.0 18.5ms 3 -1.118356214170 -7.85 -4.37 0.80 1.0 16.5ms 4 -1.118356214417 -9.61 -5.29 0.80 1.0 74.0ms 5 -1.118356214425 -11.10 -6.06 0.80 1.0 394ms 6 -1.118356214425 -12.45 -6.55 0.80 1.0 17.8ms 7 -1.118356214425 -13.53 -7.65 0.80 1.0 18.1ms 8 -1.118356214425 + -14.75 -8.45 0.80 1.0 18.7ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118356214425 -9.24 0.80 1.0 12.0ms 2 -1.118356214425 -14.65 -9.82 0.80 1.0 22.8ms ┌─────┬─────────────────┬───────────┬────────────┬────────┐ │ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │ ├─────┼─────────────────┼───────────┼────────────┼────────┤ │ 2 │ -1.118356214425 │ -4.79 │ 4.47263e-5 │ 949ms │ └─────┴─────────────────┴───────────┴────────────┴────────┘ n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118356215887 -4.94 0.80 1.0 11.6ms 2 -1.118356215917 -10.54 -5.96 0.80 1.0 23.1ms 3 -1.118356215917 -12.52 -6.70 0.80 1.0 19.8ms 4 -1.118356215917 -14.31 -7.62 0.80 1.0 18.7ms 5 -1.118356215917 -14.81 -8.39 0.80 1.0 22.3ms 6 -1.118356215917 -15.05 -8.89 0.80 1.0 25.9ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118356217810 -4.40 0.80 1.0 10.4ms 2 -1.118356218156 -9.46 -5.42 0.80 1.0 17.0ms 3 -1.118356218159 -11.45 -6.17 0.80 1.0 16.7ms 4 -1.118356218159 -13.21 -7.09 0.80 1.0 17.6ms 5 -1.118356218159 -14.65 -7.86 0.80 1.0 17.6ms 6 -1.118356218159 -14.75 -8.35 0.80 1.0 17.8ms 7 -1.118356218159 + -15.65 -9.45 0.80 1.0 21.8ms n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -1.118356218159 -10.08 0.80 1.0 10.4ms 2 -1.118356218159 + -15.35 -11.02 0.80 1.0 16.5ms ┌─────┬─────────────────┬───────────┬────────────┬────────┐ │ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │ ├─────┼─────────────────┼───────────┼────────────┼────────┤ │ 3 │ -1.118356218159 │ -8.43 │ 1.04318e-8 │ 452ms │ └─────┴─────────────────┴───────────┴────────────┴────────┘
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.0431828, -7.14517e-11, -3.1471e-10]u"a₀")
Atom(H, [ 1.44318, -6.98109e-11, -3.10227e-10]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.4863655847396426 a₀
Our results (1.486 Bohr) agrees with the equivalent tutorial from ABINIT.