Line data Source code
1 : MODULE m_kp_perturbation
2 : USE m_types_hybdat
3 :
4 : CONTAINS
5 :
6 0 : SUBROUTINE ibs_correction(nk, atoms, input, jsp, hybdat, mpdata, hybinp, &
7 : lapw, kpts, cell, mnobd, sym, noco,nococonv, &
8 0 : proj_ibsc, olap_ibsc)
9 :
10 : USE m_sphbes
11 : USE m_dsphbs
12 : USE m_constants
13 : USE m_ylm
14 : USE m_gaunt
15 : USE m_util
16 : use m_intgrf
17 : USE m_types
18 : USE m_io_hybrid
19 : IMPLICIT NONE
20 : TYPE(t_hybdat), INTENT(IN) :: hybdat
21 : TYPE(t_mpdata), intent(inout) :: mpdata
22 : TYPE(t_hybinp), INTENT(IN) :: hybinp
23 : TYPE(t_input), INTENT(IN) :: input
24 : TYPE(t_sym), INTENT(IN) :: sym
25 : TYPE(t_noco), INTENT(IN) :: noco
26 : type(t_nococonv), intent(in) :: nococonv
27 : TYPE(t_cell), INTENT(IN) :: cell
28 : TYPE(t_kpts), INTENT(IN) :: kpts
29 : TYPE(t_atoms), INTENT(IN) :: atoms
30 : TYPE(t_lapw), INTENT(IN) :: lapw
31 :
32 : ! - scalars -
33 : INTEGER, INTENT(IN) :: jsp
34 : INTEGER, INTENT(IN) :: mnobd
35 : INTEGER, INTENT(IN) :: nk
36 :
37 : ! - arrays -
38 :
39 : COMPLEX, INTENT(INOUT):: olap_ibsc(:, :, :, :)
40 : COMPLEX, INTENT(INOUT):: proj_ibsc(:, :, :)!(3,mnobd,hybdat%nbands(nk,jsp))
41 : ! - local scalars -
42 : INTEGER :: i, itype, ieq, iatom, iatom1, iband, iband1
43 : INTEGER :: iband2, ilo, ibas, ic, ikpt, ikvec, invsfct
44 : INTEGER :: irecl_cmt, irecl_z
45 : INTEGER :: j, m
46 : INTEGER :: l1, m1, p1, l2, m2, p2, l, p, lm, &
47 : lmp, lmp1, lmp2, lm1, lm2
48 : INTEGER :: ok, ig
49 : INTEGER :: idum
50 : REAL :: const
51 : REAL :: ka, kb
52 : REAL :: kvecn
53 : REAL :: olap_udot, olap_uulo, olap_udotulo
54 : REAL :: rdum
55 : REAL :: ws
56 : COMPLEX :: phase
57 : COMPLEX :: cj, cdj
58 : COMPLEX :: denom, var_enum
59 : COMPLEX :: cdum, cdum1, cdum2
60 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
61 : ! - local arrays -
62 0 : INTEGER :: lmp_start(atoms%ntype)
63 0 : REAL :: alo(atoms%nlod, atoms%ntype), blo(atoms%nlod, atoms%ntype), &
64 0 : clo(atoms%nlod, atoms%ntype)
65 0 : REAL :: u1_lo(atoms%jmtd, atoms%nlod, atoms%ntype), &
66 0 : u2_lo(atoms%jmtd, atoms%nlod, atoms%ntype)
67 : REAL :: kvec(3), qvec(3)
68 0 : REAL :: sbes(0:atoms%lmaxd + 1), dsbes(0:atoms%lmaxd + 1)
69 0 : REAL :: bas1_tmp(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd + 1, atoms%ntype), &
70 0 : bas2_tmp(atoms%jmtd, maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd + 1, atoms%ntype)
71 0 : REAL :: bas1_MT_tmp(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd + 1, atoms%ntype), &
72 0 : drbas1_MT_tmp(maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd + 1, atoms%ntype)
73 0 : REAL :: ru1(atoms%jmtd, 3, mnobd), ru2(atoms%jmtd, 3, mnobd)
74 0 : REAL :: iu1(atoms%jmtd, 3, mnobd), iu2(atoms%jmtd, 3, mnobd)
75 0 : REAL :: rintegrand(atoms%jmtd), iintegrand(atoms%jmtd), &
76 0 : integrand(atoms%jmtd)
77 :
78 : COMPLEX :: f(atoms%jmtd, mnobd)
79 0 : COMPLEX :: carr(3), carr2(3, hybdat%nbands(nk,jsp))
80 0 : COMPLEX :: ylm((atoms%lmaxd + 2)**2)
81 0 : COMPLEX, ALLOCATABLE :: u1(:, :, :, :, :), u2(:, :, :, :, :)
82 0 : COMPLEX, ALLOCATABLE :: cmt_lo(:, :, :, :)
83 0 : COMPLEX, ALLOCATABLE :: cmt_apw(:, :, :)
84 0 : TYPE(t_mat) :: z
85 0 : REAL :: work_r(input%neig)
86 0 : COMPLEX :: work_c(input%neig)
87 :
88 : !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
89 :
90 0 : bas1_tmp(:, :, 0:atoms%lmaxd, :) = hybdat%bas1(:, :, 0:atoms%lmaxd, :)
91 0 : bas2_tmp(:, :, 0:atoms%lmaxd, :) = hybdat%bas2(:, :, 0:atoms%lmaxd, :)
92 :
93 0 : bas1_MT_tmp(:, 0:atoms%lmaxd, :) = hybdat%bas1_MT(:, 0:atoms%lmaxd, :)
94 0 : drbas1_MT_tmp(:, 0:atoms%lmaxd, :) = hybdat%drbas1_MT(:, 0:atoms%lmaxd, :)
95 :
96 0 : bas1_tmp(:, :, atoms%lmaxd + 1, :) = hybdat%bas1(:, :, atoms%lmaxd, :)
97 0 : bas2_tmp(:, :, atoms%lmaxd + 1, :) = hybdat%bas2(:, :, atoms%lmaxd, :)
98 :
99 0 : bas1_MT_tmp(:, atoms%lmaxd + 1, :) = hybdat%bas1_MT(:, atoms%lmaxd, :)
100 0 : drbas1_MT_tmp(:, atoms%lmaxd + 1, :) = hybdat%drbas1_MT(:, atoms%lmaxd, :)
101 :
102 : ! read in z coefficient from direct access file z at k-point nk
103 :
104 0 : call read_z(atoms, cell, hybdat, kpts, sym, noco, nococonv, input, nk, jsp, z)
105 :
106 : ! construct local orbital consisting of radial function times spherical harmonic
107 : ! where the radial function vanishes on the MT sphere boundary
108 : ! with this the local orbitals have a trivial k-dependence
109 :
110 : ! compute radial lo matching coefficients
111 0 : mpdata%num_radfun_per_l = 2
112 0 : DO itype = 1, atoms%ntype
113 0 : DO ilo = 1, atoms%nlo(itype)
114 0 : l = atoms%llo(ilo, itype)
115 0 : mpdata%num_radfun_per_l(l, itype) = mpdata%num_radfun_per_l(l, itype) + 1
116 0 : p = mpdata%num_radfun_per_l(l, itype)
117 :
118 0 : ws = -wronskian(hybdat%bas1_MT(1, l, itype), hybdat%drbas1_MT(1, l, itype), hybdat%bas1_MT(2, l, itype), hybdat%drbas1_MT(2, l, itype))
119 :
120 0 : ka = 1.0/ws*wronskian(hybdat%bas1_MT(p, l, itype), hybdat%drbas1_MT(p, l, itype), hybdat%bas1_MT(2, l, itype), hybdat%drbas1_MT(2, l, itype))
121 :
122 0 : kb = 1.0/ws*wronskian(hybdat%bas1_MT(1, l, itype), hybdat%drbas1_MT(1, l, itype), hybdat%bas1_MT(p, l, itype), hybdat%drbas1_MT(p, l, itype))
123 :
124 0 : integrand = hybdat%bas1(:, 2, l, itype)*hybdat%bas1(:, 2, l, itype) + hybdat%bas2(:, 2, l, itype)*hybdat%bas2(:, 2, l, itype)
125 0 : olap_udot = intgrf(integrand, atoms, itype, hybdat%gridf)
126 :
127 0 : integrand = hybdat%bas1(:, 1, l, itype)*hybdat%bas1(:, p, l, itype) + hybdat%bas2(:, 1, l, itype)*hybdat%bas2(:, p, l, itype)
128 0 : olap_uulo = intgrf(integrand, atoms, itype, hybdat%gridf)
129 :
130 0 : integrand = hybdat%bas1(:, 2, l, itype)*hybdat%bas1(:, p, l, itype) + hybdat%bas2(:, 2, l, itype)*hybdat%bas2(:, p, l, itype)
131 0 : olap_udotulo = intgrf(integrand, atoms, itype, hybdat%gridf)
132 :
133 0 : rdum = ka**2 + (kb**2)*olap_udot + 1.0 + 2.0*ka*olap_uulo + 2.0*kb*olap_udotulo
134 0 : clo(ilo, itype) = 1.0/sqrt(rdum)
135 0 : alo(ilo, itype) = ka*clo(ilo, itype)
136 0 : blo(ilo, itype) = kb*clo(ilo, itype)
137 :
138 0 : u1_lo(:, ilo, itype) = alo(ilo, itype)*hybdat%bas1(:, 1, l, itype) + blo(ilo, itype)*hybdat%bas1(:, 2, l, itype) + clo(ilo, itype)*hybdat%bas1(:, p, l, itype)
139 :
140 0 : u2_lo(:, ilo, itype) = alo(ilo, itype)*hybdat%bas2(:, 1, l, itype) + blo(ilo, itype)*hybdat%bas2(:, 2, l, itype) + clo(ilo, itype)*hybdat%bas2(:, p, l, itype)
141 : END DO
142 : END DO
143 :
144 : ! calculate lo wavefunction coefficients
145 0 : allocate(cmt_lo(input%neig, -atoms%llod:atoms%llod, atoms%nlod, atoms%nat))
146 0 : cmt_lo = 0
147 0 : iatom = 0
148 0 : ic = 0
149 0 : ibas = lapw%nv(jsp)
150 0 : DO itype = 1, atoms%ntype
151 : ! the program is in hartree units, therefore 1/wronskian is
152 : ! (rmt**2)/2.
153 0 : const = fpi_const*(atoms%rmt(itype)**2)/2/sqrt(cell%omtil)
154 0 : DO ieq = 1, atoms%neq(itype)
155 0 : iatom = iatom + 1
156 0 : IF((sym%invsat(iatom) == 0) .or. (sym%invsat(iatom) == 1)) THEN
157 0 : IF(sym%invsat(iatom) == 0) invsfct = 1
158 0 : IF(sym%invsat(iatom) == 1) THEN
159 0 : invsfct = 2
160 0 : iatom1 = sym%invsatnr(iatom)
161 : END IF
162 :
163 0 : DO ilo = 1, atoms%nlo(itype)
164 0 : l = atoms%llo(ilo, itype)
165 0 : cdum = img**l*const
166 0 : DO ikvec = 1, invsfct*(2*l + 1)
167 0 : ic = ic + 1
168 0 : ibas = ibas + 1
169 0 : call judft_error("this line has to be re-enabled, but kveclo_eig is never set")
170 : !kvec = kpts%bk(:, nk) + lapw%gvec(:, hybdat%kveclo_eig(ic, nk), jsp)
171 :
172 0 : phase = exp(img*tpi_const*dot_product(atoms%taual(:, iatom), kvec))
173 0 : cdum1 = cdum*phase
174 :
175 0 : CALL ylm4(l, matmul(kvec, cell%bmat), ylm)
176 :
177 0 : lm = l**2
178 0 : DO M = -l, l
179 0 : lm = lm + 1
180 0 : cdum2 = cdum1*conjg(ylm(lm))
181 0 : if(z%l_real) THEN
182 0 : work_r = z%data_r(ibas, :)
183 0 : DO iband = 1, hybdat%nbands(nk,jsp)
184 0 : cmt_lo(iband, M, ilo, iatom) = cmt_lo(iband, M, ilo, iatom) + cdum2*work_r(iband)
185 0 : IF(invsfct == 2) THEN
186 : ! the factor (-1)**l is necessary as we do not calculate
187 : ! the cmt_lo in the local coordinate system of the atom
188 0 : cmt_lo(iband, -M, ilo, iatom1) = cmt_lo(iband, -M, ilo, iatom1) + (-1)**(l + M)*conjg(cdum2)*work_r(iband)
189 : END IF
190 : END DO
191 : else
192 0 : work_c = z%data_c(ibas, :)
193 0 : DO iband = 1, hybdat%nbands(nk,jsp)
194 0 : cmt_lo(iband, M, ilo, iatom) = cmt_lo(iband, M, ilo, iatom) + cdum2*work_c(iband)
195 0 : IF(invsfct == 2) THEN
196 : ! the factor (-1)**l is necessary as we do not calculate
197 : ! the cmt_lo in the local coordinate system of the atom
198 0 : cmt_lo(iband, -M, ilo, iatom1) = cmt_lo(iband, -M, ilo, iatom1) + (-1)**(l + M)*conjg(cdum2)*work_c(iband)
199 : END IF
200 : END DO
201 : end if
202 :
203 : END DO
204 :
205 : END DO !ikvec
206 : END DO ! ilo
207 :
208 : END IF
209 :
210 : END DO !ieq
211 : END DO !itype
212 :
213 : !
214 : ! calculate apw wavefunction coefficients up to lmax + 1
215 : ! note that the lo contribution is separated in cmt_lo
216 : !
217 :
218 0 : DO itype = 1, atoms%ntype
219 0 : lmp_start(itype) = sum([(2*(2*l + 1), l=0, atoms%lmax(itype) + 1)])
220 : END DO
221 0 : idum = maxval(lmp_start)
222 :
223 0 : allocate(cmt_apw(input%neig, idum, atoms%nat))
224 0 : cmt_apw = 0
225 0 : DO i = 1, lapw%nv(jsp)
226 0 : kvec = kpts%bk(:, nk) + lapw%gvec(:, i, jsp)
227 0 : kvecn = sqrt(dot_product(matmul(kvec, cell%bmat), matmul(kvec, cell%bmat)))
228 :
229 0 : iatom = 0
230 0 : DO itype = 1, atoms%ntype
231 : !calculate spherical sperical harmonics
232 0 : CALL ylm4(atoms%lmax(itype) + 1, matmul(kvec, cell%bmat), ylm)
233 :
234 : !calculate spherical bessel function at |kvec|*R_MT(itype)
235 0 : CALL sphbes(atoms%lmax(itype) + 1, kvecn*atoms%rmt(itype), sbes)
236 :
237 : !calculate radial derivative of spherical bessel function at |kvec|*R_MT(itype)
238 0 : CALL dsphbs(atoms%lmax(itype) + 1, kvecn*atoms%rmt(itype), sbes, dsbes)
239 0 : dsbes = kvecn*dsbes
240 :
241 0 : DO ieq = 1, atoms%neq(itype)
242 0 : iatom = iatom + 1
243 :
244 0 : phase = exp(img*tpi_const*dot_product(kvec, atoms%taual(:, iatom)))
245 :
246 0 : lm = 0
247 0 : lmp = 0
248 0 : DO l = 0, atoms%lmax(itype) + 1
249 : denom = wronskian(bas1_MT_tmp(2, l, itype), drbas1_MT_tmp(2, l, itype), &
250 0 : bas1_MT_tmp(1, l, itype), drbas1_MT_tmp(1, l, itype))
251 0 : cdum1 = fpi_const*img**l*sbes(l)*phase/sqrt(cell%omtil)
252 0 : cdum2 = fpi_const*img**l*dsbes(l)*phase/sqrt(cell%omtil)
253 0 : DO M = -l, l
254 0 : lm = lm + 1
255 0 : cj = cdum1*conjg(ylm(lm))
256 0 : cdj = cdum2*conjg(ylm(lm))
257 0 : DO p = 1, 2
258 0 : lmp = lmp + 1
259 0 : p1 = p + (-1)**(p - 1)
260 :
261 : var_enum = CMPLX(wronskian(bas1_MT_tmp(p1, l, itype), drbas1_MT_tmp(p1, l, itype), REAL(cj), REAL(cdj)), &
262 0 : wronskian(bas1_MT_tmp(p1, l, itype), drbas1_MT_tmp(p1, l, itype), AIMAG(cj), AIMAG(cdj)))
263 :
264 0 : cdum = (-1)**(p + 1)*var_enum/denom
265 0 : if(z%l_real) THEN
266 0 : work_r = z%data_r(i, :)
267 0 : DO iband = 1, hybdat%nbands(nk,jsp)
268 0 : cmt_apw(iband, lmp, iatom) = cmt_apw(iband, lmp, iatom) + cdum*work_r(iband)
269 : END DO
270 : else
271 0 : work_c = z%data_c(i, :)
272 0 : DO iband = 1, hybdat%nbands(nk,jsp)
273 0 : cmt_apw(iband, lmp, iatom) = cmt_apw(iband, lmp, iatom) + cdum*work_c(iband)
274 : END DO
275 : end if
276 : END DO !p
277 : END DO !M
278 : END DO !l
279 :
280 : END DO !iatom
281 : END DO ! itype
282 : END DO ! i
283 :
284 : ! construct radial functions (complex) for the first order
285 : ! incomplete basis set correction
286 :
287 0 : allocate(u1(atoms%jmtd, 3, mnobd,(atoms%lmaxd + 1)**2, atoms%nat), stat=ok)!hybdat%nbands
288 0 : IF(ok /= 0) call judft_error('kp_perturbation: failure allocation u1')
289 0 : allocate(u2(atoms%jmtd, 3, mnobd,(atoms%lmaxd + 1)**2, atoms%nat), stat=ok)!hybdat%nbands
290 0 : IF(ok /= 0) call judft_error('kp_perturbation: failure allocation u2')
291 0 : u1 = 0; u2 = 0
292 :
293 0 : iatom = 0
294 0 : DO itype = 1, atoms%ntype
295 0 : DO ieq = 1, atoms%neq(itype)
296 0 : iatom = iatom + 1
297 :
298 0 : lm1 = 0
299 0 : DO l1 = 0, atoms%lmax(itype)! + 1
300 0 : DO m1 = -l1, l1
301 0 : lm1 = lm1 + 1
302 :
303 0 : DO p1 = 1, 2
304 :
305 0 : carr2 = 0
306 0 : l2 = l1 + 1
307 0 : lmp2 = 2*l2**2 + p1
308 0 : DO m2 = -l2, l2
309 0 : carr = gauntvec(l1, m1, l2, m2, atoms)
310 0 : DO iband = 1, mnobd! hybdat%nbands
311 0 : carr2(1:3, iband) = carr2(1:3, iband) + carr*cmt_apw(iband, lmp2, iatom)
312 : END DO
313 0 : lmp2 = lmp2 + 2
314 : END DO
315 :
316 0 : DO iband = 1, mnobd! hybdat%nbands
317 0 : DO i = 1, 3
318 0 : DO ig = 1, atoms%jri(itype)
319 : ! the r factor is already included in bas1
320 0 : u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) - img*bas1_tmp(ig, p1, l2, itype)*carr2(i, iband)
321 0 : u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) - img*bas2_tmp(ig, p1, l2, itype)*carr2(i, iband)
322 : END DO
323 : END DO
324 : END DO
325 :
326 0 : l2 = l1 - 1
327 0 : IF(l2 >= 0) THEN
328 0 : carr2 = 0
329 0 : lmp2 = 2*l2**2 + p1
330 0 : DO m2 = -l2, l2
331 0 : carr = gauntvec(l1, m1, l2, m2, atoms)
332 0 : DO iband = 1, mnobd! hybdat%nbands
333 0 : carr2(1:3, iband) = carr2(1:3, iband) + carr*cmt_apw(iband, lmp2, iatom)
334 : END DO
335 0 : lmp2 = lmp2 + 2
336 : END DO
337 :
338 0 : DO iband = 1, mnobd! hybdat%nbands
339 0 : DO i = 1, 3
340 0 : DO ig = 1, atoms%jri(itype)
341 : ! the r factor is already included in bas1
342 0 : u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) - img*bas1_tmp(ig, p1, l2, itype)*carr2(i, iband)
343 0 : u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) - img*bas2_tmp(ig, p1, l2, itype)*carr2(i, iband)
344 : END DO
345 : END DO
346 : END DO
347 :
348 : END IF
349 :
350 0 : carr2 = 0
351 0 : l2 = l1 + 1
352 0 : lmp2 = 2*l2**2
353 0 : DO m2 = -l2, l2
354 0 : carr = gauntvec(l1, m1, l2, m2, atoms)
355 0 : DO p2 = 1, 2
356 0 : lmp2 = lmp2 + 1
357 0 : rdum = w(p1, l1, p2, l2, itype, bas1_MT_tmp, drbas1_MT_tmp, atoms%rmt)
358 0 : DO iband = 1, mnobd! hybdat%nbands
359 0 : carr2(1:3, iband) = carr2(1:3, iband) + img*carr*rdum*cmt_apw(iband, lmp2, iatom)
360 : END DO
361 : END DO
362 : END DO
363 :
364 0 : DO iband = 1, mnobd! hybdat%nbands
365 0 : DO i = 1, 3
366 0 : DO ig = 1, atoms%jri(itype)
367 0 : u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) + bas1_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
368 0 : u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) + bas2_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
369 : END DO
370 : END DO
371 : END DO
372 :
373 0 : l2 = l1 - 1
374 0 : IF(l2 >= 0) THEN
375 0 : carr2 = 0
376 0 : lmp2 = 2*l2**2
377 0 : DO m2 = -l2, l2
378 0 : carr = gauntvec(l1, m1, l2, m2, atoms)
379 0 : DO p2 = 1, 2
380 0 : lmp2 = lmp2 + 1
381 0 : rdum = w(p1, l1, p2, l2, itype, bas1_MT_tmp, drbas1_MT_tmp, atoms%rmt)
382 0 : DO iband = 1, mnobd! hybdat%nbands
383 0 : carr2(1:3, iband) = carr2(1:3, iband) + img*carr*rdum*cmt_apw(iband, lmp2, iatom)
384 : END DO
385 : END DO
386 : END DO
387 :
388 0 : DO iband = 1, mnobd! hybdat%nbands
389 0 : DO i = 1, 3
390 0 : DO ig = 1, atoms%jri(itype)
391 : u1(ig, i, iband, lm1, iatom) = u1(ig, i, iband, lm1, iatom) &
392 0 : + bas1_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
393 : u2(ig, i, iband, lm1, iatom) = u2(ig, i, iband, lm1, iatom) &
394 0 : + bas2_tmp(ig, p1, l1, itype)*carr2(i, iband)/atoms%rmsh(ig, itype)
395 : END DO
396 : END DO
397 : END DO
398 : END IF
399 :
400 : END DO ! p1
401 : END DO !m1
402 : END DO !l1
403 : END DO !ieq
404 : END DO !iatom
405 :
406 : ! construct lo contribtution
407 0 : IF(any(atoms%llo == atoms%lmaxd)) call judft_error('ibs_correction: atoms%llo=atoms%lmaxd is not implemented')
408 :
409 0 : iatom = 0
410 0 : DO itype = 1, atoms%ntype
411 0 : DO ieq = 1, atoms%neq(itype)
412 0 : mpdata%num_radfun_per_l = 2
413 0 : iatom = iatom + 1
414 0 : DO ilo = 1, atoms%nlo(itype)
415 0 : l1 = atoms%llo(ilo, itype)
416 0 : mpdata%num_radfun_per_l(l1, itype) = mpdata%num_radfun_per_l(l1, itype) + 1
417 0 : p1 = mpdata%num_radfun_per_l(l1, itype)
418 :
419 0 : l2 = l1 + 1
420 0 : lm2 = l2**2
421 0 : DO m2 = -l2, l2
422 0 : lm2 = lm2 + 1
423 0 : carr2 = 0
424 :
425 0 : DO m1 = -l1, l1
426 0 : carr = gauntvec(l2, m2, l1, m1, atoms)
427 0 : DO iband = 1, mnobd
428 0 : carr2(1:3, iband) = carr2(1:3, iband) + cmt_lo(iband, m1, ilo, iatom)*carr
429 : END DO
430 : END DO
431 :
432 0 : DO iband = 1, mnobd
433 0 : DO i = 1, 3
434 0 : DO ig = 1, atoms%jri(itype)
435 : ! the r factor is already included in
436 0 : u1(ig, i, iband, lm2, iatom) = u1(ig, i, iband, lm2, iatom) - img*u1_lo(ig, ilo, itype)*carr2(i, iband)
437 0 : u2(ig, i, iband, lm2, iatom) = u2(ig, i, iband, lm2, iatom) - img*u2_lo(ig, ilo, itype)*carr2(i, iband)
438 : END DO
439 : END DO
440 : END DO
441 :
442 : END DO
443 :
444 0 : l2 = l1 - 1
445 0 : IF(l2 >= 0) THEN
446 0 : lm2 = l2**2
447 0 : DO m2 = -l2, l2
448 0 : lm2 = lm2 + 1
449 0 : carr2 = 0
450 :
451 0 : DO m1 = -l1, l1
452 0 : carr = gauntvec(l2, m2, l1, m1, atoms)
453 0 : DO iband = 1, mnobd
454 0 : carr2(1:3, iband) = carr2(1:3, iband) + cmt_lo(iband, m1, ilo, iatom)*carr
455 : END DO
456 : END DO
457 :
458 0 : DO iband = 1, mnobd
459 0 : DO i = 1, 3
460 0 : DO ig = 1, atoms%jri(itype)
461 : ! the r factor is already included in
462 0 : u1(ig, i, iband, lm2, iatom) = u1(ig, i, iband, lm2, iatom) - img*u1_lo(ig, ilo, itype)*carr2(i, iband)
463 0 : u2(ig, i, iband, lm2, iatom) = u2(ig, i, iband, lm2, iatom) - img*u2_lo(ig, ilo, itype)*carr2(i, iband)
464 : END DO
465 : END DO
466 : END DO
467 :
468 : END DO
469 : END IF
470 :
471 : END DO
472 : END DO
473 : END DO
474 :
475 : !
476 : ! calculate projection < phi(n',k)|phi^1(n,k)>
477 : ! n' = iband1 , n= iband2
478 : !
479 0 : iatom = 0
480 0 : proj_ibsc = 0
481 0 : DO itype = 1, atoms%ntype
482 0 : DO ieq = 1, atoms%neq(itype)
483 0 : iatom = iatom + 1
484 0 : lm = 0
485 0 : lmp = 0
486 0 : DO l = 0, atoms%lmax(itype)
487 0 : DO M = -l, l
488 0 : lm = lm + 1
489 0 : ru1 = real(u1(:, :, :, lm, iatom))
490 0 : iu1 = aimag(u1(:, :, :, lm, iatom))
491 0 : ru2 = real(u2(:, :, :, lm, iatom))
492 0 : iu2 = aimag(u2(:, :, :, lm, iatom))
493 0 : DO p = 1, 2
494 0 : lmp = lmp + 1
495 :
496 0 : DO iband = 1, mnobd! hybdat%nbands
497 0 : DO i = 1, 3
498 :
499 0 : rintegrand = atoms%rmsh(:, itype)*(hybdat%bas1(:, p, l, itype)*ru1(:, i, iband) + hybdat%bas2(:, p, l, itype)*ru2(:, i, iband))
500 :
501 0 : iintegrand = atoms%rmsh(:, itype)*(hybdat%bas1(:, p, l, itype)*iu1(:, i, iband) + hybdat%bas2(:, p, l, itype)*iu2(:, i, iband))
502 :
503 : carr2(i, iband) = intgrf(rintegrand, atoms, itype, hybdat%gridf) &
504 0 : + img*intgrf(iintegrand, atoms, itype, hybdat%gridf)
505 :
506 : END DO
507 : END DO
508 :
509 0 : DO iband1 = 1, hybdat%nbands(nk,jsp)
510 0 : cdum = conjg(cmt_apw(iband1, lmp, iatom))
511 0 : DO iband2 = 1, mnobd! hybdat%nbands
512 0 : proj_ibsc(1:3, iband2, iband1) = proj_ibsc(1:3, iband2, iband1) + cdum*carr2(1:3, iband2)
513 : END DO
514 : END DO
515 :
516 : END DO!p
517 : END DO!M
518 : END DO!l
519 :
520 : END DO!ieq
521 : END DO!itype
522 :
523 0 : iatom = 0
524 0 : DO itype = 1, atoms%ntype
525 0 : DO ieq = 1, atoms%neq(itype)
526 0 : iatom = iatom + 1
527 0 : DO ilo = 1, atoms%nlo(itype)
528 0 : l = atoms%llo(ilo, itype)
529 0 : lm = l**2
530 0 : DO M = -l, l
531 0 : lm = lm + 1
532 0 : ru1 = real(u1(:, :, :, lm, iatom))
533 0 : iu1 = aimag(u1(:, :, :, lm, iatom))
534 0 : ru2 = real(u2(:, :, :, lm, iatom))
535 0 : iu2 = aimag(u2(:, :, :, lm, iatom))
536 :
537 0 : DO iband = 1, mnobd! hybdat%nbands
538 0 : DO i = 1, 3
539 :
540 0 : rintegrand = atoms%rmsh(:, itype)*(u1_lo(:, ilo, itype)*ru1(:, i, iband) + u2_lo(:, ilo, itype)*ru2(:, i, iband))
541 :
542 0 : iintegrand = atoms%rmsh(:, itype)*(u1_lo(:, ilo, itype)*iu1(:, i, iband) + u2_lo(:, ilo, itype)*iu2(:, i, iband))
543 :
544 : carr2(i, iband) = intgrf(rintegrand, atoms, itype, hybdat%gridf) &
545 0 : + img*intgrf(iintegrand, atoms, itype, hybdat%gridf)
546 :
547 : END DO
548 : END DO
549 :
550 0 : DO iband1 = 1, hybdat%nbands(nk,jsp)
551 0 : cdum = conjg(cmt_lo(iband1, M, ilo, iatom))
552 0 : DO iband2 = 1, mnobd! hybdat%nbands
553 0 : proj_ibsc(1:3, iband2, iband1) = proj_ibsc(1:3, iband2, iband1) + cdum*carr2(1:3, iband2)
554 : END DO
555 : END DO
556 :
557 : END DO
558 :
559 : END DO
560 : END DO
561 : END DO
562 :
563 : !
564 : ! calculate <phi^1(n1,k)|phi^1(n2,k)>
565 : ! n1 and n2 occupied
566 : !
567 0 : iatom = 0
568 0 : olap_ibsc = 0
569 0 : DO itype = 1, atoms%ntype
570 0 : DO ieq = 1, atoms%neq(itype)
571 0 : iatom = iatom + 1
572 0 : lm = 0
573 0 : DO l = 0, atoms%lmax(itype)!+1
574 0 : DO M = -l, l
575 0 : lm = lm + 1
576 0 : ru1 = real(u1(:, :, :, lm, iatom))
577 0 : iu1 = aimag(u1(:, :, :, lm, iatom))
578 0 : ru2 = real(u2(:, :, :, lm, iatom))
579 0 : iu2 = aimag(u2(:, :, :, lm, iatom))
580 :
581 0 : DO iband1 = 1, mnobd ! hybdat%nbands
582 0 : DO iband2 = 1, mnobd!iband1
583 0 : DO i = 1, 3
584 0 : DO j = 1, 3
585 :
586 : rintegrand = atoms%rmsh(:, itype)**2*(ru1(:, i, iband1)*ru1(:, j, iband2) + ru2(:, i, iband1)*ru2(:, j, iband2) &
587 0 : + iu1(:, i, iband1)*iu1(:, j, iband2) + iu2(:, i, iband1)*iu2(:, j, iband2))
588 :
589 : iintegrand = atoms%rmsh(:, itype)**2*(ru1(:, i, iband1)*iu1(:, j, iband2) + ru2(:, i, iband1)*iu2(:, j, iband2) &
590 0 : - iu1(:, i, iband1)*ru1(:, j, iband2) - iu2(:, i, iband1)*ru2(:, j, iband2))
591 :
592 : olap_ibsc(i, j, iband2, iband1) = olap_ibsc(i, j, iband2, iband1) &
593 : + intgrf(rintegrand, atoms, itype, hybdat%gridf) &
594 0 : + img*intgrf(iintegrand, atoms, itype, hybdat%gridf)
595 :
596 : END DO
597 : END DO
598 :
599 : END DO
600 : END DO
601 :
602 : END DO
603 : END DO
604 : END DO
605 : END DO
606 :
607 0 : END SUBROUTINE ibs_correction
608 :
609 0 : FUNCTION gauntvec(l1, m1, l2, m2, atoms)
610 :
611 : USE m_constants
612 : USE m_gaunt
613 : USE m_juDFT
614 :
615 : USE m_types
616 : IMPLICIT NONE
617 : TYPE(t_atoms), INTENT(IN) :: atoms
618 :
619 : INTEGER, INTENT(IN) :: l1, m1, l2, m2
620 :
621 : COMPLEX :: gauntvec(-1:1)
622 :
623 : INTEGER :: j, mj
624 : INTEGER :: point(-1:1)
625 : REAL :: rfac
626 : COMPLEX :: cfac
627 : COMPLEX, PARAMETER :: img = (0.0, 1.0)
628 :
629 0 : rfac = sqrt(tpi_const/3)
630 0 : DO j = -1, 1
631 : ! j = -1 corresponds to x-direction
632 : ! j = 0 corresponds to z-direction
633 : ! j = 1 corresponds to y-direction
634 0 : mj = abs(j)
635 0 : cfac = img**((j + mj)/2.)*sqrt(2.)**(1 - mj)*rfac
636 0 : gauntvec(j) = cfac*(gaunt1(1, l1, l2, -mj, m1, m2, atoms%lmaxd + 1) + j*gaunt1(1, l1, l2, mj, m1, m2, atoms%lmaxd + 1))
637 : END DO
638 :
639 : ! transform onto cartesian coordinates
640 0 : point(-1) = -1
641 0 : point(0) = 1
642 0 : point(1) = 0
643 :
644 0 : gauntvec = gauntvec(point)
645 :
646 : END FUNCTION gauntvec
647 :
648 0 : FUNCTION w(p1, l1, p2, l2, itype, bas1_mt, drbas1_mt, &
649 0 : rmt)
650 : USE m_types
651 : USE m_juDFT
652 : IMPLICIT NONE
653 :
654 : INTEGER, INTENT(IN) :: p1, l1, p2, l2
655 : INTEGER, INTENT(IN) :: itype
656 : REAL, INTENT(IN) :: rmt(:), bas1_mt(:, 0:, :), drbas1_mt(:, 0:, :)
657 :
658 : REAL :: w
659 :
660 : INTEGER :: p
661 : REAL :: denom, var_enum
662 :
663 0 : IF(p1 > 2 .or. p2 > 2) call judft_error('w: the formalism is only valid for p<=2')
664 :
665 0 : denom = wronskian(bas1_MT(2, l1, itype), drbas1_MT(2, l1, itype), bas1_MT(1, l1, itype), drbas1_MT(1, l1, itype))
666 :
667 0 : p = p1 + (-1)**(p1 - 1)
668 :
669 : var_enum = bas1_MT(p, l1, itype)*bas1_MT(p2, l2, itype) + rmt(itype)*wronskian(bas1_MT(p, l1, itype), &
670 0 : drbas1_MT(p, l1, itype), bas1_MT(p2, l2, itype), drbas1_MT(p2, l2, itype))
671 :
672 0 : w = (-1)**(p1 + 1)*var_enum/denom
673 :
674 0 : END FUNCTION
675 :
676 0 : PURE FUNCTION wronskian(f, df, g, dg)
677 : IMPLICIT NONE
678 : REAL, INTENT(IN) :: f, df, g, dg
679 : REAL :: wronskian
680 :
681 0 : wronskian = f*dg - df*g
682 :
683 0 : END FUNCTION
684 :
685 : ! Calculates the derivative
686 : ! ikr s s
687 : ! dcprod(n',n,k,xyz) = d < e phi | phi > / sqrt(vol)
688 : ! xyz qn q+k,n'
689 : !
690 : ! s s
691 : ! < phi | d | phi >
692 : ! nq xyz n'q
693 : ! = -i ------------------------ / sqrt(vol) , s = ispin
694 : ! s s , n = occ. , n' = unocc.
695 : ! e - e , bandi1 <= n <= bandf1 , bandi2 <= n' <= bandf2
696 : ! n'q nq
697 : !
698 : ! with kp perturbation theory and
699 : !
700 : ! d d d
701 : ! d = ---, ---, --- .
702 : ! xyz dk dk dk
703 : ! x y z
704 : !
705 0 : SUBROUTINE dwavefproducts( &
706 0 : dcprod, nk, bandi1, bandf1, bandi2, bandf2, lwrite, &
707 : input, atoms, mpdata, hybinp, &
708 : cell, &
709 : hybdat, kpts, sym, noco, nococonv, lapw, &
710 : jsp, &
711 0 : eig_irr)
712 :
713 : USE m_wrapper
714 : USE m_types
715 : use m_constants, only: cmplx_0
716 : IMPLICIT NONE
717 :
718 : TYPE(t_hybdat), INTENT(IN) :: hybdat
719 : TYPE(t_input), INTENT(IN) ::input
720 : TYPE(t_mpdata), intent(in) :: mpdata
721 : TYPE(t_hybinp), INTENT(IN) :: hybinp
722 : TYPE(t_cell), INTENT(IN) :: cell
723 : TYPE(t_kpts), INTENT(IN) :: kpts
724 : type(t_sym), intent(in) :: sym
725 : type(t_noco), intent(in) :: noco
726 : type(t_nococonv), intent(in) :: nococonv
727 : TYPE(t_atoms), INTENT(IN) :: atoms
728 : TYPE(t_lapw), INTENT(IN) :: lapw
729 :
730 :
731 : ! - scalars -
732 : INTEGER, INTENT(IN) :: nk, bandi1, bandf1, bandi2, bandf2
733 : INTEGER, INTENT(IN) :: jsp
734 :
735 : ! - arrays -
736 :
737 : REAL, INTENT(IN) :: eig_irr(:, :)
738 : COMPLEX, INTENT(INOUT) :: dcprod(bandi2:bandf2, bandi1:bandf1, 3)
739 :
740 : ! - local scalars -
741 : INTEGER :: ikpt, ikpt1, iband1, iband2
742 : REAL :: rdum
743 : LOGICAL :: lwrite
744 :
745 : ! __
746 : ! Get momentum-matrix elements -i < uj | \/ | ui >
747 : !
748 0 : dcprod = cmplx_0
749 : CALL momentum_matrix(dcprod, nk, bandi1, bandf1, bandi2, bandf2, &
750 : input, atoms, mpdata, hybinp, cell, hybdat, kpts, sym, noco,nococonv, lapw, &
751 0 : jsp)
752 :
753 : ! __
754 : ! Calculate expansion coefficients -i < uj | \/ | ui > / ( ei - ej ) for periodic function ui
755 : !
756 0 : DO iband1 = bandi1, bandf1
757 0 : DO iband2 = bandi2, bandf2
758 0 : rdum = eig_irr(iband2, nk) - eig_irr(iband1, nk)
759 0 : IF(abs(rdum) > 1e-6) THEN !10.0**-6
760 0 : dcprod(iband2, iband1, :) = dcprod(iband2, iband1, :)/rdum
761 : ELSE
762 0 : dcprod(iband2, iband1, :) = 0.0
763 : END IF
764 : END DO
765 : END DO
766 :
767 0 : dcprod = dcprod/sqrt(cell%omtil)
768 :
769 0 : END SUBROUTINE dwavefproducts
770 :
771 : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
772 :
773 : ! Calculates the momentum matrix elements
774 : !
775 : ! s s
776 : ! momentum(n',n,q,xyz) = -i < phi | d | phi > , s = ispin
777 : ! nq xyz n'q , n = occ. ( bandi1 <= n <= bandf1 ) ,
778 : ! n' = unocc. ( bandi2 <= n' <= bandf2 )
779 : !
780 : !
781 0 : SUBROUTINE momentum_matrix(momentum, nk, bandi1, bandf1, bandi2, bandf2, &
782 : input, atoms, mpdata, hybinp, &
783 : cell, hybdat, kpts, sym, noco,nococonv, lapw, jsp)
784 : USE m_olap
785 : USE m_wrapper
786 : USE m_util, only: derivative
787 : use m_intgrf, only: intgrf_init, intgrf
788 : USE m_dr2fdr
789 : USE m_constants
790 : USE m_types
791 : USE m_io_hybrid
792 : use m_calc_cmt
793 : IMPLICIT NONE
794 : TYPE(t_input), INTENT(IN) :: input
795 : TYPE(t_hybdat), INTENT(IN) :: hybdat
796 : TYPE(t_mpdata), intent(in) :: mpdata
797 : TYPE(t_hybinp), INTENT(IN) :: hybinp
798 : TYPE(t_cell), INTENT(IN) :: cell
799 : TYPE(t_kpts), INTENT(IN) :: kpts
800 : type(t_sym), intent(in) :: sym
801 : type(t_noco), intent(in) :: noco
802 : type(t_nococonv), intent(in)::nococonv
803 : TYPE(t_atoms), INTENT(IN) :: atoms
804 : TYPE(t_lapw), INTENT(IN) :: lapw
805 :
806 :
807 : ! - scalars -
808 : INTEGER, INTENT(IN) :: bandi1, bandf1, bandi2, bandf2
809 : INTEGER, INTENT(IN) :: nk
810 : INTEGER, INTENT(IN) :: jsp
811 :
812 : ! - arrays -
813 0 : TYPE(t_mat):: z
814 : COMPLEX, INTENT(INOUT) :: momentum(bandi2:bandf2, bandi1:bandf1, 3)
815 :
816 : ! - local scalars -
817 : INTEGER :: itype, ieq, ic, i, j, l, lm, n1, n2, ikpt, iband1, iband2, ll, mm
818 : INTEGER :: lm_0, lm_1, lm0, lm1, lm2, lm3, n0, nn, n, l1, l2, m1, m2, ikpt1
819 : INTEGER :: irecl_cmt, irecl_z, m
820 : COMPLEX :: cdum
821 : COMPLEX :: img = (0.0, 1.0)
822 :
823 : ! - local arrays -
824 0 : INTEGER :: gpt(3, lapw%nv(jsp))
825 :
826 0 : REAL :: fcoeff((atoms%lmaxd + 1)**2, -1:1), gcoeff((atoms%lmaxd + 1)**2, -1:1)
827 0 : REAL :: qmat1(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype), dbas1(atoms%jmtd)
828 0 : REAL :: qmat2(maxval(mpdata%num_radfun_per_l), maxval(mpdata%num_radfun_per_l), 0:atoms%lmaxd, atoms%ntype), dbas2(atoms%jmtd)
829 0 : REAL :: qg(lapw%nv(jsp), 3)
830 :
831 : COMPLEX :: hlp(3, 3)
832 0 : COMPLEX :: cvec1(hybdat%maxlmindx), cvec2(hybdat%maxlmindx), cvec3(hybdat%maxlmindx)
833 0 : COMPLEX :: cmt1(hybdat%maxlmindx, bandi1:bandf1), cmt2(hybdat%maxlmindx, bandi2:bandf2)
834 : COMPLEX :: carr1(3), carr2(3)
835 0 : COMPLEX :: cmt(input%neig, hybdat%maxlmindx, atoms%nat)
836 0 : REAL :: olap_r(lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2)
837 0 : COMPLEX :: olap_c(lapw%nv(jsp)*(lapw%nv(jsp) + 1)/2)
838 0 : REAL :: vec1_r(lapw%nv(jsp)), vec2_r(lapw%nv(jsp)), vec3_r(lapw%nv(jsp))
839 0 : COMPLEX :: vec1_c(lapw%nv(jsp)), vec2_c(lapw%nv(jsp)), vec3_c(lapw%nv(jsp))
840 0 : COMPLEX :: c_phase(hybdat%nbands(nk,jsp))
841 :
842 : ! read in cmt coefficients from direct access file cmt at kpoint nk
843 0 : momentum = cmplx_0
844 :
845 0 : if(nk /= kpts%bkp(nk)) call juDFT_error("We should be reading the parent z-mat here!")
846 0 : call read_z(atoms, cell, hybdat, kpts, sym, noco, nococonv, input, nk, jsp, z, c_phase=c_phase)
847 : call calc_cmt(atoms, cell, input, noco, nococonv, hybinp, hybdat, mpdata, kpts, &
848 0 : sym, z, jsp, nk, c_phase, cmt)
849 :
850 : ! read in z coefficients from direct access file z at kpoint nk
851 :
852 :
853 : !CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,hybdat%gridf)
854 0 : gpt(:, 1:lapw%nv(jsp)) = lapw%gvec(:, 1:lapw%nv(jsp), jsp)
855 :
856 : ! Define coefficients F and G
857 : lm = 0
858 0 : DO l = 0, atoms%lmaxd
859 0 : DO M = -l, l
860 0 : lm = lm + 1
861 0 : fcoeff(lm, -1) = -sqrt(1.0*(l + M + 1)*(l + M + 2)/(2*(2*l + 1)*(2*l + 3)))
862 0 : fcoeff(lm, 0) = sqrt(1.0*(l - M + 1)*(l + M + 1)/((2*l + 1)*(2*l + 3)))
863 0 : fcoeff(lm, 1) = -sqrt(1.0*(l - M + 1)*(l - M + 2)/(2*(2*l + 1)*(2*l + 3)))
864 0 : gcoeff(lm, -1) = sqrt(1.0*(l - M)*(l - M - 1)/(2*(2*l - 1)*(2*l + 1)))
865 0 : gcoeff(lm, 0) = sqrt(1.0*(l - M)*(l + M)/((2*l - 1)*(2*l + 1)))
866 0 : gcoeff(lm, 1) = sqrt(1.0*(l + M)*(l + M - 1)/(2*(2*l - 1)*(2*l + 1)))
867 : END DO
868 : END DO
869 :
870 : ! Calculate olap int r**2*u*u' + w * int r*u*u, w = -l,l+1 ( -> qmat1/2 )
871 0 : qmat1 = 0
872 0 : qmat2 = 0
873 0 : ic = 0
874 0 : DO itype = 1, atoms%ntype
875 0 : DO l = 0, atoms%lmax(itype)
876 0 : DO n2 = 1, mpdata%num_radfun_per_l(l, itype)
877 : !ic = ic + 1
878 0 : CALL derivative(dbas1, hybdat%bas1(:, n2, l, itype), atoms, itype)
879 0 : dbas1 = dbas1 - hybdat%bas1(:, n2, l, itype)/atoms%rmsh(:, itype)
880 :
881 0 : CALL derivative(dbas2, hybdat%bas2(:, n2, l, itype), atoms, itype)
882 0 : dbas2 = dbas2 - hybdat%bas2(:, n2, l, itype)/atoms%rmsh(:, itype)
883 :
884 0 : IF(l /= 0) THEN
885 0 : DO n1 = 1, mpdata%num_radfun_per_l(l - 1, itype)
886 0 : ic = ic + 1
887 : qmat1(n1, n2, l, itype) = intgrf(dbas1(:)*hybdat%bas1(:, n1, l - 1, itype) + &
888 : dbas2(:)*hybdat%bas2(:, n1, l - 1, itype), atoms, itype, hybdat%gridf) &
889 : + intgrf((hybdat%bas1(:, n2, l, itype)*hybdat%bas1(:, n1, l - 1, itype) + hybdat%bas2(:, n2, l, itype)*hybdat%bas2(:, n1, l - 1, itype)) &
890 0 : /atoms%rmsh(:, itype), atoms, itype, hybdat%gridf)*(l + 1)
891 :
892 : END DO
893 : END IF
894 0 : IF(l /= atoms%lmax(itype)) THEN
895 0 : DO n1 = 1, mpdata%num_radfun_per_l(l + 1, itype)
896 :
897 : qmat2(n1, n2, l, itype) = intgrf(dbas1(:)*hybdat%bas1(:, n1, l + 1, itype) + dbas2(:)*hybdat%bas2(:, n1, l + 1, itype), &
898 : atoms, itype, hybdat%gridf) &
899 : - intgrf((hybdat%bas1(:, n2, l, itype)*hybdat%bas1(:, n1, l + 1, itype) + hybdat%bas2(:, n2, l, itype)*hybdat%bas2(:, n1, l + 1, itype)) &
900 0 : /atoms%rmsh(:, itype), atoms, itype, hybdat%gridf)*l
901 :
902 : END DO
903 : END IF
904 :
905 : END DO
906 : END DO
907 : END DO
908 :
909 : ! __
910 : ! Calculate momentum matrix elements -i < uj | \/ | ui > wrt wave functions u (->momentum)
911 : !
912 :
913 0 : momentum = 0
914 :
915 : ! MT contribution
916 :
917 0 : ic = 0
918 0 : DO itype = 1, atoms%ntype
919 0 : DO ieq = 1, atoms%neq(itype)
920 0 : ic = ic + 1
921 0 : nn = sum([((2*l + 1)*mpdata%num_radfun_per_l(l, itype), l=0, atoms%lmax(itype))])
922 0 : DO iband1 = bandi1, bandf1
923 0 : cmt1(:nn, iband1) = cmt(iband1, :nn, ic)
924 : ENDDO
925 0 : DO iband2 = bandi2, bandf2
926 0 : cmt2(:nn, iband2) = cmt(iband2, :nn, ic)
927 : ENDDO
928 0 : DO iband1 = bandi1, bandf1
929 :
930 0 : cvec1 = 0; cvec2 = 0; cvec3 = 0
931 : ! build up left vector(s) ( -> cvec1/2/3 )
932 0 : lm_0 = 0 ! we start with s-functions (l=0)
933 0 : lm_1 = mpdata%num_radfun_per_l(0, itype) ! we start with p-functions (l=0+1)
934 0 : lm = 0
935 0 : DO l = 0, atoms%lmax(itype) - 1
936 0 : n0 = mpdata%num_radfun_per_l(l, itype)
937 0 : n1 = mpdata%num_radfun_per_l(l + 1, itype)
938 0 : DO M = -l, l
939 0 : lm = lm + 1
940 0 : lm0 = lm_0 + (M + l)*n0
941 0 : lm1 = lm_1 + (M + 1 + l + 1)*n1
942 0 : lm2 = lm_1 + (M + l + 1)*n1
943 0 : lm3 = lm_1 + (M - 1 + l + 1)*n1
944 0 : cvec1(lm0 + 1:lm0 + n0) = fcoeff(lm, -1)*matmul(cmt1(lm1 + 1:lm1 + n1, iband1), qmat2(:n1, :n0, l, itype))
945 0 : cvec2(lm0 + 1:lm0 + n0) = fcoeff(lm, 0)*matmul(cmt1(lm2 + 1:lm2 + n1, iband1), qmat2(:n1, :n0, l, itype))
946 0 : cvec3(lm0 + 1:lm0 + n0) = fcoeff(lm, 1)*matmul(cmt1(lm3 + 1:lm3 + n1, iband1), qmat2(:n1, :n0, l, itype))
947 : END DO
948 0 : lm_0 = lm_0 + (2*l + 1)*n0
949 0 : lm_1 = lm_1 + (2*l + 3)*n1
950 : END DO
951 :
952 0 : lm_0 = mpdata%num_radfun_per_l(0, itype) ! we start with p-functions (l=1)
953 0 : lm_1 = 0 ! we start with s-functions (l=1-1)
954 0 : lm = 1
955 0 : DO l = 1, atoms%lmax(itype)
956 0 : n0 = mpdata%num_radfun_per_l(l, itype)
957 0 : n1 = mpdata%num_radfun_per_l(l - 1, itype)
958 0 : DO M = -l, l
959 0 : lm = lm + 1
960 0 : lm0 = lm_0 + (M + l)*n0
961 0 : lm1 = lm_1 + (M + 1 + l - 1)*n1
962 0 : lm2 = lm_1 + (M + l - 1)*n1
963 0 : lm3 = lm_1 + (M - 1 + l - 1)*n1
964 0 : IF(abs(M + 1) <= l - 1) cvec1(lm0 + 1:lm0 + n0) = cvec1(lm0 + 1:lm0 + n0) + gcoeff(lm, -1)*matmul(cmt1(lm1 + 1:lm1 + n1, iband1), qmat1(:n1, :n0, l, itype))
965 0 : IF(abs(M) <= l - 1) cvec2(lm0 + 1:lm0 + n0) = cvec2(lm0 + 1:lm0 + n0) + gcoeff(lm, 0)*matmul(cmt1(lm2 + 1:lm2 + n1, iband1), qmat1(:n1, :n0, l, itype))
966 0 : IF(abs(M - 1) <= l - 1) cvec3(lm0 + 1:lm0 + n0) = cvec3(lm0 + 1:lm0 + n0) + gcoeff(lm, 1)*matmul(cmt1(lm3 + 1:lm3 + n1, iband1), qmat1(:n1, :n0, l, itype))
967 : END DO
968 0 : lm_0 = lm_0 + (2*l + 1)*n0
969 0 : lm_1 = lm_1 + (2*l - 1)*n1
970 : END DO
971 : ! multiply with right vector
972 0 : DO iband2 = bandi2, bandf2
973 0 : momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dot_product(cvec1(:nn), cmt2(:nn, iband2))
974 0 : momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dot_product(cvec2(:nn), cmt2(:nn, iband2))
975 0 : momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dot_product(cvec3(:nn), cmt2(:nn, iband2))
976 : END DO ! iband2
977 : END DO ! iband1
978 :
979 : END DO ! ieq
980 : END DO ! itype
981 :
982 : ! Transform to cartesian coordinates
983 0 : hlp = 0
984 0 : hlp(1, 1) = 1.0/sqrt(2.0)
985 0 : hlp(1, 3) = -1.0/sqrt(2.0)
986 0 : hlp(2, 1) = -img/sqrt(2.0)
987 0 : hlp(2, 3) = -img/sqrt(2.0)
988 0 : hlp(3, 2) = 1.0
989 0 : DO iband1 = bandi1, bandf1
990 0 : DO iband2 = bandi2, bandf2
991 0 : momentum(iband2, iband1, :) = -img*matmul(momentum(iband2, iband1, :), transpose(hlp))
992 : END DO
993 : END DO
994 :
995 : ! plane-wave contribution
996 :
997 0 : CALL olap_pwp(z%l_real, olap_r, olap_c, gpt, lapw%nv(jsp), atoms, cell)
998 :
999 0 : DO nn = 1, lapw%nv(jsp)
1000 0 : qg(nn, :) = matmul(kpts%bk(:, nk) + gpt(:, nn), cell%bmat)
1001 : END DO
1002 :
1003 0 : if(z%l_real) THEN
1004 0 : DO iband2 = bandi2, bandf2
1005 0 : vec1_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 1))
1006 0 : vec2_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 2))
1007 0 : vec3_r = matvec(olap_r, z%data_r(:lapw%nv(jsp), iband2)*qg(:, 3))
1008 0 : DO iband1 = bandi1, bandf1
1009 0 : momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dot_product(z%data_r(:lapw%nv(jsp), iband1), vec1_r)
1010 0 : momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dot_product(z%data_r(:lapw%nv(jsp), iband1), vec2_r)
1011 0 : momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dot_product(z%data_r(:lapw%nv(jsp), iband1), vec3_r)
1012 : END DO
1013 : END DO
1014 : else
1015 0 : DO iband2 = bandi2, bandf2
1016 0 : vec1_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 1))
1017 0 : vec2_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 2))
1018 0 : vec3_c = matvec(olap_c, z%data_c(:lapw%nv(jsp), iband2)*qg(:, 3))
1019 0 : DO iband1 = bandi1, bandf1
1020 0 : momentum(iband2, iband1, 1) = momentum(iband2, iband1, 1) + dot_product(z%data_c(:lapw%nv(jsp), iband1), vec1_c)
1021 0 : momentum(iband2, iband1, 2) = momentum(iband2, iband1, 2) + dot_product(z%data_c(:lapw%nv(jsp), iband1), vec2_c)
1022 0 : momentum(iband2, iband1, 3) = momentum(iband2, iband1, 3) + dot_product(z%data_c(:lapw%nv(jsp), iband1), vec3_c)
1023 : END DO
1024 : END DO
1025 : end if
1026 :
1027 0 : END SUBROUTINE momentum_matrix
1028 :
1029 0 : END MODULE m_kp_perturbation
|
