Bug Using Charged systems

[umbriquse]

Found a bug relating to an assertion of charges in energy_ewald

System Version

(Dino) pkg> st DFTK
     Status `~/git/Dino.jl/Project.toml`
 [acf6eb54] DFTK v0.3.7

Code

using DFTK
using Unitful, UnitfulAtomic

Si = ElementPsp(:Si, psp = load_psp("hgh/lda/Si-q4"))
atoms = [
    Si => [ones(3)/8, -ones(3)/8]
]

lattice = austrip(5.431u"angstrom")/2 * [
    0.0 1.0 1.0;
    1.0 0.0 1.0;
    1.0 1.0 0.0
]
Ecut = 7

nElectrons = sum(n_elec_valence(element) * length(positions) for (element, positions) in atoms) - 2

model = model_LDA(lattice, atoms; n_electrons = nElectrons)

basis = PlaneWaveBasis(model, Ecut; kgrid = kgrid_size_from_minimal_spacing(lattice, π/10))

ERROR Message

ERROR: LoadError: AssertionError: sum(charges) == model.n_electrons
Stacktrace:
  [1] energy_ewald(model::Model{Float64}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DFTK ~/.julia/packages/DFTK/ulqQF/src/terms/ewald.jl:54
  [2] energy_ewald
    @ ~/.julia/packages/DFTK/ulqQF/src/terms/ewald.jl:47 [inlined]
  [3] TermEwald
    @ ~/.julia/packages/DFTK/ulqQF/src/terms/ewald.jl:17 [inlined]
  [4] (::Ewald)(basis::PlaneWaveBasis{Float64})
    @ DFTK ~/.julia/packages/DFTK/ulqQF/src/terms/ewald.jl:9
  [5] macro expansion
    @ ~/.julia/packages/TimerOutputs/ZQ0rt/src/TimerOutput.jl:220 [inlined]
  [6] macro expansion
    @ ~/.julia/packages/DFTK/ulqQF/src/PlaneWaveBasis.jl:243 [inlined]
  [7] PlaneWaveBasis(model::Model{Float64}, Ecut::Int64, kcoords::Vector{StaticArrays.SVector{3, Rational{Int64}}}, ksymops::Vector{Vector{Tuple{StaticArrays.SMatrix{3, 3, Int64, 9}, StaticArrays.SVector{3, Float64}}}}, symmetries::Vector{Tuple{StaticArrays.SMatrix{3, 3, Int64, 9}, StaticArrays.SVector{3, Float64}}}; fft_size::Nothing, variational::Bool, optimize_fft_size::Bool, supersampling::Int64, kgrid::Vector{Int64}, kshift::Vector{Float64}, comm_kpts::MPI.Comm)
    @ DFTK ~/.julia/packages/TimerOutputs/ZQ0rt/src/TimerOutput.jl:220
  [8] PlaneWaveBasis(model::Model{Float64}, Ecut::Int64; kgrid::Vector{Int64}, kshift::Vector{Float64}, use_symmetry::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DFTK ~/.julia/packages/DFTK/ulqQF/src/PlaneWaveBasis.jl:283
  [9] top-level scope
    @ ~/git/Dino.jl/test/issue.jl:21

Version Info

Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Apple M1
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, westmere)

[antoine-levitt]

We have not thought at all about charged systems; the functionality to have a separate number of electrons is only there for "artificial" models (eg without atoms). You're packing together positive charges, so the energy is infinite here, no?

[umbriquse]

True, but I have run self_consistent_field with a charged system in the past.

[antoine-levitt]

If it yielded a finite answer that was probably a bug... What's your use case here?

[umbriquse]

The above example was just to show the ERROR. Generally the use case is an atom bonding to a surface in a charged system.

EDIT: Spelling mistake

[antoine-levitt]

We can certainly disable the error, but I'm more interested in doing the right thing here. Don't you get an infinite energy in this case also?

[umbriquse]

I agree that disabling the error wouldn't be the right thing. We did not get infinite energy in our systems. The systems gave some numbers that cardinally agreed with the literature.

[antoine-levitt]

Physically, a 2D array of charges does not blow up? I've been meaning to look at these things for a while... The first step would be to be right for an isolated ion.

[umbriquse]

I'm not sure what you mean by your question. By blow up do you mean produce an infinite energy?

[antoine-levitt]

Not only energy, potential also, because it's 1/r and there are r of them. Maybe forces are conditionally summable, though...

[antoine-levitt]

Ie physically I'd expect the ions to start moving away from each other, pulling the surface apart?

[umbriquse]

I believe the forces cancelled in the repeating lattice. Also, the surface could've been neutrally charged and the adsorbed atom lost all of the electrons.

[umbriquse]

Regardless, thank you for the quick responses. This issue is resolved, but I'm curious what would need to be done in order to use charged systems?

[antoine-levitt]

That really depends on the system. Isolated charged system (eg an ion on a molecule) is rather easy, one just needs a correction scheme for the potential (I've never seen this done in the literature but it's not hard), but we need to implement it. Charged defects in semiconductors is more involved, there's a bunch of schemes, I need to look into it more. Charged ion on top of a surface I don't know. I have a PhD student starting who's going to look among other things into this kind of issues, so we'll see.

[mfherbst]

But in @antoine-levitt's comments is clearly the mathematician speaking. Many codes support charged unit cells and I don't think they do fancy corrections by default. I tested this with QE at some point and the difference between our definitely inconsistent implementation and QE (without corrections) was small, similar to what @umbriquse has pointed out. I think to a first approximation if we ensure that the compensating background is consistently taken into account in all terms, we can disable the assertion. We just have to communicate to the user that the compensating background is implicit ... and therefore that it's not a charged ion in vacuum you are modelling, but a charged ion with compensating uniform background charge.

[antoine-levitt]

I don't know what QE does but VASP has correction schemes that appear to be enabled by default: https://www.vasp.at/wiki/index.php/Monopole_Dipole_and_Quadrupole_corrections . I'm fine with doing something very naive here if we have a strong warning

[mfherbst]

Same. And getting the compensating backgrounds to agree is not hard, one just has to go through the relevant terms and be careful.

[mfherbst]

For QE you can switch on Markov-Payne, but you have to do so explicitly.

[antoine-levitt]

I thought Makov-Payne was for defects, but maybe it's also for isolated molecules. Makov-Payne is a bit of a hack : it only corrects the total energy, not the density or the potential. There's simpler and better ways to do that for isolated systems: you have to add and subtract a reference density that compensates the extra charge, and compute the potential for that analytically. I do that eg in the anyon code.