Skip to content
Snippets Groups Projects
Commit c29f83bf authored by Luca Gravina's avatar Luca Gravina
Browse files

Added convergence check

parent 876bdb40
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Strucutral codes/.DS_Store
+ 0
− 0
View file @ c29f83bf
No preview for this file type
Strucutral codes/Kerr resonator/kerrdpt_order1.py 0 → 100644
+ 153
− 0
View file @ c29f83bf
import sys
sys.path.insert(1, '/Users/lucagravina/dptqm/Strucutral codes/Liouvillian block diagonalization')
import numpy as np
import matplotlib.pyplot as plt
from qutip import *
import scipy.sparse as spr
import block_diagonalize
import time
import os
import pickle
class Parameters():
"""Parameter class for Kerr resonator (used to study 1st order dissipative PT)
"""
def __init__(self, Δ, N, Grat, Utilde, ηtilde=1):
"""Initializing parameters
Args:
N (Int): Scaling factor for TDL. Morally equivalent to the number of resonators or the scaling of the non-linearity.
Delta (Float): Detuning.
Grat (Float): Two-photon drive in units of the semiclassical critical drive.
Utilde (Float): Normalized Kerr-coupling.
etatilde (Float): Normalized two-photon dissipation. Defaults to 1.
"""
self.Δ = Δ
self.N = N
self.ηtilde = ηtilde
self.Utilde = Utilde
self.U = self.Utilde/self.N
self.η = self.ηtilde/self.N
self.Gc = np.abs(self.η*self.Δ/np.sqrt(self.η**2+self.U**2))
self.G = Grat*self.Gc
self.n_sc = lambda g : self.U*self.Δ/(self.U**2+self.η**2)*(1+np.sqrt((g**2-1)*np.heaviside(g-1,1)))*np.heaviside(g-1,1) #Semiclassical approximation of the photon number
class Kerr_2γ():
def __init__(self, Nfock, p):
"""Dissipative Kerr resonator class
Args:
Nfock (Int): Hilbert space truncation.
p (Parameters): Configuration of the Parameters class.
"""
self.Nfock = Nfock
self.a = destroy(Nfock)
self.parity = 1.j*np.pi*self.a.dag()*self.a
self.parity = self.parity.expm()
self.p = p
self.c_ops = [np.sqrt(p.η)*self.a**2]
self.H = -p.Δ*self.a.dag()*self.a + p.G/2 *(self.a.dag()**2 +self.a**2) + p.U/2 *self.a.dag()**2*self.a**2
self.LL = liouvillian(self.H, self.c_ops)
def bd(self):
"""Block diagonalization of the Liouvillian matrix
Returns:
Tuple: Tuple containing
- the photon number in the two parity sectors with allowing a steady state;
- the Liouvillian gap characterizing the aformentioned symmetry sectors;
- the even-odd coherence;
"""
P, block_bfs, block_sizes = block_diagonalize.PermMat(self.LL)
self.bd_L = np.dot(P, np.dot(self.LL.data, np.transpose(P)))
self.num_blocks, self.blocks_list, self.bl_indices = block_diagonalize.get_blocks(self.bd_L)
done=False
for i in range(int(len(self.blocks_list))):
block = self.blocks_list[i]
evals, evecs = Qobj(block).eigenstates()
ss_block_form = evecs[-1]
evec2 = spr.dok_matrix((self.LL.shape[0], 1), dtype='complex')
evec2[self.bl_indices[i]:self.bl_indices[i]+self.blocks_list[i].shape[0]] = ss_block_form.data
evec2 = evec2.tocsr()
evec2 = np.dot(np.transpose(P),evec2)
ss= Qobj(evec2, dims=[self.LL.dims[0], [1]])
ss=vector_to_operator(ss)
if np.real(ss.tr())>0.05:
ss=ss+ss.dag()
ss/=(ss.tr())
if expect(self.parity, ss)>0.5:
num_even = expect(ss, self.a.dag()*self.a)
gap_even = evals[-2]
else:
num_odd = expect(ss, self.a.dag()*self.a)
gap_odd = evals[-2]
elif done==False:
gap_even_odd = np.real(evals[-1])+1.j*np.abs(np.imag(evals[-1]))
gap_odd_even = np.real(evals[-1])-1.j*np.abs(np.imag(evals[-1]))
done=True
return num_even, num_odd, gap_even, gap_odd, gap_even_odd, gap_odd_even
def get_cutoff(p, c0=5, cmax=100, step=5, precision=0.005):
"""Estimator of the an appropriate Hilbert space cutoff for a given set of parameters of the Kerr resonator. In each iteration the algorithm increments a trial cutoff by STEP starting from C0.
It compares successive evaluations of the even photon number and states convergence when their normalized difference is below a user-defined PRECISION.
Args:
p (Parameters): Configuration of the Kerr resonator
c0 (Int, optional): Lower bound for cutoff estimation. Defaults to 5.
cmax (Int, optional): Upper bound for cutoff estimation. If the selected bound does not lead to convergence, an infinite cutoff is returned. Defaults to 100.
step (Int, optional): Step-like increment of the trial cutoff in the search. Defaults to 5.
precision (Float, optional): threshold for convergence. Defaults to 0.005.
Returns:
Int: Optimal Hilbert space cutoff
"""
metric = lambda x1,x2 : ((np.abs(x1-x2)/np.abs(np.min([x1,x2])))<precision)
k = Kerr_2γ(Nfock=c0,p=p)
nprev = k.bd()[0]
c = c0+step
while c<cmax:
k = Kerr_2γ(Nfock=c,p=p)
n = k.bd()[0]
if metric(n,nprev):
c -= step
while 1:
step = np.max([step//2,1])
k = Kerr_2γ(Nfock=c-step,p=p)
if metric(k.bd()[0],n):
c -= step
elif step==1: return c
else:
nprev = n
c += step
return np.inf
def map_cutoff(Utilderange, Δrange, Grange):
"""Create map with optimal truncations for a user-defined parameter range.
Args:
Utilderange (ndarray): Array with selected values of U
Deltarange (ndarray): Array with selected values of Δ
Grange (ndarray): Array with selected values of Grat
Returns:
ndarray: Array containing the optimal cutoffs for each combination of the aforementioned parameters.
"""
C = np.zeros((Utilderange.size, Δrange.size, Grange.size),dtype='float')
for i in range(Utilderange.size):
Utilde=Utilderange[i]
for j in range(Δrange.size):
Δ=Δrange[j]
for k in range(Grange.size):
Grat=Grange[k]
p = Parameters(Δ=Δ, N=10, Grat=Grat, ηtilde=1, Utilde=Utilde)
C[i,j,k] = get_cutoff(p, c0=5, cmax=100, step=5, precision=0.005)
return C
\ No newline at end of file
Tests/.ipynb_checkpoints/test2-checkpoint.ipynb
+ 14
− 14
View file @ c29f83bf
%% Cell type:code id:466e6e45 tags:
``` python
import sys
sys.path.insert(1, '/Users/lucagravina/dptqm/Strucutral codes/Liouvillian block diagonalization')
import numpy as np
import matplotlib.pyplot as plt
from qutip import *
import scipy.sparse as spr
import block_diagonalize
import time
import os
import pickle
import matplotlib.image as mpimg
from collections import deque
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.colors as colors
greiner = { 'red': ((0., 1, 1,), (.2, 0, 0), (.48, 0, 0), (.728, 1, 1), (0.912, 1, 1), (1, .5, .5)),
'green': ((0., 1, 1), (.2, 0, 0), (.3, 0, 0), (.5, 1, 1), (.712, 1, 1), (.928, 0, 0), (1, 0, 0)),
'blue': ((0., 1, 1), (.2, .5, .5), (.288, 1, 1), (.472, 1, 1), (.72, 0, 0), (1, 0, 0)) }
greiner = LinearSegmentedColormap('greiner', greiner)
```
%% Cell type:code id:a967f6f9 tags:
``` python
```
%% Cell type:code id:4b11dd09 tags:
%% Cell type:code id:0fad0444 tags:
``` python
class Parameters():
def __init__(self, Δ, N, Grat, ηtilde, Utilde):
self.Δ = Δ
self.N = N
self.ηtilde = ηtilde
self.Utilde = Utilde
self.U = self.Utilde/self.N
self.η = self.ηtilde/self.N
self.Gc = np.abs(self.η*self.Δ/np.sqrt(self.η**2+self.U**2))
self.G = Grat*self.Gc
self.n_sc = lambda g : self.U*self.Δ/(self.U**2+self.η**2)*(1+np.sqrt((g**2-1)*np.heaviside(g-1,1)))*np.heaviside(g-1,1)
class Kerr_2γ():
def __init__(self, Nfock, p):
self.Nfock = Nfock
self.a = destroy(Nfock)
self.parity = 1.j*np.pi*self.a.dag()*self.a
self.parity = self.parity.expm()
self.p = p
self.c_ops = [np.sqrt(p.η)*self.a**2]
self.H = -p.Δ*self.a.dag()*self.a + p.G/2 *(self.a.dag()**2 +self.a**2) + p.U/2 *self.a.dag()**2*self.a**2
self.LL = liouvillian(self.H, self.c_ops)
def bd(self):
P, block_bfs, block_sizes = block_diagonalize.PermMat(self.LL)
self.bd_L = np.dot(P, np.dot(self.LL.data, np.transpose(P)))
self.num_blocks, self.blocks_list, self.bl_indices = block_diagonalize.get_blocks(self.bd_L)
done=False
for i in range(int(len(self.blocks_list))):
block = self.blocks_list[i]
evals, evecs = Qobj(block).eigenstates()
ss_block_form = evecs[-1]
evec2 = spr.dok_matrix((self.LL.shape[0], 1), dtype='complex')
evec2[self.bl_indices[i]:self.bl_indices[i]+self.blocks_list[i].shape[0]] = ss_block_form.data
evec2 = evec2.tocsr()
evec2 = np.dot(np.transpose(P),evec2)
ss= Qobj(evec2, dims=[self.LL.dims[0], [1]])
ss=vector_to_operator(ss)
if np.real(ss.tr())>0.05:
ss=ss+ss.dag()
ss/=(ss.tr())
if expect(self.parity, ss)>0.5:
num_even = expect(ss, self.a.dag()*self.a)
gap_even = evals[-2]
else:
num_odd = expect(ss, self.a.dag()*self.a)
gap_odd = evals[-2]
elif done==False:
gap_even_odd = np.real(evals[-1])+1.j*np.abs(np.imag(evals[-1]))
gap_odd_even = np.real(evals[-1])-1.j*np.abs(np.imag(evals[-1]))
done=True
return num_even, num_odd, gap_even, gap_odd, gap_even_odd, gap_odd_even
def get_cutoff(p, c0=5, cmax=100, step=5, precision=0.005):
metric = lambda x1,x2 : ((np.abs(x1-x2)/np.abs(np.min([x1,x2])))<precision)
k = Kerr_2γ(Nfock=c0,p=p)
nprev = k.bd()[0]
c = c0+step
while c<cmax:
k = Kerr_2γ(Nfock=c,p=p)
n = k.bd()[0]
if metric(n,nprev):
c -= step
while 1:
step = np.max([step//2,1])
k = Kerr_2γ(Nfock=c-step,p=p)
if metric(k.bd()[0],n):
c -= step
elif step==1: return c
else:
nprev = n
c += step
return np.inf
def map_cutoff(Utilderange, Δrange, Grange):
C = np.zeros((Utilderange.size, Δrange.size, Grange.size),dtype='float')
for i in range(Utilderange.size):
Utilde=Utilderange[i]
for j in range(Δrange.size):
Δ=Δrange[j]
for k in range(Grange.size):
Grat=Grange[k]
p = Parameters(Δ=Δ, N=10, Grat=Grat, ηtilde=1, Utilde=Utilde)
C[i,j,k] = get_cutoff(p, c0=5, cmax=100, step=5, precision=0.005)
return C
```
%% Cell type:code id:d73b1b75 tags:
%% Cell type:code id:ead02f7c tags:
``` python
```
%% Cell type:markdown id:8f886292 tags:
# Symmetry-preserving 1st order PT
In this case we look at the model
$$
\hat{H}=-\Delta \hat{a}^{\dagger} \hat{a}+\frac{U}{2}\hat{a}^{\dagger^2} \hat{a}^2 + \frac{G}{2} \left[\hat{a}^{{\dagger}^2} + \hat{a}^{\dagger^2}\right]
$$
with
$$
U=\tilde{U} / N, \quad \eta=\tilde{\eta} / N
$$
and two-photon dissipation $\eta\mathrm{D}\left[a^2\right]$.
Here we have a strong-parity symmetry which is preserved in the 1st order PT. In particular we have a 1st order PT for each of the two NESS with opposite parity: within each parity sector we have a level touching with subsequent PT.
%% Cell type:markdown id:375b8dd0 tags:
### General example
We are interested in the PT in $G$. In the following we observe for $\tilde\eta = 1$, $\tilde U/\tilde\eta = \Delta/\tilde\eta = 10$ and $N=10$ a typical 1st order PT in both NESS.
The level touching is observable in both symmetry sectors.
The question is: how does the hysteresys area change in the remaining parameters ($\Delta$, $U$)?
%% Cell type:code id:3004e4b2 tags:
``` python
Grange = np.linspace(0,2,20)
num_even = np.zeros(Grange.size)
num_odd = np.zeros(Grange.size)
gap_even = np.zeros(Grange.size, dtype='complex')
gap_odd = np.zeros(Grange.size, dtype='complex')
gap_even_odd = np.zeros(Grange.size, dtype='complex')
gap_odd_even = np.zeros(Grange.size, dtype='complex')
for i in range(Grange.size):
Grat = Grange[i]
p = Parameters(Δ=1.5, N=10, Grat=Grat, ηtilde=1, Utilde=1)
k = Kerr_2γ(Nfock=40, p=p)
num_even[i], num_odd[i], gap_even[i], gap_odd[i], gap_even_odd[i], gap_odd_even[i] = k.bd()
print("\rCompleted: %d/%d"%(i+1,Grange.size),end='')
```
%% Output
Completed: 20/20
%% Cell type:code id:bfe6d8e7 tags:
``` python
fig, ax = plt.subplots(1,3, figsize=(30,9))
g=np.linspace(Grange.min(),Grange.max(),1000)
ax[0].plot(Grange, num_even/p.N, 'bo-', ms=5, label="even")
ax[0].plot(g, p.n_sc(g)/p.N, 'g-', lw=3, label="SC")
ax[0].plot(Grange, num_odd/p.N, 'rs-', ms=5, label="odd")
ax[0].axvline(1,c='grey',alpha=0.5, ls='--', lw=3)
#ax[0].set_xscale("log")
ax[0].set_xlabel(r"$G\,/\,G_c$")
ax[0].set_ylabel(r"$\langle a^{\dagger} a\rangle\,/\,N$")
ax[0].legend()
ax[1].plot(Grange, -np.real(gap_even), 'bo-', ms=5, label="even")
ax[1].plot(Grange, -np.real(gap_odd), 'rs-', ms=5, label="odd")
ax[1].axhline(-np.real(gap_even)[0], c='k', ls='--')
ax[1].axvline(1,c='grey',alpha=0.5, ls='--', lw=3)
#ax[1].set_xscale("log")
ax[1].set_yscale("log")
ax[1].set_ylabel(r"$-\Re(\lambda_1)\,/\,\tilde\eta$")
ax[1].set_xlabel(r"$G\,/\,G_c$")
ax[1].legend()
ax[2].plot(Grange, -np.real(gap_even_odd), 'g^-')
ax[2].plot(Grange, -np.real(gap_odd_even), 'yh-')
ax[2].set_xlabel(r"$G\,/\,G_c$")
ax[2].axvline(1,c='grey',alpha=0.5, ls='--', lw=3)
#ax[2].set_xscale("log")
ax[2].set_yscale('log')
fig.tight_layout()
#fig.savefig("PL_Delta=Utilde=10_N=10.png", dpi=500)
```
%% Output
%% Cell type:code id:2d909281 tags:
``` python
```
%% Cell type:markdown id:a199be0e tags:
%% Cell type:markdown id:b1c6bee4 tags:
## Estimating convergence of Hilbert space truncation
%% Cell type:code id:971b6bfa tags:
``` python
step = 1
cuts = np.arange(10,50,step)
num_even = np.zeros(cuts.size)
num_odd = np.zeros(cuts.size)
gap_even = np.zeros(cuts.size, dtype='complex')
gap_odd = np.zeros(cuts.size, dtype='complex')
gap_even_odd = np.zeros(cuts.size, dtype='complex')
gap_odd_even = np.zeros(cuts.size, dtype='complex')
Grat=1.8
p = Parameters(Δ=1.5, N=10, Grat=Grat, ηtilde=1, Utilde=1)
t1=time.time()
for i in range(cuts.size):
k = Kerr_2γ(Nfock=cuts[i], p=p)
num_even[i], num_odd[i], gap_even[i], gap_odd[i], gap_even_odd[i], gap_odd_even[i] = k.bd()
time.time()-t1
```
%% Output
33.93870401382446
%% Cell type:code id:1b5e0401 tags:
%% Cell type:code id:ad64be81 tags:
``` python
t1 = time.time()
c = get_cutoff(p,c0=36,cmax=60,step=6,precision=0.01)
print("c=",c," time=",time.time()-t1)
```
%% Output
c= 35 time= 4.905393362045288
%% Cell type:code id:6876c0dd tags:
%% Cell type:code id:b26a8310 tags:
``` python
fig,ax = plt.subplots(1,3,figsize=(18,5))
ax[0].plot(cuts,num_even,'o-')
ax[0].plot(cuts,num_odd,'o-')
ax[0].axvline(c,c='r',ls='--',lw=2)
ax[0].set_xlabel(r"$N_{cut}$")
ax[0].set_ylabel(r"$n\,(\,N_{cut}\,)$")
ax[1].set_yscale("log")
ax[1].plot(cuts[1:],(np.gradient(num_even,cuts)*100)[1:],'o-')
ax[1].plot(cuts[1:],(np.gradient(num_odd,cuts)*100)[1:],'o-')
ax[1].axvline(c,c='r',ls='--',lw=2)
ax[1].set_xlabel(r"$N_{cut}$")
ax[1].set_ylabel(r"$\partial_{N_cut}\,n\,(\,N_{cut}\,)$")
ax[2].plot(cuts[2:],(np.gradient(np.gradient(num_odd,cuts),cuts)*100)[2:],'o-')
ax[2].axvline(c,c='r',ls='--',lw=2)
ax[2].set_xlabel(r"$N_{cut}$")
ax[2].set_ylabel(r"$\partial^2_{N_cut}\,n\,(\,N_{cut}\,)$");
```
%% Output
%% Cell type:code id:7aa3b8c2 tags:
``` python
```
%% Cell type:code id:2edda78b tags:
``` python
```
%% Cell type:code id:fbccad86 tags:
``` python
```
%% Cell type:code id:76680780 tags:
``` python
Δrange = np.linspace(0.,3.5,11)
Utilderange = np.linspace(1,1,1)
Grange = np.linspace(0.1,2,10)
```
%% Cell type:code id:2b123bbf tags:
%% Cell type:code id:f87c725e tags:
``` python
```
%% Cell type:code id:9b041c2e tags:
%% Cell type:code id:c4a80ec7 tags:
``` python
```
%% Cell type:code id:4b45cdb9 tags:
%% Cell type:code id:20246b30 tags:
``` python
```
%% Cell type:code id:44386e9c tags:
%% Cell type:code id:18c2baf9 tags:
``` python
```
%% Cell type:code id:22ab118b tags:
%% Cell type:code id:f3eca601 tags:
``` python
```
%% Cell type:code id:4d0839e6 tags:
%% Cell type:code id:1bce8098 tags:
``` python
```
%% Cell type:code id:93df735b tags:
%% Cell type:code id:d1cda65b tags:
``` python
```
%% Cell type:code id:b7171447 tags:
%% Cell type:code id:e30d772c tags:
``` python
```
%% Cell type:code id:da108646 tags:
%% Cell type:code id:ad5e5e58 tags:
``` python
```
%% Cell type:code id:53d5cb37 tags:
``` python
def get_cutoff_MK1(cuts, p, precision=0.01):
metric = lambda x1,x2 : ((np.abs(x1-x2)/np.abs(np.min([x1,x2])))<precision)
q = deque()
k = Kerr_2γ(Nfock=cuts[0], p=p)
q.append(0)
q.append(k.bd()[0])
for i in range(1,cuts.size):
k = Kerr_2γ(Nfock=cuts[i], p=p)
ni = k.bd()[0]
μi = metric(ni,q.pop())
q.append(μi)
print(q[0],q[1],i)
if (q.popleft())*(q[0])==True :
return i-1
q.append(ni)
if q[0]==True:
print("Carefull: Convergence on border")
return i
else: return np.inf
def get_cutoff_MK2(p, c0=5, cmax=100, step=5, precision=0.005):
metric = lambda x1,x2 : ((np.abs(x1-x2)/np.abs(np.min([x1,x2])))<precision)
k = Kerr_2γ(Nfock=c0,p=p)
nprev = k.bd()[0]
c = c0+step
while c<cmax:
k = Kerr_2γ(Nfock=c,p=p)
n = k.bd()[0]
if metric(n,nprev):
c -= step
while 1:
k = Kerr_2γ(Nfock=c-1,p=p)
if metric(k.bd()[0],n): c -= 1
else: return c
else:
nprev = n
c += step
return np.inf
```
%% Cell type:code id:b5385b39 tags:
``` python
```
......
Tests/test2.ipynb
+ 5
− 14
View file @ c29f83bf
Source diff could not be displayed: it is too large. Options to address this: view the blob.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment