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
)

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)

using GeometryOptimization
results = minimize_energy!(h2_crude, calc; tol_forces=2e-6, verbosity=2)
nothing  # hide

results.system

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