sortedindices.union is really union |> sort i thin, but otherwise i don't see why these methods should be in the repo.

For example

@define FockBasis(10)

should expand into

b_fock = FockBasis(10)
a = destroy(b_fock)
at = create(b_fock)
n = number(b_fock)

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?

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!

__

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.

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.

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.

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

As mentioned in #8 a general way to create a diagonal operator with something like diagonaloperator(b, energies) would be useful.

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 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)

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 ?

The tag name "v0.6" is not of the appropriate SemVer form (vX.Y.Z).
cc: @david-pl

In Firefox 55 and Chrome 60 on Debian the plots from the benchmarks do not load.

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).

Would be interesting to see benchmarks against QuTiP 4.1. QuTiP 3.1 is a bit dated.

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.

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.

Some of the errors thrown are quite cryptic. It would be nice to catch some of these errors and provide more useful messages.

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...

currently the coherent_state function operates on a Ket in-place, and should be named with an exclamation point accordingly.

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).

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.

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.

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 :)

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.

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.

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).

E.g. variance(index, op, state) is missing.

The following lines in the tracenorm function change the given dense operator rho:

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?

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 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.

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.

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?

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.

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?

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:

1. the documentation link from github is (as of today) throwing a 404 error
2. Some of the basic functions don't seem to be available, even though they are part of test scripts that do complete successfully. For example: everything from 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

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?

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.