diff --git a/src/Integrals/Engines/GaussianOrbitals.jl b/src/Integrals/Engines/GaussianOrbitals.jl index 9d1fb6d5..13ab0962 100644 --- a/src/Integrals/Engines/GaussianOrbitals.jl +++ b/src/Integrals/Engines/GaussianOrbitals.jl @@ -1,96 +1,184 @@ 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 = BigInt(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) - T(πPowers[:p0d5]) / sqrt(2α) * factor +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} + 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 -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) +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 + +function gaussProdCore3(dxLR::T, xpnProdOverSum::T) where {T<:Real} + exp(- xpnProdOverSum * dxLR^2) +end + +function gaussProdCore4(xpnSum::T, xpnRatio::T, dxLR::T) where {T<:Real} + gaussProdCore1(xpnSum, 0) * gaussProdCore3(dxLR, xpnRatio) +end + + +struct PrimGaussOrbInfo{T, D} + cen::NTuple{D, T} + xpn::T + ang::NTuple{D, Int} + + 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 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)) +struct GaussProductInfo{T, D} + lhs::PrimGaussOrbInfo{T, D} + rhs::PrimGaussOrbInfo{T, D} + cen::NTuple{D, T} + xpn::T + + 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}(oData1, oData2, cenNew, xpnSum) + end end -function (cPair::CenPair{T})(xPair::XpnPair{T}) where {T} - (xPair.left .* cPair.left .+ xPair.right .* cPair.right) ./ xPair.sum +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 AngPair{D} - left::NTuple{D, Int} - right::NTuple{D, Int} - sum::NTuple{D, Int} +struct NullCache{T} <: CustomCache{T} end + +const TupleOf5T2Int{T} = Tuple{T, T, T, T, Int, Int} + +struct ZeroAngMomCache{T} <: CustomCache{T} end - AngPair(l::NonEmptyTuple{Int, D}, r::NonEmptyTuple{Int, D}) where {D} = - new{D+1}(l, r, l .+ r) +const CoreIntCacheBox{T} = Union{CustomCache{T}, LRU{<:Any, T}} + +abstract type OrbIntegralComputeCache{T, D, S<:MultiBodyIntegral{D} + } <: IntegralProcessCache{T, D} end + +const Orb1BIntegralComputeCache{T, D} = OrbIntegralComputeCache{T, D, OneBodyIntegral{D}} + +struct AxialOneBodyIntCompCache{T, D, F<:NTuple{D, CoreIntCacheBox{T}} + } <: Orb1BIntegralComputeCache{T, D} + axis::F + + AxialOneBodyIntCompCache(axis::NonEmptyTuple{CoreIntCacheBox{T}, D}) where {T, D} = + new{T, D+1, typeof(axis)}(axis) end +AxialOneBodyIntCompCache(::Type{T}, ::Val{D}) where {T, D} = +AxialOneBodyIntCompCache(ntuple( _->NullCache{T}(), Val(D) )) + -function gaussProdCore1(cPair::CenPair{T}, xPair::XpnPair{T}) where {T<:Real} - if iszero(cPair.dist) - T(1) +function overlapPGTOcore(input::TupleOf5T2Int{T}) where {T} + xpnSum, xpnRatio, xML, xMR, iL, iR = input + dxLR = xMR - xML + jRange = 0:((iL + iR) ÷ 2) + + if isequal(dxLR, zero(T)) + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(iL, iR, 2j) + end else - exp(- xPair.prod / xPair.sum * cPair.dist^2) + mapreduce(+, jRange) do j + gaussProdCore1(xpnSum, j) * gaussProdCore2(xML, xMR, iL, iR, 2j) + end * gaussProdCore3(dxLR, xpnRatio) end end -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 overlapPGTOcore(xpn::T, ang::NonEmptyTuple{Int}) where {T} + mapreduce(*, ang) do i + gaussProdCore1(2xpn, i) + end +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) + setindex!(cache, res, input) end res end +lazyOverlapPGTO!(::NullCache{T}, input::TupleOf5T2Int{T}) where {T} = overlapPGTOcore(input) -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) - end - end * gaussProdCore1(cPair, xPair) -end +lazyOverlapPGTO!(::ZeroAngMomCache{T}, input::TupleOf5T2Int{T}) where {T} = +gaussProdCore4(input[begin], input[begin+1], input[begin+3]-input[begin+2]) -function overlapPGTO(xpn::T, ang::NTuple{D, Int}) where {T, D} - mapreduce(*, ang) do i - polyGaussFuncSquaredNorm(xpn, i) +const GTOrbOverlapAxialCache{T} = + Union{NullCache{T}, ZeroAngMomCache{T}, LRU{TupleOf5T2Int{T}, T}} + +const GTOrbOverlapCache{T, D} = + AxialOneBodyIntCompCache{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 + + xpnS = data.xpn + xpnRatio = data.lhs.xpn * data.rhs.xpn / xpnS + + angL = data.lhs.ang + 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)) end 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) @@ -102,49 +190,145 @@ function getCenCoordPtr(orb::PrimitiveOrbCore) orb.f.dress.select end + function preparePGTOconfig(orb::PrimGTOcore) cenIds = getCenCoordPtr(orb) xpnIdx = getExponentPtr(orb) - ang = getAngularNums(orb) - cenIds, xpnIdx, ang + ang = getAngMomentum(orb).tuple + PrimGaussOrbInfo(cenIds, xpnIdx, ang) end -struct OverlapGTOrbSelf{T, D} <: OrbitalIntegrator{T, D} - xpn::FlatPSetInnerPtr{T} - ang::NTuple{D, Int} +struct PrimGTOrbOverlap{T, D, B<:N12Tuple{PrimGaussOrbConfig{T, D}} + } <: OrbitalIntegrator{T, D} + basis::B end +const OverlapGTOrbSelf{T, D} = PrimGTOrbOverlap{T, D, Tuple{ PrimGaussOrbConfig{T, D} }} +const OverlapGTOrbPair{T, D} = PrimGTOrbOverlap{T, D, NTuple{ 2, PrimGaussOrbConfig{T, D} }} + function (f::OverlapGTOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T} - xpnVal = getField(pars, f.xpn) - overlapPGTO(xpnVal, f.ang) + sector = first(f.basis) + xpnVal = getField(pars, sector.xpn) + overlapPGTOcore(xpnVal, sector.ang) end -function genGTOrbOverlapFunc((orb,)::Tuple{PrimGTOcore{T, D}}) where {T, D} - _, xpnIdx, ang = preparePGTOconfig(orb) - OverlapGTOrbSelf(xpnIdx, ang) +function (f::OverlapGTOrbPair{T, D})(pars1::FilteredVecOfArr{T}, + pars2::FilteredVecOfArr{T}; + cache!Self::GTOrbOverlapCache{T, D}= + AxialOneBodyIntCompCache(T, Val(D))) where {T, D} + data = GaussProductInfo(f.basis, (pars1, pars2)) + overlapPGTO!(cache!Self, data) end -struct OverlapGTOrbPair{T, D} <: OrbitalIntegrator{T, D} - cen::NTuple{2, NTuple{ D, FlatPSetInnerPtr{T} }} - xpn::NTuple{2, FlatPSetInnerPtr{T}} - ang::NTuple{2, NTuple{D, Int}} + +function genGTOrbOverlapFunc(orbCoreData::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + PrimGTOrbOverlap(preparePGTOconfig.(orbCoreData)) 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...)) + +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 genGTOrbOverlapFunc(orbs::NTuple{2, PrimGTOcore{T, D}}) where {T, D} - configs = preparePGTOconfig.(orbs) - OverlapGTOrbPair(getindex.(configs, 1), getindex.(configs, 2), getindex.(configs, 3)) +function genPrimGTOrbOverlapCache(::Type{T}, ::Val{D}, ::Val{0}) where {T, D} + ntuple(_->ZeroAngMomCache{T}(), Val(D)) 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)) |> AxialOneBodyIntCompCache + +genGTOrbIntCompCache(::Identity, + ::Tuple{TypedPrimGTOcore{T, D, L1}, TypedPrimGTOcore{T, D, L2}}) where + {T, D, L1, L2} = +genPrimGTOrbOverlapCache(T, Val(D), Val(L1+L2)) |> AxialOneBodyIntCompCache + 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}, cache::GTOrbOverlapCache{T, D}, + data::GaussProductInfo{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}, cache::GTOrbOverlapCache{T, D}, + data::GaussProductInfo{T, D}) where {T, D} = +overlapPGTO!(cache, data) + + +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!Self::GTOrbOverlapCache{T, D}= + AxialOneBodyIntCompCache(T, Val(D))) where {T, D} + data = GaussProductInfo(f.basis, (pars1, pars2)) + computeMultiMomentGTO!(f.op, cache!Self, data) +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/Engines/Numerical.jl b/src/Integrals/Engines/Numerical.jl index 1b9710e8..b3df99ec 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(StableBinary(*, 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{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{D}(), 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,)) + buildCoreIntegrator(OneBodyIntegral{D}(), Identity(), (o,)) +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, D})(pVal::FilteredVecOfArr{T}) where {T, D} + numericalIntegrate(OneBodyIntegral{D}(), f.op, f.basis, (pVal,)) +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 +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/Framework.jl b/src/Integrals/Framework.jl index 9a43b6ba..14e34d36 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}} @@ -33,9 +30,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 @@ -62,86 +59,288 @@ 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, S<:MultiBodyIntegral{D}, F<:DirectOperator, + I<:IntegralData{T, S}, B<:PrimitiveOrbCore{T, D} + } <: SpatialProcessCache{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 POrb1BCoreICache{T, D, F<:DirectOperator, I<:IntegralData{T, OneBodyIntegral{D}}, + B<:PrimitiveOrbCore{T, D}} = + PrimOrbCoreIntegralCache{T, D, OneBodyIntegral{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}} = + POrb1BCoreICache{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 -genOneBodyCoreIntegrator(::Identity, orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} = -genGTOrbOverlapFunc(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} = + 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} + +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}} + + 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 + + +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} + operator::F + cache::C + + 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 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 + +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 -isHermitian(::PrimitiveOrbCore{T, D}, ::DirectOperator, - ::PrimitiveOrbCore{T, D}) where {T, D} = +function buildCoreIntegrator(::OneBodyIntegral{D}, op::MonomialMul{T, D}, + orbs::N12Tuple{PrimGTOcore{T, D}}) where {T, D} + genGTOrbMultiMomentFunc(op, orbs) +end + + +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 + +const OrbCoreIntegralComputeCacheTypes{T, D} = Union{ + OrbCoreIntConfig{T, D, <:MultiBodyIntegral{D}, TetraTupleUnion{PrimGTOcore{T, D}}} +} + +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, 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 = buildCoreIntegrator(S(), f.operator, orbData) + intSectorCache = prepareOneBodyIntCompCache!(f, orbData) + if intSectorCache === nothing + fCore(parData...) + else + fCore(parData..., cache!Self=intSectorCache) + end +end + + +# One-Body (i|O|j) hermiticity across O +getHermiticity(::PrimitiveOrbCore{T, D}, ::DirectOperator, + ::PrimitiveOrbCore{T, D}) where {T, D} = false -isHermitian(::PrimitiveOrbCore{T, D}, ::Identity, - ::PrimitiveOrbCore{T, D}) where {T, D} = +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 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) +function genCoreIntTuple(integrator::Orb1BCCIntComputeConfigUnion{T, D}, + oDataTuple::Tuple{OrbCoreData{T, D}}) where {T, D} + (integrator(oDataTuple),) end -function evalOneBodyPrimCoreIntegral(op::DirectOperator, - (orb, pars)::OrbCoreData{T, D}) where {T, D} - f = ReturnTyped(genOneBodyCoreIntegrator(op, (orb,)), T) - f(pars) +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 @@ -164,61 +363,42 @@ end 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,) -end - -function genOneBodyIntDataPairsCore(op::DirectOperator, - oDataPair::NTuple{2, OrbCoreDataSeq{T}}, - oneBasedIdxPair::NTuple{2, Int}) where {T} - orbPars1, orbPars2 = map(oDataPair, oneBasedIdxPair) do data, idx - data[begin+idx-1] - end - ijVal = evalOneBodyPrimCoreIntegral(op, orbPars1, orbPars2) - jiVal = if isHermitian(first(orbPars1), op, first(orbPars2)) - ijVal' - else - evalOneBodyPrimCoreIntegral(op, orbPars2, orbPars1) - end - (ijVal, jiVal) -end - - -function genOneBodyIntDataPairs(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T, D}}, +function genOneBodyIntDataPairs(integrator::CachedComputeOrbCoreIntegral{T, D}, + (oDataSeq,)::Tuple{OrbCoreDataSeq{T, D}}, (indexOffset,)::Tuple{Int}=(0,)) where {T, D} - nOrbs = length(oData) - offset = indexOffset + firstindex(oData) - 1 + nOrbs = length(oDataSeq) + firstIdx = firstindex(oDataSeq) + offset = indexOffset + firstIdx - 1 pairs1 = map(1:nOrbs) do i - iiVal = genOneBodyIntDataPairsCore(op, (oData,), (i,)) + 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) + 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}}, +function genOneBodyIntDataPairs(integrator::CachedComputeOrbCoreIntegral{T, D}, + oDataSeqPair::NTuple{2, OrbCoreDataSeq{T, D}}, indexOffsets::NTuple{2, Int}) where {T, D} - mnOrbs = length.(oDataPair) - offsets = indexOffsets .+ firstindex.(oDataPair) .- 1 + 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) + oDataPair = getindex.(oDataSeqPair, firstIdx .+ ijPair) + ijValPair = genCoreIntTuple(integrator, oDataPair) idxPairNew => ijValPair end |> vec end @@ -238,39 +418,41 @@ function setTensorEntries!(tensor::AbstractMatrix{T}, valPair::NTuple{2, T}, end -function computeOneBodyPrimCoreIntTensor(op::DirectOperator, - (oData,)::Tuple{OrbCoreDataSeq{T}}) where {T} - 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,)) + 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)) + 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}}) 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 - ijVal = evalOneBodyPrimCoreIntegral(op, oData1[begin+i-1], oData2[begin+j-1]) + oDataPair = (oData1[begin+i-1], oData2[begin+j-1]) + ijVal = integrator(oDataPair) res[begin+i-1, begin+j-1] = ijVal end res end -struct BasisIndexList +struct BasisIndexList <: QueryBox{Int} index::Memory{Int} endpoint::Memory{Int} @@ -324,7 +506,6 @@ function updateOrbCache!(orbCache::PrimOrbCoreCache{T, D}, end end - function updateOrbCache!(orbCache::PrimOrbCoreCache{T, D}, orbMarkerCache::OrbCoreMarkerDict, orbData::OrbCoreData{T, D}) where {T, D} @@ -380,32 +561,22 @@ function indexCacheOrbData!(targetCache::PrimOrbCoreCache{T, D}, end -function cachePrimCoreIntegrals!(intCache::IntegralCache{T, D}, - paramCache::DimSpanDataCacheBox{T}, - orbMarkerCache::OrbCoreMarkerDict, - orbs::OrbitalCollection{T, D}) where {T, D} - orbCache = intCache.basis - oldMaxIdx = lastindex(orbCache.list) - orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbs) - updatePrimCoreIntCache!(intCache, oldMaxIdx+1) - orbIdxList -end - -function cachePrimCoreIntegrals!(intCache::IntegralCache{T, D}, +function cachePrimCoreIntegrals!(intCache::PrimOrbCoreIntegralCache{T, D}, paramCache::DimSpanDataCacheBox{T}, orbMarkerCache::OrbCoreMarkerDict, - orb::FrameworkOrb{T, D}) where {T, D} + orbData::Union{OrbitalCollection{T, D}, FrameworkOrb{T, D}} + ) where {T, D} orbCache = intCache.basis oldMaxIdx = lastindex(orbCache.list) - orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orb) + orbIdxList = indexCacheOrbData!(orbCache, paramCache, orbMarkerCache, orbData) updatePrimCoreIntCache!(intCache, oldMaxIdx+1) orbIdxList end -function cachePrimCoreIntegrals!(targetIntCache::IntegralCache{T, D}, - sourceIntCache::IntegralCache{T, D}, +function cachePrimCoreIntegrals!(targetIntCache::C, sourceIntCache::C, orbMarkerCache::OrbCoreMarkerDict, - sourceOrbList::BasisIndexList) where {T, D} + sourceOrbList::BasisIndexList) where + {C<:PrimOrbCoreIntegralCache} tOrbCache = targetIntCache.basis oldMaxIdx = lastindex(tOrbCache.list) sOrbCache = sourceIntCache.basis @@ -415,45 +586,48 @@ function cachePrimCoreIntegrals!(targetIntCache::IntegralCache{T, D}, end -function updateIntCacheCore!(op::DirectOperator, ints::CompleteOneBodyIntegrals{T}, +function updateIntCacheCore!(integrator::CachedComputeOrbCoreIntegral{T, D}, + ints::OneBodyFullCoreIntegrals{T, D}, basis::Tuple{OrbCoreDataSeq{T, D}}, offset::Tuple{Int}) where {T, D} - pairs1, pairs2 = genOneBodyIntDataPairs(op, basis, offset) + 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::CompleteOneBodyIntegrals{T}, +function updateIntCacheCore!(integrator::CachedComputeOrbCoreIntegral{T, D}, + ints::OneBodyFullCoreIntegrals{T, D}, basis::NTuple{2, OrbCoreDataSeq{T, D}}, offset::NTuple{2, Int}) where {T, D} - pairs2 = genOneBodyIntDataPairs(op, basis, offset) + pairs2 = genOneBodyIntDataPairs(integrator, basis, offset) foreach(p->setIntegralData!(ints, p), pairs2) ints end -function updatePrimCoreIntCache!(cache::IntegralCache, startIdx::Int) - 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,)) + 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,)) - updateIntCacheCore!(op, ints, (oldBasis, newBasis,), (0, boundary)) + updateIntCacheCore!(integrator, ints, (newBasis,), (boundary,)) + updateIntCacheCore!(integrator, ints, (oldBasis, newBasis,), (0, boundary)) 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) @@ -461,13 +635,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 @@ -478,7 +652,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) @@ -508,7 +683,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) @@ -532,7 +707,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 @@ -543,10 +718,10 @@ end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCache{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)) @@ -558,17 +733,17 @@ function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, end function buildIndexedOrbWeights!(paramCache::DimSpanDataCacheBox{T}, - normCache::OverlapCache{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 -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} @@ -582,7 +757,7 @@ function prepareIntegralConfig!(intCache::IntegralCache{T, D}, end -function buildIntegralEntries(intCache::OneBodyIntCache{T}, +function buildIntegralEntries(intCache::POrb1BCoreICache{T}, (intWeights,)::Tuple{IndexedWeights{T}}) where {T} idxList = intWeights.list len = length(idxList) @@ -600,7 +775,7 @@ function buildIntegralEntries(intCache::OneBodyIntCache{T}, (res,) # ([1|O|1],) end -function buildIntegralEntries(intCache::OneBodyIntCache{T}, +function buildIntegralEntries(intCache::POrb1BCoreICache{T}, intWeightPair::NTuple{2, IndexedWeights{T}}) where {T} list1, list2 = getfield.(intWeightPair, :list) intValCache = intCache.data @@ -618,7 +793,7 @@ function buildIntegralEntries(intCache::OneBodyIntCache{T}, end -function buildIntegralTensor(intCache::OneBodyIntCache{T}, +function buildIntegralTensor(intCache::POrb1BCoreICache{T}, intWeights::AbstractVector{IndexedWeights{T}}) where {T} nOrbs = length(intWeights) res = ShapedMemory{T}(undef, (nOrbs, nOrbs)) @@ -657,30 +832,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 @@ -689,7 +865,8 @@ 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 @@ -704,7 +881,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) @@ -720,8 +897,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 @@ -730,7 +907,8 @@ 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(), @@ -740,57 +918,62 @@ function computeIntegral(::OneBodyIntegral, op::DirectOperator, computeIntegral!(iCache, nCache, (bf1,); paramCache, basisMarkerCache) else coreData, w = decomposeOrb(bf1; paramCache, markerCache=basisMarkerCache) - tensor = computeOneBodyPrimCoreIntTensor(op, (coreData,)) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) + tensor = genOneBodyPrimCoreIntTensor(integrator, (coreData,)) 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) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) 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 = genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) if isempty(coreDataM) block4 else - block1 = computeOneBodyPrimCoreIntTensor(op, (coreDataM,)) - block2 = computeOneBodyPrimCoreIntTensor(op, (coreData1, coreDataM)) - block3 = computeOneBodyPrimCoreIntTensor(op, (coreDataM, coreData2)) + 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, (coreData1, coreData2)) + genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) 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(), lazyCompute::Bool=false) where {T, D} + op = Identity() bf1, bf2 = bfPair + if lazyCompute if bf1 === bf2 - computeIntegral(OneBodyIntegral(), Identity(), (bf1,); paramCache, + computeIntegral(OneBodyIntegral{D}(), op, (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) + integrator = CachedComputeOrbCoreIntegral(Val(true), OneBodyIntegral{D}(), op, T) + tensor = genOneBodyPrimCoreIntTensor(integrator, (oData1, oData2)) dot(w1, tensor, w2) end end diff --git a/src/Integrals/Interface.jl b/src/Integrals/Interface.jl index 8c3b051b..02453098 100644 --- a/src/Integrals/Interface.jl +++ b/src/Integrals/Interface.jl @@ -4,11 +4,29 @@ 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) \ No newline at end of file +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}, + 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{D}(), 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{D}(), mmOp, basisSet; paramCache) +end \ No newline at end of file diff --git a/src/Integrals/Operators.jl b/src/Integrals/Operators.jl new file mode 100644 index 00000000..0799dbd1 --- /dev/null +++ b/src/Integrals/Operators.jl @@ -0,0 +1,30 @@ +struct Identity <: DirectOperator end + +(::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 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 diff --git a/src/Mapping.jl b/src/Mapping.jl index 74315188..22fe0671 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} @@ -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 @@ -218,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 @@ -340,11 +343,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") 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..5b4ee978 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}} @@ -74,29 +74,28 @@ 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 #? 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) @@ -116,6 +115,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} @@ -212,7 +214,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 +272,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 +299,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 @@ -431,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}, @@ -465,51 +467,86 @@ 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} 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) - - diagFuncs = Memory{ReturnTyped{T}}(undef, nOrbs) - pOrbCores = Memory{PrimitiveOrbCore{T, D}}(undef, nOrbs) - pOrbScopes = Memory{AwaitFilter{FlatParamSetFilter{T}}}(undef, nOrbs) - for i in eachindex(pOrbs) - o = pOrbs[i] - oCore = o.f.apply.left - oScope = first(o.f.scope) - diagFuncs[i] = if isRenormalized(o) + oCorePtr = ChainPointer((:left, :f, :apply, :left)) + oScopePtr = ChainPointer((:left, :f, :scope)) + nIntStyle = OneBodyIntegral{D}() + + hasUnit = false + diagFuncCoreType = Union{} + utriFuncCoreType = Union{} + + diagFuncs = map(weightedOrbs) do weightedOrb + if isRenormalized(weightedOrb.left) + hasUnit = true ReturnTyped(Unit(T), T) else - nfInner = genOneBodyCoreIntegrator(Identity(), (oCore,)) - ReturnTyped(EncodeApply((Retrieve(oScope),), nfInner), T) + oCore = getField(weightedOrb, oCorePtr) + oSelect = getField(weightedOrb, oScopePtr) |> first |> Retrieve + nfCore = buildCoreIntegrator(nIntStyle, Identity(), (oCore,)) + diagFuncCoreType = typejoin(diagFuncCoreType, typeof(nfCore)) + ReturnTyped(EncodeApply((oSelect,), nfCore), T) end - pOrbCores[i] = oCore - pOrbScopes[i] = oScope end - 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) - 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 + + 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) + 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), nfCore), T) + end + + utriFuncType = eltype(utriFuncs) + if !isconcretetype(utriFuncType) + utriFuncType = getNormalizerCoreTypeUnion(T, D, Val(2), utriFuncCoreType) + end + utriFuncs = Memory{utriFuncType}(utriFuncs) + else + oCore = getField(weightedOrbs[begin], oCorePtr) + utriFuncCoreType = (typeof∘buildCoreIntegrator)(OneBodyIntegral{D}(), Identity(), (oCore, oCore)) + utriFuncs = Memory{VariedNormCore{T, D, 2, utriFuncCoreType}}(undef, 0) end - coreTransform = HermitianContract(map(itself, diagFuncs), map(itself, utriFuncs)) + coreTransform = HermitianContract(diagFuncs, utriFuncs) getWeights = map(b->ReturnTyped(b.right.f, T), weightedOrbs) NormalizeCompOrb(coreTransform, getWeights) 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 375ebfa5..41ba7b1e 100644 --- a/src/Types.jl +++ b/src/Types.jl @@ -22,9 +22,15 @@ 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 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 abstract type AbstractAmpTensor{T, O} <: JaggedOperator{T, 0, O} end abstract type AbstractAmplitude{T} <: JaggedOperator{T, 0, 0} end @@ -40,6 +46,8 @@ abstract type JoinedOperator{J} <: FunctionComposer end abstract type CompositePointer <: ConfigBox end abstract type StructuredType <: ConfigBox end +abstract type IntegralProcessCache{T, D} <: SpatialProcessCache{T, D} end + abstract type ActivePointer <: CompositePointer end abstract type StaticPointer{P<:ActivePointer} <: CompositePointer end @@ -113,6 +121,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}