In this example we will use DFTK to solve the Gross-Pitaevskii equation, and use this opportunity to explore a few internals.
The Gross-Pitaevskii equation (GPE) is a simple non-linear equation used to model bosonic systems in a mean-field approach. Denoting by ψ the effective one-particle bosonic wave function, the time-independent GPE reads in atomic units: Hψ=(−12Δ+V+2C|ψ|2)ψ=μψ‖ψ‖L2=1
We wish to model this equation in 1D using DFTK. First we set up the lattice. For a 1D case we supply two zero lattice vectors,
a = 10
lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];
which is special cased in DFTK to support 1D models.
For the potential term V we pick a harmonic
potential. We use the function ExternalFromReal which uses
Cartesian coordinates ( see Lattices and lattice vectors).
pot(x) = (x - a/2)^2;
We setup each energy term in sequence: kinetic, potential and nonlinear term.
For the non-linearity we use the LocalNonlinearity(f) term of DFTK, with f(ρ) = C ρ^α.
This object introduces an energy term C∫ρ(r)αdr
to the total energy functional, thus a potential term αCρα−1.
In our case we thus need the parameters
C = 1.0
α = 2;
… and with this build the model
using DFTK
using LinearAlgebra
n_electrons = 1 # Increase this for fun
terms = [Kinetic(),
ExternalFromReal(r -> pot(r[1])),
LocalNonlinearity(ρ -> C * ρ^α),
]
model = Model(lattice; n_electrons, terms, spin_polarization=:spinless); # spinless electrons
We discretize using a moderate Ecut (For 1D values up to 5000 are completely fine)
and run a direct minimization algorithm:
basis = PlaneWaveBasis(model, Ecut=500, kgrid=(1, 1, 1))
scfres = direct_minimization(basis; tol=1e-8) # This is a constrained preconditioned LBFGS
scfres.energies
n Energy log10(ΔE) log10(Δρ) Δtime --- --------------- --------- --------- ------ 1 +169.2606726480 -1.52 361ms 2 +126.0086591069 1.64 -0.82 1.50ms 3 +72.03783082690 1.73 -0.54 1.53ms 4 +13.64379652598 1.77 -0.36 1.50ms 5 +11.28133065189 0.37 -0.31 996μs 6 +11.05645401970 -0.65 -0.64 829μs 7 +10.55002381403 -0.30 -1.06 500μs 8 +9.450164406687 0.04 -0.54 818μs 9 +7.354132105128 0.32 -0.37 1.00ms 10 +4.731470225247 0.42 -0.21 996μs 11 +1.741687359323 0.48 -0.04 983μs 12 +1.492955275769 -0.60 -0.22 817μs 13 +1.323463775139 -0.77 -0.44 673μs 14 +1.269640343253 -1.27 -0.54 664μs 15 +1.227237297763 -1.37 -0.54 637μs 16 +1.176962610400 -1.30 -0.96 654μs 17 +1.150536737251 -1.58 -1.16 666μs 18 +1.146711574614 -2.42 -1.51 633μs 19 +1.145380388805 -2.88 -1.53 654μs 20 +1.144422364102 -3.02 -1.86 658μs 21 +1.144251253364 -3.77 -2.11 636μs 22 +1.144099904777 -3.82 -2.20 660μs 23 +1.144080967071 -4.72 -2.59 676μs 24 +1.144053337500 -4.56 -2.74 677μs 25 +1.144047865766 -5.26 -2.91 699μs 26 +1.144042505909 -5.27 -2.43 679μs 27 +1.144038676988 -5.42 -2.65 639μs 28 +1.144037765972 -6.04 -3.63 664μs 29 +1.144037387767 -6.42 -4.14 483μs 30 +1.144037052483 -6.47 -4.24 499μs 31 +1.144036955532 -7.01 -4.37 658μs 32 +1.144036882112 -7.13 -4.22 640μs 33 +1.144036868946 -7.88 -4.33 661μs 34 +1.144036862487 -8.19 -4.49 639μs 35 +1.144036857093 -8.27 -4.54 679μs 36 +1.144036854705 -8.62 -4.82 664μs 37 +1.144036854215 -9.31 -4.56 635μs 38 +1.144036853822 -9.41 -5.13 662μs 39 +1.144036853032 -9.10 -4.92 660μs 40 +1.144036852876 -9.81 -5.49 641μs 41 +1.144036852797 -10.10 -5.41 653μs 42 +1.144036852774 -10.63 -6.41 487μs 43 +1.144036852766 -11.10 -6.21 636μs 44 +1.144036852759 -11.14 -6.10 660μs 45 +1.144036852757 -11.68 -6.01 657μs 46 +1.144036852756 -12.39 -6.35 638μs 47 +1.144036852756 -12.31 -6.54 651μs 48 +1.144036852755 -12.76 -7.05 638μs 49 +1.144036852755 -13.34 -6.78 656μs 50 +1.144036852755 -13.42 -7.27 664μs 51 +1.144036852755 -13.45 -7.27 637μs 52 +1.144036852755 -14.18 -7.48 660μs 53 +1.144036852755 -14.29 -7.67 667μs 54 +1.144036852755 -15.35 -8.02 627μs
Energy breakdown (in Ha):
Kinetic 0.2682057
ExternalFromReal 0.4707475
LocalNonlinearity 0.4050836
total 1.144036852755
We use the opportunity to explore some of DFTK internals.
Extract the converged density and the obtained wave function:
ρ = real(scfres.ρ)[:, 1, 1, 1] # converged density, first spin component
ψ_fourier = scfres.ψ[1][:, 1]; # first k-point, all G components, first eigenvector
Transform the wave function to real space and fix the phase:
ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));
Check whether ψ is normalised:
x = a * vec(first.(DFTK.r_vectors(basis)))
N = length(x)
dx = a / N # real-space grid spacing
@assert sum(abs2.(ψ)) * dx ≈ 1.0
The density is simply built from ψ:
norm(scfres.ρ - abs2.(ψ))
8.31875022839291e-16
We summarize the ground state in a nice plot:
using Plots
p = plot(x, real.(ψ), label="real(ψ)")
plot!(p, x, imag.(ψ), label="imag(ψ)")
plot!(p, x, ρ, label="ρ")
The energy_hamiltonian function can be used to get the energy and
effective Hamiltonian (derivative of the energy with respect to the density matrix)
of a particular state (ψ, occupation).
The density ρ associated to this state is precomputed
and passed to the routine as an optimization.
E, ham = energy_hamiltonian(basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ)
@assert E.total == scfres.energies.total
Now the Hamiltonian contains all the blocks corresponding to k-points. Here, we just have one k-point:
H = ham.blocks[1];
H can be used as a linear operator (efficiently using FFTs), or converted to a dense matrix:
ψ11 = scfres.ψ[1][:, 1] # first k-point, first eigenvector
Hmat = Array(H) # This is now just a plain Julia matrix,
# which we can compute and store in this simple 1D example
@assert norm(Hmat * ψ11 - H * ψ11) < 1e-10
Let's check that ψ11 is indeed an eigenstate:
norm(H * ψ11 - dot(ψ11, H * ψ11) * ψ11)
1.0563325319244079e-7
Build a finite-differences version of the GPE operator H, as a sanity check:
A = Array(Tridiagonal(-ones(N - 1), 2ones(N), -ones(N - 1)))
A[1, end] = A[end, 1] = -1
K = A / dx^2 / 2
V = Diagonal(pot.(x) + C .* α .* (ρ.^(α-1)))
H_findiff = K + V;
maximum(abs.(H_findiff*ψ - (dot(ψ, H_findiff*ψ) / dot(ψ, ψ)) * ψ))
0.00022340550154961012