Compatibility with QuantumOpticsBase splitting

src/independent.jl

@@ -1,17 +1,17 @@
 module independent
 
-using QuantumOptics
+using QuantumOpticsBase
 using ..interaction, ..system
 
 import ..integrate
 
 # Define Spin 1/2 operators
 spinbasis = SpinBasis(1//2)
-sigmax = spin.sigmax(spinbasis)
-sigmay = spin.sigmay(spinbasis)
-sigmaz = spin.sigmaz(spinbasis)
-sigmap = spin.sigmap(spinbasis)
-sigmam = spin.sigmam(spinbasis)
+sigmax_ = sigmax(spinbasis)
+sigmay_ = sigmay(spinbasis)
+sigmaz_ = sigmaz(spinbasis)
+sigmap_ = sigmap(spinbasis)
+sigmam_ = sigmam(spinbasis)
 I_spin = identityoperator(spinbasis)
 
 
@@ -68,7 +68,7 @@ end
 Create density operator from independent sigma expectation values.
 """
 function densityoperator(sx::Number, sy::Number, sz::Number)
-    return 0.5*(identityoperator(spinbasis) + sx*sigmax + sy*sigmay + sz*sigmaz)
+    return 0.5*(identityoperator(spinbasis) + sx*sigmax_ + sy*sigmay_ + sz*sigmaz_)
 end
 function densityoperator(state::Vector{Float64})
     N = dim(state)
@@ -126,7 +126,7 @@ function timeevolution(T, gamma::Number, state0::Vector{Float64}; kwargs...)
     end
 
     fout_(t::Float64, u::Vector{Float64}) = deepcopy(u)
-    
+
     return integrate(T, f, state0, fout_; kwargs...)
 end
 

src/meanfield.jl

@@ -4,17 +4,17 @@ export ProductState, densityoperator
 
 import ..integrate
 
-using QuantumOptics, LinearAlgebra
+using QuantumOpticsBase, LinearAlgebra
 using ..interaction, ..system
 
 # Define Spin 1/2 operators
 spinbasis = SpinBasis(1//2)
 I = dense(identityoperator(spinbasis))
-sigmax = dense(spin.sigmax(spinbasis))
-sigmay = dense(spin.sigmay(spinbasis))
-sigmaz = dense(spin.sigmaz(spinbasis))
-sigmap = dense(spin.sigmap(spinbasis))
-sigmam = dense(spin.sigmam(spinbasis))
+sigmax_ = dense(sigmax(spinbasis))
+sigmay_ = dense(sigmay(spinbasis))
+sigmaz_ = dense(sigmaz(spinbasis))
+sigmap_ = dense(sigmap(spinbasis))
+sigmam_ = dense(sigmam(spinbasis))
 
 """
 Class describing a Meanfield state (Product state).
@@ -56,9 +56,9 @@ function ProductState(rho::DenseOperator)
     sx, sy, sz = splitstate(s)
     f(ind, op) = real(expect(embed(basis, ind, op), rho))
     for k=1:N
-        sx[k] = f(k, sigmax)
-        sy[k] = f(k, sigmay)
-        sz[k] = f(k, sigmaz)
+        sx[k] = f(k, sigmax_)
+        sy[k] = f(k, sigmay_)
+        sz[k] = f(k, sigmaz_)
     end
     return state
 end
@@ -122,7 +122,7 @@ splitstate(state::ProductState) = splitstate(state.N, state.data)
 Create density operator from independent sigma expectation values.
 """
 function densityoperator(sx::Real, sy::Real, sz::Real)
-    return 0.5*(I + sx*sigmax + sy*sigmay + sz*sigmaz)
+    return 0.5*(I + sx*sigmax_ + sy*sigmay_ + sz*sigmaz_)
 end
 function densityoperator(state::ProductState)
     sx, sy, sz = splitstate(state)
@@ -172,7 +172,7 @@ function timeevolution(T, S::system.SpinCollection, state0::ProductState; fout=n
     @assert N==state0.N
     Ω = interaction.OmegaMatrix(S)
     Γ = interaction.GammaMatrix(S)
-    
+
     function f(dy::Vector{Float64}, y::Vector{Float64}, p, t)
         sx, sy, sz = splitstate(N, y)
         dsx, dsy, dsz = splitstate(N, dy)
@@ -196,7 +196,7 @@ function timeevolution(T, S::system.SpinCollection, state0::ProductState; fout=n
     else
         fout_ = fout
     end
-    
+
     return integrate(T, f, state0, fout_; kwargs...)
 end
 
@@ -225,13 +225,13 @@ function timeevolution_symmetric(T, state0::ProductState, Ωeff::Real, Γeff::Re
         dsy[1] = δ0*sx[1] - Ωeff*sx[1]*sz[1] - 0.5*γ*sy[1] + 0.5*Γeff*sy[1]*sz[1]
         dsz[1] = -γ*(1+sz[1]) - 0.5*Γeff*(sx[1]^2+sy[1]^2)
     end
-    
+
     if isa(fout, Nothing)
         fout_(t::Float64, state::ProductState) = deepcopy(state)
     else
         fout_ = fout
     end
-    
+
     return integrate(T, f, state0, fout_; kwargs...)
 
 end

src/mpc.jl

@@ -1,6 +1,6 @@
 module mpc
 
-using QuantumOptics, LinearAlgebra
+using QuantumOpticsBase, LinearAlgebra
 using ..interaction, ..system, ..quantum
 
 import ..integrate
@@ -22,11 +22,11 @@ export MPCState, densityoperator
 
 
 spinbasis = SpinBasis(1//2)
-sigmax = dense(spin.sigmax(spinbasis))
-sigmay = dense(spin.sigmay(spinbasis))
-sigmaz = dense(spin.sigmaz(spinbasis))
-sigmap = dense(spin.sigmap(spinbasis))
-sigmam = dense(spin.sigmam(spinbasis))
+sigmax_ = dense(sigmax(spinbasis))
+sigmay_ = dense(sigmay(spinbasis))
+sigmaz_ = dense(sigmaz(spinbasis))
+sigmap_ = dense(sigmap(spinbasis))
+sigmam_ = dense(sigmam(spinbasis))
 
 """
 Class describing a MPC state (Product state + Correlations).
@@ -76,19 +76,19 @@ function MPCState(rho::AbstractOperator)
     sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state)
     f(ind, op) = real(expect(embed(basis, ind, op), rho))
     for k=1:N
