Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_wann_mmkb_sph
8 : use m_juDFT
9 : contains
10 32 : subroutine wann_mmkb_sph(
11 : c***********************************************************************
12 : c computes the Mmn(k,b) matrix elements which are the overlaps
13 : c between the Bloch wavefunctions at the state m and n
14 : c at the k-point and its nearest neighbor b (bpt), in the spheres
15 : c Y.Mokrousov 15.6.06
16 : c***********************************************************************
17 : c the gaunt coefficients impose certain conditions (which can be
18 : c seen in the gaunt.F) on the l,lp,lpp angular momenta, which can be
19 : c used in order to speedup the procedure.
20 : c The conditions are the following:
21 : c 1. |l - l''| <= l' <= l + l''
22 : c 2. m' = m + m''
23 : c 3. l + l' + l'' - is even
24 : c***********************************************************************
25 : c modified for use in conjunction with wann_ujugaunt, which
26 : c calculates certain matrix elements that have to be treated only
27 : c once per calculation
28 : c Frank Freimuth, October 2006
29 : c***********************************************************************
30 : c use BLAS matrix-matrix multiplication in apw-apw part to be faster
31 : c J.-P. Hanke, Dec. 2015
32 : c***********************************************************************
33 : > nbnd,llod,nslibd,nslibd_b,nlod,natd,ntypd,lmd,jmtd,
34 128 : > taual,nop,lmax,
35 96 : > ntype,neq,nlo,llo,acof,bcof,ccof,bbpt,
36 160 : > acof_b,bcof_b,ccof_b,gb,bkpt,ujug,ujdg,djug,djdg,ujulog,
37 96 : > djulog,ulojug,ulojdg,ulojulog,kdiff,nntot,chi,
38 32 : = mmn)
39 :
40 : use m_juDFT
41 : use m_constants, only : pimach
42 : use m_sphbes
43 : use m_ylm
44 : use m_intgr, only : intgr3
45 : use m_gaunt, only: gaunt1
46 :
47 : implicit none
48 :
49 : c .. scalar arguments ..
50 : integer, intent (in) :: llod,nlod,natd,ntypd,lmd,nbnd
51 : integer, intent (in) :: nntot
52 : integer, intent (in) :: ntype,nslibd,nslibd_b,nop,jmtd
53 :
54 : c .. array arguments ..
55 : integer, intent (in) :: neq(:) !(ntypd)
56 : integer, intent (in) :: lmax(:) !(ntypd)
57 : integer, intent (in) :: nlo(:) !(ntypd)
58 : integer, intent (in) :: llo(:,:) !(nlod,ntypd)
59 : complex, intent (in) :: chi(ntypd)
60 : real, intent (in) :: bbpt(:) !(3)
61 : real, intent (in) :: taual(:,:) !(3,natd)
62 : real, intent (in) :: bkpt(:) !(3)
63 : integer, intent (in) :: gb(:) !(3)
64 : complex, intent (in) :: ccof(-llod:,:,:,:) !(-llod:llod,nslibd,nlod,natd)
65 : complex, intent (in) :: acof(:,0:,:) !(nslibd,0:lmd,natd)
66 : complex, intent (in) :: bcof(:,0:,:) !(nslibd,0:lmd,natd)
67 : complex, intent (in) :: ccof_b(-llod:,:,:,:) !(-llod:llod,nslibd_b,nlod,natd)
68 : complex, intent (in) :: acof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
69 : complex, intent (in) :: bcof_b(:,0:,:) !(nslibd_b,0:lmd,natd)
70 :
71 : complex, intent (in) :: ujug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
72 : complex, intent (in) :: ujdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
73 : complex, intent (in) :: djug(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
74 : complex, intent (in) :: djdg(0:,0:,:,:) !(0:lmd,0:lmd,1:ntype,1:nntot)
75 : complex, intent (in) :: ujulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
76 : complex, intent (in) :: djulog(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
77 : complex, intent (in) :: ulojug(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
78 : complex, intent (in) :: ulojdg(0:,:,-llod:,:,:) !(0:lmd,nlod,-llod:llod,1:ntype,1:nntot)
79 : complex, intent (in) :: ulojulog(:,-llod:,:,-llod:,:,:) !(1:nlod,-llod:llod,1:nlod,-llod:llod,1:ntype,1:nntot)
80 :
81 : real, intent (in) :: kdiff(:,:) !(3,nntot)
82 : complex,intent(inout) :: mmn(:,:)!nbnd,nbnd) !(nbnd,nbnd)
83 :
84 :
85 : c .. local scalars ..
86 : integer i,lm,nn,n,na,j,lmp,l,lp,m,mp,lwn,lo,lop
87 : integer ll,llp
88 : real rph,cph,tpi,th,t1nn,t2nn,t3nn,dummy
89 : complex ic
90 : integer nene,lmd2
91 : complex :: fac1,fac2,fac3,fac4
92 : complex :: fac5,fac6,fac7,fac8
93 32 : complex :: mat(0:lmd,nslibd_b)
94 : C ..
95 : C .. local arrays ..
96 : real bpt(3)
97 :
98 : C ..
99 : C .. intrinsic functions ..
100 : intrinsic conjg,cmplx,sqrt,cos,sin
101 :
102 32 : call timestart("wann_mmkb_sph")
103 32 : ic = cmplx(0.,1.)
104 32 : tpi = 2* pimach()
105 32 : lmd2 = lmd+1
106 :
107 32 : na = 0
108 :
109 128 : bpt(:) = bbpt(:) + gb(:) - bkpt(:)
110 144 : do nene=1,nntot
111 276 : if(all(abs(bpt(:)-kdiff(:,nene)).lt.1e-4)) exit
112 : enddo
113 32 : IF(nene==nntot+1) CALL juDFT_error
114 : + ("cannot find matching nearest neighbor k",calledby
115 0 : + ="wann_mmkb_sph")
116 :
117 64 : do n=1,ntype
118 32 : lwn = lmax(n)
119 128 : do nn = 1,neq(n) ! cycle by the atoms within the atom type
120 64 : na = na + 1
121 : c...set up phase factors ( e^{ib\tau} )
122 :
123 64 : t1nn = tpi*taual(1,na)
124 64 : t2nn = tpi*taual(2,na)
125 64 : t3nn = tpi*taual(3,na)
126 :
127 64 : th = bpt(1)*t1nn + bpt(2)*t2nn + bpt(3)*t3nn
128 64 : rph = 2*tpi*cos(th)
129 64 : cph = 2*tpi*sin(th)
130 :
131 : c...contributions from the apws
132 : call zgemm(
133 : > 'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
134 : > ujug(0,0,n,nene),lmd2,acof_b(1,0,na),nslibd_b,
135 64 : > cmplx(0.0),mat(0,1),lmd2)
136 : call zgemm(
137 : > 'N','N',nslibd,nslibd_b,lmd2,chi(n),
138 : > acof(1,0,na),nslibd,mat(0,1),lmd2,
139 64 : > cmplx(1.0),mmn,nbnd)
140 :
141 : call zgemm(
142 : > 'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
143 : > ujdg(0,0,n,nene),lmd2,bcof_b(1,0,na),nslibd_b,
144 64 : > cmplx(0.0),mat(0,1),lmd2)
145 : call zgemm(
146 : > 'N','N',nslibd,nslibd_b,lmd2,chi(n),
147 : > acof(1,0,na),nslibd,mat(0,1),lmd2,
148 64 : > cmplx(1.0),mmn,nbnd)
149 :
150 : call zgemm(
151 : > 'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
152 : > djug(0,0,n,nene),lmd2,acof_b(1,0,na),nslibd_b,
153 64 : > cmplx(0.0),mat(0,1),lmd2)
154 : call zgemm(
155 : > 'N','N',nslibd,nslibd_b,lmd2,chi(n),
156 : > bcof(1,0,na),nslibd,mat(0,1),lmd2,
157 64 : > cmplx(1.0),mmn,nbnd)
158 :
159 : call zgemm(
160 : > 'T','C',lmd2,nslibd_b,lmd2,cmplx(rph,cph),
161 : > djdg(0,0,n,nene),lmd2,bcof_b(1,0,na),nslibd_b,
162 64 : > cmplx(0.0),mat(0,1),lmd2)
163 : call zgemm(
164 : > 'N','N',nslibd,nslibd_b,lmd2,chi(n),
165 : > bcof(1,0,na),nslibd,mat(0,1),lmd2,
166 64 : > cmplx(1.0),mmn,nbnd)
167 :
168 :
169 :
170 96 : if (nlo(n).ge.1) then
171 :
172 0 : do lo = 1,nlo(n)
173 0 : l = llo(lo,n)
174 0 : ll = l*(l+1)
175 0 : do m = -l,l
176 : lm = ll + m
177 :
178 0 : do lp = 0,lwn
179 0 : llp = lp*(lp+1)
180 0 : do mp = -lp,lp
181 0 : lmp = llp + mp
182 :
183 0 : fac1=cmplx(rph,cph)*ujulog(lmp,lo,m,n,nene)*chi(n)
184 0 : fac2=cmplx(rph,cph)*djulog(lmp,lo,m,n,nene)*chi(n)
185 0 : fac3=cmplx(rph,cph)*ulojdg(lmp,lo,m,n,nene)*chi(n)
186 0 : fac4=cmplx(rph,cph)*ulojug(lmp,lo,m,n,nene)*chi(n)
187 :
188 : c if(fac1.ne.0)write(*,*)'fac1',lmp,lm
189 : c if(fac2.ne.0)write(*,*)'fac2',lmp,lm
190 : c if(fac3.ne.0)write(*,*)'fac3',lmp,lm
191 : c if(fac4.ne.0)write(*,*)'fac4',lmp,lm
192 :
193 0 : do j = 1,nslibd_b
194 0 : do i = 1,nslibd
195 : mmn(i,j) = mmn(i,j) +
196 : + ccof(m,i,lo,na)*
197 : * (conjg(acof_b(j,lmp,na))*fac1+
198 : + conjg(bcof_b(j,lmp,na))*fac2 )+
199 :
200 : + conjg(ccof_b(m,j,lo,na))*
201 : * ( bcof(i,lmp,na)*fac3 +
202 0 : + acof(i,lmp,na)*fac4 )
203 : enddo
204 : enddo
205 :
206 :
207 : enddo ! mp
208 : enddo ! lp
209 :
210 : enddo ! lo
211 : enddo ! m lo
212 :
213 : c...contributions from lo*lo
214 :
215 0 : do lo = 1,nlo(n)
216 0 : l = llo(lo,n)
217 0 : ll = l*(l+1)
218 0 : do m = -l,l
219 : lm = ll + m
220 :
221 0 : do lop = 1,nlo(n)
222 0 : lp = llo(lop,n)
223 0 : llp = lp*(lp+1)
224 0 : do mp = -lp,lp
225 0 : lmp = llp + mp
226 :
227 0 : fac1=cmplx(rph,cph)*ulojulog(lop,mp,lo,m,n,nene)*chi(n)
228 :
229 : c if(fac1.ne.0)write(*,*)'fac1',lop,lo,mp,m
230 :
231 0 : do j = 1,nslibd_b
232 0 : do i = 1,nslibd
233 :
234 : mmn(i,j) = mmn(i,j) +
235 : + ccof(m,i,lo,na)*conjg(ccof_b(mp,j,lop,na))*
236 0 : * fac1
237 :
238 : enddo
239 : enddo
240 :
241 : enddo ! mp lop
242 : enddo ! lop
243 :
244 : enddo ! m lo
245 : enddo ! lo
246 :
247 : endif ! local orbitals on this atom
248 :
249 : enddo ! atoms in the type
250 :
251 : enddo ! atom type
252 :
253 :
254 32 : call timestop("wann_mmkb_sph")
255 32 : end subroutine wann_mmkb_sph
256 : end module m_wann_mmkb_sph
|
