Line data Source code
1 : c*****************************************c
2 : c Interstitial contribution to uHu c
3 : c < u_{k+b1} | H_{k}^{int} | u_{k+b2} > c
4 : c*****************************************c
5 : c k1_b, k2_b, k3_b : G-vectors at k+b1 c
6 : c k1_b2, k2_b2, k3_b2: G-vectors at k+b2 c
7 : c c
8 : c gb : reciprocal latt. vector taking c
9 : c k+b1 back to first Brillouin zone c
10 : c gb2: reciprocal latt. vector taking c
11 : c k+b2 back to first Brillouin zone c
12 : c c
13 : c z_b : eigenvectors at k+b1 c
14 : c z_b2: eigenvectors at k+b2 c
15 : c*****************************************c
16 : c J.-P. Hanke, Dec. 2015 c
17 : c*****************************************c
18 : module m_wann_uHu_int
19 : use m_juDFT
20 : contains
21 0 : subroutine wann_uHu_int(
22 : > chi,nvd,k1d,k2d,k3d,
23 : > n3d,nv_b,nv_b2,nbnd,neigd,
24 : > nslibd_b,nslibd_b2,nbasfcn_b,nbasfcn_b2,
25 : > addnoco,addnoco2,
26 0 : > k1_b ,k2_b ,k3_b, gb,
27 0 : > k1_b2,k2_b2,k3_b2,gb2,
28 0 : > bkpt,bbmat,vpw,zMat_b,zMat_b2,
29 0 : > rgphs,ustep,ig,l_kin,sign,uHu)
30 :
31 : USE m_types
32 :
33 : implicit none
34 :
35 : TYPE(t_mat), INTENT(IN) :: zMat_b, zMat_b2
36 :
37 : c ..arguments..
38 : logical, intent(in) :: l_kin
39 : integer, intent(in) :: gb(3),gb2(3)
40 : integer, intent(in) :: k1d,k2d,k3d
41 : integer, intent(in) :: nvd,n3d,nv_b,nv_b2,nbnd,neigd
42 : integer, intent(in) :: nslibd_b,nslibd_b2,nbasfcn_b,nbasfcn_b2
43 : integer, intent(in) :: addnoco,addnoco2,sign
44 : integer, intent(in) :: k1_b(nvd) ,k2_b(nvd) ,k3_b(nvd)
45 : integer, intent(in) :: k1_b2(nvd),k2_b2(nvd),k3_b2(nvd)
46 : integer, intent(in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
47 : complex, intent(in) :: rgphs(-k1d:k1d,-k2d:k2d,-k3d:k3d)
48 : real, intent(in) :: bkpt(3),bbmat(3,3)
49 : complex, intent(in) :: ustep(n3d),vpw(n3d)
50 : complex, intent(in) :: chi
51 :
52 : complex, intent(inout) :: uHu(nbnd,nbnd)
53 :
54 : c ..local variables..
55 0 : complex,allocatable :: vstep_c(:,:)
56 0 : complex,allocatable :: mat_c(:,:)
57 : complex :: th_c
58 0 : real,allocatable :: vstep_r(:,:)
59 0 : real,allocatable :: mat_r(:,:)
60 0 : real,allocatable :: uHu_tmp(:,:)
61 : real :: th_r
62 :
63 0 : real,allocatable :: rk_b(:),rk_b2(:)
64 : real :: phase,phase2,ekk,s(3)
65 : integer :: i,j,j1,j2,j3,i1,i2,i3,in,ind
66 :
67 0 : allocate( rk_b (nv_b) )
68 0 : allocate( rk_b2(nv_b2) )
69 :
70 0 : IF(zMat_b%l_real) THEN
71 0 : allocate( uHu_tmp(nslibd_b,nslibd_b2) )
72 0 : allocate( vstep_r(nv_b2,nv_b) )
73 0 : allocate( mat_r(nv_b,nslibd_b2) )
74 0 : vstep_r(:,:) = 0.0
75 : ELSE
76 0 : allocate( vstep_c(nv_b2,nv_b) )
77 0 : allocate( mat_c(nv_b,nslibd_b2) )
78 0 : vstep_c(:,:) = CMPLX(0.0,0.0)
79 : END IF
80 :
81 : ! set up |k+G-G(k+b1)|^2
82 0 : do i=1,nv_b
83 0 : s(1) = bkpt(1) + k1_b(i) - gb(1)
84 0 : s(2) = bkpt(2) + k2_b(i) - gb(2)
85 0 : s(3) = bkpt(3) + k3_b(i) - gb(3)
86 0 : rk_b(i) = dot_product(s,matmul(bbmat,s))
87 : ! rk_b(i) = dotirp(s,s,bbmat)
88 : enddo
89 :
90 : ! set up |k+G'-G(k+b2)|^2
91 0 : do i=1,nv_b2
92 0 : s(1) = bkpt(1) + k1_b2(i) - gb2(1)
93 0 : s(2) = bkpt(2) + k2_b2(i) - gb2(2)
94 0 : s(3) = bkpt(3) + k3_b2(i) - gb2(3)
95 0 : rk_b2(i) = dot_product(s,matmul(bbmat,s))
96 : ! rk_b2(i) = dotirp(s,s,bbmat)
97 : enddo
98 :
99 : ! construct vstep(g,g') ~ V(g-g')
100 : ! + Theta(g-g')*[rk_b+rk_b2]
101 0 : do i=1,nv_b
102 0 : j1 = -k1_b(i) + gb(1) - gb2(1)
103 0 : j2 = -k2_b(i) + gb(2) - gb2(2)
104 0 : j3 = -k3_b(i) + gb(3) - gb2(3)
105 0 : do j=1,nv_b2
106 0 : i1 = j1 + k1_b2(j)
107 0 : i2 = j2 + k2_b2(j)
108 0 : i3 = j3 + k3_b2(j)
109 0 : in = ig(sign*i1,sign*i2,sign*i3)
110 0 : if(in.eq.0) cycle
111 0 : phase = rgphs(i1,i2,i3) ! TODO: sign also here?
112 0 : phase2= rgphs(sign*i1,sign*i2,sign*i3)
113 0 : if(phase.ne.phase2) then
114 : call juDFT_error("rgphs in wann_uHu_int",
115 0 : & calledby="wann_uHu_int")
116 : endif
117 :
118 0 : ekk = rk_b(i) + rk_b2(j)
119 :
120 0 : IF(zMat_b%l_real) THEN
121 0 : th_r = phase*real(vpw(in))
122 0 : if(l_kin) th_r = th_r + phase*0.25*ekk*real(ustep(in))
123 0 : vstep_r(j,i) = th_r
124 : ELSE
125 0 : th_c = phase*conjg(vpw(in))
126 0 : if(l_kin) th_c = th_c + phase*0.25*ekk*conjg(ustep(in))
127 0 : vstep_c(j,i) = th_c
128 : END IF
129 : enddo
130 : enddo
131 :
132 : ! complex conjugate of (z(k+b1,g))^* vstep(g,g') z(k+b2,g')
133 0 : IF(zMat_b%l_real) THEN
134 : call dgemm('T','N',nv_b,nslibd_b2,nv_b2,real(1.0),
135 : > vstep_r,nv_b2,zMat_b2%data_r(1+addnoco2,1),nbasfcn_b2,
136 0 : > real(0.0),mat_r,nv_b)
137 : call dgemm('T','N',nslibd_b,nslibd_b2,nv_b,
138 : > real(1.0),zMat_b%data_r(1+addnoco,1),nbasfcn_b,
139 0 : > mat_r,nv_b,real(0.0),uHu_tmp,nslibd_b)
140 : uHu(1:nslibd_b,1:nslibd_b2) = uHu(1:nslibd_b,1:nslibd_b2)
141 0 : > + uHu_tmp(1:nslibd_b,1:nslibd_b2)*chi
142 : ELSE
143 : call zgemm('T','N',nv_b,nslibd_b2,nv_b2,cmplx(1.0),
144 : > vstep_c,nv_b2,zMat_b2%data_c(1+addnoco2,1),nbasfcn_b2,
145 0 : > cmplx(0.0),mat_c,nv_b)
146 0 : mat_c = conjg(mat_c)
147 : call zgemm('T','N',nslibd_b,nslibd_b2,nv_b,
148 : > chi,zMat_b%data_c(1+addnoco,1),nbasfcn_b,
149 0 : > mat_c,nv_b,cmplx(1.0),uHu,nbnd)
150 : END IF
151 :
152 0 : deallocate( rk_b, rk_b2 )
153 :
154 0 : IF(zMat_b%l_real) THEN
155 0 : deallocate( uHu_tmp )
156 0 : deallocate( vstep_r, mat_r )
157 : ELSE
158 0 : deallocate( vstep_c, mat_c )
159 : END IF
160 :
161 0 : end subroutine
162 : end module m_wann_uHu_int
|