-        sx[k] = f(k, sigmax)
-        sy[k] = f(k, sigmay)
-        sz[k] = f(k, sigmaz)
+        sx[k] = f(k, sigmax_)
+        sy[k] = f(k, sigmay_)
+        sz[k] = f(k, sigmaz_)
         for l=1:N
             if k==l
                 continue
             end
-            Cxx[k,l] = f([k,l], [sigmax, sigmax])
-            Cyy[k,l] = f([k,l], [sigmay, sigmay])
-            Czz[k,l] = f([k,l], [sigmaz, sigmaz])
-            Cxy[k,l] = f([k,l], [sigmax, sigmay])
-            Cxz[k,l] = f([k,l], [sigmax, sigmaz])
-            Cyz[k,l] = f([k,l], [sigmay, sigmaz])
+            Cxx[k,l] = f([k,l], [sigmax_, sigmax_])
+            Cyy[k,l] = f([k,l], [sigmay_, sigmay_])
+            Czz[k,l] = f([k,l], [sigmaz_, sigmaz_])
+            Cxy[k,l] = f([k,l], [sigmax_, sigmay_])
+            Cxz[k,l] = f([k,l], [sigmax_, sigmaz_])
+            Cyz[k,l] = f([k,l], [sigmay_, sigmaz_])
         end
     end
     return state
@@ -270,9 +270,9 @@ function densityoperator(state::MPCState)
     ρ = reduce(tensor, productstate)
     for k=1:N, l=k+1:N
         ρ += 0.25*(
-              Cxx[k,l]*C(sigmax,sigmax,k,l) + Cxy[l,k]*C(sigmay,sigmax,k,l) + Cxz[l,k]*C(sigmaz,sigmax,k,l)
-            + Cxy[k,l]*C(sigmax,sigmay,k,l) + Cyy[k,l]*C(sigmay,sigmay,k,l) + Cyz[l,k]*C(sigmaz,sigmay,k,l)
-            + Cxz[k,l]*C(sigmax,sigmaz,k,l) + Cyz[k,l]*C(sigmay,sigmaz,k,l) + Czz[k,l]*C(sigmaz,sigmaz,k,l))
+              Cxx[k,l]*C(sigmax_,sigmax_,k,l) + Cxy[l,k]*C(sigmay_,sigmax_,k,l) + Cxz[l,k]*C(sigmaz_,sigmax_,k,l)
+            + Cxy[k,l]*C(sigmax_,sigmay_,k,l) + Cyy[k,l]*C(sigmay_,sigmay_,k,l) + Cyz[l,k]*C(sigmaz_,sigmay_,k,l)
+            + Cxz[k,l]*C(sigmax_,sigmaz_,k,l) + Cyz[k,l]*C(sigmay_,sigmaz_,k,l) + Czz[k,l]*C(sigmaz_,sigmaz_,k,l))
     end
     return ρ
 end
