-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathspace_q_fc.f90
More file actions
327 lines (227 loc) · 8.27 KB
/
space_q_fc.f90
File metadata and controls
327 lines (227 loc) · 8.27 KB
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
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
subroutine Qspace_fc()
use global
use operator_radial
use density
use wfunction
use operator_spatial_space
use operator_3d
use twoe_basis_set
implicit none
integer :: ij,il,ik,ii,ip,ix,iy,iip,icoloum,i_spatial
complex(kind=k2) :: cs,cff,csum
complex(kind=k2),dimension(fedvr3d%nb_r*fedvr3d%nb_angle,system%nptot) :: cbox1,cbox2
integer :: idicp,jdicp
integer :: indexp,k1p,k2p,index_here,iy_here
do i_spatial =1,system%nptot
!
! not consider the stiffness problem in the calculations
!
if(prop%stiffness==0) then
call act(phi(:,i_spatial),cbox1(:,i_spatial))
else
call act2(phi(:,i_spatial),cbox1(:,i_spatial))
endif
enddo
do ix =1,fedvr3d%nb_r*fedvr3d%nb_angle
do ii =1,np_tot_fc
cs = zzero
do ij =1,system%nptot
cs = cs + cden1(ii+n_offset,ij)*cbox1(ix,ij)
enddo
cbox2(ix,ii) = cs
enddo
enddo
do ix=1,fedvr3d%nb_r*fedvr3d%nb_angle
do ii =1,nptot_fc
cs=zzero
do il=1,system%nptot
do ik=1,system%nptot
do ij=1,system%nptot
do iy = 1,fedvr3d%nb_angle
iy_here = index_two(2,ix,iy)
cs=cs+cden2(ii,ik,il,ij)*cww(ix,iy,ik,il)*phi(iy_here,ij)
enddo
enddo
enddo
enddo
cbox3(ix,i_spatial)= cs + cbox2(ix,i_spatial)
enddo
enddo
!
do ix=1,fedvr3d%nb_r*fedvr3d%nb_angle
do ii=1,nptot_fc
cs = zzero
do ij =1,nptot_fc
cs = cs + cden1_reg_inv_fc(ii,ij)*cbox3(ix,ij)
enddo
cbox1(ix,ii) = cs
enddo
enddo
!
! the last step, projcetor acting
!
call projector_Q_hf(cbox1)
return
end subroutine Qspace_fc
!#######################################################################################
! projector operator in Q space, i is spatial orbital, n is the num of spatial orbital
! _N___ __N____
! \ \ _N___
! \ \\
! 1 - / | i > <i| == 1 - //___ |i > [O]_ij <j|
! /___ /____
! i=1 i,j=1
! O matrix is the inverse of <i|j> matrix, the reference can be found in
! M.H. Beck, Physics Report, 324, 2000, Page 25
!
!########################################################################################
subroutine projector_Q_hf(cbox1)
use global
use wfunction
use solver
implicit none
integer :: idicp,jdicp,kdicp,ldicp
integer,dimension(system%nptot) :: ipiv2
complex(8),dimension(system%nptot,system%nptot) :: phij,inv_phij
complex(8),dimension(fedvr3d%nb_r*fedvr3d%nb_angle,fedvr3d%nb_r*fedvr3d%nb_angle) :: proj_Q
complex(8),dimension(fedvr3d%nb_r*fedvr3d%nb_angle,system%nptot) :: cbox1,cbox3
complex(8) :: sumtemp
complex(8),dimension(fedvr3d%nb_r*fedvr3d%nb_angle) :: temp
integer :: info
!
! construct the matrix phij = < i | j >, i,j spatial orbital
!
do idicp =1,system%nptot
do jdicp =1, system%nptot
sumtemp = zzero
do kdicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
sumtemp = sumtemp + dconjg(phi(kdicp,idicp))*phi(kdicp,jdicp)
enddo
phij(idicp,jdicp) = sumtemp
inv_phij(idicp,jdicp) = zzero
enddo
inv_phij(idicp,idicp) = zone
enddo
!
! find the inverse matrix of phij
!
call zgesv(system%nptot,system%nptot,phij,system%nptot,ipiv2,inv_phij,system%nptot,info)
if(info.ne.0) then
write(*,*) '____________________________________________________________'
write(*,*) 'fail to find the inverse of phij <phi | phj> in projector_Q'
write(*,*) '------------------------------------------------------------'
stop
endif
!
! construct the Q.space Projector operator
!
proj_Q = zzero
do idicp =1,system%nptot
do jdicp =1,system%nptot
do kdicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
do ldicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
proj_Q(kdicp,ldicp) = proj_Q(kdicp,ldicp) + phi(kdicp,idicp)*inv_phij(idicp,jdicp)*dconjg(phi(ldicp,jdicp))
enddo
enddo
enddo
enddo
!
! projector operator act on the intermediate wavefunction
!
do idicp =1,nptot_fc
do jdicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
sumtemp = zzero
do kdicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
sumtemp = sumtemp + proj_Q(jdicp,kdicp)*cbox1(kdicp,idicp)
enddo
temp(jdicp) = sumtemp
enddo
!
! 1 - Projector
!
do ldicp =1,fedvr3d%nb_r*fedvr3d%nb_angle
cbox3(ldicp,idicp) = (-ci*(cbox1(ldicp,idicp) - temp(ldicp)))
enddo
enddo
do ip=1,nptot_fc
do ix=1,fedvr3d%nb_r*fedvr3d%nb_angle
kphi(ix,ip+n_offset) = cbox3(ix,ip)
enddo
enddo
return
end subroutine projector_Q_hf
!===============================================================================================
! not considered the stiffness effect, the matrix is hightly sparse
!===============================================================================================
subroutine act_fc(phi_in,phi_out)
use global
use operator_radial
use operator_3d
implicit none
integer :: idicp,jdicp,kdicp,ldicp
complex(kind=k2) :: phi_in(fedvr3d%nb_r,fedvr3d%nb_angle),phi_out(fedvr3d%nb_r,fedvr3d%nb_angle)
complex(kind=k2) :: temp1,temp2
integer :: mdicp,ndicp
do idicp =1,fedvr3d%nb_angle
do jdicp =1,fedvr3d%nb_r
temp1 = phi_in(jdicp,idicp)*(vmat_radial(jdicp) + tmat_3d(jdicp,idicp))
temp2 = zzero
do kdicp =index_ham_act(jdicp,1),index_ham_act(jdicp,2)
ldicp = index_kinetic_operator(max(jdicp,kdicp), min(jdicp,kdicp))
temp2 = temp2 + tmat_radial(ldicp)*phi_in(kdicp,idicp)
enddo
phi_out(jdicp,idicp) = temp1 + temp2
enddo
enddo
!=================================================================================================
! the laser interaction term
!=================================================================================================
if(laser%tdornot) then
do idicp=1,fedvr3d%nb_angle
do jdicp =1,fedvr3d%nb_r
temp1 = zzero
do ldicp =1,fedvr3d%nb_angle
temp1 = temp1 + laser%ez_t*zmat_3d(idicp,ldicp)*fedvrx_global(jdicp)*phi_in(jdicp,ldicp)
enddo
phi_out(jdicp,idicp) = phi_out(jdicp,idicp) + temp1
enddo
enddo
endif
return
end subroutine act_fc
!===============================================================================================
! considered the stiffness effect, the matrix is full with element
!===============================================================================================
subroutine act2_fc(phi_in,phi_out)
use global
use operator_radial
use operator_3d
implicit none
integer :: idicp,jdicp,kdicp,ldicp
complex(kind=k2) :: phi_in(fedvr3d%nb_r,fedvr3d%nb_angle),phi_out(fedvr3d%nb_r,fedvr3d%nb_angle)
complex(kind=k2) :: temp1
do idicp =1,fedvr3d%nb_angle
do jdicp =1,fedvr3d%nb_r
temp1 = zzero
do kdicp =1,fedvr3d%nb_r
temp1 = temp1 + h_stiffness(ia(max(jdicp,kdicp))+min(jdicp,kdicp),idicp)*phi_in(kdicp,idicp)
enddo
phi_out(jdicp,idicp) = temp1
enddo
enddo
!=================================================================================================
! the laser interaction term
!=================================================================================================
if(laser%tdornot) then
do idicp=1,fedvr3d%nb_angle
do jdicp =1,fedvr3d%nb_r
temp1 = zzero
do ldicp =1,fedvr3d%nb_angle
temp1 = temp1 + laser%ez_t*zmat_3d(idicp,ldicp)*fedvrx_global(jdicp)*phi_in(jdicp,ldicp)
enddo
phi_out(jdicp,idicp) = phi_out(jdicp,idicp) + temp1
enddo
enddo
endif
return
end subroutine act2_fc