qojulia / quantumoptics.jl Goto Github PK
View Code? Open in Web Editor NEWLibrary for the numerical simulation of closed as well as open quantum systems.
Home Page: http://qojulia.org
License: Other
Library for the numerical simulation of closed as well as open quantum systems.
Home Page: http://qojulia.org
License: Other
I'm trying to use QuantumOptics.jl to study an array of coupled resonators, with coherent pumping and dissipation.
I've always used QuTip, and I'm curious as to wherever Julia is faster or not...
Normally I build my local operators in the manybody basis by calling tensor() on a list of local operators. For example, if I want to write operator my_op
in the many body basis, for N cavities
with a Hilbert space of dimension M
I normally do:
op_list = []
for i in range(N):
op_list.append(identity(M))
op_list[0] = my_op
O_local = tensor(op_list)
Unfortunately, it seems that QuantumOptics.tensor() does not support a list as an input. How would you go about setting up something similar?
For reference:
xmin = -30
xmax = 30
b_position = PositionBasis(xmin, xmax, 100)
b_momentum = MomentumBasis(b_position)
b_positiony = PositionBasis(xmin, xmax, 100)
b_momentumy = MomentumBasis(b_positiony)
Txp = transform(b_position, b_momentum)
Tpx = transform(b_momentum, b_position)
Typ = transform(b_positiony, b_momentumy)
Tpy = transform(b_momentumy, b_positiony)
Hkinx = LazyProduct(Txp, momentum(b_momentum)^2/2, Tpx)
Hkiny = LazyProduct(Typ, momentum(b_momentumy)^2/2, Tpy)
b_comp = b_position ⊗ b_positiony
Hkinx_comp = LazyTensor(b_comp, [1, 2], [Hkinx, one(b_positiony)])
Hkiny_comp = LazyTensor(b_comp, [1, 2], [one(b_position), Hkiny])
H = LazySum(Hkinx_comp, Hkiny_comp)
Now doing
ψx = gaussianstate(b_position, -10.0, 1.0, 2)
ψy = gaussianstate(b_positiony, 0, 0.5, 2)
ψ = ψx ⊗ ψy
T = collect(linspace(0., 10.0, 50))
tout, C = timeevolution.schroedinger(T, ψ, H)
gives the error that gemm!
is not defined for LazyTensor and LazyProduct. Instead one has to work-around by passing Hfull = sparse(full(H))
instead of H
to timeevolution
.
Your test code uses x .*= y
, so you should know that in Julia 0.5 this has changed meaning to be equivalent to broadcast!(identity, x, x .* y)
, so that it mutates the x
array (see JuliaLang/julia#17510 … in Julia 0.6 the whole operation will occur in-place without temporaries). So .*
should only be used if the left-hand side is a mutable array, and you don't mind mutating it.
At first glance, this looks like it is okay for you, because you use it in psi.data[i] .*=
, where psi.data
seems like an array that intend to mutate. But if it were a problem you could always change it to *=
. For a scalar multiplication like this, in any case, *=
is probably more efficient.
Note, however, that there is currently a bug in Julia 0.5, and psi.data[i] .*=
won't work correctly until JuliaLang/julia#17546 is merged.
A more general interface to pass quantum systems in combination with others to an ODE solver (for example DifferentialEquations.jl) would be nice. This is currently a bit tedious since it requires ugly modifications of modules (see e.g. a semi-classical system where the classical differential equations are stochastic here). It's probably better to this after #92.
In terms of consistency it would be good to add this, both for a vector of rates and a matrix, using diagonaljumps
in the latter case.
Dear authors of the package, your work is fantasitc. If you have time: I don't understand very well how to use embed() and and BosonicNParticleBasis(), can you give a simple example?. Thanks in advance!
The following lines in the tracenorm function change the given dense operator rho:
QuantumOptics.jl/src/metrics.jl
Lines 22 to 25 in ca191f5
This definitely should not happen. Also just making the diagonal real might be the wrong approach. Maybe we should rethink how we want to treat non-hermitian and nearly hermitian operators. Any thoughts?
Would be interesting to see benchmarks against QuTiP 4.1. QuTiP 3.1 is a bit dated.
In the nlevelsystems branch I already implemented a first draft of a n-level basis and one function to create transition operators from one level to another. What other operators and states are useful for such systems?
For example
@define FockBasis(10)
should expand into
b_fock = FockBasis(10)
a = destroy(b_fock)
at = create(b_fock)
n = number(b_fock)
When simulating quantum trajectories via mcwf
(or similar), is there a recommended way using the API to access 1) the times at which a jump occurs, and 2) which jump operator caused the jump? (That is, the equivalent of Result.col_times
and Result.col_which
in QuTiP, respectively.)
I imagine it might be possible to obtain the jump times by passing fout(t,psi) = t
and then using display_beforeevent=true
, but would there be any ambiguity to doing the same thing but with display_afterevent=true
? Additionally, I still don't see a way to obtain which operator caused the jump without diving a little deeper into the DiffEq backend.
Any tips or workarounds would be appreciated.
The conjugate of operators
conj(H)
does not have an effect. It also does not give an error. I think this is dangerous.
In other words,
conj(H).data == conj(H.data)
it not necessarily true.
Calculating eigenstates or eigenenergies of Hermitian SparseOperators is very slow. This is due to the Julia eigs function, which seems to have a problem with the type Hermitian. One should not convert the sparse matrix, I guess.
For a Hermitian SparseOperator H
eigs(H.data, nev = 100, which = :SR)
is fast,
eigs(Hermitian(H.data), nev = 100, which = :SR)
is about 50 times slower.
Also, in the API it says that it returns the 6 lowest eigenvalues by default. But it defaults to length(basis(op)) eigenvalues. This also doesn't make to much sense for the Lanczos algorithm, which is good for calculating only "some" eigenvalues. So maybe one should rather use the eigs default.
The tag name "v0.6" is not of the appropriate SemVer form (vX.Y.Z).
cc: @david-pl
The structures Ket
and Bra
are type-unstable. The field basis
is not strictly typed. This:
basis::Basis
is ambiguous, since any Basis
can be part of a Ket
, meaning that its fields are not type-stable.
The parametrisation must happen on the basis level, e.g.
abstract type StateVector{B} where {B<:Basis} end
"""
Bra(b::Basis[, data])
Bra state defined by coefficients in respect to the basis.
"""
type Bra{B} <: StateVector{B} where {B<:Basis}
basis::B
data::Vector{Complex128}
function Bra(b::Basis, data)
if length(b) != length(data)
throw(DimensionMismatch())
end
new(b, data)
end
end
In general this may have from very big to very miniscule impact on performance (it has to be measured), However, it hinders multiple dispatch greatly, as one cannot specialize methods on the type of basis (which is required for other changes) and has to use if
branches with typeof()
checks.
This will be the first change I will do, since it will allow everything else to be done. It is a pretty huge change though so it will take time to be done properly...
Unfortunately the story with CompositeBasis
which has bases::Vector{Basis}
is quite different. The container is abstractly typed which forbids inference, but fixing this problem is quite hard. I had the same issue here: JuliaDynamics/DynamicalBilliards.jl#30 and some answers are here: https://discourse.julialang.org/t/is-there-a-way-to-forward-getindex-for-tuples-in-a-type-stable-way/2889/22
This will be a much later fix though. First type-parametrization is in order.
sortedindices.union is really union |> sort i thin, but otherwise i don't see why these methods should be in the repo.
Use 1im
and 0im
or 1.0im
and 0.0im
if you need to enforce a floating-point type (note that e.g. 2.4 + 0im
will already promote to Complex{Float64}
). The reason is that the 0.im
syntax may soon be disallowed in Julia 0.6 (JuliaLang/julia#16339).
I have installed ODE and Roots for using the Package,
I added QuantumOptics by Pkg.add("QuantumOptics")
but when I put -
julia> using QuantumOptics
INFO: Precompiling module QuantumOptics.
WARNING: deprecated syntax "[a=>b for (a,b) in c]".
Use "Dict(a=>b for (a,b) in c)" instead.
WARNING: Module ODE with uuid 16016791654596 is missing from the cache.
This may mean module ODE does not support precompilation but is imported by a module that does.
ERROR: LoadError: LoadError: Declaring __precompile__(false) is not allowed in files that are being precompiled.
in require(::Symbol) at ./loading.jl:385
in require(::Symbol) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_from_node1(::String) at ./loading.jl:488
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in include_from_node1(::String) at ./loading.jl:488
in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in macro expansion; at ./none:2 [inlined]
in anonymous at ./<missing>:?
in eval(::Module, ::Any) at ./boot.jl:234
in eval(::Module, ::Any) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
in process_options(::Base.JLOptions) at ./client.jl:239
in _start() at ./client.jl:318
in _start() at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
while loading /Users/Rakshit/.julia/v0.5/QuantumOptics/src/timeevolution_simple.jl, in expression starting on line 5
while loading /Users/Rakshit/.julia/v0.5/QuantumOptics/src/QuantumOptics.jl, in expression starting on line 45
ERROR: Failed to precompile QuantumOptics to /Users/Rakshit/.julia/lib/v0.5/QuantumOptics.ji.
in compilecache(::String) at ./loading.jl:593
in require(::Symbol) at ./loading.jl:422
in require(::Symbol) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
It shows the error. Any solutions ?
In Firefox 55 and Chrome 60 on Debian the plots from the benchmarks do not load.
I ran into an interesting issue. When creating a specific operator Julia stops and returns
Julia has stopped: null, SIGKILL
Here is a minimal code that produces this error for me:
using QuantumOptics
systemsize = 5
b = NLevelBasis(systemsize)
b_mb = ManyBodyBasis(b, bosonstates(b,systemsize))
b_total = tensor(fill(b_mb,systemsize)...)
KetSite(j) = basisstate(b_mb,j)
ProjectionOperator(i,j) = embed(b_total, i, projector(KetSite(j)))
ProjectionOperator(1,1)
Changing "systemsize" to '1' in the definition of b_mb seems to resolve the problem (which is enough for my calculations)
It seems like multiplication of sparse superoperators on sparse operators is not defined. For example,
using QuantumOptics
A = GenericBasis(2)
so = spre(identityoperator(A))
ρ = basisstate(A,1) |> dm |> sparse
so * ρ
results in
MethodError: no method matching *(::QuantumOptics.superoperators.SparseSuperOperator, ::QuantumOptics.operators_sparse.SparseOperator)
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...) at operators.jl:424
*(::QuantumOptics.states.Bra, ::QuantumOptics.operators.Operator) at C:\Users\Eric-d2\.julia\v0.6\QuantumOptics\src\operators_dense.jl:74
*(::QuantumOptics.operators_dense.DenseOperator, ::QuantumOptics.operators.Operator) at C:\Users\Eric-d2\.julia\v0.6\QuantumOptics\src\operators_dense.jl:62
...
This seems strange!
Edit: bases.multiplicable(so,ρ)
returns true
in this example as well.
Hi,
I just wanted to bring you guys to the attention a small package authored by a colleague of mine to compute the steady state of an lindblad markovian system. It basically solves the system L\rho = 0 by also enforcing the trace normalization.
It's a very efficient method, allowing for finding the steady state of quite big systems.
In our group we've been using the method for quite some time and it works very well. Maybe you'd be interested in getting in touch to integrate that in QuantumOptics...
As of now the Monte-Carlo wave function method only supports vectors of decay rates. We could change the function to accept decay rate matrices and internally diagonalize them to find a set of jump operators corresponding to a diagonal Lindblad term and then calculate the time evolution.
The method of timecorrelations.correlation
, where tspan
is omitted can later not be used in timecorrelations.correlation2spectrum
.
This is because the call to steadystate.master
with the option save_everystep=true
does not return an equidistant list of times (which is required for the discrete FFT).
currently the coherent_state function operates on a Ket in-place, and should be named with an exclamation point accordingly.
When treating correlated decay with a MCWF approach, one needs to use diagonaljumps
. Since MCWF methods are used when memory efficiency is important, it makes sense to use it in combination with lazy operators. For correlated decay, one then needs a diagonaljumps
function that can handle lazy operators. Something like:
function diagonaljumps_lazy(rates, J)
@assert length(J) == size(rates)[1] == size(rates)[2]
d, v = eig(rates)
d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)]
end
Some of the errors thrown are quite cryptic. It would be nice to catch some of these errors and provide more useful messages.
It is often practical to normalize the state vector at every time step during a time evolution with a stochastic Schrödinger equation. Since this requires switching to schroedinger_dynamic
even for otherwise linear problems, it would be nice to make this an option for the solver.
I would like to suggest a change to the way SMEs are called in QuantumOptics.jl.
Consider the H object
which takes X and \rho as arguments. The first equation is what I am calling the "linear version of H". Perhaps the default setting in "stochastic.master" and "stochastic.master_dynamic" should be with the expectation part calculated and one could imagine using flag to switch to the linear version.
I noticed that there is a factor of 0.5 in the definition of the tracenorm function. Is this on purpose? To me it seems like this was a typo from the beginning, since also the docstring states that it uses this relation. However, a more common definition would be without the 0.5. This would also be more intuitive, since with the current implementation a density matrix (or any matrix with trace 1) has a tracenorm of 0.5 instead of 1 (which is what at least I would have expected).
As mentioned in #8 a general way to create a diagonal operator with something like diagonaloperator(b, energies)
would be useful.
using QuantumOptics
bx = PositionBasis(-1,1,21)
bp = MomentumBasis(bx)
Tpx = transform(bp,bx;ket_only=true)
psi = Ket(bx)
Tpx*dm(psi) # unexpected that this works
@which Tpx*dm(psi)
@which operators.gemm!(Complex(1.), Tpx, dm(psi), Complex(0.), Tpx*dm(psi))
Currently the default setting for the type of eigenvalues to compute with the eigs function is :LM, i.e. largest magnitude (default setting of Julia eigs function). To find the low energy states it would be more apt to change it to get the smallest eigenvalues with smallest real part, i.e. option :SR.
I was thinking it might make sense to have operators subtype abstract array and implement the corresponding interface. I think that would make sense conceptually, and would make it easier to use other Julia functions with them. For example, one could then use QuantumOptics Operators in LazyArrays.jl functions immediately, and not have to worry about implementing lazy functions in QuantumOptics too.
Does that seem like a good idea?
I have read through most of the docs, at least those relevant in what I would like to achieve.
It seems like there is no support for 2 spatial directions, at least for wavepacket transport. I have tried the extremely naive:
xmin = -30
xmax = 30
Npoints = 200
b_position = PositionBasis(xmin, xmax, Npoints)
b_momentum = MomentumBasis(b_position)
ymin = -20
ymax = 20
b_positiony = PositionBasis(ymin, ymax, Npoints)
b_momentumy = MomentumBasis(b_positiony)
Txp = transform(b_position, b_momentum)
Tpx = transform(b_momentum, b_position)
Typ = transform(b_positiony, b_momentumy)
Tpy = transform(b_momentumy, b_positiony)
Hkinx = LazyProduct(Txp, momentum(b_momentum)^2/2, Tpx)
Hkiny = LazyProduct(Typ, momentum(b_momentumy)^2/2, Tpy)
H = LazySum(Hkinx, Hkiny)
but I cannot express/create a gaussianstate
in two spatial dimensions. Is this true?
P.S.: I was wondering if you have a gitter room, or some other means for users to contact you if they have questions, etc. Because github issues are not very suitable for questions.
E.g. variance(index, op, state)
is missing.
Hi there, firstly: superb work on this toolbox!
I am trying to implement some notebooks and find my way around the function structure. A couple of things:
nlevel.jl
is absent. For NLevelBasis
, (called from test_bases.jl
, apparently successfully?)If I do
julia> using QuantumOptics
julia> b=NLevelBasis(3)
ERROR: UndefVarError: NLevelBasis not defined
Perhaps I am missing something about the structure of the package, but I can't find out as documentation is not displaying. For example, it would be nice to call transition
in order to build a particular Hamiltonian... is there another way to do it?
Thanks for any help!
I was trying to implement the partial transpose and ran into some problems. It took me a while to realize that this was actually due to a convention used throughout the toolbox. Namely, the tensor product is defined in a way that is exactly opposite to what one is used to from standard notations, i.e.
b = SpinBasis(1//2)
tensor(spinup(b), spindown(b))
Ket(dim=4)
basis: [Spin(1/2) ⊗ Spin(1/2)]
0.0+0.0im
0.0+0.0im
1.0+0.0im
0.0+0.0im
even though
kron(spinup(b).data, spindown(b).data)
4-element Array{Complex{Float64},1}:
0.0+0.0im
1.0+0.0im
0.0+0.0im
0.0+0.0im
which is the usual way to write a tensor product. This is due the definition of the tensor function in the toolbox, as e.g. for states in
QuantumOptics.jl/src/states.jl
Line 90 in 9ab5d4b
Since this "switched" tensor product is consistent throughout the toolbox, no issues arise. Is there a reason, though, why it is defined like this? Shouldn't we use the more common convention in order to avoid confusions?
This bit of code
using QuantumOptics
bf = FockBasis(5)
bs = SpinBasis(1//2)
σ = sigmam(bs)
embed(bf⊗bs,1,σ)
results in
SparseOperator(dim=4x4)
basis: [Spin(1/2) ⊗ Spin(1/2)]
[2, 1] = 1.0+0.0im
[4, 3] = 1.0+0.0im
which clearly overwrites the FockBasis
. This should throw an error.
Hi, I'm trying to run a simulation which will give steady-state polarization of a two-atom system driven by classical field, in presence of correlated decay. I was trying to run steadystate.master
but it seems to be broken..
The time-evolution using timeevolution.master
seems to give proper result:
tspan = 0.0:1e-10:2e-6
rho0 = tensor(basisstate(H.basis_l, 1), dagger(basisstate(H.basis_r, 1)))
@time tout, rho_t = timeevolution.master(tspan, rho0, H, σge; rates=γ_mat);
This means that tracedistance from the initial state converges, approaching the steady state.
However, when I try using steadystate, the first error I encounter is:
I tried lowering hmin
values to 1e-10
by
@time rho_ss = steadystate.master(H, σge; rates=γ_mat, hmin=1e-10)
it runs without error, but the job seems to be running forever. Would you mind looking at steadystate.master
to check if there's any bug?
Thanks in advance.
Users as stupid as myself, may not find the documentation that you have here https://qojulia.org/documentation/index.html if they only look at the github page and not at the organiation page.
So please link the docs page in the github repo page as well :)
__
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.