Update to Wigner_sym, now wigner_d_matrix and relative test

Author: fsicignano

src/common/spherical_harmonics.jl

@@ -47,3 +47,40 @@ function ylm_real(l::Integer, m::Integer, rvec::AbstractVector{T}) where {T}
 
     error("The case l = $l and m = $m is not implemented")
 end
+
+"""
+ This function returns the Wigner D matrix for real spherical harmonics, 
+    for a given l and symmetry operation, solving a randomized linear system.
+    Such matrix gives the decomposition of a spherical harmonic after application
+    of a symmetry operation back into the basis of spherical harmonics.
+    
+                       Yₗₘ₁(R̂r) = Σₘ₂ D(l,R̂)ₘ₁ₘ₂ * Yₗₘ₂(r)
+
+ The lattice is needed to convert reduced symmetries to Cartesian space.
+"""
+function wigner_d_matrix(l::Integer, Wcart::AbstractMatrix{T}) where {T}
+    D = Matrix{T}(undef, 2*l+1, 2*l+1)
+    if l == 0
+        return D .= 1
+    end
+    rng = Xoshiro(1234)
+    neq = 4*(2*l+1)
+    for m1 in -l:l
+        b = Vector{T}(undef, neq)
+        A = Matrix{T}(undef, neq, 2*l+1)
+        for n in 1:neq
+            r = rand(rng, T, 3)
+            r = r / norm(r)
+            r0 =  Wcart * r
+            b[n] = DFTK.ylm_real(l, m1, r0)
+            for m2 in -l:l
+                A[n,m2+l+1] = DFTK.ylm_real(l, m2, r)
+            end
+        end
+        κ = cond(A)
+        @assert κ < 10.0 "The Wigner matrix computation is badly conditioned. κ(A)=$(κ)"
+        D[m1+l+1,:] = A\b
+    end
+
+    return D
+end

src/postprocess/band_structure.jl

@@ -68,8 +68,9 @@ function compute_bands(scfres::NamedTuple,
                        n_bands=default_n_bands_bandstructure(scfres), kwargs...)
     τ = haskey(scfres, :τ) ? scfres.τ : nothing
     nhubbard = haskey(scfres, :nhubbard) ? scfres.nhubbard : nothing
-    occupation = scfres.occupation
-    compute_bands(scfres.basis, kgrid; scfres.ρ, τ, nhubbard, occupation, scfres.εF, n_bands, kwargs...)
+    compute_bands(scfres.basis, kgrid; 
+                  scfres.ρ, τ, nhubbard, scfres.occupation, 
+                  scfres.εF, n_bands, kwargs...)
 end
 
 """

src/postprocess/dos.jl

@@ -166,7 +166,7 @@ function atomic_orbital_projectors(basis::PlaneWaveBasis{T};
     labels = []
     for (iatom, atom) in enumerate(basis.model.atoms)
         psp = atom.psp
-        @assert !iszero(size(psp.r2_pswfcs[1], 1)) "FATAL ERROR: No Atomic projector found within the provided PseudoPotential."
+        count_n_pswfc(psp)
         for l in 0:psp.lmax, n in 1:DFTK.count_n_pswfc_radial(psp, l)
             label = DFTK.pswfc_label(psp, n, l)
             if !isonmanifold((; iatom, atom.species, label))

src/scf/self_consistent_field.jl

@@ -179,8 +179,8 @@ Overview of parameters:
         if any(needs_τ, basis.terms)
             τ = compute_kinetic_energy_density(basis, ψ, occupation)
         end
-        for (iterm, term) in enumerate(basis.terms)
-            if typeof(term)==DFTK.TermHubbard{Vector{Matrix{ComplexF64}}, Vector{Any}}
+        for term in basis.terms
+            if isa(term, DFTK.TermHubbard)
                 nhubbard = compute_nhubbard(term.manifold, basis, ψ, occupation; projectors=term.P, labels=term.labels).nhubbard
             end
         end
@@ -190,7 +190,6 @@ Overview of parameters:
                        ρin, τ, nhubbard, α=damping, n_iter, nbandsalg.occupation_threshold,
                        runtime_ns=time_ns() - start_ns, nextstate...,
                        diagonalization=[nextstate.diagonalization])
-                       @show nhubbard
 
         # Compute the energy of the new state
         if compute_consistent_energies

src/terms/hubbard.jl

@@ -59,15 +59,14 @@ end
 """
 Symmetrize the Hubbard occupation matrix according to the l quantum number of the manifold.
 """
-function symmetrize_nhub(nhubbard::Array{Matrix{Complex{T}}}, lattice, symmetry, 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!
 
     nspins = size(nhubbard, 1)
     natoms = size(nhubbard, 2)
-    nsym = length(symmetry)
+    nsym = length(symmetries)
     l = Int64((size(nhubbard[1, 1, 1], 1)-1)/2)
-    WigD = Wigner_sym(l, lattice, symmetry)
 
      # Initialize the nhubbard matrix
     ns = Array{Matrix{Complex{T}}}(undef, nspins, natoms, natoms)
@@ -77,52 +76,22 @@ function symmetrize_nhub(nhubbard::Array{Matrix{Complex{T}}}, lattice, symmetry,
                                     size(nhubbard[σ, iatom, jatom], 2))
     end
 
-    for σ in 1:nspins, iatom in 1:natoms, isym in 1:nsym
-        for m1 in 1:size(ns[σ, iatom, iatom], 1), m2 in 1:size(ns[σ, iatom, iatom], 2)
-            sym_atom = find_symmetry_preimage(positions, positions[iatom], symmetry[isym])
-            # 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)
-                ns[σ, iatom, iatom][m1, m2] += WigD[m0, m1, isym] *
-                                               nhubbard[σ, sym_atom, sym_atom][m0, m00] *
-                                               WigD[m00, m2, isym]
-            end
-        end
-    end
-    ns .= ns / nsym
-end
-
-"""
- This function returns the Wigner matrix for a given l and symmetry operation
-    solving a randomized linear system.
- The lattice L is needed to convert reduced symmetries to Cartesian space.
-"""
-function Wigner_sym(l::Int64, L, symmetries::Vector{SymOp{T}}) where {T}
-    nsym = length(symmetries)
-    D = Array{Float64}(undef, 2*l+1, 2*l+1, nsym)
-    if l == 0
-        return D .= 1
-    end
-    Random.seed!(1234)
-    for (isym, symmetry) in enumerate(symmetries)
-        W = symmetry.W
-        for m1 in -l:l
-            b = Vector{Float64}(undef, 2*l+1)
-            A = Matrix{Float64}(undef, 2*l+1, 2*l+1)
-            for n in 1:2*l+1
-                r = rand(Float64, 3)
-                r = r / norm(r)
-                r0 = L * W * inv(L) * r
-                b[n] = DFTK.ylm_real(l, m1, r0)
-                for m2 in -l:l
-                    A[n,m2+l+1] = DFTK.ylm_real(l, m2, r)
+    for symmetry in symmetries
+        Wcart = lattice * symmetry.W * inv(lattice)
+        WigD = wigner_d_matrix(l, Wcart)
+        for σ in 1:nspins, iatom in 1:natoms
+            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)
+                    ns[σ, iatom, iatom][m1, m2] += WigD[m0, m1] *
+                                                   nhubbard[σ, sym_atom, sym_atom][m0, m00] *
+                                                   WigD[m00, m2]
                 end
             end
-            @assert cond(A) > 1/2 "The Wigner matrix computaton is badly conditioned."
-            D[m1+l+1,:,isym] = A\b
         end
     end
-
-    return D
+    ns .= ns / nsym
 end
 
 """
@@ -243,7 +212,7 @@ struct Hubbard
     U        
 end
 function (hubbard::Hubbard)(basis::AbstractBasis)
-    isempty(hubbard.U) && return TermNoop()
+    iszero(hubbard.U) && return TermNoop()
     projs, labs = atomic_orbital_projectors(basis)
     labels, projectors = extract_manifold(basis, projs, labs, hubbard.manifold)
     U = austrip(hubbard.U)

test/hubbard.jl

@@ -1,3 +1,38 @@
+@testitem "Test Wigner matrices on Silicon symmetries" setup=[TestCases] begin
+   using DFTK
+   using PseudoPotentialData
+   using LinearAlgebra
+   
+   # Identity
+   Id = Float64[1 0 0; 0 1 0; 0 0 1]
+   D = DFTK.wigner_d_matrix(1, Id)
+   @test norm(D - I) < 1e-8
+   D = DFTK.wigner_d_matrix(2, Id)
+   @test norm(D - I) < 1e-8
+   # This reverts all p orbitals, sends all d orbitals in themselves
+   Inv = -Id
+   D = DFTK.wigner_d_matrix(1, Inv)
+   @test norm(D + I) < 1e-8
+   D = DFTK.wigner_d_matrix(2, Inv)
+   @test norm(D - I) < 1e-8
+   # This keeps pz, dz2, dx2-y2 and dxy unchanged, changes sign to all others
+   A3  = Float64[1 0 0; 0 -1 0; 0 0 -1] 
+   D3p = Float64[-1 0 0; 0 -1 0; 0 0 1]
+   D3d = Float64[-1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 -1 0; 0 0 0 0 1]
+   D = DFTK.wigner_d_matrix(1, A3)
+   @test norm(D - D3p) < 1e-8
+   D = DFTK.wigner_d_matrix(2, A3)
+   @test norm(D - D3d) < 1e-8
+   # This sends: px <-> py, dxz <-> dyz, dx2-y2 -> -(dx2-y2) and keeps the other fixed
+   A3  = Float64[0 1 0; 1 0 0; 0 0 1] 
+   D3p = Float64[0 0 1; 0 1 0; 1 0 0]
+   D3d = Float64[1 0 0 0 0; 0 0 0 1 0; 0 0 1 0 0; 0 1 0 0 0; 0 0 0 0 -1]
+   D = DFTK.wigner_d_matrix(1, A3)
+   @test norm(D - D3p) < 1e-8
+   D = DFTK.wigner_d_matrix(2, A3)
+   @test norm(D - D3d) < 1e-8
+end 
+
 @testitem "Test Hubbard U term in Nickel Oxide" setup=[TestCases] begin
    using DFTK
    using PseudoPotentialData
@@ -50,25 +85,3 @@
    @test norm(n_hub .- scfres_nosym.nhubbard) < 1e-8
 
 end
-
-@testitem "Test Wigner matrices on Silicon symmetries" setup=[TestCases] begin
-   using DFTK
-   using PseudoPotentialData
-   using LinearAlgebra
-
-   lattice =  [[0 1 1.];
-               [1 0 1.];
-               [1 1 0.]]
-   positions = [ones(3)/8, -ones(3)/8]
-   pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
-   Si = ElementPsp(:Si, pseudopotentials)
-   atoms = [Si, Si]
-   model = model_DFT(lattice, atoms, positions; functionals=PBE())
-   basis = PlaneWaveBasis(model; Ecut=32, kgrid=[2, 2, 2])
-
-   D = DFTK.Wigner_sym(1, lattice, basis.symmetries)
-   D5 = [-1 0 0; 0 -1 0; 0 0 1]
-   @test norm(D[:,:,1] - I) < 1e-8
-   @test norm(D[:,:,25] + I) < 1e-8
-   @test norm(D[:,:,5] - D5) < 1e-8
-end