From b763fce0ac223940fb68afc2649ea33dd2e768fb Mon Sep 17 00:00:00 2001 From: Shuhei Ohno Date: Fri, 14 Aug 2026 13:12:23 +0900 Subject: [PATCH] Simplify Rayleigh-Ritz calculations --- docs/src/Free-Complement.md | 2 +- docs/src/Rayleigh-Ritz.md | 18 +- src/Rayleigh-Ritz.jl | 477 ++++++++++++------------------------ test/FC.jl | 2 +- test/Rayleigh-Ritz.jl | 38 ++- 5 files changed, 200 insertions(+), 337 deletions(-) diff --git a/docs/src/Free-Complement.md b/docs/src/Free-Complement.md index 0946849..44e0faf 100644 --- a/docs/src/Free-Complement.md +++ b/docs/src/Free-Complement.md @@ -142,7 +142,7 @@ reference_energies = [ @printf("%5s %-18s %s\n", "M_n", "This work", "Ref.") println("----- ------------------ ------------") for reference_energy in reference_energies - result = solve(H, basisset, info=-1) + result = solve(H, basisset) @printf( "%5d %18.15f %s\n", length(basisset), diff --git a/docs/src/Rayleigh-Ritz.md b/docs/src/Rayleigh-Ritz.md index dd70dbb..245c04a 100644 --- a/docs/src/Rayleigh-Ritz.md +++ b/docs/src/Rayleigh-Ritz.md @@ -136,14 +136,15 @@ BS = BasisSet( nothing # hide ``` -Solve the eigenvalue problem. You should find +Configure the method and solve the eigenvalue problem. You should find ```math E_{n=1} = -0.499278~E_\mathrm{h}, ``` which is amazingly good for only four basis functions according to [Thijssen(2007)](https://doi.org/10.1017/CBO9781139171397). The exact ground-state energy is ``-0.5~E_\mathrm{h}``. ```@repl example -solve(H, BS) +method = RayleighRitzMethod(BS) +solve(H, method) ``` ## Example of Hydrogen Atom @@ -297,12 +298,13 @@ The results agree with those in Table 1 of the Supporting Information. The small ### Solver ```@docs; canonical=false -solve(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4) -solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4) -solve(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4) -optimize(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) -optimize(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=1, progress=true, optimizer=Optim.NelderMead(), options...) -optimize(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) +RayleighRitzMethod +ResultRayleighRitz +ResultOptimizationRayleighRitz +solve(hamiltonian::Hamiltonian, method::RayleighRitzMethod; perturbation=Hamiltonian()) +optimize +expectation +verification ``` ### Basis Set diff --git a/src/Rayleigh-Ritz.jl b/src/Rayleigh-Ritz.jl index 4ac8c17..b1d3f1d 100644 --- a/src/Rayleigh-Ritz.jl +++ b/src/Rayleigh-Ritz.jl @@ -1,343 +1,189 @@ -export solve, optimize +export RayleighRitzMethod, ResultRayleighRitz, ResultOptimizationRayleighRitz, + solve, optimize, expectation, verification import LinearAlgebra import Optim -import Printf import QuadGK import SpecialFunctions import Subscripts -# result +# method and results + +struct RayleighRitzMethod{B<:Union{BasisSet,GeometricBasisSet}} + basisset::B +end + +RayleighRitzMethod(basis::Basis) = RayleighRitzMethod(BasisSet(basis)) + +Base.string(method::RayleighRitzMethod) = "RayleighRitzMethod($(method.basisset))" +Base.show(io::IO, method::RayleighRitzMethod) = print(io, Base.string(method)) struct ResultRayleighRitz data::Any ResultRayleighRitz(; args...) = new(NamedTuple(Dict(args))) end +struct ResultOptimizationRayleighRitz + data::Any + ResultOptimizationRayleighRitz(; args...) = new(NamedTuple(Dict(args))) +end + function ψ(result::ResultRayleighRitz, r; n::Int=1) return sum(result.C[i,n] * TwoBody.φ(result.basisset[i], r) for i in 1:result.nₘₐₓ) end -Base.getproperty(result::ResultRayleighRitz, symbol::Symbol) = Base.getproperty(getfield(result,:data), symbol) -Base.haskey(result::ResultRayleighRitz, symbol::Symbol) = Base.haskey(getfield(result,:data), symbol) +function Base.getproperty(result::ResultRayleighRitz, symbol::Symbol) + symbol === :data && return getfield(result, :data) + return Base.getproperty(getfield(result, :data), symbol) +end + +Base.haskey(result::ResultRayleighRitz, symbol::Symbol) = + Base.haskey(getfield(result, :data), symbol) + +function Base.getproperty(result::ResultOptimizationRayleighRitz, symbol::Symbol) + symbol === :data && return getfield(result, :data) + data = getfield(result, :data) + return Base.haskey(data, symbol) ? Base.getproperty(data, symbol) : Base.getproperty(data.result, symbol) +end + +Base.haskey(result::ResultOptimizationRayleighRitz, symbol::Symbol) = + Base.haskey(getfield(result, :data), symbol) || Base.haskey(getfield(result, :data).result, symbol) + Base.show(io::IO, result::ResultRayleighRitz) = print(io, Base.string(result)) +Base.show(io::IO, result::ResultOptimizationRayleighRitz) = print(io, Base.string(result)) -function Base.string(result::ResultRayleighRitz) - if 0 ≤ result.info - # method - text = "# method\n\n" - text *= "Rayleigh-Ritz method with $(typeof(result.basisset[1]))\n" - text *= "J. Thijssen, Computational Physics 2nd Edition (2013)\n" - text *= "https://doi.org/10.1017/CBO9781139171397\n" - # optimization - if haskey(result, :optimizer) - text *= "\n# optimizer\n\n" - text *= string(result.optimizer) * "\n" - text *= "P. K. Mogensen, A. N. Riseth, J. Open Source Softw., 3(24), 615 (2018)\n" - text *= "https://doi.org/10.21105/joss.00615\n" - # optimization log - if result.progress - text *= "\n# optimization log\n\n" - for (E, x) in result.history - if isinf(E) - text *= Printf.@sprintf("%+.9e %s\n", E, "[" * join([Printf.@sprintf("%+.3e", x[i]) for i in keys(x)], ", ") *"]") - else - text *= Printf.@sprintf("%+.9e %s\n", E, "[" * join([Printf.@sprintf("%+.3e", x[i]) for i in keys(x)], ", ") *"]") - end - end - end - # initial basis set - if haskey(result, :geometricbasisset) - # initial geometric progression - text *= "\n# initial geometric progression\n\n" - text *= "type \t$(result.initialgeometricbasisset.basistype)\n" - text *= string("range\tr", Subscripts.sub("$(result.initialgeometricbasisset.nₘᵢₙ)"), " - r", Subscripts.sub("$(result.initialgeometricbasisset.nₘₐₓ)"), "\n") - text *= string("r", Subscripts.sub("$(result.initialgeometricbasisset.nₘᵢₙ)"), " \t", result.initialgeometricbasisset.r₁, "\n") - text *= string("r", Subscripts.sub("$(result.initialgeometricbasisset.n)" ), " \t", result.initialgeometricbasisset.rₙ, "\n") - else - # initial basis function - text *= "\n# initial basis function\n\n" - for n in 1:result.nₘₐₓ - text *= Printf.@sprintf("φ%s(r) = TwoBody.φ(%s, r)\n", Subscripts.sub("$n"), result.initialbasisset[n]) - end - end - end - # geometric progression - if haskey(result, :geometricbasisset) && 0 < result.info - if haskey(result, :optimizer) - text *= "\n# optimized geometric progression\n\n" - else - text *= "\n# geometric progression\n\n" - end - text *= "type \t$(result.geometricbasisset.basistype)\n" - text *= string("range\tr", Subscripts.sub("$(result.geometricbasisset.nₘᵢₙ)"), " - r", Subscripts.sub("$(result.geometricbasisset.nₘₐₓ)"), "\n") - text *= string("r", Subscripts.sub("$(result.geometricbasisset.nₘᵢₙ)"), " \t", result.geometricbasisset.r₁, "\n") - text *= string("r", Subscripts.sub("$(result.geometricbasisset.n)" ), " \t", result.geometricbasisset.rₙ, "\n") - end - # basis set - if haskey(result, :optimizer) && !haskey(result, :geometricbasisset) - text *= "\n# optimized basis function\n\n" - else - text *= "\n# basis function\n\n" - end - for n in 1:result.nₘₐₓ - text *= Printf.@sprintf("φ%s(r) = TwoBody.φ(%s, r)\n", Subscripts.sub("$n"), result.basisset[n]) - end - # eigenfunction - text *= "\n# eigenfunction\n\n" - for n in 1:result.nₘₐₓ - text *= Printf.@sprintf("ψ%s(r) = ", Subscripts.sub("$n")) - if result.nₘₐₓ > 5 - text *= "( " - end - ncol = 5 - if mod(result.nₘₐₓ,5) == 0 - ncol = 5 - elseif mod(result.nₘₐₓ,4) == 0 - ncol = 4 - elseif mod(result.nₘₐₓ,3) == 0 - ncol = 3 - end - for i in 1:result.nₘₐₓ - text *= Printf.@sprintf("%s %.6fφ%s(r) ", result.C[i,n]<0 ? "-" : "+", abs(result.C[i,n]), Subscripts.sub("$i")) - if result.nₘₐₓ > 5 - if result.nₘₐₓ == i - text *= ")" - elseif mod(i,ncol) == 0 - text *= "\n " - end - end - end - text *= "\n" - end - # eigenvalue - text *= "\n# eigenvalue\n\n" - for n in 1:result.nₘₐₓ - text *= Printf.@sprintf("E%s = %s\n", Subscripts.sub("$n"), "$(result.E[n])") - end - if 0 < result.info - # verification - text *= "\n# verification\n" - # norm - text *= "\nn \tnorm, <ψₙ|ψₙ> = cₙ' * S * cₙ = 1\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", result.expectation[:S][n], "\n") - end - # ill-conditioned - text *= "\nn \till-conditioned, |<ψₙ|H|ψₙ> - E| = |cₙ' * H * cₙ - E| = 0\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", abs(result.expectation[:0][n]), "\n") - end - # expectation value - text *= "\n# expectation value\n" - # hamiltonian - text *= "\nn \thamiltonian, <ψₙ|H|ψₙ> = cₙ' * H * cₙ\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", result.expectation[:H][n], "\n") - end - # perturbation - if !isempty(result.perturbation.terms) - text *= "\nn \tperturbation\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", result.expectation[:perturbation][n], "\n") - end - text *= "\nn \thamiltonian + perturbation\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", result.expectation[:H][n] + result.expectation[:perturbation][n], "\n") - end - end - # each term - for term in [result.hamiltonian.terms..., result.perturbation.terms...] - text *= "\nn \texpectation value of $(term)\n" - for n in 1:min(result.nₘₐₓ, result.info) - text *= string("$n\t", result.expectation[term][n], "\n") - end - end - end - # end - return text - else - # info < 0 - text = "ResultRayleighRitz:\n" - text *= " E: $(typeof(result.E))\n" - for e in result.E - text *= " $(e)\n" - end - return text - end - return "error" +ψ(result::ResultOptimizationRayleighRitz, r; n::Int=1) = ψ(result.result, r; n=n) + +function expectation(result::ResultRayleighRitz, M::AbstractMatrix) + size(M) == size(result.H) || throw(DimensionMismatch("the matrix dimensions must match the result")) + return [LinearAlgebra.dot(result.C[:,n], M * result.C[:,n]) for n in eachindex(result.E)] end -# solver +expectation(result::ResultRayleighRitz, operator::Operator) = + expectation(result, matrix(operator, result.basisset)) -function solve(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4) +expectation(result::ResultRayleighRitz, hamiltonian::Hamiltonian) = + expectation(result, matrix(hamiltonian, result.basisset)) - # initialization - nₘₐₓ = length(basisset.basis) +expectation(result::ResultOptimizationRayleighRitz, operator) = + expectation(result.result, operator) - # matrix element - S = matrix(basisset) - H = matrix(hamiltonian, basisset) +function verification(result::ResultRayleighRitz) + norm = expectation(result, result.S) + residual = abs.(expectation(result, result.H) .- result.E) + return (; norm, residual) +end - # calculation - E, C = LinearAlgebra.eigen(H, S) +verification(result::ResultOptimizationRayleighRitz) = verification(result.result) - # expectation value - expectation = Dict() - if 0 < info - expectation[:S] = [C[:,n]' * S * C[:,n] for n in 1:min(nₘₐₓ, info)] - expectation[:H] = [C[:,n]' * H * C[:,n] for n in 1:min(nₘₐₓ, info)] - expectation[:0] = [expectation[:H][n] - E[n] for n in 1:min(nₘₐₓ, info)] - if !isempty(perturbation.terms) - M = matrix(perturbation, basisset) - expectation[:perturbation] = [C[:,n]' * M * C[:,n] for n in 1:min(nₘₐₓ, info)] - end - for term in [hamiltonian.terms..., perturbation.terms...] - M = matrix(term, basisset) - expectation[term] = [C[:,n]' * M * C[:,n] for n in 1:min(nₘₐₓ, info)] - end +function Base.string(result::ResultRayleighRitz) + check = verification(result) + text = "# method\n\n$(result.method)\n\n# eigenvalue\n\n" + for n in eachindex(result.E) + text *= "E$(Subscripts.sub(string(n))) = $(result.E[n])\n" end - - # return - if 0 < info - # full information - return ResultRayleighRitz(; - info = info, - hamiltonian = hamiltonian, - perturbation = perturbation, - basisset = basisset, - nₘₐₓ = nₘₐₓ, - H = H, - S = S, - E = E, - C = C, - expectation = expectation, - ) - else - # for optimization - return ResultRayleighRitz(; - info = info, - E = E, - ) + text *= "\n# verification\n\nn\tnorm\tresidual\n" + for n in eachindex(result.E) + text *= "$n\t$(check.norm[n])\t$(check.residual[n])\n" end + return text +end +function Base.string(result::ResultOptimizationRayleighRitz) + text = "# optimizer\n\n$(result.optimizer)\n" + if result.progress + text *= "\n# optimization log\n\nenergy\tparameters\n" + for entry in result.history + text *= "$(entry.energy)\t$(entry.parameters)\n" + end + end + return text * "\n" * Base.string(result.result) end -function optimize(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) +# solver - # optimize +function solve(hamiltonian::Hamiltonian, method::RayleighRitzMethod; perturbation=Hamiltonian()) + source = method.basisset + basisset = source isa GeometricBasisSet ? BasisSet(source.basis...) : source + nₘₐₓ = length(basisset.basis) + S = matrix(basisset) + H = matrix(hamiltonian, basisset) + E, C = LinearAlgebra.eigen(H, S) + common = (method=method, hamiltonian=hamiltonian, perturbation=perturbation, + basisset=basisset, nₘₐₓ=nₘₐₓ, H=H, S=S, E=E, C=C) + return source isa GeometricBasisSet ? + ResultRayleighRitz(; geometricbasisset=source, common...) : + ResultRayleighRitz(; common...) +end + +function optimize(hamiltonian::Hamiltonian, method::RayleighRitzMethod{<:BasisSet}; + perturbation=Hamiltonian(), progress=true, + optimizer=Optim.NelderMead(), options...) + basisset = method.basisset history = [] - res = Optim.optimize( + optimization = Optim.optimize( x -> try - E = solve( - hamiltonian, - BasisSet([_replace_exponent(basisset.basis[i], x[i]) for i in keys(basisset.basis)]...), - perturbation = perturbation, - info = -1 - ).E[1] - if 0 ≤ info - push!(history, (energy=E, parameters=x)) - end + varied = BasisSet([_replace_exponent(basisset.basis[i], x[i]) for i in keys(basisset.basis)]...) + E = solve(hamiltonian, RayleighRitzMethod(varied); perturbation=perturbation).E[1] + progress && push!(history, (energy=E, parameters=copy(x))) E catch - if 0 ≤ info - push!(history, (energy=Inf, parameters=x)) - end + progress && push!(history, (energy=Inf, parameters=copy(x))) Inf end, - [basisset.basis[i].a for i in keys(basisset.basis)], - method = optimizer, - options... + [float(basisset.basis[i].a) for i in keys(basisset.basis)], + optimizer, + Optim.Options(; options...) ) - # result - res = solve(hamiltonian, BasisSet([_replace_exponent(basisset.basis[i], res.minimizer[i]) for i in keys(basisset.basis)]...), perturbation=perturbation, info=info) - if 0 ≤ info - return ResultRayleighRitz(; - optimizer = optimizer, - initialbasisset = basisset, - options = options, - progress = progress, - history = history, - getfield(res, :data)..., - ) - else - return ResultRayleighRitz(; - getfield(res, :data)..., - ) - end - -end - -function solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4) - return solve(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info) -end - -function optimize(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=1, progress=true, optimizer=Optim.NelderMead(), options...) - return optimize(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info, progress=progress, optimizer=optimizer, options...) -end - -function solve(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4) - res = solve(hamiltonian, BasisSet(basisset.basis...); perturbation=perturbation, info=info) - if 0 ≤ info - return ResultRayleighRitz(; - geometricbasisset = basisset, - getfield(res, :data)..., - ) - else - return ResultRayleighRitz(; - getfield(res, :data)..., - ) - end + optimized = BasisSet([_replace_exponent(basisset.basis[i], optimization.minimizer[i]) for i in keys(basisset.basis)]...) + result = solve(hamiltonian, RayleighRitzMethod(optimized); perturbation=perturbation) + return ResultOptimizationRayleighRitz(; + result, initialmethod=method, optimizer, options, progress, history, optimization, + ) end -function optimize(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) - - # optimize +function optimize(hamiltonian::Hamiltonian, method::RayleighRitzMethod{<:GeometricBasisSet}; + perturbation=Hamiltonian(), progress=true, + optimizer=Optim.NelderMead(), options...) + basisset = method.basisset history = [] - res = Optim.optimize( + optimization = Optim.optimize( x -> try - E = solve( - hamiltonian, - GeometricBasisSet(basisset.basistype, x..., basisset.n, nₘₐₓ=basisset.nₘₐₓ, nₘᵢₙ=basisset.nₘᵢₙ), - perturbation = perturbation, - info = -1 - ).E[1] - if 0 ≤ info - push!(history, (energy=E, parameters=x)) - end + varied = GeometricBasisSet(basisset.basistype, x..., basisset.n, + nₘₐₓ=basisset.nₘₐₓ, nₘᵢₙ=basisset.nₘᵢₙ) + E = solve(hamiltonian, RayleighRitzMethod(varied); perturbation=perturbation).E[1] + progress && push!(history, (energy=E, parameters=copy(x))) E catch - if 0 ≤ info - push!(history, (energy=Inf, parameters=x)) - end + progress && push!(history, (energy=Inf, parameters=copy(x))) Inf end, - [basisset.r₁, basisset.rₙ], - method = optimizer, - options... + [float(basisset.r₁), float(basisset.rₙ)], + optimizer, + Optim.Options(; options...) ) - # result - res = solve(hamiltonian, GeometricBasisSet(basisset.basistype, res.minimizer..., basisset.n, nₘₐₓ=basisset.nₘₐₓ, nₘᵢₙ=basisset.nₘᵢₙ); perturbation=perturbation, info=info) - if 0 ≤ info - return ResultRayleighRitz(; - optimizer = optimizer, - initialbasisset = BasisSet(basisset.basis...), - initialgeometricbasisset = basisset, - options = options, - progress = progress, - history = history, - getfield(res, :data)..., - ) - else - return ResultRayleighRitz(; - getfield(res, :data)..., - ) - end - + optimized = GeometricBasisSet(basisset.basistype, optimization.minimizer..., basisset.n, + nₘₐₓ=basisset.nₘₐₓ, nₘᵢₙ=basisset.nₘᵢₙ) + result = solve(hamiltonian, RayleighRitzMethod(optimized); perturbation=perturbation) + return ResultOptimizationRayleighRitz(; + result, initialmethod=method, optimizer, options, progress, history, optimization, + ) end +solve(hamiltonian::Hamiltonian, basisset::Union{BasisSet,GeometricBasisSet}; perturbation=Hamiltonian()) = + solve(hamiltonian, RayleighRitzMethod(basisset); perturbation=perturbation) + +solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian()) = + solve(hamiltonian, RayleighRitzMethod(basis); perturbation=perturbation) + +optimize(hamiltonian::Hamiltonian, basisset::Union{BasisSet,GeometricBasisSet}; kwargs...) = + optimize(hamiltonian, RayleighRitzMethod(basisset); kwargs...) + +optimize(hamiltonian::Hamiltonian, basis::Basis; kwargs...) = + optimize(hamiltonian, RayleighRitzMethod(basis); kwargs...) + # matrix function matrix(basisset::BasisSet) @@ -513,50 +359,49 @@ end # docstring @doc raw""" -`solve(hamiltonian::Hamiltonian, basisset::BasisSet)` +`RayleighRitzMethod(basisset)` -This function returns the eigenvalues ``E`` and eigenvectors ``\pmb{c}`` for -```math -\pmb{H} \pmb{c} = E \pmb{S} \pmb{c}. -``` -The Hamiltonian matrix is defined as ``H_{ij} = \langle \phi_{i} | \hat{H} | \phi_{j} \rangle``. The overlap matrix is defined as ``S_{ij} = \langle \phi_{i} | \phi_{j} \rangle``. -""" solve(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4) +Configure a Rayleigh--Ritz calculation with a `BasisSet`, `GeometricBasisSet`, +or single `Basis`. +""" RayleighRitzMethod @doc raw""" -`function optimize(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...)` +`ResultRayleighRitz` -This function minimizes the energy by changing the exponents of the basis functions using Optim.jl. -```math -\frac{\partial E}{\partial a_i} = 0 -``` -""" optimize(hamiltonian::Hamiltonian, basisset::BasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) +Result of a Rayleigh--Ritz calculation. +""" ResultRayleighRitz @doc raw""" -`solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4)` +`ResultOptimizationRayleighRitz` -This a solver for 1-basis calculations. This function returns `solve(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info)`. -""" solve(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4) +Result of nonlinear basis optimization. Solver fields are forwarded from `result`. +""" ResultOptimizationRayleighRitz @doc raw""" -`optimize(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=4, optimizer=Optim.NelderMead())` +`solve(hamiltonian, method::RayleighRitzMethod; perturbation=Hamiltonian())` -This a optimizer for 1-basis calculations. This function returns `optimize(hamiltonian, BasisSet(basis); perturbation=perturbation, info=info, progress=progress, optimizer=optimizer, options...)`. -""" optimize(hamiltonian::Hamiltonian, basis::Basis; perturbation=Hamiltonian(), info=1, progress=true, optimizer=Optim.NelderMead(), options...) +Solve ``\boldsymbol{H}\boldsymbol{c}=E\boldsymbol{S}\boldsymbol{c}`` and return a +`ResultRayleighRitz`. A basis or basis set may be passed in place of `method`. +""" solve(hamiltonian::Hamiltonian, method::RayleighRitzMethod; perturbation=Hamiltonian()) @doc raw""" -`solve(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4)` +`optimize(hamiltonian, method::RayleighRitzMethod; kwargs...)` -This function is a wrapper for `solve(hamiltonian::Hamiltonian, basisset::BasisSet, ...)`. -""" solve(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4) +Minimize the ground-state energy with Optim.jl. Individual exponents are varied +for a `BasisSet`; ``r_1`` and ``r_n`` are varied for a `GeometricBasisSet`. +""" optimize @doc raw""" -`optimize(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4, optimizer=Optim.NelderMead())` +`expectation(result, operator)` -This function minimizes the energy by optimizing $r_1$ and $r_n$ using Optim.jl. -```math -\frac{\partial E}{\partial r_1} = \frac{\partial E}{\partial r_n} = 0 -``` -""" optimize(hamiltonian::Hamiltonian, basisset::GeometricBasisSet; perturbation=Hamiltonian(), info=4, progress=true, optimizer=Optim.NelderMead(), options...) +Return the expectation value of a matrix, operator, or Hamiltonian for each state. +""" expectation + +@doc raw""" +`verification(result)` + +Return the norms and absolute eigenvalue residuals for a result. +""" verification @doc raw""" `matrix(basisset::BasisSet)` diff --git a/test/FC.jl b/test/FC.jl index be62a42..0c7641a 100644 --- a/test/FC.jl +++ b/test/FC.jl @@ -72,7 +72,7 @@ ] energies = Float64[] for _ in eachindex(expected) - push!(energies, solve(hamiltonian, basisset, info=-1).E[1]) + push!(energies, solve(hamiltonian, basisset).E[1]) basisset = FC(hamiltonian, basisset) end diff --git a/test/Rayleigh-Ritz.jl b/test/Rayleigh-Ritz.jl index 6839826..28d4bed 100644 --- a/test/Rayleigh-Ritz.jl +++ b/test/Rayleigh-Ritz.jl @@ -12,7 +12,14 @@ SimpleGaussianBasis(0.444529), SimpleGaussianBasis(0.1219492), ) - res = solve(H, BS) + method = RayleighRitzMethod(BS) + res = solve(H, method) + + @test res.method === method + @test verification(res).norm ≈ ones(res.nₘₐₓ) + @test verification(res).residual ≈ zeros(res.nₘₐₓ) atol=1e-12 + @test expectation(res, H) ≈ res.E atol=1e-12 + @test occursin("# verification", string(res)) println("4π×∫|ψ(r)|²r²dr = 1") println(" i\tnumerical \tanalytical") @@ -109,15 +116,15 @@ @test uncontracted_result.E[1] ≈ -0.495010 atol=1e-6 @test uncontracted_result.E[1] < contracted_result.E[1] - @test contracted_result.expectation[radius][1] ≈ 1.500392 atol=2e-6 - @test contracted_result.expectation[radius_squared][1] ≈ 2.996124 atol=5e-6 - @test contracted_result.expectation[inverse_radius][1] ≈ 0.989206 atol=2e-6 - @test 2 * contracted_result.expectation[H.terms[1]][1] ≈ 0.988596 atol=2e-6 + @test expectation(contracted_result, radius)[1] ≈ 1.500392 atol=2e-6 + @test expectation(contracted_result, radius_squared)[1] ≈ 2.996124 atol=5e-6 + @test expectation(contracted_result, inverse_radius)[1] ≈ 0.989206 atol=2e-6 + @test 2 * expectation(contracted_result, H.terms[1])[1] ≈ 0.988596 atol=2e-6 - @test uncontracted_result.expectation[radius][1] ≈ 1.514699 atol=1e-6 - @test uncontracted_result.expectation[radius_squared][1] ≈ 3.046195 atol=1e-6 - @test uncontracted_result.expectation[inverse_radius][1] ≈ 0.977070 atol=1e-6 - @test 2 * uncontracted_result.expectation[H.terms[1]][1] ≈ 0.964119 atol=1e-6 + @test expectation(uncontracted_result, radius)[1] ≈ 1.514699 atol=1e-6 + @test expectation(uncontracted_result, radius_squared)[1] ≈ 3.046195 atol=1e-6 + @test expectation(uncontracted_result, inverse_radius)[1] ≈ 0.977070 atol=1e-6 + @test 2 * expectation(uncontracted_result, H.terms[1])[1] ≈ 0.964119 atol=1e-6 @test contracted_result.S[1,1] ≈ 1.0 atol=2e-6 @test 4π * quadgk(r -> r^2 * abs2(TwoBody.ψ(contracted_result, r)), 0, Inf)[1] ≈ 1.0 atol=1e-9 @@ -236,9 +243,18 @@ Gaussian(coefficient=-1, exponent=1), ) - custom_result = solve(custom_H, BS, info=1) - gaussian_result = solve(gaussian_H, BS, info=1) + custom_result = solve(custom_H, BS) + gaussian_result = solve(gaussian_H, BS) @test custom_result.H ≈ gaussian_result.H rtol=1e-9 @test custom_result.E ≈ gaussian_result.E rtol=1e-9 end + + @testset "Optimization result" begin + result = optimize(H, SimpleGaussianBasis(1); progress=false, iterations=5) + @test result isa ResultOptimizationRayleighRitz + @test result.result isa ResultRayleighRitz + @test result.E == result.result.E + @test isempty(result.history) + @test occursin("# optimizer", string(result)) + end end