Line data Source code
1 : !
2 : ! Calculates the Coulomb matrix
3 : !
4 : ! v = < M | v | M >
5 : ! k,IJ k,I k,J
6 : !
7 : ! with the mixed-basis functions M (indices I and J).
8 : !
9 : ! Note that
10 : ! *
11 : ! v = v .
12 : ! k,JI k,IJ
13 : !
14 : ! In the code: coulomb(IJ,k) = v where only the upper triangle (I<=J) is stored.
15 : ! k,IJ
16 : !
17 : ! The Coulomb matrix v(IJ,k) diverges at the Gamma-point. Here, we apply the decomposition
18 : !
19 : ! (0) (1) * 2-l (0)* (0) (1)* m (1)
20 : ! v = v + SUM v * Y (k) / k with v = v , v = (-1) v
21 : ! k,IJ IJ lm IJ lm JI IJ JI,lm IJ,l,-m
22 : !
23 : ! where a = atom index, R = position vector, T = Wigner-Seitz radius (scalar).
24 : ! a 0
25 : ! (0)
26 : ! In the code: coulomb(IJ,1) = v where only the upper triangle (I<=J) is stored,
27 : ! IJ
28 : ! (1)
29 : ! coulfac(IJ,lm) = v IJ,lm
30 : !
31 : ! For the PW contribution we have to construct plane waves within the MT spheres with the help
32 : ! of spherical Bessel functions. The value lexp (LEXP in gwinp) is the corresponding cutoff.
33 : !
34 : MODULE m_coulombmatrix
35 : #ifdef CPP_MPI
36 : use mpi
37 : #endif
38 : use m_types
39 : USE m_intgrf, ONLY: intgrf, intgrf_init
40 : use m_sphbes, only: sphbes
41 : use m_glob_tofrom_loc
42 : USE m_trafo, ONLY: symmetrize_mpimat, symmetrize, bramat_trafo
43 : use m_gamma_double_gpt_loop
44 : CONTAINS
45 :
46 12 : SUBROUTINE coulombmatrix(fmpi, fi, mpdata, hybdat, xcpot)
47 : use m_work_package
48 : use m_structureconstant
49 : USE m_types
50 : USE m_types_mpimat
51 : use m_types_mat
52 : USE m_types_hybdat
53 : USE m_juDFT
54 : USE m_constants
55 : use m_util, only: primitivef
56 : USE m_hsefunctional, ONLY: change_coulombmatrix
57 : USE m_wrapper
58 : USE m_io_hybrid
59 : use m_ylm
60 : use m_calc_l_m_from_lm
61 : use m_calc_mpsmat
62 : use m_copy_coul
63 : use m_apply_inverse_olap
64 : IMPLICIT NONE
65 :
66 : TYPE(t_xcpot_inbuild), INTENT(IN) :: xcpot
67 : TYPE(t_mpi), INTENT(IN) :: fmpi
68 : type(t_fleurinput), intent(in) :: fi
69 : TYPE(t_mpdata), intent(in) :: mpdata
70 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
71 :
72 : ! - local scalars -
73 : INTEGER :: inviop
74 : INTEGER :: nqnrm, iqnrm, iqnrm1, iqnrm2
75 : INTEGER :: itype, l, ix, iy, iy0, i, j, lm, l1, l2, m1, m2, ineq, ikpt
76 : INTEGER :: lm1, lm2, itype1, itype2, ineq1, ineq2, n1, n2, iat2
77 : INTEGER :: ic, ic1, ic2
78 : INTEGER :: igpt, igpt1, igpt2, igptp, igptp1, igptp2
79 : INTEGER :: isym, isym1, isym2, igpt0
80 : INTEGER :: iatm1, iatm2
81 : INTEGER :: m, im
82 : INTEGER :: maxfac, ix_loc, pe, pe_ix
83 :
84 : LOGICAL :: lsym
85 :
86 : REAL :: rdum, rdum1, rdum2
87 : REAL :: svol, qnorm1, qnorm2
88 : REAL :: fcoulfac
89 :
90 : COMPLEX :: cdum, cexp, csum
91 :
92 : ! - local arrays -
93 : INTEGER :: g(3)
94 12 : INTEGER, ALLOCATABLE :: pqnrm(:, :)
95 24 : INTEGER :: rrot(3, 3, fi%sym%nsym), invrrot(3, 3, fi%sym%nsym)
96 24 : INTEGER, ALLOCATABLE :: iarr(:), POINTER(:, :, :, :)!,pointer(:,:,:)
97 24 : INTEGER, ALLOCATABLE :: nsym_gpt(:, :), sym_gpt(:, :, :)
98 24 : INTEGER :: nsym1(fi%kpts%nkpt + 1), sym1(fi%sym%nsym, fi%kpts%nkpt + 1)
99 :
100 12 : INTEGER, ALLOCATABLE :: ngptm1(:)
101 12 : INTEGER, ALLOCATABLE :: pgptm1(:, :)
102 :
103 : REAL :: q(3), q1(3), q2(3), mtmt_term
104 24 : REAL :: integrand(fi%atoms%jmtd), primf1(fi%atoms%jmtd), primf2(fi%atoms%jmtd)
105 152 : REAL :: moment(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype), &
106 132 : moment2(maxval(mpdata%num_radbasfn), fi%atoms%ntype)
107 24 : REAL, ALLOCATABLE :: gmat(:, :), qnrm(:)
108 12 : REAL, ALLOCATABLE :: sphbesmoment(:, :, :)
109 12 : REAL, ALLOCATABLE :: sphbes0(:, :, :)
110 24 : REAL, ALLOCATABLE :: olap(:, :, :, :), integral(:, :, :, :)
111 12 : REAL, ALLOCATABLE :: gridf(:, :)
112 64 : REAL :: facA(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
113 64 : REAL :: facB(0:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
114 64 : REAL :: facC(-1:MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1))
115 32 : REAL :: sphbes_var(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
116 32 : REAL :: sphbesmoment1(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
117 :
118 12 : COMPLEX :: y((fi%hybinp%lexp + 1)**2), smat
119 112 : COMPLEX :: dwgn(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), 0:maxval(fi%hybinp%lcutm1), fi%sym%nsym)
120 12 : COMPLEX, ALLOCATABLE :: carr2(:, :), carr2a(:, :), carr2b(:, :)
121 12 : COMPLEX, ALLOCATABLE :: structconst(:,:,:,:)
122 :
123 : INTEGER :: ierr
124 : INTEGER :: iatom, mtmt_idx
125 12 : TYPE(t_mat) :: mat
126 12 : type(t_mat), allocatable :: mtmt_repl(:)
127 12 : class(t_mat), allocatable :: coul(:)
128 :
129 12 : CALL timestart("Coulomb matrix setup")
130 12 : call timestart("prep in coulomb")
131 12 : if (fmpi%is_root()) write (*, *) "start of coulomb calculation"
132 :
133 72 : allocate(structconst((2*fi%hybinp%lexp + 1)**2, fi%atoms%nat, fi%atoms%nat, fi%kpts%nkpt), stat=ierr)
134 12 : if(ierr /= 0) call judft_error("can't allocate structconst. error: " // int2str(ierr))
135 :
136 12 : svol = SQRT(fi%cell%vol)
137 12 : fcoulfac = fpi_const/fi%cell%vol
138 32 : maxfac = MAX(2*fi%atoms%lmaxd + maxval(fi%hybinp%lcutm1) + 1, 4*MAX(maxval(fi%hybinp%lcutm1), fi%hybinp%lexp) + 1)
139 :
140 12 : facA(0) = 1 !
141 12 : facB(0) = 1 ! Define:
142 36 : facC(-1:0) = 1 ! facA(i) = i!
143 792 : DO i = 1, maxfac ! facB(i) = sqrt(i!)
144 780 : facA(i) = facA(i - 1)*i ! facC(i) = (2i+1)!!
145 780 : facB(i) = facB(i - 1)*SQRT(i*1.0) !
146 792 : facC(i) = facC(i - 1)*(2*i + 1) !
147 : END DO
148 :
149 12 : CALL intgrf_init(fi%atoms%ntype, fi%atoms%jmtd, fi%atoms%jri, fi%atoms%dx, fi%atoms%rmsh, gridf)
150 :
151 : ! Calculate the structure constant
152 12 : CALL structureconstant(structconst, fi%cell, fi%hybinp, fi%atoms, fi%kpts, fmpi)
153 :
154 12 : IF (fmpi%irank == 0) WRITE (oUnit, '(//A)') '### subroutine: coulombmatrix ###'
155 :
156 : !
157 : ! Matrix allocation
158 : !
159 :
160 12 : call timestart("coulomb allocation")
161 :
162 12 : if(fmpi%n_size == 1) then
163 0 : allocate(t_mat::coul(fi%kpts%nkpt))
164 : else
165 72 : allocate(t_mpimat::coul(fi%kpts%nkpt))
166 : endif
167 48 : do ikpt = 1, fi%kpts%nkpt
168 84 : if(any(ikpt == fmpi%k_list))then
169 36 : call coul(ikpt)%init(.False., hybdat%nbasm(ikpt), hybdat%nbasm(ikpt), fmpi%sub_comm, .false.)
170 : else
171 0 : call coul(ikpt)%init(.False., 1, 1, fmpi%sub_comm, .false.)
172 : endif
173 : enddo
174 12 : call timestop("coulomb allocation")
175 :
176 12 : IF (fmpi%irank == 0) then
177 : write (oUnit,*) "Size of coulomb matrix: " //&
178 42 : float2str(sum([(coul(fmpi%k_list(i))%size_mb(), i=1,size(fmpi%k_list))])) // " MB"
179 : endif
180 :
181 : ! Generate Symmetry:
182 : ! Reduce list of g-Points so that only one of each symm-equivalent is calculated
183 : ! calculate rotations in reciprocal space
184 588 : DO isym = 1, fi%sym%nsym
185 588 : IF (isym <= fi%sym%nop) THEN
186 480 : inviop = fi%sym%invtab(isym)
187 6240 : rrot(:, :, isym) = TRANSPOSE(fi%sym%mrot(:, :, inviop))
188 4128 : DO l = 0, maxval(fi%hybinp%lcutm1)
189 : dwgn(:, :, l, isym) = TRANSPOSE(fi%hybinp%d_wgn2(-maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), &
190 228960 : -maxval(fi%hybinp%lcutm1):maxval(fi%hybinp%lcutm1), l, isym))
191 : END DO
192 : ELSE
193 96 : inviop = isym - fi%sym%nop
194 1248 : rrot(:, :, isym) = -rrot(:, :, inviop)
195 43776 : dwgn(:, :, :, isym) = dwgn(:, :, :, inviop)
196 864 : DO l = 0, maxval(fi%hybinp%lcutm1)
197 2976 : DO m1 = -l, l
198 9600 : DO m2 = -l, -1
199 6720 : cdum = dwgn(m1, m2, l, isym)
200 6720 : dwgn(m1, m2, l, isym) = dwgn(m1, -m2, l, isym)*(-1)**m2
201 9120 : dwgn(m1, -m2, l, isym) = cdum*(-1)**m2
202 : END DO
203 : END DO
204 : END DO
205 : END IF
206 : END DO
207 6744 : invrrot(:, :, :fi%sym%nop) = rrot(:, :, fi%sym%invtab)
208 12 : IF (fi%sym%nsym > fi%sym%nop) THEN
209 1348 : invrrot(:, :, fi%sym%nop + 1:) = rrot(:, :, fi%sym%invtab + fi%sym%nop)
210 : END IF
211 :
212 : ! Get symmetry operations that leave bk(:,ikpt) invariant -> sym1
213 60 : nsym1 = 0
214 48 : DO ikpt = 1, fi%kpts%nkpt
215 36 : isym1 = 0
216 1764 : DO isym = 1, fi%sym%nsym
217 : ! temporary fix until bramat_trafo is correct
218 : ! for systems with symmetries including translations
219 1728 : IF (isym > fi%sym%nop) THEN
220 288 : isym2 = isym - fi%sym%nop
221 : ELSE
222 : isym2 = isym
223 : END IF
224 6048 : IF (ANY(abs(fi%sym%tau(:, isym2)) > 1e-12)) CYCLE
225 :
226 34404 : IF (ALL(ABS(MATMUL(rrot(:, :, isym), fi%kpts%bk(:, ikpt)) - fi%kpts%bk(:, ikpt)) < 1e-12)) THEN
227 544 : isym1 = isym1 + 1
228 544 : sym1(isym1, ikpt) = isym
229 : END IF
230 : END DO
231 48 : nsym1(ikpt) = isym1
232 : END DO
233 : ! Define reduced lists of G points -> pgptm1(:,ikpt), ikpt=1,..,nkpt
234 : !if(allocated(pgptm1)) deallocate(fi%hybinp%pgptm1)
235 14512 : allocate (pgptm1(maxval(mpdata%n_g), fi%kpts%nkptf), source=0) !in mixedbasis
236 1916 : allocate (iarr(maxval(mpdata%n_g)), source=0)
237 : allocate (POINTER(fi%kpts%nkpt, &
238 : MINVAL(mpdata%g(1, :)) - 1:MAXVAL(mpdata%g(1, :)) + 1, &
239 : MINVAL(mpdata%g(2, :)) - 1:MAXVAL(mpdata%g(2, :)) + 1, &
240 : MINVAL(mpdata%g(3, :)) - 1:MAXVAL(mpdata%g(3, :)) + 1), &
241 58292 : source=0)
242 36 : allocate (ngptm1, mold=mpdata%n_g)
243 108 : ngptm1 = 0
244 :
245 48 : DO ikpt = 1, fi%kpts%nkpt
246 4952 : DO igpt = 1, mpdata%n_g(ikpt)
247 19664 : g = mpdata%g(:, mpdata%gptm_ptr(igpt, ikpt))
248 4952 : POINTER(ikpt, g(1), g(2), g(3)) = igpt
249 : END DO
250 5388 : iarr = 0
251 : j = 0
252 4952 : DO igpt = mpdata%n_g(ikpt), 1, -1
253 4952 : IF (iarr(igpt) == 0) THEN
254 1432 : j = j + 1
255 1432 : pgptm1(j, ikpt) = igpt
256 10936 : DO isym1 = 1, nsym1(ikpt)
257 123552 : g = MATMUL(rrot(:, :, sym1(isym1, ikpt)), mpdata%g(:, mpdata%gptm_ptr(igpt, ikpt)))
258 9504 : i = POINTER(ikpt, g(1), g(2), g(3))
259 9504 : IF (i == 0) call judft_error('coulombmatrix: zero pointer (bug?)')
260 10936 : iarr(i) = 1
261 : END DO
262 : END IF
263 : END DO
264 48 : ngptm1(ikpt) = j
265 : END DO
266 12 : deallocate (iarr)
267 :
268 : ! Distribute the work as equally as possible over the processes
269 12 : call timestop("prep in coulomb")
270 :
271 12 : call timestart("define gmat")
272 : ! Define gmat (symmetric)
273 48 : allocate (gmat((fi%hybinp%lexp + 1)**2, (fi%hybinp%lexp + 1)**2))
274 1005732 : gmat = 0
275 : lm1 = 0
276 216 : DO l1 = 0, fi%hybinp%lexp
277 3684 : DO m1 = -l1, l1
278 3468 : lm1 = lm1 + 1
279 3468 : lm2 = 0
280 41412 : lp1: DO l2 = 0, l1
281 547332 : DO m2 = -l2, l2
282 506124 : lm2 = lm2 + 1
283 506124 : IF (lm2 > lm1) EXIT lp1 ! Don't cross the diagonal!
284 : gmat(lm1, lm2) = facB(l1 + l2 + m2 - m1)*facB(l1 + l2 + m1 - m2)/ &
285 : (facB(l1 + m1)*facB(l1 - m1)*facB(l2 + m2)*facB(l2 - m2))/ &
286 502860 : SQRT(1.0*(2*l1 + 1)*(2*l2 + 1)*(2*(l1 + l2) + 1))*(fpi_const)**1.5
287 540600 : gmat(lm2, lm1) = gmat(lm1, lm2)
288 : END DO
289 : END DO LP1
290 : END DO
291 : END DO
292 12 : call timestop("define gmat")
293 :
294 : ! Calculate moments of MT functions
295 12 : call timestart("calc moments of MT")
296 32 : DO itype = 1, fi%atoms%ntype
297 120 : DO l = 0, fi%hybinp%lcutm1(itype)
298 880 : DO i = 1, mpdata%num_radbasfn(l, itype)
299 : ! note that mpdata%radbasfn_mt already contains the factor rgrid
300 : moment(i, l, itype) = intgrf(fi%atoms%rmsh(:, itype)**(l + 1)*mpdata%radbasfn_mt(:, i, l, itype), &
301 618020 : fi%atoms, itype, gridf)
302 : END DO
303 : END DO
304 212 : DO i = 1, mpdata%num_radbasfn(0, itype)
305 : moment2(i, itype) = intgrf(fi%atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
306 146512 : fi%atoms, itype, gridf)
307 : END DO
308 : END DO
309 12 : call timestop("calc moments of MT")
310 :
311 12 : call timestart("getnorm")
312 : ! Look for different qnorm = |k+G|, definition of qnrm and pqnrm.
313 12 : CALL getnorm(fi%kpts, mpdata%g, mpdata%n_g, mpdata%gptm_ptr, qnrm, nqnrm, pqnrm, fi%cell)
314 0 : allocate (sphbesmoment(0:fi%hybinp%lexp, fi%atoms%ntype, nqnrm), &
315 0 : olap(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm), &
316 460 : integral(maxval(mpdata%num_radbasfn), 0:maxval(fi%hybinp%lcutm1), fi%atoms%ntype, nqnrm))
317 48332 : sphbes_var = 0
318 13956 : sphbesmoment = 0
319 48332 : sphbesmoment1 = 0
320 39092 : olap = 0
321 39092 : integral = 0
322 :
323 : ! Calculate moments of spherical Bessel functions (for (2) and (3)) (->sphbesmoment)
324 : ! Calculate overlap of spherical Bessel functions with basis functions (for (2)) (->olap)
325 : ! Calculate overlap of sphbesmoment1(r,l) with basis functions (for (2)) (->integral)
326 : ! We use sphbes(r,l) = j_l(qr)
327 : ! and sphbesmoment1(r,l) = 1/r**(l-1) * INT(0..r) r'**(l+2) * j_l(qr') dr'
328 : ! + r**(l+2) * INT(r..S) r'**(1-l) * j_l(qr') dr' .
329 :
330 12 : call timestop("getnorm")
331 :
332 12 : call bessel_calculation(fi, fmpi, mpdata, nqnrm, gridf, qnrm, sphbesmoment, olap, integral)
333 :
334 : !
335 : ! (1) Case < MT | v | MT >
336 : !
337 :
338 :
339 : ! (1a) r,r' in same MT
340 12 : call timestart("loop 1")
341 12 : ix = 0
342 12 : iy = 0
343 12 : iy0 = 0
344 :
345 :
346 252 : call mat%alloc(.True., maxval(mpdata%num_radbasfn), maxval(mpdata%num_radbasfn))
347 :
348 536 : allocate(mtmt_repl(calc_num_mtmts(fi)))
349 :
350 12 : mtmt_idx = 0
351 32 : DO itype = 1, fi%atoms%ntype
352 52 : DO ineq = 1, fi%atoms%neq(itype)
353 : ! Here the diagonal block matrices do not depend on ineq. In (1b) they do depend on ineq, though,
354 140 : DO l = 0, fi%hybinp%lcutm1(itype)
355 100 : mat%matsize1=mpdata%num_radbasfn(l, itype)
356 100 : mat%matsize2=mpdata%num_radbasfn(l, itype)
357 860 : DO n2 = 1, mpdata%num_radbasfn(l, itype)
358 : ! note that mpdata%radbasfn_mt already contains the factor rgrid
359 : CALL primitivef(primf1, mpdata%radbasfn_mt(:, n2, l, itype) &
360 : *fi%atoms%rmsh(:, itype)**(l + 1), fi%atoms%rmsh, fi%atoms%dx, &
361 617920 : fi%atoms%jri, fi%atoms%jmtd, itype, fi%atoms%ntype)
362 : ! -itype is to enforce inward integration
363 : CALL primitivef(primf2, mpdata%radbasfn_mt(:fi%atoms%jri(itype), n2, l, itype) &
364 : /fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l, fi%atoms%rmsh, fi%atoms%dx, &
365 606832 : fi%atoms%jri, fi%atoms%jmtd, -itype, fi%atoms%ntype)
366 :
367 606072 : primf1(:fi%atoms%jri(itype)) = primf1(:fi%atoms%jri(itype))/fi%atoms%rmsh(:fi%atoms%jri(itype), itype)**l
368 617160 : primf2 = primf2*fi%atoms%rmsh(:, itype)**(l + 1)
369 :
370 4212 : DO n1 = 1, n2
371 2708224 : integrand = mpdata%radbasfn_mt(:, n1, l, itype)*(primf1 + primf2)
372 4112 : mat%data_r(n1, n2) = fpi_const/(2*l + 1) * intgrf(integrand, fi%atoms, itype, gridf)
373 : END DO
374 : END DO
375 :
376 : ! distribute mat for m=-l,l on coulomb in block-matrix form
377 620 : DO M = -l, l
378 500 : mtmt_idx = mtmt_idx + 1
379 500 : call mtmt_repl(mtmt_idx)%alloc(.True., mpdata%num_radbasfn(l, itype), mpdata%num_radbasfn(l, itype))
380 :
381 3980 : DO n2 = 1, mpdata%num_radbasfn(l, itype)
382 18180 : DO n1 = 1, n2
383 17680 : mtmt_repl(mtmt_idx)%data_r(n1, n2) = mat%data_r(n1, n2)
384 : END DO
385 : END DO
386 600 : call mtmt_repl(mtmt_idx)%u2l()
387 : END DO
388 :
389 : END DO
390 : END DO
391 : END DO
392 12 : call mat%free()
393 12 : call timestop("loop 1")
394 :
395 48 : DO im = 1, size(fmpi%k_list)
396 36 : ikpt = fmpi%k_list(im)
397 :
398 : ! only the first rank handles the MT-MT part
399 36 : call timestart("MT-MT part")
400 36 : ix = 0
401 36 : ic2 = 0
402 36 : mtmt_idx = 0
403 96 : DO itype2 = 1, fi%atoms%ntype
404 156 : DO ineq2 = 1, fi%atoms%neq(itype2)
405 60 : ic2 = ic2 + 1
406 60 : lm2 = 0
407 420 : DO l2 = 0, fi%hybinp%lcutm1(itype2)
408 1860 : DO m2 = -l2, l2
409 1500 : mtmt_idx = mtmt_idx + 1
410 1500 : lm2 = lm2 + 1
411 12240 : DO n2 = 1, mpdata%num_radbasfn(l2, itype2)
412 10440 : ix = ix + 1
413 10440 : call glob_to_loc(fmpi, ix, pe, ix_loc)
414 11940 : if(fmpi%n_rank == pe) then
415 : iy = 0
416 : ic1 = 0
417 12420 : lp2: DO itype1 = 1, itype2
418 19620 : DO ineq1 = 1, fi%atoms%neq(itype1)
419 7200 : ic1 = ic1 + 1
420 7200 : lm1 = 0
421 50400 : DO l1 = 0, fi%hybinp%lcutm1(itype1)
422 223200 : DO m1 = -l1, l1
423 180000 : lm1 = lm1 + 1
424 1480206 : DO n1 = 1, mpdata%num_radbasfn(l1, itype1)
425 1264206 : iy = iy + 1
426 1264206 : rdum = (-1)**(l2 + m2)*moment(n1, l1, itype1)*moment(n2, l2, itype2)*gmat(lm1, lm2)
427 1264206 : l = l1 + l2
428 1264206 : lm = l**2 + l + m1 - m2 + 1
429 :
430 1264206 : if(itype2 /= itype1 .or. ineq2 /= ineq1 .or. l2 /= l1 .or. m2 /= m1) then
431 : mtmt_term = 0.0
432 : else
433 37380 : mtmt_term = mtmt_repl(mtmt_idx)%data_r(n1, n2)
434 : endif
435 :
436 : coul(ikpt)%data_c(iy, ix_loc) = mtmt_term + EXP(CMPLX(0.0, 1.0)*tpi_const* &
437 : dot_PRODUCT(fi%kpts%bk(:, ikpt), &
438 : fi%atoms%taual(:, ic2) - fi%atoms%taual(:, ic1))) &
439 5236824 : *rdum*structconst(lm, ic1, ic2, ikpt)
440 : END DO
441 : END DO
442 : END DO
443 : END DO
444 : END DO lp2
445 : endif ! pe = n_rank
446 : END DO
447 : END DO
448 : END DO
449 : END DO
450 : END DO
451 :
452 36 : IF (fi%sym%invs) THEN
453 : !symmetrize makes the Coulomb matrix real symmetric
454 : CALL symmetrize_mpimat(fi, fmpi, coul(ikpt)%data_c, [1,1],[hybdat%nbasp, hybdat%nbasp], &
455 72 : 3, .FALSE., mpdata%num_radbasfn)
456 : ENDIF
457 48 : call timestop("MT-MT part")
458 : END DO
459 :
460 108 : IF (maxval(mpdata%n_g) /= 0) THEN ! skip calculation of plane-wave contribution if mixed basis does not contain plane waves
461 :
462 : !
463 : ! (2) Case < MT | v | PW >
464 : !
465 :
466 : ! (2a) r in MT, r' everywhere
467 : ! (2b) r,r' in same MT
468 : ! (2c) r,r' in different MT
469 :
470 12 : call timestart("loop over interst.")
471 48 : DO im = 1, size(fmpi%k_list)
472 36 : ikpt = fmpi%k_list(im)
473 : call loop_over_interst(fi, hybdat, mpdata, fmpi, structconst, sphbesmoment, moment, moment2, &
474 48 : qnrm, facc, gmat, integral, olap, pqnrm, pgptm1, ngptm1, ikpt, coul(ikpt))
475 :
476 : END DO
477 :
478 12 : call timestop("loop over interst.")
479 12 : deallocate (olap, integral)
480 :
481 : !
482 : ! (3) Case < PW | v | PW >
483 : !
484 : ! (3a) r,r' everywhere; r everywhere, r' in MT; r in MT, r' everywhere
485 :
486 : ! Coulomb matrix, contribution (3a)
487 12 : call timestart("coulomb matrix 3a")
488 48 : DO im = 1, size(fmpi%k_list)
489 36 : ikpt = fmpi%k_list(im)
490 :
491 1480 : DO igpt0 = 1, ngptm1(ikpt)
492 1432 : igpt2 = pgptm1(igpt0, ikpt)
493 1432 : igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
494 1432 : ix = hybdat%nbasp + igpt2
495 1432 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
496 1468 : if(fmpi%n_rank == pe_ix) then
497 11456 : q2 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp2), fi%cell%bmat)
498 2864 : rdum2 = SUM(q2**2)
499 716 : IF (abs(rdum2) > 1e-12) rdum2 = fpi_const/rdum2
500 :
501 57126 : DO igpt1 = 1, igpt2
502 56410 : igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
503 56410 : iy = hybdat%nbasp + igpt1
504 902560 : q1 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp1), fi%cell%bmat)
505 225640 : rdum1 = SUM(q1**2)
506 56410 : IF (abs(rdum1) > 1e-12) rdum1 = fpi_const/rdum1
507 56410 : smat = calc_smat_elem(fi, mpdata, igptp1, igptp2)
508 :
509 57126 : IF (ikpt == 1) THEN
510 6542 : IF (igpt1 /= 1) THEN
511 6456 : coul(1)%data_c(iy,ix_loc) = -smat*rdum1/fi%cell%vol
512 : END IF
513 6542 : IF (igpt2 /= 1) THEN
514 : coul(1)%data_c(iy,ix_loc) &
515 6536 : = coul(1)%data_c(iy,ix_loc) - smat*rdum2/fi%cell%vol
516 : END IF
517 : ELSE
518 49868 : coul(ikpt)%data_c(iy,ix_loc) = -smat*(rdum1 + rdum2)/fi%cell%vol
519 : END IF
520 : END DO
521 716 : IF (ikpt /= 1 .OR. igpt2 /= 1) THEN
522 710 : coul(ikpt)%data_c(iy,ix_loc) = coul(ikpt)%data_c(iy,ix_loc) + rdum2
523 : END IF
524 : endif
525 : END DO
526 : END DO
527 12 : call timestop("coulomb matrix 3a")
528 : ! (3b) r,r' in different MT
529 :
530 12 : call timestart("coulomb matrix 3b")
531 48 : DO im = 1, size(fmpi%k_list)
532 36 : ikpt = fmpi%k_list(im)
533 36 : if (fmpi%is_root()) write (*, *) "coulomb pw-loop nk: ("//int2str(ikpt)//"/"//int2str(fi%kpts%nkpt)//")"
534 : ! group together quantities which depend only on l,m and igpt -> carr2a
535 828 : allocate (carr2a((fi%hybinp%lexp + 1)**2, maxval(mpdata%n_g)), carr2b(fi%atoms%nat, maxval(mpdata%n_g)))
536 1567512 : carr2a = 0; carr2b = 0
537 4952 : DO igpt = 1, mpdata%n_g(ikpt)
538 4916 : igptp = mpdata%gptm_ptr(igpt, ikpt)
539 4916 : iqnrm = pqnrm(igpt, ikpt)
540 78656 : q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
541 :
542 4916 : call ylm4(fi%hybinp%lexp, q, y)
543 :
544 1425640 : y = CONJG(y)
545 : lm = 0
546 88488 : DO l = 0, fi%hybinp%lexp
547 1509212 : DO M = -l, l
548 1420724 : lm = lm + 1
549 1504296 : carr2a(lm, igpt) = fpi_const*CMPLX(0.0, 1.0)**(l)*y(lm)
550 : END DO
551 : END DO
552 14220 : DO ic = 1, fi%atoms%nat
553 : carr2b(ic, igpt) = EXP(-CMPLX(0.0, 1.0)*tpi_const* &
554 41988 : dot_PRODUCT(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%atoms%taual(:, ic)))
555 : END DO
556 : END DO
557 :
558 : !finally we can loop over the plane waves (G: igpt1,igpt2)
559 36 : call timestart("loop over plane waves")
560 :
561 : !$OMP parallel default(none) private(carr2, igpt0, igpt1, igpt2, ix, iy, ix_loc, pe_ix, iatom, lm1, l1, m1)&
562 : !$OMP private(lm2, l2, m2, cdum, ic, lm, l, m, itype, itype2, igptp1, csum, igptp2, iqnrm1, iqnrm2, cexp)&
563 : !$OMP shared(fmpi, fi, ikpt, ngptm1, pgptm1, pqnrm, coul, gmat, structconst, sphbesmoment, hybdat, mpdata)&
564 36 : !$OMP shared(carr2a, carr2b)
565 : allocate(carr2(fi%atoms%nat, (fi%hybinp%lexp + 1)**2))
566 : !$OMP do schedule(dynamic)
567 : DO igpt0 = 1, ngptm1(ikpt)
568 : igpt2 = pgptm1(igpt0, ikpt)
569 : ix = hybdat%nbasp + igpt2
570 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
571 : if(fmpi%n_rank == pe_ix) then
572 : igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
573 : iqnrm2 = pqnrm(igpt2, ikpt)
574 :
575 : carr2 = 0
576 :
577 : ! call timestart("itype loops")
578 : do iatom = 1,fi%atoms%nat
579 : itype2 = fi%atoms%itype(iatom)
580 : cexp = CONJG(carr2b(iatom, igpt2))
581 :
582 : DO lm1 = 1, (fi%hybinp%lexp+1)**2
583 : call calc_l_m_from_lm(lm1, l1, m1)
584 : do lm2 = 1, (fi%hybinp%lexp+1)**2
585 : call calc_l_m_from_lm(lm2, l2, m2)
586 : cdum = (-1)**(l2 + m2)*sphbesmoment(l2, itype2, iqnrm2)*cexp*carr2a(lm2, igpt2)*gmat(lm1, lm2)
587 : l = l1 + l2
588 : lm = l**2 + l - l1 - m2 + (m1 + l1) + 1
589 : do iat2 =1,fi%atoms%nat
590 : carr2(iat2, lm1) = carr2(iat2,lm1) + cdum*structconst(lm, iat2, iatom, ikpt)
591 : enddo
592 : enddo
593 : enddo
594 : end do ! iatom
595 :
596 : ! call timestop("itype loops")
597 :
598 : ! call timestart("igpt1")
599 : DO igpt1 = 1, igpt2
600 : iy = hybdat%nbasp + igpt1
601 : igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
602 : iqnrm1 = pqnrm(igpt1, ikpt)
603 : csum = 0
604 : do ic = 1, fi%atoms%nat
605 : do lm = 1, (fi%hybinp%lexp+1)**2
606 : itype = fi%atoms%itype(ic)
607 : call calc_l_m_from_lm(lm, l, m)
608 : cdum = carr2b(ic, igpt1)*sphbesmoment(l, itype, iqnrm1)
609 : csum = csum + cdum*carr2(ic, lm)*CONJG(carr2a(lm, igpt1)) ! for coulomb
610 : END DO
611 : END DO
612 : coul(ikpt)%data_c(iy,ix_loc) = coul(ikpt)%data_c(iy,ix_loc) + csum/fi%cell%vol
613 : END DO
614 : ! call timestop("igpt1")
615 : endif ! pe_ix
616 : END DO !igpt0
617 : !$omp end do
618 : deallocate (carr2)
619 : !$OMP end parallel
620 36 : deallocate(carr2a, carr2b)
621 :
622 48 : call timestop("loop over plane waves")
623 : END DO !ikpt
624 12 : call timestop("coulomb matrix 3b")
625 :
626 : ! check if I own the gamma point
627 12 : if(any(fmpi%k_list == 1)) then
628 : ! Add corrections from higher orders in (3b) to coulomb(:,1)
629 : ! (1) igpt1 > 1 , igpt2 > 1 (finite G vectors)
630 12 : call timestart("add corrections from higher orders")
631 12 : call gamma_double_gpt_loop(fi, fmpi, hybdat, mpdata, sphbesmoment, gmat, ngptm1, pgptm1, pqnrm, coul(1)%data_c)
632 :
633 12 : rdum = (fpi_const)**(1.5)/fi%cell%vol**2*gmat(1, 1)
634 :
635 : ! (2) igpt1 = 1 , igpt2 > 1 (first G vector vanishes, second finite)
636 12 : call timestart("igpt1=1 loop")
637 12 : iy = hybdat%nbasp + 1
638 184 : DO igpt0 = 1, ngptm1(1)
639 172 : igpt2 = pgptm1(igpt0, 1)
640 184 : IF (igpt2 /= 1) then
641 160 : ix = hybdat%nbasp + igpt2
642 160 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
643 160 : if(fmpi%n_rank == pe_ix) then
644 80 : iqnrm2 = pqnrm(igpt2, 1)
645 80 : igptp2 = mpdata%gptm_ptr(igpt2, 1)
646 80 : qnorm2 = qnrm(iqnrm2)
647 232 : DO itype1 = 1, fi%atoms%ntype
648 384 : DO ineq1 = 1, fi%atoms%neq(itype1)
649 : ic2 = 0
650 600 : DO itype2 = 1, fi%atoms%ntype
651 744 : DO ineq2 = 1, fi%atoms%neq(itype2)
652 296 : ic2 = ic2 + 1
653 1184 : cdum = EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT(mpdata%g(:, igptp2), fi%atoms%taual(:, ic2)))
654 : coul(1)%data_c(iy, ix_loc) = coul(1)%data_c(iy, ix_loc) &
655 : + rdum*cdum*fi%atoms%rmt(itype1)**3*( &
656 : +sphbesmoment(0, itype2, iqnrm2)/30*fi%atoms%rmt(itype1)**2 &
657 : - sphbesmoment(2, itype2, iqnrm2)/18 &
658 592 : + sphbesmoment(1, itype2, iqnrm2)/6/qnorm2)
659 : END DO
660 : END DO
661 : END DO
662 : END DO
663 : endif !pe_ix
664 : endif
665 : END DO
666 12 : call timestop("igpt1=1 loop")
667 :
668 : ! (2) igpt1 = 1 , igpt2 = 1 (vanishing G vectors)
669 12 : call timestart("igpt1=igpt2=1 loop")
670 12 : iy = hybdat%nbasp + 1
671 12 : ix = hybdat%nbasp + 1
672 12 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
673 12 : if(pe_ix == fmpi%n_rank) then
674 16 : DO itype1 = 1, fi%atoms%ntype
675 26 : DO ineq1 = 1, fi%atoms%neq(itype1)
676 38 : DO itype2 = 1, fi%atoms%ntype
677 46 : DO ineq2 = 1, fi%atoms%neq(itype2)
678 : coul(1)%data_c(iy, ix_loc) = coul(1)%data_c(iy, ix_loc) &
679 : + rdum*fi%atoms%rmt(itype1)**3*fi%atoms%rmt(itype2)**3* &
680 36 : (fi%atoms%rmt(itype1)**2 + fi%atoms%rmt(itype2)**2)/90
681 : END DO
682 : END DO
683 : END DO
684 : END DO
685 : endif ! pe_ix
686 12 : call timestop("igpt1=igpt2=1 loop")
687 12 : call timestop("add corrections from higher orders")
688 : endif
689 :
690 : ! (3c) r,r' in same MT
691 :
692 : ! Calculate sphbesintegral
693 12 : call timestart("sphbesintegral")
694 0 : allocate (sphbes0(-1:fi%hybinp%lexp + 2, fi%atoms%ntype, nqnrm),&
695 192 : & carr2((fi%hybinp%lexp + 1)**2, maxval(mpdata%n_g)))
696 533572 : sphbes0 = 0; carr2 = 0
697 420 : DO iqnrm = 1, nqnrm
698 1172 : DO itype = 1, fi%atoms%ntype
699 752 : rdum = qnrm(iqnrm)*fi%atoms%rmt(itype)
700 752 : call sphbes(fi%hybinp%lexp + 2, rdum, sphbes0(0, itype, iqnrm))
701 1160 : IF (abs(rdum) > 1e-12) sphbes0(-1, itype, iqnrm) = COS(rdum)/rdum
702 : END DO
703 : END DO
704 12 : call timestop("sphbesintegral")
705 :
706 12 : call timestart("loop 2")
707 48 : DO im = 1, size(fmpi%k_list)
708 36 : ikpt = fmpi%k_list(im)
709 36 : call timestart("harmonics setup")
710 4952 : DO igpt = 1, mpdata%n_g(ikpt)
711 4916 : igptp = mpdata%gptm_ptr(igpt, ikpt)
712 78656 : q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
713 4952 : call ylm4(fi%hybinp%lexp, q, carr2(:, igpt))
714 : END DO
715 36 : call timestop("harmonics setup")
716 : call perform_double_g_loop(fi, hybdat, fmpi, mpdata, sphbes0, carr2, ngptm1,pgptm1,&
717 36 : pqnrm,qnrm, nqnrm, ikpt, coul(ikpt))
718 : ! this one is needed
719 : #ifdef CPP_MPI
720 36 : call timestart("post dblgloop barrier")
721 36 : call MPI_Barrier(fmpi%sub_comm, ierr)
722 36 : call timestop("post dblgloop barrier")
723 : #endif
724 84 : call coul(ikpt)%u2l()
725 : END DO
726 12 : call timestop("loop 2")
727 12 : deallocate (carr2)
728 :
729 : !
730 : ! Symmetry-equivalent G vectors
731 : !
732 : ! All elements are needed so send all data to all processes treating the
733 : ! respective k-points
734 :
735 252 : allocate (carr2(maxval(hybdat%nbasm), 2), iarr(maxval(mpdata%n_g)))
736 : allocate (nsym_gpt(mpdata%num_gpts(), fi%kpts%nkpt), &
737 144 : sym_gpt(MAXVAL(nsym1), mpdata%num_gpts(), fi%kpts%nkpt))
738 258948 : nsym_gpt = 0; sym_gpt = 0
739 12 : call timestart("loop 3")
740 48 : DO im = 1, size(fmpi%k_list)
741 36 : ikpt = fmpi%k_list(im)
742 37044 : carr2 = 0; iarr = 0
743 1468 : iarr(pgptm1(:ngptm1(ikpt), ikpt)) = 1
744 1480 : DO igpt0 = 1, ngptm1(ikpt)
745 1432 : lsym = (1 <= igpt0) .AND. (ngptm1(ikpt) >= igpt0)
746 1432 : igpt2 = pgptm1(igpt0, ikpt)
747 1432 : ix = hybdat%nbasp + igpt2
748 1432 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
749 1432 : if(pe_ix == fmpi%n_rank) then
750 355344 : carr2(:hybdat%nbasm(ikpt),2) = coul(ikpt)%data_c(:hybdat%nbasm(ikpt),ix_loc)
751 : endif
752 : #ifdef CPP_MPI
753 1432 : call timestart("bcast carr2")
754 1432 : call MPI_Bcast(carr2(1,2), hybdat%nbasm(ikpt), MPI_DOUBLE_COMPLEX, pe_ix, fmpi%sub_comm, ierr)
755 1432 : call timestop("bcast carr2")
756 : #endif
757 :
758 1432 : IF (lsym) THEN
759 1432 : ic = 1
760 1432 : sym_gpt(ic, igpt0, ikpt) = igpt2
761 : END IF
762 9504 : DO isym1 = 2, nsym1(ikpt)
763 8072 : isym = sym1(isym1, ikpt)
764 : CALL bramat_trafo(carr2(:, 2), igpt2, ikpt, isym, .FALSE., POINTER(ikpt, :, :, :), &
765 : fi%sym, rrot(:, :, isym), invrrot(:, :, isym), mpdata, fi%hybinp, &
766 : fi%kpts, maxval(fi%hybinp%lcutm1), fi%atoms, fi%hybinp%lcutm1, &
767 : mpdata%num_radbasfn, maxval(mpdata%num_radbasfn), dwgn(:, :, :, isym), &
768 9830032 : hybdat%nbasp, hybdat%nbasm, carr2(:, 1), igpt1)
769 9504 : IF (iarr(igpt1) == 0) THEN
770 : CALL bramat_trafo(carr2(:, 2), igpt2, ikpt, isym, .TRUE., POINTER(ikpt, :, :, :), &
771 : fi%sym, rrot(:, :, isym), invrrot(:, :, isym), mpdata, fi%hybinp, &
772 : fi%kpts, maxval(fi%hybinp%lcutm1), fi%atoms, fi%hybinp%lcutm1, &
773 : mpdata%num_radbasfn, maxval(mpdata%num_radbasfn), &
774 4311088 : dwgn(:, :, :, isym), hybdat%nbasp, hybdat%nbasm, carr2(:, 1), igpt1)
775 3484 : l = (hybdat%nbasp + igpt1 - 1)*(hybdat%nbasp + igpt1)/2
776 3484 : ix = hybdat%nbasp + igpt1
777 3484 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
778 3484 : if(pe_ix == fmpi%n_rank) then
779 710330 : coul(ikpt)%data_c(:hybdat%nbasp + igpt1,ix_loc) = carr2(:hybdat%nbasp + igpt1, 1)
780 : endif
781 :
782 1420660 : do ix = 1,hybdat%nbasp + igpt1
783 1417176 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
784 1420660 : if(pe_ix == fmpi%n_rank) coul(ikpt)%data_c(hybdat%nbasp + igpt1, ix_loc) = conjg(carr2(ix, 1))
785 : enddo
786 :
787 3484 : iarr(igpt1) = 1
788 3484 : IF (lsym) THEN
789 3484 : ic = ic + 1
790 3484 : sym_gpt(ic, igpt0, ikpt) = igpt1
791 : END IF
792 : END IF
793 : END DO
794 1468 : nsym_gpt(igpt0, ikpt) = ic
795 : END DO ! igpt0
796 : END DO ! ikpt
797 12 : call timestop("loop 3")
798 12 : call timestart("gap 1:")
799 12 : deallocate (carr2, iarr, pgptm1)
800 : END IF
801 12 : deallocate (qnrm, pqnrm)
802 :
803 12 : IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN
804 : !
805 : ! The HSE functional is realized subtracting erf/r from
806 : ! the normal Coulomb matrix
807 : !
808 : ELSE
809 : ! check for gamma
810 12 : if(any(fmpi%k_list == 1)) then
811 : CALL subtract_sphaverage(fi%sym, fi%cell, fi%atoms, mpdata, &
812 12 : fi%hybinp, hybdat, fmpi, hybdat%nbasm, gridf, coul(1))
813 : endif
814 : END IF
815 :
816 : ! transform Coulomb matrix to the biorthogonal set
817 12 : call timestop("gap 1:")
818 48 : DO im = 1, size(fmpi%k_list)
819 36 : ikpt = fmpi%k_list(im)
820 36 : call apply_inverse_olaps(mpdata, fi%atoms, fi%cell, hybdat, fmpi, fi%sym, ikpt, coul(ikpt))
821 : ! lower to upper, because the lower half is better in memory
822 : #ifdef CPP_MPI
823 36 : call timestart("post inverse barrier")
824 36 : call MPI_BARRIER(fmpi%sub_comm, ierr)
825 36 : call timestop("post inverse barrier")
826 : #endif
827 84 : call coul(ikpt)%l2u()
828 : enddo
829 :
830 : !call plot_coulombmatrix() -> code was shifted to plot_coulombmatrix.F90
831 : !
832 : ! rearrange coulomb matrix
833 : !
834 12 : if(.not. allocated(hybdat%coul)) allocate(hybdat%coul(fi%kpts%nkpt))
835 :
836 48 : DO im = 1, size(fmpi%k_list)
837 36 : ikpt = fmpi%k_list(im)
838 : ! unpack coulomb into coulomb(ikpt)
839 48 : call copy_from_dense_to_sparse(fi, fmpi, mpdata, coul, ikpt, hybdat)
840 : END DO ! ikpt
841 12 : CALL timestop("Coulomb matrix setup")
842 :
843 524 : END SUBROUTINE coulombmatrix
844 :
845 : ! Calculate body of Coulomb matrix at Gamma point: v_IJ = SUM(G) c^*_IG c_JG 4*pi/G**2 .
846 : ! For this we must subtract from coulomb(:,1) the spherical average of a term that comes
847 : ! from the fact that MT functions have k-dependent Fourier coefficients (see script).
848 12 : SUBROUTINE subtract_sphaverage(sym, cell, atoms, mpdata, hybinp, hybdat, fmpi, nbasm1, gridf, coulomb)
849 :
850 : USE m_types
851 : USE m_constants
852 : USE m_wrapper
853 : USE m_trafo
854 : USE m_util
855 : use m_intgrf
856 : USE m_olap
857 : IMPLICIT NONE
858 :
859 : TYPE(t_sym), INTENT(IN) :: sym
860 : TYPE(t_cell), INTENT(IN) :: cell
861 : TYPE(t_atoms), INTENT(IN) :: atoms
862 : TYPE(t_mpdata), intent(in) :: mpdata
863 : TYPE(t_hybinp), INTENT(IN) :: hybinp
864 : TYPE(t_hybdat), INTENT(IN) :: hybdat
865 : type(t_mpi), intent(in) :: fmpi
866 :
867 : INTEGER, INTENT(IN) :: nbasm1(:)
868 : REAL, INTENT(IN) :: gridf(:, :)
869 : class(t_mat), intent(inout) :: coulomb
870 :
871 : ! - local scalars -
872 : INTEGER :: l, ix, iy, ix_loc, pe_ix, i, j, n, nn, itype, ieq, M, ierr
873 :
874 : ! - local arrays -
875 12 : TYPE(t_mat) :: olap
876 : !COMPLEX , ALLOCATABLE :: constfunc(:) !can also be real in inversion case
877 12 : COMPLEX :: coeff(1,nbasm1(1)), cderiv(-1:1, nbasm1(1)), claplace(1,nbasm1(1))
878 :
879 12 : call timestart("subtract_sphaverage")
880 12 : CALL olap%alloc(sym%invs, mpdata%n_g(1), mpdata%n_g(1), 0.)
881 :
882 12 : n = nbasm1(1)
883 12 : nn = n*(n + 1)/2
884 8124 : CALL olap_pw(olap, mpdata%g(:, mpdata%gptm_ptr(:mpdata%n_g(1), 1)), mpdata%n_g(1), atoms, cell, fmpi)
885 :
886 : ! Define coefficients (coeff) and their derivatives (cderiv,claplace)
887 10212 : coeff = 0
888 20412 : cderiv = 0
889 10212 : claplace = 0
890 12 : j = 0
891 32 : DO itype = 1, atoms%ntype
892 52 : DO ieq = 1, atoms%neq(itype)
893 140 : DO l = 0, hybinp%lcutm1(itype)
894 620 : DO M = -l, l
895 4080 : DO i = 1, mpdata%num_radbasfn(l, itype)
896 3480 : j = j + 1
897 3980 : IF (l == 0) THEN
898 : coeff(1,j) = SQRT(fpi_const) &
899 : *intgrf(atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, i, 0, itype), &
900 : atoms, itype, gridf) &
901 146492 : /SQRT(cell%vol)
902 :
903 : claplace(1,j) = -SQRT(fpi_const) &
904 : *intgrf(atoms%rmsh(:, itype)**3*mpdata%radbasfn_mt(:, i, 0, itype), &
905 : atoms, itype, gridf) &
906 146492 : /SQRT(cell%vol)
907 :
908 3300 : ELSE IF (l == 1) THEN
909 : cderiv(M,j) = -SQRT(fpi_const/3)*CMPLX(0.0, 1.0) &
910 : *intgrf(atoms%rmsh(:, itype)**2*mpdata%radbasfn_mt(:, i, 1, itype), &
911 : atoms, itype, gridf) &
912 408240 : /SQRT(cell%vol)
913 : END IF
914 : END DO
915 : END DO
916 : END DO
917 : END DO
918 : END DO
919 12 : IF (olap%l_real) THEN
920 952 : coeff(1,hybdat%nbasp + 1:n) = olap%data_r(1, 1:n - hybdat%nbasp)
921 : else
922 680 : coeff(1,hybdat%nbasp + 1:n) = olap%data_c(1, 1:n - hybdat%nbasp)
923 : END IF
924 12 : IF (sym%invs) THEN
925 : CALL symmetrize(coeff, 1, nbasm1(1), 2, &
926 : atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
927 20 : mpdata%num_radbasfn, sym)
928 : CALL symmetrize(claplace, 1, nbasm1(1), 2, &
929 : atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
930 20 : mpdata%num_radbasfn, sym)
931 : CALL symmetrize(cderiv(-1:-1,:), 1, nbasm1(1), 2, &
932 : atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
933 20 : mpdata%num_radbasfn, sym)
934 : CALL symmetrize(cderiv(0:0,:), 1, nbasm1(1), 2, &
935 : atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
936 20 : mpdata%num_radbasfn, sym)
937 : CALL symmetrize(cderiv(1:1,:), 1, nbasm1(1), 2, &
938 : atoms, hybinp%lcutm1, maxval(hybinp%lcutm1), &
939 20 : mpdata%num_radbasfn, sym)
940 : ENDIF
941 : ! Subtract head contributions from coulomb(:nn,1) to obtain the body
942 12 : l = 0
943 5112 : DO ix = 1, n
944 5100 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
945 5112 : if(fmpi%n_rank == pe_ix) then
946 597218 : DO iy = 1, ix
947 594668 : l = l + 1
948 : coulomb%data_c(iy,ix_loc) = coulomb%data_c(iy,ix_loc) - fpi_const/3 &
949 : *(dot_PRODUCT(cderiv(:,iy), cderiv(:,ix)) &
950 : + (CONJG(coeff(1,iy))*claplace(1,ix) &
951 2381222 : + CONJG(claplace(1,iy))*coeff(1,ix))/2)
952 : END DO
953 : endif
954 : END DO
955 :
956 : !needed bc apply inverse uses lower half
957 : #ifdef CPP_MPI
958 12 : call timestart("post subtr avg barrier")
959 12 : call MPI_Barrier(fmpi%sub_comm, ierr)
960 12 : call timestop("post subtr avg barrier")
961 : #endif
962 12 : call coulomb%u2l()
963 12 : call timestop("subtract_sphaverage")
964 12 : END SUBROUTINE subtract_sphaverage
965 :
966 :
967 :
968 : ! ---------
969 :
970 : ! Returns a list of (k+G) vector lengths in qnrm(1:nqnrm) and the corresponding pointer pqnrm(1:ngpt(ikpt),ikpt)
971 12 : SUBROUTINE getnorm(kpts, gpt, ngpt, pgpt, qnrm, nqnrm, pqnrm, cell)
972 : USE m_types
973 : USE m_juDFT
974 : IMPLICIT NONE
975 : TYPE(t_cell), INTENT(IN) :: cell
976 : TYPE(t_kpts), INTENT(IN) :: kpts
977 :
978 : INTEGER, INTENT(IN) :: ngpt(:), gpt(:, :), pgpt(:, :)!(dim,kpts%nkpt)
979 12 : REAL, ALLOCATABLE :: qnrm(:), help(:)
980 : INTEGER, INTENT(INOUT) :: nqnrm
981 : INTEGER, ALLOCATABLE :: pqnrm(:, :)
982 : INTEGER :: i, j, ikpt, igpt, igptp
983 : REAL :: q(3), qnorm
984 :
985 5484 : allocate (qnrm(MAXVAL(ngpt)*kpts%nkpt), source=0.0)
986 5532 : allocate (pqnrm(MAXVAL(ngpt), kpts%nkpt), source=0)
987 : i = 0
988 48 : DO ikpt = 1, kpts%nkpt
989 4964 : igptloop: DO igpt = 1, ngpt(ikpt)
990 4916 : igptp = pgpt(igpt, ikpt)
991 4916 : IF (igptp == 0) call judft_error('getnorm: zero pointer (bug?)')
992 78656 : q = MATMUL(kpts%bk(:, ikpt) + gpt(:, igptp), cell%bmat)
993 19664 : qnorm = norm2(q)
994 104060 : DO j = 1, i
995 104060 : IF (ABS(qnrm(j) - qnorm) < 1e-12) THEN
996 4508 : pqnrm(igpt, ikpt) = j
997 4508 : CYCLE igptloop
998 : END IF
999 : END DO
1000 408 : i = i + 1
1001 408 : qnrm(i) = qnorm
1002 444 : pqnrm(igpt, ikpt) = i
1003 : END DO igptloop
1004 : END DO
1005 12 : nqnrm = i
1006 :
1007 36 : allocate (help(nqnrm))
1008 420 : help(1:nqnrm) = qnrm(1:nqnrm)
1009 12 : deallocate (qnrm)
1010 36 : allocate (qnrm(1:nqnrm))
1011 432 : qnrm = help
1012 :
1013 12 : END SUBROUTINE getnorm
1014 :
1015 324 : subroutine loop_over_interst(fi, hybdat, mpdata, fmpi, structconst, sphbesmoment, moment, moment2, &
1016 36 : qnrm, facc, gmat, integral, olap, pqnrm, pgptm1, ngptm1, ikpt, coul)
1017 : use m_types
1018 : use m_juDFT
1019 : use m_ylm, only: ylm4
1020 : use m_constants, only: fpi_const, tpi_const
1021 : USE m_trafo, ONLY: symmetrize
1022 : use m_calc_l_m_from_lm
1023 : implicit none
1024 :
1025 : type(t_fleurinput), intent(in) :: fi
1026 : type(t_hybdat), intent(in) :: hybdat
1027 : type(t_mpdata), intent(in) :: mpdata
1028 : type(t_mpi), intent(in) :: fmpi
1029 : REAL, intent(in) :: sphbesmoment(0:, :, :), qnrm(:), facC(-1:), gmat(:, :), moment(:, 0:, :), moment2(:, :)
1030 : real, intent(in) :: integral(:, 0:, :, :), olap(:, 0:, :, :)
1031 : integer, intent(in) :: ikpt, ngptm1(:), pqnrm(:, :), pgptm1(:, :)
1032 : complex, intent(in) :: structconst(:, :, :, :)
1033 : class(t_mat), intent(inout) :: coul
1034 :
1035 : integer :: igpt0, igpt, igptp, iqnrm, niter
1036 : integer :: ix, iy, ic, itype, lm, l, m, itype1, ic1, l1, m1, lm1, loc_from
1037 : integer :: l2, m2, lm2, n, i, iatm, j_type, j_l, iy_start, j_m, j_lm, pe_ix, ix_loc
1038 : real :: q(3), qnorm, svol, tmp_vec(3)
1039 36 : COMPLEX :: y((fi%hybinp%lexp + 1)**2), y1((fi%hybinp%lexp + 1)**2), y2((fi%hybinp%lexp + 1)**2)
1040 : complex :: csum, csumf(9), cdum, cexp
1041 36 : integer, allocatable :: lm_arr(:), ic_arr(:)
1042 :
1043 :
1044 36 : call range_from_glob_to_loc(fmpi, hybdat%nbasp+1, loc_from)
1045 795666 : coul%data_c(:hybdat%nbasp,loc_from:) = 0
1046 :
1047 36 : svol = SQRT(fi%cell%vol)
1048 : ! start to loop over interstitial plane waves
1049 : !DO igpt0 = 1, ngptm1(ikpt)
1050 1468 : do igpt0 = 1, ngptm1(ikpt)
1051 1432 : igpt = pgptm1(igpt0, ikpt)
1052 1432 : igptp = mpdata%gptm_ptr(igpt, ikpt)
1053 1432 : ix = hybdat%nbasp + igpt
1054 1432 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
1055 1468 : if(pe_ix == fmpi%n_rank) then
1056 11456 : q = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%cell%bmat)
1057 2864 : qnorm = norm2(q)
1058 716 : iqnrm = pqnrm(igpt, ikpt)
1059 716 : IF (ABS(qnrm(iqnrm) - qnorm) > 1e-12) then
1060 0 : call judft_error('coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)') ! We shouldn't st op here!
1061 : endif
1062 :
1063 9308 : tmp_vec = MATMUL(fi%kpts%bk(:, fi%kpts%nkpt), fi%cell%bmat)
1064 716 : call ylm4(2, tmp_vec, y1)
1065 11456 : tmp_vec = MATMUL(mpdata%g(:, igptp), fi%cell%bmat)
1066 716 : call ylm4(2, tmp_vec, y2)
1067 716 : call ylm4(fi%hybinp%lexp, q, y)
1068 621488 : y1 = CONJG(y1); y2 = CONJG(y2); y = CONJG(y)
1069 :
1070 : ! this unrolls the do ic=1,atoms%nat{do lm=1,..{}}
1071 716 : call collapse_ic_and_lm_loop(fi%atoms, fi%hybinp%lcutm1, niter, ic_arr, lm_arr)
1072 :
1073 : !$OMP PARALLEL DO default(none) &
1074 : !$OMP private(ic, lm, itype, l, m, csum, csumf, ic1, itype1, cexp, lm1, l2, cdum, m2, lm2, iy) &
1075 : !$OMP private(j_m, j_type, iy_start, l1, m1) &
1076 : !$OMP shared(ic_arr, lm_arr, fi, mpdata, olap, qnorm, moment, integral, hybdat, svol) &
1077 : !$OMP shared(moment2, ix, igpt, facc, structconst, y, y1, y2, gmat, iqnrm, sphbesmoment, ikpt) &
1078 : !$OMP shared(igptp, niter, fmpi, pe_ix, coul, ix_loc) &
1079 716 : !$OMP schedule(dynamic)
1080 : do i = 1,niter
1081 : ic = ic_arr(i)
1082 : lm = lm_arr(i)
1083 :
1084 : itype = fi%atoms%itype(ic)
1085 : call calc_l_m_from_lm(lm, l, m)
1086 :
1087 : ! calculate sum over lm and centers for (2c) -> csum, csumf
1088 : csum = 0
1089 : csumf = 0
1090 : do ic1 = 1, fi%atoms%nat
1091 : itype1 = fi%atoms%itype(ic1)
1092 : cexp = fpi_const*EXP(CMPLX(0.0, 1.0)*tpi_const &
1093 : *(dot_PRODUCT(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp), fi%atoms%taual(:, ic1)) &
1094 : - dot_PRODUCT(fi%kpts%bk(:, ikpt), fi%atoms%taual(:, ic))))
1095 :
1096 : do lm1 = 1, (fi%hybinp%lexp+1)**2
1097 : call calc_l_m_from_lm(lm1, l1, m1)
1098 : l2 = l + l1 ! for structconst
1099 : cdum = sphbesmoment(l1, itype1, iqnrm)*CMPLX(0.0, 1.0)**(l1)*cexp
1100 : m2 = M - m1 ! for structconst
1101 : lm2 = l2**2 + l2 + m2 + 1 !
1102 : csum = csum - (-1)**(m1 + l1)*gmat(lm1, lm)*y(lm1)*cdum*structconst(lm2, ic, ic1, ikpt)
1103 : END DO
1104 :
1105 : ! add contribution of (2c) to csum and csumf coming from linear and quadratic orders of Y_lm*(G) / G * j_(l+1)(GS)
1106 : IF (ikpt == 1 .AND. l <= 2) THEN
1107 : cexp = EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic1))) &
1108 : *gmat(lm, 1)*fpi_const/fi%cell%vol
1109 : csumf(lm) = csumf(lm) - cexp*SQRT(fpi_const)* &
1110 : CMPLX(0.0, 1.0)**l*sphbesmoment(0, itype1, iqnrm)/facC(l - 1)
1111 : IF (l == 0) THEN
1112 : IF (igpt /= 1) THEN
1113 : csum = csum - cexp*(sphbesmoment(0, itype1, iqnrm)*fi%atoms%rmt(itype1)**2 - &
1114 : sphbesmoment(2, itype1, iqnrm)*2.0/3)/10
1115 : ELSE
1116 : csum = csum - cexp*fi%atoms%rmt(itype1)**5/30
1117 : END IF
1118 : ELSE IF (l == 1) THEN
1119 : csum = csum + cexp*CMPLX(0.0, 1.0)*SQRT(fpi_const) &
1120 : *sphbesmoment(1, itype1, iqnrm)*y(lm)/3
1121 : END IF
1122 : END IF
1123 : END DO
1124 :
1125 : ! add contribution of (2a) to csumf
1126 : IF (ikpt == 1 .AND. igpt == 1 .AND. l <= 2) THEN
1127 : csumf(lm) = csumf(lm) + (fpi_const)**2*CMPLX(0.0, 1.0)**l/facC(l)
1128 : END IF
1129 :
1130 : ! finally define coulomb
1131 : cdum = (fpi_const)**2*CMPLX(0.0, 1.0)**(l)*y(lm) &
1132 : *EXP(CMPLX(0.0, 1.0)*tpi_const &
1133 : *dot_PRODUCT(mpdata%g(:, igptp), fi%atoms%taual(:, ic)))
1134 :
1135 : !calclate iy_start on the fly for OpenMP
1136 : iy_start = 0
1137 : do iatm = 1, ic-1
1138 : j_type = fi%atoms%itype(iatm)
1139 : do j_l = 0,fi%hybinp%lcutm1(j_type)
1140 : iy_start = iy_start + mpdata%num_radbasfn(j_l, j_type) * (2*j_l+1)
1141 : end do
1142 : end do
1143 : do j_lm = 1,lm-1
1144 : call calc_l_m_from_lm(j_lm, j_l, j_m)
1145 : iy_start = iy_start + mpdata%num_radbasfn(j_l, itype)
1146 : enddo
1147 :
1148 : DO n = 1, mpdata%num_radbasfn(l, itype)
1149 : iy = iy_start + n
1150 :
1151 : IF (ikpt == 1 .AND. igpt == 1) THEN
1152 : IF (l == 0) coul%data_c(iy, ix_loc) = -cdum*moment2(n, itype)/6/svol
1153 : coul%data_c(iy, ix_loc) = coul%data_c(iy, ix_loc) &
1154 : + (-cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
1155 : + csum*moment(n, l, itype))/svol ! (2c)
1156 : ELSE
1157 : coul%data_c(iy, ix_loc) = (cdum*olap(n, l, itype, iqnrm)/qnorm**2 & ! (2a)&
1158 : - cdum/(2*l + 1)*integral(n, l, itype, iqnrm) & ! (2b)&
1159 : + csum*moment(n, l, itype))/svol ! (2c)
1160 : endif
1161 : END DO
1162 : END DO ! collapsed atom & lm loop (ic)
1163 : !$OMP END PARALLEL DO
1164 : endif !pe_ix
1165 : END DO
1166 :
1167 36 : IF (fi%sym%invs) THEN
1168 : call symmetrize_mpimat(fi, fmpi, coul%data_c, [1,hybdat%nbasp+1], [hybdat%nbasp, hybdat%nbasp+mpdata%n_g(ikpt)],&
1169 120 : 1, .false., mpdata%num_radbasfn)
1170 : ENDIF
1171 36 : endsubroutine loop_over_interst
1172 :
1173 36 : subroutine perform_double_g_loop(fi, hybdat, fmpi, mpdata, sphbes0, carr2, ngptm1,pgptm1,pqnrm,qnrm, nqnrm, ikpt, coulomb)
1174 : use m_juDFT
1175 : use m_constants, only: tpi_const,fpi_const
1176 : use m_sphbessel_integral
1177 : implicit none
1178 : type(t_fleurinput), intent(in) :: fi
1179 : TYPE(t_mpdata), intent(in) :: mpdata
1180 : TYPE(t_hybdat), INTENT(IN) :: hybdat
1181 : type(t_mpi), intent(in) :: fmpi
1182 : integer, intent(in) :: ikpt, ngptm1(:), pqnrm(:,:),pgptm1(:, :), nqnrm
1183 : real, intent(in) :: qnrm(:), sphbes0(:, :, :)
1184 : complex, intent(in) :: carr2(:, :)
1185 : !complex, intent(inout) :: coulomb(:) ! only at ikpt
1186 : class(t_mat), intent(inout) :: coulomb
1187 :
1188 : integer :: igpt0, igpt1, igpt2, ix, iy, igptp1, igptp2, iqnrm1, iqnrm2
1189 : integer :: ic, itype, lm, m, l, pe_ix, ix_loc
1190 : real :: q1(3), q2(3)
1191 36 : complex :: y1((fi%hybinp%lexp + 1)**2), y2((fi%hybinp%lexp + 1)**2)
1192 36 : COMPLEX :: cexp1(fi%atoms%ntype)
1193 : complex :: cdum, cdum1
1194 : logical :: ldum
1195 :
1196 36 : call timestart("double g-loop")
1197 :
1198 1468 : DO igpt0 = 1, ngptm1(ikpt)
1199 1432 : igpt2 = pgptm1(igpt0, ikpt)
1200 1432 : ix = hybdat%nbasp + igpt2
1201 1432 : call glob_to_loc(fmpi, ix, pe_ix, ix_loc)
1202 1468 : if(pe_ix == fmpi%n_rank) then
1203 716 : igptp2 = mpdata%gptm_ptr(igpt2, ikpt)
1204 716 : iqnrm2 = pqnrm(igpt2, ikpt)
1205 11456 : q2 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp2), fi%cell%bmat)
1206 207640 : y2 = CONJG(carr2(:, igpt2))
1207 :
1208 : !$OMP PARALLEL DO default(none) &
1209 : !$OMP private(igpt1, iy, igptp1, iqnrm1, q1, y1, cexp1, ic, itype, lm) &
1210 : !$OMP private(cdum, l, cdum1, m, ldum) &
1211 : !$OMP shared(igpt2, coulomb, hybdat, mpdata, ikpt, fi, carr2, pqnrm, igptp2)&
1212 716 : !$OMP shared(qnrm, sphbes0, iqnrm2, nqnrm, y2, ix_loc)
1213 : DO igpt1 = 1, igpt2
1214 : iy = hybdat%nbasp + igpt1
1215 : igptp1 = mpdata%gptm_ptr(igpt1, ikpt)
1216 : iqnrm1 = pqnrm(igpt1, ikpt)
1217 : q1 = MATMUL(fi%kpts%bk(:, ikpt) + mpdata%g(:, igptp1), fi%cell%bmat)
1218 : y1 = carr2(:, igpt1)
1219 : cexp1 = 0
1220 : do ic = 1,fi%atoms%nat
1221 : itype = fi%atoms%itype(ic)
1222 : cexp1(itype) = cexp1(itype) + &
1223 : EXP(CMPLX(0.0, 1.0)*tpi_const*dot_PRODUCT( &
1224 : (mpdata%g(:, igptp2) - mpdata%g(:, igptp1)), fi%atoms%taual(:, ic)))
1225 : ENDDO
1226 : lm = 0
1227 : cdum = 0
1228 : DO l = 0, fi%hybinp%lexp
1229 : cdum1 = 0
1230 : DO itype = 1, fi%atoms%ntype
1231 : cdum1 = cdum1 + cexp1(itype)*sphbessel_integral( &
1232 : fi%atoms, itype, qnrm, nqnrm, &
1233 : iqnrm1, iqnrm2, l, fi%hybinp, &
1234 : sphbes0, .False., ldum) &
1235 : /(2*l + 1)
1236 : END DO
1237 : DO M = -l, l
1238 : lm = lm + 1
1239 : cdum = cdum + cdum1*y1(lm)*y2(lm)
1240 : ENDDO
1241 : ENDDO
1242 : coulomb%data_c(iy,ix_loc) = coulomb%data_c(iy,ix_loc) + (fpi_const)**3*cdum/fi%cell%vol
1243 : END DO
1244 : !$OMP end parallel do
1245 : endif !pe_ix
1246 : END DO
1247 36 : call timestop("double g-loop")
1248 36 : end subroutine perform_double_g_loop
1249 :
1250 716 : subroutine collapse_ic_and_lm_loop(atoms, lcutm1, niter, ic_arr, lm_arr)
1251 : use m_types
1252 : implicit none
1253 : type(t_atoms), intent(in) :: atoms
1254 : integer, intent(in) :: lcutm1(:)
1255 : integer, intent(out) :: niter
1256 : integer, intent(inout), allocatable :: ic_arr(:), lm_arr(:)
1257 :
1258 : integer :: ic, lm, itype
1259 :
1260 716 : if(allocated(ic_arr)) deallocate(ic_arr)
1261 716 : if(allocated(lm_arr)) deallocate(lm_arr)
1262 :
1263 716 : niter = 0
1264 2082 : do ic = 1, atoms%nat
1265 1366 : itype = atoms%itype(ic)
1266 36232 : do lm = 1,(lcutm1(itype)+1)**2
1267 35516 : niter = niter + 1
1268 : enddo
1269 : enddo
1270 :
1271 2864 : allocate( lm_arr(niter), ic_arr(niter))
1272 716 : niter = 0
1273 2082 : do ic = 1, atoms%nat
1274 1366 : itype = atoms%itype(ic)
1275 36232 : do lm = 1,(lcutm1(itype)+1)**2
1276 34150 : niter = niter + 1
1277 34150 : lm_arr(niter) = lm
1278 35516 : ic_arr(niter) = ic
1279 : enddo
1280 : enddo
1281 716 : end subroutine collapse_ic_and_lm_loop
1282 :
1283 12 : subroutine bessel_calculation(fi, fmpi, mpdata, nqnrm, gridf, qnrm, sphbesmoment, olap, integral)
1284 : implicit NONE
1285 : type(t_fleurinput), intent(in) :: fi
1286 : type(t_mpi), intent(in) :: fmpi
1287 : type(t_mpdata), intent(in) :: mpdata
1288 : integer, intent(in) :: nqnrm
1289 : real, intent(in) :: gridf(:,:), qnrm(:)
1290 : real, intent(inout) :: sphbesmoment(0:,:,:), olap(:,0:,:,:), integral(:,0:,:,:)
1291 :
1292 : integer :: iqnrm, itype, i, l, n, ng, buf_sz, root, ierr
1293 32 : REAL :: sphbes_var(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
1294 32 : REAL :: sphbesmoment1(fi%atoms%jmtd, 0:maxval(fi%hybinp%lcutm1))
1295 32 : REAL :: rarr(0:fi%hybinp%lexp + 1), rarr1(0:maxval(fi%hybinp%lcutm1))
1296 : real :: rdum, qnorm, rdum1
1297 :
1298 12 : call timestart("Bessel calculation")
1299 :
1300 12 : do iqnrm = 1+fmpi%irank, nqnrm, fmpi%isize
1301 204 : qnorm = qnrm(iqnrm)
1302 : !$OMP parallel do default(none) &
1303 : !$OMP shared(olap, integral, sphbesmoment, fi,qnorm, iqnrm, mpdata, gridf) &
1304 204 : !$OMP private(itype, rdum, sphbes_var, sphbesmoment1, ng, rarr, rarr1, rdum1, i, l)
1305 : DO itype = 1, fi%atoms%ntype
1306 : ng = fi%atoms%jri(itype)
1307 : rdum = fi%atoms%rmt(itype)
1308 : sphbes_var = 0
1309 : sphbesmoment1 = 0
1310 : IF (abs(qnorm) < 1e-12) THEN
1311 : sphbesmoment(0, itype, iqnrm) = rdum**3/3
1312 : DO i = 1, ng
1313 : sphbes_var(i, 0) = 1
1314 : sphbesmoment1(i, 0) = fi%atoms%rmsh(i, itype)**2/3 + (rdum**2 - fi%atoms%rmsh(i, itype)**2)/2
1315 : END DO
1316 : ELSE
1317 : call sphbes(fi%hybinp%lexp + 1, qnorm*rdum, rarr)
1318 : DO l = 0, fi%hybinp%lexp
1319 : sphbesmoment(l, itype, iqnrm) = rdum**(l + 2)*rarr(l + 1)/qnorm
1320 : END DO
1321 : DO i = ng, 1, -1
1322 : rdum = fi%atoms%rmsh(i, itype)
1323 : call sphbes(fi%hybinp%lcutm1(itype) + 1, qnorm*rdum, rarr)
1324 : DO l = 0, fi%hybinp%lcutm1(itype)
1325 : sphbes_var(i, l) = rarr(l)
1326 : IF (l /= 0) THEN; rdum1 = -rdum**(1 - l)*rarr(l - 1)
1327 : ELSE; rdum1 = -COS(qnorm*rdum)/qnorm
1328 : ENDIF
1329 : IF (i == ng) rarr1(l) = rdum1
1330 : sphbesmoment1(i, l) = (rdum**(l + 2)*rarr(l + 1)/rdum**(l + 1) &
1331 : + (rarr1(l) - rdum1)*rdum**l)/qnorm
1332 : END DO
1333 : END DO
1334 : END IF
1335 : DO l = 0, fi%hybinp%lcutm1(itype)
1336 : DO n = 1, mpdata%num_radbasfn(l, itype)
1337 : ! note that mpdata%radbasfn_mt already contains one factor rgrid
1338 : olap(n, l, itype, iqnrm) = &
1339 : intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbes_var(:, l), &
1340 : fi%atoms, itype, gridf)
1341 :
1342 : integral(n, l, itype, iqnrm) = &
1343 : intgrf(fi%atoms%rmsh(:, itype)*mpdata%radbasfn_mt(:, n, l, itype)*sphbesmoment1(:, l), &
1344 : fi%atoms, itype, gridf)
1345 :
1346 : END DO
1347 : END DO
1348 : END DO
1349 : END DO
1350 :
1351 : #ifdef CPP_MPI
1352 12 : call timestart("bcast bessel")
1353 420 : do iqnrm = 1, nqnrm
1354 408 : root = mod(iqnrm - 1,fmpi%isize)
1355 408 : buf_sz = size(olap,1) * size(olap,2) * size(olap,3)
1356 408 : call MPI_Bcast(olap(1,0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
1357 :
1358 408 : buf_sz = size(integral,1) * size(integral,2) * size(integral,3)
1359 408 : call MPI_Bcast(integral(1,0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
1360 :
1361 408 : buf_sz = size(sphbesmoment,1) * size(sphbesmoment,2)
1362 420 : call MPI_Bcast(sphbesmoment(0,1,iqnrm), buf_sz, MPI_DOUBLE_PRECISION, root, fmpi%mpi_comm, ierr)
1363 : enddo
1364 12 : call timestop("bcast bessel")
1365 : #endif
1366 12 : call timestop("Bessel calculation")
1367 12 : end subroutine bessel_calculation
1368 :
1369 12 : function calc_num_mtmts(fi) result(num_mtmt)
1370 : implicit none
1371 : type(t_fleurinput), intent(in) :: fi
1372 : integer :: num_mtmt
1373 :
1374 : integer :: itype, ineq, l, m
1375 :
1376 12 : num_mtmt = 0
1377 32 : DO itype = 1, fi%atoms%ntype
1378 52 : DO ineq = 1, fi%atoms%neq(itype)
1379 140 : DO l = 0, fi%hybinp%lcutm1(itype)
1380 620 : DO M = -l, l
1381 600 : num_mtmt = num_mtmt + 1
1382 : enddo
1383 : enddo
1384 : enddo
1385 : enddo
1386 12 : end function calc_num_mtmts
1387 :
1388 : END MODULE m_coulombmatrix
|