Commit fc703d73 authored by YesunHuang's avatar YesunHuang
Browse files

Update 3D

parent 17cf66b2
......@@ -66,8 +66,8 @@ def SearchForMax(array:list):
#parameters
kappa_a=2;kappa_b=2;kappa_c=2
psi0_l=[0,3,0]
Ea=1
g=np.linspace(0.1,4,10)
Ea=10
g=np.linspace(0.1,4,40)
Ntraj=100
DeltaB=0
tlist=np.linspace(0,20,2000)
......
'''
Name: SFG_search.py
Desriptption:
Email: yesunhuang@mail.ustc.edu.cn
OpenSource: https://github.com/yesunhuang
Msg: For super computer
Author: YesunHuang
Date: 2021-06-05 21:15:28
'''
#import area
from qutip import*
import numpy as np
from scipy import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib.pyplot import MultipleLocator
import time
'''
@name: BuildOperator_Exact
@fuction: Helper function for mesolve
@param {float} Ea
@param {float} deltaB
@param {float} g
@return {qobj} operator
'''
def BuildOperator_Exact(Ea:float,DeltaB:float,g:float):
Na=int(max(math.ceil(Ea*Ea+6*Ea),4)); Nb=2*psi0_l[1]; Nc=Nb
psi0=tensor(basis(Na,psi0_l[0]),basis(Nb,psi0_l[1]),basis(Nc,psi0_l[2]))
a=tensor(destroy(Na),qeye(Nb),qeye(Nc))
b=tensor(qeye(Na),destroy(Nb),qeye(Nc))
c=tensor(qeye(Na),qeye(Nb),destroy(Nc))
H=g*(a.dag()*b*c.dag()+a*b.dag()*c)+Ea*(a.dag()+a)+DeltaB*b.dag()*b
c_ops=[]
c_ops.append(np.sqrt(kappa_a)*a)
c_ops.append(np.sqrt(kappa_b)*b)
c_ops.append(np.sqrt(kappa_c)*c)
psiIdeal=tensor(basis(Nb,psi0_l[2]),basis(Nc,psi0_l[1]))
idealState=psiIdeal*psiIdeal.dag()
operator={'Hamilton':H,'Collapse':c_ops,'Initial_state':psi0,'Ideal_state':idealState,'track':[b.dag()*b,\
c.dag()*c,\
a.dag()*a,\
b.dag()*b.dag()*b*b,\
c.dag()*c.dag()*c*c,\
b.dag()*c.dag()*b*c]}
return operator
'''
@name: SearchForMax
@fuction: Search for max
@param {list} array
@return {int} maxIndex
'''
def SearchForMax(array:list):
max=0
maxIndex=0
for i in range(len(array)):
if array[i]>max:
max=array[i]
maxIndex=i
return maxIndex
#parameters
kappa_a=2;kappa_b=2;kappa_c=2
psi0_l=[0,3,0]
Ea=np.linspace(0.1,1,3)
g=np.linspace(0.1,4,3)
DeltaB=0
Ntraj=100
tlist=np.linspace(0,20,2000)
opts=Options()
opts.store_states=True
#opts.rhs_reuse=True
#data storage
data=np.zeros([np.size(Ea),np.size(g),9])
#solve for data
if __name__=="__main__":
for i in range(0,np.size(Ea)):
for j in range(0,np.size(g)):
op=BuildOperator_Exact(Ea[i],DeltaB,g[j])
output=mcsolve(op['Hamilton'],op['Initial_state'],tlist,op['Collapse'],op['track'],ntraj=Ntraj,options=opts)
maxIndex=SearchForMax(output.expect[1])
data[i][j][0]=tlist[maxIndex] #t
data[i][j][1]=output.expect[0][maxIndex] #Nb
data[i][j][2]=output.expect[1][maxIndex] #Nc
data[i][j][3]=output.expect[2][maxIndex] #Na
rho=output.states[0][maxIndex]*output.states[0][maxIndex].dag()
for k in range(1,Ntraj):
rho+=output.states[k][maxIndex]*output.states[k][maxIndex].dag()
rho=1/Ntraj*rho
rho_bc=rho.ptrace([1,2])
data[i][j][4]=entropy_vn(rho_bc) #(A,BC) Entanglement
data[i][j][5]=expect(op['Ideal_state'],rho_bc) #Fidelity
data[i][j][6]=output.expect[3][maxIndex]/(data[i][j][1]*data[i][j][1]) #g2b
data[i][j][7]=output.expect[4][maxIndex]/(data[i][j][2]*data[i][j][2]) #g2c
data[i][j][8]=output.expect[5][maxIndex]/(data[i][j][1]*data[i][j][2]) #g2bc
#save data
#np.savetxt('Data/population_g_'+str(g)+'_Ea_'+str(Ea[0])+'-'+str(Ea[-1])+'_DeltaB_'+str(DeltaB[0])+'-'+str(DeltaB[-1])+'.txt',data)
np.save('Data/EgData3D_DeltaB_'+str(DeltaB)+'_Ea_'+str(Ea[0])+'-'+str(Ea[-1])+'_g_'+str(g[0])+'-'+str(g[-1])+'.npy',data)
import numpy as np
from scipy import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib.pyplot import MultipleLocator
import time
g_s=0.1;g_e=4.0
DeltaB=0
Ea_s=0.1;Ea_e=1.0
data=np.load('Data/EgData3D_DeltaB_'+str(DeltaB)+'_Ea_'+str(Ea_s)+'-'+str(Ea_e)+'_g_'+str(g_s)+'-'+str(g_e)+'.npy')
g=np.linspace(g_s,g_e,np.size(data,1))
Ea=np.linspace(Ea_s,Ea_e,np.size(data,0))
#print population
fig = plt.figure(figsize=(6,18))
axes=[]
axes.append(fig.add_subplot(3, 1, 1))
axes.append(fig.add_subplot(3, 1, 2))
axes.append(fig.add_subplot(3, 1, 3))
ax0=axes[0].contourf(g,Ea,data[:,:,0],cmap=cm.coolwarm)
ax1=axes[1].contourf(g,Ea,data[:,:,1],cmap=cm.coolwarm)
ax3=axes[2].contourf(g,Ea,data[:,:,2],cmap=cm.coolwarm)
ax00=axes[0].contour(g,Ea,data[:,:,0],ax0.levels)
ax11=axes[1].contour(g,Ea,data[:,:,1],ax1.levels)
ax11=axes[2].contour(g,Ea,data[:,:,2],ax1.levels)
#fmt = '%1.2f '
axes[0].clabel(ax00, ax00.levels, inline=True)
axes[1].clabel(ax11, ax11.levels, inline=True)
axes[2].clabel(ax11, ax11.levels, inline=True)
axes[0].set_xlabel(r'$g$');axes[1].set_xlabel(r'$g$');axes[2].set_xlabel(r'$g$')
axes[0].set_ylabel(r'$E_a$');axes[1].set_ylabel(r'$E_a$');axes[2].set_ylabel(r'$E_a$')
bar1=fig.colorbar(ax0,ax=axes[0],pad=0.1)
bar2=fig.colorbar(ax1,ax=axes[1],pad=0.1)
bar2=fig.colorbar(ax1,ax=axes[2],pad=0.1)
fig.savefig('imgs/EgPopulation3D_DeltaB_'+str(DeltaB)+'_Ea_'+str(Ea_s)+'-'+str(Ea_e)+'_g_'+str(g_s)+'-'+str(g_e)+'.svg',dpi=600,format='svg',bbox_inches='tight')
#print other data
fig = plt.figure(figsize=(12,18))
axes=[]
for i in range(0,6):
axes.append(fig.add_subplot(3, 2, i+1))
ax_s=axes[i].contourf(g,Ea,data[:,:,(i+4)%9],cmap=cm.coolwarm)
ax_ss=axes[i].contour(g,Ea,data[:,:,(i+4)%9],ax_s.levels)
axes[i].clabel(ax_ss, ax_ss.levels, inline=True)
axes[i].set_xlabel(r'$g$')
axes[i].set_ylabel(r'$E_a$')
fig.colorbar(ax_s,ax=axes[i],pad=0.1)
axes[0].set_ylabel(r'$[Entanglement]E_a$')
axes[1].set_ylabel(r'$[Fidelity]E_a$')
axes[2].set_ylabel(r'$[g_2b]E_a$')
axes[3].set_ylabel(r'$[g_2c]E_a$')
axes[4].set_ylabel(r'$[g_2bc]E_a$')
axes[5].set_ylabel(r'$[t_max]E_a$')
fig.savefig('imgs/EgOtherData3D_DeltaB_'+str(DeltaB)+'_Ea_'+str(Ea_s)+'-'+str(Ea_e)+'_g_'+str(g_s)+'-'+str(g_e)+'.svg',dpi=600,format='svg',bbox_inches='tight')
......@@ -44,8 +44,8 @@ fig = plt.figure(figsize=(12,18))
axes=[]
for i in range(0,6):
axes.append(fig.add_subplot(3, 2, i+1))
ax_s=axes[i].contourf(DeltaB,Ea,data[:,:,(i+3)%9],cmap=cm.coolwarm)
ax_ss=axes[i].contour(DeltaB,Ea,data[:,:,(i+3)%9],ax_s.levels)
ax_s=axes[i].contourf(DeltaB,Ea,data[:,:,(i+4)%9],cmap=cm.coolwarm)
ax_ss=axes[i].contour(DeltaB,Ea,data[:,:,(i+4)%9],ax_s.levels)
axes[i].clabel(ax_ss, ax_ss.levels, inline=True)
axes[i].set_xlabel(r'$\Delta_b$')
axes[0].set_ylabel(r'$E_a$')
......
Markdown is supported
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