SFG_search_Eg_3D.py 4.11 KB
 YesunHuang committed Jun 14, 2021 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 ``````''' 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) `````` YesunHuang committed Jun 16, 2021 42 `````` psiIdeal=tensor(qeye(Na),basis(Nb,psi0_l[2]),basis(Nc,psi0_l[1])) `````` YesunHuang committed Jun 14, 2021 43 44 45 46 47 48 `````` 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,\ `````` YesunHuang committed Jun 16, 2021 49 50 `````` b.dag()*c.dag()*b*c,\ idealState]} `````` YesunHuang committed Jun 14, 2021 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 `````` 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] `````` YesunHuang committed Jun 16, 2021 70 71 ``````Ea=np.linspace(1,10,20) g=np.linspace(0.1,4,20) `````` YesunHuang committed Jun 14, 2021 72 ``````DeltaB=0 `````` YesunHuang committed Jun 14, 2021 73 74 ``````Ntraj=2000 tlist=np.linspace(0,10,1000) `````` YesunHuang committed Jun 14, 2021 75 ``````opts=Options() `````` YesunHuang committed Jun 16, 2021 76 ``````opts.store_states=False `````` YesunHuang committed Jun 14, 2021 77 78 79 80 81 82 83 84 85 86 87 ``````#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) `````` YesunHuang committed Jun 17, 2021 88 `````` print(str(Ea[i])+'-'+str(g[j])) `````` YesunHuang committed Jun 14, 2021 89 90 91 92 93 `````` 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 `````` YesunHuang committed Jun 16, 2021 94 `````` ''' `````` YesunHuang committed Jun 14, 2021 95 96 97 98 99 100 `````` 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 `````` YesunHuang committed Jun 16, 2021 101 102 `````` ''' data[i][j][5]=output.expect[6][maxIndex] #Fidelity `````` YesunHuang committed Jun 14, 2021 103 104 105 106 107 108 109 110 `````` 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) ``````