0.5-Optimization Dev

Project.toml

@@ -1,13 +1,15 @@
 name = "Quiqbox"
 uuid = "7cb8c394-fae1-4ab9-92f2-30189d7746cd"
-version = "0.5.8"
+version = "0.5.9"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6"
+LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -18,15 +20,38 @@ TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
 Combinatorics = "1"
+Documenter = "1"
 DoubleFloats = "^1.2"
 FastGaussQuadrature = "1"
 ForwardDiff = "^0.10.25"
 LBFGSB = "0.4"
+LRUCache = "^1.6"
+LazyArrays = "^1.8"
 LineSearches = "^7.1.1"
 LinearAlgebra = "1"
+Optim = "1"
 Printf = "1"
+QuadGK = "2"
+Random = "1"
 SPGBox = "0.7"
 SpecialFunctions = "2"
 Statistics = "1"
+Suppressor = "0.2"
 TensorOperations = "^4.0.5"
+Test = "1"
 julia = "^1.6"
+libcint_jll = "5"
+
+[extras]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+libcint_jll = "574b78ca-bebd-517c-801d-4735c93a9686"
+
+[targets]
+test = ["Documenter", "ForwardDiff", "LinearAlgebra", "Optim", "QuadGK", "Random", "Suppressor", "Test", "libcint_jll"]

src/Basis.jl

@@ -2070,7 +2070,7 @@ hasNormFactor(b::FloatingGTBasisFuncs) = b.normalizeGTO
 
 
 getNijk(::Type{T}, i::Integer, j::Integer, k::Integer) where {T} = 
