We solve the 1D Gross-Pitaevskii equation with a custom potential.
This is similar to Gross-Pitaevskii equation in 1D example
and we show how to define local potentials attached to atoms, which allows for
instance to compute forces.
The custom potential is actually already defined as ElementGaussian
in DFTK, and could
be used as is.
using DFTK
using LinearAlgebra
First, we define a new element which represents a nucleus generating a Gaussian potential.
struct CustomPotential <: DFTK.Element
α # Prefactor
L # Width of the Gaussian nucleus
end
Some default values
CustomPotential() = CustomPotential(1.0, 0.5);
We extend the two methods providing access to the real and Fourier representation of the potential to DFTK.
function DFTK.local_potential_real(el::CustomPotential, r::Real)
-el.α / (√(2π) * el.L) * exp(- (r / el.L)^2 / 2)
end
function DFTK.local_potential_fourier(el::CustomPotential, p::Real)
# = ∫ V(r) exp(-ix⋅p) dx
-el.α * exp(- (p * el.L)^2 / 2)
end
!!! tip "Gaussian potentials and DFTK"
DFTK already implements CustomPotential
in form of the DFTK.ElementGaussian
,
so this explicit re-implementation is only provided for demonstration purposes.
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]];
In this example, we want to generate two Gaussian potentials generated by two "nuclei" localized at positions $x_1$ and $x_2$, that are expressed in $[0,1)$ in fractional coordinates. $|x_1 - x_2|$ should be different from $0.5$ to break symmetry and get nonzero forces.
x1 = 0.2
x2 = 0.8
positions = [[x1, 0, 0], [x2, 0, 0]]
gauss = CustomPotential()
atoms = [gauss, gauss];
We setup a Gross-Pitaevskii model
C = 1.0
α = 2;
n_electrons = 1 # Increase this for fun
terms = [Kinetic(),
AtomicLocal(),
LocalNonlinearity(ρ -> C * ρ^α)]
model = Model(lattice, atoms, positions; n_electrons, terms,
spin_polarization=:spinless); # use "spinless electrons"
We discretize using a moderate Ecut and run a SCF algorithm to compute forces
afterwards. As there is no ionic charge associated to gauss
we have to specify
a starting density and we choose to start from a zero density.
basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))
ρ = zeros(eltype(basis), basis.fft_size..., 1)
scfres = self_consistent_field(basis; tol=1e-5, ρ)
scfres.energies
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -0.143642470793 -0.42 8.0 214ms 2 -0.156043550161 -1.91 -1.10 1.0 64.3ms 3 -0.156750350217 -3.15 -1.56 2.0 1.25ms 4 -0.156882954442 -3.88 -1.83 1.0 869μs 5 -0.156868671005 + -4.85 -1.82 1.0 847μs 6 -0.157056333936 -3.73 -3.58 1.0 801μs 7 -0.157056399754 -7.18 -3.80 1.0 799μs 8 -0.157056402506 -8.56 -3.95 1.0 800μs 9 -0.157056406788 -8.37 -4.52 1.0 799μs 10 -0.157056406871 -10.08 -4.83 1.0 833μs 11 -0.157056406905 -10.46 -5.07 1.0 823μs
Energy breakdown (in Ha): Kinetic 0.0380288 AtomicLocal -0.3163455 LocalNonlinearity 0.1212603 total -0.157056406905
Computing the forces can then be done as usual:
compute_forces(scfres)
2-element Vector{StaticArraysCore.SVector{3, Float64}}: [-0.05568361387364795, -0.0, -0.0] [0.0556832693863373, -0.0, -0.0]
Extract the converged total local potential
tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1
Extract other quantities before plotting them
ρ = scfres.ρ[:, 1, 1, 1] # converged density, first spin component
ψ_fourier = scfres.ψ[1][:, 1] # first k-point, all G components, first eigenvector
101-element Vector{ComplexF64}: -0.8400051425469964 + 0.49040839725185287im -0.08632740111373867 + 0.05040011297229624im 0.10228681179269153 - 0.05971593489446094im 0.044434076717601675 - 0.025941870869170266im -0.007754205911498586 + 0.0045263640560960685im -0.011799368017153937 + 0.006888746374955894im -0.001702150674911634 + 0.0009939156403862104im 0.0018938631629404924 - 0.001105654388096641im 0.0008156454450108033 - 0.00047624995440458345im -0.00012417151240156553 + 7.24831514772523e-5im ⋮ -0.0001241631204976044 + 7.249891411333487e-5im 0.0008157037705721767 - 0.00047615661631790035im 0.0018938541065647658 - 0.0011056728965533572im -0.0017022981011222556 + 0.0009936565698363657im -0.011799444812762731 + 0.006888621618192733im -0.007753641809124232 + 0.004527355650359131im 0.044434559038039965 - 0.02594105878631225im 0.10228615399162452 - 0.059717024773320394im -0.08632809528487465 + 0.050398918821539254im
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)]));
using Plots
x = a * vec(first.(DFTK.r_vectors(basis)))
p = plot(x, real.(ψ), label="real(ψ)")
plot!(p, x, imag.(ψ), label="imag(ψ)")
plot!(p, x, ρ, label="ρ")
plot!(p, x, tot_local_pot, label="tot local pot")