@@ -429,7 +429,7 @@ function timeevolution(T, S::system.SpinCollection, state0::MPCState; fout=nothi
     else
     fout_ = fout
     end
-    
+
     return integrate(T, f, state0, fout_; kwargs...)
 end
 

src/quantum.jl

@@ -18,11 +18,11 @@ export Hamiltonian, JumpOperators
 
 # Define Spin 1/2 operators
 spinbasis = SpinBasis(1//2)
-sigmax = spin.sigmax(spinbasis)
-sigmay = spin.sigmay(spinbasis)
-sigmaz = spin.sigmaz(spinbasis)
-sigmap = spin.sigmap(spinbasis)
-sigmam = spin.sigmam(spinbasis)
+sigmax_ = sigmax(spinbasis)
+sigmay_ = sigmay(spinbasis)
+sigmaz_ = sigmaz(spinbasis)
+sigmap_ = sigmap(spinbasis)
+sigmam_ = sigmam(spinbasis)
 I_spin = identityoperator(spinbasis)
 
 """
@@ -90,15 +90,15 @@ function Hamiltonian(S::system.SpinCollection)
     H = SparseOperator(b)
     for i=1:N
         if S.spins[i].delta != 0.
-            H += 0.5*S.spins[i].delta * embed(b, i, sigmaz)
+            H += 0.5*S.spins[i].delta * embed(b, i, sigmaz_)
         end
     end
     for i=1:N, j=1:N
         if i==j
             continue
         end
-        sigmap_i = embed(b, i, sigmap)
-        sigmam_j = embed(b, j, sigmam)
+        sigmap_i = embed(b, i, sigmap_)
+        sigmam_j = embed(b, j, sigmam_)
         H += interaction.Omega(spins[i].position, spins[j].position, S.polarizations[i], S.polarizations[j], S.gammas[i], S.gammas[j])*sigmap_i*sigmam_j
     end
     return H
@@ -122,7 +122,7 @@ function Hamiltonian(S::system.CavitySpinCollection)
     at = SparseOperator(create(bc))
     for i=1:length(S.spincollection.spins)
         if S.g[i] != 0.
-            H += S.g[i]*(tensor(a, embed(bs, i, sigmap)) + tensor(at, embed(bs, i, sigmam)))
+            H += S.g[i]*(tensor(a, embed(bs, i, sigmap_)) + tensor(at, embed(bs, i, sigmam_)))
         end
     end
     return H
@@ -134,7 +134,7 @@ end
 Jump operators of the given system.
 """
 function JumpOperators(S::system.SpinCollection)
-    J = SparseOperator[embed(basis(S), i, sigmam) for i=1:length(S.spins)]
+    J = SparseOperator[embed(basis(S), i, sigmam_) for i=1:length(S.spins)]
     Γ = interaction.GammaMatrix(S)
     return Γ, J
 end
@@ -179,7 +179,7 @@ function JumpOperators_diagonal(S::system.SpinCollection)
     for i=1:N
         op = Operator(b)
         for j=1:N
-            op += M[j,i]*embed(b, j, sigmam)
+            op += M[j,i]*embed(b, j, sigmam_)
         end
         push!(J, sqrt(λ[i])*op)
     end
@@ -248,7 +248,7 @@ function rotate(axis::Vector{T1}, angles::Vector{T2}, ρ::DenseOperator) where {
     basis = ρ.basis_l
     n = axis/norm(axis)
     for i=1:N
-        nσ = n[1]*sigmax + n[2]*sigmay + n[3]*sigmaz
+        nσ = n[1]*sigmax_ + n[2]*sigmay_ + n[3]*sigmaz_
         α = angles[i]
         R = I_spin*cos(α/2) - 1im*nσ*sin(α/2)
         R_ = embed(basis, i, R)
@@ -272,7 +272,7 @@ function squeeze_sx(χT::Real, ρ₀::DenseOperator)
     N = dim(ρ₀)
     basis = ρ₀.basis_l
     totaloperator(op::SparseOperator) = sum([embed(basis, i, op) for i=1:N])/N
-    sigmax_total = totaloperator(sigmax)
+    sigmax_total = totaloperator(sigmax_)
     H = χT*sigmax_total^2
     T = [0.,1.]
     t, states = QuantumOptics.timeevolution.master(T, ρ₀, H, [])
@@ -295,7 +295,7 @@ function squeeze(axis::Vector{T}, χT::Real, ρ₀::DenseOperator) where T<:Real
     N = dim(ρ₀)
     basis = ρ₀.basis_l
     totaloperator(op::SparseOperator) = sum([embed(basis, i, op) for i=1:N])/N
-    σ = map(totaloperator, [sigmax, sigmay, sigmaz])
+    σ = map(totaloperator, [sigmax_, sigmay_, sigmaz_])
     σn = sum([axis[i]*σ[i] for i=1:3])
     H = χT*σn^2
     tout, states = QuantumOptics.timeevolution.master([0,1], ρ₀, H, [])
@@ -334,7 +334,7 @@ function squeezingparameter(ρ::DenseOperator)
     N = dim(ρ)
     basis = ρ.basis_l
     totaloperator(op::SparseOperator) = sum([embed(basis, i, op) for i=1:N])/N
-    S = map(totaloperator, [sigmax, sigmay, sigmaz])
+    S = map(totaloperator, [sigmax_, sigmay_, sigmaz_])
     n = real([expect(s, ρ) for s=S])
     e1, e2 = orthogonal_vectors(n)
     function f(phi)

src/reducedspin.jl

@@ -6,7 +6,6 @@ export ReducedSpinBasis, reducedspintransition, reducedspinstate, reducedsigmap,
 
 import Base: ==
 
-using .bases, .states, .operators, .operators_sparse
 using ..interaction, ..system