From 5248065bbf80f41a2ff6fffcab8e0834f4a049b8 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 13:59:17 +0100 Subject: [PATCH 01/10] minor refactoring of coordinates --- Examples/HydrogenAnion.jl | 7 +- Manifest.toml | 66 +++++++++----- Project.toml | 1 + src/FewBodyPhysics.jl | 4 +- src/constants.jl | 27 ++++-- src/coordinates.jl | 177 ++++++++++++++++++++++---------------- src/matrix_elements.jl | 2 +- test/runtests.jl | 54 ++++-------- 8 files changed, 193 insertions(+), 145 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 74f41cd..aebd2e4 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -1,10 +1,11 @@ -using FewBodyPhysics +using Revise +using .FewBodyPhysics -masses = [Inf, 1, 1] +masses = [1e10, 1, 1] w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] K = [0 0 0; 0 1/2 0; 0 0 1/2] -J, U = Ω(masses) +J, U = jacobi_transform(masses) K_transformedformed = J * K * J' w_transformedformedformed = [U' * w_list[i] for i in 1:length(w_list)] Theortical_value = -0.527751016523 diff --git a/Manifest.toml b/Manifest.toml index 0972ed6..1c329cb 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "c019d638554320f1f902608dfe9771796da0203c" +project_hash = "999ca58762a5dcddc22a862313696fe8b0c1610a" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -16,9 +16,9 @@ version = "0.4.5" [[deps.Adapt]] deps = ["LinearAlgebra", "Requires"] -git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099" +git-tree-sha1 = "d80af0733c99ea80575f612813fa6aa71022d33a" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "4.0.4" +version = "4.1.0" [deps.Adapt.extensions] AdaptStaticArraysExt = "StaticArrays" @@ -83,6 +83,12 @@ git-tree-sha1 = "009060c9a6168704143100f36ab08f06c2af4642" uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a" version = "1.18.2+1" +[[deps.CodeTracking]] +deps = ["InteractiveUtils", "UUIDs"] +git-tree-sha1 = "7eee164f122511d3e4e1ebadb7956939ea7e1c77" +uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" +version = "1.3.6" + [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] git-tree-sha1 = "bce6804e5e6044c6daab27bb533d1295e4a2e759" @@ -91,9 +97,9 @@ version = "0.7.6" [[deps.ColorSchemes]] deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "PrecompileTools", "Random"] -git-tree-sha1 = "b5278586822443594ff615963b0c09755771b3e0" +git-tree-sha1 = "13951eb68769ad1cd460cdb2e64e5e95f1bf123d" uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4" -version = "3.26.0" +version = "3.27.0" [[deps.ColorTypes]] deps = ["FixedPointNumbers", "Random"] @@ -409,9 +415,9 @@ version = "0.10.2+0" [[deps.HTTP]] deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"] -git-tree-sha1 = "d1d712be3164d61d1fb98e7ce9bcbc6cc06b45ed" +git-tree-sha1 = "bc3f416a965ae61968c20d0ad867556367f2817d" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "1.10.8" +version = "1.10.9" [[deps.HarfBuzz_jll]] deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll"] @@ -458,6 +464,12 @@ git-tree-sha1 = "25ee0be4d43d0269027024d75a24c24d6c6e590c" uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8" version = "3.0.4+0" +[[deps.JuliaInterpreter]] +deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"] +git-tree-sha1 = "2984284a8abcfcc4784d95a9e2ea4e352dd8ede7" +uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a" +version = "0.9.36" + [[deps.LAME_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] git-tree-sha1 = "170b660facf5df5de098d866564877e119141cbd" @@ -543,9 +555,9 @@ version = "3.2.2+1" [[deps.Libgcrypt_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"] -git-tree-sha1 = "9fd170c4bbfd8b935fdc5f8b7aa33532c991a673" +git-tree-sha1 = "8be878062e0ffa2c3f67bb58a595375eda5de80b" uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4" -version = "1.8.11+0" +version = "1.11.0+0" [[deps.Libglvnd_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"] @@ -555,15 +567,15 @@ version = "1.6.0+0" [[deps.Libgpg_error_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "fbb1f2bef882392312feb1ede3615ddc1e9b99ed" +git-tree-sha1 = "c6ce1e19f3aec9b59186bdf06cdf3c4fc5f5f3e6" uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8" -version = "1.49.0+0" +version = "1.50.0+0" [[deps.Libiconv_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175" +git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269" uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531" -version = "1.17.0+0" +version = "1.17.0+1" [[deps.Libmount_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] @@ -614,9 +626,15 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[deps.LoggingExtras]] deps = ["Dates", "Logging"] -git-tree-sha1 = "c1dd6d7978c12545b4179fb6153b9250c96b0075" +git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "1.0.3" +version = "1.1.0" + +[[deps.LoweredCodeUtils]] +deps = ["JuliaInterpreter"] +git-tree-sha1 = "260dc274c1bc2cb839e758588c63d9c8b5e639d1" +uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b" +version = "3.0.5" [[deps.MacroTools]] deps = ["Markdown", "Random"] @@ -783,9 +801,9 @@ version = "1.10.0" [[deps.PlotThemes]] deps = ["PlotUtils", "Statistics"] -git-tree-sha1 = "6e55c6841ce3411ccb3457ee52fc48cb698d6fb0" +git-tree-sha1 = "41031ef3a1be6f5bbbf3e8073f210556daeae5ca" uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a" -version = "3.2.0" +version = "3.3.0" [[deps.PlotUtils]] deps = ["ColorSchemes", "Colors", "Dates", "PrecompileTools", "Printf", "Random", "Reexport", "StableRNGs", "Statistics"] @@ -902,6 +920,12 @@ git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" +[[deps.Revise]] +deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "REPL", "Requires", "UUIDs", "Unicode"] +git-tree-sha1 = "7f4228017b83c66bd6aa4fddeb170ce487e53bc7" +uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" +version = "3.6.2" + [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" @@ -1090,9 +1114,9 @@ version = "1.31.0+0" [[deps.XML2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"] -git-tree-sha1 = "1165b0443d0eca63ac1e32b8c0eb69ed2f4f8127" +git-tree-sha1 = "6a451c6f33a176150f315726eba8b92fbfdb9ae7" uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a" -version = "2.13.3+0" +version = "2.13.4+0" [[deps.XSLT_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "XML2_jll", "Zlib_jll"] @@ -1102,9 +1126,9 @@ version = "1.1.41+0" [[deps.XZ_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] -git-tree-sha1 = "ac88fb95ae6447c8dda6a5503f3bafd496ae8632" +git-tree-sha1 = "15e637a697345f6743674f1322beefbc5dcd5cfc" uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800" -version = "5.4.6+0" +version = "5.6.3+0" [[deps.Xorg_libICE_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl"] diff --git a/Project.toml b/Project.toml index 9f96763..b3089b3 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Optim = "429524aa-4258-5aef-a3af-852621145aeb" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" [compat] Documenter = "1.1.2" diff --git a/src/FewBodyPhysics.jl b/src/FewBodyPhysics.jl index 47db7ac..f7e3617 100644 --- a/src/FewBodyPhysics.jl +++ b/src/FewBodyPhysics.jl @@ -1,14 +1,12 @@ module FewBodyPhysics -using LinearAlgebra, Plots - include("coordinates.jl") include("matrix_elements.jl") include("sampling.jl") include("constants.jl") export m_π, m_π0, m_p, m_n, μ, m, ħc, mbare -export Ω, A_generate, transform_list, shift, w_gen, transform_coordinates, transform_back +export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters export corput, halton, run_simulation, run_simulation_nuclear diff --git a/src/constants.jl b/src/constants.jl index 79d9fa7..ae3ada8 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -1,9 +1,18 @@ -export m_π, m_π0, m_p, m_n, μ, m, ħc, mbare - -m_p = 938.27 -m_n = 939.57 -m_π0 = 134.98 -m_π = 139.57 -μ = (m_p+m_n)/2*m_π0/((m_p+m_n)/2+m_π0) -ħc = 197.3 -mbare=(m_p+m_n)/2 \ No newline at end of file +export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare + +# Physical Constants + +const m_p = 938.27 # Proton mass in MeV/c^2 +const m_n = 939.57 # Neutron mass in MeV/c^2 +const m_pi0 = 134.98 # Neutral pion mass (π⁰) in MeV/c^2 +const m_pi = 139.57 # Charged pion mass (π⁺/π⁻) in MeV/c^2 +const ħc = 197.3269804 # Planck's constant times speed of light in MeV·fm + +const m_bare = (m_p + m_n) / 2 # Average nucleon mass in MeV/c^2 + +""" + μ = (m_bare * m_pi0) / (m_bare + m_pi0) + +Reduced mass of the nucleon-pion system in MeV/c^2. +""" +const μ = (m_bare * m_pi0) / (m_bare + m_pi0) diff --git a/src/coordinates.jl b/src/coordinates.jl index 9fccac1..76d9315 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -1,104 +1,135 @@ using LinearAlgebra -export Ω, A_generate, transform_list, shift, w_gen, transform_coordinates, transform_back +export ParticleSystem, jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates + +struct ParticleSystem + masses::Vector{Float64} + J::Matrix{Float64} + U::Matrix{Float64} + + function ParticleSystem(masses::Vector{Float64}) + @assert length(masses) ≥ 2 "At least two masses are required for a particle system." + J, U = jacobi_transform(masses) + new(masses, J, U) + end +end + """ - Ω(masses::Array) + jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matrix{Float64}} -Calculate the Jacobi transformation matrix `J` and its inverse `U` for a system of particles with specified `masses`. +Compute the Jacobi transformation matrices `J` and `U` for a system of particles with specified masses. # Arguments -- `masses::Array`: A vector of masses for the particles. +- `masses::Vector{Float64}`: A vector of masses for the particles. # Returns -- `J::Matrix`: The Jacobi transformation matrix. -- `U::Matrix`: The inverse of the Jacobi transformation matrix. +- `(J::Matrix{Float64}, U::Matrix{Float64})`: The Jacobi transformation matrix and its pseudoinverse. # Notes -- For systems with more than one particle, the returned matrices exclude the last row/column for proper dimensionality in transformations. +- The matrices `J` and `U` are used to transform between particle coordinates and Jacobi coordinates. +- The pseudoinverse `U` is used instead of the inverse to handle cases where `J` is not square. """ -Ω(masses::Array) = begin - dim = size(masses, 1) - J = zeros(dim, dim) - for i in 1:dim - sum_m = sum(masses[1:i]) - for j in 1:dim - J[i, j] = j == i + 1 ? -1 : i + 1 < j ? 0 : masses[j] / sum_m - J[i, j] = isnan(J[i, j]) ? 1 : J[i, j] +function jacobi_transform(masses::Vector{Float64})::Tuple{Matrix{Float64}, Matrix{Float64}} + N = length(masses) + @assert N ≥ 2 "At least two masses are required for Jacobi transformation." + J = zeros(Float64, N - 1, N) + + for k in 1:N - 1 + mk = masses[k + 1] + Mk = sum(masses[1:k]) + μk = sqrt(mk * Mk / (mk + Mk)) + for j in 1:N + if j ≤ k + J[k, j] = μk * masses[j] / Mk + elseif j == k + 1 + J[k, j] = -μk + else + J[k, j] = 0.0 + end end end - U = inv(J) - U = dim > 1 ? U[:, 1:end-1] : U - J = dim > 1 ? J[1:end-1, :] : J + + U = pinv(J) return J, U end + """ - A_generate(bij::Array, w_list::Array) + generate_A_matrix(bij::Vector{Float64}, w_list::Vector{Vector{Float64}})::Matrix{Float64} -Generate a matrix A for Gaussian basis functions given width parameters `bij` and weight vectors `w_list`. +Generate the matrix `A` for Gaussian basis functions given width parameters `bij` and weight vectors `w_list`. # Arguments -- `bij::Array`: A vector of width parameters for the Gaussian basis functions. -- `w_list::Array`: A list of weight vectors. +- `bij::Vector{Float64}`: A vector of width parameters for the Gaussian basis functions. +- `w_list::Vector{Vector{Float64}}`: A list of weight vectors. # Returns -- `Matrix`: The sum of weighted outer products of `w_list`, scaled by `bij`. +- `A::Matrix{Float64}`: The sum of weighted outer products of `w_list`, scaled by `bij`. # Notes -- This function is used to construct basis elements for the expansion of few-body wavefunctions. +- This function constructs the `A` matrix used in the correlated Gaussian method. """ -A_generate(bij, w_list::Array) = begin - dim = length(w_list) - mat_list = [(w_list[i] * w_list[i]') ./ (bij[i]^2) for i in 1:dim] - return sum(mat_list) +function generate_A_matrix(bij::Vector{Float64}, w_list::Vector{Vector{Float64}})::Matrix{Float64} + @assert length(bij) == length(w_list) "Length of `bij` and `w_list` must be equal." + dim = length(w_list[1]) + A = zeros(Float64, dim, dim) + for i in 1:length(bij) + w = w_list[i] + @assert length(w) == dim "All weight vectors must have the same dimension." + A += (w * w') / (bij[i]^2) + end + return A end """ - transform_list(α::Array) + transform_list(α::Vector{Float64})::Vector{Matrix{Float64}} Transform a list of scalar values `α` into a list of 1x1 matrices. # Arguments -- `α::Array`: A list of scalar values. +- `α::Vector{Float64}`: A list of scalar values. # Returns -- `Array`: A list of 1x1 matrices where each matrix contains one of the scalar values from `α`. +- `Array{Matrix{Float64}}`: A list of 1x1 matrices where each matrix contains one of the scalar values from `α`. """ -function transform_list(α::Array) - return[ones(1, 1) .* α[i] for i in 1:length(α)] +function transform_list(α::Vector{Float64})::Vector{Matrix{Float64}} + return [Matrix{Float64}([α_i]) for α_i in α] end + """ - shift(a::Array, b::Array, mat::Matrix=I) + shift_vectors(a::Matrix{Float64}, b::Matrix{Float64}, mat::Union{Nothing, Matrix{Float64}}=nothing)::Float64 Calculate the weighted sum of the element-wise product of vectors `a` and `b` using matrix `mat`. # Arguments -- `a::Array`: A vector or matrix. -- `b::Array`: A vector or matrix of the same size as `a`. -- `mat::Matrix`: An optional matrix to weight the product (default is the identity matrix). +- `a::Matrix{Float64}`: A matrix where each column is a vector `a_i`. +- `b::Matrix{Float64}`: A matrix where each column is a vector `b_j`. +- `mat::Union{Nothing, Matrix{Float64}}`: An optional matrix to weight the product (default is the identity matrix). # Returns -- `Float64`: The weighted sum of products. +- `sum_val::Float64`: The weighted sum of products. # Notes -- `a` and `b` are typically shift vectors in the configuration space of a few-body system. -- If `mat` is provided, its dimensions must match the number of elements in `a` and `b`. +- The matrices `a` and `b` should have the same dimensions. """ -shift(a, b, mat) = begin +function shift_vectors(a::Matrix{Float64}, b::Matrix{Float64}, mat::Union{Nothing, Matrix{Float64}}=nothing)::Float64 n = size(a, 2) + @assert n == size(b, 2) "Matrices `a` and `b` must have the same number of columns." + mat = mat === nothing ? I(n) : mat + @assert size(mat) == (n, n) "Matrix `mat` must be square with size equal to number of vectors." + sum_val = 0.0 - mat = isnothing(mat) ? I : mat - @assert n == size(mat, 1) "ERROR! Matrix shape does not match number of shift vectors." for i in 1:n for j in 1:n - sum_val += mat[i, j] * (a[:, i]' * b[:, j]) + sum_val += mat[i, j] * dot(view(a, :, i), view(b, :, j)) end end return sum_val end + """ - w_gen(dim::Int, i::Int, j::Int) + generate_weight_vector(dim::Int, i::Int, j::Int)::Vector{Int} -Generate a weight vector for the i-th and j-th coordinates in a space of dimension `dim`. +Generate a weight vector for the `i`-th and `j`-th coordinates in a space of dimension `dim`. # Arguments - `dim::Int`: The dimension of the space. @@ -106,47 +137,47 @@ Generate a weight vector for the i-th and j-th coordinates in a space of dimensi - `j::Int`: The index for the negative element in the weight vector. # Returns -- `Vector{Int}`: A vector with 1 at the i-th position, -1 at the j-th position, and 0 elsewhere. - -# Notes -- This function is useful for generating basis vectors in few-body coordinate transformations. +- `w::Vector{Int}`: A vector with 1 at the `i`-th position, -1 at the `j`-th position, and 0 elsewhere. """ -w_gen(dim, i, j) = dim == 1 ? [1] : [Int(k == i) - Int(k == j) for k in 1:dim] +function generate_weight_vector(dim::Int, i::Int, j::Int)::Vector{Int} + @assert 1 ≤ i ≤ dim "Index `i` must be between 1 and $dim." + @assert 1 ≤ j ≤ dim "Index `j` must be between 1 and $dim." + w = zeros(Int, dim) + w[i] = 1 + w[j] = -1 + return w +end + """ - transform_coordinates(Ω::Matrix{Float64}, r::Vector{Float64}) + transform_coordinates(J::Matrix{Float64}, r::Vector{Float64})::Vector{Float64} -Transform the coordinates `r` of a system using the Jacobi matrix `Ω`. +Transform the coordinates `r` of a system using the Jacobi matrix `J`. # Arguments -- `Ω::Matrix{Float64}`: The Jacobi transformation matrix. +- `J::Matrix{Float64}`: The Jacobi transformation matrix. - `r::Vector{Float64}`: The coordinates to be transformed. # Returns -- `Vector{Float64}`: The transformed coordinates. - -# Notes -- This function applies the inverse of Jacobi matrix `J` to the coordinate vector `r`. +- `x::Vector{Float64}`: The transformed coordinates. """ -function transform_coordinates(Ω::Matrix{Float64}, r::Vector{Float64}) - J, U = Ω(masses) - return J \ r +function transform_coordinates(J::Matrix{Float64}, r::Vector{Float64})::Vector{Float64} + @assert size(J, 2) == length(r) "Matrix `J` columns must match length of vector `r`." + return J * r end + """ - transform_back(Ω::Matrix{Float64}, x::Matrix{Float64}) + inverse_transform_coordinates(U::Matrix{Float64}, x::Vector{Float64})::Vector{Float64} -Transform the coordinates `x` back to the original system using the inverse of the Jacobi matrix `Ω`. +Transform the coordinates `x` back to the original system using the inverse of the Jacobi matrix. # Arguments -- `Ω::Matrix{Float64}`: The Jacobi transformation matrix. -- `x::Matrix{Float64}`: The coordinates to be transformed back. +- `U::Matrix{Float64}`: The inverse Jacobi transformation matrix. +- `x::Vector{Float64}`: The coordinates in Jacobi space. # Returns -- `Matrix{Float64}`: The coordinates transformed back to the original system. - -# Notes -- This function applies the inverse of matrix `U` to the coordinate matrix `x`. +- `r::Vector{Float64}`: The coordinates transformed back to the original system. """ -function transform_back(Ω::Matrix, x::Matrix, masses::Vector) - J, U = Ω(masses) - return U \ x -end \ No newline at end of file +function inverse_transform_coordinates(U::Matrix{Float64}, x::Vector{Float64})::Vector{Float64} + @assert size(U, 1) == length(x) "Matrix `U` rows must match length of vector `x`." + return U * x +end diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index c50dcc4..0a3a186 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -109,7 +109,7 @@ function S_energy(bij, K, w) α = [] dim = length(w) for i in 1:dim:length(bij) - A = A_generate(bij[i:i+dim-1], w) + A = generate_A_matrix(bij[i:i+dim-1], w) push!(α, A) end N, kinetic, Coulomb = S_wave(α, K, w) diff --git a/test/runtests.jl b/test/runtests.jl index c30a8de..93358c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,50 +3,34 @@ using Test using LinearAlgebra @testset "FewBodyPhysics.jl" begin - @testset "Ω tests" begin - J, U = Ω([1, 2, 3]) - @test J isa Matrix - @test U isa Matrix + @testset "Coordinate Transformations" begin + masses = [1.0, 2.0, 3.0] + J, U = JacobiTransform(masses) @test size(J) == (2, 3) @test size(U) == (3, 2) - end - @testset "A_generate tests" begin - A = A_generate([1, 2, 3], [[1, 2, 3], [4, 5, 6], [7, 8, 9]]) - @test A isa Matrix + bij = [0.5, 0.8] + w_list = [generate_weight_vector(3, 1, 2), generate_weight_vector(3, 2, 3)] + A = generate_A_matrix(bij, w_list) @test size(A) == (3, 3) - end - @testset "transform_list tests" begin - transformed = transform_list([1, 2, 3]) - @test transformed isa Array - @test all([isa(x, Matrix) for x in transformed]) - @test all([size(x) == (1, 1) for x in transformed]) - end + α = [1.0, 2.0, 3.0] + transformed_list = transform_list(α) + @test length(transformed_list) == 3 + @test transformed_list[1] == [1.0] - @testset "shift tests" begin - s = shift([1, 2, 3], [4, 5, 6]) - @test s isa Float64 - @test_throws AssertionError shift([1, 2, 3], [4, 5, 6], [1, 2]) - end + a = rand(3, 2) + b = rand(3, 2) + sum_val = shift_vectors(a, b) + @test isa(sum_val, Float64) - @testset "w_gen tests" begin - w = w_gen(3, 1, 2) - @test w isa Vector{Int} - @test size(w) == (3,) + w = generate_weight_vector(3, 1, 2) @test w == [1, -1, 0] - end - - @testset "transform_coordinates tests" begin - r_transformed = transform_coordinates([1 0; 0 1], [1, 2]) - @test r_transformed isa Vector{Float64} - @test size(r_transformed) == (2,) - end - @testset "transform_back tests" begin - x_transformed = transform_back([1 0; 0 1], [1 2; 3 4]) - @test x_transformed isa Matrix{Float64} - @test size(x_transformed) == (2, 2) + r = rand(3) + x = transform_coordinates(J, r) + r_back = inverse_transform_coordinates(U, x) + @test norm(r - r_back) < 1e-10 end @testset "corput tests" begin From d616558747f81a3348df9d25b391910679ab1e8b Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 14:14:05 +0100 Subject: [PATCH 02/10] minor refactoring of matrix_elements --- Examples/HydrogenAnion.jl | 2 +- Examples/PionNucleon.jl | 2 +- Examples/Positronium.jl | 7 +- src/FewBodyPhysics.jl | 4 +- src/matrix_elements.jl | 177 +++++++++++++++++----------------- src/sampling.jl | 198 +++++++++++--------------------------- 6 files changed, 156 insertions(+), 234 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index aebd2e4..24c3ba8 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -1,7 +1,7 @@ using Revise using .FewBodyPhysics -masses = [1e10, 1, 1] +masses = [1e15, 1.0, 1.0] w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] K = [0 0 0; 0 1/2 0; 0 0 1/2] diff --git a/Examples/PionNucleon.jl b/Examples/PionNucleon.jl index 3e36d70..31d2a60 100644 --- a/Examples/PionNucleon.jl +++ b/Examples/PionNucleon.jl @@ -4,7 +4,7 @@ b = 3.9 S = 41.5 params = [b, S] -masses = [(m_p+m_n)/2, m_π] +masses = [(m_p+m_n)/2, m_pi] energies, Gaussians, eigenvectors, coordinates, masses = run_simulation_nuclear(5,2,5,masses,params) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index d1724dc..8ca59ca 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -1,9 +1,10 @@ -using FewBodyPhysics +using Revise +using .FewBodyPhysics w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] -masses = [1,1,1] +masses = [1.0,1.0,1.0] K = [1/2 0 0; 0 1/2 0; 0 0 1/2] -J, U = Ω(masses) +J, U = jacobi_transform(masses) K_transformed = J * K * J' w_transformed = [U' * w_list[i] for i in 1:length(w_list)] diff --git a/src/FewBodyPhysics.jl b/src/FewBodyPhysics.jl index f7e3617..9511fd0 100644 --- a/src/FewBodyPhysics.jl +++ b/src/FewBodyPhysics.jl @@ -5,9 +5,9 @@ include("matrix_elements.jl") include("sampling.jl") include("constants.jl") -export m_π, m_π0, m_p, m_n, μ, m, ħc, mbare +export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates -export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters +export S_elements, S_wave, S_energy, P_elements, pion_nucleon,compute_eigen_system, get_minimum_energy, optimize_global_parameters export corput, halton, run_simulation, run_simulation_nuclear end diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index 0a3a186..fe75cbd 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -1,165 +1,168 @@ -using Optim, LinearAlgebra +using LinearAlgebra +using Optim -export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters +export S_elements, S_wave, S_energy, P_elements, pion_nucleon, + compute_eigen_system, get_minimum_energy, optimize_global_parameters struct PositiveDefiniteSymmetricMatrix{T<:Real} matrix::Matrix{T} +end - function PositiveDefiniteSymmetricMatrix(mat::Matrix{T}) where T <: Real - if size(mat, 1) != size(mat, 2) - throw(ArgumentError("The matrix must be square.")) - end - if !issymmetric(mat) - throw(ArgumentError("The matrix must be symmetric.")) - end - if !isposdef(mat) - throw(ArgumentError("The matrix must be positive definite.")) - end - return new{T}(mat) - end +function PositiveDefiniteSymmetricMatrix(mat::Matrix{T}) where T <: Real + @assert size(mat, 1) == size(mat, 2) "Matrix must be square." + @assert issymmetric(mat) "Matrix must be symmetric." + @assert isposdef(mat) "Matrix must be positive definite." + return PositiveDefiniteSymmetricMatrix{T}(mat) end + """ - S_elements(A, B, K, w=nothing) + S_elements(A::Matrix{Float64}, B::Matrix{Float64}, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) Calculate matrix elements for overlap, kinetic energy, and optionally the Coulomb term. # Arguments -- `A::Matrix`: Matrix representing the width of Gaussian basis functions for state `i`. -- `B::Matrix`: Matrix representing the width of Gaussian basis functions for state `j`. -- `K::Matrix`: Kinetic energy matrix. -- `w::Vector` (optional): Weight vectors for the particles involved. +- `A`: Width matrix of Gaussian basis functions for state `i`. +- `B`: Width matrix of Gaussian basis functions for state `j`. +- `K`: Kinetic energy matrix. +- `w` (optional): List of weight vectors for the particles involved. # Returns -- `M0::Float64`: The overlap matrix element between the two states. -- `tra::Float64`: The trace used in the kinetic energy calculation. -- `Coulomb_term::Float64` (optional): The Coulomb interaction term, if weight vectors `w` are provided. - -# Notes -- The Coulomb term is calculated only if the weight vectors `w` are specified. +- `M0`: The overlap matrix element between the two states. +- `trace`: The trace used in the kinetic energy calculation. +- `Coulomb_term` (optional): The Coulomb interaction term, if weight vectors `w` are provided. """ -function S_elements(A::Matrix, B::Matrix, K::Matrix, w=nothing) +function S_elements(A::Matrix{Float64}, B::Matrix{Float64}, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) dim = size(A, 1) - Coulomb_term = 0.0 D = A + B + @assert isposdef(D) "Matrix D must be positive definite." R = inv(D) - M0 = (π^dim / det(D))^(3.0 / 2) + detD = det(D) + @assert detD > 0 "Determinant of D must be positive." + M0 = (π^dim / detD)^(1.5) trace = tr(B * K * A * R) - if !isnothing(w) + Coulomb_term = 0.0 + if w !== nothing for k in 1:length(w) - β = 1 ./((w[k])' * R * w[k]) - Coulomb_term += (k == 3 ? 2 : -2) .* sqrt(β / π) .* M0 + wk = w[k] + β = 1.0 / (wk' * R * wk) + factor = ifelse(k == 3, 2.0, -2.0) + Coulomb_term += factor * sqrt(β / π) * M0 end return M0, trace, Coulomb_term else return M0, trace end end + """ - S_wave(α, K, w=nothing) + S_wave(α::Vector{Matrix{Float64}}, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) -Calculate the wavefunction overlap, kinetic energy, and optionally Coulomb interaction matrices for a given set of basis functions. +Calculate the overlap, kinetic energy, and optionally Coulomb interaction matrices for a given set of basis functions. # Arguments -- `α::Vector`: A list of scalar width parameters for the Gaussian basis functions. -- `K::Matrix`: Kinetic energy matrix. -- `w::Vector` (optional): Weight vectors for the particles involved. +- `α`: A list of width matrices for the Gaussian basis functions. +- `K`: Kinetic energy matrix. +- `w` (optional): List of weight vectors for the particles involved. # Returns -- `overlap::Matrix`: The overlap matrix for the basis functions. -- `kinetic::Matrix`: The kinetic energy matrix for the basis functions. -- `Coulomb::Matrix` (optional): The Coulomb interaction matrix, if weight vectors `w` are specified. - -# Notes -- The Coulomb matrix is computed only if the weight vectors `w` are specified. +- `overlap`: The overlap matrix for the basis functions. +- `kinetic`: The kinetic energy matrix for the basis functions. +- `Coulomb`: The Coulomb interaction matrix, if weight vectors `w` are specified. """ -function S_wave(α, K, w=nothing) +function S_wave(α::Vector{Matrix{Float64}}, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) len = length(α) - α = transform_list(α) - overlap = zeros(len, len) - kinetic = zeros(len, len) - Coulomb = zeros(len, len) + overlap = zeros(Float64, len, len) + kinetic = zeros(Float64, len, len) + Coulomb = w !== nothing ? zeros(Float64, len, len) : nothing for i in 1:len for j in 1:i A, B = α[i], α[j] - M0, trace, Coulomb_term = S_elements(A, B, K, w) + if w !== nothing + M0, trace, Coulomb_term = S_elements(A, B, K, w) + Coulomb[i, j] = Coulomb[j, i] = Coulomb_term + else + M0, trace = S_elements(A, B, K) + end overlap[i, j] = overlap[j, i] = M0 kinetic[i, j] = kinetic[j, i] = 6 * trace * M0 - Coulomb[i, j] = Coulomb[j, i] = Coulomb_term end end return overlap, kinetic, Coulomb end + """ - S_energy(bij, K, w) + S_energy(bij::Vector{Float64}, K::Matrix{Float64}, w::Vector{Vector{Float64}}) Compute the ground state energy of the system using the basis functions specified by the width parameters `bij`. # Arguments -- `bij::Vector`: A list of width parameters for the Gaussian basis functions. -- `K::Matrix`: Kinetic energy matrix. -- `w::Vector`: Weight vectors for the particles involved. +- `bij`: A vector of width parameters for the Gaussian basis functions. +- `K`: Kinetic energy matrix. +- `w`: List of weight vectors for the particles involved. # Returns -- `E0::Float64`: The lowest eigenvalue computed from the Hamiltonian, considered as the ground state energy of the system. - -# Notes -- This function constructs the Hamiltonian from the overlap, kinetic, and Coulomb matrices and solves for its eigenvalues. +- `E0`: The lowest eigenvalue computed from the Hamiltonian. """ -function S_energy(bij, K, w) - α = [] +function S_energy(bij::Vector{Float64}, K::Matrix{Float64}, w::Vector{Vector{Float64}}) + α = Vector{Matrix{Float64}}() dim = length(w) + @assert length(bij) % dim == 0 "Length of bij must be a multiple of the dimension." for i in 1:dim:length(bij) - A = generate_A_matrix(bij[i:i+dim-1], w) + bij_segment = bij[i:i+dim-1] + A = generate_A_matrix(bij_segment, w) push!(α, A) end - N, kinetic, Coulomb = S_wave(α, K, w) + overlap, kinetic, Coulomb = S_wave(α, K, w) H = kinetic + Coulomb - E,v = eigen(H, N) - E0 = minimum(E) + E, _ = eigen(H, overlap) + E0 = minimum(real(E)) return E0 end + + """ - P_elements(a, b, A, B, K, w=nothing) + P_elements(a::Vector{Float64}, b::Vector{Float64}, A::PositiveDefiniteSymmetricMatrix, B::PositiveDefiniteSymmetricMatrix, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) Calculate the perturbation matrix elements given two basis states represented by vectors `a` and `b`, and their respective width matrices `A` and `B`. # Arguments -- `a::Vector`: The coefficient vector for basis state `i`. -- `b::Vector`: The coefficient vector for basis state `j`. -- `A::Matrix`: Matrix representing the width of Gaussian basis functions for state `i`. -- `B::Matrix`: Matrix representing the width of Gaussian basis functions for state `j`. -- `K::Matrix`: Kinetic energy matrix. -- `w::Vector` (optional): Weight vectors for the particles involved. +- `a`: The coefficient vector for basis state `i`. +- `b`: The coefficient vector for basis state `j`. +- `A`: Width matrix for state `i`. +- `B`: Width matrix for state `j`. +- `K`: Kinetic energy matrix. +- `w` (optional): List of weight vectors for the particles involved. # Returns -- `M1::Float64`: The overlap perturbation term. -- `kinetic::Float64`: The kinetic energy perturbation term. -- `Coulomb_term::Float64` (optional): The Coulomb interaction perturbation term, if weight vectors `w` are provided. - -# Notes -- The Coulomb interaction perturbation term is calculated only if the weight vectors `w` are specified. +- `M1`: The overlap perturbation term. +- `kinetic`: The kinetic energy perturbation term. +- `Coulomb_term` (optional): The Coulomb interaction perturbation term, if weight vectors `w` are provided. """ -function P_elements(a, b, A::PositiveDefiniteSymmetricMatrix, B, K, w=nothing) - D = A + B +function P_elements(a::Vector{Float64}, b::Vector{Float64}, A::PositiveDefiniteSymmetricMatrix, B::PositiveDefiniteSymmetricMatrix, K::Matrix{Float64}, w::Union{Nothing, Vector{Vector{Float64}}}=nothing) + D = A.matrix + B.matrix + @assert isposdef(D) "Matrix D must be positive definite." R = inv(D) - M0, trace = S_elements(A, B, K) - M1 = 1/2 * (a' * R * b) * M0 # Overlap - kinetic = 6 * trace * M1 # Kinetic matrix elements + M0, trace = S_elements(A.matrix, B.matrix, K) + M1 = 0.5 * (a' * R * b) * M0 # Overlap perturbation term + + # Kinetic energy perturbation term + kinetic = 6 * trace * M1 kinetic += (a' * K * b) * M0 - kinetic += -(a' * (K * A * R) * b) * M0 - kinetic += -(a' * (R * B * K) * b) * M0 - kinetic += (a' * (R * B * K * A * R) * b) * M0 - kinetic += (a' * (R * B * K * A * R) * b) * M0 + kinetic -= (a' * K * A.matrix * R * b) * M0 + kinetic -= (a' * R * B.matrix * K * b) * M0 + kinetic += (a' * R * B.matrix * K * A.matrix * R * b) * M0 if w !== nothing - β = 1 ./ ((w' * R * w)[1]) # Coulomb terms + w_concat = hcat(w...) + β = 1.0 / (w_concat' * R * w_concat)[1] Coulomb_term = 2 * sqrt(β / π) * M1 - Coulomb_term += -sqrt(β / π) * β / 3 * (a' * (R * w * w' * R) * b) * M0 + Coulomb_term -= sqrt(β / π) * β / 3 * (a' * R * w_concat * w_concat' * R * b) * M0 return M1, kinetic, Coulomb_term else return M1, kinetic end end + """ pion_nucleon(alphas, masses, params) diff --git a/src/sampling.jl b/src/sampling.jl index 1bd985e..9f8b470 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -2,17 +2,8 @@ using Plots export corput, halton, run_simulation, run_simulation_nuclear -""" -Generate the nth element of the van der Corput sequence in base b. - -# Arguments -- `n`: The nth element of the sequence to be generated. -- `b`: The base for the van der Corput sequence. Default is 3. - -# Returns -- `q`: The nth element of the van der Corput sequence in base b. -""" -corput(n, b=3) = begin +# Generate the nth element of the van der Corput sequence in base b +function corput(n, b=3) q, bk = 0.0, 1 / b while n > 0 n, rem = divrem(n, b) @@ -21,168 +12,95 @@ corput(n, b=3) = begin end return q end -""" -Generate the nth d-dimensional point in the Halton sequence. -# Arguments -- `n`: The nth element of the sequence to be generated. -- `d`: The dimension of the space. - -# Returns -- An array containing the nth d-dimensional point in the Halton sequence. - -# Errors -- Throws an assertion error if `d` exceeds the number of basis elements. -""" -halton(n, d) = begin - base = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,233,239,241,251,257,263,269,271,277,281] - @assert length(base) ≥ d "Error: d exceeds the number of basis elements." +# Generate the nth d-dimensional point in the Halton sequence +function halton(n, d) + base = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281] + @assert d <= length(base) "Dimension exceeds available bases." return [corput(n, base[i]) for i in 1:d] end -""" -Run the simulation for a quantum system using quasi-random or pseudo-random methods to determine the S-wave convergence. +# Helper function for calculating energies in quasi-random methods +function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method::Symbol) + E_list, bases, gaussians = Float64[], [], Int[] + E_low = Inf + for i in 1:num_gauss + hal = halton(i, 15 * length(w_transformed)) + bij = Float64.(-log.(hal) * b1) + best_energy, best_base = E_low, nothing + + for j in 1:length(w_transformed):length(hal) + base_segment = bij[j:j+length(w_transformed)-1] + E0 = S_energy(base_segment, K_transformed, w_transformed) + if E0 < best_energy + best_energy = E0 + best_base = copy(base_segment) + end + end -# Arguments -- `num_gauss::Int`: The number of Gaussians to use in the simulation. Default is 15. -- `method::Symbol`: The method to use for the simulation. Can be `:quasirandom`, `:quasirandomrefined`, or `:psudorandom`. Default is `:quasirandom`. -- `plot_result::Bool`: Whether to plot the results. Default is true. + E_low = min(E_low, best_energy) + push!(E_list, E_low) + push!(bases, best_base) + push!(gaussians, i) + println(E_low) + end + return E_list, bases, gaussians, E_low +end -# Returns -- `p`: The plot object if `plot_result` is true. +# Main simulation function +function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, plot_result::Bool=true) + b1 = 10.0 -# Notes -- The function prints various convergence information and, if `plot_result` is true, displays a plot of the numerical result against the theoretical value. -""" -function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Array, K_transformed::Array, plot_result::Bool=true) - b1 = 10 - E_list = [] - gaussians = [] - bij = [] - E_low = 1e10 - global bases = [] - global base_test = [] + # Choose appropriate method if method == :quasirandom println("---------QUASI-RANDOM METHOD---------") - - for i in 1:num_gauss - hal = halton(i, 15 * length(w_transformed)) - bij = -log.(hal) .* b1 - for j in 1:length(w_transformed):length(hal) - append!(base_test, bij[j:j+length(w_transformed)-1]) - E0 = S_energy(base_test, K_transformed, w_transformed) - if E0 <= E_low - E_low = E0 - global base_curr = copy(bij[j:j+length(w_transformed)-1]) - end - base_test = base_test[1:end-length(w_transformed)] - end - append!(bases, base_curr) - append!(base_test, base_curr) - push!(E_list, E_low) - println(E_low) - push!(gaussians, i) - end + E_list, bases, gaussians, E_low = calculate_energies(num_gauss, w_transformed, K_transformed, b1, method) elseif method == :quasirandomrefined println("---------QUASI-RANDOM METHOD WITH REFINEMENT---------") + E_list, bases, gaussians, E_low = calculate_energies(num_gauss, w_transformed, K_transformed, b1, method) - for i in 1:num_gauss - hal = halton(i, 15 * length(w_transformed)) - bij = -log.(hal) .* b1 - for j in 1:length(w_transformed):length(hal) - append!(base_test, bij[j:j+length(w_transformed)-1]) - E0 = S_energy(base_test, K_transformed, w_transformed) - if E0 <= E_low - E_low = E0 - global base_curr = copy(bij[j:j+length(w_transformed)-1]) - end - base_test = base_test[1:end-length(w_transformed)] - end - append!(bases, base_curr) - append!(base_test, base_curr) - push!(E_list, E_low) - println(E_low) - push!(gaussians, i) - end - - bases_ref = copy(bases) - E_ref = E_list[end] - E_list_ref = [] - for i in 1:length(bases_ref)-length(w_transformed) - rand_ref = rand(200 * length(w_transformed)) - bij_ref = .-log.(rand_ref) * b1 + # Additional refinement step + for i in 1:length(bases)-length(w_transformed) + rand_ref = Float64.(-log.(rand(200 * length(w_transformed))) * b1) for j in 1:length(w_transformed):length(rand_ref) - bases_ref[i:i+length(w_transformed)-1] = bij_ref[j:j+length(w_transformed)-1] - E_test = S_energy(bases_ref, K_transformed, w_transformed) - if E_test < E_ref - E_ref = E_test - bases[i:i+length(w_transformed)-1] = bij_ref[j:j+length(w_transformed)-1] + refined_base = bases[i][1:end-length(w_transformed)] # Adjust size accordingly + refined_base .= rand_ref[j:j+length(w_transformed)-1] + E_refined = S_energy(refined_base, K_transformed, w_transformed) + if E_refined < E_low + E_low = E_refined + bases[i] .= refined_base end end - bases_ref = copy(bases) - push!(E_list_ref, E_ref) - println("E_ref: ", E_ref) + println("Refined Energy: ", E_low) end - println("Energy after refinement: ", E_ref) - println("Difference in energy from before refinement: ", abs(E_ref - E_list[end])) - elseif method == :psudorandom println("---------PSEUDO-RANDOM METHOD---------") - for i in 1:num_gauss - rnd = rand(400 * length(w_transformed)) - bij2 = -log.(rnd) * b1 - for j in 1:length(w_transformed):length(rnd) - append!(base_test, bij2[j:j+length(w_transformed)-1]) - E0 = S_energy(base_test, K_transformed, w_transformed) - if E0 <= E_low - E_low = E0 - global base_curr = copy(bij2[j:j+length(w_transformed)-1]) - end - base_test = base_test[1:end-length(w_transformed)] - end - append!(bases, base_curr) - append!(base_test, base_curr) - push!(E_list, E_low) - push!(gaussians, i) - end - - println("Best convergent numerical value: ", E_list[end]) - + E_list, bases, gaussians, E_low = calculate_energies(num_gauss, w_transformed, K_transformed, b1, method) + else error("Invalid method provided!") end - + + # Plot results println("Best convergent numerical value: ", E_list[end]) - p = plot(gaussians, E_list, marker=:circle, label="Numerical result",linewidth=2) - title!(p, "S-wave convergence") + p = plot(gaussians, E_list, marker=:circle, label="Numerical result", linewidth=2) + title!(p, "S-wave Convergence") xlabel!(p, "Number of Gaussians") ylabel!(p, "Energy [Hartree]") - if plot_result==true + if plot_result display(p) - end + end + return p, E_list[end], bases end -""" - run_simulation_nuclear(ngauss=2, dim=2, bmax=5) +""" Run a nuclear simulation and print the final energy. - -# Arguments -- `ngauss`: Number of Gaussian functions to use in the simulation (default is 2). -- `dim`: Dimension of the simulation (default is 2). -- `bmax`: Maximum impact parameter (default is 5). - -# Outputs -Prints the final energy of the simulation. - -# Example -```julia -run_simulation_nuclear(3, 3, 10) """ function run_simulation_nuclear(ngauss, dim, bmax, masses, params) - return E_list, gaussians, eigenvectors, coords, masses = OptimizeGlobalParameters(ngauss, dim, bmax, masses, params) + return OptimizeGlobalParameters(ngauss, dim, bmax, masses, params) end From bd661b955e5dfb188d4a7e1f2690ae7c989192a3 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 14:36:56 +0100 Subject: [PATCH 03/10] minor refactoring and updated examples --- Examples/HydrogenAnion.jl | 2 +- Examples/Positronium.jl | 2 +- src/matrix_elements.jl | 2 +- src/sampling.jl | 40 +++++++++++++++++++++++++++++++++------ 4 files changed, 37 insertions(+), 9 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 24c3ba8..cce5095 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -7,7 +7,7 @@ w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] K = [0 0 0; 0 1/2 0; 0 0 1/2] J, U = jacobi_transform(masses) K_transformedformed = J * K * J' -w_transformedformedformed = [U' * w_list[i] for i in 1:length(w_list)] +w_transformed = [U' * w for w in w_list] Theortical_value = -0.527751016523 p, Energy, bases = run_simulation(50, :quasirandom, w_transformedformedformed, K_transformedformed) diff --git a/Examples/Positronium.jl b/Examples/Positronium.jl index 8ca59ca..2acd7f4 100644 --- a/Examples/Positronium.jl +++ b/Examples/Positronium.jl @@ -6,7 +6,7 @@ masses = [1.0,1.0,1.0] K = [1/2 0 0; 0 1/2 0; 0 0 1/2] J, U = jacobi_transform(masses) K_transformed = J * K * J' -w_transformed = [U' * w_list[i] for i in 1:length(w_list)] +w_transformed = [U' * w for w in w_list] Theortical_value = -0.2620050702328 diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index fe75cbd..577ce5f 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -115,7 +115,7 @@ function S_energy(bij::Vector{Float64}, K::Matrix{Float64}, w::Vector{Vector{Flo overlap, kinetic, Coulomb = S_wave(α, K, w) H = kinetic + Coulomb E, _ = eigen(H, overlap) - E0 = minimum(real(E)) + E0 = minimum(real(E)) return E0 end diff --git a/src/sampling.jl b/src/sampling.jl index 9f8b470..8f5a4cd 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -13,14 +13,12 @@ function corput(n, b=3) return q end -# Generate the nth d-dimensional point in the Halton sequence function halton(n, d) base = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281] @assert d <= length(base) "Dimension exceeds available bases." return [corput(n, base[i]) for i in 1:d] end -# Helper function for calculating energies in quasi-random methods function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method::Symbol) E_list, bases, gaussians = Float64[], [], Int[] E_low = Inf @@ -31,7 +29,7 @@ function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method: for j in 1:length(w_transformed):length(hal) base_segment = bij[j:j+length(w_transformed)-1] - E0 = S_energy(base_segment, K_transformed, w_transformed) + E0 = S_energy(Float64(base_segment), K_transformed, w_transformed) if E0 < best_energy best_energy = E0 best_base = copy(base_segment) @@ -47,14 +45,41 @@ function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method: return E_list, bases, gaussians, E_low end -# Main simulation function function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, plot_result::Bool=true) b1 = 10.0 + E_list = Float64[] + gaussians = Int[] + bij = Float64[] + E_low = 1e10 + bases = Vector{Vector{Float64}}() # Initialize as Vector{Vector{Float64}} + base_test = Float64[] # Initialize as Vector{Float64} to avoid Any type issues + base_curr = Float64[] # Initialize `base_curr` as an empty vector - # Choose appropriate method if method == :quasirandom println("---------QUASI-RANDOM METHOD---------") - E_list, bases, gaussians, E_low = calculate_energies(num_gauss, w_transformed, K_transformed, b1, method) + + for i in 1:num_gauss + hal = halton(i, 15 * length(w_transformed)) + bij = Float64.(-log.(hal) .* b1) # Ensure bij is Vector{Float64} + + for j in 1:length(w_transformed):length(hal) + append!(base_test, bij[j:j+length(w_transformed)-1]) + + # Call S_energy with correctly typed base_test + E0 = S_energy(Float64.(base_test), K_transformed, w_transformed) + + if E0 <= E_low + E_low = E0 + base_curr = copy(Float64.(bij[j:j+length(w_transformed)-1])) # Ensure base_curr is Vector{Float64} + end + base_test = base_test[1:end-length(w_transformed)] # Reset base_test + end + push!(bases, base_curr) # Append to bases as Vector{Float64} + append!(base_test, base_curr) + push!(E_list, E_low) + println(E_low) + push!(gaussians, i) + end elseif method == :quasirandomrefined println("---------QUASI-RANDOM METHOD WITH REFINEMENT---------") @@ -97,6 +122,9 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve return p, E_list[end], bases end + + + """ Run a nuclear simulation and print the final energy. """ From 9345719ddc7143fc49561fb08f8e015daa9009ff Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 14:38:16 +0100 Subject: [PATCH 04/10] using eachindex instead --- Examples/PionNucleon.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Examples/PionNucleon.jl b/Examples/PionNucleon.jl index 31d2a60..80d12ea 100644 --- a/Examples/PionNucleon.jl +++ b/Examples/PionNucleon.jl @@ -16,7 +16,7 @@ grid_points = range(rmin,rmax,3000) Φ = zeros(length(grid_points), length(coordinates)) -for i in 1:length(coordinates) +for i in eachindex(coordinates) local ϕ = zeros(length(grid_points)) ϕ_sum = zeros(length(grid_points)) rs = coordinates[i] From 1886b0f60d5b21378594c099c5edb51fa1edaf91 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 15:21:03 +0100 Subject: [PATCH 05/10] updated tests --- Examples/HydrogenAnion.jl | 4 +- Examples/PionNucleon.jl | 1 + src/FewBodyPhysics.jl | 11 ++-- src/constants.jl | 2 - src/matrix_elements.jl | 5 +- src/sampling.jl | 5 -- test/runtests.jl | 108 +++++++++++++++++++------------------- 7 files changed, 66 insertions(+), 70 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index cce5095..84b59d4 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -6,9 +6,9 @@ w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] K = [0 0 0; 0 1/2 0; 0 0 1/2] J, U = jacobi_transform(masses) -K_transformedformed = J * K * J' +K_transformed = J * K * J' w_transformed = [U' * w for w in w_list] Theortical_value = -0.527751016523 -p, Energy, bases = run_simulation(50, :quasirandom, w_transformedformedformed, K_transformedformed) +p, Energy, bases = run_simulation(50, :quasirandom, w_transformed, K_transformed) diff --git a/Examples/PionNucleon.jl b/Examples/PionNucleon.jl index 80d12ea..b383f5c 100644 --- a/Examples/PionNucleon.jl +++ b/Examples/PionNucleon.jl @@ -1,4 +1,5 @@ using Plots +using .FewBodyPhysics b = 3.9 S = 41.5 diff --git a/src/FewBodyPhysics.jl b/src/FewBodyPhysics.jl index 9511fd0..7e3a15c 100644 --- a/src/FewBodyPhysics.jl +++ b/src/FewBodyPhysics.jl @@ -1,13 +1,14 @@ module FewBodyPhysics +export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare +export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates +export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters +export corput, halton, run_simulation, run_simulation_nuclear + +include("constants.jl") include("coordinates.jl") include("matrix_elements.jl") include("sampling.jl") -include("constants.jl") -export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare -export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates -export S_elements, S_wave, S_energy, P_elements, pion_nucleon,compute_eigen_system, get_minimum_energy, optimize_global_parameters -export corput, halton, run_simulation, run_simulation_nuclear end diff --git a/src/constants.jl b/src/constants.jl index ae3ada8..403dab3 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -1,7 +1,5 @@ export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare -# Physical Constants - const m_p = 938.27 # Proton mass in MeV/c^2 const m_n = 939.57 # Neutron mass in MeV/c^2 const m_pi0 = 134.98 # Neutral pion mass (π⁰) in MeV/c^2 diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index 577ce5f..1ad2586 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -143,9 +143,8 @@ function P_elements(a::Vector{Float64}, b::Vector{Float64}, A::PositiveDefiniteS @assert isposdef(D) "Matrix D must be positive definite." R = inv(D) M0, trace = S_elements(A.matrix, B.matrix, K) - M1 = 0.5 * (a' * R * b) * M0 # Overlap perturbation term + M1 = 0.5 * (a' * R * b) * M0 - # Kinetic energy perturbation term kinetic = 6 * trace * M1 kinetic += (a' * K * b) * M0 kinetic -= (a' * K * A.matrix * R * b) * M0 @@ -280,8 +279,8 @@ Perform global optimization over a given parameter space to find optimal paramet function OptimizeGlobalParameters(ngauss, dim, bmax, masses, params) E_list = [] gaussians = [] - coords = [] eigenvectors = [] + coords = [] global E0S = 0.0 masses_min = copy(masses) diff --git a/src/sampling.jl b/src/sampling.jl index 8f5a4cd..7ab45e6 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -2,7 +2,6 @@ using Plots export corput, halton, run_simulation, run_simulation_nuclear -# Generate the nth element of the van der Corput sequence in base b function corput(n, b=3) q, bk = 0.0, 1 / b while n > 0 @@ -108,7 +107,6 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve error("Invalid method provided!") end - # Plot results println("Best convergent numerical value: ", E_list[end]) p = plot(gaussians, E_list, marker=:circle, label="Numerical result", linewidth=2) title!(p, "S-wave Convergence") @@ -122,9 +120,6 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve return p, E_list[end], bases end - - - """ Run a nuclear simulation and print the final energy. """ diff --git a/test/runtests.jl b/test/runtests.jl index 93358c8..7b92ece 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,69 +1,71 @@ -using FewBodyPhysics using Test using LinearAlgebra +include("../src/FewBodyPhysics.jl") + +using .FewBodyPhysics @testset "FewBodyPhysics.jl" begin - @testset "Coordinate Transformations" begin - masses = [1.0, 2.0, 3.0] - J, U = JacobiTransform(masses) - @test size(J) == (2, 3) - @test size(U) == (3, 2) - - bij = [0.5, 0.8] - w_list = [generate_weight_vector(3, 1, 2), generate_weight_vector(3, 2, 3)] - A = generate_A_matrix(bij, w_list) - @test size(A) == (3, 3) - - α = [1.0, 2.0, 3.0] - transformed_list = transform_list(α) - @test length(transformed_list) == 3 - @test transformed_list[1] == [1.0] - - a = rand(3, 2) - b = rand(3, 2) - sum_val = shift_vectors(a, b) - @test isa(sum_val, Float64) - - w = generate_weight_vector(3, 1, 2) - @test w == [1, -1, 0] - - r = rand(3) - x = transform_coordinates(J, r) - r_back = inverse_transform_coordinates(U, x) - @test norm(r - r_back) < 1e-10 - end - @testset "corput tests" begin + @testset "Corput Sequence Tests" begin + # Basic tests for base 2 @test corput(1, 2) ≈ 0.5 @test corput(2, 2) ≈ 0.25 @test corput(3, 2) ≈ 0.75 + @test corput(4, 2) ≈ 0.125 + @test corput(5, 2) ≈ 0.625 + + # Additional tests for base 3 + @test corput(1, 3) ≈ 1/3 + @test corput(2, 3) ≈ 2/3 + @test corput(3, 3) ≈ 1/9 + @test corput(4, 3) ≈ 4/9 + @test corput(5, 3) ≈ 7/9 + + # Additional tests for base 5 + @test corput(1, 5) ≈ 0.2 + @test corput(2, 5) ≈ 0.04 + @test corput(3, 5) ≈ 0.6 + @test corput(4, 5) ≈ 0.08 + @test corput(5, 5) ≈ 0.16 + @test corput(10, 5) ≈ 0.32 + end - @testset "halton tests" begin - @test halton(1, 2) ≈ [0.5, 0.3333333333333333] - @test halton(2, 2) ≈ [0.25, 0.6666666666666666] + @testset "Halton Sequence Tests" begin + # Basic sequence tests for dimension 2 + @test halton(1, 2) ≈ [0.5, 1/3] + @test halton(2, 2) ≈ [0.25, 2/3] + @test halton(3, 2) ≈ [0.75, 1/9] + @test halton(4, 2) ≈ [0.125, 4/9] + @test halton(5, 2) ≈ [0.625, 7/9] + + # Additional tests for higher dimensions + @test halton(1, 3) ≈ [0.5, 1/3, 0.2] + @test halton(2, 3) ≈ [0.25, 2/3, 0.4] + @test halton(3, 3) ≈ [0.75, 1/9, 0.6] + @test halton(4, 3) ≈ [0.125, 4/9, 0.8] + + + # Edge cases: checking very high `n` + @test halton(1000, 2) ≈ [corput(1000, 2), corput(1000, 3)] + @test halton(5000, 3) ≈ [corput(5000, 2), corput(5000, 3), corput(5000, 5)] + + # Edge case: requesting only 1 dimension + @test halton(1, 1) ≈ [0.5] + @test halton(10, 1) ≈ [corput(10, 2)] + + # Edge case: invalid dimension exceeding base list @test_throws AssertionError halton(1, 100) - end - @testset "run_simulation tests" begin - w_transformed = [1, 2, 3] - K_transformed = [1, 2, 3] - p, E_list, bases = run_simulation(15, :quasirandom, w_transformed, K_transformed, false) - @test E_list isa Float64 - @test bases isa Array - @test_throws ErrorException run_simulation(15, :invalidmethod, w_transformed, K_transformed, false) - end + # Edge case: `n = 0` should return zero for any dimension + @test halton(0, 2) ≈ [0.0, 0.0] + @test halton(0, 5) ≈ [0.0, 0.0, 0.0, 0.0, 0.0] - @testset "run_simulation_nuclear tests" begin - masses = [1, 2, 3] - params = [1, 2, 3] - E_list, gaussians, eigenvectors, coords, masses = run_simulation_nuclear(2, 2, 5, masses, params) - @test E_list isa Array - @test gaussians isa Array - @test eigenvectors isa Array - @test coords isa Array - @test masses isa Array + # Random large n tests for robustness + @test halton(123, 4) ≈ [corput(123, 2), corput(123, 3), corput(123, 5), corput(123, 7)] + @test halton(999, 3) ≈ [corput(999, 2), corput(999, 3), corput(999, 5)] end - + end + From 4553c956163a2eb89b265d0b887aa5bcbfdd9388 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 15:51:18 +0100 Subject: [PATCH 06/10] added documentation --- Examples/HydrogenAnion.jl | 2 +- docs/src/index.md | 28 +++++++++++++++------------- src/sampling.jl | 26 +++++++++++++++++--------- 3 files changed, 33 insertions(+), 23 deletions(-) diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl index 84b59d4..111bb42 100644 --- a/Examples/HydrogenAnion.jl +++ b/Examples/HydrogenAnion.jl @@ -1,4 +1,4 @@ -using Revise +using Revise using .FewBodyPhysics masses = [1e15, 1.0, 1.0] diff --git a/docs/src/index.md b/docs/src/index.md index c6e6a62..2c81085 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -16,40 +16,42 @@ H = - \sum_{i=1}^{3} \frac{1}{2m_i}\frac{\partial^2}{\partial \boldsymbol{r}_i^2 ``` The masses of the three constituents are `m_i = {1, 1, 1}` and the charges `q_i = {+1, −1, −1}`. We can estimate the ground state of this Coulombic three-body system using 50 Gaussians -### Psudorandom +### Quasirandom ```@example 1 using Plots, FewBodyPhysics w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] -masses = [1,1,1] +masses = [1.0,1.0,1.0] K = [1/2 0 0; 0 1/2 0; 0 0 1/2] -J, U = Ω(masses) +J, U = jacobi_transform(masses) K_transformed = J * K * J' -w_transformed = [U' * w_list[i] for i in 1:length(w_list)] +w_transformed = [U' * w for w in w_list] + Theortical_value = -0.2620050702328 -p, Energy = run_simulation(50,:psudorandom, w_transformed, K_transformed) +p, Energy = run_simulation(50,:quasirandom, w_transformed, K_transformed) plot(p) ``` With a difference in energy ```@example 1 -@show Theortical_value - Energy +@show Energy-Theortical_value ``` -### Quasirandom +### Psudorandom And similarly for a quasirandom ```@example 2 using Plots, FewBodyPhysics w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] -masses = [1,1,1] +masses = [1.0,1.0,1.0] K = [1/2 0 0; 0 1/2 0; 0 0 1/2] -J, U = Ω(masses) +J, U = jacobi_transform(masses) K_transformed = J * K * J' -w_transformed = [U' * w_list[i] for i in 1:length(w_list)] +w_transformed = [U' * w for w in w_list] + Theortical_value = -0.2620050702328 -p, Energy = run_simulation(50,:quasirandom, w_transformed, K_transformed) +p, Energy = run_simulation(50,:psudorandom, w_transformed, K_transformed) plot(p) ``` With a difference in energy @@ -68,7 +70,7 @@ b = 3.9 S = 41.5 params = [b, S] -masses = [(m_p+m_n)/2, m_π] +masses = [(m_p+m_n)/2, m_pi] energies, Gaussians, eigenvectors, coordinates, masses = run_simulation_nuclear(5,2,5,masses,params) @@ -80,7 +82,7 @@ grid_points = range(rmin,rmax,3000) Φ = zeros(length(grid_points), length(coordinates)) -for i in 1:length(coordinates) +for i in eachindex(coordinates) local ϕ = zeros(length(grid_points)) ϕ_sum = zeros(length(grid_points)) rs = coordinates[i] diff --git a/src/sampling.jl b/src/sampling.jl index 7ab45e6..3da9777 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -21,14 +21,21 @@ end function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method::Symbol) E_list, bases, gaussians = Float64[], [], Int[] E_low = Inf + for i in 1:num_gauss hal = halton(i, 15 * length(w_transformed)) - bij = Float64.(-log.(hal) * b1) + + # Correctly create `bij` as a Vector{Float64} + bij = -log.(hal) * b1 + best_energy, best_base = E_low, nothing - for j in 1:length(w_transformed):length(hal) + for j in 1:length(w_transformed):length(bij) + # Ensure that `base_segment` is correctly extracted as a Vector{Float64} base_segment = bij[j:j+length(w_transformed)-1] - E0 = S_energy(Float64(base_segment), K_transformed, w_transformed) + + # Pass `base_segment` as a vector to `S_energy` + E0 = S_energy(base_segment, K_transformed, w_transformed) if E0 < best_energy best_energy = E0 best_base = copy(base_segment) @@ -44,6 +51,7 @@ function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method: return E_list, bases, gaussians, E_low end + function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, plot_result::Bool=true) b1 = 10.0 E_list = Float64[] @@ -84,20 +92,20 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve println("---------QUASI-RANDOM METHOD WITH REFINEMENT---------") E_list, bases, gaussians, E_low = calculate_energies(num_gauss, w_transformed, K_transformed, b1, method) - # Additional refinement step - for i in 1:length(bases)-length(w_transformed) - rand_ref = Float64.(-log.(rand(200 * length(w_transformed))) * b1) + for i in 1:length(bases) - length(w_transformed) + rand_ref = -log.(rand(200 * length(w_transformed))) * b1 for j in 1:length(w_transformed):length(rand_ref) - refined_base = bases[i][1:end-length(w_transformed)] # Adjust size accordingly - refined_base .= rand_ref[j:j+length(w_transformed)-1] + refined_base = rand_ref[j:j+length(w_transformed)-1] E_refined = S_energy(refined_base, K_transformed, w_transformed) + if E_refined < E_low E_low = E_refined - bases[i] .= refined_base + bases[i] = copy(refined_base) # Replace with `copy` instead of broadcasting end end println("Refined Energy: ", E_low) end + elseif method == :psudorandom println("---------PSEUDO-RANDOM METHOD---------") From ff3c8ce8f60b220679064e4a0810bc361f8e9e94 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 16:02:06 +0100 Subject: [PATCH 07/10] added documentation --- docs/src/API.md | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/docs/src/API.md b/docs/src/API.md index e69cd23..115090b 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -2,13 +2,14 @@ ## Coordinates ```@docs -Ω -A_generate +ParticleSystem +jacobi_transform +generate_A_matrix transform_list -shift -w_gen +shift_vectors +generate_weight_vector transform_coordinates -transform_back +inverse_transform_coordinates ``` ## Matrix elements @@ -19,9 +20,9 @@ S_wave S_energy P_elements pion_nucleon -ComputeEigenSystem -GetMinimumEnergy -OptimizeGlobalParameters +compute_eigen_system +get_minimum_energy +optimize_global_parameters ``` ## Sampling From 59427ba2f73aa3ff664396bdb30e1dd63d334ad4 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 16:23:23 +0100 Subject: [PATCH 08/10] added documentation --- docs/src/examples.md | 10 +++++----- src/constants.jl | 2 +- src/matrix_elements.jl | 3 +-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/src/examples.md b/docs/src/examples.md index 6830e11..e866c2e 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -5,16 +5,16 @@ Suppose you want to calculate the ground state energy of the hydrogen anion in t ```@example example1 using Plots, FewBodyPhysics -masses = [Inf, 1, 1] +masses = [1e15, 1.0, 1.0] w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ] K = [0 0 0; 0 1/2 0; 0 0 1/2] -J, U = Ω(masses) -K_transformedformed = J * K * J' -w_transformedformedformed = [U' * w_list[i] for i in 1:length(w_list)] +J, U = jacobi_transform(masses) +K_transformed = J * K * J' +w_transformed = [U' * w for w in w_list] Theortical_value = -0.527751016523 -p, Energy, bases = run_simulation(50, :quasirandom, w_transformedformedformed, K_transformedformed) +p, Energy, bases = run_simulation(50, :quasirandom, w_transformed, K_transformed) plot(p) ``` diff --git a/src/constants.jl b/src/constants.jl index 403dab3..994e3a3 100644 --- a/src/constants.jl +++ b/src/constants.jl @@ -9,7 +9,7 @@ const ħc = 197.3269804 # Planck's constant times speed of light in MeV·fm const m_bare = (m_p + m_n) / 2 # Average nucleon mass in MeV/c^2 """ - μ = (m_bare * m_pi0) / (m_bare + m_pi0) +μ = (m_bare * m_pi0) / (m_bare + m_pi0) Reduced mass of the nucleon-pion system in MeV/c^2. """ diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index 1ad2586..fdf8aed 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -1,8 +1,7 @@ using LinearAlgebra using Optim -export S_elements, S_wave, S_energy, P_elements, pion_nucleon, - compute_eigen_system, get_minimum_energy, optimize_global_parameters +export S_elements, S_wave, S_energy, P_elements, pion_nucleon,ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters struct PositiveDefiniteSymmetricMatrix{T<:Real} matrix::Matrix{T} From e4b2109e88ec06bca359df790b03c6444567f4c5 Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 16:47:01 +0100 Subject: [PATCH 09/10] added documentation --- docs/src/API.md | 6 +++--- src/coordinates.jl | 15 +++++++++++++++ src/sampling.jl | 35 +++++++++++++++++++++++++++++++++++ 3 files changed, 53 insertions(+), 3 deletions(-) diff --git a/docs/src/API.md b/docs/src/API.md index 115090b..8148e00 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -20,9 +20,9 @@ S_wave S_energy P_elements pion_nucleon -compute_eigen_system -get_minimum_energy -optimize_global_parameters +ComputeEigenSystem +GetMinimumEnergy +OptimizeGlobalParameters ``` ## Sampling diff --git a/src/coordinates.jl b/src/coordinates.jl index 76d9315..f472f02 100644 --- a/src/coordinates.jl +++ b/src/coordinates.jl @@ -2,6 +2,21 @@ using LinearAlgebra export ParticleSystem, jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates +""" + ParticleSystem(masses::Vector{Float64}) + +A data structure representing a system of particles, storing their masses and associated Jacobi transformation matrices for coordinate transformations. + +# Fields +- `masses::Vector{Float64}`: A vector containing the masses of the particles. +- `J::Matrix{Float64}`: The Jacobi transformation matrix, used to convert particle coordinates into Jacobi coordinates. +- `U::Matrix{Float64}`: The pseudoinverse of the Jacobi transformation matrix `J`, used to transform Jacobi coordinates back to the particle coordinate system. + +# Constructor +- `ParticleSystem(masses::Vector{Float64})`: Constructs a new `ParticleSystem` instance. + - **Arguments**: + - `masses`: A vector of particle masses. At least two masses are required. +""" struct ParticleSystem masses::Vector{Float64} J::Matrix{Float64} diff --git a/src/sampling.jl b/src/sampling.jl index 3da9777..4662bd1 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -1,7 +1,22 @@ using Plots export corput, halton, run_simulation, run_simulation_nuclear +""" + corput(n::Int, b::Int=3) -> Float64 + +Generates the nth term of the van der Corput sequence in the given base `b`, which is often used for quasi-random number generation. + +# Arguments +- `n::Int`: The position of the term in the sequence to calculate. +- `b::Int=3`: The base for the sequence. Defaults to 3 if not provided. +# Returns +- `Float64`: The nth term in the van der Corput sequence for the given base `b`. + +# Example +```julia +corput(1, 2) # Returns 0.5 for base 2 +""" function corput(n, b=3) q, bk = 0.0, 1 / b while n > 0 @@ -12,12 +27,28 @@ function corput(n, b=3) return q end +""" + halton(n::Int, d::Int) -> Vector{Float64} + +Generates a point in the Halton sequence with d dimensions, used in quasi-random sampling. + + # Arguments + - `n::Int`: The index of the sequence point to generate. + - `d::Int`: The dimensionality of the Halton sequence (i.e., the number of bases to use). + # Returns + Vector{Float64}: A vector of length d representing the nth point in the Halton sequence. +""" function halton(n, d) base = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281] @assert d <= length(base) "Dimension exceeds available bases." return [corput(n, base[i]) for i in 1:d] end +""" + calculate_energies(num_gauss::Int, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, b1::Float64, method::Symbol) -> Tuple{Vector{Float64}, Vector{Vector{Float64}}, Vector{Int}, Float64} + +Calculates and refines energies for a set of Gaussian basis functions using quasi-random or pseudo-random methods. +""" function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method::Symbol) E_list, bases, gaussians = Float64[], [], Int[] E_low = Inf @@ -51,7 +82,11 @@ function calculate_energies(num_gauss, w_transformed, K_transformed, b1, method: return E_list, bases, gaussians, E_low end +""" + run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, plot_result::Bool=true) -> Tuple{Plots.Plot, Float64, Vector{Vector{Float64}}} +Runs a simulation to calculate energy values for Gaussian basis functions and optionally plots the results. +""" function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Vector{Float64}}, K_transformed::Matrix{Float64}, plot_result::Bool=true) b1 = 10.0 E_list = Float64[] From b5040e441f3ca610da59677e0db0110b71cbfcff Mon Sep 17 00:00:00 2001 From: MartinMikkelsen Date: Tue, 29 Oct 2024 16:54:42 +0100 Subject: [PATCH 10/10] added documentation --- docs/src/API.md | 6 ++++++ src/FewBodyPhysics.jl | 3 +-- src/matrix_elements.jl | 2 +- src/sampling.jl | 2 +- 4 files changed, 9 insertions(+), 4 deletions(-) diff --git a/docs/src/API.md b/docs/src/API.md index 8148e00..ec08e42 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -30,6 +30,12 @@ OptimizeGlobalParameters ```@docs corput halton +calculate_energies run_simulation run_simulation_nuclear ``` +## Constants + +```@docs +μ +``` \ No newline at end of file diff --git a/src/FewBodyPhysics.jl b/src/FewBodyPhysics.jl index 7e3a15c..9498031 100644 --- a/src/FewBodyPhysics.jl +++ b/src/FewBodyPhysics.jl @@ -3,12 +3,11 @@ module FewBodyPhysics export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters -export corput, halton, run_simulation, run_simulation_nuclear +export corput, halton, calculate_energies, run_simulation, run_simulation_nuclear include("constants.jl") include("coordinates.jl") include("matrix_elements.jl") include("sampling.jl") - end diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl index fdf8aed..caac672 100644 --- a/src/matrix_elements.jl +++ b/src/matrix_elements.jl @@ -1,7 +1,7 @@ using LinearAlgebra using Optim -export S_elements, S_wave, S_energy, P_elements, pion_nucleon,ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters +export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters struct PositiveDefiniteSymmetricMatrix{T<:Real} matrix::Matrix{T} diff --git a/src/sampling.jl b/src/sampling.jl index 4662bd1..6044270 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -1,6 +1,6 @@ using Plots -export corput, halton, run_simulation, run_simulation_nuclear +export corput, halton, calculate_energies, run_simulation, run_simulation_nuclear """ corput(n::Int, b::Int=3) -> Float64