-T(πvals[-0.75]) * 2^T(1.5*(i+j+k) + 0.75) * T(sqrt( factorial(i) * factorial(j) * 
+T(πPowers[:n0d75]) * 2^T(1.5*(i+j+k) + 0.75) * T(sqrt( factorial(i) * factorial(j) * 
 factorial(k) / (factorial(2i) * factorial(2j) * factorial(2k)) ))
 
 getNα(i::Integer, j::Integer, k::Integer, α::T) where {T} = 

src/CI.jl

@@ -0,0 +1,129 @@
+struct RSD{𝑚ₛ, ON} <:RCI{𝑚ₛ, ON}
+    N::Int
+    spinUpOccu::NTuple{ON, Bool}
+    spinDnOccu::NTuple{ON, Bool}
+    ID::UInt
+
+    function RSD(spinUpOccu::NTuple{ON, Bool}, spinDnOccu::NTuple{ON, Bool}, 
+                 basis::NTuple{ON, <:GTBasisFuncs{T, D, 1}})
+        Nup = sum(spinUpOccu)
+        Ndn = sum(spinDnOccu)
+        new{Nup-Ndn, ON}(Nup+Ndn, spinUpOccu, spinDnOccu, objectid(basis))
+    end
+end
+
+struct CSF{𝑠, 𝑚ₛ, ON} <: RCI{𝑚ₛ, ON}
+    N::Int
+    spatialOccu::NTuple{ON, Int}
+    ID::UInt
+
+    function RSD(spatialOccu::NTuple{ON, Bool}, basis::NTuple{ON, <:GTBasisFuncs{T, D, 1}})
+
+        # for 
+        new{Nup-Ndn, ON}(Nup+Ndn, spinUpOccu, spinDnOccu, objectid(basis))
+    end
+end
+
+
+
+
+function showbitStr(num, nDigits=4)
+    str = num|>bitstring
+    str[end-nDigits+1:end]
+end
+
+function computeSpinPerm(nUp, nDn) # potential improvement: only iterate over first half.
+    u = (1 << nUp) - 1
+    siz = binomial(nUp+nDn, nUp)
+    v = Array{Int}(undef, siz)
+    v[begin] = u
+    for i in 2:siz
+        t = u | (u - 1)
+        t1 = t + 1
+        t2 = ((~t & t1) - 1) >> (trailing_zeros(u) + 1)
+        u = t1 | t2
+        # @show u
+        # @show showbitStr.([t, t1, t2], nUp+nDn)
+        v[i] = u
+    end
+    # @show showbitStr.(v, nUp+nDn)
+    v
+end
+
+## Input: [1, 0, 1, 2]
+function genSpinStrings(spatialOrbConfig::Vector{Int})
+    ids1 = findall(isone, spatialOrbConfig)
+    ids2 = findall(isequal(2), spatialOrbConfig)
+    nSingleOccu = length(ids1)
+    # nDoubleOccu = sum(ids2)
+    nUp = nSingleOccu ÷ 2
+    nDn = nSingleOccu - nUp
+    permStrs = showbitStr.(computeSpinPerm(nUp, nDn), nSingleOccu)
+    spinUpConfigs = collect(zero(spatialOrbConfig) for _ in eachindex(permStrs))
+    spinDnConfigs = deepcopy(spinUpConfigs)
+    for (perm, SUconfig, SDconfig) in zip(permStrs, spinUpConfigs, spinDnConfigs)
+        SUconfig[ids2] .= 1
+        SDconfig[ids2] .= 1
+        for (spin, idx1) in zip(perm, ids1)
+            config = spin=='1' ? SUconfig : SDconfig
+            config[idx1] = 1
+        end
+    end
+    spinUpConfigs, spinDnConfigs
+end
+
+
+
+# RHF state
+function formatSpatialOrbConfig((RHFSconfig,)::Tuple{NTuple{ON, String}})
+    res = (Int[], Int[], Int[]) # [DoubleOccuIds, SingleOccuIds, UnOccuIds]
+    for (i, e) in enumerate(RHFSconfig)
+        if e == spinOccupations[end]
+            push!(res[1], i)
+        elseif e == spinOccupations[begin]
+            push!(res[3], i)
+        else
+            push!(res[2], i)
+        end
+    end
+    res
+end
+
+function genSpatialOrbConfigCore(n::Int, refConfig::NTuple{3, Vector{Int}})
+    N = 2*length(refConfig[begin]) + length(refConfig[begin+1])
+    for i in min(n, N)
+
+    end
+end
+
+function cc()
+
+end
+
+function promoteElec!(spoConifg, aᵢ, cⱼ)
+    spoConifg[aᵢ] -= 1
+    spoConifg[cⱼ] += 1
+    spoConifg
+end
+
+[2,2,2,0,0,0]
+
+
+
+
+function restrictedSDconfig()
+    getEhf()
+    sum(diag(Hc)[1:n]) + 0.5*sum([(eeI[i,i,j,j] - eeI[i,j,j,i]) for j in 1:n, i in 1:n])
+end
+
+function genFCI()
+
+end
+
+function genXCI()
+
+end
+
+function genSCI()
+
+end

src/Exception.jl

@@ -0,0 +1,16 @@
+function checkFuncReturn(f::F, fName::Symbol, ::Type{TArg}, ::Type{TReturn}) where 
+                        {F<:Function, TArg, TReturn}
+    Base.return_types(f, (TArg,))[] == TReturn || 
+    throw(ArgumentError("`$fName`: `$f` should return a `$TReturn` given an input "*
+                        "argument of type `$TArg`."))
+end
+
+
+function checkCollectionMinLen(data, dataSym::Symbol, minLen::Int)
+    minLen = abs(minLen)
+    str = minlen > 1 ? "s." : "."
+    dataLen = length(data)
+    dataLen < minLen && 
+    throw(ArgumentError("`$dataSym`: $data should contain at least $minLen element"* str))
+    dataLen == minLen
+end

src/HartreeFock.jl

@@ -888,7 +888,7 @@ function runHFcore(bs::GTBasis{T, D, BN, BFT},
     Ntot = (N isa Int) ? N : (N[begin] + N[end])
     Ntot > Nlow || throw(DomainError(N, "$(HFT) requires more than $(Nlow) electrons."))
     Ns = splitSpins(Val(HFT), N)
-    leastNb = max(Ns...)
+    leastNb = maximum(Ns)
     BN < leastNb &&  throw(DomainError(BN, "The number of basis functions should be no "*
                            "less than $(leastNb)."))
     nuc = arrayToTuple(nuc)