Homec4science
Diffusion BosonicCodes (master) Browse

R10106/qutip_cat_state.py
22796502c783master

Last Change
Raw File

qutip_cat_state.py
View Options

import qutip
import numpy as np
from math import floor
from scipy.optimize import minimize
from math import factorial as fac
"""----------------"""
""" initialization """
"""----------------"""
def adjacency_mat(Lx, Ly):
""""return upper triangle of 2d adjacencies without periodic boundaries"""
I_a, I_b = np.eye(Lx), np.eye(Ly)
offdi_a, offdi_b = np.diag(np.ones(Lx-1), k=1), np.diag(np.ones(Ly-1), k=1)
return np.kron(offdi_a, I_b) + np.kron(I_a, offdi_b)
def set_default_init_params():
""""global variables:"""
global N_max, N_modes, Lx, Ly
global onsite, hop, kerr, pump
print("setting default initial parameters...")
N_max, N_modes = 7, 2 # choose even number of modes!
print("N_max:", N_max, "\t N_modes:", N_modes)
# Dimensions of the reservoir lattice
Lx = 1#floor(np.sqrt(N_modes))
Ly = N_modes #floor(N_modes / Lx)
#onsite energies
onsite = np.ones(N_modes) * 0
print("onsite energies\n", onsite)
#hopping matrix with shape (N_modes, N_modes)
hop = np.zeros((N_modes, N_modes))
# adjacency_mat(Lx, Ly) \
# * np.random.uniform(-1, 1, (Lx*Ly)**2).reshape(Lx*Ly, Ly*Lx)
print("hopping matrix elements\n", hop)
#cross-Kerr matrix with shape (N_modes, N_modes)
# x_kerr = adjacency_mat(Lx, Ly) \
# * np.random.uniform(0.5, 1.5, (Lx*Ly)**2).reshape(Lx*Ly, Ly*Lx) * 10
#onsite-Kerr matrix with shape (N_modes, N_modes)
onsite_kerr = np.eye((N_modes))
kerr = onsite_kerr * 0.1
print("kerr nonlinearity elements\n", kerr)
#single-particle pumping amplitudes
pump = np.ones(N_modes) # np.random.uniform(- 1, 1, N_modes) * 1
print("pump\n", pump)
print("done.")
"""----------------"""
def spectral_rad(lx, ly, onsite, hop):
"calculate spectral radius of single particle hamiltonian"
mat = np.diag(onsite) + hop
print("coupling matrix:")
print(mat)
rad = min((np.linalg.eigvals(mat)))
print("spectral radius:", rad)
return abs(rad)
set_default_init_params()
print("constructing operators...")
# define single-particle operators(1.7358564808935595+2.1726903683902553e-18j)
a = qutip.destroy(N_max)
a_dag = qutip.create(N_max)
n = qutip.num(N_max)
# Define annihilation and creation operator for entire Hilbert space
# initialize two empty lists for creat. and annih. & numb ops for each mode
aa, aa_dag, nn = [], [], []
for i in np.arange(N_modes):
a_tensor = qutip.tensor([qutip.qeye(N_max)] * i
+ [a]
+ [qutip.qeye(N_max)] * (N_modes - i - 1))
nn.append(qutip.tensor([qutip.qeye(N_max)] * i
+ [n]
+ [qutip.qeye(N_max)] * (N_modes - i - 1)))
aa.append(a_tensor)
aa_dag.append(a_tensor.dag())
del a_tensor # free memory
print("done.")
def hamiltonian(n_max, n_modes, lx, ly, onsite, hop, kerr, pump):
# initialize empty Hamilton operatornp.linspace(0.1, 10, 100)]
h = qutip.Qobj(dims=[[n_max] * n_modes, [n_max] * n_modes])
for i in range(n_modes):
h += onsite[i] * nn[i] # onsite-energies
# pump
h += pump[i] * aa_dag[i] + np.conj(pump[i]) * aa[i]
for j in np.arange(i, n_modes): # ensure that j > i
print("i,j = ", i,j)
if hop[i,j] != 0: # avoid cases with zero element
h += hop[i,j] * ( aa_dag[i] * aa[j] + aa_dag[j] * aa[i])
if kerr[i,j] !=0: # avoid cases with zero element
h += kerr[i,j] * (aa_dag[i] * aa[i] * aa_dag[j] * aa[j]) #* nn[i] * nn[j]
return h
def liouvillian(h, rho, gam):
l = -1j * (h * rho - rho * h)
for i in range(N_modes):
l += gam * (aa[i] * rho * aa_dag[i] - nn[i] * rho - rho * nn[i])
return l
def rho_out(rho, w):
aa_out = aa_mixed(w)
rho_final = qutip.Qobj(dims=[[N_max**N_modes], [N_max**N_modes]])
vac = qutip.tensor(qutip.fock(N_max, 0), qutip.fock(N_max, 0))
mat = np.empty((N_max**N_modes, N_max**N_modes), dtype=complex)
for n1 in range(N_max**N_modes):
for n2 in range(N_max**N_modes):
rho_op = (aa_out**n1) * rho * (aa_out.dag() **n2)
facul = fac(n1)*fac(n2)
if facul > 1e16:
mat[n1, n2] = 0
else:
mat[n1, n2] = qutip.expect(rho_op, vac) / np.sqrt(facul)
rho_final.set_data(qutip.fastsparse.csr2fast(qutip.fastsparse.csr_matrix(mat)))
rho_final /= rho_final.tr()
return rho_final
# create Schrödinge^r's cat density matrix state
# zeta = 1.4
# rho_cat = 0.5 * qutip.coherent_dm(N_max, zeta) \
# + 0.5 * qutip.coherent_dm(N_max, -zeta)
print("constructing target quantum state...")
w1 = 1/np.sqrt(2)* ( qutip.fock(N_max, 0) + qutip.fock(N_max, 4))
w1_dm = qutip.ket2dm(w1)
w2_dm = qutip.fock_dm(N_max**N_modes, 1)
rho_target = w2_dm #w1_dm/2 + w2_dm/2
print("done.")
# rho_single_ph = qutip.fock_dm(N_max, 1)
# # qutip.plot_wigner(rho_cat)
def aa_mixed(w):
aa_out = qutip.Qobj(dims=[[N_max] * N_modes, [N_max] * N_modes])
w /= np.sqrt(sum(abs(w)**2)) # normalize weights to unit abs value
c = w.reshape((N_modes, 2))
coeffs = c[:,0] + 1j * c[:,1]
for i in range(N_modes):
aa_out += coeffs[i] * aa[i]
return aa_out
def error(w, rho1, rho_target):
"""cost function to minimize. currently: tracedistance"""
rho = rho_out(rho1, w) # mixing the output modes
res = 1 - qutip.fidelity(rho, rho_target)# qutip.tracedist(rho, rho_target) # calculate the tracedistance
return res
vals = []
def inter_res(arg):
"""Save intermediate results of minimization routine to a list"""
c = arg.reshape((N_modes, 2))
coeffs = c[:,0] + 1j * c[:,1]
err = error(arg, rho_steady, rho_target)
val = coeffs.tolist() + [err]
vals.append(val)
print("current error:", err)
def constraint(w):
"define constraints for the absolut values of the mixing coeffs"
return sum(abs(w)**2) - 1 # returns zero if constraint is satisfied
def occupation_const(w):
aa_out = aa_mixed(w)
n = (aa_out.dag() * aa_out * rho_steady).tr()
return abs(n)-1
"""minimazation"""
print("Constructing Hamiltonian...")
set_default_init_params()
mat = np.zeros((100,100))
for i, U in enumerate(np.linspace(0.1, 10, 100)):
for j, F in enumerate(np.linspace(0.1, 10, 100)):
print("U:", U, "F:", F)
h = hamiltonian(N_max, N_modes, Lx, Ly, onsite, hop, kerr*U, pump*F)
print("done.")
print("Solving for the steady-state...")
rho_steady = qutip.steadystate(h, aa, method='direct')
print("done.")
print("starting optimization routine...")
print("initialize random mixing of modes.")
x0 = np.random.uniform(0, 1, (N_modes * 2))
x0/= np.sqrt(sum(abs(x0)**2))
res = minimize(error, x0 = x0, args=(rho_steady, rho_target), method='Nelder-Mead',
options={'disp': True})#, callback=inter_res)
# constraints = {"type":"eq", "fun": occupation_const})
mat[i,j] = error(res.x, rho_steady, rho_target)
# vals = np.array(vals)
# print("kerr: ", new_kerr)
print("all done.")
c4science · Help