Support isolated systems

[epolack]

Supersede PR #650.

[antoine-levitt]

Test script (to debug):

diff --git a/src/elements.jl b/src/elements.jl
index 839ad131e..5949354ec 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -205,3 +205,12 @@ end
 function local_potential_fourier(el::ElementGaussian, q::Real)
     -el.α * exp(- (q * el.L)^2 / 2)  # = ∫_ℝ³ V(x) exp(-ix⋅q) dx
 end
+
+
+
+function charge_real(el::Element, r)
+    zero(r)
+ end
+function charge_fourier(el::Element, q)
+    zero(q)
+ end
diff --git a/src/terms/hartree.jl b/src/terms/hartree.jl
index b2a25759d..d66a11c74 100644
--- a/src/terms/hartree.jl
+++ b/src/terms/hartree.jl
@@ -80,7 +80,11 @@ end
                                             ψ, occupation; ρ, kwargs...) where {T}
     model = basis.model
     ρtot_real = total_density(ρ)
+    ρtot_real .-= [sum(DFTK.charge_real(el, norm(r - vector_red_to_cart(model, pos)))
+                       for (el, pos) in zip(model.atoms, model.positions))
+                   for r in r_vectors_cart(basis)]
     ρtot_fourier = fft(basis, ρtot_real)
+    # println(sum(ρtot_real) * basis.dvol)
     pot_fourier = term.poisson_green_coeffs .* ρtot_fourier
     pot_real = irfft(basis, pot_fourier)
 
@@ -113,6 +117,10 @@ end
         coeffs_ders = map(1:3) do α
             get_dipole(α, center, basis, ρtot_real) / get_dipole(α, center, basis, ρders[α])
         end
+        # println(coeffs_ders)
+        if basis.model.n_dim == 1
+            coeffs_ders[2:3] .= 0
+        end
         ρref = total_charge * ρrad + sum([coeffs_ders[α]*ρders[α] for α=1:3])
 
         # compute corresponding solution of -ΔVref = 4π ρref

using DFTK
using LinearAlgebra

x = 1.234
M = 2.345
α = 3.456
d = 1
ε = 1e-4
f(x) = DFTK.Vref_real(x, M, α, d)
fpp(x) = -(f(x+ε)+f(x-ε)-2f(x))/(ε^2) / (4π)
xs = range(0, 5, 100)
using PyPlot
plot(xs, fpp.(xs))
plot(xs, DFTK.ρref_real.(xs, M, α, d))
println(-(f(x+ε)+f(x-ε)-2f(x))/(ε^2) - 4π*DFTK.ρref_real(x, M, α, d))


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

    α_charge
    L_charge
end

# Some default values
CustomPotential() = CustomPotential(1.0, 0.5, 1.0, 0.5);

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)
    -el.α * exp(- (q * el.L)^2 / 2)
end

function DFTK.charge_real(el::CustomPotential, r::Real)
    el.α_charge / (√(2π) * el.L_charge) * exp(- (r / el.L_charge)^2 / 2)
end
function DFTK.charge_fourier(el::CustomPotential, q::Real)
    el.α_charge * exp(- (q * el.L_charge)^2 / 2)
end


a = 10
lattice = a .* [[1 0 0.]; [0 0 0]; [0 0 0]];

# x1 = 0.2
# x2 = 0.8
# positions = [[x1, 0, 0], [x2, 0, 0]]
gauss     = CustomPotential(1, 0.5, 1, .5)
positions = [[0.5, 0, 0]]
atoms     = [gauss];


gauss     = CustomPotential(1, 0.5, .333333, .5)
δ = 2
positions = [[1/2+δ/a, 0, 0],
             [1/2+δ/a, 0, 0],
             [1/2-δ/a, 0, 0]]
atoms     = [gauss, gauss, gauss];

# We setup a Gross-Pitaevskii model
n_electrons = 1  # Increase this for fun
terms = [Kinetic(),
         AtomicLocal(),
         Hartree(scaling_factor=.4),
         # LocalNonlinearity(ρ -> C * ρ^α)
]
per = false
periodic = fill(per, 3)
model = Model(lattice, atoms, positions; n_electrons, terms,
              spin_polarization=:spinless, periodic=periodic);  # 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))
ρ = ones(eltype(basis), basis.fft_size..., 1)
scfres = self_consistent_field(basis; tol=1e-8, ρ, maxiter=1000)
scfres.energies

# Extract the converged total local potential

# 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
ψ = ifft(basis, basis.kpoints[1], ψ_fourier)[:, 1, 1]
ψ /= (ψ[div(end, 2)] / abs(ψ[div(end, 2)]));

using PyPlot
x = vec(first.(DFTK.r_vectors_cart(basis)))
atomic_pot = (scfres.ham).blocks[1].operators[2].potential[:, 1, 1]; # use only dimension 1
hartree_pot = (scfres.ham).blocks[1].operators[3].potential[:, 1, 1]; # use only dimension 1
# plot(x, atomic_pot, label="Vat")
dash = per ? "--" : "-"
plot(x.-a/2, ρ, dash, label="ρ per=$per a=$a")
plot(x.-a/2, hartree_pot, dash, label="Vh per=$per a=$a")
xlim(-6, 6)

legend()