Skip to content

Commit

Permalink
On-the-fly transfer matrix generation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@amartyabose
amartyabose committed Sep 12, 2024
1 parent 26de7de commit 574fa34
Show file tree
Hide file tree
Showing 7 changed files with 78 additions and 46 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "QuantumDynamics"
uuid = "5a40c529-53c2-4483-a223-e00c1cee8134"
authors = ["Amartya Bose <[email protected]> and contributors"]
version = "0.2.19"
version = "0.2.20"

[deps]
AtomsIO = "1692102d-eeb4-4df9-807b-c9517f998d44"
Expand Down
65 changes: 36 additions & 29 deletions src/DynamicMap_MasterEquation/TTM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ module TTM

using HDF5
using FLoops
using ..EtaCoefficients, ..SpectralDensities, ..Blip, ..Utilities
using ..SpectralDensities, ..Utilities

const references = """
- Cerrillo, J.; Cao, J. Non-Markovian Dynamical Maps: Numerical Processing of Open Quantum Trajectories. Phys. Rev. Lett. 2014, 112 (11), 110401. https://doi.org/10.1103/PhysRevLett.112.110401."""
Expand Down Expand Up @@ -55,35 +55,42 @@ function get_propagators_QuAPI(; fbU::Array{ComplexF64,3}, Jw::Vector{T}, β, dt
U0e, T0e
end

function get_improved_Ts(; fbU::Array{ComplexF64,3}, Jw::Vector{T}, β, dt, ntimes, rmax, kmax::Union{Int,Nothing}=nothing, path_integral_routine, extraargs::Utilities.ExtraArgs, svec=[1.0 -1.0], verbose::Bool=false, reference_prop=false) where {T<:SpectralDensities.SpectralDensity}
ηs = [EtaCoefficients.calculate_η(jw; β, dt, kmax=ntimes, imaginary_only=reference_prop) for jw in Jw]
_, _, _, sbar, Δs = Blip.setup_simulation(svec)
U0e_within_r = path_integral_routine(; fbU, Jw, β, dt, ntimes=rmax, kmax, extraargs, svec, verbose, reference_prop)
sdim2 = size(fbU, 2)
T0e = zeros(ComplexF64, ntimes, sdim2, sdim2)
U0e = zero(T0e)
for n = 1:rmax
U0e[n, :, :] .= U0e_within_r[n, :, :]
T0e[n, :, :] .= U0e[n, :, :]
for j = 1:n-1
T0e[n, :, :] .-= T0e[j, :, :] * U0e[n-j, :, :]
end
end
for n = rmax+1:ntimes
U0e[n, :, :] .= sum([T0e[j, :, :] * U0e[n-j, :, :] for j = 1:n-1])
for s1 = 1:sdim2, s2 = 1:sdim2
val = 0
for j = 1:length(Jw)
val -= Δs[j, s1] * (real(ηs[j].η0e[n]) * Δs[j, s2] + 2im * imag(ηs[j].η0e[n]) * sbar[j, s2])
end
U0e[n, s1, s2] *= exp(val)
end
T0e[n, :, :] .= U0e[n, :, :]
for j = 1:n-1
T0e[n, :, :] .-= T0e[j, :, :] * U0e[n-j, :, :]
end
# function get_improved_Ts(; fbU::Array{ComplexF64,3}, Jw::Vector{T}, β, dt, ntimes, rmax, kmax::Union{Int,Nothing}=nothing, path_integral_routine, extraargs::Utilities.ExtraArgs, svec=[1.0 -1.0], verbose::Bool=false, reference_prop=false) where {T<:SpectralDensities.SpectralDensity}
# ηs = [EtaCoefficients.calculate_η(jw; β, dt, kmax=ntimes, imaginary_only=reference_prop) for jw in Jw]
# _, _, _, sbar, Δs = Blip.setup_simulation(svec)
# U0e_within_r = path_integral_routine(; fbU, Jw, β, dt, ntimes=rmax, kmax, extraargs, svec, verbose, reference_prop)
# sdim2 = size(fbU, 2)
# T0e = zeros(ComplexF64, ntimes, sdim2, sdim2)
# U0e = zero(T0e)
# for n = 1:rmax
# U0e[n, :, :] .= U0e_within_r[n, :, :]
# T0e[n, :, :] .= U0e[n, :, :]
# for j = 1:n-1
# T0e[n, :, :] .-= T0e[j, :, :] * U0e[n-j, :, :]
# end
# end
# for n = rmax+1:ntimes
# U0e[n, :, :] .= sum([T0e[j, :, :] * U0e[n-j, :, :] for j = 1:n-1])
# for s1 = 1:sdim2, s2 = 1:sdim2
# val = 0
# for j = 1:length(Jw)
# val -= Δs[j, s1] * (real(ηs[j].η0e[n]) * Δs[j, s2] + 2im * imag(ηs[j].η0e[n]) * sbar[j, s2])
# end
# U0e[n, s1, s2] *= exp(val)
# end
# T0e[n, :, :] .= U0e[n, :, :]
# for j = 1:n-1
# T0e[n, :, :] .-= T0e[j, :, :] * U0e[n-j, :, :]
# end
# end
# U0e, T0e
# end

function update_Ts!(T0e::Array{<:Complex, 3}, U0e::Array{<:Complex, 3}, n::Int64)
T0e[n, :, :] .= U0e[n, :, :]
for j = 1:n-1
T0e[n, :, :] .-= T0e[j, :, :] * U0e[n-j, :, :]
end
U0e, T0e
end

"""
Expand Down
11 changes: 6 additions & 5 deletions src/PathIntegral/Blip.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ using FLoops

using LinearAlgebra

using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Utilities
using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Utilities, ..TTM

const references = """
- Makri, N. Blip decomposition of the path integral: Exponential acceleration of real-time calculations on quantum dissipative systems. The Journal of Chemical Physics 2014, 141 (13), 134117. https://doi.org/10.1063/1.4896736."""
Expand Down Expand Up @@ -145,15 +145,14 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Vect

ndim = length(group_states)
U0e = zeros(ComplexF64, ntimes, sdim2, sdim2)
T0e = zero(U0e)
if !isnothing(output)
if !from_TTM
Utilities.check_or_insert_value(output, "U0e", U0e)
end
Utilities.check_or_insert_value(output, "T0e", T0e)
Utilities.check_or_insert_value(output, "time_taken", zeros(Float64, ntimes))
Utilities.check_or_insert_value(output, "num_paths", zeros(Int64, ntimes))
for (j, ηj) in enumerate(η)
Utilities.check_or_insert_value(output, "eta_$(j)", ηj.ηmn)
end
end
@inbounds begin
for i = 1:ntimes
Expand Down Expand Up @@ -181,15 +180,17 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Vect
@reduce tmpU0e = zeros(ComplexF64, sdim2, sdim2) .+ get_total_amplitude(; tmpprops=propagators, path, group_Δs, sbar, η, propagator_type="0e", nsteps=i, sdim2, val1, valend, valjkp)
end
@inbounds U0e[i, :, :] .+= tmpU0e
TTM.update_Ts!(T0e, U0e, i)
end
if !isnothing(output)
output["U0e"][i, :, :] = U0e[i, :, :]
output["T0e"][i, :, :] = T0e[i, :, :]
output["time_taken"][i] = time_taken
output["num_paths"][i] = num_paths
flush(output)
end
if verbose
@info "Done time step $(i); # paths = $(sum(num_paths)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
@info "Done time step $(i); # paths = $(sum(num_paths)); trace(Tmat) = $(round(tr(T0e[i, :, :]); digits=3)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
end
end
end
Expand Down
8 changes: 7 additions & 1 deletion src/PathIntegral/PCTNPI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ using HDF5
using LinearAlgebra
using ITensors
using FLoops
using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Blip, ..Utilities
using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Blip, ..Utilities, ..TTM

const references = """
- Bose, A. Pairwise Connected Tensor Network Representation of Path Integrals. Physical Review B 2022, 105 (2), 024309. https://doi.org/10.1103/PhysRevB.105.024309."""
Expand Down Expand Up @@ -149,6 +149,7 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Abst
_, state_to_blip, group_Δs, sbar, Δs = Blip.setup_simulation(svec)

U0e = zeros(ComplexF64, ntimes, sdim2, sdim2)
T0e = zero(U0e)

# calculate for the first time step separately, s1->s2
U0e[1, :, :] = fbU[1, :, :]
Expand All @@ -161,11 +162,14 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Abst
U0e[1, t, :] .*= infl1
end
end
TTM.update_Ts!(T0e, U0e, 1)
if !isnothing(output)
if !from_TTM
Utilities.check_or_insert_value(output, "U0e", U0e)
end
Utilities.check_or_insert_value(output, "T0e", T0e)
output["U0e"][1, :, :] = U0e[1, :, :]
output["T0e"][1, :, :] = T0e[1, :, :]
flush(output)
end

Expand Down Expand Up @@ -199,8 +203,10 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Abst
prop *= tensor
end
U0e[j, :, :] = Utilities.convert_ITensor_to_matrix(prop, sites[1], sites[j+1])
TTM.update_Ts!(T0e, U0e, j)
if !isnothing(output)
output["U0e"][j, :, :] = U0e[j, :, :]
output["T0e"][j, :, :] = T0e[j, :, :]
flush(output)
end
end
Expand Down
16 changes: 12 additions & 4 deletions src/PathIntegral/QuAPI.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ module QuAPI

using HDF5
using FLoops
using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Utilities
using ..EtaCoefficients, ..Propagators, ..SpectralDensities, ..Utilities, ..TTM

const references = """
- Makri, N.; Makarov, D. E. Tensor Propagator for Iterative Quantum Time Evolution of Reduced Density Matrices. I. Theory. The Journal of Chemical Physics 1995, 102 (11), 4600–4610. https://doi.org/10.1063/1.469509.
Expand Down Expand Up @@ -450,10 +450,12 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Vect
@info "Starting propagation within memory"
end
U0e = Array{ComplexF64}(undef, ntimes, sdim2, sdim2)
T0e = zero(U0e)
if !isnothing(output)
if !from_TTM
Utilities.check_or_insert_value(output, "U0e", U0e)
end
Utilities.check_or_insert_value(output, "T0e", T0e)
Utilities.check_or_insert_value(output, "time_taken", zeros(Float64, ntimes))
Utilities.check_or_insert_value(output, "num_paths", zeros(Int64, ntimes))
end
Expand Down Expand Up @@ -482,15 +484,17 @@ function build_augmented_propagator(; fbU::AbstractArray{ComplexF64,3}, Jw::Vect
@reduce tmpU0e = zeros(ComplexF64, sdim2, sdim2) .+ tmpval
end
@inbounds U0e[i, :, :] .= tmpU0e
TTM.update_Ts!(T0e, U0e, i)
end
if !isnothing(output)
output["U0e"][i, :, :] = U0e[i, :, :]
output["T0e"][i, :, :] = T0e[i, :, :]
output["time_taken"][i] = time_taken
output["num_paths"][i] = num_paths
flush(output)
end
if verbose
@info "Done time step $(i); # paths = $(sum(num_paths)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
@info "Done time step $(i); # paths = $(sum(num_paths)); trace(Tmat) = $(round(tr(T0e[i, :, :]); digits=3)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
end
end
U0e
Expand All @@ -507,10 +511,12 @@ function build_augmented_propagator_kink(; fbU::AbstractArray{ComplexF64,3}, Jw:
@info "Starting propagation within memory"
end
U0e = Array{ComplexF64}(undef, ntimes, sdim2, sdim2)
T0e = zero(U0e)
if !isnothing(output)
if !from_TTM
Utilities.check_or_insert_value(output, "U0e", U0e)
end
Utilities.check_or_insert_value(output, "T0e", T0e)
Utilities.check_or_insert_value(output, "time_taken", zeros(Float64, ntimes))
Utilities.check_or_insert_value(output, "num_paths", zeros(Int64, ntimes))
end
Expand All @@ -525,7 +531,7 @@ function build_augmented_propagator_kink(; fbU::AbstractArray{ComplexF64,3}, Jw:
_, time_taken, memory_allocated, gc_time, _ = @timed begin
forward_paths, forward_amplitudes = Utilities.generate_paths_kink_limit(forward_paths, forward_amplitudes, nkinks, sdim, fbU, extraargs.prop_cutoff, sqrt(extraargs.cutoff))
if verbose
@info "Path generation done"
@info "Path generation done. $(length(forward_paths)) forward paths generated. Expect up to $(length(forward_paths)^2) forward-backward paths based on filtration."
end
@floop exec for ((fp, fa), (bp, ba)) in Iterators.product(zip(forward_paths, forward_amplitudes), zip(forward_paths, forward_amplitudes))
num_blips = count(fp .!= bp)
Expand All @@ -548,15 +554,17 @@ function build_augmented_propagator_kink(; fbU::AbstractArray{ComplexF64,3}, Jw:
@reduce tmpU0e = zeros(ComplexF64, sdim2, sdim2) .+ tmpval
end
@inbounds U0e[i, :, :] .= tmpU0e
TTM.update_Ts!(T0e, U0e, i)
end
if !isnothing(output)
output["U0e"][i, :, :] = U0e[i, :, :]
output["T0e"][i, :, :] = T0e[i, :, :]
output["time_taken"][i] = time_taken
output["num_paths"][i] = num_paths
flush(output)
end
if verbose
@info "Done time step $(i); # paths = $(sum(num_paths)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
@info "Done time step $(i); # paths = $(sum(num_paths)); trace(Tmat) = $(round(tr(T0e[i, :, :]); digits=3)); time = $(round(time_taken; digits=3)) sec; memory allocated = $(round(memory_allocated / 1e6; digits=3)) GB; gc time = $(round(gc_time; digits=3)) sec"
end
end
U0e
Expand Down
Loading

2 comments on commit 574fa34

@amartyabose
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/115051

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.20 -m "<description of version>" 574fa348ffd77f099b7b38b9a27a73892d760f70
git push origin v0.2.20

Please sign in to comment.