-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathorb_eqs_p_zet.f90
84 lines (77 loc) · 2.39 KB
/
orb_eqs_p_zet.f90
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
!
! same with haru
!
!==================================================================
subroutine orb_eqs_p_zet
!==================================================================
use global
implicit none
call orb_eqs_p_zet2(system%NP1,system%NP2,system%NP0,system%NP0+system%NP1)
return
end subroutine orb_eqs_p_zet
!==================================================================
subroutine orb_eqs_p_zet2(NNP1,NNP2,NNA1,NNA2)
!==================================================================
use global
use auxiliary
use operator_spatial_space
implicit none
integer(4) :: NNP1,NNP2,NNP12,NNA1,NNA2
integer(4) :: iijj,kkll,i,j,ii,jjp,k,l
integer(4) :: info,ipiv2(NNP1*NNP2)
complex(kind=k2) :: cs,ce2(NNP1*NNP2),ca2(NNP1*NNP2,NNP1*NNP2)
complex(kind=k2) :: ca2_reg(NNP1*NNP2,NNP1*NNP2)
complex(kind=k2) :: here_temp,here_temp1,here_temp2
NNP12=NNP1*NNP2
!============================================================================
! calculation of ce1 and ce1
!============================================================================
iijj=0
do j=1,NNP2
do i=1,NNP1
iijj=iijj+1
kkll=0
do k=1,NNP2
do l=1,NNP1
kkll=kkll+1
ca2(iijj,kkll)=czet1(iijj,kkll)
enddo
enddo
enddo
enddo
!==============================================================================
!==============================================================================
!---- calculataion of ce1 -------------------
iijj=0
do j=1,NNP2
do i=1,NNP1
iijj=iijj+1
ce2(iijj)=cv_zet2(iijj)
enddo
enddo
call reg_mat(ca2,ca2_reg,NNP12,eps_Q)
if (NNP2.gt.0) call zgesv(NNP12,1,ca2_reg,NNP12,ipiv2,ce2,NNP12,info)
!----------------
iijj=0
do j=1,NNP2
jjp=j+NNA2
do i=1,NNP1
ii=i+NNA1
iijj=iijj+1
cs=-ce2(iijj)
ch_dummy(jjp,ii)=cs
! CE(ii,jjp)=dconjg(cs)
here_temp = ham_onebody(ia( max(ii,jjp) ) +min(ii,jjp) )
if(jjp >= ii) then
here_temp1 = here_temp
here_temp2 = dconjg(here_temp)
else
here_temp1 = dconjg(here_temp)
here_temp2 = here_temp
endif
! CE(ii,jjp)=dconjg(cs-CH(jjp,ii))+CH(ii,jjp) ! <-- Be careful.
ch_dummy(ii,jjp)=dconjg(cs-here_temp1)+ here_temp2 !! ch_dummy is h_i^j-i eta_i^j
enddo
enddo
return
end subroutine orb_eqs_p_zet2