We use the DFTK and Optim packages in this example to find the minimal-energy bond length of the $H_2$ molecule. We setup $H_2$ in an LDA model just like in the Tutorial for silicon.
using DFTK using Optim using LinearAlgebra using Printf kgrid = [1, 1, 1] # k-point grid Ecut = 5 # kinetic energy cutoff in Hartree tol = 1e-8 # tolerance for the optimization routine a = 10 # lattice constant in Bohr lattice = a * I(3) H = ElementPsp(:H, psp=load_psp("hgh/lda/h-q1")); atoms = [H, H];
We define a Bloch wave and a density to be used as global variables so that we can transfer the solution from one iteration to another and therefore reduce the optimization time.
ψ = nothing ρ = nothing
First, we create a function that computes the solution associated to the
position $x ∈ ℝ^6$ of the atoms in reduced coordinates
(cf. Reduced and cartesian coordinates for more
details on the coordinates system).
They are stored as a vector:
x[1:3] represents the position of the
first atom and
x[4:6] the position of the second.
We also update
ρ for the next iteration.
function compute_scfres(x) model = model_LDA(lattice, atoms, [x[1:3], x[4:6]]) basis = PlaneWaveBasis(model; Ecut, kgrid) global ψ, ρ if isnothing(ρ) ρ = guess_density(basis) end scfres = self_consistent_field(basis; ψ, ρ, tol=tol / 10, callback=identity) ψ = scfres.ψ ρ = scfres.ρ scfres end;
Then, we create the function we want to optimize:
fg! is used to update the
value of the objective function
F, namely the energy, and its gradient
here computed with the forces (which are, by definition, the negative gradient
of the energy).
function fg!(F, G, x) scfres = compute_scfres(x) if G != nothing grad = compute_forces(scfres) G .= -[grad; grad] end scfres.energies.total end;
Now, we can optimize on the 6 parameters
x = [x1, y1, z1, x2, y2, z2] in
reduced coordinates, using
LBFGS(), the default minimization algorithm
in Optim. We start from
x0, which is a first guess for the coordinates. By
optimize traces the output of the optimization algorithm during the
iterations. Once we have the minimizer
xmin, we compute the bond length in
x0 = vcat(lattice \ [0., 0., 0.], lattice \ [1.4, 0., 0.]) xres = optimize(Optim.only_fg!(fg!), x0, LBFGS(), Optim.Options(show_trace=true, f_tol=tol)) xmin = Optim.minimizer(xres) dmin = norm(lattice*xmin[1:3] - lattice*xmin[4:6]) @printf "\nOptimal bond length for Ecut=%.2f: %.3f Bohr\n" Ecut dmin
Iter Function value Gradient norm 0 -1.061651e+00 6.219595e-01 * time: 5.1975250244140625e-5 1 -1.064076e+00 2.917716e-01 * time: 2.589946985244751 2 -1.066008e+00 4.814326e-02 * time: 3.065577983856201 3 -1.066049e+00 4.063030e-04 * time: 3.3044850826263428 4 -1.066049e+00 4.521245e-06 * time: 3.511681079864502 5 -1.066049e+00 1.667006e-07 * time: 3.7725119590759277 Optimal bond length for Ecut=5.00: 1.556 Bohr
We used here very rough parameters to generate the example and
Ecut to 10 Ha yields a bond length of 1.523 Bohr,
which agrees with ABINIT.
!!! note "Degrees of freedom" We used here a very general setting where we optimized on the 6 variables representing the position of the 2 atoms and it can be easily extended to molecules with more atoms (such as $H_2O$). In the particular case of $H_2$, we could use only the internal degree of freedom which, in this case, is just the bond length.