diff --git a/.github/workflows/ci-short.yml b/.github/workflows/ci-short.yml index fce60ad..c19e9bd 100644 --- a/.github/workflows/ci-short.yml +++ b/.github/workflows/ci-short.yml @@ -16,14 +16,11 @@ jobs: - '1.7' os: - ubuntu-latest - arch: - - x64 steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} - arch: ${{ matrix.arch }} - uses: actions/cache@v1 env: cache-name: cache-artifacts diff --git a/benchmark/benchmark.jl b/benchmark/benchmark.jl index 10aa0b4..4c6061e 100644 --- a/benchmark/benchmark.jl +++ b/benchmark/benchmark.jl @@ -21,7 +21,7 @@ import QuantumOpticsBase BenchmarkTools.DEFAULT_PARAMETERS.samples = 100 # dic_commits = JSON.parsefile("test/benchmark_commits.json") -open("benchmark/benchmark_commits.json","r") do f +open("benchmark/benchmark_commits.json", "r") do f global dic_commits = JSON.parse(f) end @@ -34,18 +34,18 @@ println ## Aggregate building ### agg_dimer_small -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5.0 t1_1 = @benchmark begin HR = 0.01 shift = (2.0 * HR) - modes = [Mode(180., shift)] + modes = [Mode(180.0, shift)] mols = [ - Molecule([Mode(300., shift)], 3, [12500., 12710.]), - Molecule([Mode(300., shift)], 3, [12500., 12690.]) + Molecule([Mode(300.0, shift)], 3, [12500.0, 12710.0]), + Molecule([Mode(300.0, shift)], 3, [12500.0, 12690.0]), ] agg = Aggregate(mols) - for mol_i in 2:length(agg.molecules) + for mol_i = 2:length(agg.molecules) agg.coupling[mol_i, mol_i+1] = 200 agg.coupling[mol_i+1, mol_i] = 200 end @@ -59,7 +59,7 @@ t1_1 = @benchmark begin Ham = getAggHamiltonian(agg, aggIndices, FCFact) Ham_bath = getAggHamiltonianBath(agg) - Ham_sys = getAggHamiltonianSystem(agg; groundEnergy=false) + Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false) b = GenericBasis([aggIndLen]) b_sys = GenericBasis([size(Ham_sys, 1)]) b_bath = GenericBasis([size(Ham_bath, 1)]) @@ -77,18 +77,18 @@ t1_1 = @benchmark begin end ### agg_dimer_medium -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0 t1_2 = @benchmark begin HR = 0.01 shift = (2.0 * HR) - modes = [Mode(180., shift)] + modes = [Mode(180.0, shift)] mols = [ - Molecule([Mode(300., shift)], 7, [12500., 12710.]), - Molecule([Mode(300., shift)], 7, [12500., 12690.]) + Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]), + Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]), ] agg = Aggregate(mols) - for mol_i in 2:length(agg.molecules) + for mol_i = 2:length(agg.molecules) agg.coupling[mol_i, mol_i+1] = 200 agg.coupling[mol_i+1, mol_i] = 200 end @@ -102,7 +102,7 @@ t1_2 = @benchmark begin Ham = getAggHamiltonian(agg, aggIndices, FCFact) Ham_bath = getAggHamiltonianBath(agg) - Ham_sys = getAggHamiltonianSystem(agg; groundEnergy=false) + Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false) b = GenericBasis([aggIndLen]) b_sys = GenericBasis([size(Ham_sys, 1)]) b_bath = GenericBasis([size(Ham_bath, 1)]) @@ -120,28 +120,28 @@ t1_2 = @benchmark begin end ### agg_building -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0 t1_4 = @benchmark begin HR = 0.01 shift = (2.0 * HR) - modes = [Mode(180., shift)] + modes = [Mode(180.0, shift)] mols = [ - Molecule([Mode(300., shift)], 7, [12500., 12710.]), - Molecule([Mode(300., shift)], 7, [12500., 12690.]) + Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]), + Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]), ] agg = Aggregate(mols) end HR = 0.01 shift = (2.0 * HR) -modes = [Mode(180., shift)] +modes = [Mode(180.0, shift)] mols = [ - Molecule([Mode(300., shift)], 7, [12500., 12710.]), - Molecule([Mode(300., shift)], 7, [12500., 12690.]) + Molecule([Mode(300.0, shift)], 7, [12500.0, 12710.0]), + Molecule([Mode(300.0, shift)], 7, [12500.0, 12690.0]), ] agg = Aggregate(mols) -for mol_i in 2:length(agg.molecules) +for mol_i = 2:length(agg.molecules) agg.coupling[mol_i, mol_i+1] = 200 agg.coupling[mol_i+1, mol_i] = 200 end @@ -152,57 +152,57 @@ aggIndLen = length(aggIndices) base = GenericBasis([aggIndLen]) ### getFranckCondonFactors -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 t1_5 = @benchmark begin FCFact = getFranckCondonFactors(agg, aggIndices) end FCFact = getFranckCondonFactors(agg, aggIndices) ### getFCProd -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 100. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 100.0 t1_6 = @benchmark begin FCProd = getFCProd(agg, FCFact, aggIndices, vibindices) end FCProd = getFCProd(agg, FCFact, aggIndices, vibindices) ### getFCProd -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 t1_7 = @benchmark begin Ham = getAggHamiltonian(agg, aggIndices, FCFact) end Ham = getAggHamiltonian(agg, aggIndices, FCFact) Ham_bath = getAggHamiltonianBath(agg) -Ham_sys = getAggHamiltonianSystem(agg; groundEnergy=false) +Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false) b = GenericBasis([aggIndLen]) b_sys = GenericBasis([size(Ham_sys, 1)]) b_bath = GenericBasis([size(Ham_bath, 1)]) ### getAggHamiltonianInteraction -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 t1_8 = @benchmark begin Ham_int = getAggHamiltonianInteraction(agg, aggIndices, FCFact) end dic_commit["agg_dimer_small"] = t1_1.times # 27 dic_commit["agg_dimer_medium"] = t1_2.times # 147 -dic_commit["agg_building"] = t1_4.times -dic_commit["getFranckCondonFactors"] = t1_5.times +dic_commit["agg_building"] = t1_4.times +dic_commit["getFranckCondonFactors"] = t1_5.times dic_commit["getFCProd"] = t1_6.times -dic_commit["getAggHamiltonian"] = t1_7.times -dic_commit["getAggHamiltonianInteraction"] = t1_8.times +dic_commit["getAggHamiltonian"] = t1_7.times +dic_commit["getAggHamiltonianInteraction"] = t1_8.times ## Simulations HR = 0.01 shift = (2.0 * HR) -modes = [Mode(180., shift)] +modes = [Mode(180.0, shift)] mols = [ - Molecule([Mode(300., shift)], 3, [12500., 12710.]), - Molecule([Mode(300., shift)], 3, [12500., 12690.]) + Molecule([Mode(300.0, shift)], 3, [12500.0, 12710.0]), + Molecule([Mode(300.0, shift)], 3, [12500.0, 12690.0]), ] agg = Aggregate(mols) -for mol_i in 2:length(agg.molecules) +for mol_i = 2:length(agg.molecules) agg.coupling[mol_i, mol_i+1] = 200 agg.coupling[mol_i+1, mol_i] = 200 end @@ -217,7 +217,7 @@ Ham = getAggHamiltonian(agg, aggIndices, FCFact) basis = GenericBasis([length(aggIndices)]) Ham_bath = getAggHamiltonianBath(agg) -Ham_sys = getAggHamiltonianSystem(agg; groundEnergy=false) +Ham_sys = getAggHamiltonianSystem(agg; groundEnergy = false) b = GenericBasis([aggIndLen]) b_sys = GenericBasis([size(Ham_sys, 1)]) b_bath = GenericBasis([size(Ham_bath, 1)]) @@ -237,7 +237,7 @@ H_lambda = diagm(H_lambda) t_max = 0.01 t_count = 200 -t0 = 0. +t0 = 0.0 t_step = (t_max - t0) / (t_count) tspan = [t0:t_step:t_max;] T = 300 @@ -245,7 +245,7 @@ mu_array = [[2, 1]] W0_1 = thermal_state(T, mu_array, Ham, vibindices, aggIndices; diagonalize = false) mu_array = [[1, 2]] W0_2 = thermal_state(T, mu_array, Ham, vibindices, aggIndices; diagonalize = false) -W0 = 0.8*W0_1 + 0.2*W0_2 +W0 = 0.8 * W0_1 + 0.2 * W0_2 W0 = DenseOperator(W0.basis_l, W0.basis_r, complex(W0.data)) W0_bath = get_rho_bath(W0, agg, FCProd, aggIndices, vibindices) W0_bath = DenseOperator(W0_bath.basis_l, W0_bath.basis_r, complex(W0_bath.data)) @@ -258,19 +258,19 @@ t2_0 = @benchmark begin end ### evolutionExact -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 t2_1 = @benchmark begin op_array = evolutionExact(W0, tspan, Ham) end ### evolutionApproximate -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 2.0 t2_2 = @benchmark begin op_array = evolutionApproximate(W0, tspan, Ham) end ### schroedinger -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1.0 psi0 = randstate(basis) t2_3 = @benchmark begin Tspan, psi_t = schroedinger( @@ -284,7 +284,7 @@ t2_3 = @benchmark begin end ### schroedinger -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 t2_4 = @benchmark begin Tspan, rho_t = liouvilleVonNeumann( W0, @@ -297,12 +297,14 @@ t2_4 = @benchmark begin end ### evolutionOperatorIterator -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10. +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10.0 elLen = length(agg.molecules) -rho_t_exact = zeros(ComplexF64, length(tspan), elLen+1, elLen+1) +rho_t_exact = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1) t2_5 = @benchmark begin t_i = 0 - foreach(evolutionOperatorIterator(Ham, tspan; diagonalize = false, approximate=true)) do U_op + foreach( + evolutionOperatorIterator(Ham, tspan; diagonalize = false, approximate = true), + ) do U_op t_i = t_i + 1 # println(t_i / 150. * 100) W = U_op * W0 * U_op' @@ -312,8 +314,21 @@ t2_5 = @benchmark begin end ### master_int -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20. -p = (Ham_S, Ham_int, H_lambda, H_S, H_Sinv, W0, W0_bath, agg, FCProd, aggIndices, vibindices, ComplexF64) +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0 +p = ( + Ham_S, + Ham_int, + H_lambda, + H_S, + H_Sinv, + W0, + W0_bath, + agg, + FCProd, + aggIndices, + vibindices, + ComplexF64, +) t2_6 = @benchmark begin Tspan, W_t_master_int = master_int( W0, @@ -322,12 +337,12 @@ t2_6 = @benchmark begin Ham_int; reltol = 1.0e-4, abstol = 1.0e-4, - alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()) + alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()), ) elLen = length(agg.molecules) - rho_t_master = zeros(ComplexF64, length(tspan), elLen+1, elLen+1) - for t_i in 1:length(tspan) + rho_t_master = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1) + for t_i = 1:length(tspan) t = tspan[t_i] U_op_S = evolutionOperator(Ham_sys, t) rho = trace_bath(W_t_master_int[t_i], agg, FCProd, aggIndices, vibindices) @@ -337,8 +352,22 @@ t2_6 = @benchmark begin end ### evolutionOperatorIterator -BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20. -p = (Ham_S, Ham_int, H_lambda, H_S, H_Sinv, Ham_B, W0, W0_bath, agg, FCProd, aggIndices, vibindices, ComplexF64) +BenchmarkTools.DEFAULT_PARAMETERS.seconds = 20.0 +p = ( + Ham_S, + Ham_int, + H_lambda, + H_S, + H_Sinv, + Ham_B, + W0, + W0_bath, + agg, + FCProd, + aggIndices, + vibindices, + ComplexF64, +) t2_7 = @benchmark begin Tspan, rho_t_ansatz_int = master_ansatz( rho0, @@ -346,12 +375,12 @@ t2_7 = @benchmark begin p; reltol = 1.0e-4, abstol = 1.0e-4, - alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()) + alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5()), ) - + elLen = length(agg.molecules) - rho_t_ansatz = zeros(ComplexF64, length(tspan), elLen+1, elLen+1) - for t_i in 1:length(tspan) + rho_t_ansatz = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1) + for t_i = 1:length(tspan) t = tspan[t_i] U_op_S = evolutionOperator(Ham_sys, t) rho = U_op_S * rho_t_ansatz_int[t_i] * U_op_S' @@ -359,18 +388,18 @@ t2_7 = @benchmark begin end end -dic_commit["trace_bath"] = t2_0.times -dic_commit["evolutionExact"] = t2_1.times -dic_commit["evolutionApproximate"] = t2_2.times -dic_commit["schroedinger"] = t2_3.times -dic_commit["liouvilleVonNeumann"] = t2_4.times -dic_commit["evolutionOperatorIterator"] = t2_5.times -dic_commit["master_int"] = t2_6.times -dic_commit["master_ansatz"] = t2_7.times +dic_commit["trace_bath"] = t2_0.times +dic_commit["evolutionExact"] = t2_1.times +dic_commit["evolutionApproximate"] = t2_2.times +dic_commit["schroedinger"] = t2_3.times +dic_commit["liouvilleVonNeumann"] = t2_4.times +dic_commit["evolutionOperatorIterator"] = t2_5.times +dic_commit["master_int"] = t2_6.times +dic_commit["master_ansatz"] = t2_7.times dic_commits[commit_hash] = dic_commit json_string = JSON.json(dic_commits) -open("benchmark/benchmark_commits.json","w") do f - write(f, json_string) -end \ No newline at end of file +open("benchmark/benchmark_commits.json", "w") do f + write(f, json_string) +end diff --git a/docs/make.jl b/docs/make.jl index e419b9a..77eb62e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,9 +16,7 @@ makedocs( pages = [ "Home" => "index.md", "Documentation" => "documentation.md", - "Tutorials" => [ - "Dimer" => "tutorials/dimer.md" - ] + "Tutorials" => ["Dimer" => "tutorials/dimer.md"], ], ) diff --git a/src/OpenQuantumSystems.jl b/src/OpenQuantumSystems.jl index 18c4065..a971f72 100644 --- a/src/OpenQuantumSystems.jl +++ b/src/OpenQuantumSystems.jl @@ -178,6 +178,7 @@ export # aggregate.jl Aggregate, setupAggregate, + setupAggregate!, # evolution.jl get_tspan, @@ -252,6 +253,8 @@ export QME_sI_ansatz_upart2, # postprocessing.jl + OperatorVector, + OperatorVectorArray, operator_recast, interaction_pic_to_schroedinger_pic, schroedinger_pic_to_interaction_pic, diff --git a/src/aggregate.jl b/src/aggregate.jl index 0be969b..73ee6fb 100644 --- a/src/aggregate.jl +++ b/src/aggregate.jl @@ -4,23 +4,25 @@ abstract type AbstractAggregate end mutable struct Aggregate <: AbstractAggregate - core::Union{AggregateCore, Nothing} - tools::Union{AggregateTools, Nothing} - operators::Union{AggregateOperators, Nothing} + core::Union{AggregateCore,Nothing} + tools::Union{AggregateTools,Nothing} + operators::Union{AggregateOperators,Nothing} function Aggregate( - core::Union{AggregateCore, Nothing}, - tools::Union{AggregateTools, Nothing}, - operators::Union{AggregateOperators, Nothing} + core::Union{AggregateCore,Nothing}, + tools::Union{AggregateTools,Nothing}, + operators::Union{AggregateOperators,Nothing}, )::Aggregate new(core, tools, operators) end end -Aggregate(aggCore::AggregateCore)::Aggregate = Aggregate(aggCore, nothing, nothing)::Aggregate +Aggregate(aggCore::AggregateCore)::Aggregate = + Aggregate(aggCore, nothing, nothing)::Aggregate AggregateTools(agg::Aggregate)::AggregateTools = AggregateTools(agg.core) -AggregateOperators(agg::Aggregate)::AggregateOperators = AggregateOperators(agg.core, agg.tools) +AggregateOperators(agg::Aggregate)::AggregateOperators = + AggregateOperators(agg.core, agg.tools) """ setupAggregate(agg; groundEnergy=true, verbose=false) @@ -31,15 +33,15 @@ Generate all basic data from the [`Aggregate`](@ref). Returns """ function setupAggregate(aggCore::AggregateCore; groundEnergy::Bool = true)::Aggregate aggTools = AggregateTools(aggCore) - aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy=groundEnergy) + aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy = groundEnergy) return Aggregate(aggCore, aggTools, aggOperators) end function setupAggregate!(agg::Aggregate; groundEnergy::Bool = true)::Aggregate agg.tools = AggregateTools(agg.core) - agg.operators = AggregateOperators(agg.core, agg.tools; groundEnergy=groundEnergy) + agg.operators = AggregateOperators(agg.core, agg.tools; groundEnergy = groundEnergy) return agg end -Base.:(==)(x::Aggregate, y::Aggregate) = +Base.:(==)(x::Aggregate, y::Aggregate) = x.core == y.core && x.tools == y.tools && x.operators == y.operators diff --git a/src/aggregateCore.jl b/src/aggregateCore.jl index c5fe656..7cf8a82 100644 --- a/src/aggregateCore.jl +++ b/src/aggregateCore.jl @@ -11,14 +11,15 @@ Mutable stuct which purpose is to store important information about the whole ag * `molecules`: Vector of molecules ([`Molecule`](@ref)). * `coupling`: Matrix of couplings ``J_{nm}`` between molecules. """ -struct AggregateCore{T<:Integer,C1<:ComputableType,C2<:ComputableType} <: AbstractAggregateCore +struct AggregateCore{T<:Integer,C1<:ComputableType,C2<:ComputableType} <: + AbstractAggregateCore molecules::Vector{Molecule{T,C1,C2}} coupling::Matrix{C1} molCount::Int64 function AggregateCore{T,C1,C2}( molecules::Vector{Molecule{T,C1,C2}}, coupling::Matrix{C1}, - molCount::Int64 + molCount::Int64, ) where {T<:Integer,C1<:ComputableType,C2<:ComputableType} new(molecules, coupling, molCount) end @@ -27,10 +28,14 @@ end AggregateCore(molecules::Vector{Molecule{T,C1,C2}}, coupling::Matrix{C1}) where {T,C1,C2} = AggregateCore{T,C1,C2}(molecules, coupling, length(molecules)) AggregateCore(molecules::Vector{Molecule{T,C1,C2}}) where {T,C1,C2} = - AggregateCore{T,C1,C2}(molecules, zeros(C1, (length(molecules) + 1, length(molecules) + 1)), length(molecules)) + AggregateCore{T,C1,C2}( + molecules, + zeros(C1, (length(molecules) + 1, length(molecules) + 1)), + length(molecules), + ) -Base.:(==)(x::AggregateCore, y::AggregateCore) = +Base.:(==)(x::AggregateCore, y::AggregateCore) = x.molecules == y.molecules && x.coupling == y.coupling && x.molCount == y.molCount """ @@ -110,4 +115,4 @@ function elIndOrder(elInd::Vector{T})::T where {T<:Integer} ind += (elInd[i] - 1) * i end return ind -end \ No newline at end of file +end diff --git a/src/aggregateOperators.jl b/src/aggregateOperators.jl index d06faea..15270c1 100644 --- a/src/aggregateOperators.jl +++ b/src/aggregateOperators.jl @@ -13,7 +13,7 @@ function getAggHamSystemSmall( aggCore::AggregateCore, aggTools::AggregateTools; groundEnergy::Bool = true, -) +) Ham_sys = zeros(Float64, (aggCore.molCount + 1, aggCore.molCount + 1)) agg_shifts = getShifts(aggCore) @@ -21,9 +21,10 @@ function getAggHamSystemSmall( reorganisation_energies = [] for mol_i = 1:aggCore.molCount mol = aggCore.molecules[mol_i] - reorganisation_energy = 0. + reorganisation_energy = 0.0 for mode_i = 1:length(mol.modes) - reorganisation_energy += agg_frequencies[mol_i][mode_i] * agg_shifts[mol_i][mode_i]^2 / 2. + reorganisation_energy += + agg_frequencies[mol_i][mode_i] * agg_shifts[mol_i][mode_i]^2 / 2.0 end push!(reorganisation_energies, reorganisation_energy) end @@ -41,13 +42,13 @@ function getAggHamSystemSmall( if ind > 1 E_state += reorganisation_energies[ind-1] end - Ham_sys[ind, ind] = E_state + Ham_sys[ind, ind] = E_state end Ham_sys[:, :] += aggCore.coupling[:, :] if !groundEnergy E0 = Ham_sys[1, 1] - for i = 1:aggTool.bSystemSize + for i = 1:aggTools.bSystemSize Ham_sys[i, i] -= E0 end end @@ -67,7 +68,7 @@ function getAggHamSystemBig( aggCore::AggregateCore, aggTools::AggregateTools; groundEnergy::Bool = true, -) +) Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy = groundEnergy) Ham = zeros(Float64, (aggTools.bSize, aggTools.bSize)) for I = 1:aggTools.bSize @@ -99,12 +100,12 @@ Get Hamiltonian of the bath, ``H_B`` only for ``\\mathcal{H}_B``. function getAggHamBathSmall( aggCore::AggregateCore, aggTools::AggregateTools; - groundEnergy::Bool = true + groundEnergy::Bool = true, ) Ham_bath = zeros(Float64, (aggTools.bBathSize, aggTools.bBathSize)) elind = fill(1, (aggCore.molCount + 1)) - Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy=true) + Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy = true) for I = 1:aggTools.bBathSize vibind = aggTools.vibIndices[I] Ham_bath[I, I] = getAggStateEnergy(aggCore, elind, vibind) - Ham_sys.data[1, 1] @@ -128,7 +129,7 @@ Get Hamiltonian of the bath, ``H_B`` only for ``\\mathcal{H}_S \\otimes \\mathca function getAggHamBathBig( aggCore::AggregateCore, aggTools::AggregateTools; - groundEnergy::Bool = true + groundEnergy::Bool = true, ) Ham_bath = getAggHamBathSmall(aggCore, aggTools) @@ -163,8 +164,8 @@ function getAggHamSystemBath( aggTools::AggregateTools; groundEnergy::Bool = true, ) - Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy=true) - Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy=true) + Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy = true) + Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy = true) Ham_0 = Ham_B + Ham_S if !groundEnergy E0 = Ham_0.data[1, 1] @@ -188,10 +189,7 @@ Get interation Hamiltonian of the [`Aggregate`](@ref). * `aggTools.indices`: Aggregate indices generated by [`getIndices`](@ref). * `franckCondonFactors`: Franck-Condon factors generated by [`getFranckCondonFactors`](@ref). """ -function getAggHamInteraction( - aggCore::AggregateCore, - aggTools::AggregateTools -) +function getAggHamInteraction(aggCore::AggregateCore, aggTools::AggregateTools) Ham_I = zeros(Float64, (aggTools.bSize, aggTools.bSize)) agg_shifts = getShifts(aggCore) agg_frequencies = getShifts(aggCore) @@ -215,25 +213,28 @@ function getAggHamInteraction( if elOrder2 == 1 || elOrder1 != elOrder2 continue end - diff_num = 0; mol_j = 0; mode_j = 0 + diff_num = 0 + mol_j = 0 + mode_j = 0 for mol_i = 1:aggCore.molCount mol = aggCore.molecules[mol_i] for mode_i = 1:length(mol.modes) diff_vib = abs(vibind1[mol_i][mode_i] - vibind2[mol_i][mode_i]) diff_num += diff_vib if diff_vib == 1 && diff_num == 1 - mol_j = mol_i; mode_j = mode_i; + mol_j = mol_i + mode_j = mode_i end end end - if diff_num != 1 || mol_j != elOrder1-1 + if diff_num != 1 || mol_j != elOrder1 - 1 continue end vib_n = vibind1[mol_j][mode_j] vib_m = vibind2[mol_j][mode_j] coeff = agg_coeffs[mol_j][mode_j] - Ham_I[I, J] = - coeff * sqrt(min(vib_n, vib_m)) + Ham_I[I, J] = -coeff * sqrt(min(vib_n, vib_m)) end end @@ -256,11 +257,11 @@ Get Hamiltonian of the [`Aggregate`](@ref). function getAggHamiltonian( aggCore::AggregateCore, aggTools::AggregateTools; - groundEnergy::Bool = true + groundEnergy::Bool = true, ) Ham_I = getAggHamInteraction(aggCore, aggTools) Ham_0 = getAggHamSystemBath(aggCore, aggTools; groundEnergy = groundEnergy) - return Ham_0 + Ham_I + return Ham_0 + Ham_I end ### Sparse versions @@ -334,7 +335,7 @@ nothing =# # TODO: add docs and proper types -struct AggregateOperators{O_sys, O_bath, O} <: AbstractAggregateOperators +struct AggregateOperators{O_sys,O_bath,O} <: AbstractAggregateOperators Ham_sys::O_sys Ham_bath::O_bath Ham_S::O @@ -343,7 +344,7 @@ struct AggregateOperators{O_sys, O_bath, O} <: AbstractAggregateOperators Ham_I::O Ham::O groundEnergy::Bool - function AggregateOperators{O_sys, O_bath, O}( + function AggregateOperators{O_sys,O_bath,O}( Ham_sys::O_sys, Ham_bath::O_bath, Ham_S::O, @@ -352,42 +353,47 @@ struct AggregateOperators{O_sys, O_bath, O} <: AbstractAggregateOperators Ham_I::O, Ham::O, groundEnergy::Bool, - ) where {O_sys<:Operator, O_bath<:Operator, O<:Operator} - new( - Ham_sys, Ham_bath, - Ham_S, Ham_B, Ham_0, Ham_I, Ham, - groundEnergy - ) + ) where {O_sys<:Operator,O_bath<:Operator,O<:Operator} + new(Ham_sys, Ham_bath, Ham_S, Ham_B, Ham_0, Ham_I, Ham, groundEnergy) end end function AggregateOperators( - aggCore::AggregateCore, + aggCore::AggregateCore, aggTools::AggregateTools; - groundEnergy::Bool = true -) + groundEnergy::Bool = true, +) Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy = true) Ham_bath = getAggHamBathSmall(aggCore, aggTools; groundEnergy = true) - Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy=groundEnergy) - Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy=groundEnergy) - Ham_0 = getAggHamSystemBath(aggCore, aggTools; groundEnergy=groundEnergy) + Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy = groundEnergy) + Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy = groundEnergy) + Ham_0 = getAggHamSystemBath(aggCore, aggTools; groundEnergy = groundEnergy) Ham_I = getAggHamInteraction(aggCore, aggTools) - Ham = getAggHamiltonian(aggCore, aggTools; groundEnergy=groundEnergy) + Ham = getAggHamiltonian(aggCore, aggTools; groundEnergy = groundEnergy) O_sys = typeof(Ham_sys) O_bath = typeof(Ham_bath) O = typeof(Ham_S) - return AggregateOperators{O_sys, O_bath, O}( - Ham_sys, Ham_bath, - Ham_S, Ham_B, Ham_0, Ham_I, Ham, - groundEnergy + return AggregateOperators{O_sys,O_bath,O}( + Ham_sys, + Ham_bath, + Ham_S, + Ham_B, + Ham_0, + Ham_I, + Ham, + groundEnergy, ) end -Base.:(==)(x::AggregateOperators, y::AggregateOperators) = - x.Ham_sys == y.Ham_sys && x.Ham_bath == y.Ham_bath && - x.Ham_S == y.Ham_S && x.Ham_B == y.Ham_B && x.Ham_0 == y.Ham_0 && - x.Ham_I == y.Ham_I && x.Ham == y.Ham && - x.groundEnergy == y.groundEnergy \ No newline at end of file +Base.:(==)(x::AggregateOperators, y::AggregateOperators) = + x.Ham_sys == y.Ham_sys && + x.Ham_bath == y.Ham_bath && + x.Ham_S == y.Ham_S && + x.Ham_B == y.Ham_B && + x.Ham_0 == y.Ham_0 && + x.Ham_I == y.Ham_I && + x.Ham == y.Ham && + x.groundEnergy == y.groundEnergy diff --git a/src/aggregateTools.jl b/src/aggregateTools.jl index 63e599d..558e9d4 100644 --- a/src/aggregateTools.jl +++ b/src/aggregateTools.jl @@ -1,13 +1,13 @@ abstract type AbstractAggregateTools end -const ElIndices = Vector{Vector{T}} where T <: Integer -const VibIndices = Vector{Vector{Vector{Int64}}} where T <: Integer +const ElIndices = Vector{Vector{T}} where {T<:Integer} +const VibIndices = Vector{Vector{Vector{Int64}}} where {T<:Integer} # TODO: be more specific for T const Indices = Vector{Vector{Vector{T} where T}} -const IndicesMap = Vector{Vector{T}} where T <: Integer -const FCfactorsT = Matrix{T} where T <: AbstractFloat -const FCproductT = Matrix{T} where T <: AbstractFloat +const IndicesMap = Vector{Vector{T}} where {T<:Integer} +const FCfactorsT = Matrix{T} where {T<:AbstractFloat} +const FCproductT = Matrix{T} where {T<:AbstractFloat} """ vibrationalIndices(aggCore) @@ -115,10 +115,7 @@ Get Frack-Condon factors of the [`Aggregate`](@ref) in a form of matrix. * `aggCore`: Instance of [`Aggregate`](@ref). * `aggIndices`: Aggregate indices generated by [`getIndices`](@ref). """ -function getFranckCondonFactors( - aggCore::AggregateCore, - aggIndices::Indices -)::FCfactorsT +function getFranckCondonFactors(aggCore::AggregateCore, aggIndices::Indices)::FCfactorsT aggIndLen = length(aggIndices) molLen = length(aggCore.molecules) FC = zeros(Float64, (aggIndLen, aggIndLen)) @@ -144,12 +141,9 @@ function getFranckCondonFactors( return FC end -getFranckCondonFactors( - aggCore::AggregateCore - )::FCfactorsT = ( - aggIndices = getIndices(aggCore); - getFranckCondonFactors(aggCore, aggIndices) -) +getFranckCondonFactors(aggCore::AggregateCore)::FCfactorsT = + (aggIndices = getIndices(aggCore); + getFranckCondonFactors(aggCore, aggIndices)) # TODO: resurrect sparse version #= @@ -209,11 +203,11 @@ Get product of Franck-Condon factors. This way the trace over bath will be faste """ function getFCproduct( - aggCore::AggregateCore, - indices::Indices, - indicesMap::IndicesMap, - FCfactors::FCfactorsT - )::FCproductT + aggCore::AggregateCore, + indices::Indices, + indicesMap::IndicesMap, + FCfactors::FCfactorsT, +)::FCproductT elLen = length(aggCore.molecules) aggIndLen = length(indices) vibLen = length(indicesMap[2]) @@ -274,9 +268,18 @@ struct AggregateTools <: AbstractAggregateTools bSize::Int64, )::AggregateTools new( - elIndices, vibIndices, indices, indicesMap, - FCfactors, FCproduct, - basisSystem, basisBath, basis, bSystemSize, bBathSize, bSize + elIndices, + vibIndices, + indices, + indicesMap, + FCfactors, + FCproduct, + basisSystem, + basisBath, + basis, + bSystemSize, + bBathSize, + bSize, ) end end @@ -288,7 +291,7 @@ function AggregateTools(aggCore::AggregateCore)::AggregateTools indicesMap = getIndicesMap(aggCore, indices) FCfactors = getFranckCondonFactors(aggCore, indices) FCproduct = getFCproduct(aggCore, indices, indicesMap, FCfactors) - + bSystemSize = length(elIndices) bBathSize = length(vibIndices) bSize = length(indices) @@ -296,20 +299,36 @@ function AggregateTools(aggCore::AggregateCore)::AggregateTools basisSystem = GenericBasis([bSystemSize]) basisBath = GenericBasis([bBathSize]) basis = GenericBasis([bSize]) - + return AggregateTools( - elIndices, vibIndices, indices, indicesMap, - FCfactors, FCproduct, - basisSystem, basisBath, basis, bSystemSize, bBathSize, bSize + elIndices, + vibIndices, + indices, + indicesMap, + FCfactors, + FCproduct, + basisSystem, + basisBath, + basis, + bSystemSize, + bBathSize, + bSize, ) end -Base.:(==)(x::AggregateTools, y::AggregateTools) = - x.elIndices == y.elIndices && x.vibIndices == y.vibIndices && x.indices == y.indices && - x.FCfactors == y.FCfactors && x.FCproduct == y.FCproduct && - x.bSystemSize == y.bSystemSize && x.bBathSize == y.bBathSize && x.bSize == y.bSize && - x.basisSystem == y.basisSystem && x.basisBath == y.basisBath && x.basis == y.basis +Base.:(==)(x::AggregateTools, y::AggregateTools) = + x.elIndices == y.elIndices && + x.vibIndices == y.vibIndices && + x.indices == y.indices && + x.FCfactors == y.FCfactors && + x.FCproduct == y.FCproduct && + x.bSystemSize == y.bSystemSize && + x.bBathSize == y.bBathSize && + x.bSize == y.bSize && + x.basisSystem == y.basisSystem && + x.basisBath == y.basisBath && + x.basis == y.basis """ take_el_part(A, a, b, indicesMap) @@ -324,4 +343,4 @@ function take_el_part(A::Array, a, b, indicesMap) b2 = indicesMap[b][end] return A[a1:a2, b1:b2] -end \ No newline at end of file +end diff --git a/src/evolution.jl b/src/evolution.jl index 0db7c17..bad8035 100644 --- a/src/evolution.jl +++ b/src/evolution.jl @@ -147,11 +147,11 @@ function evolutionOperatorIterator end U_diagonal = map(lambda -> exp(-1im * lambda * tspan[1]), Ham_lambda) U = S * diagm(U_diagonal) * inv(S) U_op = DenseOperator(basis, basis, U) - elseif !approximate + elseif approximate U_op = evolutionOperator(Hamiltonian, tspan[1]) + U_step = evolutionOperator(Hamiltonian, t_step) else U_op = evolutionOperator(Hamiltonian, tspan[1]) - U_step = evolutionOperator(Hamiltonian, t_step) end @yield U_op @@ -160,10 +160,10 @@ function evolutionOperatorIterator end if diagonalize U_diagonal .= map(lambda -> exp(-1im * lambda * t), Ham_lambda) U_op.data .= S * diagm(U_diagonal) * inv(S) - elseif !approximate - U_op = evolutionOperator(Hamiltonian, t) - else + elseif approximate U_op = U_step * U_op + else + U_op = evolutionOperator(Hamiltonian, t) end @yield U_op end @@ -199,11 +199,11 @@ function evolutionSuperOperatorIterator end U_diagonal = map(lambda -> exp(-1im * lambda * tspan[1]), Ham_lambda) U = S * diagm(U_diagonal) * Sinv U_op = DenseOperator(basis, basis, U) - elseif !approximate + elseif approximate U_op = evolutionOperator(Hamiltonian, tspan[1]) + U_step = evolutionOperator(Hamiltonian, t_step) else U_op = evolutionOperator(Hamiltonian, tspan[1]) - U_step = evolutionOperator(Hamiltonian, t_step) end U_supop = spre(U_op) * spost(U_op') @yield U_supop @@ -214,11 +214,11 @@ function evolutionSuperOperatorIterator end U_diagonal .= map(lambda -> exp(-1im * lambda * tspan[t_i]), Ham_lambda) U_op.data .= S * diagm(U_diagonal) * Sinv U_supop.data .= (spre(U_op)).data * (spost(U_op')).data - elseif !approximate - U_op = evolutionOperator(Hamiltonian, t) + elseif approximate + U_op = U_step * U_op U_supop = spre(U_op) * spost(U_op') else - U_op = U_step * U_op + U_op = evolutionOperator(Hamiltonian, t) U_supop = spre(U_op) * spost(U_op') end @yield U_supop @@ -585,48 +585,50 @@ function evolution_approximate( return tspan, rho_t end -function evolutionOperatorExp( - Ham::Array, t::AbstractFloat, n::Integer -)::Array +#= +function evolutionOperatorExp(Ham::Array, t::AbstractFloat, n::Integer)::Array # c = ones(size(Ham)) c = one(Ham) Len = size(Ham, 1) - s = c + s = c if n > 0 - for k in 1:n - c = c * (-1.0im/k * Ham * t) + for k = 1:n + c = c * (-1.0im / k * Ham * t) s = s + c # s[:, :] = s + (-1.0im * Ham * t)^k/factorial(big(k)) end end - s /= abs(det(s))^(1.0/Len) + s /= abs(det(s))^(1.0 / Len) s end function evolutionOperatorExp( - Ham::T, t::AbstractFloat, n::Integer + Ham::T, + t::AbstractFloat, + n::Integer, ) where {B<:Basis,T<:Operator{B,B}} data = evolutionOperatorExp(Ham.data, t, n) DenseOperator(Ham.basis_l, Ham.basis_r, data) end +=# function evolution_el_part( - Ham::Array, - t::AbstractFloat, - a::Integer, - b::Integer, - indicesMap::IndicesMap + Ham::Array, + t::AbstractFloat, + a::Integer, + b::Integer, + indicesMap::IndicesMap, ) data = take_el_part(Ham, a, b, indicesMap) return evolutionOperator(data, t) end function evolution_el_part( - Ham::T, - t::AbstractFloat, - a::Integer, - b::Integer, - indicesMap::IndicesMap + Ham::T, + t::AbstractFloat, + a::Integer, + b::Integer, + indicesMap::IndicesMap, ) where {B<:Basis,T<:Operator{B,B}} data = evolution_el_part(Ham.data, t, a, b, indicesMap) b = GenericBasis([size(data, 1)]) @@ -636,12 +638,12 @@ end function Evolution_SI_exact( W0::T, tspan::Array, - agg::Aggregate + agg::Aggregate, ) where {B<:Basis,T<:Operator{B,B}} W_t_exact = zeros(ComplexF64, length(tspan), agg.tools.bSize, agg.tools.bSize) Ham = agg.operators.Ham Ham_0 = agg.operators.Ham_0 - for t_i in 1:length(tspan) + for t_i = 1:length(tspan) t = tspan[t_i] U_op = evolutionOperator(Ham, t) W = U_op * W0 * U_op' @@ -655,14 +657,14 @@ end function Evolution_sI_exact( W0::T, tspan::Array, - agg::Aggregate + agg::Aggregate, ) where {B<:Basis,T<:Operator{B,B}} elLen = agg.core.molCount - rho_t_exact = zeros(ComplexF64, length(tspan), elLen+1, elLen+1) + rho_t_exact = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1) Ham = agg.operators.Ham Ham_0 = agg.operators.Ham_0 - for t_i in 1:length(tspan) + for t_i = 1:length(tspan) t = tspan[t_i] U_op = evolutionOperator(Ham, t) W = U_op * W0 * U_op' @@ -676,12 +678,12 @@ end function Evolution_SS_exact( W0::T, tspan::Array, - agg::Aggregate + agg::Aggregate, ) where {B<:Basis,T<:Operator{B,B}} W_t_exact = zeros(ComplexF64, length(tspan), agg.tools.bSize, agg.tools.bSize) Ham = agg.operators.Ham Ham_0 = agg.operators.Ham_0 - for t_i in 1:length(tspan) + for t_i = 1:length(tspan) t = tspan[t_i] U_op = evolutionOperator(Ham, t) W = U_op * W0 * U_op' @@ -693,18 +695,18 @@ end function Evolution_sS_exact( W0::T, tspan::Array, - agg::Aggregate + agg::Aggregate, ) where {B<:Basis,T<:Operator{B,B}} elLen = agg.core.molCount - rho_t_exact = zeros(ComplexF64, length(tspan), elLen+1, elLen+1) + rho_t_exact = zeros(ComplexF64, length(tspan), elLen + 1, elLen + 1) Ham = agg.operators.Ham Ham_0 = agg.operators.Ham_0 - for t_i in 1:length(tspan) + for t_i = 1:length(tspan) t = tspan[t_i] U_op = evolutionOperator(Ham, t) W = U_op * W0 * U_op' rho_t_exact[t_i, :, :] = trace_bath(W.data[:, :], agg.core, agg.tools) end return tspan, rho_t_exact -end \ No newline at end of file +end diff --git a/src/initial_state.jl b/src/initial_state.jl index 8b0dbdb..6cd3608 100644 --- a/src/initial_state.jl +++ b/src/initial_state.jl @@ -81,51 +81,17 @@ function thermal_state( return W0 end -function thermal_state_old( +thermal_state(T, mu_array, agg::Aggregate; boltzmann_const::Float64 = 0.69503476, diagonalize::Bool = false, diagonal = false +) = thermal_state( T, mu_array, - aggCore::AggregateCore, - aggTools::AggregateTools, - aggOperators::AggregateOperators; - boltzmann_const::Float64 = 0.69503476, - diagonalize::Bool = false, - diagonal = false + agg.core, + agg.tools, + agg.operators; + boltzmann_const=boltzmann_const, + diagonalize=diagonalize, + diagonal=diagonal ) - vibindices = aggTools.indicesMap - aggIndices = aggTools.indices - Ham = aggOperators.Ham - aggIndsLen = length(aggIndices) - data = -Ham.data / (T * boltzmann_const) - if diagonal - data = Diagonal(data) - end - if !diagonalize - data = exp(data) - else - H_lambda, H_S = eigen(data) - H_lambda = diagm(H_lambda) - data = H_S * exp(-H_lambda / (T * boltzmann_const)) * inv(H_S) - end - W0 = DenseOperator(Ham.basis_r, Ham.basis_l, data) - - for I = 1:aggIndsLen - elind1, vibind1 = aggIndices[I] - elOrder1 = OpenQuantumSystems.elIndOrder(elind1) - if elind1 in mu_array - continue - end - for J = 1:aggIndsLen - elind2, vibind2 = aggIndices[J] - elOrder2 = OpenQuantumSystems.elIndOrder(elind2) - if elind2 in mu_array - continue - end - W0.data[I, J] = 0.0 - end - end - normalize!(W0) - return W0 -end """ thermal_state_composite(T, mu_weighted, Ham, aggIndices; @@ -190,49 +156,17 @@ function thermal_state_composite( return W0 end -function thermal_state_composite_old( +thermal_state_composite(T, mu_weighted, agg::Aggregate; boltzmann_const::Float64 = 0.69503476, diagonalize::Bool = false, diagonal = false +) = thermal_state_composite( T, mu_weighted, - aggCore::AggregateCore, - aggTools::AggregateTools, - aggOperators::AggregateOperators; - boltzmann_const::Float64 = 0.69503476, - diagonalize::Bool = false, - diagonal = false, + agg.core, + agg.tools, + agg.operators; + boltzmann_const=boltzmann_const, + diagonalize=diagonalize, + diagonal=diagonal ) - vibindices = aggTools.indicesMap - aggIndices = aggTools.indices - Ham = aggOperators.Ham - aggIndsLen = length(aggIndices) - data = -Ham.data / (T * boltzmann_const) - if diagonal - data = Diagonal(data) - end - if !diagonalize - data = exp(data) - else - H_lambda, H_S = eigen(data) - H_lambda = diagm(H_lambda) - data = H_S * exp(-H_lambda / (T * boltzmann_const)) * inv(H_S) - end - W0 = DenseOperator(Ham.basis_r, Ham.basis_l, zero(data)) - - for I = 1:aggIndsLen - elind1, vibind1 = aggIndices[I] - elOrder1 = OpenQuantumSystems.elIndOrder(elind1) - mu_w = mu_weighted[elOrder1] - for J = 1:aggIndsLen - elind2, vibind2 = aggIndices[J] - elOrder2 = OpenQuantumSystems.elIndOrder(elind2) - if elOrder1 != elOrder2 - continue - end - W0.data[I, J] = data[I, J] * mu_w - end - end - normalize!(W0) - return W0 -end function ultrafast_laser_excitation(T::AbstractFloat, weights::Array, agg::Aggregate; diagonalize = true) molCount = agg.core.molCount diff --git a/src/metrics.jl b/src/metrics.jl index cae0eb4..d9f7e54 100644 --- a/src/metrics.jl +++ b/src/metrics.jl @@ -17,6 +17,7 @@ function tracenorm(rho::DenseSuperOpType) ishermitian(rho) ? tracenorm_h(rho) : tracenorm_nh(rho) end +#= function tracenorm!(rho::DenseSuperOpType) if ishermitian(rho) rho.data ./= tracenorm_h(rho) @@ -24,6 +25,7 @@ function tracenorm!(rho::DenseSuperOpType) rho.data ./= tracenorm_nh(rho) end end +=# """ tracenorm_h(rho) diff --git a/src/postprocessing.jl b/src/postprocessing.jl index dc94b90..4d14e7f 100644 --- a/src/postprocessing.jl +++ b/src/postprocessing.jl @@ -1,6 +1,6 @@ const OperatorVector = Vector{Operator{GenericBasis{Vector{Int64}}, GenericBasis{Vector{Int64}}, Matrix{ComplexF64}}} -# const OperatorVectorArray = Array{ComplexF64, 3} +const OperatorVectorArray = Array{ComplexF64, 3} function operator_recast(rho::OperatorVector)::Array N = length(rho) @@ -12,13 +12,13 @@ function operator_recast(rho::OperatorVector)::Array return rho_array end -function operator_recast(rho::Array)::Array +function operator_recast(rho::OperatorVectorArray)::Array return rho end ###### -function interaction_pic_to_schroedinger_pic(rho_int::Array, tspan::Array, agg::Aggregate) +function interaction_pic_to_schroedinger_pic(rho_int::OperatorVectorArray, tspan::Array, agg::Aggregate) rho_sch = deepcopy(rho_int) Ham_sys = agg.operators.Ham_sys for t_i in 1:length(tspan) @@ -30,11 +30,11 @@ function interaction_pic_to_schroedinger_pic(rho_int::Array, tspan::Array, agg:: end function interaction_pic_to_schroedinger_pic(rho_int::OperatorVector, tspan::Array, agg::Aggregate) - rho_array = operator_recast(rho) + rho_array = operator_recast(rho_int) return interaction_pic_to_schroedinger_pic(rho_array, tspan, agg) end -function schroedinger_pic_to_interaction_pic(rho_sch::Array, tspan::Array, agg::Aggregate) +function schroedinger_pic_to_interaction_pic(rho_sch::OperatorVectorArray, tspan::Array, agg::Aggregate) rho_int = deepcopy(rho_sch) Ham_sys = agg.operators.Ham_sys for t_i in 1:length(tspan) @@ -45,14 +45,14 @@ function schroedinger_pic_to_interaction_pic(rho_sch::Array, tspan::Array, agg:: return rho_int end -function schroedinger_pic_to_interaction_pic(rho_int::OperatorVector, tspan::Array, agg::Aggregate) - rho_array = operator_recast(rho) +function schroedinger_pic_to_interaction_pic(rho_sch::OperatorVector, tspan::Array, agg::Aggregate) + rho_array = operator_recast(rho_sch) return schroedinger_pic_to_interaction_pic(rho_array, tspan, agg) end ##### -function local_st_to_exciton_st(rho_local::Array, agg::Aggregate) +function local_st_to_exciton_st(rho_local::OperatorVectorArray, agg::Aggregate) N = size(rho_local, 1) rho_exciton = deepcopy(rho_local) Ham_sys = agg.operators.Ham_sys @@ -66,8 +66,13 @@ function local_st_to_exciton_st(rho_local::Array, agg::Aggregate) return rho_exciton end -function exciton_st_to_local_st(rho_exciton::Array, agg::Aggregate) - N = size(rho_exciton, N) +function local_st_to_exciton_st(rho_local::OperatorVector, agg::Aggregate) + rho_array = operator_recast(rho_local) + return local_st_to_exciton_st(rho_array, agg) +end + +function exciton_st_to_local_st(rho_exciton::OperatorVectorArray, agg::Aggregate) + N = size(rho_exciton, 1) rho_local = deepcopy(rho_exciton) Ham_sys = agg.operators.Ham_sys Ham_sys_lambda, Ham_sys_S = eigen(Ham_sys.data) @@ -80,3 +85,7 @@ function exciton_st_to_local_st(rho_exciton::Array, agg::Aggregate) return rho_local end +function exciton_st_to_local_st(rho_exciton::OperatorVector, agg::Aggregate) + rho_array = operator_recast(rho_exciton) + return exciton_st_to_local_st(rho_array, agg) +end diff --git a/src/scoring.jl b/src/scoring.jl index c0af617..501f1bf 100644 --- a/src/scoring.jl +++ b/src/scoring.jl @@ -1,34 +1,29 @@ ######## -function get_rmse_in_time(operator_vec1::OperatorVector, operator_vec2::OperatorVector) - t_count = length(operator_vec1) +function get_rmse_in_time(operator_data_vec1::OperatorVectorArray, operator_data_vec2::OperatorVectorArray) + t_count, M, K = size(operator_data_vec1) c = 0 for t_i in 1:t_count - op1 = operator_vec1[t_i] - op2 = operator_vec2[t_i] - c += norm(op1.data - op2.data)^2 + c += norm(operator_data_vec1[t_i, :, :] - operator_data_vec2[t_i, :, :])^2 end return sqrt(c / t_count) end -function get_rmse_in_time(operator_vec1::OperatorVector, operator_data_vec2::Array) - t_count = length(operator_vec1) - c = 0 - for t_i in 1:t_count - op1 = operator_vec1[t_i] - c += norm(op1.data - operator_data_vec2[t_i, :, :])^2 - end - return sqrt(c / t_count) +function get_rmse_in_time(operator_data_vec1::OperatorVectorArray, operator_vec2::OperatorVector) + operator_data_vec2 = operator_recast(operator_vec2) + return get_rmse_in_time(operator_data_vec1, operator_data_vec2) end -function get_rmse_in_time(operator_data_vec1::Array, operator_data_vec2::Array) - t_count = length(operator_vec1) - c = 0 - for t_i in 1:t_count - c += norm(operator_data_vec1[t_i, :, :] - operator_data_vec2[t_i, :, :])^2 - end - return sqrt(c / t_count) +function get_rmse_in_time(operator_vec1::OperatorVector, operator_data_vec2::OperatorVectorArray) + operator_data_vec1 = operator_recast(operator_vec1) + return get_rmse_in_time(operator_data_vec1, operator_data_vec2) +end + +function get_rmse_in_time(operator_vec1::OperatorVector, operator_vec2::OperatorVector) + operator_data_vec1 = operator_recast(operator_vec1) + operator_data_vec2 = operator_recast(operator_vec2) + return get_rmse_in_time(operator_data_vec1, operator_data_vec2) end ####### @@ -46,12 +41,12 @@ function compare_rho(rho::Array, rho_ref::Array; smooth_const=1e-9) return rho_sum end -function compare_rho(rho::Array, rho_ref::OperatorVector; smooth_const=1e-9) +function compare_rho(rho::OperatorVectorArray, rho_ref::OperatorVector; smooth_const=1e-9) rho_ref_array = operator_recast(rho_ref) return compare_rho(rho, rho_ref_array, smooth_const=smooth_const) end -function compare_rho(rho::OperatorVector, rho_ref::Array; smooth_const=1e-9) +function compare_rho(rho::OperatorVector, rho_ref::OperatorVectorArray; smooth_const=1e-9) rho_array = operator_recast(rho) return compare_rho(rho_array, rho_ref, smooth_const=smooth_const) end @@ -64,7 +59,7 @@ end ###### -function compare_rho_in_time(rho::Array, rho_ref::Array; smooth_const=1e-9) +function compare_rho_in_time(rho::OperatorVectorArray, rho_ref::OperatorVectorArray; smooth_const=1e-9) N, M, K = size(rho) rho_rel = zeros(Float64, N, M, K) @@ -76,12 +71,12 @@ function compare_rho_in_time(rho::Array, rho_ref::Array; smooth_const=1e-9) return rho_rel end -function compare_rho_in_time(rho::Array, rho_ref::OperatorVector; smooth_const=1e-9) +function compare_rho_in_time(rho::OperatorVectorArray, rho_ref::OperatorVector; smooth_const=1e-9) rho_ref_array = operator_recast(rho_ref) return compare_rho_in_time(rho, rho_ref_array; smooth_const=smooth_const) end -function compare_rho_in_time(rho::OperatorVector, rho_ref::Array; smooth_const=1e-9) +function compare_rho_in_time(rho::OperatorVector, rho_ref::OperatorVectorArray; smooth_const=1e-9) rho_array = operator_recast(rho) return compare_rho_in_time(rho_array, rho_ref; smooth_const=smooth_const) end diff --git a/test/local_test.ipynb b/test/local_test.ipynb new file mode 100644 index 0000000..d6528c6 --- /dev/null +++ b/test/local_test.ipynb @@ -0,0 +1,156 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, + "outputs": [], + "source": [ + "using Revise\n", + "using LinearAlgebra\n", + "using SparseArrays\n", + "using Pkg\n", + "using DelayDiffEq \n", + "using Plots\n", + "using BenchmarkTools\n", + "using LaTeXStrings" + ] + }, + { + "cell_type": "code", + "execution_count": 61, + "metadata": {}, + "outputs": [], + "source": [ + "using OpenQuantumSystems" + ] + }, + { + "cell_type": "code", + "execution_count": 94, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "postprocessing | \u001b[32m 8 \u001b[39m\u001b[36m 8\u001b[39m\n" + ] + }, + { + "data": { + "text/plain": [ + "Test.DefaultTestSet(\"postprocessing\", Any[], 8, false, false)" + ] + }, + "execution_count": 94, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "include(\"test_postprocessing.jl\")" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "operators dense | \u001b[32m 6 \u001b[39m\u001b[36m 6\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "superoperators | \u001b[32m 2 \u001b[39m\u001b[36m 2\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "metrics | \u001b[32m 19 \u001b[39m\u001b[36m 19\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "molecules | \u001b[32m 20 \u001b[39m\u001b[36m 20\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "aggregateCore | \u001b[32m 15 \u001b[39m\u001b[36m 15\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "aggregateTools | \u001b[32m 13 \u001b[39m\u001b[36m 13\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "aggregateOperators | \u001b[32m 14 \u001b[39m\u001b[36m 14\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "aggregate | \u001b[32m 7 \u001b[39m\u001b[36m 7\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "evolution | \u001b[32m 104 \u001b[39m\u001b[36m 104\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "schroedinger | \u001b[32m 25 \u001b[39m\u001b[36m 25\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "liouville | \u001b[32m 44 \u001b[39m\u001b[36m 44\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "interaction picture | \u001b[32m 12 \u001b[39m\u001b[36m 12\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: |\u001b[22m\u001b[39m\n", + "master | \u001b[36mNo tests\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "trace | \u001b[32m 41 \u001b[39m\u001b[36m 41\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "initial state | \u001b[32m 13 \u001b[39m\u001b[36m 13\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "memory kernel | \u001b[32m 51 \u001b[39m\u001b[36m 51\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: |\u001b[22m\u001b[39m\n", + "master ansatz | \u001b[36mNo tests\u001b[39m\n", + "\u001b[37m\u001b[1mTest Summary: | \u001b[22m\u001b[39m\u001b[32m\u001b[1mPass \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n", + "scoring | \u001b[32m 1 \u001b[39m\u001b[36m 1\u001b[39m\n" + ] + } + ], + "source": [ + "names = [\n", + " \"test_operators_dense.jl\",\n", + " \"test_superoperators.jl\",\n", + " \"test_metrics.jl\",\n", + " \"test_molecules.jl\",\n", + " \"test_aggregateCore.jl\",\n", + " \"test_aggregateTools.jl\",\n", + " \"test_aggregateOperators.jl\",\n", + " \"test_aggregate.jl\",\n", + " \"test_evolution.jl\",\n", + " \"test_schroedinger.jl\",\n", + " \"test_liouville.jl\",\n", + " \"test_interaction_picture.jl\",\n", + " \"test_master_exact.jl\",\n", + " \"test_trace.jl\",\n", + " \"test_initial_state.jl\",\n", + " \"test_memory_kernel.jl\",\n", + " \"test_master_ansatz.jl\",\n", + " \"test_postprocessing.jl\",\n", + " \"test_scoring.jl\"\n", + "]\n", + "\n", + "for name in names\n", + " include(name)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.3", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.3" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/test/runtests.jl b/test/runtests.jl index c6507db..138f554 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,6 +18,8 @@ names = [ "test_initial_state.jl", "test_memory_kernel.jl", "test_master_ansatz.jl", + "test_postprocessing.jl", + "test_scoring.jl", ] #= diff --git a/test/test_aggregate.jl b/test/test_aggregate.jl index f35d25a..cac6566 100644 --- a/test/test_aggregate.jl +++ b/test/test_aggregate.jl @@ -63,4 +63,22 @@ end @test agg.core == aggCore @test agg.tools == aggTools @test agg.operators == aggOperators + + agg = Aggregate(aggCore, nothing, nothing) + agg_ = Aggregate(aggCore) + @test agg_ == agg + + agg = setupAggregate(aggCore) + aggTools_ = AggregateTools(agg) + @test aggTools_ == aggTools + + aggOperators_ = AggregateOperators(agg) + @test aggOperators_ == aggOperators + + agg = setupAggregate(aggCore) + agg_ = Aggregate(aggCore, nothing, nothing) + setupAggregate!(agg_) + @test agg_ == agg + + end diff --git a/test/test_aggregateCore.jl b/test/test_aggregateCore.jl index c8667ac..c27e0e3 100644 --- a/test/test_aggregateCore.jl +++ b/test/test_aggregateCore.jl @@ -71,4 +71,7 @@ end @test OpenQuantumSystems.elIndOrder([1, 1, 1]) == 1 @test OpenQuantumSystems.elIndOrder([2, 1, 1]) == 2 @test OpenQuantumSystems.elIndOrder([1, 2, 1]) == 3 + + aggCore_ = AggregateCore([mol1, mol2], [0.0 0.0 0.0; 0.0 0.0 50.0; 0.0 50.0 0.0]) + @test aggCore_ == aggCore end diff --git a/test/test_aggregateOperators.jl b/test/test_aggregateOperators.jl index cae7cc3..4932838 100644 --- a/test/test_aggregateOperators.jl +++ b/test/test_aggregateOperators.jl @@ -78,4 +78,47 @@ end ref = aggOperators.Ham_0.data + aggOperators.Ham_I.data @test 1e-14 > D(aggOperators.Ham.data, ref) + # groundEnergy = false + aggOperators = AggregateOperators(aggCore, aggTools; groundEnergy=false) + @test 1e-14 > D(aggOperators.Ham_sys.data, [5.0 0.0 0.0; 0.0 203.5 50.0; 0.0 50.0 306.0]) + + ref = [0.25 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.5500000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.8499999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.4499999999999993 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.75 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0499999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.6500000000000004 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.9500000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.25] + @test 1e-14 > D(aggOperators.Ham_bath.data, ref) + + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -20.258948470231473 -14.325239843009513 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -14.325239843009516 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 7.162619921504748 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -7.1626199215047555 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -7.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 1.7906549803761975 2.5323685587789493 2.5323685587789493; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -1.790654980376199 2.1057475427305525e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -1.7906549803761995; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 14.325239843009498 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 20.25894847023148 -14.325239843009513 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 0.0 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 20.25894847023148 1.6845980341844322e-14 -14.325239843009516 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 0.0 301.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -10.12947423511573 14.325239843009507 -14.325239843009502 7.162619921504748 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 0.0 0.0 301.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -7.1626199215047555 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 0.0 0.0 0.0 0.0 301.0 0.0 0.0 0.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -7.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 0.0 0.0 0.0 0.0 301.0 0.0 0.0 0.0; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 1.7906549803761975 -2.532368558778949 2.532368558778948 0.0 0.0 0.0 0.0 0.0 0.0 301.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -1.790654980376199 -1.111293314979136e-15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -1.7906549803761995 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.0] + @test 1e-14 > D(aggOperators.Ham_S.data, ref) + + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0] + @test 1e-14 > D(aggOperators.Ham_B.data, ref) + + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -20.258948470231473 -14.325239843009513 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.1 0.0 0.0 0.0 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -14.325239843009516 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.7 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 7.162619921504748 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.0 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -7.1626199215047555 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.3 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -7.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.9 0.0 0.0 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 1.7906549803761975 2.5323685587789493 2.5323685587789493; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.2 0.0 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -1.790654980376199 2.1057475427305525e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.5 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -1.7906549803761995; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 -14.325239843009513 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 0.0 301.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 1.6845980341844322e-14 -14.325239843009516 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 0.0 301.6 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.12947423511573 14.325239843009507 -14.325239843009502 7.162619921504748 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 0.0 0.0 301.2 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -7.1626199215047555 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 0.0 0.0 0.0 0.0 301.5 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -7.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 0.0 0.0 0.0 0.0 301.8 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 1.7906549803761975 -2.532368558778949 2.532368558778948 0.0 0.0 0.0 0.0 0.0 0.0 301.4 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -1.790654980376199 -1.111293314979136e-15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.7 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -1.7906549803761995 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.0] + @test 1e-14 > D(aggOperators.Ham_0.data, ref) + + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0 -2.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0 -2.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 0.0 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 0.0] + @test 1e-9 > D(aggOperators.Ham_I.data, ref) + + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 0.0 14.325239843009498 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.8 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 0.0 -20.258948470231473 -14.325239843009513 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.1 0.0 0.0 -0.7071067811865475 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -14.325239843009516 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 198.7 0.0 0.0 -1.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 7.162619921504748 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 199.0 0.0 0.0 -1.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -7.1626199215047555 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.7071067811865475 0.0 0.0 199.3 0.0 0.0 -1.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -7.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 198.9 0.0 0.0 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 1.7906549803761975 2.5323685587789493 2.5323685587789493; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 199.2 0.0 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -1.790654980376199 2.1057475427305525e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 199.5 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -1.7906549803761995; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 301.0 -1.414213562373095 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 -14.325239843009513 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 -1.414213562373095 301.3 -2.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 1.6845980341844322e-14 -14.325239843009516 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 -2.0 301.6 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.12947423511573 14.325239843009507 -14.325239843009502 7.162619921504748 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 0.0 0.0 301.2 -1.414213562373095 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -7.1626199215047555 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 0.0 0.0 0.0 -1.414213562373095 301.5 -2.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -7.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 0.0 0.0 0.0 -2.0 301.8 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 1.7906549803761975 -2.532368558778949 2.532368558778948 0.0 0.0 0.0 0.0 0.0 0.0 301.4 -1.414213562373095 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -1.790654980376199 -1.111293314979136e-15 0.0 0.0 0.0 0.0 0.0 0.0 -1.414213562373095 301.7 -2.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -1.7906549803761995 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -2.0 302.0] + @test 1e-14 > D(aggOperators.Ham.data, ref) + + # groundEnergy = false + Ham_sys = getAggHamSystemSmall(aggCore, aggTools; groundEnergy = false) + ref = [0.0 0.0 0.0; 0.0 198.5 50.0; 0.0 50.0 301.0] + @test 1e-14 > D(Ham_sys.data, ref) + + Ham_S = getAggHamSystemBig(aggCore, aggTools; groundEnergy = false) + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -20.258948470231473 -14.325239843009513 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -14.325239843009516 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 7.162619921504748 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -7.1626199215047555 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -7.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 1.7906549803761975 2.5323685587789493 2.5323685587789493; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -1.790654980376199 2.1057475427305525e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -1.7906549803761995; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 14.325239843009498 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 20.25894847023148 -14.325239843009513 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 0.0 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 20.25894847023148 1.6845980341844322e-14 -14.325239843009516 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 0.0 301.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 -10.12947423511573 14.325239843009507 -14.325239843009502 7.162619921504748 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 0.0 0.0 301.0 0.0 0.0 0.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -7.1626199215047555 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 0.0 0.0 0.0 0.0 301.0 0.0 0.0 0.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -7.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 0.0 0.0 0.0 0.0 301.0 0.0 0.0 0.0; 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 1.7906549803761975 -2.532368558778949 2.532368558778948 0.0 0.0 0.0 0.0 0.0 0.0 301.0 0.0 0.0; 0.0 -0.0 -0.0 0.0 -0.0 -0.0 0.0 -0.0 -0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -1.790654980376199 -1.111293314979136e-15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.0 0.0; 0.0 0.0 -0.0 0.0 0.0 -0.0 0.0 0.0 -0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -1.7906549803761995 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.0] + @test 1e-14 > D(Ham_S.data, ref) + + Ham_bath = getAggHamBathSmall(aggCore, aggTools; groundEnergy = false) + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0] + @test 1e-14 > D(Ham_bath.data, ref) + + Ham_B = getAggHamBathBig(aggCore, aggTools; groundEnergy = false) + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0] + @test 1e-14 > D(Ham_B.data, ref) + + Ham_0 = getAggHamSystemBath(aggCore, aggTools; groundEnergy = false) + ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.3000000000000007 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.5999999999999996 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.1999999999999993 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.7999999999999998 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.40000000000000036 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7000000000000002 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 20.25894847023148 20.25894847023148 -10.12947423511573 -14.325239843009511 -14.325239843009511 5.064737117557862 7.162619921504751 7.162619921504751; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -20.258948470231473 -14.325239843009513 1.6845980341844322e-14 14.325239843009507 10.129474235115739 -1.1911906935453394e-14 -7.162619921504749 -5.064737117557866 5.955953467726693e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.1 0.0 0.0 0.0 0.0 0.0 0.0 20.258948470231466 -8.890346519833037e-15 -14.325239843009516 -14.325239843009502 6.286424311272163e-15 10.12947423511574 7.162619921504747 -3.1432121556360795e-15 -5.064737117557867; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.7 0.0 0.0 0.0 0.0 0.0 10.12947423511573 14.325239843009513 14.325239843009513 7.162619921504748 10.129474235115739 10.129474235115739 -10.743929882257122 -15.194211352673607 -15.194211352673607; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.0 0.0 0.0 0.0 0.0 -14.325239843009511 -10.12947423511574 1.1911906935453396e-14 -10.129474235115737 -7.1626199215047555 8.422990170922161e-15 15.194211352673603 10.743929882257133 -1.2634485256383239e-14; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.3 0.0 0.0 0.0 14.325239843009504 -6.286424311272164e-15 -10.129474235115742 10.129474235115731 -4.445173259916518e-15 -7.162619921504756 -15.194211352673598 6.667759889874776e-15 10.743929882257133; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 198.9 0.0 0.0 5.064737117557865 7.1626199215047555 7.1626199215047555 10.743929882257127 15.194211352673614 15.194211352673614 1.7906549803761975 2.5323685587789493 2.5323685587789493; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.2 0.0 -7.162619921504754 -5.064737117557869 5.955953467726697e-15 -15.194211352673612 -10.743929882257138 1.2634485256383247e-14 -2.532368558778949 -1.790654980376199 2.1057475427305525e-15; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 199.5 7.162619921504751 -3.1432121556360815e-15 -5.06473711755787 15.194211352673607 -6.66775988987478e-15 -10.743929882257138 2.532368558778948 -1.111293314979136e-15 -1.7906549803761995; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 14.325239843009498 -20.258948470231473 20.258948470231466 10.12947423511573 -14.325239843009511 14.325239843009504 5.064737117557865 -7.162619921504754 7.162619921504751 301.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 -14.325239843009513 -8.890346519833037e-15 14.325239843009513 -10.12947423511574 -6.286424311272164e-15 7.1626199215047555 -5.064737117557869 -3.1432121556360815e-15 0.0 301.3 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20.25894847023148 1.6845980341844322e-14 -14.325239843009516 14.325239843009513 1.1911906935453396e-14 -10.129474235115742 7.1626199215047555 5.955953467726697e-15 -5.06473711755787 0.0 0.0 301.6 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -10.12947423511573 14.325239843009507 -14.325239843009502 7.162619921504748 -10.129474235115737 10.129474235115731 10.743929882257127 -15.194211352673612 15.194211352673607 0.0 0.0 0.0 301.2 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 10.129474235115739 6.286424311272163e-15 10.129474235115739 -7.1626199215047555 -4.445173259916518e-15 15.194211352673614 -10.743929882257138 -6.66775988987478e-15 0.0 0.0 0.0 0.0 301.5 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -14.325239843009511 -1.1911906935453394e-14 10.12947423511574 10.129474235115739 8.422990170922161e-15 -7.162619921504756 15.194211352673614 1.2634485256383247e-14 -10.743929882257138 0.0 0.0 0.0 0.0 0.0 301.8 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.064737117557862 -7.162619921504749 7.162619921504747 -10.743929882257122 15.194211352673603 -15.194211352673598 1.7906549803761975 -2.532368558778949 2.532368558778948 0.0 0.0 0.0 0.0 0.0 0.0 301.4 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 -5.064737117557866 -3.1432121556360795e-15 -15.194211352673607 10.743929882257133 6.667759889874776e-15 2.5323685587789493 -1.790654980376199 -1.111293314979136e-15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 301.7 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.162619921504751 5.955953467726693e-15 -5.064737117557867 -15.194211352673607 -1.2634485256383239e-14 10.743929882257133 2.5323685587789493 2.1057475427305525e-15 -1.7906549803761995 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 302.0] + @test 1e-14 > D(Ham_0.data, ref) + end diff --git a/test/test_aggregateTools.jl b/test/test_aggregateTools.jl index 85c777c..603b663 100644 --- a/test/test_aggregateTools.jl +++ b/test/test_aggregateTools.jl @@ -84,4 +84,7 @@ end @test aggTools.basisBath == GenericBasis([4]) @test aggTools.basis == GenericBasis([12]) + FCfactors = getFranckCondonFactors(aggCore) + @test 1e-14 > D(aggTools.FCfactors, FCfactors) + end \ No newline at end of file diff --git a/test/test_evolution.jl b/test/test_evolution.jl index 3d496f0..a7f22da 100644 --- a/test/test_evolution.jl +++ b/test/test_evolution.jl @@ -13,6 +13,10 @@ using Random, SparseArrays, LinearAlgebra, StableRNGs D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) = abs(tracedistance_nh(dense(op1), dense(op2))) + tspan = get_tspan(0., 1., 10) + ref = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + @test tspan == ref + mode1 = Mode(0.2, 1.0) mode2 = Mode(0.3, 2.0) Energy = [0.0, 200.0] @@ -68,7 +72,8 @@ using Random, SparseArrays, LinearAlgebra, StableRNGs @test 1e-10 > D(U_sop_array[3], U_sop3) t = 0.0 - foreach(evolutionOperatorIterator(Ham, tspan)) do U_op + foreach(evolutionOperatorIterator(Ham, tspan; diagonalize = true, approximate = false) + ) do U_op U_op_ref = evolutionOperator(Ham, t) @test 1e-11 > D(U_op, U_op_ref) t += 0.5 @@ -93,7 +98,9 @@ using Random, SparseArrays, LinearAlgebra, StableRNGs end t = 0.0 - foreach(evolutionSuperOperatorIterator(Ham, tspan)) do U_sop + foreach( + evolutionSuperOperatorIterator(Ham, tspan; diagonalize = false, approximate = true) + ) do U_sop U_sop_ref = evolutionSuperOperator(Ham, t) @test 1e-10 > D(U_sop, U_sop_ref) t += 0.5 @@ -101,12 +108,7 @@ using Random, SparseArrays, LinearAlgebra, StableRNGs t = 0.0 foreach( - evolutionSuperOperatorIterator( - Ham, - tspan; - diagonalize = false, - approximate = false, - ), + evolutionSuperOperatorIterator(Ham, tspan; diagonalize = false, approximate = false), ) do U_op U_op_ref = evolutionSuperOperator(Ham, t) @test 1e-12 > D(U_op, U_op_ref) @@ -186,28 +188,87 @@ using Random, SparseArrays, LinearAlgebra, StableRNGs for t_i = 1:length(tspan) U_op = U_op_array[t_i] rho_ref = U_op * rho0 * U_op' - 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) + @test 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) end T, rho_t = evolution_exact(rho0, tspan, Ham; diagonalize = true) for t_i = 1:length(tspan) U_op = U_op_array[t_i] rho_ref = U_op * rho0 * U_op' - 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) + @test 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) end T, rho_t = evolution_approximate(rho0, tspan, Ham; diagonalize = false) for t_i = 1:length(tspan) U_op = U_op_array[t_i] rho_ref = U_op * rho0 * U_op' - 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) + @test 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) end T, rho_t = evolution_approximate(rho0, tspan, Ham; diagonalize = true) for t_i = 1:length(tspan) U_op = U_op_array[t_i] rho_ref = U_op * rho0 * U_op' - 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) + @test 1e-12 > D(rho_t_ref[t_i], rho_t[t_i]) end + t = 1. + indicesMap = agg.tools.indicesMap + U_part = evolution_el_part(Ham, t, 2, 2, indicesMap) + data = take_el_part(Ham.data, 2, 2, indicesMap) + b = GenericBasis([size(data, 1)]) + data = evolutionOperator(data, t) + ref = DenseOperator(b, b, data) + @test 1e-12 > D(ref, U_part) + + + mode1 = Mode(0.2, 1.0) + mode2 = Mode(0.3, 2.0) + Energy = [0.0, 200.0] + mol1 = Molecule([mode1], 2, [2.0, 200.0]) + mol2 = Molecule([mode2], 2, [3.0, 300.0]) + aggCore = AggregateCore([mol1, mol2]) + aggCore.coupling[2, 3] = 50 + aggCore.coupling[3, 2] = 50 + agg = setupAggregate(aggCore) + aggTools = agg.tools + aggOperators = agg.operators + + Ham_I = agg.operators.Ham_I + Ham_0 = agg.operators.Ham_0 + Ham = agg.operators.Ham + + basis = agg.tools.basis + indicesLen = agg.tools.bSize + indices = agg.tools.indices + indicesMap = agg.tools.indicesMap + FCFact = agg.tools.FCfactors + FCProd = agg.tools.FCproduct + + W0 = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25000143755694526 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2499997124870238 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2500002875090083 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24999856244702248] + W0 = DenseOperator(basis, basis, complex(W0)) + + t = 1. + tspan = get_tspan(0., t, 1) + _, W_int_t_ = Evolution_SI_exact(W0, tspan, agg) + U_op = evolutionOperator(Ham, t) + W = U_op * W0 * U_op' + U_0_op = evolutionOperator(Ham_0, t) + W = U_0_op' * W * U_0_op + W_int_t_ref = W.data + @test 1e-12 > D(W_int_t_ref, W_int_t_[2, :, :]) + + _, rho_int_t_ = Evolution_sI_exact(W0, tspan, agg) + rho_int_t_ref = trace_bath(W_int_t_ref, agg.core, agg.tools) + @test 1e-12 > D(rho_int_t_ref, rho_int_t_[2, :, :]) + + _, W_t_ = Evolution_SS_exact(W0, tspan, agg) + U_op = evolutionOperator(Ham, t) + W = U_op * W0 * U_op' + W_t_ref = W.data + @test 1e-12 > D(W_t_ref, W_t_[2, :, :]) + + _, rho_t_ = Evolution_sS_exact(W0, tspan, agg) + rho_t_ref = trace_bath(W_t_ref, agg.core, agg.tools) + @test 1e-12 > D(rho_t_ref, rho_t_[2, :, :]) end diff --git a/test/test_initial_state.jl b/test/test_initial_state.jl index 019d2a1..93b6a3a 100644 --- a/test/test_initial_state.jl +++ b/test/test_initial_state.jl @@ -43,6 +43,9 @@ import QuantumOpticsBase W0 = thermal_state(T, mu_array, aggCore, aggTools, aggOperators; diagonalize = true) @test 1e-15 > D(W0_ref, W0.data) + W0 = thermal_state(T, mu_array, agg; diagonalize = true) + @test 1e-15 > D(W0_ref, W0.data) + W01 = thermal_state(T, [[1, 2, 1]], aggCore, aggTools, aggOperators; diagonalize = true) W02 = thermal_state(T, [[1, 1, 2]], aggCore, aggTools, aggOperators; diagonalize = true) W0_composite_ref = 0.8 * W01 + 0.2 * W02 @@ -50,6 +53,9 @@ import QuantumOpticsBase W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], aggCore, aggTools, aggOperators; diagonalize = true) @test 1e-15 > D(W0_composite_ref, W0) + W0 = thermal_state_composite(T, [0.0, 0.8, 0.2], agg; diagonalize = true) + @test 1e-15 > D(W0_composite_ref, W0) + W0_ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25029983141333245 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.24993996472850133 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25005986276494624 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.2497003410932199] W0 = thermal_state(T, mu_array, aggCore, aggTools, aggOperators; diagonalize = false, diagonal = true) @test 1e-15 > D(W0_ref, W0.data) @@ -95,4 +101,31 @@ import QuantumOpticsBase ) @test 1e-15 > D(W0_composite_ref, W0) + T = 0.5 + mu_array = [[1, 1, 2]] + W0_ref = [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.7750674418239181 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.06464140010302448 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1479518276234691 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.012339330449588293] + W0 = thermal_state(T, mu_array, aggCore, aggTools, aggOperators; diagonalize = true) + @test 1e-15 > D(W0_ref, W0.data) + + W0, rho0, W0_bath = ultrafast_laser_excitation(T, [0.0, 0.8, 0.2], agg; diagonalize = true) + W0_ref = thermal_state_composite( + T, + [0.0, 0.8, 0.2], + aggCore, + aggTools, + aggOperators; + diagonalize = true + ) + normalize!(W0_ref) + W0_ref = DenseOperator(W0_ref.basis_l, W0_ref.basis_r, complex(W0_ref.data)) + + rho0_ref = trace_bath(W0_ref, agg.core, agg.tools) + rho0_ref = DenseOperator(rho0_ref.basis_l, rho0_ref.basis_r, complex(rho0_ref.data)) + + W0_bath_ref = get_rho_bath(W0_ref, agg.core, agg.tools) + W0_bath_ref = DenseOperator(W0_bath_ref.basis_l, W0_bath_ref.basis_r, complex(W0_bath_ref.data)) + @test 1e-15 > D(W0_ref, W0) + @test 1e-15 > D(rho0_ref, rho0) + @test 1e-15 > D(W0_bath_ref, W0_bath) + end # testset diff --git a/test/test_postprocessing.jl b/test/test_postprocessing.jl new file mode 100644 index 0000000..8b9e0f1 --- /dev/null +++ b/test/test_postprocessing.jl @@ -0,0 +1,103 @@ +using Test +using OpenQuantumSystems +using Random, SparseArrays, LinearAlgebra, StableRNGs + +@testset "postprocessing" begin + + Random.seed!(StableRNG(0), 1) + + D(op1::Array, op2::Array) = abs(norm(op1 - op2)) + D(x1::StateVector, x2::StateVector) = norm(x2 - x1) + D(op1::AbstractOperator, op2::AbstractOperator) = + abs(tracedistance_nh(dense(op1), dense(op2))) + D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) = + abs(tracedistance_nh(dense(op1), dense(op2))) + + + mode1 = Mode(0.2, 1.0) + mode2 = Mode(0.3, 2.0) + Energy = [0.0, 200.0] + mol1 = Molecule([mode1], 3, [2.0, 200.0]) + mol2 = Molecule([mode2], 3, [3.0, 300.0]) + aggCore = AggregateCore([mol1, mol2]) + aggCore.coupling[2, 3] = 50 + aggCore.coupling[3, 2] = 50 + agg = setupAggregate(aggCore) + aggTools = agg.tools + aggOperators = agg.operators + + Ham_B = agg.operators.Ham_B + Ham_I = agg.operators.Ham_I + Ham_0 = agg.operators.Ham_0 + Ham = agg.operators.Ham + + Ham_0_lambda, Ham_0_S = eigen(Ham_0.data) + Ham_0_Sinv = inv(Ham_0_S) + Ham_0_lambda = diagm(Ham_0_lambda) + + basis = agg.tools.basis + indicesLen = agg.tools.bSize + indices = agg.tools.indices + indicesMap = agg.tools.indicesMap + FCFact = agg.tools.FCfactors + FCProd = agg.tools.FCproduct + + # recast of OperatorVectorArray + tspan = get_tspan(0., 0.02, 100) + W0, rho0, W0_bath = ultrafast_laser_excitation(10., [0.0, 0.3, 0.7], agg) + _, rho_t_1 = Evolution_sS_exact(W0, tspan, agg) + _, rho_t_2 = LvN_sS( + W0, + tspan, + agg; + reltol = 1e-12, + abstol = 1e-12, + alg = OrdinaryDiffEq.Tsit5(), + ) + + rho_t_1_recasted = operator_recast(rho_t_1) + @test typeof(rho_t_1) == OperatorVectorArray + @test typeof(rho_t_1_recasted) == OperatorVectorArray + + rho_t_2_recasted = operator_recast(rho_t_2) + @test typeof(rho_t_2) == OperatorVector + @test typeof(rho_t_2_recasted) == OperatorVectorArray + + # interaction and schrodinger picture + _, rho_int_t_ref = Evolution_sI_exact(W0, tspan, agg) + rho_t = interaction_pic_to_schroedinger_pic(rho_int_t_ref, tspan, agg) + rho_int_t = schroedinger_pic_to_interaction_pic(rho_t, tspan, agg) + @test D(rho_int_t_ref, rho_int_t) < 1e-14 + + _, rho_int_t_ref = LvN_sI( + W0, + tspan, + agg; + reltol = 1e-12, + abstol = 1e-12, + alg = OrdinaryDiffEq.Tsit5(), + ) + rho_t = interaction_pic_to_schroedinger_pic(rho_int_t_ref, tspan, agg) + rho_int_t = schroedinger_pic_to_interaction_pic(rho_t, tspan, agg) + rho_int_t_ref = operator_recast(rho_int_t_ref) + @test D(rho_int_t_ref, rho_int_t) < 1e-14 + + # local and exciton basis + _, rho_t_ref = Evolution_sS_exact(W0, tspan, agg) + rho_t_exc = local_st_to_exciton_st(rho_t_ref, agg) + rho_t = exciton_st_to_local_st(rho_t_exc, agg) + @test D(rho_t_ref, rho_t) < 1e-14 + + _, rho_t_ref = LvN_sS( + W0, + tspan, + agg; + reltol = 1e-12, + abstol = 1e-12, + alg = OrdinaryDiffEq.Tsit5(), + ) + rho_t_exc = local_st_to_exciton_st(rho_t_ref, agg) + rho_t = exciton_st_to_local_st(rho_t_exc, agg) + rho_t_ref = operator_recast(rho_t_ref) + @test D(rho_t_ref, rho_t) < 1e-14 +end \ No newline at end of file diff --git a/test/test_scoring.jl b/test/test_scoring.jl new file mode 100644 index 0000000..916720e --- /dev/null +++ b/test/test_scoring.jl @@ -0,0 +1,138 @@ +using Test +using OpenQuantumSystems +using Random, SparseArrays, LinearAlgebra, StableRNGs + +@testset "scoring" begin + + Random.seed!(StableRNG(0), 1) + + D(op1::Array, op2::Array) = abs(norm(op1 - op2)) + D(x1::StateVector, x2::StateVector) = norm(x2 - x1) + D(op1::AbstractOperator, op2::AbstractOperator) = + abs(tracedistance_nh(dense(op1), dense(op2))) + D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) = + abs(tracedistance_nh(dense(op1), dense(op2))) + + + mode1 = Mode(0.2, 1.0) + mode2 = Mode(0.3, 2.0) + Energy = [0.0, 200.0] + mol1 = Molecule([mode1], 3, [2.0, 200.0]) + mol2 = Molecule([mode2], 3, [3.0, 300.0]) + aggCore = AggregateCore([mol1, mol2]) + aggCore.coupling[2, 3] = 50 + aggCore.coupling[3, 2] = 50 + agg = setupAggregate(aggCore) + aggTools = agg.tools + aggOperators = agg.operators + + Ham_B = agg.operators.Ham_B + Ham_I = agg.operators.Ham_I + Ham_0 = agg.operators.Ham_0 + Ham = agg.operators.Ham + + Ham_0_lambda, Ham_0_S = eigen(Ham_0.data) + Ham_0_Sinv = inv(Ham_0_S) + Ham_0_lambda = diagm(Ham_0_lambda) + + basis = agg.tools.basis + indicesLen = agg.tools.bSize + indices = agg.tools.indices + indicesMap = agg.tools.indicesMap + FCFact = agg.tools.FCfactors + FCProd = agg.tools.FCproduct + + tspan = get_tspan(0., 0.02, 100) + W0, rho0, W0_bath = ultrafast_laser_excitation(10., [0.0, 0.3, 0.7], agg) + _, W_t_1 = Evolution_SS_exact(W0, tspan, agg) + _, W_t_2 = Evolution_SS_exact(W0, tspan, agg) + rmse = get_rmse_in_time(W_t_1, W_t_2) + @test rmse == 0. + + _, W_t_3 = LvN_SS( + W0, + tspan, + agg; + reltol = 1e-10, + abstol = 1e-10, + alg = OrdinaryDiffEq.Tsit5(), + ) + rmse = get_rmse_in_time(W_t_3, W_t_3) + @test rmse == 0. + rmse = get_rmse_in_time(W_t_1, W_t_3) + @test rmse < 1e-9 + rmse = get_rmse_in_time(W_t_3, W_t_1) + @test rmse < 1e-9 + + + tspan = get_tspan(0., 0.02, 100) + W0, rho0, W0_bath = ultrafast_laser_excitation(10., [0.0, 0.3, 0.7], agg) + _, rho_t_1 = Evolution_sS_exact(W0, tspan, agg) + _, rho_t_2 = Evolution_sS_exact(W0, tspan, agg) + score = compare_rho(rho_t_1, rho_t_2) + @test sum(score) == 0. + + _, rho_t_3 = LvN_sS( + W0, + tspan, + agg; + reltol = 1e-12, + abstol = 1e-12, + alg = OrdinaryDiffEq.Tsit5(), + ) + score = compare_rho(rho_t_3, rho_t_3) + @test sum(score) == 0. + + score = compare_rho(rho_t_1, rho_t_3) + @test sum(score) < 1e-9 + + score = compare_rho(rho_t_3, rho_t_1) + @test sum(score) < 1e-9 + + tspan = get_tspan(0., 0.02, 100) + _, rho_t_1 = Evolution_sS_exact(W0, tspan, agg) + _, rho_t_2 = LvN_sS( + W0, + tspan, + agg; + reltol = 1e-12, + abstol = 1e-12, + alg = OrdinaryDiffEq.Tsit5(), + ) + + score_ref = compare_rho(rho_t_2, rho_t_2) + score_t = compare_rho_in_time(rho_t_2, rho_t_2; smooth_const=1e-9) + N, M, K = size(rho_t_1) + score = zeros(Float64, M, K) + for t_i in 1:N + score[:, :] += score_t[t_i, :, :] + end + @test D(score_ref, score) == 0. + + score_ref = compare_rho(rho_t_1, rho_t_2) + score_t = compare_rho_in_time(rho_t_1, rho_t_2; smooth_const=1e-9) + N, M, K = size(rho_t_1) + score = zeros(Float64, M, K) + for t_i in 1:N + score[:, :] += score_t[t_i, :, :] + end + @test D(score_ref, score) < 1e-8 + + score_ref = compare_rho(rho_t_2, rho_t_1) + score_t = compare_rho_in_time(rho_t_2, rho_t_1; smooth_const=1e-9) + N, M, K = size(rho_t_1) + score = zeros(Float64, M, K) + for t_i in 1:N + score[:, :] += score_t[t_i, :, :] + end + @test D(score_ref, score) < 1e-8 + + score_ref = compare_rho(rho_t_1, rho_t_1) + score_t = compare_rho_in_time(rho_t_1, rho_t_1; smooth_const=1e-9) + N, M, K = size(rho_t_1) + score = zeros(Float64, M, K) + for t_i in 1:N + score[:, :] += score_t[t_i, :, :] + end + @test D(score_ref, score) == 0. +end \ No newline at end of file