SFG_search_Eg_3D.py 4.11 KB
Newer Older
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
42
    psiIdeal=tensor(qeye(Na),basis(Nb,psi0_l[2]),basis(Nc,psi0_l[1]))
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
49
50
                                                                                                   b.dag()*c.dag()*b*c,\
                                                                                                   idealState]}
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
70
71
Ea=np.linspace(1,10,20)
g=np.linspace(0.1,4,20)
YesunHuang's avatar
YesunHuang committed
72
DeltaB=0
YesunHuang's avatar
YesunHuang committed
73
74
Ntraj=2000
tlist=np.linspace(0,10,1000)
YesunHuang's avatar
YesunHuang committed
75
opts=Options()
YesunHuang's avatar
YesunHuang committed
76
opts.store_states=False
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
88
            print(str(Ea[i])+'-'+str(g[j]))
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
94
            '''
YesunHuang's avatar
YesunHuang committed
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's avatar
YesunHuang committed
101
102
            '''
            data[i][j][5]=output.expect[6][maxIndex] #Fidelity
YesunHuang's avatar
YesunHuang committed
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)