diff --git a/.gitignore b/.gitignore index 29126e4..c43b7a0 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,9 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml + +# workstation files +*.nfs* + +# claude +.claude/ \ No newline at end of file diff --git a/Manifest.toml.bak-20260829 b/Manifest.toml.bak-20260829 new file mode 100644 index 0000000..bb51ad3 --- /dev/null +++ b/Manifest.toml.bak-20260829 @@ -0,0 +1,679 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.10.4" +manifest_format = "2.0" +project_hash = "ad36f52fceeec06b44451769f3349de416042b62" + +[[deps.ADTypes]] +git-tree-sha1 = "bbc22a9a08a0ef6460041086d8a7b27940ed4ffd" +uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b" +version = "1.22.0" + + [deps.ADTypes.extensions] + ADTypesChainRulesCoreExt = "ChainRulesCore" + ADTypesConstructionBaseExt = "ConstructionBase" + ADTypesEnzymeCoreExt = "EnzymeCore" + + [deps.ADTypes.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + +[[deps.Accessors]] +deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"] +git-tree-sha1 = "2eeb2c9bef11013efc6f8f97f32ee59b146b09fb" +uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" +version = "0.1.44" + + [deps.Accessors.extensions] + AxisKeysExt = "AxisKeys" + IntervalSetsExt = "IntervalSets" + LinearAlgebraExt = "LinearAlgebra" + StaticArraysExt = "StaticArrays" + StructArraysExt = "StructArrays" + TestExt = "Test" + UnitfulExt = "Unitful" + + [deps.Accessors.weakdeps] + AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5" + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" + Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[[deps.Adapt]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "28e1637322d4019ed2577cbec9268fab9b7da117" +uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +version = "4.6.0" +weakdeps = ["SparseArrays", "StaticArrays"] + + [deps.Adapt.extensions] + AdaptSparseArraysExt = "SparseArrays" + AdaptStaticArraysExt = "StaticArrays" + +[[deps.ArrayInterface]] +deps = ["Adapt", "LinearAlgebra"] +git-tree-sha1 = "3d0cabd25fab32390e3bcb82cd67e700aebd9816" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "7.25.0" + + [deps.ArrayInterface.extensions] + ArrayInterfaceAMDGPUExt = "AMDGPU" + ArrayInterfaceBandedMatricesExt = "BandedMatrices" + ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices" + ArrayInterfaceCUDAExt = "CUDA" + ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"] + ArrayInterfaceChainRulesCoreExt = "ChainRulesCore" + ArrayInterfaceChainRulesExt = "ChainRules" + ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore" + ArrayInterfaceMetalExt = "Metal" + ArrayInterfaceReverseDiffExt = "ReverseDiff" + ArrayInterfaceSparseArraysExt = "SparseArrays" + ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore" + ArrayInterfaceTrackerExt = "Tracker" + + [deps.ArrayInterface.weakdeps] + AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e" + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" + CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e" + ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2" + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + Metal = "dde4c033-4e86-420c-a63e-0dd931031962" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + +[[deps.Artifacts]] +uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" + +[[deps.AxisAlgorithms]] +deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] +git-tree-sha1 = "01b8ccb13d68535d73d2b0c23e39bd23155fb712" +uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950" +version = "1.1.0" + +[[deps.Base64]] +uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" + +[[deps.ChainRulesCore]] +deps = ["Compat", "LinearAlgebra"] +git-tree-sha1 = "12177ad6b3cad7fd50c8b3825ce24a99ad61c18f" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "1.26.1" +weakdeps = ["SparseArrays"] + + [deps.ChainRulesCore.extensions] + ChainRulesCoreSparseArraysExt = "SparseArrays" + +[[deps.CommonSolve]] +git-tree-sha1 = "dd91a10d8b8ae06e15706158eaf1a3e87e97b5f5" +uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" +version = "0.2.7" + + [deps.CommonSolve.extensions] + CommonSolveEnzymeCoreExt = "EnzymeCore" + + [deps.CommonSolve.weakdeps] + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + +[[deps.Compat]] +deps = ["TOML", "UUIDs"] +git-tree-sha1 = "9d8a54ce4b17aa5bdce0ea5c34bc5e7c340d16ad" +uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" +version = "4.18.1" +weakdeps = ["Dates", "LinearAlgebra"] + + [deps.Compat.extensions] + CompatLinearAlgebraExt = "LinearAlgebra" + +[[deps.CompilerSupportLibraries_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" +version = "1.1.1+0" + +[[deps.CompositionsBase]] +git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" +uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b" +version = "0.1.2" +weakdeps = ["InverseFunctions"] + + [deps.CompositionsBase.extensions] + CompositionsBaseInverseFunctionsExt = "InverseFunctions" + +[[deps.ConstructionBase]] +git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb" +uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9" +version = "1.6.0" + + [deps.ConstructionBase.extensions] + ConstructionBaseIntervalSetsExt = "IntervalSets" + ConstructionBaseLinearAlgebraExt = "LinearAlgebra" + ConstructionBaseStaticArraysExt = "StaticArrays" + + [deps.ConstructionBase.weakdeps] + IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" + LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.DataStructures]] +deps = ["OrderedCollections"] +git-tree-sha1 = "e86f4a2805f7f19bec5129bc9150c38208e5dc23" +uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +version = "0.19.4" + +[[deps.Dates]] +deps = ["Printf"] +uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" + +[[deps.DifferentiationInterface]] +deps = ["ADTypes", "LinearAlgebra"] +git-tree-sha1 = "2147a95a217cc8a78ec96ee03581adf129468e49" +uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" +version = "0.7.18" + + [deps.DifferentiationInterface.extensions] + DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore" + DifferentiationInterfaceDiffractorExt = "Diffractor" + DifferentiationInterfaceEnzymeExt = ["EnzymeCore", "Enzyme"] + DifferentiationInterfaceFastDifferentiationExt = "FastDifferentiation" + DifferentiationInterfaceFiniteDiffExt = "FiniteDiff" + DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences" + DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"] + DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"] + DifferentiationInterfaceGTPSAExt = "GTPSA" + DifferentiationInterfaceHyperHessiansExt = "HyperHessians" + DifferentiationInterfaceMooncakeExt = "Mooncake" + DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"] + DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"] + DifferentiationInterfaceSparseArraysExt = "SparseArrays" + DifferentiationInterfaceSparseConnectivityTracerExt = "SparseConnectivityTracer" + DifferentiationInterfaceSparseMatrixColoringsExt = "SparseMatrixColorings" + DifferentiationInterfaceStaticArraysExt = "StaticArrays" + DifferentiationInterfaceSymbolicsExt = "Symbolics" + DifferentiationInterfaceTrackerExt = "Tracker" + DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"] + + [deps.DifferentiationInterface.weakdeps] + Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5" + Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c" + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869" + FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be" + FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" + FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" + GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8" + HyperHessians = "06b494a0-c8e0-40cc-ad32-d99506a00a6c" + Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6" + PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b" + ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5" + SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" + Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" + Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[[deps.Distributed]] +deps = ["Random", "Serialization", "Sockets"] +uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" + +[[deps.DocStringExtensions]] +git-tree-sha1 = "7442a5dfe1ebb773c29cc2962a8980f47221d76c" +uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +version = "0.9.5" + +[[deps.EnumX]] +git-tree-sha1 = "c49898e8438c828577f04b92fc9368c388ac783c" +uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56" +version = "1.0.7" + +[[deps.FillArrays]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "2f979084d1e13948a3352cf64a25df6bd3b4dca3" +uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" +version = "1.16.0" + + [deps.FillArrays.extensions] + FillArraysPDMatsExt = "PDMats" + FillArraysSparseArraysExt = "SparseArrays" + FillArraysStaticArraysExt = "StaticArrays" + FillArraysStatisticsExt = "Statistics" + + [deps.FillArrays.weakdeps] + PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" + +[[deps.FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Setfield"] +git-tree-sha1 = "f7017a4f337f8df189fcce98e32b67a1298a2115" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.31.0" + + [deps.FiniteDiff.extensions] + FiniteDiffBandedMatricesExt = "BandedMatrices" + FiniteDiffBlockBandedMatricesExt = "BlockBandedMatrices" + FiniteDiffSparseArraysExt = "SparseArrays" + FiniteDiffStaticArraysExt = "StaticArrays" + + [deps.FiniteDiff.weakdeps] + BandedMatrices = "aae01518-5342-5314-be14-df237901396f" + BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0" + SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[[deps.Future]] +deps = ["Random"] +uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820" + +[[deps.HalfIntegers]] +git-tree-sha1 = "9c3149243abb5bc0bad0431d6c4fcac0f4443c7c" +uuid = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" +version = "1.6.0" + +[[deps.IOCapture]] +deps = ["Logging", "Random"] +git-tree-sha1 = "0ee181ec08df7d7c911901ea38baf16f755114dc" +uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89" +version = "1.0.0" + +[[deps.IntegerMathUtils]] +git-tree-sha1 = "4c1acff2dc6b6967e7e750633c50bc3b8d83e617" +uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235" +version = "0.1.3" + +[[deps.InteractiveUtils]] +deps = ["Markdown"] +uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" + +[[deps.Interpolations]] +deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"] +git-tree-sha1 = "65d505fa4c0d7072990d659ef3fc086eb6da8208" +uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +version = "0.16.2" + + [deps.Interpolations.extensions] + InterpolationsForwardDiffExt = "ForwardDiff" + InterpolationsUnitfulExt = "Unitful" + + [deps.Interpolations.weakdeps] + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[[deps.InverseFunctions]] +git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb" +uuid = "3587e190-3f89-42d0-90ee-14403ec27112" +version = "0.1.17" +weakdeps = ["Dates", "Test"] + + [deps.InverseFunctions.extensions] + InverseFunctionsDatesExt = "Dates" + InverseFunctionsTestExt = "Test" + +[[deps.IrrationalConstants]] +git-tree-sha1 = "b2d91fe939cae05960e760110b328288867b5758" +uuid = "92d709cd-6900-40b7-9082-c6be49f344b6" +version = "0.2.6" + +[[deps.JLLWrappers]] +deps = ["Artifacts", "Preferences"] +git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e" +uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210" +version = "1.8.0" + +[[deps.JSON]] +deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"] +git-tree-sha1 = "f76f7560267b840e492180f9899b472f30b88450" +uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +version = "1.6.0" + + [deps.JSON.extensions] + JSONArrowExt = ["ArrowTypes"] + + [deps.JSON.weakdeps] + ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd" + +[[deps.LRUCache]] +git-tree-sha1 = "5519b95a490ff5fe629c4a7aa3b3dfc9160498b3" +uuid = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +version = "1.6.2" +weakdeps = ["Serialization"] + + [deps.LRUCache.extensions] + SerializationExt = ["Serialization"] + +[[deps.Libdl]] +uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb" + +[[deps.LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Printf"] +git-tree-sha1 = "cef1ba655e8c1f65af9d96c4fffe18bf1a3a3291" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.7.1" + +[[deps.LinearAlgebra]] +deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"] +uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" + +[[deps.Literate]] +deps = ["Base64", "IOCapture", "JSON", "REPL"] +git-tree-sha1 = "bb26d8b8ed0fa451ce3511e99c950653a2f31fe1" +uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +version = "2.21.0" + +[[deps.LogExpFunctions]] +deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"] +git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f" +uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688" +version = "0.3.29" + + [deps.LogExpFunctions.extensions] + LogExpFunctionsChainRulesCoreExt = "ChainRulesCore" + LogExpFunctionsChangesOfVariablesExt = "ChangesOfVariables" + LogExpFunctionsInverseFunctionsExt = "InverseFunctions" + + [deps.LogExpFunctions.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ChangesOfVariables = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + +[[deps.Logging]] +uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" + +[[deps.MacroTools]] +git-tree-sha1 = "1e0228a030642014fe5cfe68c2c0a818f9e3f522" +uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" +version = "0.5.16" + +[[deps.Markdown]] +deps = ["Base64"] +uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" + +[[deps.Mmap]] +uuid = "a63ad114-7e13-5084-954f-fe012c677804" + +[[deps.NLSolversBase]] +deps = ["ADTypes", "DifferentiationInterface", "FiniteDiff", "LinearAlgebra"] +git-tree-sha1 = "b3f76b463c7998473062992b246045e6961a074e" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "8.0.0" + +[[deps.NaNMath]] +deps = ["OpenLibm_jll"] +git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae" +uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +version = "1.1.3" + +[[deps.OffsetArrays]] +git-tree-sha1 = "117432e406b5c023f665fa73dc26e79ec3630151" +uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +version = "1.17.0" +weakdeps = ["Adapt"] + + [deps.OffsetArrays.extensions] + OffsetArraysAdaptExt = "Adapt" + +[[deps.OpenBLAS_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] +uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" +version = "0.3.23+4" + +[[deps.OpenLibm_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "05823500-19ac-5b8b-9628-191a04bc5112" +version = "0.8.1+2" + +[[deps.OpenSpecFun_jll]] +deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"] +git-tree-sha1 = "1346c9208249809840c91b26703912dff463d335" +uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" +version = "0.5.6+0" + +[[deps.Optim]] +deps = ["ADTypes", "EnumX", "FillArrays", "LineSearches", "LinearAlgebra", "NLSolversBase", "NaNMath", "PositiveFactorizations", "Printf", "SparseArrays", "Statistics"] +git-tree-sha1 = "7ecee58b8bd88cc8bbdb90d056d1c6546322eebd" +uuid = "429524aa-4258-5aef-a3af-852621145aeb" +version = "2.1.0" + + [deps.Optim.extensions] + OptimMOIExt = "MathOptInterface" + + [deps.Optim.weakdeps] + MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" + +[[deps.OrderedCollections]] +git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee" +uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +version = "1.8.1" + +[[deps.Parsers]] +deps = ["Dates", "PrecompileTools", "UUIDs"] +git-tree-sha1 = "5d5e0a78e971354b1c7bff0655d11fdc1b0e12c8" +uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" +version = "2.8.4" + +[[deps.PartialWaveFunctions]] +git-tree-sha1 = "9e94fb961f02c4eaeb8e83268cb34245395030dd" +uuid = "793d2195-304b-438e-bbb1-bc33c872ac39" +version = "0.2.0" + +[[deps.PositiveFactorizations]] +deps = ["LinearAlgebra"] +git-tree-sha1 = "17275485f373e6673f7e7f97051f703ed5b15b20" +uuid = "85a6dd25-e78a-55b7-8502-1745935b8125" +version = "0.2.4" + +[[deps.PrecompileTools]] +deps = ["Preferences"] +git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f" +uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a" +version = "1.2.1" + +[[deps.Preferences]] +deps = ["TOML"] +git-tree-sha1 = "8b770b60760d4451834fe79dd483e318eee709c4" +uuid = "21216c6a-2e73-6563-6e65-726566657250" +version = "1.5.2" + +[[deps.Primes]] +deps = ["IntegerMathUtils"] +git-tree-sha1 = "25cdd1d20cd005b52fc12cb6be3f75faaf59bb9b" +uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" +version = "0.5.7" + +[[deps.Printf]] +deps = ["Unicode"] +uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" + +[[deps.QuadGK]] +deps = ["DataStructures", "LinearAlgebra"] +git-tree-sha1 = "5e8e8b0ab68215d7a2b14b9921a946fee794749e" +uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +version = "2.11.3" + + [deps.QuadGK.extensions] + QuadGKEnzymeExt = "Enzyme" + + [deps.QuadGK.weakdeps] + Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" + +[[deps.REPL]] +deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"] +uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" + +[[deps.Random]] +deps = ["SHA"] +uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" + +[[deps.RationalRoots]] +git-tree-sha1 = "e5f5db699187a4810fda9181b34250deeedafd81" +uuid = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681" +version = "0.2.1" + +[[deps.Ratios]] +deps = ["Requires"] +git-tree-sha1 = "1342a47bf3260ee108163042310d26f2be5ec90b" +uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439" +version = "0.4.5" + + [deps.Ratios.extensions] + RatiosFixedPointNumbersExt = "FixedPointNumbers" + + [deps.Ratios.weakdeps] + FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93" + +[[deps.Reexport]] +git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "1.2.2" + +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "62389eeff14780bfe55195b7204c0d8738436d64" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.1" + +[[deps.Roots]] +deps = ["Accessors", "CommonSolve", "Printf"] +git-tree-sha1 = "91cfb1cb4f6e27557cc2df798a31eff6089a41eb" +uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +version = "3.0.0" + + [deps.Roots.extensions] + RootsChainRulesCoreExt = "ChainRulesCore" + RootsForwardDiffExt = "ForwardDiff" + RootsIntervalRootFindingExt = "IntervalRootFinding" + RootsSymPyExt = "SymPy" + RootsSymPyPythonCallExt = "SymPyPythonCall" + RootsUnitfulExt = "Unitful" + + [deps.Roots.weakdeps] + ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" + IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807" + SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6" + SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c" + Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" + +[[deps.SHA]] +uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" +version = "0.7.0" + +[[deps.Serialization]] +uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" + +[[deps.Setfield]] +deps = ["ConstructionBase", "Future", "MacroTools", "StaticArraysCore"] +git-tree-sha1 = "c5391c6ace3bc430ca630251d02ea9687169ca68" +uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46" +version = "1.1.2" + +[[deps.SharedArrays]] +deps = ["Distributed", "Mmap", "Random", "Serialization"] +uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383" + +[[deps.Sockets]] +uuid = "6462fe0b-24de-5631-8697-dd941f90decc" + +[[deps.SparseArrays]] +deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"] +uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +version = "1.10.0" + +[[deps.SpecialFunctions]] +deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"] +git-tree-sha1 = "2700b235561b0335d5bef7097a111dc513b8655e" +uuid = "276daf66-3868-5448-9aa4-cd146d93841b" +version = "2.7.2" +weakdeps = ["ChainRulesCore"] + + [deps.SpecialFunctions.extensions] + SpecialFunctionsChainRulesCoreExt = "ChainRulesCore" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] +git-tree-sha1 = "246a8bb2e6667f832eea063c3a56aef96429a3db" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.9.18" +weakdeps = ["ChainRulesCore", "Statistics"] + + [deps.StaticArrays.extensions] + StaticArraysChainRulesCoreExt = "ChainRulesCore" + StaticArraysStatisticsExt = "Statistics" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6ab403037779dae8c514bad259f32a447262455a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.4" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +version = "1.10.0" + +[[deps.StructUtils]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "82bee338d650aa515f31866c460cb7e3bcef90b8" +uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42" +version = "2.8.2" + + [deps.StructUtils.extensions] + StructUtilsMeasurementsExt = ["Measurements"] + StructUtilsStaticArraysCoreExt = ["StaticArraysCore"] + StructUtilsTablesExt = ["Tables"] + + [deps.StructUtils.weakdeps] + Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7" + StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" + Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" + +[[deps.SuiteSparse_jll]] +deps = ["Artifacts", "Libdl", "libblastrampoline_jll"] +uuid = "bea87d4a-7f5b-5778-9afe-8cc45184846c" +version = "7.2.1+1" + +[[deps.TOML]] +deps = ["Dates"] +uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76" +version = "1.0.3" + +[[deps.Test]] +deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] +uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[[deps.UUIDs]] +deps = ["Random", "SHA"] +uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" + +[[deps.Unicode]] +uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" + +[[deps.WignerSymbols]] +deps = ["HalfIntegers", "LRUCache", "Primes", "RationalRoots"] +git-tree-sha1 = "960e5f708871c1d9a28a7f1dbcaf4e0ee34ee960" +uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" +version = "2.0.0" + +[[deps.WoodburyMatrices]] +deps = ["LinearAlgebra", "SparseArrays"] +git-tree-sha1 = "248a7031b3da79a127f14e5dc5f417e26f9f6db7" +uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6" +version = "1.1.0" + +[[deps.libblastrampoline_jll]] +deps = ["Artifacts", "Libdl"] +uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" +version = "5.8.0+1" diff --git a/Project.toml b/Project.toml index 9d53f5a..a693745 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ version = "0.5.0" [deps] Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Optim = "429524aa-4258-5aef-a3af-852621145aeb" PartialWaveFunctions = "793d2195-304b-438e-bbb1-bc33c872ac39" @@ -16,14 +15,12 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" [compat] Interpolations = "0.16" -Literate = "2.20" OffsetArrays = "1.17" Optim = "2.0.1" PartialWaveFunctions = "0.2" diff --git a/src/FewBodyToolkit.jl b/src/FewBodyToolkit.jl index d800650..21df75a 100644 --- a/src/FewBodyToolkit.jl +++ b/src/FewBodyToolkit.jl @@ -6,7 +6,7 @@ using LinearAlgebra, SpecialFunctions, QuadGK, Optim, Roots, StaticArrays, Inter include("common/potentialtypes.jl") include("common/eigen2step.jl") include("common/auxiliary.jl") -export PotentialFunction, CentralPotential, GaussianPotential, ContactPotential1D, comparison#, SpinOrbitPotential +export PotentialFunction, CentralPotential, GaussianPotential, PowerLawPotential, ContactPotential1D, comparison#, SpinOrbitPotential ### GEM-2body include("GEM-2body/GEM2B.jl") diff --git a/src/GEM-2body/GEM2B.jl b/src/GEM-2body/GEM2B.jl index 6da6523..cee97dc 100644 --- a/src/GEM-2body/GEM2B.jl +++ b/src/GEM-2body/GEM2B.jl @@ -116,6 +116,14 @@ function GEM2B_solve!(prealloc_arrs,phys_params,num_params,return_wavefunctions: ## 1. Preliminaries: + # power-law interactions: the radial integral int |r|^(2*lmax+p) exp(-..) |r|^(dim-1) dr + # converges only for p > -2*lmax-dim (e.g. p > -1 for lmax=0 in 1D). + for vint in interactions + if vint isa PowerLawPotential && vint.p <= -2*lmax-dim + error("PowerLawPotential: p = $(vint.p) is too singular for lmax=$lmax in dim=$dim; requires p > $(-2*lmax-dim)") + end + end + # gamma function: Dict due to half-integer arguments. Change to array possible by multiplication of argument by 2 gamma_dict = Dict{Float64, Float64}() for i = 0.5:0.5:max(nmax,2*lmax+1)+0.5 diff --git a/src/GEM-2body/MatrixElements123D.jl b/src/GEM-2body/MatrixElements123D.jl index 218c7ea..1e4dc80 100644 --- a/src/GEM-2body/MatrixElements123D.jl +++ b/src/GEM-2body/MatrixElements123D.jl @@ -29,6 +29,17 @@ function element_V(vint::GaussianPotential,lmax,nu1,nu2,gamma_dict,buf,dim,domai return v0*(2*(nu1*nu2)^(1/2)/(nu1+nu2+mu_g))^(lmax+dim/2) end +# power-law interaction V(r) = v0*|r|^p; the radial integration is closed-form: +# V_{l,n,n'} = v0 * Gamma(L+p/2)/Gamma(L) * element_S(l,nu1,nu2,dim) * (nu1+nu2)^(-p/2), L = lmax+dim/2 +# Complex scaling needs no extra factor: nu1,nu2 arrive already scaled from MatrixV, and +# (nu1+nu2)^(-p/2) then reproduces exactly the rotated potential. +function element_V(vint::PowerLawPotential,lmax,nu1,nu2,gamma_dict,buf,dim,domain,dimfac) #power-law interaction + v0 = vint.v0 + p = vint.p + L = lmax + dim/2 + return v0 * gamma(L+p/2)/gamma_dict[L] * element_S(lmax,nu1,nu2,dim) / (nu1+nu2)^(p/2) +end + function element_V(vint::ContactPotential1D,lmax,nu1,nu2,gamma_dict,buf,dim,domain,dimfac) #contact interaction v0 = vint.v0 z0 = vint.z0 diff --git a/src/GEM-2body/auxiliary.jl b/src/GEM-2body/auxiliary.jl index 5b17954..ce33b66 100644 --- a/src/GEM-2body/auxiliary.jl +++ b/src/GEM-2body/auxiliary.jl @@ -25,6 +25,10 @@ make_phys_params2B(mur=0.5, interactions=[r -> -1/r], lmax=2) # 3D Coulomb ``` """ function make_phys_params2B(;hbar = 1.0, mur=1.0, interactions=[GaussianPotential(-1.0, 1.0)], lmin=0, lmax=0, dim=3) + # Abstract eltype on purpose: a concrete Vector{typeof(v)} makes the whole solver + # pipeline re-specialize for every new potential function. element_V provides the + # function barrier that keeps the quadrature itself fully specialized. + interactions = Vector{Any}(interactions) return (;hbar, mur, interactions, lmin, lmax, dim) end diff --git a/src/GEM-2body/optimV0.jl b/src/GEM-2body/optimV0.jl index d29ac4b..0700bde 100644 --- a/src/GEM-2body/optimV0.jl +++ b/src/GEM-2body/optimV0.jl @@ -45,7 +45,7 @@ function v0GEMOptim(phys_params, num_params, stateindex, target_e2; complex_rang function vint_updt(r) return v0crit*vint(r) end - phys_params = (;hbar,mur,interactions=[vint_updt],lmax,lmin,dim) + phys_params = (;hbar,mur,interactions=Any[vint_updt],lmax,lmin,dim) params[1:3] = GEM_Optim_2B(phys_params, num_params, stateindex; complex_ranged=complex_ranged, g_tol=g_tol) num_params = (;gem_params = (nmax, r1 = params[1], rnmax = params[2]), complex_range_freq, complex_scaling_angle, threshold) @@ -109,7 +109,7 @@ function find_v0crit(phys_params, num_params, stateindex, target_e2; complex_ran return v0*vint(r) end - return real(GEM2B.GEM2B_solve((hbar,mur,interactions=[vint_local],lmax,min,dim), num_params; complex_ranged)[stateindex]) - target_e2; + return real(GEM2B.GEM2B_solve((hbar,mur,interactions=Any[vint_local],lmax,lmin,dim), num_params; complex_ranged)[stateindex]) - target_e2; end v0crit = find_zeros(fun1,(0, 200))[1]; diff --git a/src/GEM-3body-1D/auxiliary.jl b/src/GEM-3body-1D/auxiliary.jl index 6cd8769..34a774f 100644 --- a/src/GEM-3body-1D/auxiliary.jl +++ b/src/GEM-3body-1D/auxiliary.jl @@ -31,6 +31,8 @@ function normalize_species(species) end function make_phys_params3B1D(;hbar = 1.0, masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], parity=+1) + # Abstract eltype on purpose: see make_phys_params2B. + interactions = Vector{Any}[Vector{Any}(vv) for vv in interactions] species = normalize_species(species) return (;hbar, masses, species, interactions, parity) end diff --git a/src/ISGL-3body/ISGL.jl b/src/ISGL-3body/ISGL.jl index 5b82745..a57be17 100644 --- a/src/ISGL-3body/ISGL.jl +++ b/src/ISGL-3body/ISGL.jl @@ -34,7 +34,7 @@ const DEFAULT_OBS = ( ) """ - ISGL_solve(phys_params, num_params; return_wavefunctions=false, complex_scaling=false, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0])) + ISGL_solve(phys_params, num_params; return_wavefunctions=false, complex_scaling=false, complex_ranged=:none, observ_params=(;stateindices=[],centobs_arr=[[],[],[]],R2_arr=[0,0,0])) Solves the 3D three-body problem using the Gaussian Expansion Method (GEM). @@ -43,6 +43,15 @@ Solves the 3D three-body problem using the Gaussian Expansion Method (GEM). - `num_params`: Numerical parameters for the GEM calculation (e.g., basis size, grid parameters, etc.). - `return_wavefunctions`: (optional) If `true`, also returns wavefunction-related observables. Default is `false`. - `complex_scaling`: (optional) If `true`, uses complex scaling method. Default is `false`. +- `complex_ranged`: (optional) Selects complex-ranged basis functions `nu -> nu*(1 + i*omega)` per Jacobi + coordinate, with `omega = complex_range_freq` from `num_params`. One of `:none` (default), `:r` + (only the internal coordinate `` r ``), `:R` (only `` R ``), or `:both`. A `Bool` is also accepted, + with `true` meaning `:both`, for consistency with `GEM2B_solve`. + Each complex-ranged coordinate doubles its number of basis functions, since the ranges enter as + conjugate pairs; choosing `:both` therefore multiplies the total basis size by four. + Only analytically treated potentials (`GaussianPotential`, `PowerLawPotential`) are supported with + complex ranges, as the interpolated path for a generic central potential requires real ranges. + Observables are likewise unsupported with complex ranges. - `observ_params`: (optional) Parameters for observable calculations. + `stateindices`: Indices of states for which observables are calculated. + `centobs_arr`: Array of central (only dependent on `` r ``; must be defined as functions) observables, for each Jacobi set (similar to `interactions` in `phys_params`). @@ -64,9 +73,12 @@ energies = ISGL_solve(phys_params, num_params) #solving with default parameters: ``` """ function ISGL_solve(phys_params, num_params; - return_wavefunctions = false, complex_scaling = false, observ_params=DEFAULT_OBS, debug = false, + return_wavefunctions = false, complex_scaling = false, complex_ranged = :none, observ_params=DEFAULT_OBS, debug = false, wf_bool=nothing, csm_bool=nothing, debug_bool=nothing) + complex_ranged_r, complex_ranged_R = parse_complex_ranged(complex_ranged) + cr_any = complex_ranged_r || complex_ranged_R + if !isnothing(wf_bool) @warn "wf_bool is deprecated, use return_wavefunctions instead" return_wavefunctions = wf_bool @@ -90,26 +102,26 @@ function ISGL_solve(phys_params, num_params; end ## 2. sanity checks: - error_code = sanity_checks(phys_params); + error_code = sanity_checks(phys_params,cr_any,observ_params); if error_code != 0 println("Erroneous inputs. Program stopped") return end - - ## 3. computations to determine sizes of arrays for allocation: - size_params = size_estimate(phys_params,num_params,observ_params,complex_scaling) - - ## 4. preallocation: #is it really necessary? and/or can it not simply be done within the precomputation? better like this for performance analysis - precomp_arrs,temp_arrs,interpol_arrs,fill_arrs,result_arrs = preallocate_data(phys_params,num_params,observ_params,size_params,complex_scaling) + + ## 3. computations to determine sizes of arrays for allocation: + size_params = size_estimate(phys_params,num_params,observ_params,complex_scaling,complex_ranged_r,complex_ranged_R) + + ## 4. preallocation: #is it really necessary? and/or can it not simply be done within the precomputation? better like this for performance analysis + precomp_arrs,temp_arrs,interpol_arrs,fill_arrs,result_arrs = preallocate_data(phys_params,num_params,observ_params,size_params,complex_scaling,complex_ranged_r,complex_ranged_R) ## 5. precomputation: - precompute_ISGL(phys_params,num_params,size_params,precomp_arrs,temp_arrs) - + precompute_ISGL(phys_params,num_params,size_params,precomp_arrs,temp_arrs,complex_ranged_r,complex_ranged_R) + ## 6. preparation of interpolation & shoulder: - interpolNshoulder(phys_params,num_params,observ_params,size_params,precomp_arrs,interpol_arrs,return_wavefunctions,complex_scaling) - + interpolNshoulder(phys_params,num_params,observ_params,size_params,precomp_arrs,interpol_arrs,return_wavefunctions,complex_scaling,cr_any) + ## 7. Calculation of matrix elements - fill_TVS(num_params,size_params,precomp_arrs,interpol_arrs,fill_arrs,complex_scaling,phys_params.hbar,debug) + fill_TVS(num_params,size_params,precomp_arrs,interpol_arrs,fill_arrs,complex_scaling,phys_params.hbar,debug,complex_ranged_r,complex_ranged_R) ## 8. Solving the generalized eigenproblem: solveHS(num_params,fill_arrs,result_arrs,return_wavefunctions) diff --git a/src/ISGL-3body/auxiliary.jl b/src/ISGL-3body/auxiliary.jl index f5bd2e5..c571717 100644 --- a/src/ISGL-3body/auxiliary.jl +++ b/src/ISGL-3body/auxiliary.jl @@ -33,6 +33,8 @@ function normalize_species(species) end function make_phys_params3B3D(;hbar = 1.0, masses=[1.0, 1.0, 1.0], species=[:x,:y,:z], interactions=[[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)],[GaussianPotential(-1.0, 1.0)]], J_tot=0, parity=1, spins=[0,0,0]) + # Abstract eltype on purpose: see make_phys_params2B. + interactions = Vector{Any}[Vector{Any}(vv) for vv in interactions] species = normalize_species(species) return (;hbar, masses, species, interactions, J_tot, parity, spins) end @@ -49,7 +51,7 @@ Create and return a named tuple containing the numerical parameters for a three- - `Lmax::Int = 0`: Maximum power `r^L` used in the basis functions of the `` R `` Jacobi coordinate. - `gem_params::NamedTuple = (nmax=5, r1=1.0, rnmax=10.0, Nmax=5, R1=1.0, RNmax=10.0)`: Parameters for the Gaussian Expansion Method (number of basis functions, smallest and largest range parameters for both Jacobi coordinates). - `complex_scaling_angle::Float64 = 0.0`: Complex scaling angle (in degrees) for the Complex Scaling Method. -- `complex_range_freq::Float64 = 0.9`: Parameter controlling the frequency for complex-ranged basis functions. Currently unsupported. +- `complex_range_freq::Float64 = 0.9`: Parameter `` \\omega `` controlling the frequency for complex-ranged basis functions, `` \\nu \\to \\nu(1 + i\\omega) ``. Only used when `ISGL_solve` is called with `complex_ranged` set to something other than `:none`. The same `` \\omega `` is used for both Jacobi coordinates. - `mu0::Float64 = 0.08`: Parameter (prefactor) for the ISGL method. - `c_shoulder::Float64 = 1.6`: Parameter (base) for the ISGL method. - `kmax_interpol::Int = 1000`: Number of numerical integration with effective Gaussian ranges used for interpolation. @@ -68,4 +70,24 @@ function make_num_params3B3D(;lmin=0, Lmin=0, lmax=0, Lmax=0, gem_params=(nmax=1 complex_scaling_angle = isnothing(theta_csm) ? complex_scaling_angle : theta_csm complex_range_freq = isnothing(omega_cr) ? complex_range_freq : omega_cr return (;lmax, Lmax, gem_params, complex_scaling_angle, complex_range_freq, mu0, c_shoulder, kmax_interpol, threshold, lmin, Lmin) -end \ No newline at end of file +end + + +""" + parse_complex_ranged(complex_ranged) -> (complex_ranged_r, complex_ranged_R) + +Translates the user-facing `complex_ranged` option into two booleans, one per Jacobi coordinate. + +Accepted values are the symbols `:none`, `:r`, `:R`, `:both`, and a `Bool` (`true` meaning `:both`, +`false` meaning `:none`) for consistency with `GEM2B_solve`, whose basis has only one coordinate. + +Kept here rather than in `common/` for now; it should move once `GEM3B1D` gains the same option. +""" +function parse_complex_ranged(complex_ranged) + complex_ranged isa Bool && return (complex_ranged, complex_ranged) + complex_ranged === :none && return (false, false) + complex_ranged === :r && return (true, false) + complex_ranged === :R && return (false, true) + complex_ranged === :both && return (true, true) + error("complex_ranged must be one of :none, :r, :R, :both (or a Bool); got $(repr(complex_ranged))") +end diff --git a/src/ISGL-3body/fillTVS.jl b/src/ISGL-3body/fillTVS.jl index 7cb3e96..ebbbce8 100644 --- a/src/ISGL-3body/fillTVS.jl +++ b/src/ISGL-3body/fillTVS.jl @@ -1,33 +1,53 @@ ## Function for calculating the matrix elements and filling the matrices T,V,S within the ISGL program -@views @inbounds function fill_TVS(num_params,size_params,precomp_arrs,interpol_arrs,fill_arrs,complex_scaling::Bool,hbar,debug::Bool) - +@views @inbounds function fill_TVS(num_params,size_params,precomp_arrs,interpol_arrs,fill_arrs,complex_scaling::Bool,hbar,debug::Bool,complex_ranged_r::Bool=false,complex_ranged_R::Bool=false) + (;gem_params,mu0,c_shoulder,complex_scaling_angle) = num_params (;nmax,Nmax,r1,rnmax,R1,RNmax) = gem_params - (;abvals_arr,cvals,gauss_indices,central_indices,so_indices,groupindex_arr,abI,factor_bf,box_size_arr,starts,ends,bvalsdiag,s_arr,JsS_arr,s_complete,JsS_complete,JlL_arr,lL_nested,maxlmax,mij_arr_dict,mijSO_arr_dict,gaussopt_arr) = size_params + (;abvals_arr,cvals,gauss_indices,central_indices,so_indices,pow_indices,groupindex_arr,abI,factor_bf,box_size_arr,starts,ends,bvalsdiag,s_arr,JsS_arr,s_complete,JsS_complete,JlL_arr,lL_nested,maxlmax,mij_arr_dict,mijSO_arr_dict,gaussopt_arr,powopt_arr) = size_params (;gamma_dict,spintrafo_dict,spinoverlap_dict,global6j_dict,facsymm_dict,jmat,murR_arr,nu_arr,NU_arr,norm_arr,NORM_arr,Clmk_arr,Dlmk_arr,S_arr,SSO_arr) = precomp_arrs - (;alpha_arr,v_arr,A_mat,w_interpol_arr,Ainv_arr_kine) = interpol_arrs + (;alpha_arr,v_arr,A_mat,w_interpol_arr,Ainv_arr_kine,w_pow_arr) = interpol_arrs (;w_arr_kine,wn_interpol_arr,kij_arr,gij_arr,T,V,S,temp_args_arr,temp_fill_mat) = fill_arrs + complex_ranged = complex_ranged_r || complex_ranged_R + # With complex ranges the bra ranges are conjugated, so S, T and V become hermitian instead of + # symmetric. Combined with complex scaling neither symmetry survives and the full matrix is filled. + fill_full = complex_ranged && complex_scaling + + # the effective numbers of ranges already include the doubling for a complex-ranged coordinate + nmax_eff = lastindex(nu_arr) + Nmax_eff = lastindex(NU_arr) + # reducing complicated many-loop structure to a single 1D loop: - flati = flattento1Dloop(temp_args_arr,groupindex_arr,bvalsdiag,abvals_arr,s_arr,JsS_arr,JlL_arr,lL_nested,nmax,Nmax,nu_arr,NU_arr,norm_arr,NORM_arr,mij_arr_dict,starts) - + flati = flattento1Dloop(temp_args_arr,groupindex_arr,bvalsdiag,abvals_arr,s_arr,JsS_arr,JlL_arr,lL_nested,nmax_eff,Nmax_eff,nu_arr,NU_arr,norm_arr,NORM_arr,mij_arr_dict,starts,fill_full) + ## Calculation of matrix elements and matrix filling via 1d loop (so we dont have redundant loops over all spin configs, lL combinations, and nmax) for index in 1:flati (;rowi,coli) = temp_args_arr[index] temp_fill_mat[rowi,coli] = sab(jmat,temp_args_arr[index],abI,factor_bf,S_arr,spintrafo_dict,facsymm_dict) end - # transpose fill: - S .= Symmetric(temp_fill_mat,:L); - + # transpose fill: the overlap is hermitian for complex ranges (also together with complex scaling, + # since the ISGL complex-scaling rotates the potential and leaves the ranges alone) + if complex_ranged + hermitian_fill!(S,temp_fill_mat); + else + S .= Symmetric(temp_fill_mat,:L); + end + #combined t and v for some speedup: reverted combination for CSM for index in 1:flati (;rowi,coli) = temp_args_arr[index] temp_fill_mat[rowi,coli] = tab(jmat,murR_arr,w_arr_kine,Ainv_arr_kine,kij_arr,mu0,c_shoulder,temp_args_arr[index],abI,factor_bf,S_arr,spintrafo_dict,facsymm_dict,hbar) end # transpose fill: - T .= Symmetric(temp_fill_mat,:L); - + if fill_full + T .= temp_fill_mat; + elseif complex_ranged + hermitian_fill!(T,temp_fill_mat); + else + T .= Symmetric(temp_fill_mat,:L); + end + if complex_scaling T .*= exp(-2*im*complex_scaling_angle*pi/180) end @@ -35,11 +55,17 @@ #v for index in 1:flati (;rowi,coli) = temp_args_arr[index] - temp_fill_mat[rowi,coli] = vab(jmat,gij_arr,mu0,c_shoulder,w_interpol_arr,wn_interpol_arr,temp_args_arr[index],abI,factor_bf,S_arr,SSO_arr,cvals,spintrafo_dict,spinoverlap_dict,facsymm_dict,gauss_indices,central_indices,so_indices,s_arr,global6j_dict,mijSO_arr_dict,gaussopt_arr,complex_scaling) # do we need hbar^2 for SO? + temp_fill_mat[rowi,coli] = vab(jmat,gij_arr,mu0,c_shoulder,w_interpol_arr,wn_interpol_arr,w_pow_arr,temp_args_arr[index],abI,factor_bf,S_arr,SSO_arr,cvals,spintrafo_dict,spinoverlap_dict,facsymm_dict,gauss_indices,central_indices,so_indices,pow_indices,s_arr,global6j_dict,mijSO_arr_dict,gaussopt_arr,powopt_arr,complex_scaling) # do we need hbar^2 for SO? end # transpose fill: - V .= Symmetric(temp_fill_mat,:L); - + if fill_full + V .= temp_fill_mat; + elseif complex_ranged + hermitian_fill!(V,temp_fill_mat); + else + V .= Symmetric(temp_fill_mat,:L); + end + if debug stp = min(9, size(T, 1)) # Adjust size_to_print as needed println("T:") @@ -55,10 +81,24 @@ end +# Reconstruct a hermitian matrix from its lower triangle. +# Hermitian(src,:L) leaves the stored diagonal untouched, so the O(eps) imaginary parts that the +# matrix elements pick up numerically would survive and make ishermitian(dest) false. eigen2step +# keys its "the eigenvalues are real" decision on exactly that predicate, so the diagonal is made +# real explicitly here. Mathematically the diagonal of S, T and V is real anyway. +function hermitian_fill!(dest,src) + dest .= Hermitian(src,:L) + for i in axes(dest,1) + dest[i,i] = real(dest[i,i]) + end + return dest +end + + ## functions to calculate one matrix element, summing over the necessary a,b,c values: function sab(jmat,temp_args_i,abI,factor_bf,S_arr,spintrafo_dict,facsymm_dict) (;rowi,coli,avals_new,bvals_new,factor_ab,ranges,norm4,mij_arr,sa,JsSa,sb,JsSb,la,La,lb,Lb,JlLa,JlLb) = temp_args_i - tempS = 0.0 + tempS = zero(norm4) # complex for complex-ranged basis functions # maybe easier: (JsSa != JsSb || JlLa != JlLb) && return tempS #immediately skip if either one is violated @@ -78,7 +118,7 @@ end function tab(jmat,murR_arr,w_arr_kine,Ainv_arr_kine,kij_arr,mu0,c_shoulder,temp_args_i,abI,factor_bf,S_arr,spintrafo_dict,facsymm_dict,hbar) (;avals_new,bvals_new,factor_ab,ranges,norm4,mij_arr,sa,JsSa,sb,JsSb,la,La,lb,Lb,Lsum,JlLa,JlLb) = temp_args_i - tempT = 0.0 + tempT = zero(promote_type(typeof(norm4),eltype(w_arr_kine))) # maybe easier: (JsSa != JsSb || JlLa != JlLb) && return tempT #immediately skip if either one is violated @@ -96,9 +136,10 @@ function tab(jmat,murR_arr,w_arr_kine,Ainv_arr_kine,kij_arr,mu0,c_shoulder,temp_ return tempT end -function vab(jmat,gij_arr,mu0,c_shoulder,w_interpol_arr,wn_interpol_arr,temp_args_i,abI,factor_bf,S_arr,SSO_arr,cvals,spintrafo_dict,spinoverlap_dict,facsymm_dict,gauss_indices,central_indices,so_indices,s_arr,global6j_dict,mijSO_arr_dict,gaussopt_arr,complex_scaling) +function vab(jmat,gij_arr,mu0,c_shoulder,w_interpol_arr,wn_interpol_arr,w_pow_arr,temp_args_i,abI,factor_bf,S_arr,SSO_arr,cvals,spintrafo_dict,spinoverlap_dict,facsymm_dict,gauss_indices,central_indices,so_indices,pow_indices,s_arr,global6j_dict,mijSO_arr_dict,gaussopt_arr,powopt_arr,complex_scaling) (;rowi,coli,avals_new,bvals_new,factor_ab,ranges,norm4,mij_arr,sa,JsSa,sb,JsSb,la,La,lb,Lb,Lsum,JlLa,JlLb) = temp_args_i - tempV = 0.0 + # complex for complex ranges (via norm4/gij_arr) and/or complex scaling (via wn_interpol_arr) + tempV = zero(promote_type(typeof(norm4),eltype(gij_arr),eltype(wn_interpol_arr))) for a in avals_new for b in bvals_new @@ -117,6 +158,12 @@ function vab(jmat,gij_arr,mu0,c_shoulder,w_interpol_arr,wn_interpol_arr,temp_arg tempV += factor_ab*factor_symm*uab*element_VGauss(c,ranges,norm4,jmat[a,c],jmat[b,c],mij_arr,S_arr,la,La,lb,Lb,gij_arr,mu0,c_shoulder,Lsum,wn_interpol_arr,v0,mu_g,JlLa) end + for ivp in pow_indices[c] #loop over the power-law interactions for this c. + (JsSa != JsSb || JlLa != JlLb) && return tempV #immediately skip if it is violated. + v0,p_pow = powopt_arr[c][ivp] + + tempV += factor_ab*factor_symm*uab*element_VPow(c,ranges,norm4,jmat[a,c],jmat[b,c],mij_arr,S_arr,la,La,lb,Lb,gij_arr,mu0,c_shoulder,Lsum,w_pow_arr,v0,p_pow,JlLa,ivp) + end for ivc in central_indices[c] #loop over the central interactions for this c. (JsSa != JsSb || JlLa != JlLb) && return tempV #immediately skip if it is violated. tempV += factor_ab*factor_symm*uab*element_V(c,ranges,norm4,jmat[a,c],jmat[b,c],mij_arr,S_arr,la,La,lb,Lb,gij_arr,mu0,c_shoulder,w_interpol_arr,Lsum,wn_interpol_arr,JlLa,ivc) @@ -154,10 +201,12 @@ function element_S(ranges,norm4,jab,mij_arr_i,S_arr,la,La,lb,Lb,JlL) f = xi/zeta q = ppr .- f*p - fij_arr = @SMatrix[2*p[i]*p[j]/zeta + 2*q[i]*q[j]/etapr for i=1:4, j=1:4] - prefac = norm4 * (pi^2/zeta/etapr)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); - - sum = 0.0 + fij_arr = @SMatrix[2*p[i]*p[j]/zeta + 2*q[i]*q[j]/etapr for i=1:4, j=1:4] + # the two 3/2-powers are kept separate on purpose: for complex ranges each factor is then + # continued from its own (small) argument, instead of from the argument of the merged product. + prefac = norm4 * (pi/zeta)^(3/2)*(pi/etapr)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); + + sum = zero(prefac) for ii = 1:lastindex(mij_arr_i) (m12,m13,m14,m23,m24,m34) = mij_arr_i[ii] sum += S_arr[la,La,lb,Lb,JlL,ii]*fij_arr[1,2]^m12*fij_arr[1,3]^m13*fij_arr[1,4]^m14*fij_arr[2,3]^m23*fij_arr[2,4]^m24*fij_arr[3,4]^m34 * prefac @@ -209,13 +258,13 @@ function element_T(ranges,norm4,jab,murR_arr,mij_arr_i,S_arr,la,La,lb,Lb,w_arr_k end end - prefac = norm4 * (pi^2/zeta/etapr)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); - - sum = 0.0 + prefac = norm4 * (pi/zeta)^(3/2)*(pi/etapr)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); + + sum = zero(promote_type(typeof(prefac),eltype(w_arr_kine))) for ii = 1:lastindex(mij_arr_i) (m12,m13,m14,m23,m24,m34) = mij_arr_i[ii] - - sum2 = 0.0 + + sum2 = zero(promote_type(eltype(w_arr_kine),eltype(kij_arr))) for n=0:Lsum sum2 +=w_arr_kine[n]*kij_arr[1,2,n]^m12*kij_arr[1,3,n]^m13*kij_arr[1,4,n]^m14*kij_arr[2,3,n]^m23*kij_arr[2,4,n]^m24*kij_arr[3,4,n]^m34 end @@ -254,18 +303,20 @@ function element_V(c,ranges,norm4,jac,jbc,mij_arr_i,S_arr,la,La,lb,Lb,gij_arr,mu end end - prefac = norm4 * (pi^2/zetac/etaprc)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); - - # here, interpolation is called (outside of ii-loop): + prefac = norm4 * (pi/zetac)^(3/2)*(pi/etaprc)^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); + + # here, interpolation is called (outside of ii-loop). + # NOTE: this path requires a real etaprc, i.e. real Gaussian ranges. Complex-ranged basis + # functions are therefore rejected for interpolated potentials in sanity_checks. for n = 0:Lsum wn_interpol_arr[n] = w_interpol_arr[c,iv,Lsum,n](log(etaprc)) end - - summe = 0.0 + + summe = zero(promote_type(typeof(prefac),eltype(wn_interpol_arr),eltype(gij_arr))) for ii = 1:lastindex(mij_arr_i) # the correct mij_arr is already selected in the flatento1Dloop function (m12,m13,m14,m23,m24,m34) = mij_arr_i[ii] - - sum2 = 0.0 + + sum2 = zero(promote_type(eltype(wn_interpol_arr),eltype(gij_arr))) for n=0:Lsum sum2 += wn_interpol_arr[n]*gij_arr[1,2,n]^m12*gij_arr[1,3,n]^m13*gij_arr[1,4,n]^m14*gij_arr[2,3,n]^m23*gij_arr[2,4,n]^m24*gij_arr[3,4,n]^m34 end @@ -276,6 +327,61 @@ function element_V(c,ranges,norm4,jac,jbc,mij_arr_i,S_arr,la,La,lb,Lb,gij_arr,mu return prefac*summe end + +# calculation of a single matrix element: power-law interaction V(r) = v0*r^p (analytic) +# Same structure as element_V, but the interpolated weights w_interpol_arr[c,iv,Lsum,n](log(etaprc)) +# are replaced by the exact w_n(etaprc) = etaprc^(-p/2) * w_pow_arr[c,iv,Lsum,n]. The common factor +# etaprc^(-p/2) is independent of n and is therefore absorbed into prefac. +function element_VPow(c,ranges,norm4,jac,jbc,mij_arr_i,S_arr,la,La,lb,Lb,gij_arr,mu0,c_shoulder,Lsum,w_pow_arr,v0,p_pow,JlL,iv) + + (;nua,nub,NUa,NUb) = ranges; + (alphaAC,gammaAC,betaAC,deltaAC) = jac + (alphaBC,gammaBC,betaBC,deltaBC) = jbc # careful with order (gamma before beta) + + etac = nua*alphaAC^2 + NUa*gammaAC^2 + nub*alphaBC^2 + NUb*gammaBC^2 + zetac = nua*betaAC^2 + NUa*deltaAC^2 + nub*betaBC^2 + NUb*deltaBC^2 + xic = nua*alphaAC*betaAC + NUa*gammaAC*deltaAC + nub*alphaBC*betaBC + NUb*gammaBC*deltaBC + etaprc = etac - xic^2/zetac + + # p, pprime, f, q: + p = SA[nua*betaAC,NUa*deltaAC,nub*betaBC,NUb*deltaBC] + ppr = SA[nua*alphaAC,NUa*gammaAC,nub*alphaBC,NUb*gammaBC] + f = xic/zetac + #fpr = xic/etac + q = ppr .- f*p + + #gij_arr: + for n = 0:Lsum + mun = mu0*c_shoulder^n + for i=1:3 + for j=(i+1):4 + gij_arr[i,j,n] = 2*p[i]*p[j]/zetac + 2*mun*q[i]*q[j]/etaprc + end + end + end + + # etaprc^(-p/2) uses the principal branch. For complex ranges arg(etaprc) is bounded by + # atan(omega) < pi/2, so the principal branch is the continuous continuation of the real-range + # result; x^y is evaluated as exp(y*log(x)) with the principal log of the *base*, so no wrapping + # occurs even for large |p|. + prefac = norm4 * (pi/zetac)^(3/2)*(pi/etaprc)^(3/2) * etaprc^(-p_pow/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); + + # no interpolation necessary here: w_pow_arr is exact and alpha-independent + summe = zero(promote_type(typeof(prefac),eltype(w_pow_arr),eltype(gij_arr))) + for ii = 1:lastindex(mij_arr_i) # the correct mij_arr is already selected in the flatento1Dloop function + (m12,m13,m14,m23,m24,m34) = mij_arr_i[ii] + + sum2 = zero(promote_type(eltype(w_pow_arr),eltype(gij_arr))) + for n=0:Lsum + sum2 += w_pow_arr[c,iv,Lsum,n]*gij_arr[1,2,n]^m12*gij_arr[1,3,n]^m13*gij_arr[1,4,n]^m14*gij_arr[2,3,n]^m23*gij_arr[2,4,n]^m24*gij_arr[3,4,n]^m34 + end + + summe += S_arr[la,La,lb,Lb,JlL,ii]*sum2 + end + + return v0*prefac*summe +end + # matrix element for Spin-Orbit interaction; postponed to future version #= function element_VSO(c,ranges,norm4,jac,jbc,mijSO_arr_dict,SSO_arr,la,La,lb,Lb,gij_arr,mu0,c_shoulder,w_interpol_arr,Lsum,wn_interpol_arr,JlLa,JlLb,ivso) @@ -379,12 +485,12 @@ function element_VGauss(c,ranges,norm4,jac,jbc,mij_arr_i,S_arr,la,La,lb,Lb,gij_a #gij_arr: ggij_arr = @SMatrix[2*p[i]*p[j]/zetac + 2*q[i]*q[j]/(etaprc + mu_g) for i=1:4, j=1:4] - prefac = norm4 * (pi^2/zetac/(etaprc + mu_g))^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); + prefac = norm4 * (pi/zetac)^(3/2)*(pi/(etaprc + mu_g))^(3/2)/(nua^la*NUa^La*nub^lb*NUb^Lb); - summe = 0.0 + summe = zero(eltype(ggij_arr)) for ii = 1:lastindex(mij_arr_i) (m12,m13,m14,m23,m24,m34) = mij_arr_i[ii] - + summe += S_arr[la,La,lb,Lb,JlL,ii]*ggij_arr[1,2]^m12*ggij_arr[1,3]^m13*ggij_arr[1,4]^m14*ggij_arr[2,3]^m23*ggij_arr[2,4]^m24*ggij_arr[3,4]^m34 end @@ -393,19 +499,22 @@ end # returns flati and fills temp_args_arr -function flattento1Dloop(temp_args_arr,groupindex_arr,bvalsdiag,abvals_arr,s_arr,JsS_arr,JlL_arr,lL_nested,nmax,Nmax,nu_arr,NU_arr,norm_arr,NORM_arr,mij_arr_dict,starts) +function flattento1Dloop(temp_args_arr,groupindex_arr,bvalsdiag,abvals_arr,s_arr,JsS_arr,JlL_arr,lL_nested,nmax,Nmax,nu_arr,NU_arr,norm_arr,NORM_arr,mij_arr_dict,starts,fill_full::Bool=false) # Keep loop-structure and write necessary functions arguments for matrix-element-calculation into 1-dim array temp_args_arr + # nmax,Nmax are the *effective* numbers of ranges here, i.e. already doubled for a complex-ranged coordinate. flati = 0 # Iterate over boxes: for boxC in groupindex_arr for boxR in groupindex_arr - boxR < boxC && continue # fill only lower-triangular (boxes, not elements!) only works for real-symmetric or hermitian matrices; NOT anymore for CSM? Also works for CSM -> complex symmetric matrices - + # fill only lower-triangular (boxes, not elements!). Valid for real-symmetric, complex-symmetric + # (CSM) and hermitian (complex ranges) matrices, but not for CSM and complex ranges together. + !fill_full && boxR < boxC && continue + # if there are some identical particles: we can ignore the sum over a-values and simply multiply by a factor which is equal to the number of a-values. ONLY ON THE BOX-DIAGONAL! (boxC = boxR) if boxC == boxR bvals_new = bvalsdiag[boxC] # for changing the role of a,b to be in line with lower triangular! factor_ab = lastindex(abvals_arr[boxR]) # factor for amount of a-values normally - diag_bool = 1 # for skipping lower-triangular calculation within each box! + diag_bool = fill_full ? 0 : 1 # for skipping lower-triangular calculation within each box! else #avals_new = abvals_arr[boxR] bvals_new = abvals_arr[boxC] @@ -440,11 +549,13 @@ function flattento1Dloop(temp_args_arr,groupindex_arr,bvalsdiag,abvals_arr,s_arr Lsum=Int64((la+La+lb+Lb)/2) mij_arr = mij_arr_dict[(la,La),(lb,Lb)] for na = 1:nmax - nua = nu_arr[na] - norma = norm_arr[la,na] + # the bra ranges enter conjugated. This is a no-op for real + # ranges and is what makes S,T,V hermitian for complex ones. + nua = conj(nu_arr[na]) + norma = conj(norm_arr[la,na]) for Na = 1:Nmax - NUa = NU_arr[Na] - NORMa = NORM_arr[La,Na] + NUa = conj(NU_arr[Na]) + NORMa = conj(NORM_arr[La,Na]) alpha += 1 diag_bool == 1 && alpha < alphab && continue # skip upper triangular only on diagonal boxes diff --git a/src/ISGL-3body/interpolationNshoulder.jl b/src/ISGL-3body/interpolationNshoulder.jl index 676bae4..ddf2ce0 100644 --- a/src/ISGL-3body/interpolationNshoulder.jl +++ b/src/ISGL-3body/interpolationNshoulder.jl @@ -1,19 +1,27 @@ # functions to precompute the w-arrays for the interaction via interpolation and the upon-the shoulder method -function interpolNshoulder(phys_params,num_params,observ_params,size_params,precomp_arrs,interpol_arrs,return_wavefunctions::Bool,complex_scaling::Bool) - +function interpolNshoulder(phys_params,num_params,observ_params,size_params,precomp_arrs,interpol_arrs,return_wavefunctions::Bool,complex_scaling::Bool,complex_ranged::Bool=false) + # Destruct Structs: (;interactions) = phys_params (;lmax,Lmax,gem_params,complex_scaling_angle,complex_range_freq,mu0,c_shoulder,kmax_interpol,complex_scaling_angle) = num_params (;nmax,Nmax,r1,rnmax,R1,RNmax) = gem_params (;stateindices,centobs_arr,R2_arr) = observ_params - (;cvals,central_indices,so_indices,maxlmax,nint_arr) = size_params + (;cvals,central_indices,so_indices,pow_indices,powopt_arr,maxlmax,nint_arr) = size_params (;gamma_dict,jmat,nu_arr,NU_arr) = precomp_arrs - (;alpha_arr,v_arr,A_mat,w_arr,w_interpol_arr,Ainv_arr_kine,v_obs_arr,w_obs_arr,w_obs_interpol_arr) = interpol_arrs + (;alpha_arr,v_arr,A_mat,w_arr,w_interpol_arr,Ainv_arr_kine,v_pow,w_pow_arr,v_obs_arr,w_obs_arr,w_obs_interpol_arr) = interpol_arrs # range interpolation: - precompute_alpha_arr(alpha_arr,r1,rnmax,R1,RNmax,nu_arr,NU_arr,jmat) - + # The interpolation mesh is only meaningful for real ranges: it is looked up at log(etaprc), + # which is complex once the ranges are. Complex ranges are therefore restricted to the + # analytically treated potentials (see sanity_checks), and the mesh is not used at all; a + # placeholder mesh is filled so that the (then empty) interpolation loops below stay well-defined. + if complex_ranged + buildnu(10.0,1.0,lastindex(alpha_arr),alpha_arr) # arguments swapped so that alpha_arr is increasing, as in precompute_alpha_arr + else + precompute_alpha_arr(alpha_arr,r1,rnmax,R1,RNmax,nu_arr,NU_arr,jmat) + end + # wrap function types: # Auto-wrap plain functions as central, to ensure compatibility wrap_potential(f::Function) = CentralPotential(f) @@ -24,7 +32,7 @@ function interpolNshoulder(phys_params,num_params,observ_params,size_params,prec # upon-the-shoulder - precompute_w(w_arr,v_arr,alpha_arr,A_mat,w_interpol_arr,Ainv_arr_kine,gamma_dict,maxlmax,mu0,c_shoulder,cvals,vint_arr_wrapped,centobs_arr_wrapped,w_obs_arr,v_obs_arr,w_obs_interpol_arr,return_wavefunctions,complex_scaling,complex_scaling_angle,central_indices,so_indices,nint_arr) + precompute_w(w_arr,v_arr,alpha_arr,A_mat,w_interpol_arr,Ainv_arr_kine,v_pow,w_pow_arr,gamma_dict,maxlmax,mu0,c_shoulder,cvals,vint_arr_wrapped,centobs_arr_wrapped,w_obs_arr,v_obs_arr,w_obs_interpol_arr,return_wavefunctions,complex_scaling,complex_scaling_angle,central_indices,so_indices,pow_indices,powopt_arr,nint_arr) end @@ -50,7 +58,7 @@ function precompute_alpha_arr(alpha_arr,r1,rnmax,R1,RNmax,nu_arr,NU_arr,jmat) Aa = SA[nua 0.0 ; 0.0 NUa] Ab = SA[nub 0.0 ; 0.0 NUb] tempA = transpose(jmat[a,c])*Aa*jmat[a,c] + transpose(jmat[b,c])*Ab*jmat[b,c] - temp = det(tempA)/tempA[2,2] + temp = abs(det(tempA)/tempA[2,2]) # abs is a no-op for real ranges (etaprc>0); it keeps this scan well-defined should complex ranges ever be routed through the interpolation path i==1 && (tempmin = temp; tempmax=temp;) temp > tempmax && (tempmax = temp) temp < tempmin && (tempmin = temp) @@ -95,6 +103,20 @@ function precompute_varr!(v_arr,alpha_arr,Lsum,gamma_dict,vcent_fun::SpinOrbitPo end end =# +# power-law interaction: the integral is analytic, so no numerical integration is needed. +# The main interaction path does not use this (see the pow_indices block in precompute_w), but +# observables (centobs_arr) go through precompute_varr! as well, so a method is provided here. +function precompute_varr!(v_arr,alpha_arr,Lsum,gamma_dict,vcent_fun::PowerLawPotential,buf,csmfac) + (;v0,p) = vcent_fun + for n = 0:Lsum + for k=1:lastindex(alpha_arr) + norm_interpol = 1/2 * gamma_dict[n+1.5]/alpha_arr[k]^(n+3/2) + vcent_analytic = v0*csmfac^(-p) * 1/2*gamma(n+(3+p)/2)/alpha_arr[k]^(n+(3+p)/2) + v_arr[k,n+1] = vcent_analytic/gamma_dict[n+1.0]/norm_interpol + end + end +end + function vcent_integration(vcent_fun,alpha,n,buf) #where {V} val = quadgk(r -> integrand(r,alpha,n,vcent_fun),0,Inf;segbuf=buf)[1] end @@ -104,7 +126,7 @@ end ### w_arr: upon-the-shoulder method -@views @inbounds function precompute_w(w_arr,v_arr,alpha_arr,A_mat,w_interpol_arr,Ainv_arr_kine,gamma_dict,maxlmax,mu0,c_shoulder,cvals,interactions,centobs_arr,w_obs_arr,v_obs_arr,w_obs_interpol_arr,return_wavefunctions::Bool,complex_scaling::Bool,complex_scaling_angle,central_indices,so_indices,nint_arr) +@views @inbounds function precompute_w(w_arr,v_arr,alpha_arr,A_mat,w_interpol_arr,Ainv_arr_kine,v_pow,w_pow_arr,gamma_dict,maxlmax,mu0,c_shoulder,cvals,interactions,centobs_arr,w_obs_arr,v_obs_arr,w_obs_interpol_arr,return_wavefunctions::Bool,complex_scaling::Bool,complex_scaling_angle,central_indices,so_indices,pow_indices,powopt_arr,nint_arr) # returns the Array w_arr[c in cvals,alpha=1:alphamax,Lsum = 1:2*maxlmax,n=1:Lsum+1] # note the +1 in the last argument: w_arr, v_arr are NOT offset-arrays due to problems with linear algebra package. @@ -142,6 +164,28 @@ end end end + # necessary for power-law interactions: + # For V(r)=v0*r^p the radial integral of the range-interpolation method is closed-form, + # Vtilde_n(alpha) = Gamma(n+(3+p)/2)/Gamma(n+3/2) * alpha^(-p/2), i.e. the alpha-dependence is a + # single power alpha^(-p/2), the same for every n. It therefore factors out of the shoulder solve + # and is applied later in element_VPow (via etaprc^(-p/2)); no interpolation over alpha is needed. + # v0 is applied in element_VPow as well (as for the Gaussian), so only p enters here. + for cc in cvals + for iv in pow_indices[cc] + p_pow = powopt_arr[cc][iv][2] + for j = 0:Lsum + v_pow[j+1] = csmfac^(-p_pow) * gamma(j+(3+p_pow)/2)/gamma_dict[j+1.5]/gamma_dict[j+1.0] + end + for n = 0:Lsum # w = Ainv*v_pow, using the Ainv already formed above (as for Ainv_arr_kine) + temp_pow = zero(eltype(v_pow)) + for j = 0:Lsum + temp_pow += Ainv[n+1,j+1]*v_pow[j+1] + end + w_pow_arr[cc,iv,Lsum,n] = temp_pow + end + end + end + for cc in cvals #performance: precompute_varr needs more time than the interpolation procedure for w_arr below. this is mostly due to the use of quadgk diff --git a/src/ISGL-3body/preallocate.jl b/src/ISGL-3body/preallocate.jl index 16e321b..cbf37b3 100644 --- a/src/ISGL-3body/preallocate.jl +++ b/src/ISGL-3body/preallocate.jl @@ -1,6 +1,9 @@ # functions to preallocate all arrays for the ISGL program -struct PrecomputeStruct +# TR is the element type of the Gaussian ranges, and of everything derived from the ranges alone +# (norms, gij/kij, the overlap matrix S). It is Float64 for the usual real ranges and ComplexF64 +# as soon as complex-ranged basis functions are used in at least one Jacobi coordinate. +struct PrecomputeStruct{TR} gamma_dict::Dict{Float64, Float64} cleb_arr::OffsetArray{Float64, 5, Array{Float64, 5}} spintrafo_dict::Dict{Tuple{Int64,Int64,Float64,Float64,Float64},Float64} @@ -9,10 +12,10 @@ struct PrecomputeStruct facsymm_dict::Dict{Tuple{Int64,Int64,Int64,Int64,Float64,Float64},Float64} jmat::Matrix{SMatrix{2, 2, Float64, 4}} murR_arr::MMatrix{2, 3, Float64, 6} - nu_arr::Vector{Float64} - NU_arr::Vector{Float64} - norm_arr::OffsetMatrix{Float64, Matrix{Float64}}#OffsetMatrix{Float64, Matrix{Float64}} - NORM_arr::OffsetMatrix{Float64, Matrix{Float64}}#OffsetMatrix{Float64, Matrix{Float64}} + nu_arr::Vector{TR} + NU_arr::Vector{TR} + norm_arr::OffsetMatrix{TR, Matrix{TR}} + NORM_arr::OffsetMatrix{TR, Matrix{TR}} Clmk_arr::OffsetArray{Float64, 3, Array{Float64, 3}} # 3,4,5 defines the dimensions of the array, not the size. Dlmk_arr::OffsetArray{ComplexF64, 4, Array{ComplexF64, 4}} S_arr::OffsetArray{Float64, 6, Array{Float64, 6}} @@ -34,18 +37,20 @@ struct InterpolationStruct{T} w_arr::Array{T, 4} w_interpol_arr::OffsetArray{Interpolations.Extrapolation{T, 1, ScaledInterpolation{T, 1, Interpolations.BSplineInterpolation{T, 1, OffsetVector{T, Vector{T}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Throw{Nothing}}, 4, Array{Interpolations.Extrapolation{T, 1, ScaledInterpolation{T, 1, Interpolations.BSplineInterpolation{T, 1, OffsetVector{T, Vector{T}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Throw{Nothing}}, 4}} Ainv_arr_kine::OffsetMatrix{Float64, Matrix{Float64}} + v_pow::Vector{T} + w_pow_arr::OffsetArray{T, 4, Array{T, 4}} v_obs_arr::Matrix{Float64} w_obs_arr::Array{Float64, 5} w_obs_interpol_arr::OffsetArray{Interpolations.Extrapolation{Float64, 1, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetVector{Float64, Vector{Float64}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Throw{Nothing}}, 4, Array{Interpolations.Extrapolation{Float64, 1, ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetVector{Float64, Vector{Float64}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{Base.OneTo{Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSpline{Cubic{Line{OnGrid}}}, Throw{Nothing}}, 4}} end # new struct for temporary arguments -struct TempArgs +struct TempArgs{TR} rowi::Int coli::Int ranges::NamedTuple{(:nua, :nub, :NUa, :NUb), - Tuple{Float64,Float64,Float64,Float64}} - norm4::Float64 + NTuple{4,TR}} + norm4::TR mij_arr::Array{SArray{Tuple{6},Int,1,6},1} sa::Float64 JsSa::Float64 @@ -65,41 +70,62 @@ struct TempArgs bvals::Vector{Int} end +# Outer constructor: takes TR from the ranges and lets the inner constructor convert the remaining +# fields (in particular the integer factor_ab). The auto-generated constructor of a *parametric* +# struct matches the declared field types exactly and would not perform that conversion. +function TempArgs(rowi,coli,ranges::NamedTuple{(:nua, :nub, :NUa, :NUb),NTuple{4,TR}},norm4,mij_arr,sa,JsSa,sb,JsSb,JlLa,JlLb,la,La,lb,Lb,Lsum,avals_new,bvals_new,factor_ab,avals,bvals) where {TR} + TempArgs{TR}(rowi,coli,ranges,norm4,mij_arr,sa,JsSa,sb,JsSb,JlLa,JlLb,la,La,lb,Lb,Lsum,avals_new,bvals_new,factor_ab,avals,bvals) +end + -struct FillStruct{T} +# T : element type of H, V and everything that can pick up a complex-scaling phase +# (complex as soon as complex_scaling OR complex ranges are active) +# TR : element type of the ranges and of the quantities derived from them alone +# (complex only for complex ranges; the ISGL complex-scaling rotates the potential, not the ranges) +struct FillStruct{T,TR} fij_arr::Matrix{Float64} Kij_arr::Matrix{Float64} - w_arr_kine::OffsetVector{Float64, Vector{Float64}} + w_arr_kine::OffsetVector{TR, Vector{TR}} wn_interpol_arr::OffsetVector{T, Vector{T}} - kij_arr::OffsetArray{Float64, 3, Array{Float64, 3}}#OffsetArray{Float64} - gij_arr::OffsetArray{Float64, 3, Array{Float64, 3}}#OffsetArray{Float64} + kij_arr::OffsetArray{TR, 3, Array{TR, 3}} + gij_arr::OffsetArray{TR, 3, Array{TR, 3}} T::Matrix{T} V::Matrix{T} - S::Matrix{Float64} - temp_args_arr::Vector{TempArgs} + S::Matrix{TR} + temp_args_arr::Vector{TempArgs{TR}} temp_fill_mat::Matrix{T} wn_obs_interpol_arr::OffsetVector{Float64, Vector{Float64}} end -struct ResultStruct{T} - energies_arr::Vector{T} - wavefun_arr::Matrix{T} +# TE: element type of the energies. Complex ranges alone leave H and S hermitian, so the +# eigenvalues stay real; only complex scaling makes them genuinely complex. +struct ResultStruct{TE,TW} + energies_arr::Vector{TE} + wavefun_arr::Matrix{TW} centobs_output::Array{Float64} R2_output::Matrix{Float64} end -function preallocate_data(phys_params,num_params,observ_params,size_params,complex_scaling::Bool) - if !complex_scaling - TT = Float64 - elseif complex_scaling - TT = ComplexF64 - end - +function preallocate_data(phys_params,num_params,observ_params,size_params,complex_scaling::Bool,complex_ranged_r::Bool=false,complex_ranged_R::Bool=false) + complex_ranged = complex_ranged_r || complex_ranged_R + + # TT: H, V, and everything that can carry the complex-scaling phase + # TR: ranges, norms, overlap S, gij/kij (the ISGL complex-scaling rotates the potential, not the ranges) + # TE: energies. With complex ranges alone, H and S stay hermitian -> real eigenvalues. + TT = (complex_scaling || complex_ranged) ? ComplexF64 : Float64 + TR = complex_ranged ? ComplexF64 : Float64 + TE = complex_scaling ? ComplexF64 : Float64 + #Destructing Structs: (;J_tot) = phys_params (;gem_params,kmax_interpol) = num_params (;nmax,Nmax) = gem_params + + # complex-ranged basis functions come in conjugate pairs, so the corresponding coordinate + # carries twice as many ranges. The two coordinates are independent. + nmax_eff = complex_ranged_r ? 2*nmax : nmax + Nmax_eff = complex_ranged_R ? 2*Nmax : Nmax (;nintmax,nbasis_total,nlL,nl,maxlmax,maximax,maxkmax,maxobs,JlL_complete) = size_params (;stateindices,centobs_arr,R2_arr) = observ_params @@ -124,13 +150,13 @@ function preallocate_data(phys_params,num_params,observ_params,size_params,compl murR_arr = MMatrix{2,3,Float64,6}(zeros(2,3)); # reduced masses # Allocate ranges - nu_arr = Vector{Float64}(undef, nmax) - NU_arr = Vector{Float64}(undef, Nmax) + nu_arr = Vector{TR}(undef, nmax_eff) + NU_arr = Vector{TR}(undef, Nmax_eff) #ranges = (nu_arr, NU_arr) - + # Allocate norms - norm_arr = OffsetMatrix{Float64}(undef, 0:maxlmax, nmax)*0.0 # one norm_arr for all l,n combinations. changed to OffsetArray for all values of l=0:maxlmax - NORM_arr = OffsetMatrix{Float64}(undef, 0:maxlmax, Nmax)*0.0 # one norm_arr for all L,N combinations + norm_arr = OffsetMatrix{TR}(undef, 0:maxlmax, nmax_eff)*zero(TR) # one norm_arr for all l,n combinations. changed to OffsetArray for all values of l=0:maxlmax + NORM_arr = OffsetMatrix{TR}(undef, 0:maxlmax, Nmax_eff)*zero(TR) # one norm_arr for all L,N combinations # ISGL arrays: Clmk_arr = OffsetArray{Float64}(undef, 0:maxlmax, -maxlmax:maxlmax, maxkmax)*0.0 @@ -149,7 +175,9 @@ function preallocate_data(phys_params,num_params,observ_params,size_params,compl w_interpol_arr=OffsetArray{interpoltypeC}(undef,3,nintmax,0:2*maxlmax,0:2*maxlmax) # penultimate dimensionality for different Lsum values!; nintmax for (maximum) number of interactions #aac=zeros(27);gac=zeros(27);abc=zeros(27);gbc=zeros(27) Ainv_arr_kine = OffsetArray{Float64}(undef,0:2*maxlmax,0:2*maxlmax)*0.0 # penultimate dimensionality for different Lsum values! - w_arr_kine = OffsetArray{Float64}(undef,0:2*maxlmax+1)*0.0 + v_pow = zeros(TT,2*maxlmax+1) # scratch for the power-law shoulder solve (alpha-independent); NOT an OffsetArray, as for v_arr + w_pow_arr = OffsetArray(zeros(TT,3,nintmax,2*maxlmax+1,2*maxlmax+1),1:3,1:nintmax,0:2*maxlmax,0:2*maxlmax) # analytic shoulder-weights for power-law interactions; same indexing as w_interpol_arr, but plain numbers (no interpolation needed) + w_arr_kine = OffsetArray{TR}(undef,0:2*maxlmax+1)*zero(TR) wn_interpol_arr = OffsetArray{TT}(undef,0:2*maxlmax)*0.0 # for the interpolated wn_values that are actually used #for observables (in range-interpolation) v_obs_arr = zeros(kmax_interpol,2*maxlmax+1) @@ -160,15 +188,15 @@ function preallocate_data(phys_params,num_params,observ_params,size_params,compl # arrays for shoulder method (used in matrix-element calculations) fij_arr = zeros(3,4) Kij_arr = zeros(3,4) - kij_arr = OffsetArray{Float64}(undef,3,4,0:2*maxlmax)*0.0 - gij_arr = OffsetArray{Float64}(undef,3,4,0:2*maxlmax)*0.0 - S = zeros(nbasis_total,nbasis_total) #Matrix{Float64}(undef, nbasis_total, nbasis_total) + kij_arr = OffsetArray{TR}(undef,3,4,0:2*maxlmax)*zero(TR) + gij_arr = OffsetArray{TR}(undef,3,4,0:2*maxlmax)*zero(TR) + S = zeros(TR,nbasis_total,nbasis_total) # hermitian (not just symmetric) for complex ranges T = zeros(TT,nbasis_total,nbasis_total) #Matrix{Float64}(undef, nbasis_total, nbasis_total) V = zeros(TT,nbasis_total,nbasis_total);#Matrix{Float64}(undef, nlL * nmax * Nmax, nlL * nmax * Nmax) # for results: - energies_arr = zeros(TT,nbasis_total);# Vector{TT}(undef, nbasis_total) # changed to zeros to avoid bad behavior due to undef values not occupied thresholding (eigen2step) + energies_arr = zeros(TE,nbasis_total);# changed to zeros to avoid bad behavior due to undef values not occupied thresholding (eigen2step) wavefun_arr = Matrix{TT}(undef, nbasis_total, nbasis_total) centobs_output = Array{Float64}(undef, 3, maxobs, lastindex(stateindices)) R2_output = zeros(3, lastindex(stateindices)) @@ -181,17 +209,18 @@ function preallocate_data(phys_params,num_params,observ_params,size_params,compl temp_D2=zeros(3)#MVector{3, ComplexF64}(zeros(3)) # temporary arrays for filling: function arguments and matrix - ntot = Int64((nbasis_total^2+nbasis_total)/2) # only for lower triangular! - # Define the type of the tuple for the function arguments - tuple_type = NamedTuple{(:rowi, :coli, :ranges, :norm4, :mij_arr, :sa, :JsSa, :sb, :JsSb, :JlLa, :JlLb, :la, :La, :lb, :Lb, :Lsum, :avals_new, :bvals_new, :factor_ab, :avals, :bvals),Tuple{Int64,Int64,NamedTuple{(:nua, :nub, :NUa, :NUb),Tuple{Float64,Float64,Float64,Float64}},Float64,Array{SArray{Tuple{6},Int64,1,6},1},Float64,Float64,Float64,Float64,Int64,Int64,Int64,Int64,Int64,Int64,Int64,Array{Int64,1},Vector{Int64},Int64,Vector{Int64},Vector{Int64}}} - temp_args_arr = Vector{TempArgs}(undef, ntot) + # Combined complex ranges and complex scaling leave H neither symmetric nor hermitian, + # so the full matrix has to be filled instead of only the lower triangle. + fill_full = complex_ranged && complex_scaling + ntot = fill_full ? Int64(nbasis_total^2) : Int64((nbasis_total^2+nbasis_total)/2) + temp_args_arr = Vector{TempArgs{TR}}(undef, ntot) temp_fill_mat = zeros(TT,nbasis_total, nbasis_total) # now constructing Structs for different steps in the program: precomp_arrs = PrecomputeStruct(gamma_dict,cleb_arr,spintrafo_dict,spinoverlap_dict,global6j_dict,facsymm_dict,jmat,murR_arr,nu_arr,NU_arr,norm_arr,NORM_arr,Clmk_arr,Dlmk_arr,S_arr,SSO_arr) temp_arrs = TempStruct(temp_clmk,temp_dlmk,temp_S,temp_D1,temp_D2) - interpol_arrs = InterpolationStruct(alpha_arr,v_arr,A_mat,w_arr,w_interpol_arr,Ainv_arr_kine,v_obs_arr,w_obs_arr,w_obs_interpol_arr) + interpol_arrs = InterpolationStruct(alpha_arr,v_arr,A_mat,w_arr,w_interpol_arr,Ainv_arr_kine,v_pow,w_pow_arr,v_obs_arr,w_obs_arr,w_obs_interpol_arr) fill_arrs = FillStruct(fij_arr,Kij_arr,w_arr_kine,wn_interpol_arr,kij_arr,gij_arr,T,V,S,temp_args_arr,temp_fill_mat,wn_obs_interpol_arr) result_arrs = ResultStruct(energies_arr,wavefun_arr,centobs_output,R2_output) diff --git a/src/ISGL-3body/precomputation.jl b/src/ISGL-3body/precomputation.jl index 7b97fea..b45a778 100644 --- a/src/ISGL-3body/precomputation.jl +++ b/src/ISGL-3body/precomputation.jl @@ -2,10 +2,10 @@ # for spins we added transformation coefficients from spins in jacobi-set a,b to c (and a to b) -function precompute_ISGL(phys_params,num_params,size_params,precomp_arrs,temp_arrs) - #Destructing Structs: +function precompute_ISGL(phys_params,num_params,size_params,precomp_arrs,temp_arrs,complex_ranged_r::Bool=false,complex_ranged_R::Bool=false) + #Destructing Structs: (;masses,J_tot,spins) = phys_params - (;lmax,Lmax,gem_params) = num_params + (;lmax,Lmax,gem_params,complex_range_freq) = num_params (;nmax,Nmax,r1,rnmax,R1,RNmax) = gem_params (;cvals,s_arr,abI,factor_bf,s_complete,JsS_arr,JsS_complete,JlL_arr,JlL_complete,lL_complete,l_complete,nl,imax_dict,mij_arr_dict,kmax_dict,imaxSO_dict,mijSO_arr_dict,so_indices) = size_params (;gamma_dict,cleb_arr,spintrafo_dict,spinoverlap_dict,global6j_dict,facsymm_dict,jmat,murR_arr,nu_arr,NU_arr,norm_arr,NORM_arr,Clmk_arr,Dlmk_arr,S_arr,SSO_arr) = precomp_arrs @@ -31,7 +31,7 @@ function precompute_ISGL(phys_params,num_params,size_params,precomp_arrs,temp_ar precompute_murR(murR_arr,masses) - precompute_ranges(nu_arr,NU_arr,r1,rnmax,nmax,R1,RNmax,Nmax) + precompute_ranges(nu_arr,NU_arr,r1,rnmax,nmax,R1,RNmax,Nmax,complex_ranged_r,complex_ranged_R,complex_range_freq) precompute_norms(norm_arr,NORM_arr,nu_arr,NU_arr,nl,l_complete,gamma_dict) precompute_spintrafo(spins,s_arr,JsS_arr,spintrafo_dict) @@ -124,10 +124,23 @@ end ### ranges ### -function precompute_ranges(nu_arr,NU_arr,r1,rnmax,nmax,R1,RNmax,Nmax) - # returns the array nu_arr[n=1:nmax] for each value of n. same for NU_arr[N=1:Nmax] - nu_arr .= buildnu(r1,rnmax,nmax,nu_arr) - NU_arr .= buildnu(R1,RNmax,Nmax,NU_arr) +function precompute_ranges(nu_arr,NU_arr,r1,rnmax,nmax,R1,RNmax,Nmax,complex_ranged_r::Bool=false,complex_ranged_R::Bool=false,complex_range_freq=0.0) + # fills nu_arr[n=1:nmax] with the geometric sequence of ranges. same for NU_arr[N=1:Nmax]. + # For a complex-ranged coordinate the array is twice as long: the first half carries + # nu*(1+i*omega), the second half its complex conjugate. + fillranges!(nu_arr,r1,rnmax,nmax,complex_ranged_r,complex_range_freq) + fillranges!(NU_arr,R1,RNmax,Nmax,complex_ranged_R,complex_range_freq) +end + +# note: nmax is always the *undoubled* number of ranges here, since the geometric sequence +# (and in particular its exponent 2*(n-1)/(nmax-1)) must not be affected by the doubling. +@views function fillranges!(nu_arr,r1,rnmax,nmax,complex_ranged::Bool,complex_range_freq) + buildnu(r1,rnmax,nmax,nu_arr) # fills nu_arr[1:nmax] in place; leaves any further entries untouched + if complex_ranged + nu_arr[1:nmax] .*= (1 + complex_range_freq*im) + nu_arr[nmax+1:2*nmax] .= conj.(nu_arr[1:nmax]) + end + return nu_arr end @views function buildnu(r1,rnmax,nmax,nu_arr) @@ -141,7 +154,11 @@ end function precompute_norms(norm_arr,NORM_arr,nu_arr,NU_arr,nl,l_complete,gamma_dict) # returns the array norm_arr[l,n] for each combination of n,l. same for NORM_arr - norm(nu,l) = (2*(2*nu)^(l+3/2)/gamma_dict[l+3/2])^(1/2); + # Equivalent to (2*(2*nu)^(l+3/2)/Gamma(l+3/2))^(1/2) for real, positive nu, but written so that + # the fractional power is taken directly on 2*nu. For a complex range, arg(2*nu) = atan(omega) is + # small, so this is the continuous continuation of the real-range norm; the original nesting could + # push arg((2*nu)^(l+3/2)) past pi for larger l and then pick the wrong branch in the outer sqrt. + norm(nu,l) = sqrt(2/gamma_dict[l+3/2]) * (2*nu)^((2*l+3)/4); for n = 1:lastindex(nu_arr) for l in l_complete #lindex = 1:nl diff --git a/src/ISGL-3body/sanitycheck.jl b/src/ISGL-3body/sanitycheck.jl index 88b9a5e..bf86a3b 100644 --- a/src/ISGL-3body/sanitycheck.jl +++ b/src/ISGL-3body/sanitycheck.jl @@ -1,8 +1,8 @@ # function to check the inputs for the ISGL program of sanity: -function sanity_checks(phys_params) +function sanity_checks(phys_params,complex_ranged::Bool=false,observ_params=nothing) (;masses,species,interactions,J_tot,parity) = phys_params - + if (lastindex(masses) !=3) || lastindex(species) !=3 println("masses and/or species have wrong size, must be 3") error_code = 1 @@ -21,6 +21,46 @@ function sanity_checks(phys_params) return error_code end + # power-law interactions: the radial integral int r^(2n+2) V(r) exp(-alpha r^2) dr + # with n=0 converges only for p > -3. + for c in 1:lastindex(interactions) + for vint in interactions[c] + if vint isa PowerLawPotential && vint.p <= -3.0 + println("PowerLawPotential: p = $(vint.p) is too singular for ISGL; requires p > -3") + error_code = 4 + return error_code + end + end + end + + # complex-ranged basis functions: + # The generic path for a central potential obtains its radial integral numerically on a real + # mesh of effective ranges and then interpolates at log(etaprc). Complex ranges make etaprc + # complex, which that lookup cannot represent. Only the analytically treated potential types + # (GaussianPotential, PowerLawPotential) are therefore supported with complex ranges; they need + # neither numerical integration nor interpolation. + if complex_ranged + for c in 1:lastindex(interactions) + for vint in interactions[c] + if !(vint isa GaussianPotential || vint isa PowerLawPotential) + println("Complex-ranged basis functions are only supported for analytically treated potentials (GaussianPotential, PowerLawPotential). Got $(typeof(vint)) in Jacobi set c=$c.") + error_code = 5 + return error_code + end + end + end + + # observables go through the same interpolation machinery and are equally unsupported + if !isnothing(observ_params) + (;centobs_arr,R2_arr) = observ_params + if any(.!isempty.(centobs_arr)) || any(R2_arr .!= 0) + println("Observables (centobs_arr, R2_arr) are not supported together with complex-ranged basis functions.") + error_code = 6 + return error_code + end + end + end + # tests required: ok, even if fasb or fasf is empty? if allequal(masses[fasb]) == false || allequal(masses[fasf]) == false println("Problem with symmetrization: species does not fit to m_arr") diff --git a/src/ISGL-3body/size_estimate.jl b/src/ISGL-3body/size_estimate.jl index b65bd3f..9141840 100644 --- a/src/ISGL-3body/size_estimate.jl +++ b/src/ISGL-3body/size_estimate.jl @@ -14,6 +14,8 @@ struct SizeParams{T<:Number} cvals::Vector{Int64} gauss_indices::Vector{Vector{Int}} gaussopt_arr::Vector{Vector{Tuple{Float64,T}}} + pow_indices::Vector{Vector{Int}} + powopt_arr::Vector{Vector{Tuple{Float64,Float64}}} central_indices::Vector{Vector{Int}} so_indices::Vector{Vector{Int}} nint_arr::Vector{Int64} @@ -50,8 +52,8 @@ struct SizeParams{T<:Number} maxobs::Int64 end -function size_estimate(phys_params,num_params,observ_params,complex_scaling::Bool) - +function size_estimate(phys_params,num_params,observ_params,complex_scaling::Bool,complex_ranged_r::Bool=false,complex_ranged_R::Bool=false) + # input interpretation: (;masses,species,interactions,J_tot,parity,spins) = phys_params (;lmax,Lmax,gem_params,complex_scaling_angle,complex_range_freq,mu0,c_shoulder,kmax_interpol,lmin,Lmin) = num_params @@ -62,7 +64,7 @@ function size_estimate(phys_params,num_params,observ_params,complex_scaling::Boo cvals = findall(isempty.(interactions) .==0 ) # consider only values for c where there are interactions (any type) # number of interactions per Jacobi-set c: - gauss_indices, gaussopt_arr, central_indices, so_indices, nint_arr, nintmax = index_interaction_types(interactions,complex_scaling, complex_scaling_angle) + gauss_indices, gaussopt_arr, pow_indices, powopt_arr, central_indices, so_indices, nint_arr, nintmax = index_interaction_types(interactions,complex_scaling, complex_scaling_angle) # box sizes, indices and factors for symmetrization abvals_arr,groupindex_arr,nboxes,abI,factor_bf = abc_size(cvals,species) @@ -73,7 +75,10 @@ function size_estimate(phys_params,num_params,observ_params,complex_scaling::Boo lL_nested,lL_complete,l_complete = lLcoupl(J_tot,parity,cvals,species,spins,s_arr,JsS_arr,JlL_arr,lmin,Lmin,lmax,Lmax) # box sizes = number of basis functions in each box; starts,ends = indices for boxes within big matrix; bvalsdiag = simplification for identical particles - box_size_arr,nbasis_total = boxsize_fun(groupindex_arr,abvals_arr,lL_nested,nmax,Nmax) + # a complex-ranged coordinate doubles its number of basis functions (nu and its conjugate) + nmax_eff = complex_ranged_r ? 2*nmax : nmax + Nmax_eff = complex_ranged_R ? 2*Nmax : Nmax + box_size_arr,nbasis_total = boxsize_fun(groupindex_arr,abvals_arr,lL_nested,nmax_eff,Nmax_eff) starts = [1; cumsum(box_size_arr[1:end-1]) .+ 1] ends = cumsum(box_size_arr) #bvalsdiag = [[abvals_arr[boxR][1]] for boxR in groupindex_arr] @@ -99,7 +104,7 @@ function size_estimate(phys_params,num_params,observ_params,complex_scaling::Boo maxobs = maximum(lastindex.(centobs_arr)) # max number of observables # Constructing Struct (collective data structure size_params with all the size parameters) - size_params = SizeParams(abvals_arr,cvals,gauss_indices,gaussopt_arr,central_indices,so_indices,nint_arr,nintmax,groupindex_arr,nboxes,abI,factor_bf,box_size_arr,nbasis_total,starts,ends,bvalsdiag,s_arr,JsS_arr,s_complete,JsS_complete,JlL_arr,JlL_complete,lL_nested,lL_complete,l_complete,nlL,nl,maxlmax,imax_dict,imaxSO_dict,maximax,mij_arr_dict,mijSO_arr_dict,kmax_dict,maxkmax,maxobs) + size_params = SizeParams(abvals_arr,cvals,gauss_indices,gaussopt_arr,pow_indices,powopt_arr,central_indices,so_indices,nint_arr,nintmax,groupindex_arr,nboxes,abI,factor_bf,box_size_arr,nbasis_total,starts,ends,bvalsdiag,s_arr,JsS_arr,s_complete,JsS_complete,JlL_arr,JlL_complete,lL_nested,lL_complete,l_complete,nlL,nl,maxlmax,imax_dict,imaxSO_dict,maximax,mij_arr_dict,mijSO_arr_dict,kmax_dict,maxkmax,maxobs) return size_params end @@ -131,22 +136,26 @@ end function index_interaction_types(interactions,complex_scaling::Bool, complex_scaling_angle) gauss_indices = [Int[] for _ in 1:3] gaussopt_arr = [Tuple{Float64,Float64}[] for _ in 1:3] + pow_indices = [Int[] for _ in 1:3] + powopt_arr = [Tuple{Float64,Float64}[] for _ in 1:3] central_indices = [Int[] for _ in 1:3] so_indices = [Int[] for _ in 1:3] nint_arr = zeros(Int64,3) - pushindexpotentialtype!(v::Function, central_indices, gauss_indices, so_indices, i) = push!(central_indices, i) # treat function as a central potential - pushindexpotentialtype!(v::CentralPotential, central_indices, gauss_indices, so_indices, i) = push!(central_indices, i) - #pushindexpotentialtype!(v::SpinOrbitPotential, central_indices, gauss_indices, so_indices, i) = push!(so_indices, i) # postponed to future version - pushindexpotentialtype!(v::GaussianPotential, central_indices, gauss_indices, so_indices, i) = push!(gauss_indices, i) + pushindexpotentialtype!(v::Function, central_indices, gauss_indices, pow_indices, so_indices, i) = push!(central_indices, i) # treat function as a central potential + pushindexpotentialtype!(v::CentralPotential, central_indices, gauss_indices, pow_indices, so_indices, i) = push!(central_indices, i) + #pushindexpotentialtype!(v::SpinOrbitPotential, central_indices, gauss_indices, pow_indices, so_indices, i) = push!(so_indices, i) # postponed to future version + pushindexpotentialtype!(v::GaussianPotential, central_indices, gauss_indices, pow_indices, so_indices, i) = push!(gauss_indices, i) + pushindexpotentialtype!(v::PowerLawPotential, central_indices, gauss_indices, pow_indices, so_indices, i) = push!(pow_indices, i) for c in 1:3 gauss_indices[c] = Int[] + pow_indices[c] = Int[] central_indices[c] = Int[] so_indices[c] = Int[] for (i, v) in enumerate(interactions[c]) - pushindexpotentialtype!(v, central_indices[c], gauss_indices[c], so_indices[c], i) + pushindexpotentialtype!(v, central_indices[c], gauss_indices[c], pow_indices[c], so_indices[c], i) if i in gauss_indices[c] # if this is a Gaussian potential push!(gaussopt_arr[c], (v.v0, v.mu_g)) # store the parameters of the Gaussian potential @@ -154,6 +163,12 @@ function index_interaction_types(interactions,complex_scaling::Bool, complex_sca push!(gaussopt_arr[c], (NaN, NaN)) # if not a Gaussian potential, store NaN. This is necessary to keep the order of indices. end + if i in pow_indices[c] # if this is a power-law potential + push!(powopt_arr[c], (v.v0, v.p)) # store the parameters of the power-law potential + else + push!(powopt_arr[c], (NaN, NaN)) # if not a power-law potential, store NaN. This is necessary to keep the order of indices. + end + end nint_arr[c] = lastindex(interactions[c]) end @@ -161,8 +176,10 @@ function index_interaction_types(interactions,complex_scaling::Bool, complex_sca nintmax = maximum(nint_arr) gaussopt_arrC = csmgaussopt(gaussopt_arr, complex_scaling, complex_scaling_angle) # adjust to complex values when complex_scaling=true + # note: powopt_arr needs no such adjustment; the complex-scaling factor for a power law is a single + # global csmfac^(-p) which is applied in precompute_w (see interpolationNshoulder.jl). - return gauss_indices, gaussopt_arrC, central_indices, so_indices, nint_arr, nintmax + return gauss_indices, gaussopt_arrC, pow_indices, powopt_arr, central_indices, so_indices, nint_arr, nintmax end diff --git a/src/common/potentialtypes.jl b/src/common/potentialtypes.jl index 4388243..cf3b999 100644 --- a/src/common/potentialtypes.jl +++ b/src/common/potentialtypes.jl @@ -54,6 +54,53 @@ function (gp::GaussianPotential)(r) end +""" + PowerLawPotential(v0::Float64, p::Float64) +A concrete implementation of `PotentialFunction` that represents a power-law potential: +```math +V(r) = v_0 |r|^{p} +``` + +Note the **absolute value**: the potential is defined via `|r|`, i.e. it is an even +function. This matters only in 1D, where the coordinate can be negative; there, +`PowerLawPotential(v0,p)` describes `v0*|x|^p`. Odd potentials are not supported. +In 2D and 3D, `r >= 0` anyway and `|r| = r`. + +Matrix elements are evaluated analytically (closed-form radial integration), so +neither numerical integration nor range-interpolation is needed. Prominent cases +are the Coulomb interaction (`p=-1`) and the harmonic oscillator (`p=2`). + +The matrix elements are finite only if `p` is not too negative. The bound depends +on the dimension and on the angular momentum, and is therefore **not** checked +here but in the individual modules: +- `GEM2B`: requires `p > -2*lmax - dim`, e.g. `p > -1` for `lmax=0` in 1D, + so a pure `1/|x|` is already divergent in 1D. +- `ISGL`: requires `p > -3`. + +# Arguments: +- `v0::Float64`: The strength of the potential. +- `p::Float64`: The exponent of the power law. +""" +struct PowerLawPotential <: PotentialFunction + v0::Float64 + p::Float64 +end + +# convenience: allow integer exponents, e.g. PowerLawPotential(-1.0,-1) +PowerLawPotential(v0::Real, p::Real) = PowerLawPotential(Float64(v0), Float64(p)) + +""" + function (pp::PowerLawPotential)(r) + +Evaluates the power-law potential at a given radial distance `r`. Uses `abs(r)`, +see the note on 1D in [`PowerLawPotential`](@ref). +""" +function (pp::PowerLawPotential)(r) + pp.v0 * abs(r)^pp.p +end + + + # postponed to future version #= """ SpinOrbitPotential(f::Function) diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..a0c58be --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,6 @@ +[deps] +FewBodyToolkit = "28c40a4d-cd08-48a0-9e99-cc2df8099be4" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +Test = "1" diff --git a/test/test2B1D.jl b/test/test2B1D.jl index 4919d84..bb2e342 100644 --- a/test/test2B1D.jl +++ b/test/test2B1D.jl @@ -38,3 +38,32 @@ energies_arr = GEM2B.GEM2B_solve(phys_params,num_paramsC;complex_scaling=true) num_paramsC = make_num_params2B(;gem_params,complex_scaling_angle=5.0) energies_arr = GEM2B.GEM2B_solve(phys_params,num_paramsC;complex_scaling=true,complex_ranged=true) @test all(isapprox.(real.(energies_arr[1:4]), exact_results; atol=1e-3)) + +# PowerLawPotential (1D): analytic treatment of V(x) = v0*|x|^p +# 1D harmonic oscillator; lmax=0 selects the even states E=(2n+1/2)*omega, +# lmax=1 the odd ones E=(2n+3/2)*omega +omega_1d = 0.6 +v_ho_pow = PowerLawPotential(0.5*mur*omega_1d^2, 2.0) +v_ho_cent(r) = 0.5*mur*omega_1d^2*r^2 +gp_1d = (;nmax=24,r1=0.2,rnmax=12.0) +np_1d = make_num_params2B(;gem_params=gp_1d) +for (l,off) in [(0,0.5),(1,1.5)] + pp_pow = make_phys_params2B(;mur,interactions=[v_ho_pow], dim=1,lmax=l,lmin=l) + pp_cent = make_phys_params2B(;mur,interactions=[v_ho_cent],dim=1,lmax=l,lmin=l) + e_pow = GEM2B.GEM2B_solve(pp_pow, np_1d) + e_cent = GEM2B.GEM2B_solve(pp_cent,np_1d) + exact_1d = [(2*n+off)*omega_1d for n=0:3] + @test all(isapprox.(e_pow[1:4], exact_1d; atol=1e-2)) + @test all(isapprox.(e_pow[1:4], e_cent[1:4]; rtol=1e-6)) +end + +# |x|^p with non-integer p: analytic vs numerical (the numerical path needs abs() explicitly) +v_pl_pow = PowerLawPotential(0.8,1.5) +v_pl_cent(r) = 0.8*abs(r)^1.5 +@test all(isapprox.(GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[v_pl_pow],dim=1),np_1d)[1:4], + GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[v_pl_cent],dim=1),np_1d)[1:4]; rtol=1e-8)) + +# validity check: in 1D with lmax=0 the bound is p > -1, so a pure 1/|x| already diverges +@test_throws ErrorException GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[PowerLawPotential(-1.0,-1.0)],dim=1),np_1d) +# for lmax=1 the bound is p > -3, so p=-1 is fine there +@test all(isfinite.(GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[PowerLawPotential(-1.0,-1.0)],dim=1,lmax=1,lmin=1),np_1d)[1:2])) diff --git a/test/test2B2D.jl b/test/test2B2D.jl index 62d5780..13f73cd 100644 --- a/test/test2B2D.jl +++ b/test/test2B2D.jl @@ -10,7 +10,7 @@ v_ho(r) = 0.5*mur*omega^2*r^2 phys_params = make_phys_params2B(;mur,interactions=[v_ho],dim=2) # numerical parameters: -gem_params = (;nmax=14,r1=0.82,rnmax=10.62) # gem_params +gem_params = (;nmax=24,r1=0.82,rnmax=10.62) # gem_params num_params = make_num_params2B(;gem_params) @@ -20,7 +20,7 @@ exact_results = ([2*i for i=0:15] .+ 1) .*omega # 1. Standard inputs: csm_bool = 0, cr_bool = 0 energies_arr = GEM2B.GEM2B_solve(phys_params,num_params) -@test all(isapprox.(energies_arr[1:4], exact_results[1:4]; atol=1e-3)) +@test all(isapprox.(energies_arr[1:4], exact_results[1:4]; atol=1e-2)) # 2. Using complex-ranged basis functions: complex_scaling = false, complex_ranged = true num_paramsCR = make_num_params2B(;gem_params=(nmax=7,r1=1.3398861224184124,rnmax=4.4150781608810705)) @@ -29,9 +29,28 @@ energies_arr = GEM2B.GEM2B_solve(phys_params,num_paramsCR;complex_ranged=true) # 3. Complex scaling (angle 0°) should have no effect: complex_scaling = true energies_arr = GEM2B.GEM2B_solve(phys_params,num_params;complex_scaling=true) -@test all(isapprox.(real.(energies_arr[1:4]), exact_results[1:4]; atol=1e-3)) +@test all(isapprox.(real.(energies_arr[1:4]), exact_results[1:4]; atol=1e-2)) # 4. Finite complex scaling angle (1°) should have very little effect on the bound states num_paramsC = make_num_params2B(;gem_params,complex_scaling_angle=1.0) energies_arr = GEM2B.GEM2B_solve(phys_params,num_paramsC;complex_scaling=true) -@test all(isapprox.(real.(energies_arr[1:4]), exact_results[1:4]; atol=1e-3)) +@test all(isapprox.(real.(energies_arr[1:4]), exact_results[1:4]; atol=1e-2)) + +# PowerLawPotential (2D): analytic treatment of V(r) = v0*|r|^p +# harmonic oscillator (p=2) against the exact 2D spectrum, and against the numerical path +v_ho_pow = PowerLawPotential(0.5*mur*omega^2, 2.0) +pp_ho_pow = make_phys_params2B(;mur,interactions=[v_ho_pow],dim=2) +e_ho_pow = GEM2B.GEM2B_solve(pp_ho_pow,num_params) +@test all(isapprox.(e_ho_pow[1:4], exact_results[1:4]; atol=1e-2)) +@test all(isapprox.(e_ho_pow[1:4], GEM2B.GEM2B_solve(phys_params,num_params)[1:4]; rtol=1e-8)) + +# attractive p=-1 in 2D: analytic vs numerical (allowed, since the bound is p > -2 for lmax=0) +v_c_pow = PowerLawPotential(-1.0,-1.0) +v_c_cent(r) = -1.0/abs(r) +gp_c = (;nmax=12,r1=0.1,rnmax=25.0) +np_c = make_num_params2B(;gem_params=gp_c) +@test all(isapprox.(GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[v_c_pow],dim=2),np_c)[1:3], + GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[v_c_cent],dim=2),np_c)[1:3]; rtol=1e-6)) + +# validity check: in 2D with lmax=0 the bound is p > -2 +@test_throws ErrorException GEM2B.GEM2B_solve(make_phys_params2B(;mur,interactions=[PowerLawPotential(1.0,-2.0)],dim=2),num_params) diff --git a/test/test2B3D.jl b/test/test2B3D.jl index 8f93e91..d9c7dd1 100644 --- a/test/test2B3D.jl +++ b/test/test2B3D.jl @@ -84,3 +84,44 @@ np_optim = make_num_params2B(;gem_params=(;nmax=8, r1=0.5, rnmax=20.0)) result_optim = GEM_Optim_2B(phys_params, np_optim, 1) @test length(result_optim) == 3 # [r1_opt, rnmax_opt, energy] @test isapprox(result_optim[3], -0.5; atol=1e-2) # hydrogen ground state E = -0.5 + +# 12. PowerLawPotential (3D): analytic treatment of V(r) = v0*|r|^p +# 12a. Coulomb (p=-1) against the exact results and against the numerical path +v_coulomb_pow = PowerLawPotential(-Z,-1.0) +pp_pow = make_phys_params2B(;interactions=[v_coulomb_pow]) +e_pow = GEM2B.GEM2B_solve(pp_pow,num_params) +@test all(isapprox.(e_pow[1:4], exact_results; atol=1e-3)) +@test all(isapprox.(e_pow[1:4], GEM2B.GEM2B_solve(phys_params,num_params)[1:4]; rtol=1e-8)) + +# 12b. harmonic oscillator (p=2) against the exact 3D spectrum E=(2n+l+3/2)*omega +omega_ho = 0.7 +v_ho_pow = PowerLawPotential(0.5*1.0*omega_ho^2, 2.0) +v_ho_cent(r) = 0.5*1.0*omega_ho^2*r^2 +gp_ho = (;nmax=24,r1=0.2,rnmax=12.0) +np_ho = make_num_params2B(;gem_params=gp_ho) +for l in [0,1,2] + pp_ho_pow = make_phys_params2B(;interactions=[v_ho_pow], lmax=l, lmin=l) + pp_ho_cent = make_phys_params2B(;interactions=[v_ho_cent], lmax=l, lmin=l) + e_ho_pow = GEM2B.GEM2B_solve(pp_ho_pow, np_ho) + e_ho_cent = GEM2B.GEM2B_solve(pp_ho_cent,np_ho) + exact_ho = [(2*n+l+1.5)*omega_ho for n=0:3] + @test all(isapprox.(e_ho_pow[1:4], exact_ho; atol=1e-2)) + @test all(isapprox.(e_ho_pow[1:4], e_ho_cent[1:4]; rtol=1e-8)) +end + +# 12c. non-integer exponent +v_pl_pow = PowerLawPotential(0.9,1.5) +v_pl_cent(r) = 0.9*abs(r)^1.5 +@test all(isapprox.(GEM2B.GEM2B_solve(make_phys_params2B(;interactions=[v_pl_pow]),num_params)[1:4], + GEM2B.GEM2B_solve(make_phys_params2B(;interactions=[v_pl_cent]),num_params)[1:4]; rtol=1e-8)) + +# 12d. complex scaling at finite angle: analytic vs numerical +np_csm = make_num_params2B(;gem_params,theta_csm=8.0) +@test all(isapprox.(GEM2B.GEM2B_solve(pp_pow,np_csm;complex_scaling=true)[1:4], + GEM2B.GEM2B_solve(phys_params,np_csm;complex_scaling=true)[1:4]; atol=1e-6)) + +# 12e. validity check: in 3D with lmax=0 the bound is p > -3 +@test_throws ErrorException GEM2B.GEM2B_solve(make_phys_params2B(;interactions=[PowerLawPotential(1.0,-3.0)]),num_params) +@test_throws ErrorException GEM2B.GEM2B_solve(make_phys_params2B(;interactions=[PowerLawPotential(1.0,-4.0)]),num_params) +# but for lmax=1 the bound is p > -5, so p=-4 is fine there +@test all(isfinite.(GEM2B.GEM2B_solve(make_phys_params2B(;interactions=[PowerLawPotential(1.0,-4.0)],lmax=1,lmin=1),num_params)[1:2])) diff --git a/test/testISGL.jl b/test/testISGL.jl index 9934a81..1112662 100644 --- a/test/testISGL.jl +++ b/test/testISGL.jl @@ -125,3 +125,194 @@ for ints in multi_cases @test all(isapprox.(e, e_double_fun_ref; atol=1e-5)) end + +# 11. PowerLawPotential: analytic treatment of V(r) = v0*r^p +# The analytic path must reproduce the numerical CentralPotential path, and for the +# harmonic oscillator (p=2) also the exact spectrum. + +# 11a. constructor: p <= -3 must be rejected (radial integral diverges) +# 11a. validity of p is checked in sanity_checks (ISGL requires p > -3), not in the +# constructor, since the bound depends on dimension/angular momentum. ISGL_solve +# returns nothing instead of throwing (as for the other sanity-check failures). +@test PowerLawPotential(1.0, -3.0).p == -3.0 # construction itself is allowed +@test isnothing(ISGL_solve(make_phys_params3B3D(;masses=[m,m,m],species=[:x,:y,:z],interactions=[[PowerLawPotential(1.0,-3.0)],[],[]]),num_params)) +@test isnothing(ISGL_solve(make_phys_params3B3D(;masses=[m,m,m],species=[:x,:y,:z],interactions=[[],[PowerLawPotential(1.0,-4.0)],[]]),num_params)) +@test PowerLawPotential(-1.0, -1).p == -1.0 # integer exponents are accepted + +# 11b. harmonic oscillator (p=2) against the exact spectrum +vho_pow = PowerLawPotential(1/a*(m)/2*omega^2, 2.0) +pp_pow = make_phys_params3B3D(;masses=[m,m,m], species=[:x,:y,:z], interactions=[[vho_pow],[vho_pow],[vho_pow]]) +e_pow = ISGL_solve(pp_pow,num_params) /omega .- a; +@test all(isapprox.(e_pow[1:12], exact_results[1:12]; atol=1e-2)) + +# 11c. analytic vs numerical path for the very same potential (vcent_ho == vho_pow) +e_cent_ho = ISGL_solve(phys_params,num_params) +e_pow_ho = ISGL_solve(pp_pow,num_params) +@test all(isapprox.(e_cent_ho[1:12], e_pow_ho[1:12]; rtol=1e-5)) # the interpolated path is only good to ~3e-6 here + +# 11d. also with higher partial waves (exercises Lsum>0 in the shoulder solve) +e_pow_22 = ISGL_solve(pp_pow,num_params22) /omega .- a; +@test all(isapprox.(e_pow_22[1:12], exact_results[1:12]; atol=1e-3)) + +# 11e. Coulomb (p=-1): analytic vs numerical, on a Ps^- -like system +vatt_pow = PowerLawPotential(-1.0,-1.0); vrep_pow = PowerLawPotential(1.0,-1.0) +vatt_cent(r) = -1.0/r; vrep_cent(r) = 1.0/r +gp_coul = (;nmax=8,Nmax=8,r1=0.5,rnmax=25.0,R1=0.5,RNmax=25.0) +np_coul = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_coul) +pp_coul_pow = make_phys_params3B3D(;masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[vrep_pow],[vatt_pow],[vatt_pow]]) +pp_coul_cent = make_phys_params3B3D(;masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[vrep_cent],[vatt_cent],[vatt_cent]]) +e_coul_pow = ISGL_solve(pp_coul_pow, np_coul) +e_coul_cent = ISGL_solve(pp_coul_cent, np_coul) +@test all(isapprox.(e_coul_pow[1:3], e_coul_cent[1:3]; atol=1e-6)) +@test e_coul_pow[1] < -0.25 # Ps^- is bound below the Ps threshold + +# 11f. non-integer exponent (exercises the general Gamma-function path) +vpl_pow = PowerLawPotential(0.7, 1.5) +vpl_cent(r) = 0.7*r^1.5 +pp_pl_pow = make_phys_params3B3D(;masses=[m,m,m], species=[:x,:y,:z], interactions=[[vpl_pow],[vpl_pow],[vpl_pow]]) +pp_pl_cent = make_phys_params3B3D(;masses=[m,m,m], species=[:x,:y,:z], interactions=[[vpl_cent],[vpl_cent],[vpl_cent]]) +@test all(isapprox.(ISGL_solve(pp_pl_pow,num_params)[1:5], ISGL_solve(pp_pl_cent,num_params)[1:5]; rtol=1e-5)) + +# 11g. mixed interaction lists: checks index bookkeeping (powopt_arr padding, nint_arr) +ints_mix_pow = [[vga,vho_pow],[vho_pow],[vho_pow]] +ints_mix_cent = [[vga,vcent_ho],[vcent_ho],[vcent_ho]] +pp_mix_pow = make_phys_params3B3D(;masses=[m,m,m], species=[:x,:y,:z], interactions=ints_mix_pow) +pp_mix_cent = make_phys_params3B3D(;masses=[m,m,m], species=[:x,:y,:z], interactions=ints_mix_cent) +e_mix_pow = ISGL_solve(pp_mix_pow,num_params) +e_mix_cent = ISGL_solve(pp_mix_cent,num_params) +@test all(isfinite.(e_mix_pow[1:5])) +@test all(isapprox.(e_mix_pow[1:5], e_mix_cent[1:5]; rtol=1e-5)) + +# 11h. complex scaling: theta=0 must reproduce the non-scaled result +np_pow_00csm = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp, theta_csm = 0.0) +@test all(isapprox.(ISGL_solve(pp_pow,np_pow_00csm,complex_scaling=false)[1:5], + ISGL_solve(pp_pow,np_pow_00csm,complex_scaling=true)[1:5]; atol=1e-3)) + +# 11i. complex scaling at finite angle: analytic path must match the numerical one. +# This is the actual check of the csmfac^(-p) factor used for power-law potentials. +np_pow_10csm = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp, theta_csm = 10.0) +e_csm_pow = ISGL_solve(pp_pow, np_pow_10csm, complex_scaling=true) +e_csm_cent = ISGL_solve(phys_params,np_pow_10csm, complex_scaling=true) +@test all(isapprox.(e_csm_pow[1:5], e_csm_cent[1:5]; rtol=1e-3)) + +# 11j. PowerLawPotential as a central observable (uses the analytic precompute_varr! method) +obs_pow = (;stateindices=[1], centobs_arr=[PotentialFunction[PowerLawPotential(1.0,2.0)],PotentialFunction[],PotentialFunction[]], R2_arr=[1,0,0]) +obs_cent = (;stateindices=[1], centobs_arr=[PotentialFunction[CentralPotential(r->r^2)],PotentialFunction[],PotentialFunction[]], R2_arr=[1,0,0]) +_,_,co_pow,_ = ISGL_solve(phys_params, num_params; return_wavefunctions=true, observ_params=obs_pow) +_,_,co_cent,_ = ISGL_solve(phys_params, num_params; return_wavefunctions=true, observ_params=obs_cent) +@test isfinite(co_pow[1,1,1]) +@test isapprox(co_pow[1,1,1], co_cent[1,1,1]; rtol=1e-6) + + +# 12. Complex-ranged basis functions: nu -> nu*(1 + i*omega), selectable per Jacobi coordinate. +# Supported only for the analytically treated potentials (GaussianPotential, PowerLawPotential), +# since the interpolated path for a generic central potential needs a real etaprc. + +# 12a. option parsing +@test FewBodyToolkit.ISGL.parse_complex_ranged(:none) == (false,false) +@test FewBodyToolkit.ISGL.parse_complex_ranged(:r) == (true,false) +@test FewBodyToolkit.ISGL.parse_complex_ranged(:R) == (false,true) +@test FewBodyToolkit.ISGL.parse_complex_ranged(:both) == (true,true) +@test FewBodyToolkit.ISGL.parse_complex_ranged(true) == (true,true) # Bool alias, as in GEM2B_solve +@test FewBodyToolkit.ISGL.parse_complex_ranged(false) == (false,false) +@test_throws ErrorException FewBodyToolkit.ISGL.parse_complex_ranged(:rR) + +# 12b. unsupported potentials and observables are rejected (sanity_checks returns nothing) +# Unit masses (not the m=1/40 of the tests above): with m=1/40 this Gaussian binds nothing, so the +# lowest eigenvalue is a discretised continuum state and cannot be compared across bases in 12f. +pp_cr_gauss = make_phys_params3B3D(;masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[vga],[vga],[vga]]) +pp_cr_cent = make_phys_params3B3D(;masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[vg],[vg],[vg]]) +gp_cr = (;nmax=6,Nmax=6,r1=0.5,rnmax=8.0,R1=0.5,RNmax=7.0) +np_cr = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_cr,omega_cr=0.9) + +@test isnothing(ISGL_solve(pp_cr_cent,np_cr;complex_ranged=:both)) # plain function -> CentralPotential +@test isnothing(ISGL_solve(pp_cr_cent,np_cr;complex_ranged=:r)) +obs_r2 = (;stateindices=[1], centobs_arr=[PotentialFunction[],PotentialFunction[],PotentialFunction[]], R2_arr=[1,0,0]) +@test isnothing(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:r,return_wavefunctions=true,observ_params=obs_r2)) +# ... while the analytic potentials are accepted: +@test !isnothing(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:r)) + +# 12c. basis-size doubling per coordinate +n_none = length(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:none)) +n_r = length(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:r)) +n_R = length(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:R)) +n_both = length(ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:both)) +@test n_r == 2*n_none +@test n_R == 2*n_none +@test n_both == 4*n_none + +# 12d. energies stay real for complex ranges alone (H and S are hermitian, not just symmetric) +e_cr_r = ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:r) +@test eltype(e_cr_r) <: Real +@test all(isfinite.(e_cr_r[1:3])) + +# 12e. omega -> 0 regression. With omega=0 the conjugate half of the basis duplicates the original +# half exactly, so S is rank-deficient by construction and the redundant directions are removed by +# the threshold in cutSmallEV. The surviving spectrum must be the real-basis one. This is the +# sharpest check of the conjugation of the bra ranges: getting it wrong changes these numbers. +np_cr0 = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_cr,omega_cr=0.0) +e_ref0 = ISGL_solve(pp_cr_gauss,np_cr0;complex_ranged=:none) +for mode in (:r,:R,:both) + e0 = ISGL_solve(pp_cr_gauss,np_cr0;complex_ranged=mode) + @test all(isapprox.(e0[1:3], e_ref0[1:3]; atol=1e-6)) +end + +# the same regression with l,L > 0. This is the check on the branch cuts: for l>0 the norm carries +# the fractional power (2*nu)^((2l+3)/4) and the prefactors carry (pi/zeta)^(3/2)*(pi/etapr)^(3/2), +# none of which is exercised at l=L=0. A smaller basis keeps the (4x doubled) :both case affordable. +gp_cr1 = (;nmax=4,Nmax=4,r1=0.5,rnmax=8.0,R1=0.5,RNmax=7.0) +np_cr1_0 = make_num_params3B3D(;lmax=1,Lmax=1,gem_params=gp_cr1,omega_cr=0.0) +e_ref1 = ISGL_solve(pp_cr_gauss,np_cr1_0;complex_ranged=:none) +for mode in (:r,:R,:both) + e1 = ISGL_solve(pp_cr_gauss,np_cr1_0;complex_ranged=mode) + @test all(isapprox.(e1[1:3], e_ref1[1:3]; atol=1e-6)) +end + +# 12f. finite omega: the complex-ranged basis is a different span of twice the size, not a superset +# of the real one (the real Gaussian is not in it), so no variational inequality is claimed here. +# Both must locate the same bound state. The tolerance is set by the basis, not by the method: with +# this gp_cr the two bases differ by ~2.5e-3, shrinking to ~8e-4 at nmax=Nmax=8 and ~2e-4 at 10. +e_ref = ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=:none) +@test e_ref[1] < 0 # a genuine bound state, so that the comparison below is meaningful +for mode in (:r,:R,:both) + ecr = ISGL_solve(pp_cr_gauss,np_cr;complex_ranged=mode) + @test isapprox(ecr[1], e_ref[1]; atol=3e-3) # same state, not a spurious one +end + +# 12g. power-law interactions with complex ranges (Ps^- -like Coulomb system) +np_coul_cr = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_coul,omega_cr=0.9) +e_coul_real = ISGL_solve(pp_coul_pow, np_coul_cr; complex_ranged=:none) +e_coul_cr = ISGL_solve(pp_coul_pow, np_coul_cr; complex_ranged=:r) +@test isapprox(e_coul_cr[1], e_coul_real[1]; atol=1e-3) +@test e_coul_cr[1] < -0.25 # still bound below the Ps threshold + +# 12h. complex ranges together with complex scaling. Neither symmetry survives, so the full matrix +# is filled; at theta=0 this must still reproduce the complex-ranged result. +np_cr_00csm = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_cr,omega_cr=0.9,theta_csm=0.0) +e_cr_csm0 = ISGL_solve(pp_cr_gauss,np_cr_00csm;complex_ranged=:r,complex_scaling=true) +@test all(isapprox.(real.(e_cr_csm0[1:3]), e_cr_r[1:3]; atol=1e-6)) +@test all(abs.(imag.(e_cr_csm0[1:3])) .< 1e-8) + +# 12i. reduction to the two-body problem: switch off two of the three interactions and give the +# spectator coordinate a single, very broad Gaussian, so that -> 0. The r-part of the basis +# then mirrors the GEM2B basis exactly, and complex_ranged=:r must reproduce GEM2B's complex-ranged +# result. Equal masses are used so that the reduced mass of the interacting pair is unambiguous. +# Complex ranges are applied to r only: with Nmax=1 the conjugate partner of a single very broad +# range is nearly linearly dependent on it, which is exactly the situation per-coordinate control +# is meant to avoid. +gp_spec = (;nmax=8,Nmax=1,r1=1.0,rnmax=10.0,R1=1.0e4,RNmax=1.0e4) +np_spec = make_num_params3B3D(;lmax=0,Lmax=0,gem_params=gp_spec,omega_cr=0.9) +pp_spec = make_phys_params3B3D(;masses=[1.0,1.0,1.0], species=[:x,:y,:z], interactions=[[vga],[],[]]) + +pp2B_spec = make_phys_params2B(;mur=0.5, interactions=[vga], dim=3) +np2B_spec = make_num_params2B(;gem_params=(;nmax=8,r1=1.0,rnmax=10.0),omega_cr=0.9) + +e2_real = GEM2B_solve(pp2B_spec,np2B_spec) +e2_cr = GEM2B_solve(pp2B_spec,np2B_spec;complex_ranged=true) +e3_real = ISGL_solve(pp_spec,np_spec;complex_ranged=:none) +e3_cr = ISGL_solve(pp_spec,np_spec;complex_ranged=:r) + +nb2 = count(<(0), real.(e2_real)) # compare the bound states only +@test nb2 >= 1 +@test all(isapprox.(e3_real[1:nb2], e2_real[1:nb2]; atol=1e-3)) +@test all(isapprox.(e3_cr[1:nb2], e2_cr[1:nb2]; atol=1e-3))