# Custom potential¶

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.

In :
using DFTK
using LinearAlgebra


First, we define a new element which represents a nucleus generating a Gaussian potential.

In :
struct CustomPotential <: DFTK.Element
α  # Prefactor
L  # Width of the Gaussian nucleus
end


Some default values

In :
CustomPotential() = CustomPotential(1.0, 0.5)

Out:
Main.var"##419".CustomPotential

We extend the two methods providing access to the real and Fourier representation of the potential to DFTK.

In :
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, q::Real)
# = ∫ V(r) exp(-ix⋅q) dx
-el.α * exp(- (q * el.L)^2 / 2)
end


We set up the lattice. For a 1D case we supply two zero lattice vectors

In :
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.

In :
x1 = 0.2
x2 = 0.8
positions = [[x1, 0, 0], [x2, 0, 0]]
gauss     = CustomPotential()
atoms     = [gauss, gauss]

Out:
2-element Vector{Main.var"##419".CustomPotential}:
Main.var"##419".CustomPotential(X)
Main.var"##419".CustomPotential(X)

We setup a Gross-Pitaevskii model

In :
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.

In :
basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))
ρ = zeros(eltype(basis), basis.fft_size..., 1)
scfres = self_consistent_field(basis; tol=1e-8, ρ=ρ)
scfres.energies

n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
1   -0.143581349776                   -0.42    8.0
2   -0.156036092192       -1.90       -1.10    1.0
3   -0.156768710397       -3.14       -1.56    2.0
4   -0.157045890563       -3.56       -2.32    2.0
5   -0.157055546819       -5.02       -2.99    2.0
6   -0.157056387622       -6.08       -3.71    1.0
7   -0.157056406098       -7.73       -4.37    1.0
8   -0.157056406912       -9.09       -5.40    2.0

Out:
Energy breakdown (in Ha):
Kinetic             0.0380295
AtomicLocal         -0.3163467
LocalNonlinearity   0.1212608

total               -0.157056406912

Computing the forces can then be done as usual:

In :
compute_forces(scfres)

Out:
2-element Vector{StaticArraysCore.SVector{3, Float64}}:
[-0.05568392329096156, 0.0, 0.0]
[0.05568743546713474, 0.0, 0.0]

Extract the converged total local potential

In :
tot_local_pot = DFTK.total_local_potential(scfres.ham)[:, 1, 1]; # use only dimension 1


Extract other quantities before plotting them

In :
ρ = scfres.ρ[:, 1, 1, 1]        # converged density, first spin component
ψ_fourier = scfres.ψ[:, 1]   # first k-point, all G components, first eigenvector
ψ = ifft(basis, basis.kpoints, ψ_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")

Out: