add Mulliken population analysis for RHF

src/Core/Options.jl

@@ -100,8 +100,9 @@ const Default = Dict{String,Union{Float64,Int,String,Bool}}(
                                   "eri_cutoff" => 10^-12,
                                   "findif_disp_size" => 5e-3,
                                   "return_ints" => false,
-                                  "deriv_type" => "findif"
-                                 )
+                                  "deriv_type" => "findif",
+                                  "population_analysis" => true
+                                  )
 """
     Fermi.Options.Current
 

src/Methods/HartreeFock/RHF/RHF.jl

@@ -47,20 +47,21 @@ Computes RHF using information from `Fermi.Options.Current`
 
 These options can be set with `@set <option> <value>`
 
-| Option         | What it does                      | Type      | choices [default]     |
-|----------------|-----------------------------------|-----------|-----------------------|
-| `rhf_alg`      | Picks RHF algorithm               | `Int`     | [1]                   |
-| `scf_max_rms`  | RMS density convergence criterion | `Float64` | [10^-9]               |
-| `scf_max_iter` | Max number of iterations          | `Int`     | [50]                  |
-| `scf_e_conv`   | Energy convergence criterion      | `Float64` | [10^-10]              |
-| `basis`        | What basis set to use             | `String`  | ["sto-3g"]            |
-| `df`           | Whether to use density fitting    | `Bool`    | `true` [`false`]      |
-| `jkfit`        | What aux. basis set to use for JK | `String`  | ["auto"]              |
-| `diis`         | Whether to use DIIS               | `Bool`    | [`true`] `false`      |
-| `oda`          | Whether to use ODA                | `Bool`    | [`true`] `false`      |
-| `oda_cutoff`   | When to turn ODA off (RMS)        | `Float64` | [1E-1]                |
-| `oda_shutoff`  | When to turn ODA off (iter)       | `Int`     | [20]                  |
-| `scf_guess`    | Which guess density to use        | `String`  | "core" ["gwh"]        |
+| Option               | What it does                      | Type      | choices [default]     |
+|----------------------|-----------------------------------|-----------|-----------------------|
+| `rhf_alg`            | Picks RHF algorithm               | `Int`     | [1]                   |
+| `scf_max_rms`        | RMS density convergence criterion | `Float64` | [10^-9]               |
+| `scf_max_iter`       | Max number of iterations          | `Int`     | [50]                  |
+| `scf_e_conv`         | Energy convergence criterion      | `Float64` | [10^-10]              |
+| `basis`              | What basis set to use             | `String`  | ["sto-3g"]            |
+| `df`                 | Whether to use density fitting    | `Bool`    | `true` [`false`]      |
+| `jkfit`              | What aux. basis set to use for JK | `String`  | ["auto"]              |
+| `diis`               | Whether to use DIIS               | `Bool`    | [`true`] `false`      |
+| `oda`                | Whether to use ODA                | `Bool`    | [`true`] `false`      |
+| `oda_cutoff`         | When to turn ODA off (RMS)        | `Float64` | [1E-1]                |
+| `oda_shutoff`        | When to turn ODA off (iter)       | `Int`     | [20]                  |
+| `scf_guess`          | Which guess density to use        | `String`  | "core" ["gwh"]        |
+| `population_analysis`| Whether to do population analysis | `Bool`    | [`true`] `false`      |
 
 # Struct tree
 

src/Methods/HartreeFock/RHF/RHFHelper.jl

@@ -205,6 +205,69 @@ function build_fock!(F::Array{Float64}, H::Array{Float64}, D::Array{Float64}, in
     return F
 end
 
+"""
+    shell_to_atom_index(bset, sh)
+
+Map a basis 'shell' object `sh` to the 1-based atom index in the given
+`bset` (GaussianBasis.BasisSet). The function accepts several common
+representations that shell objects use to indicate their owning atom:
+
+- an integer index (already 1-based or 0-based depending on provider)
+- a field named `:atom` or `:atom_index` containing a numeric index
+- a field named `:center` containing either an Atom-like object or a
+  coordinate vector
+
+"""
+function shell_to_atom_index(bset, sh)
+    fn = fieldnames(typeof(sh))
+    if :atom in fn
+        val = getfield(sh, :atom)
+    elseif :atom_index in fn
+        val = getfield(sh, :atom_index)
+    elseif :center in fn
+        val = getfield(sh, :center)
+    else
+        throw(ErrorException("Cannot determine shell->atom mapping. Shell fields: $(fn)"))
+    end
+
+    # If value is already an integer index, return it
+    if isa(val, Integer)
+        return Int(val)
+    end
+
+    # If value is an Atom-like object, try to find its index in bset.atoms
+    for (i,a) in enumerate(bset.atoms)
+        try
+            if a === val
+                return i
+            end
+            if hasfield(typeof(a), :Z) && hasfield(typeof(val), :Z)
+                if getfield(a, :Z) == getfield(val, :Z) && isapprox(getfield(a, :xyz), getfield(val, :xyz); atol=1e-8)
+                    return i
+                end
+            end
+            if hasfield(typeof(a), :AtomicSymbol) && hasfield(typeof(val), :AtomicSymbol)
+                if getfield(a, :AtomicSymbol) == getfield(val, :AtomicSymbol) && isapprox(getfield(a, :xyz), getfield(val, :xyz); atol=1e-8)
+                    return i
+                end
+            end
+        catch
+            # ignore and continue
+        end
+    end
+
+    # If the value is a coordinate vector, match by coordinates
+    if isa(val, AbstractVector)
+        for (i,a) in enumerate(bset.atoms)
+            if isapprox(getfield(a, :xyz), val; atol=1e-8)
+                return i
+            end
+        end
+    end
+
+    throw(ErrorException("Failed to map shell to atom index for shell value: $(val)"))
+end
+
 ## This is the old version of build_fock! that uses threaded accumulation into thread-specific arrays, 
 # which is no longer recommended. Left here a token of admiraton for the times where threading was simpler to implement and understand RIP.
 function old_build_fock!(F::Array{Float64}, H::Array{Float64}, D::Array{Float64}, ints::IntegralHelper{Float64,SparseERI,AtomicOrbitals})

src/Methods/HartreeFock/RHF/RHFa.jl

@@ -76,6 +76,7 @@ function RHF(ints::IntegralHelper{Float64, <:AbstractERI, AtomicOrbitals}, C::Ab
     oda = Options.get("oda")
     oda_cutoff = Options.get("oda_cutoff")
     oda_shutoff = Options.get("oda_shutoff")
+    population_analysis = Options.get("population_analysis")
 
     # Variables that will get updated iteration-to-iteration
     ite = 1
@@ -218,6 +219,49 @@ function RHF(ints::IntegralHelper{Float64, <:AbstractERI, AtomicOrbitals}, C::Ab
     end
     output(repeat("-",80))
 
+    # Do population analysis?
+    if population_analysis
+        Co = C[:,1:ndocc]
+        @tensor D[u,v] := Co[u,m]*Co[v,m]
+        S = ints["S"]
+
+        bset = ints.orbitals.basisset
+        Nbasis = GaussianBasis.num_basis.(bset.basis)
+        ao_offset = cumsum([0; Nbasis[1:end-1]])
+
+        # Build shell->atom index mapping 
+        shell_atom = [try shell_to_atom_index(bset, sh) catch e; nothing end for sh in bset.basis]
+
+        # Build atom->AO mapping 
+        atom_to_ao = [Int[] for _ = 1:length(bset.atoms)]
+        for sh in 1:bset.nshells
+            idxs = (ao_offset[sh]+1):(ao_offset[sh]+Nbasis[sh])
+            atom_idx = shell_atom[sh]
+            if atom_idx === nothing
+                throw(ErrorException("Failed to obtain atom index for shell $sh"))
+            end
+            ai = Int(atom_idx)
+            append!(atom_to_ao[ai], collect(idxs))
+        end
+        # --- Mulliken population analysis ---
+        pop=D*S
+        popsm = [sum(2*pop[μ,μ] for μ in atom_to_ao[A]) for A in 1:length(atom_to_ao)]
+        # popsm = [sum(2*D[μ,ν]*S[μ,ν] for μ in atom_to_ao[A], ν in 1:nao) for A in 1:length(atom_to_ao)]
+
+        output("\n   • Mulliken population analysis")
+        output(" {:>4} {:>4} {:>12} {:>12}", "Idx", "At", "Population", "Charge")
+        for A in 1:length(popsm)
+            Z = getfield(bset.atoms[A], :Z)
+            sym = hasfield(typeof(bset.atoms[A]), :AtomicSymbol) ? getfield(bset.atoms[A], :AtomicSymbol) : string(Z)
+            output(" {:>4} {:>4} {:12.8f} {:12.8f}", A, sym, popsm[A], Z - popsm[A])
+        end
+        total_mulliken = sum(popsm)
+        nel = molecule.Nα + molecule.Nβ
+        # output(" Mulliken total population: {:12.8f}  Electrons: {:12.8f}  Error: {:12.8f}", total_mulliken, nel, total_mulliken - nel)
+        output(" Mulliken total population: {:12.8f}  Electrons: {:12.8f}  ", total_mulliken, nel)
+
+    end ## population analysis
+
     Orbitals = RHFOrbitals(molecule, ints.basis, eps, E, C)
 
     return RHF(molecule, E, ndocc, nvir, Orbitals, ΔE, Drms)