Line data Source code
1 : c****************************************c
2 : c Muffin tin contribution to uHu c
3 : c < u_{k+b1} | H_{k}^{mt} | u_{k+b2} > c
4 : c****************************************c
5 : c Use T(lmp,lm) calculated before for c
6 : c every pair (b1,b2) and all atoms n c
7 : c c
8 : c acof ,bcof ,ccof : coefficients c
9 : c at k+b1 c
10 : c acof_b,bcof_b,ccof_b: coefficients c
11 : c at k+b2 c
12 : c c
13 : c bkpt : k-point c
14 : c bkpt_b : (k+b1)-point c
15 : c bkpt_b2: (k+b2)-point c
16 : c****************************************c
17 : c J.-P. Hanke, Dec. 2015 c
18 : c****************************************c
19 : MODULE m_wann_uHu_sph
20 : contains
21 0 : subroutine wann_uHu_sph(
22 : > chi,nbnd,llod,nslibd,nslibd_b,nlod,natd,ntypd,lmd,jmtd,
23 0 : > taual,nop,lmax,
24 0 : > ntype,neq,nlo,llo,acof,bcof,ccof,bkpt_b2,
25 0 : > acof_b,bcof_b,ccof_b,bkpt_b,bkpt,gb_b,gb_b2,
26 0 : > tuu,tud,tdu,tdd,tuulo,tulou,tdulo,tulod,tuloulo,
27 0 : > kdiff,kdiff2,nntot,nntot2,uHu)
28 :
29 : USE m_juDFT
30 : use m_constants, only : pimach
31 :
32 : implicit none
33 :
34 : c .. scalar arguments ..
35 : integer, intent (in) :: llod,nlod,natd,ntypd,lmd,nbnd
36 : integer, intent (in) :: nntot,nntot2
37 : integer, intent (in) :: ntype,nslibd,nslibd_b,nop,jmtd
38 : complex, intent (in) :: chi
39 :
40 : c .. array arguments ..
41 : integer, intent (in) :: neq(:) !(ntypd)
42 : integer, intent (in) :: lmax(:) !(ntypd)
43 : integer, intent (in) :: nlo(:) !(ntypd)
44 : integer, intent (in) :: llo(:,:) !(nlod,ntypd)
45 : real, intent (in) :: bkpt_b(:) !(3)
46 : real, intent (in) :: bkpt_b2(:) !(3)
47 : real, intent (in) :: taual(:,:) !(3,natd)
48 : real, intent (in) :: bkpt(:) !(3)
49 : integer, intent (in) :: gb_b(:) !(3)
50 : integer, intent (in) :: gb_b2(:) !(3)
51 : complex, intent (in) :: ccof(-llod:,:,:,:) !(-llod:llod,nslibd,nlod,natd)
52 : complex, intent (in) :: acof(:,0:,:) !(nslibd,0:lmd,natd)
53 : complex, intent (in) :: bcof(:,0:,:) !(nslibd,0:lmd,natd)
54 : complex, intent (in) :: ccof_b(-llod:,:,:,:) !(-llod:llod,nslibd_b,nlod,natd)
55 : complex, intent (in) :: acof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
56 : complex, intent (in) :: bcof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
57 :
58 : complex, intent (in) :: tuu(0:,0:,:,:)
59 : complex, intent (in) :: tud(0:,0:,:,:)
60 : complex, intent (in) :: tdu(0:,0:,:,:)
61 : complex, intent (in) :: tdd(0:,0:,:,:)
62 : complex, intent (in) :: tuulo(0:,:,-llod:,:,:)
63 : complex, intent (in) :: tulou(0:,:,-llod:,:,:)
64 : complex, intent (in) :: tdulo(0:,:,-llod:,:,:)
65 : complex, intent (in) :: tulod(0:,:,-llod:,:,:)
66 : complex, intent (in) :: tuloulo(:,-llod:,:,-llod:,:,:)
67 :
68 : real, intent (in) :: kdiff(:,:),kdiff2(:,:)
69 : complex,intent(inout) :: uHu(:,:)
70 :
71 :
72 : c .. local scalars ..
73 : integer i,lm,nn,n,na,j,lmp,l,lp,m,mp,lwn,lo,lop
74 : integer ll,llp,lmd2
75 : real rph,cph,tpi,sqpi16,th,t1nn,t2nn,t3nn
76 : integer nene,nene2,indexx
77 : complex :: fac1,fac2,fac3,fac4
78 0 : complex :: mat(0:lmd,nslibd_b)
79 : C ..
80 : C .. local arrays ..
81 : real bpt(3)
82 : real bpt2(3)
83 :
84 : C ..
85 : C .. intrinsic functions ..
86 : intrinsic conjg,cmplx,sqrt,cos,sin
87 :
88 0 : tpi = 2* pimach()
89 0 : sqpi16 = 4*tpi*tpi
90 0 : lmd2 = lmd+1
91 :
92 0 : na = 0
93 :
94 : ! find neighbor k+b1
95 0 : bpt(:) = bkpt_b(:) + gb_b(:) - bkpt(:)
96 0 : do nene=1,nntot
97 0 : if(all(abs(bpt(:)-kdiff(:,nene)).lt.1e-4)) exit
98 : enddo
99 0 : IF(nene==nntot+1) CALL juDFT_error
100 : + ("cannot find matching nearest neighbor k+b1",calledby
101 0 : + ="wann_uHu_sph")
102 :
103 : ! find neighbor k+b2
104 0 : bpt2(:) = bkpt_b2(:) + gb_b2(:) - bkpt(:)
105 0 : do nene2=1,nntot2
106 0 : if(all(abs(bpt2(:)-kdiff2(:,nene2)).lt.1e-4)) exit
107 : enddo
108 0 : IF(nene2==nntot2+1) CALL juDFT_error
109 : + ("cannot find matching nearest neighbor k+b2",calledby
110 0 : + ="wann_uHu_sph")
111 :
112 0 : indexx=nene2+(nene-1)*nntot2
113 :
114 : c if(nene2.ne.1) stop 'nene2.ne.1'
115 : c if(indexx.ne.nene) stop 'nene.ne.indexx'
116 : c if(ANY(bpt2.ne.0.0)) stop 'bpt2.ne.0'
117 :
118 0 : do n=1,ntype
119 0 : lwn = lmax(n)
120 0 : do nn = 1,neq(n) ! cycle by the atoms within the atom type
121 0 : na = na + 1
122 : c...set up phase factors ( 4pi e^{ib2\tau} 4pi e^{-ib1\tau} )
123 :
124 0 : t1nn = tpi*taual(1,na)
125 0 : t2nn = tpi*taual(2,na)
126 0 : t3nn = tpi*taual(3,na)
127 :
128 : th = (bpt2(1)-bpt(1))*t1nn
129 : > +(bpt2(2)-bpt(2))*t2nn
130 0 : > +(bpt2(3)-bpt(3))*t3nn
131 0 : rph = sqpi16*cos(th)
132 0 : cph = sqpi16*sin(th)
133 :
134 :
135 : c...apw-apw
136 : call zgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
137 : > tuu(0,0,n,indexx),lmd2,
138 : > acof_b(1,0,na),nslibd_b,
139 0 : > cmplx(0.0),mat(0,1),lmd2)
140 : call zgemm('N','N',nslibd,nslibd_b,lmd2,chi,
141 : > acof(1,0,na),nslibd,mat(0,1),lmd2,
142 0 : > cmplx(1.0),uHu,nbnd)
143 :
144 : call zgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
145 : > tud(0,0,n,indexx),lmd2,
146 : > bcof_b(1,0,na),nslibd_b,
147 0 : > cmplx(0.0),mat(0,1),lmd2)
148 : call zgemm('N','N',nslibd,nslibd_b,lmd2,chi,
149 : > acof(1,0,na),nslibd,mat(0,1),lmd2,
150 0 : > cmplx(1.0),uHu,nbnd)
151 :
152 : call zgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
153 : > tdu(0,0,n,indexx),lmd2,
154 : > acof_b(1,0,na),nslibd_b,
155 0 : > cmplx(0.0),mat(0,1),lmd2)
156 : call zgemm('N','N',nslibd,nslibd_b,lmd2,chi,
157 : > bcof(1,0,na),nslibd,mat(0,1),lmd2,
158 0 : > cmplx(1.0),uHu,nbnd)
159 :
160 : call zgemm('T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
161 : > tdd(0,0,n,indexx),lmd2,
162 : > bcof_b(1,0,na),nslibd_b,
163 0 : > cmplx(0.0),mat(0,1),lmd2)
164 : call zgemm('N','N',nslibd,nslibd_b,lmd2,chi,
165 : > bcof(1,0,na),nslibd,mat(0,1),lmd2,
166 0 : > cmplx(1.0),uHu,nbnd)
167 :
168 0 : if(nlo(n).ge.1) then
169 :
170 : c...apw-lo
171 0 : do lo = 1,nlo(n)
172 0 : l = llo(lo,n)
173 0 : do m = -l, l
174 :
175 0 : do lp = 0, lwn
176 0 : llp = lp*(lp+1)
177 0 : do mp = -lp, lp
178 0 : lmp = llp + mp
179 :
180 0 : fac1=cmplx(rph,cph)*tulou(lmp,lo,m,n,indexx)*chi
181 0 : fac2=cmplx(rph,cph)*tulod(lmp,lo,m,n,indexx)*chi
182 0 : fac3=cmplx(rph,cph)*tuulo(lmp,lo,m,n,indexx)*chi
183 0 : fac4=cmplx(rph,cph)*tdulo(lmp,lo,m,n,indexx)*chi
184 :
185 0 : do i = 1,nslibd
186 0 : do j = 1,nslibd_b
187 : uHu(i,j) = uHu(i,j)
188 : > + ccof(m,i,lo,na)* fac1 *conjg(acof_b(j,lmp,na))
189 : > + ccof(m,i,lo,na)* fac2 *conjg(bcof_b(j,lmp,na))
190 : > + acof(i,lmp,na) * fac3 *conjg(ccof_b(m,j,lo,na))
191 0 : > + bcof(i,lmp,na) * fac4 *conjg(ccof_b(m,j,lo,na))
192 : enddo
193 : enddo
194 :
195 : enddo !mp
196 : enddo !lp
197 :
198 : enddo ! m
199 : enddo ! lo
200 :
201 : c...lo-lo
202 0 : do lo = 1,nlo(n)
203 0 : l = llo(lo,n)
204 0 : do m = -l, l
205 :
206 0 : do lop = 1,nlo(n)
207 0 : lp = llo(lop,n)
208 0 : do mp = -lp, lp
209 :
210 0 : fac1=cmplx(rph,cph)*tuloulo(lop,mp,lo,m,n,indexx)*chi
211 :
212 0 : do i = 1,nslibd
213 0 : do j = 1,nslibd_b
214 : uHu(i,j) = uHu(i,j)
215 0 : > + ccof(m,i,lo,na)*fac1*conjg(ccof_b(mp,j,lop,na))
216 : enddo
217 : enddo
218 :
219 : enddo !mp
220 : enddo !lop
221 :
222 : enddo ! m
223 : enddo ! lo
224 :
225 : endif !(nlo(n).ge.1)
226 :
227 : enddo ! atoms in the type
228 :
229 : enddo ! atom type
230 :
231 0 : end subroutine wann_uHu_sph
232 : end module m_wann_uHu_sph
|
