Add Hubbard corrections and DFT+U #1158

fsicignano · 2025-09-25T13:03:36Z

New Hubbard U term implementation.

Key features: - New hubbard.jl program
- New operator added in operators.jl
- OrbitalManifold struct has been moved to hubbard.jl
- self_consistent_field.jl has been updated to store the Hubbard n matrix and the occupation
- Hubbard term has been added to the tests in hamiltonian_consistency
- The pseudopotential used in hamiltonian_consistency has been changed with one having atomic orbitals

test/hamiltonian_consistency.jl

Technici4n

I haven't looked at the math yet, but have some comments already

src/eigen/diag.jl

src/postprocess/band_structure.jl

src/terms/hubbard.jl

src/terms/operators.jl

Technici4n · 2025-09-26T15:45:01Z

test/hamiltonian_consistency.jl

    test_consistency_term(ExternalFromFourier(X -> abs(norm(X))))
    test_consistency_term(LocalNonlinearity(ρ -> ρ^2))
    test_consistency_term(Hartree())
+    test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3S"), 1.0))


Suggested change

test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3S"), 1.0))

test_consistency_term(Hubbard(DFTK.OrbitalManifold(; label="3D"), 1.0))

Spacing, and let's use P to test multiple m quantum numbers. And does it work if species is removed as well?

No it doesn't work without species, should it? The intended way should be (in the multiU update that we'll make after this pr is resolved) to add each of the manifolds one by one (so to get a 3d manifold you explicitly ask for all of the single atomic manifolds having 3d orbitals)

test/hubbard.jl

src/terms/hubbard.jl

- New handling of measure unit for U value; - Hubbard matrix has been renamed everywhere as "nhubbard"; - Function "compute_hubbard_nIJ" has been renamed as "compute_nhubbard"; - New documentation for "compute_nhubbard" function; - New documentation for "Hubbard" struct; - New test on Wigner matrices added in tests/hubbard.jl; - Refactoring inside apply! for HubbardUOperator; - Several lines have been shortened for compactness; - Unutilized functions have been removed from terms/hubbard.jl; - Function reshape_hubbard_proj has been updated.

test/hubbard.jl

src/postprocess/band_structure.jl

src/postprocess/dos.jl

src/scf/self_consistent_field.jl

src/terms/hubbard.jl

Technici4n · 2025-09-30T09:53:17Z

src/terms/hubbard.jl

+    U        :: Float64
+end
+function (hubbard::Hubbard)(basis::AbstractBasis)
+    isempty(hubbard.U) && return TermNoop()


I think it's weird. And it's probably not what you're looking for:

julia> isempty(0) false

- Kwarg use_nlcc=false has been added for all meta_gga tests; - Kwarg spin_polarization=:none added to anyonic test, to avoid failure from @Assert length(basis.kpoints) == 1; - Tolerance for condition number in wigner_d_matrix increased to 100, neq decreased to 2*l+2.

ext/DFTKPlotsExt.jl

src/common/spherical_harmonics.jl

src/terms/hubbard.jl

src/terms/operators.jl

test/testcases.jl

Technici4n · 2025-10-09T12:20:03Z

src/terms/hubbard.jl

+    natoms = length(manifold_atoms)  # Number of atoms of the species in the manifold
+    nhubbard = Array{Matrix{Complex{T}}}(undef, nspins, natoms, natoms)
+    p_I = [Vector{Matrix{Complex{T}}}(undef, natoms) for ik in 1:length(basis.kpoints)]
+    # Very low-level, but works


Extract l from first label, make sure that all the labels have the same l, and that each block of 2l+1 consecutive labels is the same except for m.

Technici4n · 2025-10-29T08:15:26Z

@antoine-levitt you said you wanted to have a look, PR is ready!

Technici4n · 2025-11-04T13:11:05Z

@antoine-levitt do you still want to have a look or should we go ahead with this PR?

antoine-levitt · 2025-11-04T14:00:17Z

