From ac4b99263cfe113ba63a3a0774ca05036114bb69 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 13 Jan 2025 14:58:54 -0500 Subject: [PATCH 01/29] Added constructors for `StableAdd` and `StableMul`. --- src/Integrals/Engines/Numerical.jl | 2 +- src/Mapping.jl | 3 ++- src/Spatial.jl | 4 ++-- src/SpatialBasis.jl | 8 ++++---- 4 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Integrals/Engines/Numerical.jl b/src/Integrals/Engines/Numerical.jl index 1b9710e8..224b89f4 100644 --- a/src/Integrals/Engines/Numerical.jl +++ b/src/Integrals/Engines/Numerical.jl @@ -19,7 +19,7 @@ end function composeOneBodyKernel(op::O, ::Type{T}, termL::F1, termR::F2) where {O<:Function, T, F1<:Function, F2<:Function} - PairCombine(StableBinary(*, T), adjoint∘termL, op(termR)) + PairCombine(StableMul(T), adjoint∘termL, op(termR)) end function composeOneBodyKernel(op::O, ::Type{T}, term::F) where {O<:Function, T, F<:Function} diff --git a/src/Mapping.jl b/src/Mapping.jl index 74315188..1f0b0a06 100644 --- a/src/Mapping.jl +++ b/src/Mapping.jl @@ -59,7 +59,8 @@ end const StableAdd{T} = StableBinary{T, typeof(+)} const StableMul{T} = StableBinary{T, typeof(*)} -StableBinary(f::Function) = Base.Fix1(StableBinary, f) +StableAdd(::Type{T}) where {T} = StableBinary(+, T) +StableMul(::Type{T}) where {T} = StableBinary(*, T) struct Retrieve{P<:CompositePointer} <: FunctionComposer diff --git a/src/Spatial.jl b/src/Spatial.jl index c1f95279..f71eb6e3 100644 --- a/src/Spatial.jl +++ b/src/Spatial.jl @@ -92,7 +92,7 @@ function unpackParamFuncCore!(f::AxialProdFunc{T, D}, paramSet::FlatParamSet) wh fEvalComps[i] = InsertInward(fInner, (OnlyHead∘Retrieve)(ptr)) map(x->(ChainPointer(anchor, x.first)=>x.second), axialPairs) end - fEvalCore = Tuple(fEvalComps) |> (ChainReduce∘StableBinary)(*, T) + fEvalCore = Tuple(fEvalComps) |> (ChainReduce∘StableMul)(T) EvalAxialProdFunc(fEvalCore), paramSet, pairs end @@ -122,7 +122,7 @@ function unpackParamFuncCore!(f::PolyRadialFunc{T, D}, paramSet::FlatParamSet) w end coordEncoder = InsertInward(fInner, OnlyHead(LinearAlgebra.norm)) angularFunc = (OnlyHead∘ReturnTyped)(f.angular, T) - fEvalCore = PairCombine(StableBinary(*, T), coordEncoder, angularFunc) + fEvalCore = PairCombine(StableMul(T), coordEncoder, angularFunc) EvalPolyRadialFunc(fEvalCore), paramSet, radialPairs end diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 10f71704..0d8aea5e 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -74,7 +74,7 @@ end function ScaledOrbital(orb::EvalComposedOrb{T}, scalar::Function, scope::FlatParamSetFilter{T}) where {T} - fCoreLocal = PairCombine(StableBinary(*, T), orb, scalar) + fCoreLocal = PairCombine(StableMul(T), orb, scalar) fCore = ParamFilterFunc(fCoreLocal, AwaitFilter(scope)) ScaledOrbital(fCore) end @@ -212,7 +212,7 @@ function getEffectiveWeight(o::CompositeOrb{T}, weight::FlattenedParam{T, 1}, "`CompositeOrb` with another value is prohibited.")) len = first(outputSizeOf(o.weight)) outWeight = indexParam(weight, idx) - map(i->CellParam(StableBinary(*, T), indexParam(o.weight, i), outWeight, :w), 1:len) + map(i->CellParam(StableMul(T), indexParam(o.weight, i), outWeight, :w), 1:len) end function getEffectiveWeight(o::AbstractVector{<:ComposedOrb{T, D}}, @@ -270,7 +270,7 @@ function restrainEvalOrbType(weightedFs::AbstractVector{<:WeightedPF{T, D}}) whe fInnerType = (eltype∘map)(f->getField(f, cPtr), fInnerObjs) nInnerType = (eltype∘map)(f->getField(f, ChainPointer( (:right, :f) )), fInnerObjs) V = compressWeightedPF(fInnerType, nInnerType) - ChainReduce(StableBinary(+, T), V(weightedFs)) + ChainReduce(StableAdd(T), V(weightedFs)) end struct CompOrbParamPtr{T, D, R<:FieldPtrDict{T}, @@ -297,7 +297,7 @@ function unpackComposedOrbCore!(f::CompositeOrb{T, D}, paramSet::FlatParamSet, i += 1 fInnerCore, _, basisPtr = unpackParamFunc!(b, pSetLocal, paramSetId) getWeight = Retrieve(ChainPointer( weightPtr, ChainPointer(i, TensorType(T)) )) - weightedPrimOrb = PairCombine(StableBinary(*, T), fInnerCore, OnlyBody(getWeight)) + weightedPrimOrb = PairCombine(StableMul(T), fInnerCore, OnlyBody(getWeight)) push!(weightedFields, weightedPrimOrb) basisPtr end From d19053be506fc4f988973cdd74c9707a0996e4f2 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 13 Jan 2025 15:00:02 -0500 Subject: [PATCH 02/29] Moved `Identity`. --- src/Integrals/Operators.jl | 3 +++ src/Mapping.jl | 5 ----- src/Quiqbox.jl | 1 + 3 files changed, 4 insertions(+), 5 deletions(-) create mode 100644 src/Integrals/Operators.jl diff --git a/src/Integrals/Operators.jl b/src/Integrals/Operators.jl new file mode 100644 index 00000000..2d238240 --- /dev/null +++ b/src/Integrals/Operators.jl @@ -0,0 +1,3 @@ +struct Identity <: DirectOperator end + +(::Identity)(f::Function) = itself(f) \ No newline at end of file diff --git a/src/Mapping.jl b/src/Mapping.jl index 1f0b0a06..07a4249a 100644 --- a/src/Mapping.jl +++ b/src/Mapping.jl @@ -341,11 +341,6 @@ function unpackTypedFunc!(f::ReturnTyped{T}, paramSet::AbstractVector, end -struct Identity <: DirectOperator end - -(::Identity)(f::Function) = itself(f) - - struct LeftPartial{F<:Function, A<:NonEmptyTuple{Any}} <: FunctionModifier f::F header::A diff --git a/src/Quiqbox.jl b/src/Quiqbox.jl index eb9d58b6..9e93267c 100644 --- a/src/Quiqbox.jl +++ b/src/Quiqbox.jl @@ -28,6 +28,7 @@ include("SpatialBasis.jl") # include("Matter.jl") # include("Overload.jl") +include("Integrals/Operators.jl") include("Integrals/Engines/Numerical.jl") include("Integrals/Engines/GaussianOrbitals.jl") include("Integrals/Framework.jl") From 5146c31e4fee787125d5101333741ee80c2ff7e9 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 13 Jan 2025 21:21:07 -0500 Subject: [PATCH 03/29] Code optimization. --- src/Integrals/Framework.jl | 2 +- src/Lexicons.jl | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 9a43b6ba..769ef260 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -33,9 +33,9 @@ const OneBodyIdxSymDict = let tempDict=Base.ImmutableDict(true=>:aa) end const TwoBodyIdxSymDict = let - valTemp = (:aaaa, :aabb, :abab, :aaxy, :abxx, :abxy) keyTemp = ((true, true, true), (true, true, false), (false, false, true), (true, false, false), (false, true, false), (false, false, false)) + valTemp = (:aaaa, :aabb, :abab, :aaxy, :abxx, :abxy) mapreduce(Base.ImmutableDict, keyTemp, valTemp, init=Base.ImmutableDict{NTuple{3, Bool}, Symbol}()) do key, val key=>val diff --git a/src/Lexicons.jl b/src/Lexicons.jl index 1dfcc9cf..6706c660 100644 --- a/src/Lexicons.jl +++ b/src/Lexicons.jl @@ -1,2 +1,8 @@ -const πPowers = Dict( [:n0d75, :p0d5, :p1d0, :p1d5, :p2d5] .=> big(π).^ - [ -0.75, 0.5, 1.0, 1.5, 2.5] ) \ No newline at end of file +const πPowers = let + keyTemp = (:n0d75, :p0d5, :p1d0, :p1d5, :p2d5) + valTemp = big(π).^ (-0.75, 0.5, 1.0, 1.5, 2.5) + mapreduce(Base.ImmutableDict, keyTemp, valTemp, + init=Base.ImmutableDict{Symbol, BigFloat}()) do key, val + key=>val + end +end \ No newline at end of file From 4e6456d27c5f61edd0ac46e1f87994cb348b7e6f Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 13 Jan 2025 22:09:53 -0500 Subject: [PATCH 04/29] Extended factorial computation range and added the corresponding cache. --- src/Integrals/Engines/GaussianOrbitals.jl | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 9d1fb6d5..57adb00c 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -1,12 +1,22 @@ using LinearAlgebra: norm +using LRUCache -function oddFactorialCore(a::Int) # a * (a-2) * ... * 1 - factorial(2a) ÷ (2^a * factorial(a)) +const DefaultOddFactorialCacheSizeLimit = 25 +const OddFactorialCache = LRU{Int, BigInt}(maxsize=DefaultOddFactorialCacheSizeLimit) + +function oddFactorial(a::Int) # a * (a-2) * ... * 1 + get!(OddFactorialCache, a) do + i = (a > 33 ? BigInt(1) : 1) + for j = 1:2:a + i *= j + end + i + end end function polyGaussFuncSquaredNorm(α::T, degree::Int) where {T<:Real} - factor = degree > 0 ? (oddFactorialCore(2degree - 1) / (4α)^degree) : one(T) + factor = degree > 0 ? (T(oddFactorial(2degree - 1)) / (4α)^degree) : one(T) T(πPowers[:p0d5]) / sqrt(2α) * factor end From 06dcb53494526075b3f3af5461d121df10ed9b75 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 13 Jan 2025 23:46:28 -0500 Subject: [PATCH 05/29] Removed flag in case for different integer-type inputs. --- src/Integrals/Engines/GaussianOrbitals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 57adb00c..36bf3968 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -6,7 +6,7 @@ const OddFactorialCache = LRU{Int, BigInt}(maxsize=DefaultOddFactorialCacheSizeL function oddFactorial(a::Int) # a * (a-2) * ... * 1 get!(OddFactorialCache, a) do - i = (a > 33 ? BigInt(1) : 1) + i = BigInt(1) for j = 1:2:a i *= j end From 034dd19dee3ae9a93b0180b1d2559e576fbc2388 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Wed, 15 Jan 2025 23:22:07 -0500 Subject: [PATCH 06/29] Reorganized and optimized core functions. --- src/Integrals/Engines/GaussianOrbitals.jl | 141 ++++++++++++---------- 1 file changed, 79 insertions(+), 62 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 36bf3968..80bd157f 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -21,83 +21,94 @@ function polyGaussFuncSquaredNorm(α::T, degree::Int) where {T<:Real} end -struct XpnPair{T<:Real} - left::T - right::T - sum::T - prod::T - - function XpnPair(l::T, r::T) where {T} - checkPositivity(l) - checkPositivity(r) - new{T}(l, r, l+r, l*r) - end +function gaussProdCore1(dxLR::T, xpnProdOverSum::T) where {T<:Real} + exp(- xpnProdOverSum * dxLR^2) end - -struct CenPair{T<:Real, D} - left::NTuple{D, T} - right::NTuple{D, T} - dist::T - - CenPair(l::NonEmptyTuple{T, D}, r::NonEmptyTuple{T, D}) where {T, D} = - new{T, D+1}(l, r, norm(l .- r)) +function gaussProdCore2(dxML::T, dxMR::T, lxL::Int, lxR::Int, lx::Int) where {T<:Real} + lb = max(-lx, lx - 2lxR) + ub = min( lx, 2lxL - lx ) + res = zero(T) + for q in lb:2:ub + i = (lx + q) >> 1 + j = (lx - q) >> 1 + res += binomial(lxL, i) * binomial(lxR, j) * dxML^(lxL - i) * dxMR^(lxR - j) + end + res end -function (cPair::CenPair{T})(xPair::XpnPair{T}) where {T} - (xPair.left .* cPair.left .+ xPair.right .* cPair.right) ./ xPair.sum +function gaussProdCore2(lxL::Int, lxR::Int, lx::Int) + lb = max(-lx, lx - 2lxR) + ub = min( lx, 2lxL - lx ) + res = zero(Int) + for q in lb:2:ub + res += Int(isequal(2lxL, lx+q) * isequal(2lxR, lx-q)) + end + res end -struct AngPair{D} - left::NTuple{D, Int} - right::NTuple{D, Int} - sum::NTuple{D, Int} +struct PrimGaussOrbData{T, D} + cen::NTuple{D, T} + xpn::T + ang::NTuple{D, Int} - AngPair(l::NonEmptyTuple{Int, D}, r::NonEmptyTuple{Int, D}) where {D} = - new{D+1}(l, r, l .+ r) + PrimGaussOrbData(cen::NonEmptyTuple{T, D}, xpn::T, + ang::NonEmptyTuple{Int, D}) where {T, D} = + new{T, D+1}(cen, xpn, ang) end -function gaussProdCore1(cPair::CenPair{T}, xPair::XpnPair{T}) where {T<:Real} - if iszero(cPair.dist) - T(1) - else - exp(- xPair.prod / xPair.sum * cPair.dist^2) - end -end +struct GaussProductData{T, D} + lhs::PrimGaussOrbData{T, D} + rhs::PrimGaussOrbData{T, D} + cen::NTuple{D, T} + xpn::T -function gaussProdCore2(x1::T, x2::T, x::T, lx1::Int, lx2::Int, lx::Int) where {T<:Real} - lb = max(-lx, lx - 2lx2) - ub = min( lx, 2lx1 - lx ) - res = zero(T) - for q in lb:2:ub - i = (lx + q) ÷ 2 - j = (lx - q) ÷ 2 - res += binomial(lx1, i) * binomial(lx2, j) * (x - x1)^(lx1 -i) * (x - x2)^(lx2 -j) + function GaussProductData(oData1::PrimGaussOrbData{T, D}, + oData2::PrimGaussOrbData) where {T, D} + xpnSum = oData1.xpn + oData2.xpn + cenNew = (oData1.xpn .* oData1.cen .+ oData2.xpn .* oData2.cen) ./ xpnSum + new{T, D}(oData1, oData2, cenNew, xpnSum) end - res end -function overlapPGTO(cPair::CenPair{T, D}, xPair::XpnPair{T}, - aPair::AngPair{D}) where {T, D} - α = xPair.sum - mapreduce(*, cPair.left, cPair.right, cPair(xPair), aPair.left, aPair.right, - aPair.sum) do x1, x2, x, i1, i2, i - mapreduce(+, 0:(i÷2)) do j - gaussProdCore2(x1, x2, x, i1, i2, 2j) * polyGaussFuncSquaredNorm(α/2, j) +function overlapPGTO(data::GaussProductData{T}) where {T} + cenL = data.lhs.cen + cenR = data.rhs.cen + cenM = data.cen + + xpn = data.xpn + xpnRatio = data.lhs.xpn * data.rhs.xpn / xpn + + angL = data.lhs.ang + angR = data.rhs.ang + + mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR + dxLR = xL - xR + jRange = 0:((iL + iR) ÷ 2) + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore2(iL, iR, 2j) * polyGaussFuncSquaredNorm(xpn/2, j) + end + else + mapreduce(+, jRange) do j + xML = xM - xL + xMR = xM - xR + gaussProdCore2(xML, xMR, iL, iR, 2j) * polyGaussFuncSquaredNorm(xpn/2, j) + end * gaussProdCore1(dxLR, xpnRatio) end - end * gaussProdCore1(cPair, xPair) + end end - -function overlapPGTO(xpn::T, ang::NTuple{D, Int}) where {T, D} +function overlapPGTO(xpn::T, ang::NonEmptyTuple{Int}) where {T} mapreduce(*, ang) do i polyGaussFuncSquaredNorm(xpn, i) end end + function getAngularNums(orb::PrimATOcore) angFunc = orb.f.apply.f.right.f angFunc.f.m.tuple @@ -112,6 +123,7 @@ function getCenCoordPtr(orb::PrimitiveOrbCore) orb.f.dress.select end + function preparePGTOconfig(orb::PrimGTOcore) cenIds = getCenCoordPtr(orb) xpnIdx = getExponentPtr(orb) @@ -119,6 +131,7 @@ function preparePGTOconfig(orb::PrimGTOcore) cenIds, xpnIdx, ang end + struct OverlapGTOrbSelf{T, D} <: OrbitalIntegrator{T, D} xpn::FlatPSetInnerPtr{T} ang::NTuple{D, Int} @@ -129,10 +142,6 @@ function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} overlapPGTO(xpnVal, f.ang) end -function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} - _, xpnIdx, ang = preparePGTOconfig(orb) - OverlapGTOrbSelf(xpnIdx, ang) -end struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} cen::NTuple{2, NTuple{ D, FlatPSetInnerPtr{T} }} @@ -142,18 +151,26 @@ end function (f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}) where {T} - cen1 = getField(pars1, f.cen[1]) - cen2 = getField(pars2, f.cen[2]) - xpn1 = getField(pars1, f.xpn[1]) - xpn2 = getField(pars2, f.xpn[2]) - overlapPGTO(CenPair(cen1, cen2), XpnPair(xpn1, xpn2), AngPair(f.ang...)) + c1 = getField(pars1, first(f.cen)) + x1 = getField(pars1, first(f.xpn)) + d1 = PrimGaussOrbData(c1, x1, first(f.ang)) + c2 = getField(pars2, last(f.cen)) + x2 = getField(pars2, last(f.xpn)) + d2 = PrimGaussOrbData(c2, x2, last(f.ang)) + GaussProductData(d1, d2) |> overlapPGTO end + function genGTOrbOverlapFunc(orbs::NTuple{2, PrimGTOcore{T, D}}) where {T, D} configs = preparePGTOconfig.(orbs) OverlapGTOrbPair(getindex.(configs, 1), getindex.(configs, 2), getindex.(configs, 3)) end +function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} + _, xpnIdx, ang = preparePGTOconfig(orb) + OverlapGTOrbSelf(xpnIdx, ang) +end + function buildNormalizerCore(o::PrimGTOcore{T, D}) where {T, D} genGTOrbOverlapFunc((o,)) From 8f15979d258ad862b7fac956c18cbd31c20b1cfc Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Fri, 17 Jan 2025 01:28:12 -0500 Subject: [PATCH 07/29] Code optimization. --- src/Integrals/Engines/GaussianOrbitals.jl | 61 +++++++++++++---------- 1 file changed, 34 insertions(+), 27 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 80bd157f..9bd33fa3 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -15,14 +15,9 @@ function oddFactorial(a::Int) # a * (a-2) * ... * 1 end -function polyGaussFuncSquaredNorm(α::T, degree::Int) where {T<:Real} - factor = degree > 0 ? (T(oddFactorial(2degree - 1)) / (4α)^degree) : one(T) - T(πPowers[:p0d5]) / sqrt(2α) * factor -end - - -function gaussProdCore1(dxLR::T, xpnProdOverSum::T) where {T<:Real} - exp(- xpnProdOverSum * dxLR^2) +function gaussProdCore1(xpnLRsum::T, degree::Int) where {T<:Real} + factor = degree > 0 ? (T(oddFactorial(2degree - 1)) / (2xpnLRsum)^degree) : one(T) + T(πPowers[:p0d5]) / sqrt(xpnLRsum) * factor end function gaussProdCore2(dxML::T, dxMR::T, lxL::Int, lxR::Int, lx::Int) where {T<:Real} @@ -47,6 +42,14 @@ function gaussProdCore2(lxL::Int, lxR::Int, lx::Int) res end +function gaussProdCore3(dxLR::T, xpnProdOverSum::T) where {T<:Real} + exp(- xpnProdOverSum * dxLR^2) +end + +function gaussProdCore4(xpnSum::T, dxLR::T, xpnRatio::T) where {T} + gaussProdCore1(xpnSum, 0) * gaussProdCore3(dxLR, xpnRatio) +end + struct PrimGaussOrbData{T, D} cen::NTuple{D, T} @@ -79,32 +82,36 @@ function overlapPGTO(data::GaussProductData{T}) where {T} cenR = data.rhs.cen cenM = data.cen - xpn = data.xpn - xpnRatio = data.lhs.xpn * data.rhs.xpn / xpn + xpnS = data.xpn + xpnRatio = data.lhs.xpn * data.rhs.xpn / xpnS angL = data.lhs.ang angR = data.rhs.ang mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR dxLR = xL - xR - jRange = 0:((iL + iR) ÷ 2) - if isequal(dxLR, zero(T)) - mapreduce(+, jRange) do j - gaussProdCore2(iL, iR, 2j) * polyGaussFuncSquaredNorm(xpn/2, j) - end + if iL == iR == 0 + gaussProdCore4(xpnS, dxLR, xpnRatio) else - mapreduce(+, jRange) do j - xML = xM - xL - xMR = xM - xR - gaussProdCore2(xML, xMR, iL, iR, 2j) * polyGaussFuncSquaredNorm(xpn/2, j) - end * gaussProdCore1(dxLR, xpnRatio) + jRange = 0:((iL + iR) ÷ 2) + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore1(xpnS, j) * gaussProdCore2(iL, iR, 2j) + end + else + mapreduce(+, jRange) do j + xML = xM - xL + xMR = xM - xR + gaussProdCore1(xpnS, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) + end * gaussProdCore3(dxLR, xpnRatio) + end end end end function overlapPGTO(xpn::T, ang::NonEmptyTuple{Int}) where {T} mapreduce(*, ang) do i - polyGaussFuncSquaredNorm(xpn, i) + gaussProdCore1(2xpn, i) end end @@ -151,12 +158,12 @@ end function (f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}) where {T} - c1 = getField(pars1, first(f.cen)) - x1 = getField(pars1, first(f.xpn)) - d1 = PrimGaussOrbData(c1, x1, first(f.ang)) - c2 = getField(pars2, last(f.cen)) - x2 = getField(pars2, last(f.xpn)) - d2 = PrimGaussOrbData(c2, x2, last(f.ang)) + c1 = getField(pars1, f.cen[begin]) + x1 = getField(pars1, f.xpn[begin]) + d1 = PrimGaussOrbData(c1, x1, f.ang[begin]) + c2 = getField(pars2, f.cen[end]) + x2 = getField(pars2, f.xpn[end]) + d2 = PrimGaussOrbData(c2, x2, f.ang[end]) GaussProductData(d1, d2) |> overlapPGTO end From f857dfcd047305e909e1105953be6973f3a202bb Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Fri, 17 Jan 2025 12:40:48 -0500 Subject: [PATCH 08/29] Added caching to the overlap computation of PGTO. --- src/Integrals/Engines/GaussianOrbitals.jl | 54 +++++++++++++++-------- 1 file changed, 36 insertions(+), 18 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 9bd33fa3..2e2734ad 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -46,7 +46,7 @@ function gaussProdCore3(dxLR::T, xpnProdOverSum::T) where {T<:Real} exp(- xpnProdOverSum * dxLR^2) end -function gaussProdCore4(xpnSum::T, dxLR::T, xpnRatio::T) where {T} +function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T} gaussProdCore1(xpnSum, 0) * gaussProdCore3(dxLR, xpnRatio) end @@ -77,6 +77,40 @@ struct GaussProductData{T, D} end +const DefaultPrimGTOIntCacheSizeLimit = 100 + +@generated function lazyOverlapPGTO(input::Tuple{T, T, T, T, Int, Int}) where {T} + cache = LRU{Tuple{T, T, T, T, Int, Int}, T}(maxsize=DefaultPrimGTOIntCacheSizeLimit) + return quote + if input[end-1] == input[end] == 0 + gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) + else + res = get($cache, input, nothing) + if res === nothing + res = overlapPGTOcore(input) + setindex!($cache, res, input) + end + res + end + end +end + +function overlapPGTOcore(input::Tuple{T, T, T, T, Int, Int}) where {T} + xpnSum, xpnRatio, xML, xMR, iL, iR = input + jRange = 0:((iL + iR) ÷ 2) + dxLR = xMR - xML + + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) + end + else + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) + end * gaussProdCore3(dxLR, xpnRatio) + end +end + function overlapPGTO(data::GaussProductData{T}) where {T} cenL = data.lhs.cen cenR = data.rhs.cen @@ -89,23 +123,7 @@ function overlapPGTO(data::GaussProductData{T}) where {T} angR = data.rhs.ang mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR - dxLR = xL - xR - if iL == iR == 0 - gaussProdCore4(xpnS, dxLR, xpnRatio) - else - jRange = 0:((iL + iR) ÷ 2) - if isequal(dxLR, zero(T)) - mapreduce(+, jRange) do j - gaussProdCore1(xpnS, j) * gaussProdCore2(iL, iR, 2j) - end - else - mapreduce(+, jRange) do j - xML = xM - xL - xMR = xM - xR - gaussProdCore1(xpnS, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) - end * gaussProdCore3(dxLR, xpnRatio) - end - end + lazyOverlapPGTO((xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end From 206c39cb169cc94ff681e81275e0e491645ce2f7 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Fri, 17 Jan 2025 21:53:07 -0500 Subject: [PATCH 09/29] Fixed type constraint for `NormalizeCompOrb` with inhomogeneous core integrators. --- src/SpatialBasis.jl | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 0d8aea5e..3d3922c0 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -23,7 +23,7 @@ end const GetParamSubset{T} = Retrieve{AwaitFilter{FlatParamSetFilter{T}}} const VariedNormCore{T, D, N, F<:OrbitalIntegrator{T, D}} = - ReturnTyped{T, <:EncodeApply{N, NTuple{N, GetParamSubset{T}}, F}} + ReturnTyped{T, EncodeApply{N, NTuple{N, GetParamSubset{T}}, F}} const UnitScalarCore{T} = ReturnTyped{T, Unit{T}} @@ -486,30 +486,51 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} diagFuncs = Memory{ReturnTyped{T}}(undef, nOrbs) pOrbCores = Memory{PrimitiveOrbCore{T, D}}(undef, nOrbs) pOrbScopes = Memory{AwaitFilter{FlatParamSetFilter{T}}}(undef, nOrbs) + hasUnit = false + diagFuncCoreType = Union{} + utriFuncCoreType = Union{} for i in eachindex(pOrbs) o = pOrbs[i] oCore = o.f.apply.left oScope = first(o.f.scope) diagFuncs[i] = if isRenormalized(o) + hasUnit = true ReturnTyped(Unit(T), T) else nfInner = genOneBodyCoreIntegrator(Identity(), (oCore,)) + diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfInner)) ReturnTyped(EncodeApply((Retrieve(oScope),), nfInner), T) end pOrbCores[i] = oCore pOrbScopes[i] = oScope end + diagFuncType = if isconcretetype(diagFuncCoreType) + VariedNormCore{T, D, 1, diagFuncCoreType} + else + VariedNormCore{T, D, 1, <:diagFuncCoreType} + end + hasUnit && (diagFuncType = Union{UnitScalarCore{T}, diagFuncType}) + diagFuncs = convert(Memory{diagFuncType}, diagFuncs) + utriFuncs = Memory{ReturnTyped{T}}(undef, nOrbs*(nOrbs-1)÷2) for j in 1:length(utriFuncs) n, m = convert1DidxTo2D(nOrbs-1, j) corePairs = (pOrbCores[begin+m-1], pOrbCores[begin+n-1]) nfInner = genOneBodyCoreIntegrator(Identity(), corePairs) + utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfInner)) filters = (pOrbScopes[begin+m-1], pOrbScopes[begin+n-1]) utriFuncs[begin+j-1] = ReturnTyped(EncodeApply(Retrieve.(filters), nfInner), T) end - coreTransform = HermitianContract(map(itself, diagFuncs), map(itself, utriFuncs)) + utriFuncType = if isconcretetype(utriFuncCoreType) + VariedNormCore{T, D, 2, utriFuncCoreType} + else + VariedNormCore{T, D, 2, <:utriFuncCoreType} + end + utriFuncs = convert(Memory{utriFuncType}, utriFuncs) + + coreTransform = HermitianContract(diagFuncs, utriFuncs) getWeights = map(b->ReturnTyped(b.right.f, T), weightedOrbs) NormalizeCompOrb(coreTransform, getWeights) end From 3c779e25da3e4a0d475c032920934dcf18c00c8b Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Fri, 17 Jan 2025 21:55:06 -0500 Subject: [PATCH 10/29] Added multi-cache support for different GTO subshells. --- src/Integrals/Engines/GaussianOrbitals.jl | 86 +++++++++++++---------- 1 file changed, 50 insertions(+), 36 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 2e2734ad..dce8d22a 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -51,48 +51,60 @@ function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T} end -struct PrimGaussOrbData{T, D} +struct PrimGaussOrbData{T, D, S} cen::NTuple{D, T} xpn::T - ang::NTuple{D, Int} + ang::WeakComp{D, S} - PrimGaussOrbData(cen::NonEmptyTuple{T, D}, xpn::T, - ang::NonEmptyTuple{Int, D}) where {T, D} = - new{T, D+1}(cen, xpn, ang) + PrimGaussOrbData(cen::NTuple{D, T}, xpn::T, ang::WeakComp{D, S}) where {T, D, S} = + new{T, D, S}(cen, xpn, ang) end -struct GaussProductData{T, D} - lhs::PrimGaussOrbData{T, D} - rhs::PrimGaussOrbData{T, D} +struct GaussProductData{T, D, S1, S2} + lhs::PrimGaussOrbData{T, D, S1} + rhs::PrimGaussOrbData{T, D, S2} cen::NTuple{D, T} xpn::T - function GaussProductData(oData1::PrimGaussOrbData{T, D}, - oData2::PrimGaussOrbData) where {T, D} + function GaussProductData(oData1::PrimGaussOrbData{T, D, S1}, + oData2::PrimGaussOrbData{T, D, S2}) where {T, D, S1, S2} xpnSum = oData1.xpn + oData2.xpn cenNew = (oData1.xpn .* oData1.cen .+ oData2.xpn .* oData2.cen) ./ xpnSum - new{T, D}(oData1, oData2, cenNew, xpnSum) + new{T, D, S1, S2}(oData1, oData2, cenNew, xpnSum) end end -const DefaultPrimGTOIntCacheSizeLimit = 100 +const DefaultPGTOrbOverlapCacheSizeLimit = 100 -@generated function lazyOverlapPGTO(input::Tuple{T, T, T, T, Int, Int}) where {T} - cache = LRU{Tuple{T, T, T, T, Int, Int}, T}(maxsize=DefaultPrimGTOIntCacheSizeLimit) - return quote - if input[end-1] == input[end] == 0 - gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) - else - res = get($cache, input, nothing) - if res === nothing - res = overlapPGTOcore(input) - setindex!($cache, res, input) - end - res - end +struct NullCache{T} <: QueryBox{T} end + +@generated function gen1DimPGTOrbOverlapCache(::Type{T}, ::Val{S}) where {T, S} + if S==0 + return :( NullCache{T}() ) + else + K = Tuple{T, T, T, T, Int, Int} + cache = LRU{K, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit) + return :( $cache ) + end +end + +function lazyOverlapPGTO(::NullCache{T}, input::Tuple{T, T, T, T, Int, Int}) where {T} + gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) +end + +function lazyOverlapPGTO(cache::LRU{K, T}, input::K) where + {T, K<:Tuple{T, T, T, T, Int, Int}} + # get!(cache, input) do + # overlapPGTOcore(input) + # end + res = get(cache, input, nothing) # Fewer allocations than using `get!` + if res === nothing + res = overlapPGTOcore(input) + setindex!(cache, res, input) end + res end function overlapPGTOcore(input::Tuple{T, T, T, T, Int, Int}) where {T} @@ -111,7 +123,7 @@ function overlapPGTOcore(input::Tuple{T, T, T, T, Int, Int}) where {T} end end -function overlapPGTO(data::GaussProductData{T}) where {T} +function overlapPGTO(data::GaussProductData{T, D, S1, S2}) where {T, D, S1, S2} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -122,8 +134,10 @@ function overlapPGTO(data::GaussProductData{T}) where {T} angL = data.lhs.ang angR = data.rhs.ang - mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR - lazyOverlapPGTO((xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) + intCache = gen1DimPGTOrbOverlapCache(T, Val(S1+S2)) + + mapreduce(*, cenL, cenR, cenM, angL.tuple, angR.tuple) do xL, xR, xM, iL, iR + lazyOverlapPGTO(intCache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end @@ -134,9 +148,9 @@ function overlapPGTO(xpn::T, ang::NonEmptyTuple{Int}) where {T} end -function getAngularNums(orb::PrimATOcore) +function getAngMomentum(orb::PrimATOcore) angFunc = orb.f.apply.f.right.f - angFunc.f.m.tuple + angFunc.f.m end function getExponentPtr(orb::PrimGTOcore) @@ -152,26 +166,26 @@ end function preparePGTOconfig(orb::PrimGTOcore) cenIds = getCenCoordPtr(orb) xpnIdx = getExponentPtr(orb) - ang = getAngularNums(orb) + ang = getAngMomentum(orb) cenIds, xpnIdx, ang end -struct OverlapGTOrbSelf{T, D} <: OrbitalIntegrator{T, D} +struct OverlapGTOrbSelf{T, D, S} <: OrbitalIntegrator{T, D} xpn::FlatPSetInnerPtr{T} - ang::NTuple{D, Int} + ang::WeakComp{D, S} end function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} xpnVal = getField(pars, f.xpn) - overlapPGTO(xpnVal, f.ang) + overlapPGTO(xpnVal, f.ang.tuple) end -struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} +struct OverlapGTOrbPair{T, D, S1, S2} <: OrbitalIntegrator{T, D} cen::NTuple{2, NTuple{ D, FlatPSetInnerPtr{T} }} xpn::NTuple{2, FlatPSetInnerPtr{T}} - ang::NTuple{2, NTuple{D, Int}} + ang::Tuple{WeakComp{D, S1}, WeakComp{D, S2}} end function (f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, From c48a9a171935ad15424509a30f4c691ca263e35d Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 18 Jan 2025 18:52:28 -0500 Subject: [PATCH 11/29] Improved the caching stability. --- src/Integrals/Engines/GaussianOrbitals.jl | 159 +++++++++++----------- 1 file changed, 83 insertions(+), 76 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index dce8d22a..b3165605 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -51,54 +51,67 @@ function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T} end -struct PrimGaussOrbData{T, D, S} +struct PrimGaussOrbInfo{T, D} cen::NTuple{D, T} xpn::T - ang::WeakComp{D, S} + ang::NTuple{D, Int} - PrimGaussOrbData(cen::NTuple{D, T}, xpn::T, ang::WeakComp{D, S}) where {T, D, S} = - new{T, D, S}(cen, xpn, ang) + PrimGaussOrbInfo(cen::NonEmptyTuple{T, D}, xpn::T, + ang::NonEmptyTuple{Int, D}) where {T, D} = + new{T, D+1}(cen, xpn, ang) end +const PrimGaussOrbConfig{T, D} = PrimGaussOrbInfo{FlatPSetInnerPtr{T}, D} -struct GaussProductData{T, D, S1, S2} - lhs::PrimGaussOrbData{T, D, S1} - rhs::PrimGaussOrbData{T, D, S2} + +struct GaussProductInfo{T, D} + lhs::PrimGaussOrbInfo{T, D} + rhs::PrimGaussOrbInfo{T, D} cen::NTuple{D, T} xpn::T - function GaussProductData(oData1::PrimGaussOrbData{T, D, S1}, - oData2::PrimGaussOrbData{T, D, S2}) where {T, D, S1, S2} + function GaussProductInfo(oData1::PrimGaussOrbInfo{T, D}, + oData2::PrimGaussOrbInfo{T, D}) where {T, D} xpnSum = oData1.xpn + oData2.xpn cenNew = (oData1.xpn .* oData1.cen .+ oData2.xpn .* oData2.cen) ./ xpnSum - new{T, D, S1, S2}(oData1, oData2, cenNew, xpnSum) + new{T, D}(oData1, oData2, cenNew, xpnSum) end end +struct NullCache{T} <: QueryBox{T} end + +const TupleOf5T2Int{T} = Tuple{T, T, T, T, Int, Int} -const DefaultPGTOrbOverlapCacheSizeLimit = 100 +const NullOrT5Int2LCU{T} = Union{NullCache{T}, LRU{TupleOf5T2Int{T}, T}} -struct NullCache{T} <: QueryBox{T} end +function overlapPGTOcore(input::TupleOf5T2Int{T}) where {T} + xpnSum, xpnRatio, xML, xMR, iL, iR = input + dxLR = xMR - xML -@generated function gen1DimPGTOrbOverlapCache(::Type{T}, ::Val{S}) where {T, S} - if S==0 - return :( NullCache{T}() ) + if iL == iR == 0 + gaussProdCore4(xpnSum, xpnRatio, dxLR) else - K = Tuple{T, T, T, T, Int, Int} - cache = LRU{K, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit) - return :( $cache ) + jRange = 0:((iL + iR) ÷ 2) + + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) + end + else + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) + end * gaussProdCore3(dxLR, xpnRatio) + end end end -function lazyOverlapPGTO(::NullCache{T}, input::Tuple{T, T, T, T, Int, Int}) where {T} - gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) +function overlapPGTOcore(xpn::T, ang::NonEmptyTuple{Int}) where {T} + mapreduce(*, ang) do i + gaussProdCore1(2xpn, i) + end end -function lazyOverlapPGTO(cache::LRU{K, T}, input::K) where - {T, K<:Tuple{T, T, T, T, Int, Int}} - # get!(cache, input) do - # overlapPGTOcore(input) - # end +function lazyOverlapPGTO(cache::LRU{TupleOf5T2Int{T}, T}, input::TupleOf5T2Int{T}) where {T} res = get(cache, input, nothing) # Fewer allocations than using `get!` if res === nothing res = overlapPGTOcore(input) @@ -107,23 +120,9 @@ function lazyOverlapPGTO(cache::LRU{K, T}, input::K) where res end -function overlapPGTOcore(input::Tuple{T, T, T, T, Int, Int}) where {T} - xpnSum, xpnRatio, xML, xMR, iL, iR = input - jRange = 0:((iL + iR) ÷ 2) - dxLR = xMR - xML +lazyOverlapPGTO(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) - if isequal(dxLR, zero(T)) - mapreduce(+, jRange) do j - gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) - end - else - mapreduce(+, jRange) do j - gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) - end * gaussProdCore3(dxLR, xpnRatio) - end -end - -function overlapPGTO(data::GaussProductData{T, D, S1, S2}) where {T, D, S1, S2} +function overlapPGTO!(cache::NullOrT5Int2LCU{T}, data::GaussProductInfo{T, D}) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -134,16 +133,8 @@ function overlapPGTO(data::GaussProductData{T, D, S1, S2}) where {T, D, S1, S2} angL = data.lhs.ang angR = data.rhs.ang - intCache = gen1DimPGTOrbOverlapCache(T, Val(S1+S2)) - - mapreduce(*, cenL, cenR, cenM, angL.tuple, angR.tuple) do xL, xR, xM, iL, iR - lazyOverlapPGTO(intCache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) - end -end - -function overlapPGTO(xpn::T, ang::NonEmptyTuple{Int}) where {T} - mapreduce(*, ang) do i - gaussProdCore1(2xpn, i) + mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR + lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end @@ -166,48 +157,64 @@ end function preparePGTOconfig(orb::PrimGTOcore) cenIds = getCenCoordPtr(orb) xpnIdx = getExponentPtr(orb) - ang = getAngMomentum(orb) - cenIds, xpnIdx, ang + ang = getAngMomentum(orb).tuple + PrimGaussOrbInfo(cenIds, xpnIdx, ang) end -struct OverlapGTOrbSelf{T, D, S} <: OrbitalIntegrator{T, D} - xpn::FlatPSetInnerPtr{T} - ang::WeakComp{D, S} +struct OverlapGTOrbSelf{T, D} <: OrbitalIntegrator{T, D} + core::PrimGaussOrbConfig{T, D} end function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} - xpnVal = getField(pars, f.xpn) - overlapPGTO(xpnVal, f.ang.tuple) + xpnVal = getField(pars, f.core.xpn) + overlapPGTOcore(xpnVal, f.core.ang) end +struct OverlapGTOrbPairCore{T, D} <: OrbitalIntegrator{T, D} + lhs::PrimGaussOrbConfig{T, D} + rhs::PrimGaussOrbConfig{T, D} +end -struct OverlapGTOrbPair{T, D, S1, S2} <: OrbitalIntegrator{T, D} - cen::NTuple{2, NTuple{ D, FlatPSetInnerPtr{T} }} - xpn::NTuple{2, FlatPSetInnerPtr{T}} - ang::Tuple{WeakComp{D, S1}, WeakComp{D, S2}} +struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} + core::OverlapGTOrbPairCore{T, D} + cache::LRU{TupleOf5T2Int{T}, T} end -function (f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, - pars2::FilteredVecOfArr{T}) where {T} - c1 = getField(pars1, f.cen[begin]) - x1 = getField(pars1, f.xpn[begin]) - d1 = PrimGaussOrbData(c1, x1, f.ang[begin]) - c2 = getField(pars2, f.cen[end]) - x2 = getField(pars2, f.xpn[end]) - d2 = PrimGaussOrbData(c2, x2, f.ang[end]) - GaussProductData(d1, d2) |> overlapPGTO +function (f::OverlapGTOrbPairCore{T})(pars1::FilteredVecOfArr{T}, + pars2::FilteredVecOfArr{T}, + cache::NullOrT5Int2LCU{T}=NullCache{T}()) where {T} + d1, d2 = map((:lhs, :rhs), (pars1, pars2)) do sym, pars + sector = getfield(f, sym) + cen = getField(pars, sector.cen) + xpn = getField(pars, sector.xpn) + PrimGaussOrbInfo(cen, xpn, sector.ang) + end + overlapPGTO!(cache, GaussProductInfo(d1, d2)) end +(f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}) where {T} = +f.core(pars1, pars2, f.cache) + + +const DefaultPGTOrbOverlapCacheSizeLimit = 256 + +@generated function genPGTOrbOverlapCache(::Type{T}, ::Val{L}) where {T, L} + maxsize = min(8(L+2), DefaultPGTOrbOverlapCacheSizeLimit) + cache = LRU{TupleOf5T2Int{T}, T}(;maxsize) + return :( $cache ) +end -function genGTOrbOverlapFunc(orbs::NTuple{2, PrimGTOcore{T, D}}) where {T, D} - configs = preparePGTOconfig.(orbs) - OverlapGTOrbPair(getindex.(configs, 1), getindex.(configs, 2), getindex.(configs, 3)) +function genGTOrbOverlapFunc((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, + TypedPrimGTOcore{T, D, L2}}) where + {T, D, L1, L2} + cache = genPGTOrbOverlapCache(T, Val(L1+L2)) + core = OverlapGTOrbPairCore(preparePGTOconfig(orb1), preparePGTOconfig(orb2)) + OverlapGTOrbPair(core, cache) end function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} - _, xpnIdx, ang = preparePGTOconfig(orb) - OverlapGTOrbSelf(xpnIdx, ang) + preparePGTOconfig(orb) |> OverlapGTOrbSelf end From 89ea75872015c145e566bc213be41c8c952e6d73 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 18 Jan 2025 19:01:08 -0500 Subject: [PATCH 12/29] Added `TypedPrimGTOcore`. --- src/SpatialBasis.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 3d3922c0..2a05b807 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -116,6 +116,9 @@ const PrimATOcore{T, D, B<:EvalPolyGaussProd{T, D}} = const PrimGTOcore{T, D, B<:EvalPolyGaussProd{T, D}} = PrimitiveOrbCore{T, D, B} +const TypedPrimGTOcore{T, D, L} = + PrimitiveOrbCore{T, D, EvalPolyRadialFunc{T, D, EvalGaussFunc{T}, L}} + const EvalPrimGTO{T, D, B<:EvalPolyGaussProd{T, D}, F<:EvalOrbNormalizer{T, D}} = EvalPrimOrb{T, D, B, F} From 6e28895d41dae849e1be3d7f533b6ae64a81e863 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 18 Jan 2025 22:21:14 -0500 Subject: [PATCH 13/29] Optimized `NormalizeCompOrb`. --- src/SpatialBasis.jl | 77 ++++++++++++++++++++++++++------------------- 1 file changed, 45 insertions(+), 32 deletions(-) diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 2a05b807..366a9e08 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -481,57 +481,70 @@ function buildNormalizer(::Type{T}, ::Val{D}) where {T, D} EvalOrbNormalizer(T, Val(D)) end + +function getNormalizerCoreTypeUnion(::Type{T}, ::Val{D}, ::Val{N}, ::Type{F}) where + {T, D, N, F<:OrbitalIntegrator{T, D}} + if isconcretetype(F) + VariedNormCore{T, D, N, F} + else + VariedNormCore{T, D, N, <:F} + end +end + function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} weightedOrbs = f.f.chain nOrbs = length(weightedOrbs) - pOrbs = map(x->x.left, weightedOrbs) + oCorePtr = ChainPointer((:left, :f, :apply, :left)) + oScopePtr = ChainPointer((:left, :f, :scope)) - diagFuncs = Memory{ReturnTyped{T}}(undef, nOrbs) - pOrbCores = Memory{PrimitiveOrbCore{T, D}}(undef, nOrbs) - pOrbScopes = Memory{AwaitFilter{FlatParamSetFilter{T}}}(undef, nOrbs) hasUnit = false diagFuncCoreType = Union{} utriFuncCoreType = Union{} - for i in eachindex(pOrbs) - o = pOrbs[i] - oCore = o.f.apply.left - oScope = first(o.f.scope) - diagFuncs[i] = if isRenormalized(o) + + diagFuncs = map(weightedOrbs) do weightedOrb + if isRenormalized(weightedOrb.left) hasUnit = true ReturnTyped(Unit(T), T) else + oCore = getField(weightedOrb, oCorePtr) + oSelect = getField(weightedOrb, oScopePtr) |> first |> Retrieve nfInner = genOneBodyCoreIntegrator(Identity(), (oCore,)) diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfInner)) - ReturnTyped(EncodeApply((Retrieve(oScope),), nfInner), T) + ReturnTyped(EncodeApply((oSelect,), nfInner), T) end - pOrbCores[i] = oCore - pOrbScopes[i] = oScope end - diagFuncType = if isconcretetype(diagFuncCoreType) - VariedNormCore{T, D, 1, diagFuncCoreType} - else - VariedNormCore{T, D, 1, <:diagFuncCoreType} - end - hasUnit && (diagFuncType = Union{UnitScalarCore{T}, diagFuncType}) - diagFuncs = convert(Memory{diagFuncType}, diagFuncs) - - utriFuncs = Memory{ReturnTyped{T}}(undef, nOrbs*(nOrbs-1)÷2) - for j in 1:length(utriFuncs) - n, m = convert1DidxTo2D(nOrbs-1, j) - corePairs = (pOrbCores[begin+m-1], pOrbCores[begin+n-1]) - nfInner = genOneBodyCoreIntegrator(Identity(), corePairs) - utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfInner)) - filters = (pOrbScopes[begin+m-1], pOrbScopes[begin+n-1]) - utriFuncs[begin+j-1] = ReturnTyped(EncodeApply(Retrieve.(filters), nfInner), T) + if !isconcretetype(eltype(diagFuncs)) + diagFuncType = getNormalizerCoreTypeUnion(T, D, Val(1), diagFuncCoreType) + hasUnit && (diagFuncType = Union{UnitScalarCore{T}, diagFuncType}) + diagFuncs = Memory{diagFuncType}(diagFuncs) end - utriFuncType = if isconcretetype(utriFuncCoreType) - VariedNormCore{T, D, 2, utriFuncCoreType} + if nOrbs > 1 + utriFuncNum = nOrbs * (nOrbs-1) ÷ 2 + utriFuncs = map(1:utriFuncNum) do j + n, m = convert1DidxTo2D(nOrbs-1, j) + weightedOrb1 = weightedOrbs[begin+m-1] + weightedOrb2 = weightedOrbs[begin+n-1] + oCore1 = getField(weightedOrb1, oCorePtr) + oCore2 = getField(weightedOrb2, oCorePtr) + nfInner = genOneBodyCoreIntegrator(Identity(), (oCore1, oCore2)) + utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfInner)) + oSelect1 = getField(weightedOrb1, oScopePtr) |> first |> Retrieve + oSelect2 = getField(weightedOrb2, oScopePtr) |> first |> Retrieve + ReturnTyped(EncodeApply((oSelect1, oSelect2), nfInner), T) + end + + utriFuncType = eltype(utriFuncs) + if !isconcretetype(utriFuncType) + utriFuncType = getNormalizerCoreTypeUnion(T, D, Val(2), utriFuncCoreType) + end + utriFuncs = Memory{utriFuncType}(utriFuncs) else - VariedNormCore{T, D, 2, <:utriFuncCoreType} + oCore = getField(weightedOrbs[begin], oCorePtr) + utriFuncCoreType = (typeof∘genOneBodyCoreIntegrator)(Identity(), (oCore, oCore)) + utriFuncs = Memory{VariedNormCore{T, D, 2, utriFuncCoreType}}(undef, 0) end - utriFuncs = convert(Memory{utriFuncType}, utriFuncs) coreTransform = HermitianContract(diagFuncs, utriFuncs) getWeights = map(b->ReturnTyped(b.right.f, T), weightedOrbs) From 871a9eb9b167a49262ea2e50f3de42d17f6dc937 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 19 Jan 2025 00:59:41 -0500 Subject: [PATCH 14/29] Split cache for different dimensions. --- src/Integrals/Engines/GaussianOrbitals.jl | 30 ++++++++++++----------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index b3165605..d22611eb 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -122,7 +122,8 @@ end lazyOverlapPGTO(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) -function overlapPGTO!(cache::NullOrT5Int2LCU{T}, data::GaussProductInfo{T, D}) where {T, D} +function overlapPGTO!(caches::NTuple{D, NullOrT5Int2LCU{T}}, + data::GaussProductInfo{T, D}) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -133,7 +134,7 @@ function overlapPGTO!(cache::NullOrT5Int2LCU{T}, data::GaussProductInfo{T, D}) w angL = data.lhs.ang angR = data.rhs.ang - mapreduce(*, cenL, cenR, cenM, angL, angR) do xL, xR, xM, iL, iR + mapreduce(*, caches, cenL, cenR, cenM, angL, angR) do cache, xL, xR, xM, iL, iR lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end @@ -178,39 +179,40 @@ end struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} core::OverlapGTOrbPairCore{T, D} - cache::LRU{TupleOf5T2Int{T}, T} + cache::NTuple{D, LRU{TupleOf5T2Int{T}, T}} end -function (f::OverlapGTOrbPairCore{T})(pars1::FilteredVecOfArr{T}, - pars2::FilteredVecOfArr{T}, - cache::NullOrT5Int2LCU{T}=NullCache{T}()) where {T} +function (f::OverlapGTOrbPairCore{T, D})(pars1::FilteredVecOfArr{T}, + pars2::FilteredVecOfArr{T}, + caches::NTuple{D, NullOrT5Int2LCU{T}}= + ntuple( _->NullCache{T}(), Val(D) )) where {T, D} d1, d2 = map((:lhs, :rhs), (pars1, pars2)) do sym, pars sector = getfield(f, sym) cen = getField(pars, sector.cen) xpn = getField(pars, sector.xpn) PrimGaussOrbInfo(cen, xpn, sector.ang) end - overlapPGTO!(cache, GaussProductInfo(d1, d2)) + overlapPGTO!(caches, GaussProductInfo(d1, d2)) end (f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}) where {T} = f.core(pars1, pars2, f.cache) -const DefaultPGTOrbOverlapCacheSizeLimit = 256 +const DefaultPGTOrbOverlapCacheSizeLimit = 128 -@generated function genPGTOrbOverlapCache(::Type{T}, ::Val{L}) where {T, L} - maxsize = min(8(L+2), DefaultPGTOrbOverlapCacheSizeLimit) - cache = LRU{TupleOf5T2Int{T}, T}(;maxsize) - return :( $cache ) +@generated function genPGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} + maxsize = DefaultPGTOrbOverlapCacheSizeLimit + caches = ntuple(_->LRU{TupleOf5T2Int{T}, T}(; maxsize), Val(D)) + return :( $caches ) end function genGTOrbOverlapFunc((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where {T, D, L1, L2} - cache = genPGTOrbOverlapCache(T, Val(L1+L2)) + caches = genPGTOrbOverlapCache(T, Val(D), Val(L1+L2)) core = OverlapGTOrbPairCore(preparePGTOconfig(orb1), preparePGTOconfig(orb2)) - OverlapGTOrbPair(core, cache) + OverlapGTOrbPair(core, caches) end function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} From b43dd148a89041115206546150cff722d4a20dfc Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 19 Jan 2025 14:45:03 -0500 Subject: [PATCH 15/29] Formatted code and replaced anonymous functions with typed functors. --- src/Integrals/Engines/Numerical.jl | 124 +++++++++++++++-------------- 1 file changed, 66 insertions(+), 58 deletions(-) diff --git a/src/Integrals/Engines/Numerical.jl b/src/Integrals/Engines/Numerical.jl index 224b89f4..0ac0a867 100644 --- a/src/Integrals/Engines/Numerical.jl +++ b/src/Integrals/Engines/Numerical.jl @@ -1,90 +1,98 @@ using HCubature -function changeFuncRange(f::F, ::Type{T}) where {F, T} - fCore = function (x) - x = Tuple(x) - f(x ./ (one(T) .- x.^2)) * mapreduce(*, x) do t - (1 + t^2) / (1 - t^2)^2 - end - end - ReturnTyped(fCore, T) +struct ConfineInterval{T, F<:ReturnTyped{T}} <: FunctionModifier + f::F + + ConfineInterval(f::F) where {T, F<:ReturnTyped{T}} = new{T, F}(f) end +ConfineInterval(::Type{T}, f::Function) where {T} = ConfineInterval(ReturnTyped(f, T)) -function numericalIntegration(integrand::ReturnTyped{T}, - interval::NTuple{2, NonEmptyTuple{T, D}}) where {T, D} - (first∘hcubature)(integrand, first(interval), last(interval); maxevals=typemax(Int)) +function (f::ConfineInterval{T})(x::AbstractVector{T}) where {T} + val = f.f(x ./ (one(T) .- x.^2)) + mapreduce(*, x) do t + tSquare = t * t + (1 + tSquare) / (1 - tSquare)^2 + end * val end -function composeOneBodyKernel(op::O, ::Type{T}, termL::F1, termR::F2) where - {O<:Function, T, F1<:Function, F2<:Function} - PairCombine(StableMul(T), adjoint∘termL, op(termR)) -end +struct ModSquaredMap{T, F<:ReturnTyped{T}} <: FunctionModifier + f::F -function composeOneBodyKernel(op::O, ::Type{T}, term::F) where {O<:Function, T, F<:Function} - composeOneBodyKernel(op, T, term, term) + ModSquaredMap(f::F) where {T, F<:ReturnTyped{T}} = new{T, F}(f) end -function composeOneBodyKernel(::Identity, ::Type{T}, term::F) where {T, F<:Function} - function (input::I) where {I} - val = term(input)::T - val' * val - end -end +ModSquaredMap(::Type{T}, f::Function) where {T} = ModSquaredMap(ReturnTyped(f, T)) -function numericalOneBodyInt(op::F, (orb,)::Tuple{EvalFieldFunction{T, D}}, - (pVal,)::Tuple{FilteredVecOfArr{T}}) where - {F<:DirectOperator, T, D} - term = Base.Fix2(orb, pVal) - integrand = changeFuncRange(composeOneBodyKernel(op, T, term), T) - bound = ntuple(_->one(T), Val(D)) - numericalIntegration(integrand, (.-(bound), bound)) +function (f::ModSquaredMap{T})(input::AbstractVector{T}) where {T} + val = f.f(input) + val' * val end -function numericalOneBodyInt(op::F, orbs::NTuple{2, EvalFieldFunction{T, D}}, - pValPair::NTuple{2, FilteredVecOfArr{T}}) where - {F<:DirectOperator, T, D} - termL, termR = Base.Fix2.(orbs, pValPair) - integrand = changeFuncRange(composeOneBodyKernel(op, T, termL, termR), T) - bound = ntuple(_->one(T), Val(D)) - numericalIntegration(integrand, (.-(bound), bound)) + +function numericalIntegrateCore(integrand::ConfineInterval{T}, + interval::NTuple{2, NonEmptyTuple{T, D}}) where {T, D} + fullRes = hcubature(integrand, first(interval), last(interval), maxevals=typemax(Int)) + first(fullRes) end -struct NumOverlapOrbSelf{T, D, B<:PrimitiveOrbCore{T, D}} <: OrbitalIntegrator{T, D} - orb::Tuple{B} +function composeOneBodyKernel(op::O, ::Type{T}, (termL, termR)::Tuple{F1, F2}) where + {O<:Function, T, F1<:Function, F2<:Function} + PairCombine(StableMul(T), adjoint∘termL, op(termR)) end -function (f::NumOverlapOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} - numericalOneBodyInt(Identity(), f.orb, (pars,)) +function composeOneBodyKernel(op::O, ::Type{T}, (term,)::Tuple{F}) where + {O<:Function, T, F<:Function} + composeOneBodyKernel(op, T, (term, term)) end -struct NumOverlapOrbPair{T, D, B1<:PrimitiveOrbCore{T, D}, - B2<:PrimitiveOrbCore{T, D}} <: OrbitalIntegrator{T, D} - orb::Tuple{B1, B2} +function composeOneBodyKernel(::Identity, ::Type{T}, (term,)::Tuple{F}) where + {T, F<:Function} + ModSquaredMap(T, term) end -function (f::NumOverlapOrbPair{T})(pars::Vararg{FilteredVecOfArr{T}, 2}) where {T} - numericalOneBodyInt(Identity(), f.orb, pars) + +function composeIntegralKernel(::OneBodyIntegral, op::O, ::Type{T}, + terms::N12Tuple{Function}) where {O<:Function, T} + composeOneBodyKernel(op, T, terms) end -genOneBodyCoreIntegrator(::Identity, orbs::Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = -NumOverlapOrbSelf(orbs) +function numericalIntegrate(::OneBodyIntegral, op::F, + orbs::NonEmptyTuple{EvalFieldFunction{T, D}, N}, + pVals::NonEmptyTuple{FilteredVecOfArr{T}, N}) where + {F<:DirectOperator, T, D, N} + terms = Base.Fix2.(orbs, pVals) + integralKernel = composeIntegralKernel(OneBodyIntegral(), op, T, terms) + integrand = ConfineInterval(T, integralKernel) + bound = ntuple(_->one(T), Val(D)) + numericalIntegrateCore(integrand, (.-(bound), bound)) +end -genOneBodyCoreIntegrator(::Identity, orbs::NTuple{2, PrimitiveOrbCore{T, D}}) where {T, D} = -NumOverlapOrbPair(orbs) function buildNormalizerCore(o::PrimitiveOrbCore{T, D}) where {T, D} - genOneBodyCoreIntegrator(Identity(), (o,)) + buildOneBodyCoreIntegrator(Identity(), (o,)) end -#!! Implement a functional struct (<:OrbitalIntegrator{T, D}) for the closure -function genOneBodyCoreIntegrator(op::F, - orbs::NonEmptyTuple{PrimitiveOrbCore{T, D}, N}) where - {F<:DirectOperator, N, T, D} - function (pVal::AbtVecOfAbtArr{T}, pVals::Vararg{AbtVecOfAbtArr{T}, N}) where {T} - numericalOneBodyInt(op, orbs, (pVal, pVals...)) - end +struct OneBodyNumIntegrate{T, D, F<:DirectOperator, + P<:N12Tuple{PrimitiveOrbCore{T, D}}} <: OrbitalIntegrator{T, D} + op::F + basis::P +end + +const OneBodySelfNumInt{T, D, F<:DirectOperator, P<:PrimitiveOrbCore{T, D}} = + OneBodyNumIntegrate{T, D, F, Tuple{P}} +const OnyBodyPairNumInt{T, D, F<:DirectOperator, P1<:PrimitiveOrbCore{T, D}, + P2<:PrimitiveOrbCore{T, D}} = + OneBodyNumIntegrate{T, D, F, Tuple{P1, P2}} + +function (f::OneBodySelfNumInt{T})(pVal::FilteredVecOfArr{T}) where {T} + numericalIntegrate(OneBodyIntegral(), f.op, f.basis, (pVal,)) +end + +function (f::OnyBodyPairNumInt{T})(pVal1::FilteredVecOfArr{T}, + pVal2::FilteredVecOfArr{T}) where {T} + numericalIntegrate(OneBodyIntegral(), f.op, f.basis, (pVal1, pVal2)) end \ No newline at end of file From 641718c86b234e46dc8bbfc5e96af0b37b41168c Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 19 Jan 2025 14:46:48 -0500 Subject: [PATCH 16/29] Added support of keyword arguments for `ReturnedTyped`. --- src/Mapping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Mapping.jl b/src/Mapping.jl index 07a4249a..86451c1f 100644 --- a/src/Mapping.jl +++ b/src/Mapping.jl @@ -41,7 +41,7 @@ ReturnTyped(::Type{T}) where {T} = ReturnTyped(itself, T) ReturnTyped(f::ReturnTyped{T}, ::Type{T}) where {T} = itself(f) -(f::ReturnTyped{T, F})(arg...) where {T, F} = convert(T, f.f(arg...)) +(f::ReturnTyped{T, F})(arg...; kws...) where {T, F} = convert(T, f.f(arg...; kws...)) const Return{T} = ReturnTyped{T, ItsType} From dd3933e016fdd05aea0346fafe6926136aa6b1ae Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 19 Jan 2025 14:47:45 -0500 Subject: [PATCH 17/29] Added support of vector input for `ShiftByArg`. --- src/Mapping.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Mapping.jl b/src/Mapping.jl index 86451c1f..22fe0671 100644 --- a/src/Mapping.jl +++ b/src/Mapping.jl @@ -219,10 +219,12 @@ end (f::Plus{T})(arg::T) where {T} = arg + f.val -struct ShiftByArg{T<:Real, D} <: FieldlessFunction end +struct ShiftByArg{T, D} <: FieldlessFunction end -(::ShiftByArg{T, D})(input::NTuple{D, Real}, args::Vararg{T, D}) where {T, D} = -(input .- args) +function (::ShiftByArg{T, D})(input::Union{NTuple{D, T}, AbstractVector{T}}, + args::Vararg{T, D}) where {T, D} + input .- args +end struct HermitianContract{T, F1<:ReturnTyped{T}, F2<:ReturnTyped{T}} <: FunctionComposer From 426fd26e8909b4d4dce8f52178a4b543fc309485 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 19 Jan 2025 14:49:45 -0500 Subject: [PATCH 18/29] Optimized the interface of core integral caching. --- src/Integrals/Engines/GaussianOrbitals.jl | 104 ++++++++++++++-------- src/Integrals/Framework.jl | 55 ++++++++---- src/Types.jl | 4 + 3 files changed, 105 insertions(+), 58 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index d22611eb..e489d194 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -78,30 +78,41 @@ struct GaussProductInfo{T, D} end end -struct NullCache{T} <: QueryBox{T} end + +struct NullCache{T} <: CustomCache{T} end const TupleOf5T2Int{T} = Tuple{T, T, T, T, Int, Int} -const NullOrT5Int2LCU{T} = Union{NullCache{T}, LRU{TupleOf5T2Int{T}, T}} +struct ZeroAngMomCache{T} <: CustomCache{T} end + +const CoreIntCacheBox{T} = Union{CustomCache{T}, LRU{<:Any, T}} + + +struct OneBodyCoreIntCache{T, D, + F<:NTuple{D, CoreIntCacheBox{T}}} <: SpatialIntegralCache{T, D} + axis::F + + OneBodyCoreIntCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = + new{T, D+1, typeof(axis)}(axis) +end + +OneBodyCoreIntCache(::Type{T}, ::Val{D}) where {T, D} = +OneBodyCoreIntCache(ntuple( _->NullCache{T}(), Val(D) )) + function overlapPGTOcore(input::TupleOf5T2Int{T}) where {T} xpnSum, xpnRatio, xML, xMR, iL, iR = input dxLR = xMR - xML + jRange = 0:((iL + iR) ÷ 2) - if iL == iR == 0 - gaussProdCore4(xpnSum, xpnRatio, dxLR) - else - jRange = 0:((iL + iR) ÷ 2) - - if isequal(dxLR, zero(T)) - mapreduce(+, jRange) do j - gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) - end - else - mapreduce(+, jRange) do j - gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) - end * gaussProdCore3(dxLR, xpnRatio) + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) end + else + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) + end * gaussProdCore3(dxLR, xpnRatio) end end @@ -122,8 +133,17 @@ end lazyOverlapPGTO(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) -function overlapPGTO!(caches::NTuple{D, NullOrT5Int2LCU{T}}, - data::GaussProductInfo{T, D}) where {T, D} +lazyOverlapPGTO(::ZeroAngMomCache{T}, input::TupleOf5T2Int{T}) where {T} = +gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) + + +const GTOrbOverlapAxialCache{T} = + Union{NullCache{T}, ZeroAngMomCache{T}, LRU{TupleOf5T2Int{T}, T}} + +const GTOrbOverlapCache{T, D} = + OneBodyCoreIntCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} + +function overlapPGTO!(cache::GTOrbOverlapCache{T, D}, data::GaussProductInfo{T, D}) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -134,7 +154,7 @@ function overlapPGTO!(caches::NTuple{D, NullOrT5Int2LCU{T}}, angL = data.lhs.ang angR = data.rhs.ang - mapreduce(*, caches, cenL, cenR, cenM, angL, angR) do cache, xL, xR, xM, iL, iR + mapreduce(*, cache.axis, cenL, cenR, cenM, angL, angR) do cache, xL, xR, xM, iL, iR lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end @@ -172,51 +192,57 @@ function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} overlapPGTOcore(xpnVal, f.core.ang) end -struct OverlapGTOrbPairCore{T, D} <: OrbitalIntegrator{T, D} +struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} lhs::PrimGaussOrbConfig{T, D} rhs::PrimGaussOrbConfig{T, D} end -struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} - core::OverlapGTOrbPairCore{T, D} - cache::NTuple{D, LRU{TupleOf5T2Int{T}, T}} -end -function (f::OverlapGTOrbPairCore{T, D})(pars1::FilteredVecOfArr{T}, - pars2::FilteredVecOfArr{T}, - caches::NTuple{D, NullOrT5Int2LCU{T}}= - ntuple( _->NullCache{T}(), Val(D) )) where {T, D} +function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, + pars2::FilteredVecOfArr{T}; + cache::GTOrbOverlapCache{T, D}= + OneBodyCoreIntCache(T, Val(D))) where {T, D} d1, d2 = map((:lhs, :rhs), (pars1, pars2)) do sym, pars sector = getfield(f, sym) cen = getField(pars, sector.cen) xpn = getField(pars, sector.xpn) PrimGaussOrbInfo(cen, xpn, sector.ang) end - overlapPGTO!(caches, GaussProductInfo(d1, d2)) + overlapPGTO!(cache, GaussProductInfo(d1, d2)) +end + + +function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} + preparePGTOconfig(orb) |> OverlapGTOrbSelf end -(f::OverlapGTOrbPair{T})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}) where {T} = -f.core(pars1, pars2, f.cache) + +function genGTOrbOverlapFunc((orb1, orb2)::NTuple{2, TypedPrimGTOcore{T, D}}) where {T, D} + OverlapGTOrbPair(preparePGTOconfig(orb1), preparePGTOconfig(orb2)) +end const DefaultPGTOrbOverlapCacheSizeLimit = 128 -@generated function genPGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} +@generated function getGTOrbOverlapCacheCore(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} maxsize = DefaultPGTOrbOverlapCacheSizeLimit caches = ntuple(_->LRU{TupleOf5T2Int{T}, T}(; maxsize), Val(D)) return :( $caches ) end -function genGTOrbOverlapFunc((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, - TypedPrimGTOcore{T, D, L2}}) where - {T, D, L1, L2} - caches = genPGTOrbOverlapCache(T, Val(D), Val(L1+L2)) - core = OverlapGTOrbPairCore(preparePGTOconfig(orb1), preparePGTOconfig(orb2)) - OverlapGTOrbPair(core, caches) +@generated function getGTOrbOverlapCacheCore(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} + caches = ntuple(_->ZeroAngMomCache{T}(), Val(D)) + return :( $caches ) end -function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} - preparePGTOconfig(orb) |> OverlapGTOrbSelf +function getGTOrbOverlapCache((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, + TypedPrimGTOcore{T, D, L2}}) where + {T, D, L1, L2} + getGTOrbOverlapCacheCore(T, Val(D), Val(L1+L2)) |> OneBodyCoreIntCache +end + +function getGTOrbOverlapCache(::Tuple{TypedPrimGTOcore{T}}) where {T} + NullCache{T}() end diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 769ef260..729c5b98 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -1,8 +1,5 @@ using LinearAlgebra: dot -const N12Tuple{T} = Union{Tuple{T}, NTuple{2, T}} -const N24Tuple{T} = Union{NTuple{2, T}, NTuple{4, T}} - const OrbitalBasisSet{T, D} = AbstractVector{<:OrbitalBasis{T, D}} const OrbitalCollection{T, D} = NonEmpTplOrAbtArr{FrameworkOrb{T, D}, 1} const OrbitalInput{T, D} = Union{FrameworkOrb{T, D}, OrbitalCollection{T, D}} @@ -116,10 +113,6 @@ function getPrimCoreIntData(cache::CompleteOneBodyIntegrals, idx::Tuple{Int}) end -genOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -genGTOrbOverlapFunc(orbs) - - isHermitian(::PrimitiveOrbCore{T, D}, ::DirectOperator, ::PrimitiveOrbCore{T, D}) where {T, D} = false @@ -129,17 +122,41 @@ isHermitian(::PrimitiveOrbCore{T, D}, ::Identity, true -function evalOneBodyPrimCoreIntegral(op::DirectOperator, - (orb1, pars1)::OrbCoreData{T, D}, - (orb2, pars2)::OrbCoreData{T, D}) where {T, D} - f = ReturnTyped(genOneBodyCoreIntegrator(op, (orb1, orb2)), T) - f(pars1, pars2) +buildOneBodyCoreIntegrator(op::DirectOperator, + orbs::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = +OneBodyNumIntegrate(op, orbs) + +buildOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = +genGTOrbOverlapFunc(orbs) + + +prepareOneBodyCoreIntCache(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = +getGTOrbOverlapCache(orbs) + +prepareOneBodyCoreIntCache(::Identity, ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = +NullCache{T}() + + +function lazyEvalCoreIntegral(integrator::OrbitalIntegrator{T}, + paramData::N12Tuple{FilteredVecOfArr{T}}, + cache::C) where {T, C} + integrator(paramData...; cache)::T end +function lazyEvalCoreIntegral(integrator::OrbitalIntegrator{T}, + paramData::N12Tuple{FilteredVecOfArr{T}}, + ::NullCache{T}) where {T} + integrator(paramData...)::T +end + + function evalOneBodyPrimCoreIntegral(op::DirectOperator, - (orb, pars)::OrbCoreData{T, D}) where {T, D} - f = ReturnTyped(genOneBodyCoreIntegrator(op, (orb,)), T) - f(pars) + data::N12Tuple{OrbCoreData{T, D}}) where {T, D} + orbData = first.(data) + parData = last.(data) + fCore = buildOneBodyCoreIntegrator(op, orbData) + cache = prepareOneBodyCoreIntCache(op, orbData) + lazyEvalCoreIntegral(fCore, parData, cache) end @@ -167,7 +184,7 @@ sortTensorIndex(arg::Vararg{Int}) = sortTensorIndex(arg) function genOneBodyIntDataPairsCore(op::DirectOperator, (oData,)::Tuple{OrbCoreDataSeq{T}}, (oneBasedIdx,)::Tuple{Int}) where {T} - iiVal = evalOneBodyPrimCoreIntegral(op, oData[begin+oneBasedIdx-1]) + iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],)) (iiVal,) end @@ -177,11 +194,11 @@ function genOneBodyIntDataPairsCore(op::DirectOperator, orbPars1, orbPars2 = map(oDataPair, oneBasedIdxPair) do data, idx data[begin+idx-1] end - ijVal = evalOneBodyPrimCoreIntegral(op, orbPars1, orbPars2) + ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2)) jiVal = if isHermitian(first(orbPars1), op, first(orbPars2)) ijVal' else - evalOneBodyPrimCoreIntegral(op, orbPars2, orbPars1) + evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1)) end (ijVal, jiVal) end @@ -263,7 +280,7 @@ function computeOneBodyPrimCoreIntTensor(op::DirectOperator, len1, len2 = length.(oDataPair) res = ShapedMemory{T}(undef, (len1, len2)) for j in 1:len2, i in 1:len1 - ijVal = evalOneBodyPrimCoreIntegral(op, oData1[begin+i-1], oData2[begin+j-1]) + ijVal = evalOneBodyPrimCoreIntegral(op, (oData1[begin+i-1], oData2[begin+j-1])) res[begin+i-1, begin+j-1] = ijVal end res diff --git a/src/Types.jl b/src/Types.jl index 375ebfa5..1e908441 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -25,6 +25,8 @@ abstract type FunctionalStruct <: CompositeFunction end # A struct with defined abstract type TypedEvaluator{T, F} <: Evaluator{F} end abstract type ViewedObject{T, P} <: QueryBox{T} end +abstract type CustomCache{T} <: QueryBox{T} end +abstract type SpatialIntegralCache{T, D} <: QueryBox{T} end abstract type AbstractAmpTensor{T, O} <: JaggedOperator{T, 0, O} end abstract type AbstractAmplitude{T} <: JaggedOperator{T, 0, 0} end @@ -113,6 +115,8 @@ const FlattenedParam{T, N} = Union{InnerSpanParam{T, N}, OuterSpanParam{T, N}} const AVectorOrNTuple{T, NNMO} = Union{Tuple{T, Vararg{T, NNMO}}, AbstractVector{<:T}} const NonEmptyTuple{T, NMO} = Tuple{T, Vararg{T, NMO}} +const N12Tuple{T} = Union{Tuple{T}, NTuple{2, T}} +const N24Tuple{T} = Union{NTuple{2, T}, NTuple{4, T}} const MissingOr{T} = Union{Missing, T} const AbtArray0D{T} = AbstractArray{T, 0} From fb79d1078d854b2cb9a9dca945c008a3b26dcac2 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 00:25:59 -0500 Subject: [PATCH 19/29] Formatted code and temporarily disabled lazy generations of integral cache. --- src/Integrals/Engines/GaussianOrbitals.jl | 64 +++++++++++------------ src/Integrals/Framework.jl | 11 ++-- src/SpatialBasis.jl | 6 +-- 3 files changed, 41 insertions(+), 40 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index e489d194..00834760 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -46,7 +46,7 @@ function gaussProdCore3(dxLR::T, xpnProdOverSum::T) where {T<:Real} exp(- xpnProdOverSum * dxLR^2) end -function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T} +function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T<:Real} gaussProdCore1(xpnSum, 0) * gaussProdCore3(dxLR, xpnRatio) end @@ -78,6 +78,16 @@ struct GaussProductInfo{T, D} end end +function GaussProductInfo(configs::NTuple{2, PrimGaussOrbConfig{T, D}}, + parsPair::NTuple{2, FilteredVecOfArr{T}}) where {T, D} + d1, d2 = map(configs, parsPair) do sector, pars + cen = getField(pars, sector.cen) + xpn = getField(pars, sector.xpn) + PrimGaussOrbInfo(cen, xpn, sector.ang) + end + GaussProductInfo(d1, d2) +end + struct NullCache{T} <: CustomCache{T} end @@ -143,7 +153,8 @@ const GTOrbOverlapAxialCache{T} = const GTOrbOverlapCache{T, D} = OneBodyCoreIntCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} -function overlapPGTO!(cache::GTOrbOverlapCache{T, D}, data::GaussProductInfo{T, D}) where {T, D} +function overlapPGTO(data::GaussProductInfo{T, D}, + cache::GTOrbOverlapCache{T, D}) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -182,57 +193,42 @@ function preparePGTOconfig(orb::PrimGTOcore) PrimGaussOrbInfo(cenIds, xpnIdx, ang) end - -struct OverlapGTOrbSelf{T, D} <: OrbitalIntegrator{T, D} - core::PrimGaussOrbConfig{T, D} +struct PrimGTOrbOverlap{T, D, B<:N12Tuple{PrimGaussOrbConfig{T, D}} + } <: OrbitalIntegrator{T, D} + basis::B end -function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} - xpnVal = getField(pars, f.core.xpn) - overlapPGTOcore(xpnVal, f.core.ang) -end +const OverlapGTOrbSelf{T, D} = PrimGTOrbOverlap{T, D, Tuple{ PrimGaussOrbConfig{T, D} }} +const OverlapGTOrbPair{T, D} = PrimGTOrbOverlap{T, D, NTuple{ 2, PrimGaussOrbConfig{T, D} }} -struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} - lhs::PrimGaussOrbConfig{T, D} - rhs::PrimGaussOrbConfig{T, D} +function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} + sector = first(f.basis) + xpnVal = getField(pars, sector.xpn) + overlapPGTOcore(xpnVal, sector.ang) end - function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; cache::GTOrbOverlapCache{T, D}= OneBodyCoreIntCache(T, Val(D))) where {T, D} - d1, d2 = map((:lhs, :rhs), (pars1, pars2)) do sym, pars - sector = getfield(f, sym) - cen = getField(pars, sector.cen) - xpn = getField(pars, sector.xpn) - PrimGaussOrbInfo(cen, xpn, sector.ang) - end - overlapPGTO!(cache, GaussProductInfo(d1, d2)) -end - - -function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} - preparePGTOconfig(orb) |> OverlapGTOrbSelf + data = GaussProductInfo(f.basis, (pars1, pars2)) + overlapPGTO(data, cache) end -function genGTOrbOverlapFunc((orb1, orb2)::NTuple{2, TypedPrimGTOcore{T, D}}) where {T, D} - OverlapGTOrbPair(preparePGTOconfig(orb1), preparePGTOconfig(orb2)) +function genGTOrbOverlapFunc(orbCoreData::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + PrimGTOrbOverlap(preparePGTOconfig.(orbCoreData)) end const DefaultPGTOrbOverlapCacheSizeLimit = 128 -@generated function getGTOrbOverlapCacheCore(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} - maxsize = DefaultPGTOrbOverlapCacheSizeLimit - caches = ntuple(_->LRU{TupleOf5T2Int{T}, T}(; maxsize), Val(D)) - return :( $caches ) +function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} + ntuple(_->LRU{TupleOf5T2Int{T}, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit), Val(D)) end -@generated function getGTOrbOverlapCacheCore(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} - caches = ntuple(_->ZeroAngMomCache{T}(), Val(D)) - return :( $caches ) +function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} + ntuple(_->ZeroAngMomCache{T}(), Val(D)) end function getGTOrbOverlapCache((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 729c5b98..2639ced1 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -130,12 +130,17 @@ buildOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where genGTOrbOverlapFunc(orbs) -prepareOneBodyCoreIntCache(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -getGTOrbOverlapCache(orbs) -prepareOneBodyCoreIntCache(::Identity, ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = +prepareOneBodyCoreIntCache(op::DirectOperator, + orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = +genGTO1BCoreIntCache(op, orbs) + +prepareOneBodyCoreIntCache(::DirectOperator, + ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = NullCache{T}() +# adjustOneBodyCoreIntConfig + function lazyEvalCoreIntegral(integrator::OrbitalIntegrator{T}, paramData::N12Tuple{FilteredVecOfArr{T}}, diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 366a9e08..6d644f29 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -508,7 +508,7 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} else oCore = getField(weightedOrb, oCorePtr) oSelect = getField(weightedOrb, oScopePtr) |> first |> Retrieve - nfInner = genOneBodyCoreIntegrator(Identity(), (oCore,)) + nfInner = buildOneBodyCoreIntegrator(Identity(), (oCore,)) diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfInner)) ReturnTyped(EncodeApply((oSelect,), nfInner), T) end @@ -528,7 +528,7 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} weightedOrb2 = weightedOrbs[begin+n-1] oCore1 = getField(weightedOrb1, oCorePtr) oCore2 = getField(weightedOrb2, oCorePtr) - nfInner = genOneBodyCoreIntegrator(Identity(), (oCore1, oCore2)) + nfInner = buildOneBodyCoreIntegrator(Identity(), (oCore1, oCore2)) utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfInner)) oSelect1 = getField(weightedOrb1, oScopePtr) |> first |> Retrieve oSelect2 = getField(weightedOrb2, oScopePtr) |> first |> Retrieve @@ -542,7 +542,7 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} utriFuncs = Memory{utriFuncType}(utriFuncs) else oCore = getField(weightedOrbs[begin], oCorePtr) - utriFuncCoreType = (typeof∘genOneBodyCoreIntegrator)(Identity(), (oCore, oCore)) + utriFuncCoreType = (typeof∘buildOneBodyCoreIntegrator)(Identity(), (oCore, oCore)) utriFuncs = Memory{VariedNormCore{T, D, 2, utriFuncCoreType}}(undef, 0) end From c6c689aa5e37ed340d75ee6b88479ca152d1e666 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 00:27:43 -0500 Subject: [PATCH 20/29] Added `MonomialMul`. --- src/Integrals/Operators.jl | 29 ++++++++++++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/src/Integrals/Operators.jl b/src/Integrals/Operators.jl index 2d238240..0799dbd1 100644 --- a/src/Integrals/Operators.jl +++ b/src/Integrals/Operators.jl @@ -1,3 +1,30 @@ struct Identity <: DirectOperator end -(::Identity)(f::Function) = itself(f) \ No newline at end of file +(::Identity)(f::Function) = itself(f) + + +struct MonomialMul{T, D, L} <: DirectOperator + center::NTuple{D, T} + degree::WeakComp{D, L} +end + +MonomialMul(center::NonEmptyTuple{T, D}, degree::NonEmptyTuple{Int, D}) where {T, D} = +MonomialMul(center, WeakComp(degree)) + +abstract type ConstantOperator end +abstract type ParamBoxOperator end + +function (f::MonomialMul{T, D})(arg::NTuple{D, T}) where {T, D} + mapreduce(StableMul(T), arg, f.center, f.degree) do a, c, d + (a - c)^d + end +end + +function (f::MonomialMul{T, D})(target::F) where {T, D, F<:Function} + PairCombine(StableMul(T), f, target) +end + + +# function transform!(::DimSpanDataCacheBox{T}, ::Identity, orb::FrameworkOrb{T}) where {T} +# itself(orb) +# end \ No newline at end of file From c6728efff5df8c09dbd5936fef93073017572ce4 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 00:29:20 -0500 Subject: [PATCH 21/29] Implemented integral computations of multipole moments. --- src/Integrals/Engines/GaussianOrbitals.jl | 104 ++++++++++++++++++++-- src/Integrals/Framework.jl | 3 + src/Integrals/Interface.jl | 20 ++++- 3 files changed, 118 insertions(+), 9 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 00834760..86e03afc 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -231,17 +231,105 @@ function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} ntuple(_->ZeroAngMomCache{T}(), Val(D)) end -function getGTOrbOverlapCache((orb1, orb2)::Tuple{TypedPrimGTOcore{T, D, L1}, - TypedPrimGTOcore{T, D, L2}}) where - {T, D, L1, L2} - getGTOrbOverlapCacheCore(T, Val(D), Val(L1+L2)) |> OneBodyCoreIntCache -end +genGTO1BCoreIntCache(::MonomialMul{T, D, L}, + ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where + {T, D, L, L1, L2} = +genPrimGTOrbOverlapCache(T, Val(D), Val(L+L1+L2)) |> OneBodyCoreIntCache -function getGTOrbOverlapCache(::Tuple{TypedPrimGTOcore{T}}) where {T} - NullCache{T}() -end +genGTO1BCoreIntCache(::Identity, + ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where + {T, D, L1, L2} = +genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> OneBodyCoreIntCache + +genGTO1BCoreIntCache(::DirectOperator, ::Tuple{TypedPrimGTOcore{T}}) where {T} = +NullCache{T}() + +genGTO1BCoreIntCache(op::MonomialMul{T}, (orb,)::Tuple{TypedPrimGTOcore{T}}) where {T} = +genGTO1BCoreIntCache(op, (orb, orb)) function buildNormalizerCore(o::PrimGTOcore{T, D}) where {T, D} genGTOrbOverlapFunc((o,)) +end + + +## Multiple Moment ## +function computeMultiMomentGTO(op::MonomialMul{T, D}, + data::PrimGaussOrbInfo{T, D}) where {T, D} + xpn = data.xpn + mapreduce(*, op.center, op.degree.tuple, data.center, data.ang) do xMM, n, x, i + dx = x - xMM + m = iszero(dx) ? 0 : n + mapreduce(+, 0:m) do k + l = n - k + if isodd(l) + zero(T) + else + binomial(n, k) * dx^k * gaussProdCore1(2xpn, i + (l >> 1)) + end + end + end +end + + +function computeMultiMomentGTO(op::MonomialMul{T, D}, data::GaussProductInfo{T, D}, + cache::GTOrbOverlapCache{T, D}) where {T, D} + + cenL = data.lhs.cen + cenR = data.rhs.cen + cenM = data.cen + + xpnS = data.xpn + xpnRatio = data.lhs.xpn * data.rhs.xpn / xpnS + + angL = data.lhs.ang + angR = data.rhs.ang + + mapreduce(*, cache.axis, op.center, op.degree.tuple, + cenL, cenR, cenM, angL, angR) do cache, xMM, n, xL, xR, xM, iL, iR + dx = xR - xMM #! Consider when xL == xMM + m = iszero(dx) ? 0 : n + mapreduce(+, 0:m) do k + binomial(n, k) * dx^k * + lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR+n-k)) + end + end +end + +computeMultiMomentGTO(::MonomialMul{T, D, 0}, data::GaussProductInfo{T, D}, + cache::GTOrbOverlapCache{T, D}) where {T, D} = +overlapPGTO(data, cache) + + +struct PrimGTOrbMultiMoment{T, D, L, B<:N12Tuple{PrimGaussOrbConfig{T, D}} + } <: OrbitalIntegrator{T, D} + op::MonomialMul{T, D, L} + basis::B +end + +const MultiMomentGTOrbSelf{T, D, L} = + PrimGTOrbMultiMoment{T, D, L, Tuple{ PrimGaussOrbConfig{T, D} }} +const MultiMomentGTOrbPair{T, D, L} = + PrimGTOrbMultiMoment{T, D, L, NTuple{ 2, PrimGaussOrbConfig{T, D} }} + +function (f::MultiMomentGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} + sector = first(f.basis) + cen = getField(pars, sector.cen) + xpn = getField(pars, sector.xpn) + data = PrimGaussOrbInfo(cen, xpn, sector.ang) + computeMultiMomentGTO(f.op, data) +end + +function (f::MultiMomentGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, + pars2::FilteredVecOfArr{T}; + cache::GTOrbOverlapCache{T, D}= + OneBodyCoreIntCache(T, Val(D))) where {T, D} + data = GaussProductInfo(f.basis, (pars1, pars2)) + computeMultiMomentGTO(f.op, data, cache) +end + + +function genGTOrbMultiMomentFunc(op::MonomialMul{T, D}, + orbCoreData::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + PrimGTOrbMultiMoment(op, preparePGTOconfig.(orbCoreData)) end \ No newline at end of file diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 2639ced1..37d4cbf3 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -129,6 +129,9 @@ OneBodyNumIntegrate(op, orbs) buildOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = genGTOrbOverlapFunc(orbs) +buildOneBodyCoreIntegrator(op::MonomialMul{T, D}, + orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = +genGTOrbMultiMomentFunc(op, orbs) prepareOneBodyCoreIntCache(op::DirectOperator, diff --git a/src/Integrals/Interface.jl b/src/Integrals/Interface.jl index 8c3b051b..ec837bc2 100644 --- a/src/Integrals/Interface.jl +++ b/src/Integrals/Interface.jl @@ -11,4 +11,22 @@ end overlaps(basisSet::OrbitalBasisSet{T}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T)) where {T} = -computeIntTensor(OneBodyIntegral(), Identity(), basisSet; paramCache) \ No newline at end of file +computeIntTensor(OneBodyIntegral(), Identity(), basisSet; paramCache) + + +function multipoleMoment(center::NTuple{D, Real}, degrees::NTuple{D, Int}, + orb1::OrbitalBasis{T, D}, orb2::OrbitalBasis{T, D}; + paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), + lazyCompute::Bool=false) where {T, D} + mmOp = MonomialMul(T.(center), degrees) + orbData = orb1 === orb2 ? (orb1,) : (orb1, orb2) + computeIntegral(OneBodyIntegral(), mmOp, (orb1, orb2); paramCache, lazyCompute) +end + +function multipoleMoments(center::NTuple{D, Real}, degrees::NTuple{D, Int}, + basisSet::OrbitalBasisSet{T, D}; + paramCache::DimSpanDataCacheBox{T}= + DimSpanDataCacheBox(T)) where {T, D} + mmOp = MonomialMul(T.(center), degrees) + computeIntTensor(OneBodyIntegral(), mmOp, basisSet; paramCache) +end \ No newline at end of file From 867575862f581345d7575728a318b25f03fc8813 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 18:18:23 -0500 Subject: [PATCH 22/29] Formatting code. --- src/Integrals/Engines/GaussianOrbitals.jl | 30 +++++++++-------------- src/Traits.jl | 6 ++--- src/Types.jl | 5 +++- 3 files changed, 19 insertions(+), 22 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 86e03afc..6800d4b4 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -98,16 +98,16 @@ struct ZeroAngMomCache{T} <: CustomCache{T} end const CoreIntCacheBox{T} = Union{CustomCache{T}, LRU{<:Any, T}} -struct OneBodyCoreIntCache{T, D, - F<:NTuple{D, CoreIntCacheBox{T}}} <: SpatialIntegralCache{T, D} +struct OneBodyIntProcessCache{T, D, F<:NTuple{D, CoreIntCacheBox{T}} + } <: IntegralProcessCache{T, D} axis::F - OneBodyCoreIntCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = + OneBodyIntProcessCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = new{T, D+1, typeof(axis)}(axis) end -OneBodyCoreIntCache(::Type{T}, ::Val{D}) where {T, D} = -OneBodyCoreIntCache(ntuple( _->NullCache{T}(), Val(D) )) +OneBodyIntProcessCache(::Type{T}, ::Val{D}) where {T, D} = +OneBodyIntProcessCache(ntuple( _->NullCache{T}(), Val(D) )) function overlapPGTOcore(input::TupleOf5T2Int{T}) where {T} @@ -151,7 +151,7 @@ const GTOrbOverlapAxialCache{T} = Union{NullCache{T}, ZeroAngMomCache{T}, LRU{TupleOf5T2Int{T}, T}} const GTOrbOverlapCache{T, D} = - OneBodyCoreIntCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} + OneBodyIntProcessCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} function overlapPGTO(data::GaussProductInfo{T, D}, cache::GTOrbOverlapCache{T, D}) where {T, D} @@ -210,7 +210,7 @@ end function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; cache::GTOrbOverlapCache{T, D}= - OneBodyCoreIntCache(T, Val(D))) where {T, D} + OneBodyIntProcessCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) overlapPGTO(data, cache) end @@ -231,21 +231,15 @@ function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} ntuple(_->ZeroAngMomCache{T}(), Val(D)) end -genGTO1BCoreIntCache(::MonomialMul{T, D, L}, +genGTOrbIntCompCache(::MonomialMul{T, D, L}, ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where {T, D, L, L1, L2} = -genPrimGTOrbOverlapCache(T, Val(D), Val(L+L1+L2)) |> OneBodyCoreIntCache +genPrimGTOrbOverlapCache(T, Val(D), Val(L+L1+L2)) |> OneBodyIntProcessCache -genGTO1BCoreIntCache(::Identity, +genGTOrbIntCompCache(::Identity, ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where {T, D, L1, L2} = -genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> OneBodyCoreIntCache - -genGTO1BCoreIntCache(::DirectOperator, ::Tuple{TypedPrimGTOcore{T}}) where {T} = -NullCache{T}() - -genGTO1BCoreIntCache(op::MonomialMul{T}, (orb,)::Tuple{TypedPrimGTOcore{T}}) where {T} = -genGTO1BCoreIntCache(op, (orb, orb)) +genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> OneBodyIntProcessCache function buildNormalizerCore(o::PrimGTOcore{T, D}) where {T, D} @@ -323,7 +317,7 @@ end function (f::MultiMomentGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; cache::GTOrbOverlapCache{T, D}= - OneBodyCoreIntCache(T, Val(D))) where {T, D} + OneBodyIntProcessCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) computeMultiMomentGTO(f.op, data, cache) end diff --git a/src/Traits.jl b/src/Traits.jl index dff862a0..edf88e80 100644 --- a/src/Traits.jl +++ b/src/Traits.jl @@ -55,8 +55,8 @@ returnUndefinedTraitError(InputStyle, F) abstract type IntegralStyle <: AnyInterface end -abstract type MultiBodyIntegral <: IntegralStyle end +abstract type MultiBodyIntegral{D} <: IntegralStyle end -struct OneBodyIntegral <: MultiBodyIntegral end +struct OneBodyIntegral{D} <: MultiBodyIntegral{D} end -struct TwoBodyIntegral <: MultiBodyIntegral end \ No newline at end of file +struct TwoBodyIntegral{D} <: MultiBodyIntegral{D} end \ No newline at end of file diff --git a/src/Types.jl b/src/Types.jl index 1e908441..3d288a86 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -24,9 +24,10 @@ abstract type FunctionalStruct <: CompositeFunction end # A struct with defined abstract type TypedEvaluator{T, F} <: Evaluator{F} end +abstract type SpatialIntegralCache{T, D} <: QueryBox{T} end +abstract type IntegralData{T, S} <: QueryBox{T} end abstract type ViewedObject{T, P} <: QueryBox{T} end abstract type CustomCache{T} <: QueryBox{T} end -abstract type SpatialIntegralCache{T, D} <: QueryBox{T} end abstract type AbstractAmpTensor{T, O} <: JaggedOperator{T, 0, O} end abstract type AbstractAmplitude{T} <: JaggedOperator{T, 0, 0} end @@ -42,6 +43,8 @@ abstract type JoinedOperator{J} <: FunctionComposer end abstract type CompositePointer <: ConfigBox end abstract type StructuredType <: ConfigBox end +abstract type IntegralProcessCache{T, D} <: SpatialIntegralCache{T, D} end + abstract type ActivePointer <: CompositePointer end abstract type StaticPointer{P<:ActivePointer} <: CompositePointer end From 5a1e60e79abe28a8324342845f45c8b895b19fd6 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 18:19:38 -0500 Subject: [PATCH 23/29] Implemented interfaces for integral computation cache. --- src/Integrals/Framework.jl | 292 ++++++++++++++++++++++--------------- 1 file changed, 173 insertions(+), 119 deletions(-) diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 37d4cbf3..d0646812 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -59,56 +59,54 @@ struct IndexedWeights{T} <: QueryBox{T} end -abstract type IntegralData{T, S} <: QueryBox{T} end - -struct CompleteOneBodyIntegrals{T<:Number} <: IntegralData{T, OneBodyIntegral} +struct OneBodyFullCoreIntegrals{T<:Number, D} <: IntegralData{T, OneBodyIntegral{D}} aa::Dict{ Tuple{ Int}, Tuple{ T}} ab::Dict{NTuple{2, Int}, NTuple{2, T}} -end -CompleteOneBodyIntegrals(::Type{T}) where {T} = -CompleteOneBodyIntegrals(Dict{ Tuple{ Int}, Tuple{ T}}(), - Dict{NTuple{2, Int}, NTuple{2, T}}()) - -struct OneBodySymmetricDiffData{T<:Number} <: IntegralData{T, OneBodyIntegral} - ab::Dict{NTuple{2, Int}, T} + OneBodyFullCoreIntegrals(::Type{T}, ::Val{D}) where {T, D} = + new{T, D}(Dict{ Tuple{Int}, Tuple{T}}(), Dict{NTuple{2, Int}, NTuple{2, T}}()) end +# struct OneBodySymmetricDiffData{T<:Number} <: IntegralData{T, OneBodyIntegral} +# ab::Dict{NTuple{2, Int}, T} +# end + -struct IntegralCache{T, D, F<:DirectOperator, B<:PrimitiveOrbCore{T, D}, - I<:IntegralData{T}} <: QueryBox{T} +struct PrimOrbCoreIntegralCache{T, D, F<:DirectOperator, + I<:IntegralData{T, <:MultiBodyIntegral{D}}, + B<:PrimitiveOrbCore{T, D}} <: SpatialIntegralCache{T, D} operator::F - basis::PrimOrbCoreCache{T, D, B} data::I + basis::PrimOrbCoreCache{T, D, B} end -const OneBodyIntCache{T, D, F<:DirectOperator, B<:PrimitiveOrbCore{T, D}, - I<:IntegralData{T, OneBodyIntegral}} = - IntegralCache{T, D, F, B, I} +const OneBodyCoreIntCache{T, D, F<:DirectOperator, I<:IntegralData{T, OneBodyIntegral{D}}, + B<:PrimitiveOrbCore{T, D}} = + PrimOrbCoreIntegralCache{T, D, F, I, B} -const OverlapCache{T, D, B<:PrimitiveOrbCore{T, D}, - I<:IntegralData{T, OneBodyIntegral}} = - OneBodyIntCache{T, D, Identity, B, I} +const OverlapCoreCache{T, D, I<:IntegralData{T, OneBodyIntegral{D}}, + B<:PrimitiveOrbCore{T, D}} = + OneBodyCoreIntCache{T, D, Identity, I, B} -function setIntegralData!(ints::CompleteOneBodyIntegrals{T}, +function setIntegralData!(ints::OneBodyFullCoreIntegrals{T}, pair::Pair{Tuple{Int}, Tuple{T}}) where {T} setindex!(getfield(ints, OneBodyIdxSymDict[true ]), pair.second, pair.first) ints end -function setIntegralData!(ints::CompleteOneBodyIntegrals{T}, +function setIntegralData!(ints::OneBodyFullCoreIntegrals{T}, pair::Pair{NTuple{2, Int}, NTuple{2, T}}) where {T} setindex!(getfield(ints, OneBodyIdxSymDict[false]), pair.second, pair.first) ints end -function getPrimCoreIntData(cache::CompleteOneBodyIntegrals, idx::NTuple{2, Int}) +function getPrimCoreIntData(cache::OneBodyFullCoreIntegrals, idx::NTuple{2, Int}) getfield(cache, OneBodyIdxSymDict[false])[idx] end -function getPrimCoreIntData(cache::CompleteOneBodyIntegrals, idx::Tuple{Int}) +function getPrimCoreIntData(cache::OneBodyFullCoreIntegrals, idx::Tuple{Int}) getfield(cache, OneBodyIdxSymDict[true ])[idx] end @@ -133,38 +131,61 @@ buildOneBodyCoreIntegrator(op::MonomialMul{T, D}, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = genGTOrbMultiMomentFunc(op, orbs) +const OneBodyCoreIntData{T, D} = Tuple{DirectOperator, Vararg{PrimitiveOrbCore{T, D}, 2}} -prepareOneBodyCoreIntCache(op::DirectOperator, - orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -genGTO1BCoreIntCache(op, orbs) +const OneBodyCICacheDict{T, D} = + Dict{Type{<:OneBodyCoreIntData{T, D}}, OneBodyIntProcessCache{T, D}} +#! Consider a new type union when two-body computation cache is needed +const IntCompCacheDict{T, D} = Union{OneBodyCICacheDict{T, D}, TypedEmptyDict{Val{D}, T}} + +IntCompCacheDict{T, D}() where {T, D} = TypedEmptyDict{Val{D}, T}() -prepareOneBodyCoreIntCache(::DirectOperator, - ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = -NullCache{T}() +IntCompCacheDict(::Type{T}, ::Type{OneBodyIntegral{D}}) where {T, D} = +OneBodyCICacheDict{T, D}() + +IntCompCacheDict(::Type{T}, ::Type{TwoBodyIntegral{D}}) where {T, D} = +TypedEmptyDict{Val{D}, T}() + +IntCompCacheDict(::Type{<:IntegralData{T, S}}) where {T, S<:MultiBodyIntegral} = +IntCompCacheDict(T, S) + +function prepareOneBodyIntCompCache!(cacheDict::OneBodyCICacheDict{T, D}, + op::DirectOperator, + orbs::NTuple{2, PrimGTOcore{T, D}}) where {T, D} + get!(cacheDict, typeof( (op, orbs...) )) do + genGTOrbIntCompCache(op, orbs) + end +end + +function prepareOneBodyIntCompCache!(::IntCompCacheDict{T, D}, ::DirectOperator, + ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} + NullCache{T}() +end # adjustOneBodyCoreIntConfig -function lazyEvalCoreIntegral(integrator::OrbitalIntegrator{T}, - paramData::N12Tuple{FilteredVecOfArr{T}}, - cache::C) where {T, C} +function lazyEvalCoreIntegral!(cache::IntegralProcessCache{T}, + integrator::OrbitalIntegrator{T}, + paramData::N12Tuple{FilteredVecOfArr{T}}) where {T} integrator(paramData...; cache)::T end -function lazyEvalCoreIntegral(integrator::OrbitalIntegrator{T}, - paramData::N12Tuple{FilteredVecOfArr{T}}, - ::NullCache{T}) where {T} +function lazyEvalCoreIntegral!(::NullCache{T}, integrator::OrbitalIntegrator{T}, + paramData::N12Tuple{FilteredVecOfArr{T}}) where {T} integrator(paramData...)::T end function evalOneBodyPrimCoreIntegral(op::DirectOperator, - data::N12Tuple{OrbCoreData{T, D}}) where {T, D} + data::N12Tuple{OrbCoreData{T, D}}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} orbData = first.(data) parData = last.(data) fCore = buildOneBodyCoreIntegrator(op, orbData) - cache = prepareOneBodyCoreIntCache(op, orbData) - lazyEvalCoreIntegral(fCore, parData, cache) + intCompCache = prepareOneBodyIntCompCache!(computeCache, op, orbData) + lazyEvalCoreIntegral!(intCompCache, fCore, parData) end @@ -190,23 +211,27 @@ sortTensorIndex(arg::Vararg{Int}) = sortTensorIndex(arg) function genOneBodyIntDataPairsCore(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T}}, - (oneBasedIdx,)::Tuple{Int}) where {T} - iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],)) + (oData,)::Tuple{OrbCoreDataSeq{T, D}}, + (oneBasedIdx,)::Tuple{Int}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} + iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],); computeCache) (iiVal,) end function genOneBodyIntDataPairsCore(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T}}, - oneBasedIdxPair::NTuple{2, Int}) where {T} + oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, + oneBasedIdxPair::NTuple{2, Int}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} orbPars1, orbPars2 = map(oDataPair, oneBasedIdxPair) do data, idx data[begin+idx-1] end - ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2)) + ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2); computeCache) jiVal = if isHermitian(first(orbPars1), op, first(orbPars2)) ijVal' else - evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1)) + evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1); computeCache) end (ijVal, jiVal) end @@ -214,19 +239,21 @@ end function genOneBodyIntDataPairs(op::DirectOperator, (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - (indexOffset,)::Tuple{Int}=(0,)) where {T, D} + (indexOffset,)::Tuple{Int}=(0,); + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} nOrbs = length(oData) offset = indexOffset + firstindex(oData) - 1 pairs1 = map(1:nOrbs) do i - iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,)) + iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,); computeCache) (i + offset,) => iiVal end pairs2 = map(1:triMatEleNum(nOrbs-1)) do l n, m = convert1DidxTo2D(nOrbs-1, l) ijPair = sortTensorIndex((m, n+1)) - ijValPair = genOneBodyIntDataPairsCore(op, (oData, oData), ijPair) + ijValPair = genOneBodyIntDataPairsCore(op, (oData, oData), ijPair; computeCache) (ijPair .+ offset) => ijValPair end @@ -235,7 +262,9 @@ end function genOneBodyIntDataPairs(op::DirectOperator, oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - indexOffsets::NTuple{2, Int}) where {T, D} + indexOffsets::NTuple{2, Int}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} mnOrbs = length.(oDataPair) offsets = indexOffsets .+ firstindex.(oDataPair) .- 1 @@ -243,7 +272,7 @@ function genOneBodyIntDataPairs(op::DirectOperator, idxPairOld = mnIdx .+ offsets idxPairNew = sortTensorIndex(idxPairOld) ijPair = ifelse(idxPairNew == idxPairOld, mnIdx, reverse(mnIdx)) - ijValPair = genOneBodyIntDataPairsCore(op, oDataPair, ijPair) + ijValPair = genOneBodyIntDataPairsCore(op, oDataPair, ijPair; computeCache) idxPairNew => ijValPair end |> vec end @@ -264,38 +293,44 @@ end function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T}}) where {T} + (oData,)::Tuple{OrbCoreDataSeq{T, D}}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}() + ) where {T, D} nOrbs = length(oData) res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) for i in 1:nOrbs - temp = genOneBodyIntDataPairsCore(op, (oData,), (i,)) + temp = genOneBodyIntDataPairsCore(op, (oData,), (i,); computeCache) setTensorEntries!(res, temp, (i,)) end for l in 1:triMatEleNum(nOrbs-1) n, m = convert1DidxTo2D(mBasis-1, l) - temp = genOneBodyIntDataPairsCore(op, (oData, oData), (m, n+1)) + temp = genOneBodyIntDataPairsCore(op, (oData, oData), (m, n+1); computeCache) setTensorEntries!(res, temp, (m, n+1)) end res end function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}) where - {T, D} + oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}() + ) where {T, D} oData1, oData2 = oDataPair len1, len2 = length.(oDataPair) res = ShapedMemory{T}(undef, (len1, len2)) for j in 1:len2, i in 1:len1 - ijVal = evalOneBodyPrimCoreIntegral(op, (oData1[begin+i-1], oData2[begin+j-1])) + oCorePair = (oData1[begin+i-1], oData2[begin+j-1]) + ijVal = evalOneBodyPrimCoreIntegral(op, oCorePair; computeCache) res[begin+i-1, begin+j-1] = ijVal end res end -struct BasisIndexList +struct BasisIndexList <: QueryBox{Int} index::Memory{Int} endpoint::Memory{Int} @@ -405,80 +440,94 @@ function indexCacheOrbData!(targetCache::PrimOrbCoreCache{T, D}, end -function cachePrimCoreIntegrals!(intCache::IntegralCache{T, D}, +function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, paramCache::DimSpanDataCacheBox{T}, orbMarkerCache::OrbCoreMarkerDict, - orbs::OrbitalCollection{T, D}) where {T, D} + orbs::OrbitalCollection{T, D}) where {T, D, + F<:DirectOperator, + I<:IntegralData{ T, <:MultiBodyIntegral{D} }} orbCache = intCache.basis + computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(orbCache.list) orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbs) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1) + updatePrimCoreIntCache!(intCache, oldMaxIdx+1; computeCache) orbIdxList end -function cachePrimCoreIntegrals!(intCache::IntegralCache{T, D}, +function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, paramCache::DimSpanDataCacheBox{T}, orbMarkerCache::OrbCoreMarkerDict, - orb::FrameworkOrb{T, D}) where {T, D} + orb::FrameworkOrb{T, D}) where {T, D, F<:DirectOperator, + I<:IntegralData{ T, <:MultiBodyIntegral{D} }} orbCache = intCache.basis + computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(orbCache.list) orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orb) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1) + updatePrimCoreIntCache!(intCache, oldMaxIdx+1; computeCache) orbIdxList end -function cachePrimCoreIntegrals!(targetIntCache::IntegralCache{T, D}, - sourceIntCache::IntegralCache{T, D}, +function cachePrimCoreIntegrals!(targetIntCache::PrimOrbCoreIntegralCache{T, D, F, I}, + sourceIntCache::PrimOrbCoreIntegralCache{T, D, F, I}, orbMarkerCache::OrbCoreMarkerDict, - sourceOrbList::BasisIndexList) where {T, D} + sourceOrbList::BasisIndexList) where {T, D, + F<:DirectOperator, + I<:IntegralData{ T, <:MultiBodyIntegral{D} }} tOrbCache = targetIntCache.basis + computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(tOrbCache.list) sOrbCache = sourceIntCache.basis orbIdxList = indexCacheOrbData!(tOrbCache, sOrbCache, orbMarkerCache, sourceOrbList) - updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1) + updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1; computeCache) orbIdxList end -function updateIntCacheCore!(op::DirectOperator, ints::CompleteOneBodyIntegrals{T}, +function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, basis::Tuple{OrbCoreDataSeq{T, D}}, - offset::Tuple{Int}) where {T, D} - pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset) + offset::Tuple{Int}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} + pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset; computeCache) foreach(p->setIntegralData!(ints, p), pairs1) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updateIntCacheCore!(op::DirectOperator, ints::CompleteOneBodyIntegrals{T}, +function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, basis::NTuple{2, OrbCoreDataSeq{T, D}}, - offset::NTuple{2, Int}) where {T, D} - pairs2 = genOneBodyIntDataPairs(op, basis, offset) + offset::NTuple{2, Int}; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} + pairs2 = genOneBodyIntDataPairs(op, basis, offset; computeCache) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updatePrimCoreIntCache!(cache::IntegralCache, startIdx::Int) +function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D}, startIdx::Int; + computeCache::IntCompCacheDict{T, D}= + IntCompCacheDict{T, D}()) where {T, D} op = cache.operator basis = cache.basis.list ints = cache.data firstIdx = firstindex(basis) if startIdx == firstIdx - updateIntCacheCore!(op, ints, (basis,), (0,)) + updateIntCacheCore!(op, ints, (basis,), (0,); computeCache) elseif firstIdx < startIdx <= lastindex(basis) boundary = startIdx - 1 oldBasis = @view basis[begin:boundary] newBasis = @view basis[startIdx: end] - updateIntCacheCore!(op, ints, (newBasis,), (boundary,)) - updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary)) + updateIntCacheCore!(op, ints, (newBasis,), (boundary,); computeCache) + updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary); computeCache) end cache end -function decodePrimCoreInt(cache::CompleteOneBodyIntegrals{T}, ptrPair::NTuple{2, Int}, +function decodePrimCoreInt(cache::OneBodyFullCoreIntegrals{T}, ptrPair::NTuple{2, Int}, (coeff1, coeff2)::NTuple{2, T}=( one(T), one(T) )) where {T} coeffProd = coeff1' * coeff2 ptrPairNew = sortTensorIndex(ptrPair) @@ -486,13 +535,13 @@ function decodePrimCoreInt(cache::CompleteOneBodyIntegrals{T}, ptrPair::NTuple{2 (ptrPairNew == ptrPair ? data : reverse(data)) .* (coeffProd, coeffProd') end -function decodePrimCoreInt(cache::CompleteOneBodyIntegrals{T}, ptr::Tuple{Int}, +function decodePrimCoreInt(cache::OneBodyFullCoreIntegrals{T}, ptr::Tuple{Int}, (coeff1,)::Tuple{T}=( one(T), )) where {T} getPrimCoreIntData(cache, ptr) .* coeff1' .* coeff1 end -function buildPrimOrbWeight(normCache::OverlapCache{T, D}, orb::EvalPrimOrb{T, D}, +function buildPrimOrbWeight(normCache::OverlapCoreCache{T, D}, orb::EvalPrimOrb{T, D}, idx::Int) where {T, D} if isRenormalized(orb) decodePrimCoreInt(normCache.data, (idx,)) |> first |> AbsSqrtInv @@ -503,7 +552,8 @@ end function buildNormalizedCompOrbWeight!(weight::AbstractVector{T}, - normCache::OverlapCache{T, D}, orb::FCompOrb{T, D}, + normCache::OverlapCoreCache{T, D}, + orb::FCompOrb{T, D}, idxSeq::AbstractVector{Int}) where {T, D} overlapCache = normCache.data nPrimOrbs = length(weight) @@ -533,7 +583,7 @@ end function buildOrbWeight!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCache{T, D}, orb::FrameworkOrb{T, D}, + normCache::OverlapCoreCache{T, D}, orb::FrameworkOrb{T, D}, idxSeq::AbstractVector{Int}) where {T, D} nPrimOrbs = orbSizeOf(orb) weight = Memory{T}(undef, nPrimOrbs) @@ -568,7 +618,7 @@ end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCache{T, D}, + normCache::OverlapCoreCache{T, D}, orbs::OrbitalCollection{T, D}, normIdxList::BasisIndexList, intIdxList::BasisIndexList=normIdxList) where {T, D} @@ -583,7 +633,7 @@ function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCache{T, D}, + normCache::OverlapCoreCache{T, D}, orb::FrameworkOrb{T, D}, normIdxList::BasisIndexList, intIdxList::BasisIndexList=normIdxList) where {T, D} @@ -592,8 +642,8 @@ function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, end -function prepareIntegralConfig!(intCache::IntegralCache{T, D}, - normCache::OverlapCache{T, D}, +function prepareIntegralConfig!(intCache::PrimOrbCoreIntegralCache{T, D}, + normCache::OverlapCoreCache{T, D}, paramCache::DimSpanDataCacheBox{T}, orbMarkerCache::OrbCoreMarkerDict, orbInput::OrbitalInput{T, D}) where {T, D} @@ -607,7 +657,7 @@ function prepareIntegralConfig!(intCache::IntegralCache{T, D}, end -function buildIntegralEntries(intCache::OneBodyIntCache{T}, +function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, (intWeights,)::Tuple{IndexedWeights{T}}) where {T} idxList = intWeights.list len = length(idxList) @@ -625,7 +675,7 @@ function buildIntegralEntries(intCache::OneBodyIntCache{T}, (res,) # ([1|O|1],) end -function buildIntegralEntries(intCache::OneBodyIntCache{T}, +function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, intWeightPair::NTuple{2, IndexedWeights{T}}) where {T} list1, list2 = getfield.(intWeightPair, :list) intValCache = intCache.data @@ -643,7 +693,7 @@ function buildIntegralEntries(intCache::OneBodyIntCache{T}, end -function buildIntegralTensor(intCache::OneBodyIntCache{T}, +function buildIntegralTensor(intCache::OneBodyCoreIntCache{T}, intWeights::AbstractVector{IndexedWeights{T}}) where {T} nOrbs = length(intWeights) res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) @@ -682,30 +732,31 @@ function getPrimOrbCoreTypeUnion(orbs::OrbitalCollection) end -function initializeIntegralCache!(::OneBodyIntegral, op::DirectOperator, +function initializeIntegralCache!(::OneBodyIntegral{D}, op::DirectOperator, paramCache::DimSpanDataCacheBox{T}, orbInput::OrbitalInput{T, D}) where {T, D} coreType = getPrimOrbCoreTypeUnion(orbInput) orbCache = PrimOrbCoreCache(T, Val(D), coreType) - IntegralCache(op, orbCache, CompleteOneBodyIntegrals(T)) + PrimOrbCoreIntegralCache(op, OneBodyFullCoreIntegrals(T, Val(D)), orbCache) end function initializeOverlapCache!(paramCache::DimSpanDataCacheBox{T}, - orbInput::OrbitalInput{T}) where {T} - initializeIntegralCache!(OneBodyIntegral(), Identity(), paramCache, orbInput) + orbInput::OrbitalInput{T, D}) where {T, D} + initializeIntegralCache!(OneBodyIntegral{D}(), Identity(), paramCache, orbInput) end function initializeOneBodyCachePair!(op::F, paramCache::DimSpanDataCacheBox{T}, - orbInput::OrbitalInput{T}) where {F<:DirectOperator, T} - intCache = initializeIntegralCache!(OneBodyIntegral(), op, paramCache, orbInput) + orbInput::OrbitalInput{T, D}) where + {F<:DirectOperator, T, D} + intCache = initializeIntegralCache!(OneBodyIntegral{D}(), op, paramCache, orbInput) normCache = F <: Identity ? intCache : initializeOverlapCache!(paramCache, orbInput) intCache, normCache end -function computeIntTensor!(intCache::IntegralCache{T, D}, - normCache::OverlapCache{T, D}, +function computeIntTensor!(intCache::PrimOrbCoreIntegralCache{T, D}, + normCache::OverlapCoreCache{T, D}, basisSet::AbstractVector{<:FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict()) where @@ -714,7 +765,7 @@ function computeIntTensor!(intCache::IntegralCache{T, D}, buildIntegralTensor(intCache, ws) end -function computeIntTensor(::OneBodyIntegral, op::F, +function computeIntTensor(::OneBodyIntegral{D}, op::F, basisSet::AbstractVector{<:FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict()) where @@ -729,7 +780,7 @@ computeIntTensor(style::MultiBodyIntegral, op::DirectOperator, orbs::OrbitalBasi computeIntTensor(style, op, map(FrameworkOrb, orbs); paramCache, basisMarkerCache) -function decomposeOrb!(normCache::OverlapCache{T, D}, orb::FrameworkOrb{T, D}; +function decomposeOrb!(normCache::OverlapCoreCache{T, D}, orb::FrameworkOrb{T, D}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), markerCache::OrbCoreMarkerDict=OrbCoreMarkerDict()) where {T, D} normIdxList = cachePrimCoreIntegrals!(normCache, paramCache, markerCache, orb) @@ -745,8 +796,8 @@ function decomposeOrb(orb::FrameworkOrb{T}; end -function computeIntegral!(intCache::IntegralCache{T, D}, - normCache::OverlapCache{T, D}, +function computeIntegral!(intCache::PrimOrbCoreIntegralCache{T, D}, + normCache::OverlapCoreCache{T, D}, bfs::NonEmptyTuple{FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}, basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict()) where @@ -755,7 +806,7 @@ function computeIntegral!(intCache::IntegralCache{T, D}, buildIntegralEntries(intCache, ws) |> first end -function computeIntegral(::OneBodyIntegral, op::DirectOperator, +function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, (bf1,)::Tuple{FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict(), @@ -765,40 +816,42 @@ function computeIntegral(::OneBodyIntegral, op::DirectOperator, computeIntegral!(iCache, nCache, (bf1,); paramCache, basisMarkerCache) else coreData, w = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) - tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,)) + computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) + tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,); computeCache) dot(w, tensor, w) end end -function computeIntegral(::OneBodyIntegral, op::DirectOperator, +function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, (bf1, bf2)::NTuple{2, FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict(), lazyCompute::Bool=false) where {T, D} - coreData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) - coreData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) + oData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) + oData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) + computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) tensor = if lazyCompute - coreData1 = Vector(coreData1) - coreData2 = Vector(coreData2) + oData1 = Vector(oData1) + oData2 = Vector(oData2) transformation = (b::PrimitiveOrbCore{T, D})->lazyMarkObj!(basisMarkerCache, b) - coreDataM = intersectMultisets!(coreData1, coreData2; transformation) - block4 = computeOneBodyPrimCoreIntTensor(op, (coreData1, coreData2)) + coreDataM = intersectMultisets!(oData1, oData2; transformation) + block4 = computeOneBodyPrimCoreIntTensor(op, (oData1, oData2); computeCache) if isempty(coreDataM) block4 else - block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,)) - block2 = computeOneBodyPrimCoreIntTensor(op, (coreData1, coreDataM)) - block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, coreData2)) + block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,); computeCache) + block2 = computeOneBodyPrimCoreIntTensor(op, (oData1, coreDataM); computeCache) + block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, oData2); computeCache) hvcat((2, 2), block1, block3, block2, block4) end else - computeOneBodyPrimCoreIntTensor(op, (coreData1, coreData2)) + computeOneBodyPrimCoreIntTensor(op, (oData1, oData2); computeCache) end dot(w1, tensor, w2) end -function computeIntegral(::OneBodyIntegral, ::Identity, +function computeIntegral(::OneBodyIntegral{D}, ::Identity, bfPair::NTuple{2, FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict(), @@ -806,16 +859,17 @@ function computeIntegral(::OneBodyIntegral, ::Identity, bf1, bf2 = bfPair if lazyCompute if bf1 === bf2 - computeIntegral(OneBodyIntegral(), Identity(), (bf1,); paramCache, + computeIntegral(OneBodyIntegral{D}(), Identity(), (bf1,); paramCache, basisMarkerCache, lazyCompute) else normCache = initializeOverlapCache!(paramCache, bfPair) computeIntegral!(normCache, normCache, bfPair; paramCache, basisMarkerCache) end else - coreData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) - coreData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) - tensor = computeOneBodyPrimCoreIntTensor(Identity(), (coreData1, coreData2)) + oData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) + oData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) + computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) + tensor = computeOneBodyPrimCoreIntTensor(Identity(), (oData1, oData2); computeCache) dot(w1, tensor, w2) end end From 0cee1873063707eab9df19e609bb9d1be0eb8f7b Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 20 Jan 2025 21:31:10 -0500 Subject: [PATCH 24/29] Changed the field order of `PrimitiveOrb`. --- src/SpatialBasis.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index 6d644f29..a1bfa6cd 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -81,22 +81,21 @@ end #? Allow .renormalize mutable -#! Change .center to the first field struct PrimitiveOrb{T, D, B<:FieldAmplitude{T, D}, C<:NTuple{D, ElementalParam{T}}} <: ComposedOrb{T, D, B} - body::B center::C + body::B renormalize::Bool end const PrimGTO{T, D, B<:PolyGaussProd{T, D}, C<:NTuple{D, ElementalParam{T}}} = PrimitiveOrb{T, D, B, C} -function PrimitiveOrb(body::B, center::NTuple{D, ParamOrValue{T}}; +function PrimitiveOrb(center::NTuple{D, ParamOrValue{T}}, body::B; renormalize::Bool=false) where {T, D, B<:FieldAmplitude{T, D}} length(center)!=D && throw(AssertionError("The length of `center` must match `D=$D`.")) encoder = genCellEncoder(T, :cen) - PrimitiveOrb(body, encoder.(center), renormalize) + PrimitiveOrb(encoder.(center), body, renormalize) end PrimitiveOrb(ob::PrimitiveOrb) = itself(ob) @@ -434,7 +433,7 @@ function genGaussTypeOrb(center::NonEmptyTuple{ParamOrValue{T}, D}, ijk::NonEmptyTuple{Int, D}=ntuple(_->0, Val(D+1)); renormalize::Bool=false) where {T, D} gf = GaussFunc(xpn) - PrimitiveOrb(PolyRadialFunc(gf, ijk), center; renormalize) + PrimitiveOrb(center, PolyRadialFunc(gf, ijk); renormalize) end function genGaussTypeOrb(center::NonEmptyTuple{ParamOrValue{T}, D}, From ca8376b9b7d18579367d4bac43db6d0e544a510f Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 21 Jan 2025 10:45:40 -0500 Subject: [PATCH 25/29] Added type parameters. --- src/Integrals/Engines/Numerical.jl | 14 +++++++------- src/Integrals/Interface.jl | 12 ++++++------ 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/src/Integrals/Engines/Numerical.jl b/src/Integrals/Engines/Numerical.jl index 0ac0a867..44c09754 100644 --- a/src/Integrals/Engines/Numerical.jl +++ b/src/Integrals/Engines/Numerical.jl @@ -59,12 +59,12 @@ function composeIntegralKernel(::OneBodyIntegral, op::O, ::Type{T}, composeOneBodyKernel(op, T, terms) end -function numericalIntegrate(::OneBodyIntegral, op::F, +function numericalIntegrate(::OneBodyIntegral{D}, op::F, orbs::NonEmptyTuple{EvalFieldFunction{T, D}, N}, pVals::NonEmptyTuple{FilteredVecOfArr{T}, N}) where {F<:DirectOperator, T, D, N} terms = Base.Fix2.(orbs, pVals) - integralKernel = composeIntegralKernel(OneBodyIntegral(), op, T, terms) + integralKernel = composeIntegralKernel(OneBodyIntegral{D}(), op, T, terms) integrand = ConfineInterval(T, integralKernel) bound = ntuple(_->one(T), Val(D)) numericalIntegrateCore(integrand, (.-(bound), bound)) @@ -88,11 +88,11 @@ const OnyBodyPairNumInt{T, D, F<:DirectOperator, P1<:PrimitiveOrbCore{T, D}, P2<:PrimitiveOrbCore{T, D}} = OneBodyNumIntegrate{T, D, F, Tuple{P1, P2}} -function (f::OneBodySelfNumInt{T})(pVal::FilteredVecOfArr{T}) where {T} - numericalIntegrate(OneBodyIntegral(), f.op, f.basis, (pVal,)) +function (f::OneBodySelfNumInt{T, D})(pVal::FilteredVecOfArr{T}) where {T, D} + numericalIntegrate(OneBodyIntegral{D}(), f.op, f.basis, (pVal,)) end -function (f::OnyBodyPairNumInt{T})(pVal1::FilteredVecOfArr{T}, - pVal2::FilteredVecOfArr{T}) where {T} - numericalIntegrate(OneBodyIntegral(), f.op, f.basis, (pVal1, pVal2)) +function (f::OnyBodyPairNumInt{T, D})(pVal1::FilteredVecOfArr{T}, + pVal2::FilteredVecOfArr{T}) where {T, D} + numericalIntegrate(OneBodyIntegral{D}(), f.op, f.basis, (pVal1, pVal2)) end \ No newline at end of file diff --git a/src/Integrals/Interface.jl b/src/Integrals/Interface.jl index ec837bc2..02453098 100644 --- a/src/Integrals/Interface.jl +++ b/src/Integrals/Interface.jl @@ -4,14 +4,14 @@ function overlap(orb1::OrbitalBasis{T, D}, orb2::OrbitalBasis{T, D}; if orb1 === orb2 && isRenormalized(orb1) one(T) else - computeIntegral(OneBodyIntegral(), Identity(), (orb1, orb2); + computeIntegral(OneBodyIntegral{D}(), Identity(), (orb1, orb2); paramCache, lazyCompute) end end -overlaps(basisSet::OrbitalBasisSet{T}; - paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T)) where {T} = -computeIntTensor(OneBodyIntegral(), Identity(), basisSet; paramCache) +overlaps(basisSet::OrbitalBasisSet{T, D}; + paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T)) where {T, D} = +computeIntTensor(OneBodyIntegral{D}(), Identity(), basisSet; paramCache) function multipoleMoment(center::NTuple{D, Real}, degrees::NTuple{D, Int}, @@ -20,7 +20,7 @@ function multipoleMoment(center::NTuple{D, Real}, degrees::NTuple{D, Int}, lazyCompute::Bool=false) where {T, D} mmOp = MonomialMul(T.(center), degrees) orbData = orb1 === orb2 ? (orb1,) : (orb1, orb2) - computeIntegral(OneBodyIntegral(), mmOp, (orb1, orb2); paramCache, lazyCompute) + computeIntegral(OneBodyIntegral{D}(), mmOp, (orb1, orb2); paramCache, lazyCompute) end function multipoleMoments(center::NTuple{D, Real}, degrees::NTuple{D, Int}, @@ -28,5 +28,5 @@ function multipoleMoments(center::NTuple{D, Real}, degrees::NTuple{D, Int}, paramCache::DimSpanDataCacheBox{T}= DimSpanDataCacheBox(T)) where {T, D} mmOp = MonomialMul(T.(center), degrees) - computeIntTensor(OneBodyIntegral(), mmOp, basisSet; paramCache) + computeIntTensor(OneBodyIntegral{D}(), mmOp, basisSet; paramCache) end \ No newline at end of file From e812f5c57025d1229f512325d17e38d04d1727f9 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 21 Jan 2025 13:11:31 -0500 Subject: [PATCH 26/29] Changed `computeCache` from a keyword argument to a positional argument to reduce extra overhead. --- src/Integrals/Framework.jl | 68 +++++++++++++++++++------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index d0646812..34b06221 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -178,7 +178,7 @@ end function evalOneBodyPrimCoreIntegral(op::DirectOperator, - data::N12Tuple{OrbCoreData{T, D}}; + data::N12Tuple{OrbCoreData{T, D}}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} orbData = first.(data) @@ -212,26 +212,26 @@ sortTensorIndex(arg::Vararg{Int}) = sortTensorIndex(arg) function genOneBodyIntDataPairsCore(op::DirectOperator, (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - (oneBasedIdx,)::Tuple{Int}; + (oneBasedIdx,)::Tuple{Int}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} - iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],); computeCache) + iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],), computeCache) (iiVal,) end function genOneBodyIntDataPairsCore(op::DirectOperator, oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - oneBasedIdxPair::NTuple{2, Int}; + oneBasedIdxPair::NTuple{2, Int}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} orbPars1, orbPars2 = map(oDataPair, oneBasedIdxPair) do data, idx data[begin+idx-1] end - ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2); computeCache) + ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2), computeCache) jiVal = if isHermitian(first(orbPars1), op, first(orbPars2)) ijVal' else - evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1); computeCache) + evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1), computeCache) end (ijVal, jiVal) end @@ -239,21 +239,21 @@ end function genOneBodyIntDataPairs(op::DirectOperator, (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - (indexOffset,)::Tuple{Int}=(0,); + (indexOffset,)::Tuple{Int}=(0,), computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} nOrbs = length(oData) offset = indexOffset + firstindex(oData) - 1 pairs1 = map(1:nOrbs) do i - iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,); computeCache) + iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,), computeCache) (i + offset,) => iiVal end pairs2 = map(1:triMatEleNum(nOrbs-1)) do l n, m = convert1DidxTo2D(nOrbs-1, l) ijPair = sortTensorIndex((m, n+1)) - ijValPair = genOneBodyIntDataPairsCore(op, (oData, oData), ijPair; computeCache) + ijValPair = genOneBodyIntDataPairsCore(op, (oData, oData), ijPair, computeCache) (ijPair .+ offset) => ijValPair end @@ -262,7 +262,7 @@ end function genOneBodyIntDataPairs(op::DirectOperator, oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - indexOffsets::NTuple{2, Int}; + indexOffsets::NTuple{2, Int}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} mnOrbs = length.(oDataPair) @@ -272,7 +272,7 @@ function genOneBodyIntDataPairs(op::DirectOperator, idxPairOld = mnIdx .+ offsets idxPairNew = sortTensorIndex(idxPairOld) ijPair = ifelse(idxPairNew == idxPairOld, mnIdx, reverse(mnIdx)) - ijValPair = genOneBodyIntDataPairsCore(op, oDataPair, ijPair; computeCache) + ijValPair = genOneBodyIntDataPairsCore(op, oDataPair, ijPair, computeCache) idxPairNew => ijValPair end |> vec end @@ -293,7 +293,7 @@ end function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T, D}}; + (oData,)::Tuple{OrbCoreDataSeq{T, D}}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}() ) where {T, D} @@ -301,20 +301,20 @@ function computeOneBodyPrimCoreIntTensor(op::DirectOperator, res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) for i in 1:nOrbs - temp = genOneBodyIntDataPairsCore(op, (oData,), (i,); computeCache) + temp = genOneBodyIntDataPairsCore(op, (oData,), (i,), computeCache) setTensorEntries!(res, temp, (i,)) end for l in 1:triMatEleNum(nOrbs-1) n, m = convert1DidxTo2D(mBasis-1, l) - temp = genOneBodyIntDataPairsCore(op, (oData, oData), (m, n+1); computeCache) + temp = genOneBodyIntDataPairsCore(op, (oData, oData), (m, n+1), computeCache) setTensorEntries!(res, temp, (m, n+1)) end res end function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}; + oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}() ) where {T, D} @@ -323,7 +323,7 @@ function computeOneBodyPrimCoreIntTensor(op::DirectOperator, res = ShapedMemory{T}(undef, (len1, len2)) for j in 1:len2, i in 1:len1 oCorePair = (oData1[begin+i-1], oData2[begin+j-1]) - ijVal = evalOneBodyPrimCoreIntegral(op, oCorePair; computeCache) + ijVal = evalOneBodyPrimCoreIntegral(op, oCorePair, computeCache) res[begin+i-1, begin+j-1] = ijVal end res @@ -450,7 +450,7 @@ function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(orbCache.list) orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbs) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1; computeCache) + updatePrimCoreIntCache!(intCache, oldMaxIdx+1, computeCache) orbIdxList end @@ -463,7 +463,7 @@ function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(orbCache.list) orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orb) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1; computeCache) + updatePrimCoreIntCache!(intCache, oldMaxIdx+1, computeCache) orbIdxList end @@ -478,17 +478,17 @@ function cachePrimCoreIntegrals!(targetIntCache::PrimOrbCoreIntegralCache{T, D, oldMaxIdx = lastindex(tOrbCache.list) sOrbCache = sourceIntCache.basis orbIdxList = indexCacheOrbData!(tOrbCache, sOrbCache, orbMarkerCache, sourceOrbList) - updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1; computeCache) + updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1, computeCache) orbIdxList end function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, basis::Tuple{OrbCoreDataSeq{T, D}}, - offset::Tuple{Int}; + offset::Tuple{Int}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} - pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset; computeCache) + pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset, computeCache) foreach(p->setIntegralData!(ints, p), pairs1) foreach(p->setIntegralData!(ints, p), pairs2) ints @@ -496,16 +496,16 @@ end function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, basis::NTuple{2, OrbCoreDataSeq{T, D}}, - offset::NTuple{2, Int}; + offset::NTuple{2, Int}, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} - pairs2 = genOneBodyIntDataPairs(op, basis, offset; computeCache) + pairs2 = genOneBodyIntDataPairs(op, basis, offset, computeCache) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D}, startIdx::Int; +function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D}, startIdx::Int, computeCache::IntCompCacheDict{T, D}= IntCompCacheDict{T, D}()) where {T, D} op = cache.operator @@ -514,13 +514,13 @@ function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D}, startIdx firstIdx = firstindex(basis) if startIdx == firstIdx - updateIntCacheCore!(op, ints, (basis,), (0,); computeCache) + updateIntCacheCore!(op, ints, (basis,), (0,), computeCache) elseif firstIdx < startIdx <= lastindex(basis) boundary = startIdx - 1 oldBasis = @view basis[begin:boundary] newBasis = @view basis[startIdx: end] - updateIntCacheCore!(op, ints, (newBasis,), (boundary,); computeCache) - updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary); computeCache) + updateIntCacheCore!(op, ints, (newBasis,), (boundary,), computeCache) + updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary), computeCache) end cache @@ -817,7 +817,7 @@ function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, else coreData, w = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) - tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,); computeCache) + tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,), computeCache) dot(w, tensor, w) end end @@ -836,17 +836,17 @@ function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, oData2 = Vector(oData2) transformation = (b::PrimitiveOrbCore{T, D})->lazyMarkObj!(basisMarkerCache, b) coreDataM = intersectMultisets!(oData1, oData2; transformation) - block4 = computeOneBodyPrimCoreIntTensor(op, (oData1, oData2); computeCache) + block4 = computeOneBodyPrimCoreIntTensor(op, (oData1, oData2), computeCache) if isempty(coreDataM) block4 else - block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,); computeCache) - block2 = computeOneBodyPrimCoreIntTensor(op, (oData1, coreDataM); computeCache) - block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, oData2); computeCache) + block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,), computeCache) + block2 = computeOneBodyPrimCoreIntTensor(op, (oData1, coreDataM), computeCache) + block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, oData2), computeCache) hvcat((2, 2), block1, block3, block2, block4) end else - computeOneBodyPrimCoreIntTensor(op, (oData1, oData2); computeCache) + computeOneBodyPrimCoreIntTensor(op, (oData1, oData2), computeCache) end dot(w1, tensor, w2) end @@ -869,7 +869,7 @@ function computeIntegral(::OneBodyIntegral{D}, ::Identity, oData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) oData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) - tensor = computeOneBodyPrimCoreIntTensor(Identity(), (oData1, oData2); computeCache) + tensor = computeOneBodyPrimCoreIntTensor(Identity(), (oData1, oData2), computeCache) dot(w1, tensor, w2) end end From 9608b8115bf37494a138fbd9b90ce4b2258ac9fd Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 27 Jan 2025 23:34:59 -0500 Subject: [PATCH 27/29] Optimized the interface for cached computing of integrals. --- src/Integrals/Engines/GaussianOrbitals.jl | 43 ++- src/Integrals/Engines/Numerical.jl | 2 +- src/Integrals/Framework.jl | 445 +++++++++++++--------- src/SpatialBasis.jl | 23 +- src/Types.jl | 7 +- 5 files changed, 320 insertions(+), 200 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 6800d4b4..6091f7ea 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -97,17 +97,21 @@ struct ZeroAngMomCache{T} <: CustomCache{T} end const CoreIntCacheBox{T} = Union{CustomCache{T}, LRU{<:Any, T}} +abstract type OrbIntegralComputeCache{T, D, S<:MultiBodyIntegral{D} + } <: IntegralProcessCache{T, D} end -struct OneBodyIntProcessCache{T, D, F<:NTuple{D, CoreIntCacheBox{T}} - } <: IntegralProcessCache{T, D} +const Orb1BIntegralComputeCache{T, D} = OrbIntegralComputeCache{T, D, OneBodyIntegral{D}} + +struct AxialOneBodyIntCompCache{T, D, F<:NTuple{D, CoreIntCacheBox{T}} + } <: Orb1BIntegralComputeCache{T, D} axis::F - OneBodyIntProcessCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = + AxialOneBodyIntCompCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = new{T, D+1, typeof(axis)}(axis) end -OneBodyIntProcessCache(::Type{T}, ::Val{D}) where {T, D} = -OneBodyIntProcessCache(ntuple( _->NullCache{T}(), Val(D) )) +AxialOneBodyIntCompCache(::Type{T}, ::Val{D}) where {T, D} = +AxialOneBodyIntCompCache(ntuple( _->NullCache{T}(), Val(D) )) function overlapPGTOcore(input::TupleOf5T2Int{T}) where {T} @@ -151,7 +155,7 @@ const GTOrbOverlapAxialCache{T} = Union{NullCache{T}, ZeroAngMomCache{T}, LRU{TupleOf5T2Int{T}, T}} const GTOrbOverlapCache{T, D} = - OneBodyIntProcessCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} + AxialOneBodyIntCompCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} function overlapPGTO(data::GaussProductInfo{T, D}, cache::GTOrbOverlapCache{T, D}) where {T, D} @@ -210,7 +214,7 @@ end function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; cache::GTOrbOverlapCache{T, D}= - OneBodyIntProcessCache(T, Val(D))) where {T, D} + AxialOneBodyIntCompCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) overlapPGTO(data, cache) end @@ -223,23 +227,34 @@ end const DefaultPGTOrbOverlapCacheSizeLimit = 128 -function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} - ntuple(_->LRU{TupleOf5T2Int{T}, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit), Val(D)) +# function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} +# ntuple(_->LRU{TupleOf5T2Int{T}, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit), Val(D)) +# end + +# function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} +# ntuple(_->ZeroAngMomCache{T}(), Val(D)) +# end + +@generated function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} + maxsize = DefaultPGTOrbOverlapCacheSizeLimit + caches = ntuple(_->LRU{TupleOf5T2Int{T}, T}(; maxsize), Val(D)) + return :( $caches ) end -function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} - ntuple(_->ZeroAngMomCache{T}(), Val(D)) +@generated function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} + caches = ntuple(_->ZeroAngMomCache{T}(), Val(D)) + return :( $caches ) end genGTOrbIntCompCache(::MonomialMul{T, D, L}, ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where {T, D, L, L1, L2} = -genPrimGTOrbOverlapCache(T, Val(D), Val(L+L1+L2)) |> OneBodyIntProcessCache +genPrimGTOrbOverlapCache(T, Val(D), Val(L+L1+L2)) |> AxialOneBodyIntCompCache genGTOrbIntCompCache(::Identity, ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where {T, D, L1, L2} = -genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> OneBodyIntProcessCache +genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> AxialOneBodyIntCompCache function buildNormalizerCore(o::PrimGTOcore{T, D}) where {T, D} @@ -317,7 +332,7 @@ end function (f::MultiMomentGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; cache::GTOrbOverlapCache{T, D}= - OneBodyIntProcessCache(T, Val(D))) where {T, D} + AxialOneBodyIntCompCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) computeMultiMomentGTO(f.op, data, cache) end diff --git a/src/Integrals/Engines/Numerical.jl b/src/Integrals/Engines/Numerical.jl index 44c09754..b3df99ec 100644 --- a/src/Integrals/Engines/Numerical.jl +++ b/src/Integrals/Engines/Numerical.jl @@ -72,7 +72,7 @@ end function buildNormalizerCore(o::PrimitiveOrbCore{T, D}) where {T, D} - buildOneBodyCoreIntegrator(Identity(), (o,)) + buildCoreIntegrator(OneBodyIntegral{D}(), Identity(), (o,)) end diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 34b06221..cccd2b33 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -71,23 +71,22 @@ end # ab::Dict{NTuple{2, Int}, T} # end - -struct PrimOrbCoreIntegralCache{T, D, F<:DirectOperator, - I<:IntegralData{T, <:MultiBodyIntegral{D}}, - B<:PrimitiveOrbCore{T, D}} <: SpatialIntegralCache{T, D} +#!!! Just changed the signature +struct PrimOrbCoreIntegralCache{T, D, S<:MultiBodyIntegral{D}, F<:DirectOperator, + I<:IntegralData{T, S}, B<:PrimitiveOrbCore{T, D} + } <: SpatialProcessCache{T, D} operator::F data::I basis::PrimOrbCoreCache{T, D, B} end -const OneBodyCoreIntCache{T, D, F<:DirectOperator, I<:IntegralData{T, OneBodyIntegral{D}}, - B<:PrimitiveOrbCore{T, D}} = - PrimOrbCoreIntegralCache{T, D, F, I, B} +const POrb1BCoreICache{T, D, F<:DirectOperator, I<:IntegralData{T, OneBodyIntegral{D}}, + B<:PrimitiveOrbCore{T, D}} = + PrimOrbCoreIntegralCache{T, D, OneBodyIntegral{D}, F, I, B} const OverlapCoreCache{T, D, I<:IntegralData{T, OneBodyIntegral{D}}, B<:PrimitiveOrbCore{T, D}} = - OneBodyCoreIntCache{T, D, Identity, I, B} - + POrb1BCoreICache{T, D, Identity, I, B} function setIntegralData!(ints::OneBodyFullCoreIntegrals{T}, pair::Pair{Tuple{Int}, Tuple{T}}) where {T} @@ -120,74 +119,226 @@ isHermitian(::PrimitiveOrbCore{T, D}, ::Identity, true -buildOneBodyCoreIntegrator(op::DirectOperator, - orbs::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = -OneBodyNumIntegrate(op, orbs) +struct Bar <: Any end + +const OrbCoreIntLayoutAllOrbs{T, D, N} = NonEmptyTuple{PrimitiveOrbCore{T, D}, N} +const OrbCoreIntLayoutOrbBar1{T, D} = + Tuple{PrimitiveOrbCore{T, D}, Bar, PrimitiveOrbCore{T, D}} +const OrbCoreIntLayoutOrbBar2{T, D} = + Tuple{PrimitiveOrbCore{T, D}, PrimitiveOrbCore{T, D}, Bar} +const OrbCoreIntLayoutOrbBar3{T, D} = + Tuple{PrimitiveOrbCore{T, D}, Bar, PrimitiveOrbCore{T, D}, PrimitiveOrbCore{T, D}} +const OrbCoreIntLayoutOrbBar4{T, D} = + Tuple{PrimitiveOrbCore{T, D}, PrimitiveOrbCore{T, D}, Bar, PrimitiveOrbCore{T, D}} + +const OneBodyOrbCoreIntLayout{T, D} = Union{ + OrbCoreIntLayoutAllOrbs{T, D, 0}, OrbCoreIntLayoutAllOrbs{T, D, 1} +} + +const TwoBodyOrbCoreIntLayout{T, D} = Union{ + OrbCoreIntLayoutAllOrbs{T, D, 0}, OrbCoreIntLayoutAllOrbs{T, D, 3}, + OrbCoreIntLayoutOrbBar1{T, D}, OrbCoreIntLayoutOrbBar2{T, D}, + OrbCoreIntLayoutOrbBar3{T, D}, OrbCoreIntLayoutOrbBar4{T, D}, +} + +# const OrbCoreIntLayout{T, D, U<:Union{Bar, PrimitiveOrbCore{T, D}}, +# PrimitiveOrbCore{T, D}<:V<:U, N} = +# Tuple{PrimitiveOrbCore{T, D}, Vararg{V, N}} + +const OrbCoreIntLayout{T, D} = + Union{OneBodyOrbCoreIntLayout{T, D}, TwoBodyOrbCoreIntLayout{T, D}} + +struct OrbCoreIntConfig{T, D, S<:MultiBodyIntegral{D}, + P<:OrbCoreIntLayout{T, D}} <: StructuredType + + function OrbCoreIntConfig(::OneBodyIntegral{D}, ::P) where + {D, T, P<:OneBodyOrbCoreIntLayout{T, D}} + new{T, D, OneBodyIntegral{D}, P}() + end + + function OrbCoreIntConfig(::TwoBodyIntegral{D}, ::P) where + {D, T, P<:TwoBodyOrbCoreIntLayout{T, D}} + new{T, D, TwoBodyIntegral{D}, P}() + end +end + +const OneBodyOrbCoreIntConfig{T, D, P<:OneBodyOrbCoreIntLayout{T, D}} = + OrbCoreIntConfig{T, D, OneBodyIntegral{D}, P} -buildOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -genGTOrbOverlapFunc(orbs) +const TwoBodyOrbCoreIntConfig{T, D, P<:TwoBodyOrbCoreIntLayout{T, D}} = + OrbCoreIntConfig{T, D, TwoBodyIntegral{D}, P} -buildOneBodyCoreIntegrator(op::MonomialMul{T, D}, - orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -genGTOrbMultiMomentFunc(op, orbs) +abstract type OrbIntegralCompCacheBox{T, D, S<:MultiBodyIntegral{D}} <: QueryBox{T} end -const OneBodyCoreIntData{T, D} = Tuple{DirectOperator, Vararg{PrimitiveOrbCore{T, D}, 2}} +struct OrbIntCompCacheBox{T, D, S<:MultiBodyIntegral{D}} <: OrbIntegralCompCacheBox{T, D, S} + dict::Dict{OrbCoreIntConfig{T, D, S}, OrbIntegralComputeCache{T, D, S}} -const OneBodyCICacheDict{T, D} = - Dict{Type{<:OneBodyCoreIntData{T, D}}, OneBodyIntProcessCache{T, D}} -#! Consider a new type union when two-body computation cache is needed -const IntCompCacheDict{T, D} = Union{OneBodyCICacheDict{T, D}, TypedEmptyDict{Val{D}, T}} + function OrbIntCompCacheBox(::S, ::Type{T}) where {D, S<:MultiBodyIntegral{D}, T} + dict = Dict{OrbCoreIntConfig{T, D, S}, OrbIntegralComputeCache{T, D, S}}() + new{T, D, S}(dict) + end +end -IntCompCacheDict{T, D}() where {T, D} = TypedEmptyDict{Val{D}, T}() +struct OrbIntNullCacheBox{T, D, S<:MultiBodyIntegral{D}} <: OrbIntegralCompCacheBox{T, D, S} -IntCompCacheDict(::Type{T}, ::Type{OneBodyIntegral{D}}) where {T, D} = -OneBodyCICacheDict{T, D}() + OrbIntCompCacheBox(::S, ::Type{T}) where {D, S<:MultiBodyIntegral{D}, T} = + new{T, D, S}() +end -IntCompCacheDict(::Type{T}, ::Type{TwoBodyIntegral{D}}) where {T, D} = -TypedEmptyDict{Val{D}, T}() +abstract type CachedComputeSpatialIntegral{T, D} <: StatefulFunction{T} end + +struct CachedComputeOrbCoreIntegral{T, D, F<:DirectOperator, + C<:OrbIntegralCompCacheBox{T, D} + } <: CachedComputeSpatialIntegral{T, D} + operator::F + cache::C -IntCompCacheDict(::Type{<:IntegralData{T, S}}) where {T, S<:MultiBodyIntegral} = -IntCompCacheDict(T, S) + function CachedComputeOrbCoreIntegral(::Val{true}, ::S, operator::F, ::Type{T}) where + {F<:DirectOperator, T, D, S<:MultiBodyIntegral{D}} + new{T, D, F, OrbIntCompCacheBox{T, D, S}}(operator, OrbIntCompCacheBox(S(), T)) + end -function prepareOneBodyIntCompCache!(cacheDict::OneBodyCICacheDict{T, D}, - op::DirectOperator, - orbs::NTuple{2, PrimGTOcore{T, D}}) where {T, D} - get!(cacheDict, typeof( (op, orbs...) )) do - genGTOrbIntCompCache(op, orbs) + function CachedComputeOrbCoreIntegral(::Val{false}, ::S, operator::F, ::Type{T}) where + {F<:DirectOperator, T, D, S<:MultiBodyIntegral{D}} + new{T, D, F, OrbIntNullCacheBox{T, D, S}}(operator, OrbIntNullCacheBox(S(), T)) end end -function prepareOneBodyIntCompCache!(::IntCompCacheDict{T, D}, ::DirectOperator, - ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} - NullCache{T}() +const OrbCoreIntegralComputeConfig{T, D, S<:MultiBodyIntegral{D}, F<:DirectOperator, + C<:OrbIntegralCompCacheBox{T, D, S}} = + CachedComputeOrbCoreIntegral{T, D, F, C} + +const Orb1BCCIntComputeConfigUnion{T, D, F<:DirectOperator} = + OrbCoreIntegralComputeConfig{T, D, OneBodyIntegral{D}, F} + +const Orb2BCCIntComputeConfigUnion{T, D, F<:DirectOperator} = + OrbCoreIntegralComputeConfig{T, D, TwoBodyIntegral{D}, F} + +const DirectComputeOrbCoreIntegral{T, D, S<:MultiBodyIntegral{D}, F<:DirectOperator, + C<:OrbIntNullCacheBox{T, D, S}} = + CachedComputeOrbCoreIntegral{T, D, F, C} + +const ReusedComputeOrbCoreIntegral{T, D, S<:MultiBodyIntegral{D}, F<:DirectOperator, + C<:OrbIntCompCacheBox{T, D, S}} = + CachedComputeOrbCoreIntegral{T, D, F, C} + + +function buildCoreIntegrator(::OneBodyIntegral{D}, op::DirectOperator, + orbs::N12Tuple{PrimitiveOrbCore{T, D}} + ) where {T, D} + OneBodyNumIntegrate(op, orbs) +end + +function buildCoreIntegrator(::OneBodyIntegral{D}, ::Identity, + orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + genGTOrbOverlapFunc(orbs) end -# adjustOneBodyCoreIntConfig +function buildCoreIntegrator(::OneBodyIntegral{D}, op::MonomialMul{T, D}, + orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + genGTOrbMultiMomentFunc(op, orbs) +end -function lazyEvalCoreIntegral!(cache::IntegralProcessCache{T}, - integrator::OrbitalIntegrator{T}, - paramData::N12Tuple{FilteredVecOfArr{T}}) where {T} - integrator(paramData...; cache)::T +function (f::DirectComputeOrbCoreIntegral{T, D, S} + )(data::N12Tuple{OrbCoreData{T, D}}) where {T, D, S<:MultiBodyIntegral{D}} + orbData = first.(data) + parData = last.(data) + fCore = buildCoreIntegrator(S(), f.operator, orbData) + fCore(parData...) end -function lazyEvalCoreIntegral!(::NullCache{T}, integrator::OrbitalIntegrator{T}, - paramData::N12Tuple{FilteredVecOfArr{T}}) where {T} - integrator(paramData...)::T + +function prepareOneBodyIntCompCache!(f::ReusedComputeOrbCoreIntegral{T, D}, + orbs::TetraTupleUnion{PrimGTOcore{T, D}}) where {T, D} + # get!(f.cache, typeof( (op, orbs...) )) do + genGTOrbIntCompCache(f.operator, orbs) + # end end +prepareOneBodyIntCompCache!(::ReusedComputeOrbCoreIntegral{T, D}, + ::TetraTupleUnion{PrimitiveOrbCore{T, D}}) where {T, D} = +nothing -function evalOneBodyPrimCoreIntegral(op::DirectOperator, - data::N12Tuple{OrbCoreData{T, D}}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} +prepareOneBodyIntCompCache!(::ReusedComputeOrbCoreIntegral{T, D, OneBodyIntegral{D}}, + ::Tuple{PrimGTOcore{T, D}}) where {T, D} = +nothing + +function (f::OrbCoreIntegralComputeConfig{T, D, S} + )(data::TetraTupleUnion{OrbCoreData{T, D}}) where {T, D, S<:MultiBodyIntegral{D}} orbData = first.(data) parData = last.(data) - fCore = buildOneBodyCoreIntegrator(op, orbData) - intCompCache = prepareOneBodyIntCompCache!(computeCache, op, orbData) - lazyEvalCoreIntegral!(intCompCache, fCore, parData) + fCore = buildCoreIntegrator(S(), f.operator, orbData) + intCacheSector = prepareOneBodyIntCompCache!(f, orbData) + if intCacheSector === nothing + fCore(parData...) + else + fCore(parData..., cache=intCacheSector) + end end +# One-Body (i|O|j) hermiticity across O +getHermiticity(::PrimitiveOrbCore{T, D}, ::DirectOperator, + ::PrimitiveOrbCore{T, D}) where {T, D} = +false + +getHermiticity(::PrimitiveOrbCore{T, D}, ::Identity, + ::PrimitiveOrbCore{T, D}) where {T, D} = +true + +# Two-Body (ii|O|jj) (ii|O|jk) (ij|O|kk) (ij|O|kl) hermiticity across O +getHermiticity(::N12Tuple{PrimitiveOrbCore{T, D}}, ::Identity, + ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = +true + +function genCoreIntTuple(integrator::Orb1BCCIntComputeConfigUnion{T, D}, + oDataTuple::Tuple{OrbCoreData{T, D}}) where {T, D} + (integrator(oDataTuple),) +end + +function genCoreIntTuple(integrator::Orb1BCCIntComputeConfigUnion{T, D}, + oDataTuple::NTuple{2, OrbCoreData{T, D}}) where {T, D} + d1, d2 = oDataTuple + ijVal = integrator(oDataTuple) + jiVal = if getHermiticity(first(d1), integrator.operator, first(d2)) + ijVal' + else + integrator(oDataTuple|>reverse) + end + (ijVal, jiVal) +end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::Tuple{OrbCoreDataSeq{T, D}}, +# oneBasedIds::Tuple{Int}) where {T, D} +# end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::Tuple{OrbCoreDataSeq{T, D}, Bar, OrbCoreDataSeq{T, D}}, +# oneBasedIds::NTuple{2, Int}) where {T, D} +# end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::Tuple{OrbCoreDataSeq{T, D}, OrbCoreDataSeq{T, D}, Bar}, +# oneBasedIds::NTuple{2, Int}) where {T, D} +# end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::Tuple{OrbCoreDataSeq{T, D}, Bar, OrbCoreDataSeq{T, D}, +# OrbCoreDataSeq{T, D}}, +# oneBasedIds::NTuple{3, Int}) where {T, D} +# end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::Tuple{OrbCoreDataSeq{T, D}, OrbCoreDataSeq{T, D}, Bar, +# OrbCoreDataSeq{T, D}}, +# oneBasedIds::NTuple{3, Int}) where {T, D} +# end + +# function genCoreIntTuple(integrator::Orb2BCCIntComputeConfigUnion{T, D}, +# oDataTuple::NTuple{4, OrbCoreDataSeq{T, D}}, +# oneBasedIds::NTuple{4, Int}) where {T, D} +# end function sortTensorIndex((i, j)::NTuple{2, Int}) if i > j @@ -210,69 +361,42 @@ end sortTensorIndex(arg::Vararg{Int}) = sortTensorIndex(arg) -function genOneBodyIntDataPairsCore(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - (oneBasedIdx,)::Tuple{Int}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - iiVal = evalOneBodyPrimCoreIntegral(op, (oData[begin+oneBasedIdx-1],), computeCache) - (iiVal,) -end - -function genOneBodyIntDataPairsCore(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - oneBasedIdxPair::NTuple{2, Int}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - orbPars1, orbPars2 = map(oDataPair, oneBasedIdxPair) do data, idx - data[begin+idx-1] - end - ijVal = evalOneBodyPrimCoreIntegral(op, (orbPars1, orbPars2), computeCache) - jiVal = if isHermitian(first(orbPars1), op, first(orbPars2)) - ijVal' - else - evalOneBodyPrimCoreIntegral(op, (orbPars2, orbPars1), computeCache) - end - (ijVal, jiVal) -end - - -function genOneBodyIntDataPairs(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - (indexOffset,)::Tuple{Int}=(0,), - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - nOrbs = length(oData) - offset = indexOffset + firstindex(oData) - 1 +function genOneBodyIntDataPairs(integrator::CachedComputeOrbCoreIntegral{T, D}, + (oDataSeq,)::Tuple{OrbCoreDataSeq{T, D}}, + (indexOffset,)::Tuple{Int}=(0,)) where {T, D} + nOrbs = length(oDataSeq) + firstIdx = firstindex(oDataSeq) + offset = indexOffset + firstIdx - 1 pairs1 = map(1:nOrbs) do i - iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,), computeCache) + iiVal = genCoreIntTuple(integrator, (oDataSeq[begin+i-1],)) (i + offset,) => iiVal end pairs2 = map(1:triMatEleNum(nOrbs-1)) do l n, m = convert1DidxTo2D(nOrbs-1, l) - ijPair = sortTensorIndex((m, n+1)) - ijValPair = genOneBodyIntDataPairsCore(op, (oData, oData), ijPair, computeCache) + i, j = ijPair = sortTensorIndex((m, n+1)) + oDataPair = (oDataSeq[begin+i-1], oDataSeq[begin+j-1]) + ijValPair = genCoreIntTuple(integrator, oDataPair) (ijPair .+ offset) => ijValPair end pairs1, pairs2 end -function genOneBodyIntDataPairs(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - indexOffsets::NTuple{2, Int}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - mnOrbs = length.(oDataPair) - offsets = indexOffsets .+ firstindex.(oDataPair) .- 1 +function genOneBodyIntDataPairs(integrator::CachedComputeOrbCoreIntegral{T, D}, + oDataSeqPair::NTuple{2, OrbCoreDataSeq{T, D}}, + indexOffsets::NTuple{2, Int}) where {T, D} + mnOrbs = length.(oDataSeqPair) + firstIdx = firstindex.(oDataSeqPair) + offsets = indexOffsets .+ firstindex.(oDataSeqPair) .- 1 map(Iterators.product( Base.OneTo.(mnOrbs)... )) do mnIdx idxPairOld = mnIdx .+ offsets idxPairNew = sortTensorIndex(idxPairOld) ijPair = ifelse(idxPairNew == idxPairOld, mnIdx, reverse(mnIdx)) - ijValPair = genOneBodyIntDataPairsCore(op, oDataPair, ijPair, computeCache) + oDataPair = getindex.(oDataSeqPair, firstIdx .+ ijPair) + ijValPair = genCoreIntTuple(integrator, oDataPair) idxPairNew => ijValPair end |> vec end @@ -292,38 +416,34 @@ function setTensorEntries!(tensor::AbstractMatrix{T}, valPair::NTuple{2, T}, end -function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T, D}}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}() - ) where {T, D} - nOrbs = length(oData) +function genOneBodyPrimCoreIntTensor(integrator::CachedComputeOrbCoreIntegral{T, D}, + (oDataSeq,)::Tuple{OrbCoreDataSeq{T, D}}) where {T, D} + nOrbs = length(oDataSeq) res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) for i in 1:nOrbs - temp = genOneBodyIntDataPairsCore(op, (oData,), (i,), computeCache) + temp = genCoreIntTuple(integrator, (oDataSeq[begin+i-1],)) setTensorEntries!(res, temp, (i,)) end for l in 1:triMatEleNum(nOrbs-1) n, m = convert1DidxTo2D(mBasis-1, l) - temp = genOneBodyIntDataPairsCore(op, (oData, oData), (m, n+1), computeCache) + oDataPair = (oDataSeq[begin+m-1], oDataSeq[begin+n]) + temp = genCoreIntTuple(integrator, oDataPair) setTensorEntries!(res, temp, (m, n+1)) end res end -function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T, D}}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}() - ) where {T, D} - oData1, oData2 = oDataPair - len1, len2 = length.(oDataPair) +function genOneBodyPrimCoreIntTensor(integrator::CachedComputeOrbCoreIntegral{T, D}, + oDataSeqPair::NTuple{2, OrbCoreDataSeq{T, D}}) where + {T, D} + oData1, oData2 = oDataSeqPair + len1, len2 = length.(oDataSeqPair) res = ShapedMemory{T}(undef, (len1, len2)) for j in 1:len2, i in 1:len1 - oCorePair = (oData1[begin+i-1], oData2[begin+j-1]) - ijVal = evalOneBodyPrimCoreIntegral(op, oCorePair, computeCache) + oDataPair = (oData1[begin+i-1], oData2[begin+j-1]) + ijVal = integrator(oDataPair) res[begin+i-1, begin+j-1] = ijVal end res @@ -440,87 +560,66 @@ function indexCacheOrbData!(targetCache::PrimOrbCoreCache{T, D}, end -function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, +function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D}, paramCache::DimSpanDataCacheBox{T}, orbMarkerCache::OrbCoreMarkerDict, - orbs::OrbitalCollection{T, D}) where {T, D, - F<:DirectOperator, - I<:IntegralData{ T, <:MultiBodyIntegral{D} }} + orbData::Union{OrbitalCollection{T, D}, FrameworkOrb{T, D}} + ) where {T, D} orbCache = intCache.basis - computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(orbCache.list) - orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbs) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1, computeCache) + orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbData) + updatePrimCoreIntCache!(intCache, oldMaxIdx+1) orbIdxList end -function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D, F, I}, - paramCache::DimSpanDataCacheBox{T}, +function cachePrimCoreIntegrals!(targetIntCache::C, sourceIntCache::C, orbMarkerCache::OrbCoreMarkerDict, - orb::FrameworkOrb{T, D}) where {T, D, F<:DirectOperator, - I<:IntegralData{ T, <:MultiBodyIntegral{D} }} - orbCache = intCache.basis - computeCache = IntCompCacheDict(I) - oldMaxIdx = lastindex(orbCache.list) - orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orb) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1, computeCache) - orbIdxList -end - -function cachePrimCoreIntegrals!(targetIntCache::PrimOrbCoreIntegralCache{T, D, F, I}, - sourceIntCache::PrimOrbCoreIntegralCache{T, D, F, I}, - orbMarkerCache::OrbCoreMarkerDict, - sourceOrbList::BasisIndexList) where {T, D, - F<:DirectOperator, - I<:IntegralData{ T, <:MultiBodyIntegral{D} }} + sourceOrbList::BasisIndexList) where + {C<:PrimOrbCoreIntegralCache} tOrbCache = targetIntCache.basis - computeCache = IntCompCacheDict(I) oldMaxIdx = lastindex(tOrbCache.list) sOrbCache = sourceIntCache.basis orbIdxList = indexCacheOrbData!(tOrbCache, sOrbCache, orbMarkerCache, sourceOrbList) - updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1, computeCache) + updatePrimCoreIntCache!(targetIntCache, oldMaxIdx+1) orbIdxList end -function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, +function updateIntCacheCore!(integrator::CachedComputeOrbCoreIntegral{T, D}, + ints::OneBodyFullCoreIntegrals{T, D}, basis::Tuple{OrbCoreDataSeq{T, D}}, - offset::Tuple{Int}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset, computeCache) + offset::Tuple{Int}) where {T, D} + pairs1, pairs2 = genOneBodyIntDataPairs(integrator, basis, offset) foreach(p->setIntegralData!(ints, p), pairs1) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updateIntCacheCore!(op::DirectOperator, ints::OneBodyFullCoreIntegrals{T, D}, +function updateIntCacheCore!(integrator::CachedComputeOrbCoreIntegral{T, D}, + ints::OneBodyFullCoreIntegrals{T, D}, basis::NTuple{2, OrbCoreDataSeq{T, D}}, - offset::NTuple{2, Int}, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - pairs2 = genOneBodyIntDataPairs(op, basis, offset, computeCache) + offset::NTuple{2, Int}) where {T, D} + pairs2 = genOneBodyIntDataPairs(integrator, basis, offset) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D}, startIdx::Int, - computeCache::IntCompCacheDict{T, D}= - IntCompCacheDict{T, D}()) where {T, D} - op = cache.operator +function updatePrimCoreIntCache!(cache::PrimOrbCoreIntegralCache{T, D, S}, + startIdx::Int) where {T, D, S<:MultiBodyIntegral{D}} basis = cache.basis.list ints = cache.data firstIdx = firstindex(basis) + integrator = CachedComputeOrbCoreIntegral(Val(true), S(), cache.operator, T) if startIdx == firstIdx - updateIntCacheCore!(op, ints, (basis,), (0,), computeCache) + updateIntCacheCore!(integrator, ints, (basis,), (0,)) elseif firstIdx < startIdx <= lastindex(basis) boundary = startIdx - 1 oldBasis = @view basis[begin:boundary] newBasis = @view basis[startIdx: end] - updateIntCacheCore!(op, ints, (newBasis,), (boundary,), computeCache) - updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary), computeCache) + updateIntCacheCore!(integrator, ints, (newBasis,), (boundary,)) + updateIntCacheCore!(integrator, ints, (oldBasis, newBasis,), (0, boundary)) end cache @@ -657,7 +756,7 @@ function prepareIntegralConfig!(intCache::PrimOrbCoreIntegralCache{T, D}, end -function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, +function buildIntegralEntries(intCache::POrb1BCoreICache{T}, (intWeights,)::Tuple{IndexedWeights{T}}) where {T} idxList = intWeights.list len = length(idxList) @@ -675,7 +774,7 @@ function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, (res,) # ([1|O|1],) end -function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, +function buildIntegralEntries(intCache::POrb1BCoreICache{T}, intWeightPair::NTuple{2, IndexedWeights{T}}) where {T} list1, list2 = getfield.(intWeightPair, :list) intValCache = intCache.data @@ -693,7 +792,7 @@ function buildIntegralEntries(intCache::OneBodyCoreIntCache{T}, end -function buildIntegralTensor(intCache::OneBodyCoreIntCache{T}, +function buildIntegralTensor(intCache::POrb1BCoreICache{T}, intWeights::AbstractVector{IndexedWeights{T}}) where {T} nOrbs = length(intWeights) res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) @@ -816,8 +915,8 @@ function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, computeIntegral!(iCache, nCache, (bf1,); paramCache, basisMarkerCache) else coreData, w = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) - computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) - tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,), computeCache) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) + tensor = genOneBodyPrimCoreIntTensor(integrator, (coreData,)) dot(w, tensor, w) end end @@ -829,24 +928,24 @@ function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, lazyCompute::Bool=false) where {T, D} oData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) oData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) - computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) tensor = if lazyCompute oData1 = Vector(oData1) oData2 = Vector(oData2) transformation = (b::PrimitiveOrbCore{T, D})->lazyMarkObj!(basisMarkerCache, b) coreDataM = intersectMultisets!(oData1, oData2; transformation) - block4 = computeOneBodyPrimCoreIntTensor(op, (oData1, oData2), computeCache) + block4 = genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) if isempty(coreDataM) block4 else - block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,), computeCache) - block2 = computeOneBodyPrimCoreIntTensor(op, (oData1, coreDataM), computeCache) - block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, oData2), computeCache) + block1 = genOneBodyPrimCoreIntTensor(integrator, (coreDataM,)) + block2 = genOneBodyPrimCoreIntTensor(integrator, (oData1, coreDataM)) + block3 = genOneBodyPrimCoreIntTensor(integrator, (coreDataM, oData2)) hvcat((2, 2), block1, block3, block2, block4) end else - computeOneBodyPrimCoreIntTensor(op, (oData1, oData2), computeCache) + genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) end dot(w1, tensor, w2) end @@ -856,10 +955,12 @@ function computeIntegral(::OneBodyIntegral{D}, ::Identity, paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), basisMarkerCache::OrbCoreMarkerDict=OrbCoreMarkerDict(), lazyCompute::Bool=false) where {T, D} + op = Identity() bf1, bf2 = bfPair + if lazyCompute if bf1 === bf2 - computeIntegral(OneBodyIntegral{D}(), Identity(), (bf1,); paramCache, + computeIntegral(OneBodyIntegral{D}(), op, (bf1,); paramCache, basisMarkerCache, lazyCompute) else normCache = initializeOverlapCache!(paramCache, bfPair) @@ -868,8 +969,8 @@ function computeIntegral(::OneBodyIntegral{D}, ::Identity, else oData1, w1 = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) oData2, w2 = decomposeOrb(bf2; paramCache, markerCache=basisMarkerCache) - computeCache = IntCompCacheDict(T, OneBodyIntegral{D}) - tensor = computeOneBodyPrimCoreIntTensor(Identity(), (oData1, oData2), computeCache) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) + tensor = genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) dot(w1, tensor, w2) end end diff --git a/src/SpatialBasis.jl b/src/SpatialBasis.jl index a1bfa6cd..5b4ee978 100644 --- a/src/SpatialBasis.jl +++ b/src/SpatialBasis.jl @@ -467,13 +467,13 @@ end function buildNormalizer(f::PrimitiveOrbCore{T}) where {T} - nfInner = (NormalizePrimOrb∘buildNormalizerCore)(f) - EvalOrbNormalizer(nfInner, T) + nfCore = (NormalizePrimOrb∘buildNormalizerCore)(f) + EvalOrbNormalizer(nfCore, T) end function buildNormalizer(f::CompositeOrbCore{T}) where {T} - nfInner = NormalizeCompOrb(f) - EvalOrbNormalizer(nfInner, T) + nfCore = NormalizeCompOrb(f) + EvalOrbNormalizer(nfCore, T) end function buildNormalizer(::Type{T}, ::Val{D}) where {T, D} @@ -495,6 +495,7 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} nOrbs = length(weightedOrbs) oCorePtr = ChainPointer((:left, :f, :apply, :left)) oScopePtr = ChainPointer((:left, :f, :scope)) + nIntStyle = OneBodyIntegral{D}() hasUnit = false diagFuncCoreType = Union{} @@ -507,9 +508,9 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} else oCore = getField(weightedOrb, oCorePtr) oSelect = getField(weightedOrb, oScopePtr) |> first |> Retrieve - nfInner = buildOneBodyCoreIntegrator(Identity(), (oCore,)) - diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfInner)) - ReturnTyped(EncodeApply((oSelect,), nfInner), T) + nfCore = buildCoreIntegrator(nIntStyle, Identity(), (oCore,)) + diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfCore)) + ReturnTyped(EncodeApply((oSelect,), nfCore), T) end end @@ -527,11 +528,11 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} weightedOrb2 = weightedOrbs[begin+n-1] oCore1 = getField(weightedOrb1, oCorePtr) oCore2 = getField(weightedOrb2, oCorePtr) - nfInner = buildOneBodyCoreIntegrator(Identity(), (oCore1, oCore2)) - utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfInner)) + nfCore = buildCoreIntegrator(nIntStyle, Identity(), (oCore1, oCore2)) + utriFuncCoreType = typejoin(utriFuncCoreType, typeof(nfCore)) oSelect1 = getField(weightedOrb1, oScopePtr) |> first |> Retrieve oSelect2 = getField(weightedOrb2, oScopePtr) |> first |> Retrieve - ReturnTyped(EncodeApply((oSelect1, oSelect2), nfInner), T) + ReturnTyped(EncodeApply((oSelect1, oSelect2), nfCore), T) end utriFuncType = eltype(utriFuncs) @@ -541,7 +542,7 @@ function NormalizeCompOrb(f::CompositeOrbCore{T, D}) where {T, D} utriFuncs = Memory{utriFuncType}(utriFuncs) else oCore = getField(weightedOrbs[begin], oCorePtr) - utriFuncCoreType = (typeof∘buildOneBodyCoreIntegrator)(Identity(), (oCore, oCore)) + utriFuncCoreType = (typeof∘buildCoreIntegrator)(OneBodyIntegral{D}(), Identity(), (oCore, oCore)) utriFuncs = Memory{VariedNormCore{T, D, 2, utriFuncCoreType}}(undef, 0) end diff --git a/src/Types.jl b/src/Types.jl index 3d288a86..41ba7b1e 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -22,9 +22,12 @@ abstract type FunctionModifier <: CompositeFunction end # Modify a function abstract type FunctionComposer <: CompositeFunction end # Compose only functions together abstract type FunctionalStruct <: CompositeFunction end # A struct with defined methods +abstract type StatefulFunction{T} <: CompositeFunction end +abstract type StatefulVariable{T} <: CompositeFunction end + abstract type TypedEvaluator{T, F} <: Evaluator{F} end -abstract type SpatialIntegralCache{T, D} <: QueryBox{T} end +abstract type SpatialProcessCache{T, D} <: QueryBox{T} end abstract type IntegralData{T, S} <: QueryBox{T} end abstract type ViewedObject{T, P} <: QueryBox{T} end abstract type CustomCache{T} <: QueryBox{T} end @@ -43,7 +46,7 @@ abstract type JoinedOperator{J} <: FunctionComposer end abstract type CompositePointer <: ConfigBox end abstract type StructuredType <: ConfigBox end -abstract type IntegralProcessCache{T, D} <: SpatialIntegralCache{T, D} end +abstract type IntegralProcessCache{T, D} <: SpatialProcessCache{T, D} end abstract type ActivePointer <: CompositePointer end abstract type StaticPointer{P<:ActivePointer} <: CompositePointer end From 89be2f273f76e75b02b6dd2306d89db6b71bc05a Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 28 Jan 2025 00:08:45 -0500 Subject: [PATCH 28/29] Formatting code. --- src/Integrals/Framework.jl | 46 ++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index cccd2b33..1027ec69 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -71,7 +71,7 @@ end # ab::Dict{NTuple{2, Int}, T} # end -#!!! Just changed the signature + struct PrimOrbCoreIntegralCache{T, D, S<:MultiBodyIntegral{D}, F<:DirectOperator, I<:IntegralData{T, S}, B<:PrimitiveOrbCore{T, D} } <: SpatialProcessCache{T, D} @@ -88,6 +88,7 @@ const OverlapCoreCache{T, D, I<:IntegralData{T, OneBodyIntegral{D}}, B<:PrimitiveOrbCore{T, D}} = POrb1BCoreICache{T, D, Identity, I, B} + function setIntegralData!(ints::OneBodyFullCoreIntegrals{T}, pair::Pair{Tuple{Int}, Tuple{T}}) where {T} setindex!(getfield(ints, OneBodyIdxSymDict[true ]), pair.second, pair.first) @@ -110,15 +111,6 @@ function getPrimCoreIntData(cache::OneBodyFullCoreIntegrals, idx::Tuple{Int}) end -isHermitian(::PrimitiveOrbCore{T, D}, ::DirectOperator, - ::PrimitiveOrbCore{T, D}) where {T, D} = -false - -isHermitian(::PrimitiveOrbCore{T, D}, ::Identity, - ::PrimitiveOrbCore{T, D}) where {T, D} = -true - - struct Bar <: Any end const OrbCoreIntLayoutAllOrbs{T, D, N} = NonEmptyTuple{PrimitiveOrbCore{T, D}, N} @@ -141,13 +133,10 @@ const TwoBodyOrbCoreIntLayout{T, D} = Union{ OrbCoreIntLayoutOrbBar3{T, D}, OrbCoreIntLayoutOrbBar4{T, D}, } -# const OrbCoreIntLayout{T, D, U<:Union{Bar, PrimitiveOrbCore{T, D}}, -# PrimitiveOrbCore{T, D}<:V<:U, N} = -# Tuple{PrimitiveOrbCore{T, D}, Vararg{V, N}} - const OrbCoreIntLayout{T, D} = Union{OneBodyOrbCoreIntLayout{T, D}, TwoBodyOrbCoreIntLayout{T, D}} + struct OrbCoreIntConfig{T, D, S<:MultiBodyIntegral{D}, P<:OrbCoreIntLayout{T, D}} <: StructuredType @@ -168,8 +157,10 @@ const OneBodyOrbCoreIntConfig{T, D, P<:OneBodyOrbCoreIntLayout{T, D}} = const TwoBodyOrbCoreIntConfig{T, D, P<:TwoBodyOrbCoreIntLayout{T, D}} = OrbCoreIntConfig{T, D, TwoBodyIntegral{D}, P} + abstract type OrbIntegralCompCacheBox{T, D, S<:MultiBodyIntegral{D}} <: QueryBox{T} end + struct OrbIntCompCacheBox{T, D, S<:MultiBodyIntegral{D}} <: OrbIntegralCompCacheBox{T, D, S} dict::Dict{OrbCoreIntConfig{T, D, S}, OrbIntegralComputeCache{T, D, S}} @@ -179,14 +170,17 @@ struct OrbIntCompCacheBox{T, D, S<:MultiBodyIntegral{D}} <: OrbIntegralCompCache end end + struct OrbIntNullCacheBox{T, D, S<:MultiBodyIntegral{D}} <: OrbIntegralCompCacheBox{T, D, S} OrbIntCompCacheBox(::S, ::Type{T}) where {D, S<:MultiBodyIntegral{D}, T} = new{T, D, S}() end + abstract type CachedComputeSpatialIntegral{T, D} <: StatefulFunction{T} end + struct CachedComputeOrbCoreIntegral{T, D, F<:DirectOperator, C<:OrbIntegralCompCacheBox{T, D} } <: CachedComputeSpatialIntegral{T, D} @@ -277,6 +271,7 @@ function (f::OrbCoreIntegralComputeConfig{T, D, S} end end + # One-Body (i|O|j) hermiticity across O getHermiticity(::PrimitiveOrbCore{T, D}, ::DirectOperator, ::PrimitiveOrbCore{T, D}) where {T, D} = @@ -291,6 +286,7 @@ getHermiticity(::N12Tuple{PrimitiveOrbCore{T, D}}, ::Identity, ::N12Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = true + function genCoreIntTuple(integrator::Orb1BCCIntComputeConfigUnion{T, D}, oDataTuple::Tuple{OrbCoreData{T, D}}) where {T, D} (integrator(oDataTuple),) @@ -340,6 +336,7 @@ end # oneBasedIds::NTuple{4, Int}) where {T, D} # end + function sortTensorIndex((i, j)::NTuple{2, Int}) if i > j (j, i) @@ -504,7 +501,6 @@ function updateOrbCache!(orbCache::PrimOrbCoreCache{T, D}, end end - function updateOrbCache!(orbCache::PrimOrbCoreCache{T, D}, orbMarkerCache::OrbCoreMarkerDict, orbData::OrbCoreData{T, D}) where {T, D} @@ -706,7 +702,7 @@ end function buildIndexedOrbWeightsCore!(primOrbPtrs::AbstractVector{Int}, - primOrbWeight::AbstractVector{T}) where {T} + primOrbWeight::AbstractVector{T}) where {T} nPtrs = length(primOrbWeight) list = Memory{Pair{Int, T}}(undef, nPtrs) for i in 1:nPtrs @@ -717,10 +713,10 @@ end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCoreCache{T, D}, - orbs::OrbitalCollection{T, D}, - normIdxList::BasisIndexList, - intIdxList::BasisIndexList=normIdxList) where {T, D} + normCache::OverlapCoreCache{T, D}, + orbs::OrbitalCollection{T, D}, + normIdxList::BasisIndexList, + intIdxList::BasisIndexList=normIdxList) where {T, D} i = 0 map(orbs) do orb iRange = getBasisIndexRange(intIdxList, (i+=1)) @@ -732,10 +728,10 @@ function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCoreCache{T, D}, - orb::FrameworkOrb{T, D}, - normIdxList::BasisIndexList, - intIdxList::BasisIndexList=normIdxList) where {T, D} + normCache::OverlapCoreCache{T, D}, + orb::FrameworkOrb{T, D}, + normIdxList::BasisIndexList, + intIdxList::BasisIndexList=normIdxList) where {T, D} orbWeight = buildOrbWeight!(paramCache, normCache, orb, normIdxList.index) buildIndexedOrbWeightsCore!(intIdxList.index, orbWeight) end @@ -864,6 +860,7 @@ function computeIntTensor!(intCache::PrimOrbCoreIntegralCache{T, D}, buildIntegralTensor(intCache, ws) end + function computeIntTensor(::OneBodyIntegral{D}, op::F, basisSet::AbstractVector{<:FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), @@ -905,6 +902,7 @@ function computeIntegral!(intCache::PrimOrbCoreIntegralCache{T, D}, buildIntegralEntries(intCache, ws) |> first end + function computeIntegral(::OneBodyIntegral{D}, op::DirectOperator, (bf1,)::Tuple{FrameworkOrb{T, D}}; paramCache::DimSpanDataCacheBox{T}=DimSpanDataCacheBox(T), From 66b1f1d723903b062a3790f80c064faeb9fa2a51 Mon Sep 17 00:00:00 2001 From: Weishi Wang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 28 Jan 2025 01:13:40 -0500 Subject: [PATCH 29/29] Replaced generated functions with dict getters for generating integral caches. --- src/Integrals/Engines/GaussianOrbitals.jl | 52 +++++++++-------------- src/Integrals/Framework.jl | 29 +++++++------ 2 files changed, 38 insertions(+), 43 deletions(-) diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 6091f7ea..13ab0962 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -136,7 +136,8 @@ function overlapPGTOcore(xpn::T, ang::NonEmptyTuple{Int}) where {T} end end -function lazyOverlapPGTO(cache::LRU{TupleOf5T2Int{T}, T}, input::TupleOf5T2Int{T}) where {T} +function lazyOverlapPGTO!(cache::LRU{TupleOf5T2Int{T}, T}, + input::TupleOf5T2Int{T}) where {T} res = get(cache, input, nothing) # Fewer allocations than using `get!` if res === nothing res = overlapPGTOcore(input) @@ -145,9 +146,9 @@ function lazyOverlapPGTO(cache::LRU{TupleOf5T2Int{T}, T}, input::TupleOf5T2Int{T res end -lazyOverlapPGTO(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) +lazyOverlapPGTO!(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) -lazyOverlapPGTO(::ZeroAngMomCache{T}, input::TupleOf5T2Int{T}) where {T} = +lazyOverlapPGTO!(::ZeroAngMomCache{T}, input::TupleOf5T2Int{T}) where {T} = gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) @@ -157,8 +158,8 @@ const GTOrbOverlapAxialCache{T} = const GTOrbOverlapCache{T, D} = AxialOneBodyIntCompCache{T, D, <:NTuple{D, GTOrbOverlapAxialCache{T}}} -function overlapPGTO(data::GaussProductInfo{T, D}, - cache::GTOrbOverlapCache{T, D}) where {T, D} +function overlapPGTO!(cache::GTOrbOverlapCache{T, D}, data::GaussProductInfo{T, D}, + ) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen cenM = data.cen @@ -170,7 +171,7 @@ function overlapPGTO(data::GaussProductInfo{T, D}, angR = data.rhs.ang mapreduce(*, cache.axis, cenL, cenR, cenM, angL, angR) do cache, xL, xR, xM, iL, iR - lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) + lazyOverlapPGTO!(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR)) end end @@ -213,10 +214,10 @@ end function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; - cache::GTOrbOverlapCache{T, D}= + cache!Self::GTOrbOverlapCache{T, D}= AxialOneBodyIntCompCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) - overlapPGTO(data, cache) + overlapPGTO!(cache!Self, data) end @@ -227,23 +228,12 @@ end const DefaultPGTOrbOverlapCacheSizeLimit = 128 -# function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} -# ntuple(_->LRU{TupleOf5T2Int{T}, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit), Val(D)) -# end - -# function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} -# ntuple(_->ZeroAngMomCache{T}(), Val(D)) -# end - -@generated function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} - maxsize = DefaultPGTOrbOverlapCacheSizeLimit - caches = ntuple(_->LRU{TupleOf5T2Int{T}, T}(; maxsize), Val(D)) - return :( $caches ) +function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{L}) where {T, D, L} + ntuple(_->LRU{TupleOf5T2Int{T}, T}(maxsize=DefaultPGTOrbOverlapCacheSizeLimit), Val(D)) end -@generated function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} - caches = ntuple(_->ZeroAngMomCache{T}(), Val(D)) - return :( $caches ) +function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} + ntuple(_->ZeroAngMomCache{T}(), Val(D)) end genGTOrbIntCompCache(::MonomialMul{T, D, L}, @@ -281,8 +271,8 @@ function computeMultiMomentGTO(op::MonomialMul{T, D}, end -function computeMultiMomentGTO(op::MonomialMul{T, D}, data::GaussProductInfo{T, D}, - cache::GTOrbOverlapCache{T, D}) where {T, D} +function computeMultiMomentGTO!(op::MonomialMul{T, D}, cache::GTOrbOverlapCache{T, D}, + data::GaussProductInfo{T, D}) where {T, D} cenL = data.lhs.cen cenR = data.rhs.cen @@ -300,14 +290,14 @@ function computeMultiMomentGTO(op::MonomialMul{T, D}, data::GaussProductInfo{T, m = iszero(dx) ? 0 : n mapreduce(+, 0:m) do k binomial(n, k) * dx^k * - lazyOverlapPGTO(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR+n-k)) + lazyOverlapPGTO!(cache, (xpnS, xpnRatio, xM-xL, xM-xR, iL, iR+n-k)) end end end -computeMultiMomentGTO(::MonomialMul{T, D, 0}, data::GaussProductInfo{T, D}, - cache::GTOrbOverlapCache{T, D}) where {T, D} = -overlapPGTO(data, cache) +computeMultiMomentGTO!(::MonomialMul{T, D, 0}, cache::GTOrbOverlapCache{T, D}, + data::GaussProductInfo{T, D}) where {T, D} = +overlapPGTO!(cache, data) struct PrimGTOrbMultiMoment{T, D, L, B<:N12Tuple{PrimGaussOrbConfig{T, D}} @@ -331,10 +321,10 @@ end function (f::MultiMomentGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, pars2::FilteredVecOfArr{T}; - cache::GTOrbOverlapCache{T, D}= + cache!Self::GTOrbOverlapCache{T, D}= AxialOneBodyIntCompCache(T, Val(D))) where {T, D} data = GaussProductInfo(f.basis, (pars1, pars2)) - computeMultiMomentGTO(f.op, data, cache) + computeMultiMomentGTO!(f.op, cache!Self, data) end diff --git a/src/Integrals/Framework.jl b/src/Integrals/Framework.jl index 1027ec69..14e34d36 100644 --- a/src/Integrals/Framework.jl +++ b/src/Integrals/Framework.jl @@ -242,18 +242,23 @@ function (f::DirectComputeOrbCoreIntegral{T, D, S} fCore(parData...) end +const OrbCoreIntegralComputeCacheTypes{T, D} = Union{ + OrbCoreIntConfig{T, D, <:MultiBodyIntegral{D}, TetraTupleUnion{PrimGTOcore{T, D}}} +} -function prepareOneBodyIntCompCache!(f::ReusedComputeOrbCoreIntegral{T, D}, - orbs::TetraTupleUnion{PrimGTOcore{T, D}}) where {T, D} - # get!(f.cache, typeof( (op, orbs...) )) do - genGTOrbIntCompCache(f.operator, orbs) - # end +function prepareOneBodyIntCompCache!(f::ReusedComputeOrbCoreIntegral{T, D, S}, + orbs::TetraTupleUnion{PrimitiveOrbCore{T, D}}) where + {T, D, S<:MultiBodyIntegral{D}} + key = OrbCoreIntConfig(S(), orbs) + if key isa OrbCoreIntegralComputeCacheTypes + get!(f.cache.dict, key) do + genGTOrbIntCompCache(f.operator, orbs) + end + else + nothing + end end -prepareOneBodyIntCompCache!(::ReusedComputeOrbCoreIntegral{T, D}, - ::TetraTupleUnion{PrimitiveOrbCore{T, D}}) where {T, D} = -nothing - prepareOneBodyIntCompCache!(::ReusedComputeOrbCoreIntegral{T, D, OneBodyIntegral{D}}, ::Tuple{PrimGTOcore{T, D}}) where {T, D} = nothing @@ -263,11 +268,11 @@ function (f::OrbCoreIntegralComputeConfig{T, D, S} orbData = first.(data) parData = last.(data) fCore = buildCoreIntegrator(S(), f.operator, orbData) - intCacheSector = prepareOneBodyIntCompCache!(f, orbData) - if intCacheSector === nothing + intSectorCache = prepareOneBodyIntCompCache!(f, orbData) + if intSectorCache === nothing fCore(parData...) else - fCore(parData..., cache=intCacheSector) + fCore(parData..., cache!Self=intSectorCache) end end