Simulate coherent multidimensional spectroscopy signals from quantum mechanical models.
Note that examples are currently not working, due to transition from CMDS.jl to MultidimensionalSpectroscopy.jl
Working examples:
- 01_coupled_dimer
- 02_displaced_harmonic_oscillator
The code relies primarily on qojulia/QuantumOptics.jl, which is well described here. Further helpful examples and functionalities are found the Python project QuTiP.
The module MultidimensionalSpectroscopy.jl contains the necessary functions to calculate 2D spectra from QM models and will be described below. examples/ shows example scenarios.
MultidimensionalSpectroscopy.jl requires the Julia language and qojulia/QuantumOptics.jl, which can be installed from via the standard sources:
Clone the repo ... and include the module via using MultidimensionalSpectroscopy
.
Type ?<function>
into the REPL to access the documentation for a certain function.
-
create_colormap: creates a blue-white-green-red colormap with zero values being white (
scheme="bright"
) or dark (scheme="dark"
) -
zeropad: zeropadding of time domain data before Fourier transformation into spectral domain
-
interpt: interpolate time vector after zeropadding
-
make2Dspectra: invokes cmds.correlations to calculate different signals
-
correlations: calculate time evolution and interaction with laser pulses
-
view_dm_evo: quick way to visualize evolution of density matrix
-
save_2d: saves 2D plots to folder
fn_base
-
plot2d: plots 2D data in out2d
-
plot_timeTrace: plots time traces from T evolution of 2D spectrum
-
crop2d: crops 2D data to reduce size
-
tri: select upper or lower triangular matrix, used to simplify pathways
-
create_subspace: create subspace (ground state -> single excitation sector, ground state -> double excitation sector, ...) Hamiltonian and operators
-
absorptionSpectrum: shortcut to calculate and plot absorption spectrum
-
plot2d_comps: plot components (GSB, SE, ESA) of 2D spectrum
-
vib_analysis: plot rephasing and non-rephasing parts of GSB, SE and ESA 2D spectra -> beating maps
-
plot_levels: diagonalise Hamiltonian and plot level structure
Set up your QM model of interest!
Example (see displaced harmonic oscillator, as well as documentation of QuantumOptics.jl):
b_tls = NLevelBasis(2) # Hilbert-space of system Basis: {|ui⟩}
b_vib = FockBasis(5) # Hilbert-space of oscillator Basis: {|vi⟩}
b = b_tls ⊗ b_vib # combined basis
j21 = transition(b_tls,2,1) # |e⟩⟨g|
j12 = dagger(j21) # |g⟩⟨e|
at = create(b_vib) # ...
a = dagger(at)
D = displace(b_vib,d) # diplacement operator with the displacement d
H = ...
Use plot_levels()
to get an overview of the energy level structure (Eigenstates) of the full systems and its components.
plot_levels(H,0) # 0: shift in x-direction (change to display multiple energy diagrams next to each other)
Result:
Operators can be ordered by increasing energy (diagonal elements of Hamiltonian) by
idx = sortperm(real(diag((H).data)))
H.data = H.data[idx,idx]
...
This is required to convert between eigen- and site basis (or coupled- and uncoupled basis) when working with a subspace of the full Hilbert space, by using the transform operator (transf_op
) from the output of create_subspace()
(see below):
H, transf_op, P, L, rho0, ... = create_subspace(H,"bi", L, rho0, ...)
op_site = transf_op * op * transf_op'
If the dimension of the basis is too large, create_subspace()
can be used to project all necessary operators into the ground + single excitation subspace (e.g. for linear spectra), or the ground + single + double excitation subspace (e.g. for 2D spectra). NOTE that creating a subspace requires ordering the operators into excitation sectors
H, transf_op, P, L, ... = create_subspace([H],"bi", L, ...)
...
replaces as maybe operators is desired.
To calculate 2D spectra first initialize the output array
spectra2d = Array{out2d}(undef, length(T))
where T = [0., 5., 10., ...]
is a vector containing population/evolution time steps.
Next call make2dspectra()
in a for loop
for i = 1:length(T)
spectra2d[i] = make2Dspectra(tlist,rho0,H,F,μ12,μ23,T[i],"lindblad";debug=false,zp=zp)
end
with tlist
being the coherence/detection time steps, rho0
the equilibrium/ground state density matrix, H
the system Hamiltonian, F
the jump operator (Lindblad operator) or the Redfield tensor, μ12
the transition dipole operator between the ground and single excited states, and μ23
the transition dipole operator between the single and double excited states. T[i]
is the current population time. The option "lindblad"
or "redfield"
selects which ... to use. debug=true
activates the debugging output and zp
is the zero padding value of 10zp.
Using multithreading, several population time steps can be evaluated simultaneously:
Threads.@threads for i = 1:length(T)
spectra2d[i] = make2Dspectra(...)
end
Make sure to disable all output plots within make2Dspectra()
when using multithreading, as these might crash the execution.
JLD2 can be used to conveniently save the spectra2d
structure (does not work with cloud drives, such as pCloud). Remember to round2d()
the data to save disk space.
@save "C:\\path\\to\\data\\file.jld2" out2d
Data can then be load as
@load "C:\\path\\to\\data\\file.jld2" out2d
However, the data format is not compatible with other software. save2d()
saves ASCII files for real (.re) and imaginary (.im) parts.
...
You can create a slider to flip through 2D spectra (work in progress):
using Blink, Interactive
mp = @manipulate for i in slider(1:length(spectra2d))
clf(); plot2d(spectra2d[i].ω,spectra2d[i].full2d)
end
w = Window();
body!(w, mp);
The following examples (scripts) are available:
Using view_dm_evo()
you can visualize the temporal (T) evolution of the system density matrix.
MultidimensionalSpectroscopy.jl outputs the full2d spectrum, as well as the GSB (spectra2d.gsb), SE (spectra2d.se) and ESA (spectra2d.esa) components. These can be conveniently visualized using plot2d_comps(spectra2d[1])
:
As you can see, the absorptive 2D data is composed of overlapping components, which can partially cancel out. One specific example is a hetero dimer without coupling: In this case the cross-peaks are absent from the absorptive 2D data, due to the cancellation of positive GSB with negative ESA.
In addition, also the rephasing (spectra2d.full2d_r) and non-rephasing (spectra2d.full2d_nr) parts of the signal can be plotted. Comparison with the absorptive representation shows that the latter exhibits sharper features in the 2D spectrum.