forked from NingyiLyu/TT-TFD-GQME
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSI_TT-TFD.py
239 lines (214 loc) · 8.6 KB
/
SI_TT-TFD.py
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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Oct 17 04:11:14 2021
@author: ningyi
"""
import numpy as np
import tt
import tt.ksl
import matplotlib.pyplot as plt
import os
qmodes=60 #number of quantum bath modes
nsc = 10200 # number of propagation steps
tau = 0.00150082999505279 # propagation time step
eps = 1e-14 # tt approx error
rma = 2000000 # max tt rank
dim = qmodes # number of coords
nstates=2 # number of surfaces
occ=10 #maximum occupation number; low for harmonic systems
EYE=1j
kondo=0.1 #kondo parameter
au=1.#3.00166*10**(-4)
cmn1toau=4.5563353e-6 # Conversion of wavenumbers to atomic units
au2ps=0.00002418884254# Conversion of attoseconds to atomic units
wc=1.*au #max freq for Ohmic bath discretization
wmax=5.*au
om=wc/qmodes*(1-np.exp(-wmax/wc))
J2au=2.2937104486906*10**17 #Conversion of Joules to atomic units
kB=1.308649*10**(-23) #Boltzmann's constant, in J/K
#T=300 #Temperature, in K
kT=0.2*au #kBT, in atomic units
beta=1./kT #beta
lam=20.*cmn1toau
alpha=lam/(2*wc*np.pi)
#initialize arrays for parameters
freq=np.zeros((qmodes)) #frequency
ck=np.zeros((qmodes)) #linear electron-phonon coupling constant
gk=np.zeros((qmodes)) #ck in occupation number representation
thetak=np.zeros((qmodes)) #temperature-dependent mixing parameter in TFD
sinhthetak=np.zeros((qmodes)) #sinh(theta)
coshthetak=np.zeros((qmodes)) #cosh(theta)
for i in range(qmodes):
freq[i]=-wc*np.log(1-(i+1)*om/(wc)) # Ohmic frequency
ck[i]=np.sqrt(kondo*om)*freq[i] #Ohmic coupling constant
gk[i]=-ck[i]/np.sqrt(2*freq[i]) #Transfer ck to occ. num. representation
thetak[i]=np.arctanh(np.exp(-beta*freq[i]/2)) #theta, defined for harmonic models
sinhthetak[i]=np.sinh(thetak[i]) #sinh(theta)
coshthetak[i]=np.cosh(thetak[i]) #cosh(theta)
eelec=1.0*au #electronic state energy, in a.u.
coupling=1.0*au #electronic interstate coupling, in a.u.
#Build initial ground state
su=np.array([1,0])
sd=np.array([0,1])
e1=np.sqrt(0.5)*(su+sd)
e2=np.sqrt(0.5)*(su+EYE*sd)
tt_su=tt.tensor(su)
tt_sd=tt.tensor(sd)
tt_e1=tt.tensor(e1)
tt_e2=tt.tensor(e2)
tt_Ie=tt.eye(2,1)
gs=np.zeros((occ))
gs[0]=1.
tt_gs=tt.tensor(gs)
tt_psi0=tt_su
tt_flag=tt_psi0
for k in range(2*qmodes):#double space formation
tt_psi0=tt.kron(tt_psi0,tt_gs)
#constructing Pauli operators
px=np.array([[0,1],[1,0]])
pz=np.array([[1,0],[0,-1]])
#Build electronic site energy matrix
He=eelec*pz+coupling*px
#TT-ize that energy matrix
tt_He=tt.matrix(He)
tt_He=tt.kron(tt_He,tt.eye(occ,qmodes*2))
#Build number operator, corresponds to harmonic oscillator Hamiltonian, see notesc
numoc=np.diag(np.arange(0,occ,1))
tt_numoc=tt.eye(occ,qmodes)*0.
#Build displacement operator, corresponds to x operator in real space
eneroc=np.zeros((occ,occ))
for i in range(occ-1):
eneroc[i,i+1]=np.sqrt(i+1)
eneroc[i+1,i]=eneroc[i,i+1]
#Construct number operator as TT
for k in range(qmodes):
if k==0:
tmp=tt.kron(tt.matrix(numoc)*freq[k],tt.eye(occ,qmodes-1))
elif 0<k<qmodes-1:
tmp=tt.kron(tt.eye(occ,k-1),tt.matrix(numoc)*freq[k])
tmp=tt.kron(tmp,tt.eye(occ,qmodes-k))
else:
tmp=tt.kron(tt.eye(occ,k),tt.matrix(numoc)*freq[k])
tt_numoc=tt_numoc+tmp
tt_numoc=tt_numoc.round(eps)
tt_systemnumoc=tt.kron(tt_Ie,tt_numoc)
tt_systemnumoc=tt.kron(tt_systemnumoc,tt.eye(occ,qmodes))
#create a duplicate of number operator for the ficticious system
tt_tildenumoc=tt.kron(tt_Ie,tt.eye(occ,qmodes))
tt_tildenumoc=tt.kron(tt_tildenumoc,tt_numoc)
#initialize displacement operator
tt_energy=tt.eye(occ,qmodes)*0.
tt_tilenergy=tt.eye(occ,qmodes)*0.
for k in range(qmodes):
if k==0:
#coshtheta takes account for energy flow from real to ficticious system
#thus takes account for temperature effect
tmp=tt.kron(tt.matrix(eneroc)*gk[k]*coshthetak[k],tt.eye(occ,qmodes-1))
elif 0<k<qmodes-1:
tmp=tt.kron(tt.eye(occ,k-1),tt.matrix(eneroc)*gk[k]*coshthetak[k])
tmp=tt.kron(tmp,tt.eye(occ,qmodes-k))
else:
tmp=tt.kron(tt.eye(occ,k),tt.matrix(eneroc)*gk[k]*coshthetak[k])
tt_energy=tt_energy+tmp
tt_energy=tt_energy.round(eps)
tt_systemenergy=tt.kron(tt.matrix(pz),tt_energy)
tt_systemenergy=tt.kron(tt_systemenergy,tt.eye(occ,qmodes))
for k in range(qmodes):
if k==0:
tmp=tt.kron(tt.matrix(eneroc)*gk[k]*sinhthetak[k],tt.eye(occ,qmodes-1))
elif 0<k<qmodes-1:
tmp=tt.kron(tt.eye(occ,k-1),tt.matrix(eneroc)*gk[k]*sinhthetak[k])
tmp=tt.kron(tmp,tt.eye(occ,qmodes-k))
else:
tmp=tt.kron(tt.eye(occ,k),tt.matrix(eneroc)*gk[k]*sinhthetak[k])
tt_tilenergy=tt_tilenergy+tmp
tt_tilenergy=tt_tilenergy.round(eps)
tt_tildeenergy=tt.kron(tt.matrix(pz),tt.eye(occ,qmodes))
tt_tildeenergy=tt.kron(tt_tildeenergy,tt_tilenergy)
#Note that ficticious Harmonic oscillators carry negative sign
H=tt_He+tt_systemnumoc-tt_tildenumoc+tt_systemenergy+tt_tildeenergy
H=H.round(eps)
#Construct propagation operator, d/dt psi(t0)=A psi(t0)
A=-EYE*H
y0=tt_psi0 #Initialize wavefunction
#Heaviside functions, for selecting electronic states from overall wavefunction
tt_heavu=tt.kron(tt_su,tt.ones(occ,dim*2))
tt_heavd=tt.kron(tt_sd,tt.ones(occ,dim*2))
#Propagation time step and range
t=np.arange(0,nsc*tau,tau)
t=t#*au2ps
#Add noise, for higher rank KSL propagation
radd = 19
#radd = np.array([1,9,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,19,9,1])
radd=np.array([1,2])
radd=np.append(radd,np.repeat(2,qmodes*2-3))
radd=np.append(radd,np.array([2,1]))
#if ( radd > 0 ):
tt_rand=tt.rand(occ,qmodes*2,radd)
#tt_rand=tt_rand*10**(-140)
#for i in range(radd-1):
# tt_rand=tt_rand+tt.rand(occ,qmodes*2,1)
wxw=tt_rand.to_list(tt_rand)
for i in range(qmodes*2):
wxw[i]=wxw[i]*0.005
tt_rand=tt_rand.from_list(wxw)
#tt_rand=tt_rand*tt_rand.norm()**(-1) #Renormalize noise
tt_rand=tt.kron(tt.ones(2,1),tt_rand)
#tt_rand=tt.kron(tmp,tt_rand)
y0 = y0+tt_rand*1e-10 #Ensure noise is small
ul=np.array([[1,0],[0,0]])
ur=np.array([[0,1],[0,0]])
tt_ul=tt.matrix(ul)
tt_ur=tt.matrix(ur)
for i in range(qmodes*2):
tt_ul=tt.kron(tt_ul,tt.eye(occ,1))
tt_ur=tt.kron(tt_ur,tt.eye(occ,1))
#Initalize population arrays
psu=np.zeros((nsc))
psd=np.zeros((nsc))
coh12=np.zeros((nsc),dtype=complex)
coh21=np.zeros((nsc),dtype=complex)
ptot=np.zeros((nsc))
expec_H=np.zeros((nsc))
expec_Hu=np.zeros((nsc))
#Propagation loop
for ii in range(nsc):
y0=tt.ksl.ksl(A,y0,tau)
print(t[ii])
psu[ii]=np.abs(tt.dot(tt_heavu*y0,tt_heavu*y0))
psd[ii]=np.abs(tt.dot(tt_heavd*y0,tt_heavd*y0))
coh12[ii]=tt.dot(tt.matvec(tt_ul,y0),tt.matvec(tt_ur,y0))
coh21[ii]=tt.dot(tt.matvec(tt_ur,y0),tt.matvec(tt_ul,y0))
expec_H=tt.dot(y0,tt.matvec(H,y0))
expec_Hu=tt.dot(tt_heavu*y0,tt.matvec(H,tt_heavu*y0))
#Plot population difference
#plt.figure(dpi=600)
#plt.xlim(0.,1.)
#plt.ylim(-1.,1.)
#plt.xlabel('time(ps)')
#plt.ylabel('Populations')
#plt.plot(t,psu-psd,label='Model 1, sigma_z')
#plt.legend()
np.save('t.npy',t)
docs_dir=os.path.expanduser('./Output/TT-TFD_Output/')
if tt_flag==tt_su:
np.save(os.path.join(docs_dir,'psu_initsu.npy'),psu)
np.save(os.path.join(docs_dir,'psd_initsu.npy'),psd)
np.save(os.path.join(docs_dir,'coh12_initsu.npy'),coh12)
np.save(os.path.join(docs_dir,'coh21_initsu.npy'),coh21)
elif tt_flag==tt_sd:
np.save(os.path.join(docs_dir,'psu_initsd.npy'),psu)
np.save(os.path.join(docs_dir,'psd_initsd.npy'),psd)
np.save(os.path.join(docs_dir,'coh12_initsd.npy'),coh12)
np.save(os.path.join(docs_dir,'coh21_initsd.npy'),coh21)
elif tt_flag==tt_e1:
np.save(os.path.join(docs_dir,'psu_inite1.npy'),psu)
np.save(os.path.join(docs_dir,'psd_inite1.npy'),psd)
np.save(os.path.join(docs_dir,'coh12_inite1.npy'),coh12)
np.save(os.path.join(docs_dir,'coh21_inite1.npy'),coh21)
elif tt_flag==tt_e2:
np.save(os.path.join(docs_dir,'psu_inite2.npy'),psu)
np.save(os.path.join(docs_dir,'psd_inite2.npy'),psd)
np.save(os.path.join(docs_dir,'coh12_inite2.npy'),coh12)
np.save(os.path.join(docs_dir,'coh21_inite2.npy'),coh21)