From 0a447fd6e5463bfecb7e6ae03fa842ed981bbd43 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Tue, 2 Jan 2024 15:43:04 +0100 Subject: [PATCH 1/6] Exact exchange Co-authored-by: bkaperick --- .github/workflows/CompatHelper.yaml | 1 + examples/silicon_hf.jl | 21 ++++++++ src/DFTK.jl | 3 +- src/standard_models.jl | 14 +++++- src/terms/exact_exchange.jl | 76 +++++++++++++++++++++++++++++ src/terms/operators.jl | 44 ++++++++++++++++- src/terms/terms.jl | 1 + test/exact_exchange.jl | 19 ++++++++ test/hamiltonian_consistency.jl | 23 ++++----- 9 files changed, 186 insertions(+), 16 deletions(-) create mode 100644 examples/silicon_hf.jl create mode 100644 src/terms/exact_exchange.jl create mode 100644 test/exact_exchange.jl diff --git a/.github/workflows/CompatHelper.yaml b/.github/workflows/CompatHelper.yaml index 06ec3fcc6d..0b7f2c98d7 100644 --- a/.github/workflows/CompatHelper.yaml +++ b/.github/workflows/CompatHelper.yaml @@ -6,6 +6,7 @@ on: jobs: CompatHelper: + if: github.repository_owner = "JuliaMolSim" runs-on: ubuntu-latest steps: - uses: julia-actions/setup-julia@latest diff --git a/examples/silicon_hf.jl b/examples/silicon_hf.jl new file mode 100644 index 0000000000..64b634b555 --- /dev/null +++ b/examples/silicon_hf.jl @@ -0,0 +1,21 @@ +# Very basic setup, useful for testing +using DFTK + +a = 10.26 # Silicon lattice constant in Bohr +lattice = a / 2 * [[0 1 1.]; + [1 0 1.]; + [1 1 0.]] +Si = ElementPsp(:Si; psp=load_psp("hgh/lda/Si-q4")) +atoms = [Si, Si] +positions = [ones(3)/8, -ones(3)/8] + +model = model_PBE(lattice, atoms, positions) +basis = PlaneWaveBasis(model; Ecut=30, kgrid=[1, 1, 1]) +scfres = self_consistent_field(basis; tol=1e-6) + +model = model_HF(lattice, atoms, positions) +basis = PlaneWaveBasis(model; basis.Ecut, basis.kgrid) +scfres = self_consistent_field(basis; + solver=DFTK.scf_damping_solver(1.0), + tol=1e-10, scfres.ρ, scfres.ψ, + diagtolalg=DFTK.AdaptiveDiagtol(; ratio_ρdiff=5e-4)) diff --git a/src/DFTK.jl b/src/DFTK.jl index c6fbc59e68..c64e506ab3 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -100,6 +100,7 @@ export Hamiltonian export HamiltonianBlock export energy_hamiltonian export Kinetic +export ExactExchange export ExternalFromFourier export ExternalFromReal export AtomicLocal @@ -143,7 +144,7 @@ include("eigen/preconditioners.jl") include("eigen/diag.jl") export model_atomic -export model_DFT, model_PBE, model_LDA, model_SCAN +export model_DFT, model_PBE, model_LDA, model_SCAN, model_HF include("standard_models.jl") export KerkerMixing, KerkerDosMixing, SimpleMixing, DielectricMixing diff --git a/src/standard_models.jl b/src/standard_models.jl index c8d1a69945..db0b7cb869 100644 --- a/src/standard_models.jl +++ b/src/standard_models.jl @@ -44,6 +44,18 @@ function model_DFT(lattice::AbstractMatrix, model_DFT(lattice, atoms, positions, Xc(functionals); kwargs...) end +""" +Build an Hartree-Fock model from the specified atoms. +""" +function model_HF(lattice::AbstractMatrix, atoms::Vector{<:Element}, + positions::Vector{<:AbstractVector}; + extra_terms=[], scaling_factor=1, kwargs...) + # TODO scaling_factor is a dirty trick for testing for now ... remove this argument later + @warn "Exact exchange in DFTK is hardly optimised and not yet production-ready." + model_atomic(lattice, atoms, positions; + extra_terms=[Hartree(), ExactExchange(scaling_factor), extra_terms...], + model_name="HF", kwargs...) +end """ Build an LDA model (Perdew & Wang parametrization) from the specified atoms. @@ -76,7 +88,7 @@ end # Generate equivalent functions for AtomsBase -for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN) +for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN, :model_HF) @eval function $fun(system::AbstractSystem, args...; kwargs...) parsed = parse_system(system) $fun(parsed.lattice, parsed.atoms, parsed.positions, args...; diff --git a/src/terms/exact_exchange.jl b/src/terms/exact_exchange.jl new file mode 100644 index 0000000000..9f17705c48 --- /dev/null +++ b/src/terms/exact_exchange.jl @@ -0,0 +1,76 @@ +""" +Exact exchange term: the Hartree-Exact exchange energy of the orbitals + +-1/2 ∑ ∫∫ ϕ_i^*(r)ϕ_j^*(r')ϕ_i(r')ϕ_j(r) / |r - r'| dr dr' + +""" +struct ExactExchange + scaling_factor::Real # to scale the term +end +ExactExchange(; scaling_factor=1) = ExactExchange(scaling_factor) +(exchange::ExactExchange)(basis) = TermExactExchange(basis, exchange.scaling_factor) +function Base.show(io::IO, exchange::ExactExchange) + fac = isone(exchange.scaling_factor) ? "" : "scaling_factor=$(exchange.scaling_factor)" + print(io, "ExactExchange($fac)") +end +struct TermExactExchange <: Term + scaling_factor::Real # scaling factor, absorbed into poisson_green_coeffs + poisson_green_coeffs::AbstractArray +end +function TermExactExchange(basis::PlaneWaveBasis{T}, scaling_factor) where T + poisson_green_coeffs = 4T(π) ./ [sum(abs2, basis.model.recip_lattice * G) + for G in to_cpu(G_vectors(basis))] + poisson_green_coeffs[1] = 0 + TermExactExchange(T(scaling_factor), scaling_factor .* poisson_green_coeffs) +end + +# Note: Implementing exact exchange in a scalable and numerically stable way, such that it +# rapidly converges with k-points is tricky. This implementation here is far too simple and +# slow to be useful. +# +# ABINIT/src/66_wfs/m_getghc.F90 +# ABINIT/src/66_wfs/m_fock_getghc.F90 +# ABINIT/src/66_wfs/m_prep_kgb.F90 +# ABINIT/src/66_wfs/m_bandfft_kpt.F90 +# +# For further information (in particular on regularising the Coulomb), consider the following +# https://www.vasp.at/wiki/index.php/Coulomb_singularity +# https://journals.aps.org/prb/pdf/10.1103/PhysRevB.34.4405 (QE default) +# https://journals.aps.org/prb/pdf/10.1103/PhysRevB.73.205119 +# https://journals.aps.org/prb/pdf/10.1103/PhysRevB.77.193110 +# https://docs.abinit.org/topics/documents/hybrids-2017.pdf (Abinit apparently +# uses a short-ranged Coulomb) + +@timing "ene_ops: ExactExchange" function ene_ops(term::TermExactExchange, + basis::PlaneWaveBasis{T}, ψ, occupation; + kwargs...) where {T} + if isnothing(ψ) || isnothing(occupation) + return (; E=T(0), ops=NoopOperator.(basis, basis.kpoints)) + end + + # TODO Occupation threshold + ψ, occupation = select_occupied_orbitals(basis, ψ, occupation; threshold=1e-8) + + E = zero(T) + + @assert length(ψ) == 1 # TODO: make it work for more kpoints + ik = 1 + kpt = basis.kpoints[ik] + occk = occupation[ik] + ψk = ψ[ik] + + for (n, ψn) in enumerate(eachcol(ψk)) + ψn_real = ifft(basis, kpt, ψn) + for (m, ψm) in enumerate(eachcol(ψk)) + ψm_real = ifft(basis, kpt, ψm) + ρnm_real = conj(ψn_real) .* ψm_real + ρnm_fourier = fft(basis, ρnm_real) + + fac_mn = occk[n] * occk[m] / T(2) + E -= fac_mn * real(dot(ρnm_fourier .* term.poisson_green_coeffs, ρnm_fourier)) + end + end + + ops = [ExchangeOperator(basis, kpt, term.poisson_green_coeffs, occk, ψk)] + (; E, ops) +end diff --git a/src/terms/operators.jl b/src/terms/operators.jl index 1602e0db26..941ce22a19 100644 --- a/src/terms/operators.jl +++ b/src/terms/operators.jl @@ -10,6 +10,22 @@ They also implement mul! and Matrix(op) for exploratory use. abstract type RealFourierOperator end # RealFourierOperator contain fields `basis` and `kpoint` +Base.Array(op::RealFourierOperator) = Matrix(op) + +# Slow fallback for getting operator as a dense matrix +function Matrix(op::RealFourierOperator) + T = complex(eltype(op.basis)) + n_G = length(G_vectors(op.basis, op.kpoint)) + H = zeros(T, n_G, n_G) + ψ = zeros(T, n_G) + @views for i in 1:n_G + ψ[i] = one(T) + mul!(H[:, i], op, ψ) + ψ[i] = zero(T) + end + H +end + # Unoptimized fallback, intended for exploratory use only. # For performance, call through Hamiltonian which saves FFTs. function LinearAlgebra.mul!(Hψ::AbstractVector, op::RealFourierOperator, ψ::AbstractVector) @@ -147,7 +163,6 @@ function apply!(Hψ, op::DivAgradOperator, ψ, ∂αψ_real = ifft!(ψ_scratch, op.basis, op.kpoint, im .* G_plus_k[α] .* ψ.fourier) A∇ψ = fft(op.basis, op.kpoint, ∂αψ_real .* op.A) Hψ.fourier .-= im .* G_plus_k[α] .* A∇ψ ./ 2 - end end # TODO Implement Matrix(op::DivAgrad) @@ -164,3 +179,30 @@ function optimize_operators_(ops) sum([op.potential for op in RSmults])) [nonRSmults..., combined_RSmults] end + +struct ExchangeOperator{T <: Real} <: RealFourierOperator + basis::PlaneWaveBasis{T} + kpoint::Kpoint{T} + poisson_green_coeffs # TODO This needs to be typed + occk # TODO This needs to be typed + ψk # TODO This needs to be typed +end +function apply!(Hψ, op::ExchangeOperator, ψ) + # Hψ = - ∑_n f_n ψ_n(r) ∫ (ψ_n)†(r') * ψ(r') / |r-r'| dr' + for (n, ψnk) in enumerate(eachcol(op.ψk)) + ψnk_real = ifft(op.basis, op.kpoint, ψnk) + x_real = conj(ψnk_real) .* ψ.real + + # TODO I think some symmetrisation is needed here since there are G-vectors, which are not balanced for even grids, which do not get properly converged. + + + # Compute integral by Poisson solve + x_four = fft(op.basis, x_real) + Vx_four = x_four .* op.poisson_green_coeffs + Vx_real = ifft(op.basis, Vx_four) + + # Real-space multiply and accumulate + Hψ.real .-= op.occk[n] .* ψnk_real .* Vx_real + end +end +# TODO Implement Matrix(op::ExchangeOperator) diff --git a/src/terms/terms.jl b/src/terms/terms.jl index a8bf5bf769..8e56e06110 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -51,6 +51,7 @@ include("ewald.jl") include("psp_correction.jl") include("entropy.jl") include("pairwise.jl") +include("exact_exchange.jl") include("magnetic.jl") breaks_symmetries(::Magnetic) = true diff --git a/test/exact_exchange.jl b/test/exact_exchange.jl new file mode 100644 index 0000000000..8907830bbe --- /dev/null +++ b/test/exact_exchange.jl @@ -0,0 +1,19 @@ +@testitem "Helium exchange energy" setup=[TestCases] begin + using DFTK + silicon = TestCases.silicon + + lattice = 10diagm(ones(3)) + positions = [0.5ones(3)] + atoms = [ElementCoulomb(:He)] + model = model_atomic(lattice, atoms, positions) + basis = PlaneWaveBasis(model; Ecut=15, kgrid=(1, 1, 1)) + scfres = self_consistent_field(basis) + + Eh, oph = DFTK.ene_ops(Hartree()(basis), basis, scfres.ψ, scfres.occupation; scfres...) + Ex, opx = DFTK.ene_ops(ExactExchange()(basis), basis, scfres.ψ, scfres.occupation; scfres...) + @test abs(Ex + Eh) < 1e-12 + + mat_h = real.(Matrix(only(oph))) + mat_x = real.(Matrix(only(opx))) + @show maximum(abs, mat_x - mat_x') +end diff --git a/test/hamiltonian_consistency.jl b/test/hamiltonian_consistency.jl index 11128503d8..e10f531a75 100644 --- a/test/hamiltonian_consistency.jl +++ b/test/hamiltonian_consistency.jl @@ -9,21 +9,14 @@ testcase = silicon function test_matrix_repr_operator(hamk, ψk; atol=1e-8) for operator in hamk.operators - try - operator_matrix = Matrix(operator) - @test norm(operator_matrix*ψk - operator*ψk) < atol - catch e - allowed_missing_operators = Union{DFTK.DivAgradOperator, - DFTK.MagneticFieldOperator} - @assert operator isa allowed_missing_operators - @info "Matrix of operator $(nameof(typeof(operator))) not yet supported" maxlog=1 - end + operator_matrix = Matrix(operator) + @test norm(operator_matrix*ψk - operator*ψk) < atol end end function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, 3], - kshift=[0, 1, 0]/2, lattice=testcase.lattice, Ecut=10, - spin_polarization=:none) + kshift=[0, 1, 0]/2, lattice=testcase.lattice, + Ecut=10, integer_occupation=false, spin_polarization=:none) sspol = spin_polarization != :none ? " ($spin_polarization)" : "" xc = term isa Xc ? "($(first(term.functionals)))" : "" @testset "$(typeof(term))$xc $sspol" begin @@ -43,6 +36,9 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, occupation = [filled_occ * rand(n_bands) for _ = 1:length(basis.kpoints)] occ_scaling = n_electrons / sum(sum(occupation)) occupation = [occ * occ_scaling for occ in occupation] + if integer_occupation + occupation = [filled_occ * ones(n_bands) for _ in 1:length(basis.kpoints)] + end ρ = with_logger(NullLogger()) do compute_density(basis, ψ, occupation) end @@ -64,10 +60,10 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, diff_predicted = 0.0 for ik = 1:length(basis.kpoints) - Hψk = ham.blocks[ik]*ψ[ik] + Hψk = ham.blocks[ik] * ψ[ik] test_matrix_repr_operator(ham.blocks[ik], ψ[ik]; atol) δψkHψk = sum(occupation[ik][iband] * real(dot(δψ[ik][:, iband], Hψk[:, iband])) - for iband=1:n_bands) + for iband = 1:n_bands) diff_predicted += 2 * basis.kweights[ik] * δψkHψk end diff_predicted = mpi_sum(diff_predicted, basis.comm_kpts) @@ -101,6 +97,7 @@ end test_consistency_term(Xc(:mgga_c_scan), spin_polarization=:collinear) test_consistency_term(Xc(:mgga_x_b00)) test_consistency_term(Xc(:mgga_c_b94), spin_polarization=:collinear) + test_consistency_term(ExactExchange(); kgrid=(1, 1, 1), Ecut=7) let a = 6 From 0c26b76669bb02b69c496466f839d85937a40600 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Thu, 4 Jan 2024 10:06:01 +0100 Subject: [PATCH 2/6] Hybrids --- Project.toml | 2 +- examples/{silicon_hf.jl => silicon_exx.jl} | 12 ++++++++++- src/DFTK.jl | 2 +- src/DispatchFunctional.jl | 18 +++++++++++++--- src/standard_models.jl | 25 ++++++++++++++++++++-- src/terms/exact_exchange.jl | 8 +++---- src/terms/operators.jl | 8 +++---- test/hamiltonian_consistency.jl | 23 ++++++++------------ 8 files changed, 68 insertions(+), 30 deletions(-) rename examples/{silicon_hf.jl => silicon_exx.jl} (56%) diff --git a/Project.toml b/Project.toml index 9d87d82ab8..fceea3c7e5 100644 --- a/Project.toml +++ b/Project.toml @@ -89,7 +89,7 @@ IterativeSolvers = "0.9" JLD2 = "0.4" JSON3 = "1" LazyArtifacts = "1.3" -Libxc = "0.3.17" +Libxc = "0.3.18" LineSearches = "7" LinearAlgebra = "1" LinearMaps = "3" diff --git a/examples/silicon_hf.jl b/examples/silicon_exx.jl similarity index 56% rename from examples/silicon_hf.jl rename to examples/silicon_exx.jl index 64b634b555..6bbcf3419f 100644 --- a/examples/silicon_hf.jl +++ b/examples/silicon_exx.jl @@ -13,9 +13,19 @@ model = model_PBE(lattice, atoms, positions) basis = PlaneWaveBasis(model; Ecut=30, kgrid=[1, 1, 1]) scfres = self_consistent_field(basis; tol=1e-6) +# Run hybrid-DFT tuning some DFTK defaults +# (Anderson does not work well right now as orbitals not taken into account) +model = model_PBE0(lattice, atoms, positions) +basis = PlaneWaveBasis(model; basis.Ecut, basis.kgrid) +scfres = self_consistent_field(basis; + solver=DFTK.scf_damping_solver(1.0), + tol=1e-8, scfres.ρ, scfres.ψ, + diagtolalg=DFTK.AdaptiveDiagtol(; ratio_ρdiff=5e-4)) + +# Run Hartree-Fock model = model_HF(lattice, atoms, positions) basis = PlaneWaveBasis(model; basis.Ecut, basis.kgrid) scfres = self_consistent_field(basis; solver=DFTK.scf_damping_solver(1.0), - tol=1e-10, scfres.ρ, scfres.ψ, + tol=1e-8, scfres.ρ, scfres.ψ, diagtolalg=DFTK.AdaptiveDiagtol(; ratio_ρdiff=5e-4)) diff --git a/src/DFTK.jl b/src/DFTK.jl index c64e506ab3..5d8ec42468 100644 --- a/src/DFTK.jl +++ b/src/DFTK.jl @@ -144,7 +144,7 @@ include("eigen/preconditioners.jl") include("eigen/diag.jl") export model_atomic -export model_DFT, model_PBE, model_LDA, model_SCAN, model_HF +export model_DFT, model_PBE, model_LDA, model_SCAN, model_PBE0, model_HF include("standard_models.jl") export KerkerMixing, KerkerDosMixing, SimpleMixing, DielectricMixing diff --git a/src/DispatchFunctional.jl b/src/DispatchFunctional.jl index 20373c191b..29b7ff4016 100644 --- a/src/DispatchFunctional.jl +++ b/src/DispatchFunctional.jl @@ -12,12 +12,17 @@ function LibxcFunctional(identifier::Symbol) @assert fun.kind in (:exchange, :correlation, :exchange_correlation) kind = Dict(:exchange => :x, :correlation => :c, :exchange_correlation => :xc)[fun.kind] - @assert fun.family in (:lda, :gga, :mgga) # Hybrids not supported yet. - if fun.family == :mgga && Libxc.needs_laplacian(fun) - family = :mggal + @assert fun.family in (:lda, :gga, :mgga, :hyb_lda, :hyb_gga, :hyb_mgga) + # Libxc maintains the distinction between hybrid and non-hybrid equivalents, + # but this distinction is not relevant for the functional interface + if startswith(string(fun.family), "hyb") + family = Symbol(string(fun.family)[5:end]) else family = fun.family end + if family == :mgga && Libxc.needs_laplacian(fun) + family = :mggal + end LibxcFunctional{family,kind}(identifier) end @@ -154,3 +159,10 @@ for fun in (:potential_terms, :kernel_terms) end end end + +# TODO This is hackish for now until Libxc has fully picked up the DftFunctionals.jl interface +exx_coefficient(::Functional{:lda}) = nothing +exx_coefficient(::Functional{:gga}) = nothing +exx_coefficient(::Functional{:mgga}) = nothing +exx_coefficient(fun::DispatchFunctional) = exx_coefficient(fun.inner) +exx_coefficient(fun::LibxcFunctional) = Libxc.Functional(fun.identifier).exx_coefficient diff --git a/src/standard_models.jl b/src/standard_models.jl index db0b7cb869..0b23767939 100644 --- a/src/standard_models.jl +++ b/src/standard_models.jl @@ -33,8 +33,14 @@ function model_DFT(lattice::AbstractMatrix, xc::Xc; extra_terms=[], kwargs...) model_name = isempty(xc.functionals) ? "rHF" : join(string.(xc.functionals), "+") + exx = [ExactExchange(; scaling_factor=exx_coefficient(f)) + for f in xc.functionals if !isnothing(exx_coefficient(f))] + if !isempty(exx) + @assert length(exx) == 1 + @warn "Exact exchange in DFTK is hardly optimised and not yet production-ready." + end model_atomic(lattice, atoms, positions; - extra_terms=[Hartree(), xc, extra_terms...], model_name, kwargs...) + extra_terms=[Hartree(), xc, exx..., extra_terms...], model_name, kwargs...) end function model_DFT(lattice::AbstractMatrix, atoms::Vector{<:Element}, @@ -57,6 +63,10 @@ function model_HF(lattice::AbstractMatrix, atoms::Vector{<:Element}, model_name="HF", kwargs...) end +# +# DFT shortcut models +# + """ Build an LDA model (Perdew & Wang parametrization) from the specified atoms. DOI:10.1103/PhysRevB.45.13244 @@ -87,8 +97,19 @@ function model_SCAN(lattice::AbstractMatrix, atoms::Vector{<:Element}, end +""" +Build an PBE0 model from the specified atoms. +DOI:10.1103/PhysRevLett.77.3865 +""" +function model_PBE0(lattice::AbstractMatrix, atoms::Vector{<:Element}, + positions::Vector{<:AbstractVector}; kwargs...) + model_DFT(lattice, atoms, positions, [:hyb_gga_xc_pbeh]; kwargs...) +end + + # Generate equivalent functions for AtomsBase -for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN, :model_HF) +for fun in (:model_atomic, :model_DFT, :model_HF, + :model_LDA, :model_PBE, :model_SCAN, :model_PBE0) @eval function $fun(system::AbstractSystem, args...; kwargs...) parsed = parse_system(system) $fun(parsed.lattice, parsed.atoms, parsed.positions, args...; diff --git a/src/terms/exact_exchange.jl b/src/terms/exact_exchange.jl index 9f17705c48..be25ff35f0 100644 --- a/src/terms/exact_exchange.jl +++ b/src/terms/exact_exchange.jl @@ -1,8 +1,8 @@ -""" +@doc raw""" Exact exchange term: the Hartree-Exact exchange energy of the orbitals - --1/2 ∑ ∫∫ ϕ_i^*(r)ϕ_j^*(r')ϕ_i(r')ϕ_j(r) / |r - r'| dr dr' - +```math +-1/2 ∑_{nm} f_n f_m ∫∫ \frac{ψ_n^*(r)ψ_m^*(r')ψ_n(r')ψ_m(r)}{|r - r'|} dr dr' +``` """ struct ExactExchange scaling_factor::Real # to scale the term diff --git a/src/terms/operators.jl b/src/terms/operators.jl index 941ce22a19..afeaf843f9 100644 --- a/src/terms/operators.jl +++ b/src/terms/operators.jl @@ -180,12 +180,12 @@ function optimize_operators_(ops) [nonRSmults..., combined_RSmults] end -struct ExchangeOperator{T <: Real} <: RealFourierOperator +struct ExchangeOperator{T <: Real,Tocc,Tpsi} <: RealFourierOperator basis::PlaneWaveBasis{T} kpoint::Kpoint{T} - poisson_green_coeffs # TODO This needs to be typed - occk # TODO This needs to be typed - ψk # TODO This needs to be typed + poisson_green_coeffs::Array{T, 3} + occk::Tocc + ψk::Tpsi end function apply!(Hψ, op::ExchangeOperator, ψ) # Hψ = - ∑_n f_n ψ_n(r) ∫ (ψ_n)†(r') * ψ(r') / |r-r'| dr' diff --git a/test/hamiltonian_consistency.jl b/test/hamiltonian_consistency.jl index e10f531a75..a91d73055b 100644 --- a/test/hamiltonian_consistency.jl +++ b/test/hamiltonian_consistency.jl @@ -7,16 +7,9 @@ using LinearAlgebra using ..TestCases: silicon testcase = silicon -function test_matrix_repr_operator(hamk, ψk; atol=1e-8) - for operator in hamk.operators - operator_matrix = Matrix(operator) - @test norm(operator_matrix*ψk - operator*ψk) < atol - end -end - function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, 3], kshift=[0, 1, 0]/2, lattice=testcase.lattice, - Ecut=10, integer_occupation=false, spin_polarization=:none) + Ecut=10, spin_polarization=:none) sspol = spin_polarization != :none ? " ($spin_polarization)" : "" xc = term isa Xc ? "($(first(term.functionals)))" : "" @testset "$(typeof(term))$xc $sspol" begin @@ -26,6 +19,7 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, model = Model(lattice, atoms, testcase.positions; terms=[term], spin_polarization, symmetries=true) basis = PlaneWaveBasis(model; Ecut, kgrid, kshift) + @assert length(basis.terms) == 1 n_electrons = testcase.n_electrons n_bands = div(n_electrons, 2, RoundUp) @@ -36,16 +30,19 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, occupation = [filled_occ * rand(n_bands) for _ = 1:length(basis.kpoints)] occ_scaling = n_electrons / sum(sum(occupation)) occupation = [occ * occ_scaling for occ in occupation] - if integer_occupation - occupation = [filled_occ * ones(n_bands) for _ in 1:length(basis.kpoints)] - end ρ = with_logger(NullLogger()) do compute_density(basis, ψ, occupation) end E0, ham = energy_hamiltonian(basis, ψ, occupation; ρ) - @assert length(basis.terms) == 1 + # Test operator is derivative of the energy + for ik = 1:length(basis.kpoints) + for operator in ham.blocks[ik].operators + @test norm(Matrix(operator) * ψ[ik] - operator * ψ[ik]) < atol + end + end + # Test operator is derivative of the energy δψ = [randn(ComplexF64, size(ψ[ik])) for ik = 1:length(basis.kpoints)] function compute_E(ε) ψ_trial = ψ .+ ε .* δψ @@ -55,13 +52,11 @@ function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, E = energy_hamiltonian(basis, ψ_trial, occupation; ρ=ρ_trial).energies E.total end - diff = (compute_E(ε) - compute_E(-ε)) / (2ε) diff_predicted = 0.0 for ik = 1:length(basis.kpoints) Hψk = ham.blocks[ik] * ψ[ik] - test_matrix_repr_operator(ham.blocks[ik], ψ[ik]; atol) δψkHψk = sum(occupation[ik][iband] * real(dot(δψ[ik][:, iband], Hψk[:, iband])) for iband = 1:n_bands) diff_predicted += 2 * basis.kweights[ik] * δψkHψk From 79629f21f1bcb228b470667636e82a1523b2dddd Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Thu, 4 Jan 2024 13:44:47 +0100 Subject: [PATCH 3/6] Update silicon_exx.jl --- examples/silicon_exx.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/silicon_exx.jl b/examples/silicon_exx.jl index 6bbcf3419f..672b2c67b7 100644 --- a/examples/silicon_exx.jl +++ b/examples/silicon_exx.jl @@ -10,7 +10,7 @@ atoms = [Si, Si] positions = [ones(3)/8, -ones(3)/8] model = model_PBE(lattice, atoms, positions) -basis = PlaneWaveBasis(model; Ecut=30, kgrid=[1, 1, 1]) +basis = PlaneWaveBasis(model; Ecut=20, kgrid=[1, 1, 1]) scfres = self_consistent_field(basis; tol=1e-6) # Run hybrid-DFT tuning some DFTK defaults From f71042e5830e505c10dab60cf20715af92daa817 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Thu, 4 Jan 2024 13:46:06 +0100 Subject: [PATCH 4/6] Update standard_models.jl --- src/standard_models.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/standard_models.jl b/src/standard_models.jl index 0b23767939..d9a677fb87 100644 --- a/src/standard_models.jl +++ b/src/standard_models.jl @@ -54,12 +54,10 @@ end Build an Hartree-Fock model from the specified atoms. """ function model_HF(lattice::AbstractMatrix, atoms::Vector{<:Element}, - positions::Vector{<:AbstractVector}; - extra_terms=[], scaling_factor=1, kwargs...) - # TODO scaling_factor is a dirty trick for testing for now ... remove this argument later + positions::Vector{<:AbstractVector}; extra_terms=[], kwargs...) @warn "Exact exchange in DFTK is hardly optimised and not yet production-ready." model_atomic(lattice, atoms, positions; - extra_terms=[Hartree(), ExactExchange(scaling_factor), extra_terms...], + extra_terms=[Hartree(), ExactExchange(), extra_terms...], model_name="HF", kwargs...) end From b1948b80df19e7346a22c45e7ed29131d062b4e8 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Thu, 4 Jan 2024 13:49:11 +0100 Subject: [PATCH 5/6] Update operators.jl --- src/terms/operators.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/terms/operators.jl b/src/terms/operators.jl index afeaf843f9..389f00607b 100644 --- a/src/terms/operators.jl +++ b/src/terms/operators.jl @@ -192,9 +192,7 @@ function apply!(Hψ, op::ExchangeOperator, ψ) for (n, ψnk) in enumerate(eachcol(op.ψk)) ψnk_real = ifft(op.basis, op.kpoint, ψnk) x_real = conj(ψnk_real) .* ψ.real - - # TODO I think some symmetrisation is needed here since there are G-vectors, which are not balanced for even grids, which do not get properly converged. - + # TODO Some symmetrisation of x_real might be needed here ... # Compute integral by Poisson solve x_four = fft(op.basis, x_real) From e8bee2111460da79f3d473c704348ca2fb2e348e Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Thu, 4 Jan 2024 14:28:09 +0100 Subject: [PATCH 6/6] Fix --- test/exact_exchange.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/exact_exchange.jl b/test/exact_exchange.jl index 8907830bbe..5945e7e039 100644 --- a/test/exact_exchange.jl +++ b/test/exact_exchange.jl @@ -1,5 +1,6 @@ @testitem "Helium exchange energy" setup=[TestCases] begin using DFTK + using LinearAlgebra silicon = TestCases.silicon lattice = 10diagm(ones(3))