Sorry I didn't see your previous message (got lost in the flood of notifications). I'll take a look!

antoine-levitt

Looks nice! I just have an issue with the manifold selection API. Needs to be clearer that one atomic orbital = one Hubbard term = one l. I would suggest just

struct Hubbard
U
element::Element
atomic_orbital::String
end
struct TermHubbard
U
P
end

and getting rid of the manifold thing entirely. If you want to +U only selected atoms, then just make them different element types.

src/common/spherical_harmonics.jl

src/scf/self_consistent_field.jl

antoine-levitt · 2025-11-04T14:37:10Z

src/terms/hubbard.jl

+  stating whether the orbital belongs to the manifold.
+"""
+@kwdef struct OrbitalManifold
+    iatom   ::Union{Int64,  Nothing} = nothing


Why is this just an int? If I want to put a +U on two atom types, should I have two hubbard terms?

At the moment it is not possible to have more than one +U correction, but I have a different version of the code which allows "multi-U" calculations with different atom types. I was planning to add it as a different PR once this one is closed.

I think since it very much determines the structure of the code we should do it now

src/terms/hubbard.jl

antoine-levitt · 2025-11-04T14:45:15Z

src/terms/hubbard.jl

+
+  where n or ibnd is the band index, ``weights[ik ibnd] = kweights[ik] * occupation[ik, ibnd]``
+  and ``Pᵢₘ₁`` is the pseudoatomic orbital projector for atom i and orbital m₁
+  (usually just the magnetic quantum number, since l is usually fixed).


This is not a "usually" it's enforced in the code, right? So it's actually enforced that each Hubbard term only deals with one type of atomic orbital. This needs to be enforced by the design.

In the code this is enforced, yes, because the symmetrization of the nhubbard is l-dipendent.

src/terms/hubbard.jl

…to hubbard-v2

mfherbst

Thanks @fsicignano for updating this so quickly. I think it indeed makes for a cleaner code.

examples/hubbard.jl

src/postprocess/dos.jl

src/scf/self_consistent_field.jl

src/symmetry.jl

src/terms/hubbard.jl

src/symmetry.jl

examples/hubbard.jl

src/terms/hubbard.jl

Technici4n · 2025-11-13T08:17:07Z

There is one more thing: currently the example takes a few minutes to run, maybe its discretization should be made coarser.

mfherbst · 2025-11-13T08:55:21Z

There is one more thing: currently the example takes a few minutes to run, maybe its discretization should be made coarser.

Yes, I agree. If we can, we should. But with magnetism that's sometimes not possible.

Technici4n · 2025-11-13T13:32:07Z

I think we are good now!

Technici4n · 2025-11-13T15:47:28Z

@antoine-levitt now is the time if you want to have another look :)

antoine-levitt · 2025-11-14T13:33:36Z

LGTM, thanks for doing this!

fsicignano added 7 commits September 25, 2025 14:54

Adding Hubbard

bc0e2eb

Complex Spherical Harmonics have been removed

d4dfa7b

Correction to testcases.jl

48a33f0

Added test for hubbard term

6c77122

Added explicit import of Hubbard in testcases.jl

5540b22

Speed up of the hubbard test

9e894bf

Hubbard test now functional

4a1d0e2

fsicignano commented Sep 25, 2025

View reviewed changes

test/hamiltonian_consistency.jl Outdated Show resolved Hide resolved

Technici4n reviewed Sep 26, 2025

View reviewed changes

fsicignano added 2 commits September 29, 2025 10:56

Still errors in the consistency test

db7d948

Some comments have been resolved

e0a1002

Technici4n reviewed Sep 29, 2025

View reviewed changes

src/terms/hubbard.jl Outdated Show resolved Hide resolved

fsicignano added 3 commits September 29, 2025 15:06

Some comments have been resolved

7045a6d

Technici4n reviewed Sep 30, 2025

View reviewed changes

fsicignano added 2 commits September 30, 2025 11:54

Added nhubbard field to scfres in forwarddiff_rules.

5aefec2

Update to Wigner_sym, now wigner_d_matrix and relative test

0d69edd

fsicignano force-pushed the hubbard-v2 branch 2 times, most recently from 15ea5c0 to 053450f Compare September 30, 2025 14:31

Change in testing

d0cf01a

- Kwarg use_nlcc=false has been added for all meta_gga tests; - Kwarg spin_polarization=:none added to anyonic test, to avoid failure from @Assert length(basis.kpoints) == 1; - Tolerance for condition number in wigner_d_matrix increased to 100, neq decreased to 2*l+2.

fsicignano force-pushed the hubbard-v2 branch from 053450f to d0cf01a Compare September 30, 2025 14:33

fsicignano added 2 commits October 3, 2025 11:04

New example in examples/hubbard.jl

16a3f6f

use_nlcc=false for all mgga tests

c9d8328

Technici4n requested changes Oct 6, 2025

View reviewed changes

fsicignano added 3 commits October 7, 2025 11:23

Most comments have been addressed

f905197

Several issues solved

7ea2342

Merge branch 'master' into hubbard-v2

3b72e16

mfherbst changed the title ~~Hubbard v2~~ Add Hubbard corrections and DFT+U Oct 7, 2025

Technici4n reviewed Oct 9, 2025

View reviewed changes

mfherbst approved these changes Oct 29, 2025

View reviewed changes

Announce feature

aecf0bd

antoine-levitt reviewed Nov 4, 2025

View reviewed changes

fsicignano added 5 commits November 5, 2025 16:51

Minor changes

1d4d73c

Modifications to the OrbitalManifold interface

15304e7

Merge branch 'hubbard-v2' of https://github.com/fsicignano/DFTK.jl in…

253931d

…to hubbard-v2

Name 'nhubbard' converted to 'hubbard_n'

f649481

Merge branch 'master' into hubbard-v2

df51495

mfherbst reviewed Nov 11, 2025

View reviewed changes

fsicignano and others added 2 commits November 11, 2025 17:01

New constructors added for OrbitalManifold. Testing

5c05a9f

Various changes

a021e32

mfherbst reviewed Nov 12, 2025

View reviewed changes

src/terms/hubbard.jl Outdated Show resolved Hide resolved

src/terms/hubbard.jl Outdated Show resolved Hide resolved

src/terms/hubbard.jl Outdated Show resolved Hide resolved

src/symmetry.jl Outdated Show resolved Hide resolved

examples/hubbard.jl Show resolved Hide resolved

Technici4n added 4 commits November 12, 2025 12:27

Minor tweaks

e2d97d6

Fix tests not testing anything

2da2d03

Reorder OrbitalManifold constructors

9f7dbcd

Move hubbard_n computation down

ca16a44

mfherbst reviewed Nov 12, 2025

View reviewed changes

src/terms/hubbard.jl Show resolved Hide resolved

Review comments

7aab38e

fsicignano and others added 3 commits November 13, 2025 10:51

Minor typo in documentation

3533669

Merge branch 'master' into hubbard-v2

99bddba

Reduce example running time, update manifold naming comment

872c29a

mfherbst merged commit 908e3c8 into JuliaMolSim:master Nov 14, 2025
8 of 10 checks passed

	test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3S"), 1.0))
	test_consistency_term(Hubbard(DFTK.OrbitalManifold(; label="3D"), 1.0))

Add Hubbard corrections and DFT+U #1158

Add Hubbard corrections and DFT+U #1158

Uh oh!

Conversation

fsicignano commented Sep 25, 2025

Uh oh!

Uh oh!

Technici4n left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Technici4n commented Oct 29, 2025

Uh oh!

Technici4n commented Nov 4, 2025

Uh oh!

antoine-levitt commented Nov 4, 2025

Uh oh!

antoine-levitt left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

mfherbst left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!