{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Geometry optimization\n",
"\n",
"We use DFTK and the [GeometryOptimization](https://github.com/JuliaMolSim/GeometryOptimization.jl/)\n",
"package to find the minimal-energy bond length of the $H_2$ molecule.\n",
"First we set up an appropriate `DFTKCalculator` (see AtomsCalculators integration),\n",
"for which we use the LDA model just like in the Tutorial for silicon\n",
"in combination with a pseudodojo pseudopotential (see Pseudopotentials)."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "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])"
},
"metadata": {},
"execution_count": 1
}
],
"cell_type": "code",
"source": [
"using DFTK\n",
"using PseudoPotentialData\n",
"\n",
"pseudopotentials = PseudoFamily(\"dojo.nc.sr.pbe.v0_4_1.standard.upf\")\n",
"calc = DFTKCalculator(;\n",
" model_kwargs = (; functionals=LDA(), pseudopotentials), # model_DFT keyword arguments\n",
" basis_kwargs = (; kgrid=[1, 1, 1], Ecut=10) # PlaneWaveBasis keyword arguments\n",
")"
],
"metadata": {},
"execution_count": 1
},
{
"cell_type": "markdown",
"source": [
"Next we set up an initial hydrogen molecule within a box of vacuum.\n",
"We use the parameters of the\n",
"[equivalent tutorial from ABINIT](https://docs.abinit.org/tutorial/base1/)\n",
"and DFTK's AtomsBase integration to setup the hydrogen molecule.\n",
"We employ a simulation box of 10 bohr times 10 bohr times 10 bohr and a\n",
"pseudodojo pseudopotential."
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "FlexibleSystem(H₂, periodicity = TTT):\n cell_vectors : [ 10 0 0;\n 0 10 0;\n 0 0 10]u\"a₀\"\n\n Atom(H, [ 0, 0, 0]u\"a₀\")\n Atom(H, [ 1.4, 0, 0]u\"a₀\")\n"
},
"metadata": {},
"execution_count": 2
}
],
"cell_type": "code",
"source": [
"using LinearAlgebra\n",
"using Unitful\n",
"using UnitfulAtomic\n",
"\n",
"r0 = 1.4 # Initial bond length in Bohr\n",
"a = 10.0 # Box size in Bohr\n",
"\n",
"cell_vectors = [[a, 0, 0]u\"bohr\", [0, a, 0]u\"bohr\", [0, 0, a]u\"bohr\"]\n",
"h2_crude = periodic_system([:H => [0, 0, 0.]u\"bohr\",\n",
" :H => [r0, 0, 0]u\"bohr\"],\n",
" cell_vectors)"
],
"metadata": {},
"execution_count": 2
},
{
"cell_type": "markdown",
"source": [
"Finally we call `minimize_energy!` to start the geometry optimisation.\n",
"We use `verbosity=2` to get some insight into the minimisation.\n",
"With `verbosity=1` only a summarising table would be printed and with\n",
"`verbosity=0` (default) the minimisation would be quiet."
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.110282381972 -0.82 0.80 8.0 4.08s\n",
" 2 -1.117191702987 -2.16 -1.82 0.80 1.0 2.73s\n",
" 3 -1.117250264002 -4.23 -2.66 0.80 1.0 19.3ms\n",
" 4 -1.117251216249 -6.02 -3.30 0.80 1.0 21.9ms\n",
" 5 -1.117251300360 -7.08 -3.77 0.80 1.0 33.4ms\n",
" 6 -1.117251310001 -8.02 -4.86 0.80 1.0 19.3ms\n",
" 7 -1.117251310375 -9.43 -5.14 0.80 2.0 23.8ms\n",
" 8 -1.117251310379 -11.34 -5.59 0.80 1.0 20.4ms\n",
" 9 -1.117251310381 -11.80 -6.54 0.80 1.0 23.3ms\n",
" 10 -1.117251310381 -13.30 -7.26 0.80 2.0 25.6ms\n",
" 11 -1.117251310381 -15.18 -7.62 0.80 1.0 21.8ms\n",
" 12 -1.117251310381 + -15.05 -7.90 0.80 1.0 22.2ms\n",
" 13 -1.117251310381 -14.35 -8.59 0.80 1.0 21.9ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.117251310381 -8.94 0.80 1.0 2.13s\n",
" 2 -1.117251310381 + -14.61 -9.91 0.80 1.0 18.5ms\n",
"\n",
"Geometry optimisation convergence (in atomic units)\n",
"┌─────┬─────────────────┬───────────┬────────────┬────────┐\n",
"│ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │\n",
"├─────┼─────────────────┼───────────┼────────────┼────────┤\n",
"│ 0 │ -1.117251310381 │ │ 0.0269317 │ 22.4s │\n",
"└─────┴─────────────────┴───────────┴────────────┴────────┘\n",
"\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.117518116179 -2.46 0.80 1.0 38.1ms\n",
" 2 -1.117520704471 -5.59 -3.49 0.80 1.0 16.6ms\n",
" 3 -1.117520731407 -7.57 -4.23 0.80 1.0 17.6ms\n",
" 4 -1.117520731867 -9.34 -5.14 0.80 1.0 31.7ms\n",
" 5 -1.117520731883 -10.81 -5.93 0.80 1.0 17.5ms\n",
" 6 -1.117520731883 -12.16 -6.43 0.80 1.0 17.9ms\n",
" 7 -1.117520731883 -13.37 -7.42 0.80 1.0 20.5ms\n",
" 8 -1.117520731883 + -Inf -8.25 0.80 1.0 20.9ms\n",
" 9 -1.117520731883 + -14.57 -8.88 0.80 1.0 18.9ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118247378325 -1.69 0.80 1.0 13.5ms\n",
" 2 -1.118339157117 -4.04 -2.70 0.80 1.0 16.4ms\n",
" 3 -1.118340124290 -6.01 -3.44 0.80 1.0 18.5ms\n",
" 4 -1.118340141731 -7.76 -4.36 0.80 1.0 19.0ms\n",
" 5 -1.118340142315 -9.23 -5.12 0.80 1.0 19.2ms\n",
" 6 -1.118340142342 -10.57 -5.62 0.80 1.0 19.6ms\n",
" 7 -1.118340142344 -11.73 -6.68 0.80 1.0 18.2ms\n",
" 8 -1.118340142344 -13.58 -7.58 0.80 1.0 20.7ms\n",
" 9 -1.118340142344 + -15.05 -8.12 0.80 1.0 20.9ms\n",
" 10 -1.118340142344 -14.95 -8.55 0.80 1.0 21.1ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118340142344 -9.05 0.80 1.0 10.5ms\n",
" 2 -1.118340142344 + -14.88 -10.36 0.80 1.0 22.5ms\n",
"\n",
"┌─────┬─────────────────┬───────────┬────────────┬────────┐\n",
"│ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │\n",
"├─────┼─────────────────┼───────────┼────────────┼────────┤\n",
"│ 1 │ -1.118340142344 │ -2.96 │ 0.00297053 │ 1.72s │\n",
"└─────┴─────────────────┴───────────┴────────────┴────────┘\n",
"\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118346834631 -3.09 0.80 1.0 10.5ms\n",
" 2 -1.118346977882 -6.84 -4.11 0.80 1.0 22.2ms\n",
" 3 -1.118346979359 -8.83 -4.86 0.80 1.0 16.8ms\n",
" 4 -1.118346979384 -10.59 -5.78 0.80 1.0 19.4ms\n",
" 5 -1.118346979385 -12.08 -6.55 0.80 1.0 17.4ms\n",
" 6 -1.118346979385 -13.43 -7.04 0.80 1.0 19.8ms\n",
" 7 -1.118346979385 -14.75 -8.13 0.80 1.0 20.0ms\n",
" 8 -1.118346979385 -14.57 -8.95 0.80 1.0 20.3ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118354837055 -2.61 0.80 1.0 10.3ms\n",
" 2 -1.118356200096 -5.87 -3.62 0.80 1.0 18.5ms\n",
" 3 -1.118356214170 -7.85 -4.37 0.80 1.0 16.5ms\n",
" 4 -1.118356214417 -9.61 -5.29 0.80 1.0 74.0ms\n",
" 5 -1.118356214425 -11.10 -6.06 0.80 1.0 394ms\n",
" 6 -1.118356214425 -12.45 -6.55 0.80 1.0 17.8ms\n",
" 7 -1.118356214425 -13.53 -7.65 0.80 1.0 18.1ms\n",
" 8 -1.118356214425 + -14.75 -8.45 0.80 1.0 18.7ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118356214425 -9.24 0.80 1.0 12.0ms\n",
" 2 -1.118356214425 -14.65 -9.82 0.80 1.0 22.8ms\n",
"\n",
"┌─────┬─────────────────┬───────────┬────────────┬────────┐\n",
"│ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │\n",
"├─────┼─────────────────┼───────────┼────────────┼────────┤\n",
"│ 2 │ -1.118356214425 │ -4.79 │ 4.47263e-5 │ 949ms │\n",
"└─────┴─────────────────┴───────────┴────────────┴────────┘\n",
"\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118356215887 -4.94 0.80 1.0 11.6ms\n",
" 2 -1.118356215917 -10.54 -5.96 0.80 1.0 23.1ms\n",
" 3 -1.118356215917 -12.52 -6.70 0.80 1.0 19.8ms\n",
" 4 -1.118356215917 -14.31 -7.62 0.80 1.0 18.7ms\n",
" 5 -1.118356215917 -14.81 -8.39 0.80 1.0 22.3ms\n",
" 6 -1.118356215917 -15.05 -8.89 0.80 1.0 25.9ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118356217810 -4.40 0.80 1.0 10.4ms\n",
" 2 -1.118356218156 -9.46 -5.42 0.80 1.0 17.0ms\n",
" 3 -1.118356218159 -11.45 -6.17 0.80 1.0 16.7ms\n",
" 4 -1.118356218159 -13.21 -7.09 0.80 1.0 17.6ms\n",
" 5 -1.118356218159 -14.65 -7.86 0.80 1.0 17.6ms\n",
" 6 -1.118356218159 -14.75 -8.35 0.80 1.0 17.8ms\n",
" 7 -1.118356218159 + -15.65 -9.45 0.80 1.0 21.8ms\n",
"n Energy log10(ΔE) log10(Δρ) α Diag Δtime\n",
"--- --------------- --------- --------- ---- ---- ------\n",
" 1 -1.118356218159 -10.08 0.80 1.0 10.4ms\n",
" 2 -1.118356218159 + -15.35 -11.02 0.80 1.0 16.5ms\n",
"\n",
"┌─────┬─────────────────┬───────────┬────────────┬────────┐\n",
"│ n │ Energy │ log10(ΔE) │ max(Force) │ Δtime │\n",
"├─────┼─────────────────┼───────────┼────────────┼────────┤\n",
"│ 3 │ -1.118356218159 │ -8.43 │ 1.04318e-8 │ 452ms │\n",
"└─────┴─────────────────┴───────────┴────────────┴────────┘\n",
"\n"
]
}
],
"cell_type": "code",
"source": [
"using GeometryOptimization\n",
"results = minimize_energy!(h2_crude, calc; tol_forces=2e-6, verbosity=2)\n",
"nothing # hide"
],
"metadata": {},
"execution_count": 3
},
{
"cell_type": "markdown",
"source": [
"Structure after optimisation (note that the atom has wrapped around)"
],
"metadata": {}
},
{
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": "FlexibleSystem(H₂, periodicity = TTT):\n cell_vectors : [ 10 0 0;\n 0 10 0;\n 0 0 10]u\"a₀\"\n\n Atom(H, [-0.0431828, -7.14517e-11, -3.1471e-10]u\"a₀\")\n Atom(H, [ 1.44318, -6.98109e-11, -3.10227e-10]u\"a₀\")\n"
},
"metadata": {},
"execution_count": 4
}
],
"cell_type": "code",
"source": [
"results.system"
],
"metadata": {},
"execution_count": 4
},
{
"cell_type": "markdown",
"source": [
"Compute final bond length:"
],
"metadata": {}
},
{
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Optimal bond length: 1.4863655847396426 a₀\n"
]
}
],
"cell_type": "code",
"source": [
"rmin = norm(position(results.system[1]) - position(results.system[2]))\n",
"println(\"Optimal bond length: \", rmin)"
],
"metadata": {},
"execution_count": 5
},
{
"cell_type": "markdown",
"source": [
"Our results (1.486 Bohr) agrees with the\n",
"[equivalent tutorial from ABINIT](https://docs.abinit.org/tutorial/base1/)."
],
"metadata": {}
}
],
"nbformat_minor": 3,
"metadata": {
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.4"
},
"kernelspec": {
"name": "julia-1.11",
"display_name": "Julia 1.11.4",
"language": "julia"
}
},
"nbformat": 4
}