Change in testing

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

Author: fsicignano

src/common/spherical_harmonics.jl

@@ -64,7 +64,7 @@ function wigner_d_matrix(l::Integer, Wcart::AbstractMatrix{T}) where {T}
         return D .= 1
     end
     rng = Xoshiro(1234)
-    neq = 4*(2*l+1)
+    neq = (2*l+2) # This value should work for p and d orbitals, but can be increased if needed
     for m1 in -l:l
         b = Vector{T}(undef, neq)
         A = Matrix{T}(undef, neq, 2*l+1)
@@ -78,7 +78,7 @@ function wigner_d_matrix(l::Integer, Wcart::AbstractMatrix{T}) where {T}
             end
         end
         κ = cond(A)
-        @assert κ < 10.0 "The Wigner matrix computation is badly conditioned. κ(A)=$(κ)"
+        @assert κ < 100.0 "The Wigner matrix computation is badly conditioned. κ(A)=$(κ)"
         D[m1+l+1,:] = A\b
     end
 

src/scf/self_consistent_field.jl

@@ -181,7 +181,8 @@ Overview of parameters:
         end
         for term in basis.terms
             if isa(term, DFTK.TermHubbard)
-                nhubbard = compute_nhubbard(term.manifold, basis, ψ, occupation; projectors=term.P, labels=term.labels).nhubbard
+                nhubbard = compute_nhubbard(term.manifold, basis, ψ, occupation;
+                                            projectors=term.P, labels=term.labels).nhubbard
             end
         end
 

src/terms/hubbard.jl

@@ -59,7 +59,8 @@ end
 """
 Symmetrize the Hubbard occupation matrix according to the l quantum number of the manifold.
 """
-function symmetrize_nhub(nhubbard::Array{Matrix{Complex{T}}}, lattice, symmetries, positions) where {T}
+function symmetrize_nhub(nhubbard::Array{Matrix{Complex{T}}}, 
+                         lattice, symmetries, positions) where {T}
     # For now we apply symmetries only on nII terms, not on cross-atom terms (nIJ)
     # WARNING: To implement +V this will need to be changed!
 
@@ -83,7 +84,7 @@ function symmetrize_nhub(nhubbard::Array{Matrix{Complex{T}}}, lattice, symmetrie
             sym_atom = find_symmetry_preimage(positions, positions[iatom], symmetry)
             for m1 in 1:size(ns[σ, iatom, iatom], 1), m2 in 1:size(ns[σ, iatom, iatom], 2)
                 # TODO: Here QE flips spin for time-reversal in collinear systems, should we?
-                for m0 in 1:size(nhubbard[σ, iatom, iatom], 1), m00 in 1:size(nhubbard[σ, iatom, iatom], 2)
+                for m0 in 1:size(nhubbard[σ,iatom,iatom],1), m00 in 1:size(nhubbard[σ,iatom,iatom],2)
                     ns[σ, iatom, iatom][m1, m2] += WigD[m0, m1] *
                                                    nhubbard[σ, sym_atom, sym_atom][m0, m00] *
                                                    WigD[m00, m2]
@@ -135,14 +136,15 @@ function compute_nhubbard(manifold::OrbitalManifold,
         # We divide by filled_occ to deal with the physical two spin channels separately.
         ψk, projk, nk = @views ψ[ik], projectors[ik], occupation[ik]/filled_occ  
         c = projk' * ψk      # <ϕ|ψ>
-        # The matrix product is done over the bands. In QE, basis.kweights[ik]*nk[ibnd] would be wg(ik,ibnd)
+        # The matrix product is done over the bands. 
+        # In QE, basis.kweights[ik]*nk[ibnd] would be wg(ik,ibnd)
         n_matrix[σ, :, :] .+= basis.kweights[ik] * c * diagm(nk) * c' 
     end
     n_matrix = mpi_sum(n_matrix, basis.comm_kpts)
 
     # Now I want to reshape it to match the notation used in the papers.
-    # Reshape into n[I, J, σ][m1, m2] where I, J indicate the atom in the Hubbard manifold, σ is the spin, 
-    # m1 and m2 are magnetic quantum numbers (n, l are fixed)
+    # Reshape into n[I, J, σ][m1, m2] where I, J indicate the atom in the Hubbard manifold, 
+    # σ is the spin, m1 and m2 are magnetic quantum numbers (n, l are fixed)
     manifold_atoms = findall(at -> at.species == Symbol(manifold.species), basis.model.atoms)
     natoms = length(manifold_atoms)  # Number of atoms of the species in the manifold
     nhubbard = Array{Matrix{Complex{T}}}(undef, nspins, natoms, natoms)
@@ -178,7 +180,8 @@ function compute_nhubbard(manifold::OrbitalManifold,
     return (; nhubbard, manifold_labels=labels, p_I)
 end
 
-function reshape_hubbard_proj(basis, projectors::Vector{Matrix{Complex{T}}}, labels, manifold) where {T}
+function reshape_hubbard_proj(basis, projectors::Vector{Matrix{Complex{T}}}, 
+                              labels, manifold) where {T}
     manifold_atoms = findall(at -> at.species == Symbol(manifold.species), basis.model.atoms)
     natoms = length(manifold_atoms)
     nprojs = length(labels)

test/hamiltonian_consistency.jl

@@ -5,8 +5,7 @@ using Logging
 using DFTK: mpi_sum
 using LinearAlgebra
 using ..TestCases: silicon
-using PseudoPotentialData
-testcase = silicon 
+testcase = silicon
 
 function test_matrix_repr_operator(hamk, ψk; atol=1e-8)
     for operator in hamk.operators
@@ -94,17 +93,17 @@ end
     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="3P"), 10))
+    test_consistency_term(Hubbard(DFTK.OrbitalManifold(;species=:Si, label="3P"), 0.01), )
     test_consistency_term(Ewald())
     test_consistency_term(PspCorrection())
     test_consistency_term(Xc([:lda_xc_teter93]))
     test_consistency_term(Xc([:lda_xc_teter93]), spin_polarization=:collinear)
     test_consistency_term(Xc([:gga_x_pbe]), spin_polarization=:collinear)
-    test_consistency_term(Xc([:mgga_x_tpss]))
-    test_consistency_term(Xc([:mgga_x_scan]))
-    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(Xc([:mgga_x_tpss]; use_nlcc=false))
+    test_consistency_term(Xc([:mgga_x_scan]; use_nlcc=false))
+    test_consistency_term(Xc([:mgga_c_scan]; use_nlcc=false), spin_polarization=:collinear)
+    test_consistency_term(Xc([:mgga_x_b00]; use_nlcc=false))
+    test_consistency_term(Xc([:mgga_c_b94]; use_nlcc=false), spin_polarization=:collinear)
 
     let
         a = 6
@@ -114,6 +113,6 @@ end
         test_consistency_term(Magnetic(Apot); kgrid=[1, 1, 1], kshift=[0, 0, 0],
                               lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
         test_consistency_term(DFTK.Anyonic(2, 3.2); kgrid=[1, 1, 1], kshift=[0, 0, 0],
-                              lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
+                              lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20, spin_polarization=:none)
     end
 end

test/hubbard.jl

@@ -1,4 +1,4 @@
-@testitem "Test Wigner matrices on Silicon symmetries" setup=[TestCases] begin
+@testitem "Test Wigner matrices" setup=[TestCases] begin
    using DFTK
    using PseudoPotentialData
    using LinearAlgebra