Line data Source code
1 : !!****m* ABINIT/m_paw_slater
2 : !! NAME
3 : !! m_paw_slater
4 : !!
5 : !! FUNCTION
6 : !! This module defines objects and procedures to evaluate Slater-like integrals
7 : !! using real spherical Harmonics.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! NOTES
16 : !! * Routines tagged with "@type_name" are tightly connected to the definition of the data type.
17 : !! Tightly connected means that the proper functioning of the implementation relies on the
18 : !! assumption that the tagged procedure is consistent with the type declaration.
19 : !! Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
20 : !! that all the tightly connected routines are changed accordingly to accommodate the modification of the data type.
21 : !! Typical examples of tightly connected routines are creation, destruction or reset methods.
22 : !!
23 : !! SOURCE
24 :
25 : #if defined HAVE_CONFIG_H
26 : #include "config.h"
27 : #endif
28 :
29 : #include "abi_common.h"
30 :
31 : MODULE m_paw_slater
32 :
33 : use defs_basis
34 : use m_abicore
35 : use m_errors
36 : use m_splines
37 :
38 : use m_fstrings, only : basename
39 : use m_paw_atomorb, only : atomorb_type, print_atomorb, destroy_atomorb, get_overlap
40 : use m_pawpsp, only : pawpsp_init_core
41 : use m_crystal, only : crystal_t
42 : use m_paw_io, only : pawio_print_ij
43 : use m_pawang, only : pawang_type
44 : use m_paw_sphharm, only : realgaunt
45 : use m_pawrad, only : pawrad_type, pawrad_free, pawrad_isame, &
46 : & pawrad_deducer0, simp_gen, calc_slatradl,pawrad_copy
47 : use m_pawtab, only : pawtab_type
48 : use m_pawrhoij, only : pawrhoij_type
49 : use m_paw_lmn, only : make_kln2ln, make_klm2lm, make_indln, klmn2ijlmn
50 :
51 : implicit none
52 :
53 : private
54 :
55 : public :: paw_sigxcore ! The onsite matrix elements of the Fock operator generated by (closed) core shells.
56 : public :: paw_mkdijexc_core ! Calculate the onsite matrix element of the Fock operator generated by the core.
57 : public :: paw_dijhf ! Compute the onsite D_{ij} strengths of the exchange parth of the self energy.
58 : !!***
59 :
60 : !!****t* m_paw_slater/slatrad_t
61 : !! NAME
62 : !! slatrad_t
63 : !!
64 : !! FUNCTION
65 : !! Object used to store radial integrals of the form.
66 : !!
67 : !! $ F_{ijkl}^L = \dfrac{4\pi}{2L+1} \int u_i(1) u_j(1) \dfrac{r_<^L}{r_>^{L+1}} u_k(2) u_l(2) d1d2 $
68 : !!
69 : !! for a given quadruple (i,j,k,l) as a function L \in [L_min, L_max].
70 : !! i,j,k,l are shorthand indices for (nn,ll) quantum numbers.
71 : !!
72 : !! NOTES
73 : !! Basic symmetry properties:
74 : !! 1) invariant under the exchange (i<-->j) and (k<-->l).
75 : !! 2) invariant under the exchange (i,j) <--> (k,l).
76 : !!
77 : !! Memory saving is achieved by storing the upper triangle of the (ij) (kl) matrix
78 : !! and, for each dimension, only the upper triangle of the two matrices (iln,jln) (kln,lln).
79 : !!
80 : !! Some matrix elements will never contribute to <ij|1/|1-2||kl> due to selection rules
81 : !! introduced by the integration of the angular part.
82 : !!
83 : !! SOURCE
84 :
85 : type,public :: slatrad_t
86 :
87 : integer :: iln,jln,kln,lln
88 : ! The (l,n) indices associated to the partial waves.
89 :
90 : integer :: lslat_min
91 : ! Min l+1 in the expansion of the Coulomb potential.
92 :
93 : integer :: lslat_max
94 : ! Max l+1 in the expansion of the Coulomb potential.
95 :
96 : integer :: nintgl
97 : ! The number of non-zero integrals stored in intgl.
98 :
99 : integer,allocatable :: intgl_select(:)
100 : ! intgl_select(lslat_min:lslat_max)
101 : ! Index of the non-zero integrals in intgl, 0 if intgl has not been
102 : ! calculated thanks to selection rules coming from the angular integration.
103 :
104 : real(dp),allocatable :: intgl(:)
105 : ! intgl(1:nintgl)
106 : ! The integrals:
107 : ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phi_{\nj\lj}(1) \dfrac{r_<^L}{r_>^{L+1}} \phi_{\nk\lk}(2) \phi_{\nl\ll} d1d2
108 : ! for given (i,j,k,l) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
109 :
110 : end type slatrad_t
111 :
112 : public :: slatrad_init ! Creation method
113 : public :: slatrad_free ! Free memory
114 : !!***
115 :
116 : interface slatrad_free
117 : module procedure slatrad_free_0D
118 : module procedure slatrad_free_1D
119 : end interface slatrad_free
120 :
121 : !----------------------------------------------------------------------
122 :
123 : !!****t* m_paw_slater/slatang_cshell_t
124 : !! NAME
125 : !! slatang_cshell_t
126 : !!
127 : !! FUNCTION
128 : !! Object used to store:
129 : !! $ F^{lsl,lc}_{li,lj,mi,mj} = sum_{msl mc} \<li mi|lsl msl;lc mc\> \<lsl msl;lc mc| lj mj\> $
130 : !! This (less general) type of radial integral is needed to evaluate the Exchange term generated
131 : !! by a closed-shell atom. In the equation, (lc,mc) are the set of angular quantum number associated
132 : !! to (closed) core electrons while (lsl,msl) comes from the expansion of 1/|r1-r2|.
133 : !! Since the F is invariant under exchange of i and j we use an array of structures indexed
134 : !! by kln = (iln,jln) in packed form.
135 : !!
136 : !! SOURCE
137 :
138 : type, public :: slatang_cshell_t
139 :
140 : integer :: nsggsel
141 : ! Number of non null matrix elements
142 :
143 : integer :: lslat_max
144 : ! Max l+1 in the expansion of the Coulomb potential
145 :
146 : integer :: lslat_min
147 : ! Min l+1 in the expansion of the Coulomb potential
148 :
149 : integer :: lc_max
150 : ! Max l+1 for orbitals summed over (usually core orbitals)
151 :
152 : integer,allocatable :: sggselect(:,:)
153 : ! sggselect(lslat_max,lc_max)
154 : ! Index of non null sgg, 0 if sgg is zero by symmetry.
155 :
156 : real(dp),allocatable :: sgg(:)
157 : ! sgg(nsggsel)
158 : ! Non null matrix elements in packed form. The index is given by sggselect.
159 :
160 : end type slatang_cshell_t
161 :
162 : public :: slatang_cshell_init ! Creation method for slatang_cshell_t.
163 : public :: slatang_cshell_free ! Destruction method for the slatang_cshell_t.
164 : !!***
165 :
166 : !----------------------------------------------------------------------
167 :
168 : !!****t* m_paw_slater/slatrad_cshell_t
169 : !! NAME
170 : !! slatrad_cshell_t
171 : !!
172 : !! FUNCTION
173 : !! Object used to store the set of radial integrals:
174 : !!
175 : !! $ \dfrac{4\pi}{2L+1} \times
176 : !! \int \phi_{\ni\li}(1) \phic_{\nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{\nc\lc}(2) \phi_{\nj\lj} d1d2 $
177 : !!
178 : !! for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
179 : !! This (less general) type of radial integral is needed to evaluate the Exchange term generated
180 : !! by a closed-shell atom. In the equation, (\nc,\lc) are the set of angular quantum number associated
181 : !! to core electrons while (lsl,msl) comes from the expansion of 1/|r1-r2|.
182 : !! Since the F is invariant under exchange of i and j we use an array of structures indexed
183 : !! by kln = (iln,jln) in packed form.
184 : !!
185 : !! SOURCE
186 :
187 : type,public :: slatrad_cshell_t
188 :
189 : integer :: lnc_size
190 : ! Number of (n,l) channel for core orbitals.
191 :
192 : integer :: lslat_max
193 : ! Max l+1 in the expansion of the Coulomb potential
194 :
195 : integer :: lslat_min
196 : ! Min l+1 in the expansion of the Coulomb potential
197 :
198 : integer :: nrlphic_int
199 : ! The number of non-zero integrals stored in rlphic_int.
200 :
201 : integer,allocatable :: rlphic_select(:,:)
202 : ! rlphic_select(lslat_max,lnc_size) TODO should be allocated with lslat_min:lslat_max
203 : ! Index of the non-zero integrals in rlphic_int, 0 if rlphic_int has not been
204 : ! calculated thanks to selection rules coming from the angular integration.
205 :
206 : real(dp),allocatable :: rlphic_int(:)
207 : ! rlphic_int(1:nrlphic_int)
208 : ! The integrals:
209 : ! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) \phic_{nc\lc}(1) \dfrac{r_<^L}{r_>^{L+1}} \phic_{nc\lc}(2) \phi_{\nj\lj} d1d2
210 : ! for given (in,il) and (jn,jl) as a function L = |il-jl|, |il-jl|+2, ..., |il+il| and ilnc = (lc,nc).
211 :
212 : end type slatrad_cshell_t
213 :
214 :
215 : CONTAINS !========================================================================================
216 : !!***
217 :
218 : !!****f* m_paw_slater/slatang_cshell_init
219 : !! NAME
220 : !! slatang_cshell_init
221 : !!
222 : !! FUNCTION
223 : !! Initialize the structure slatrad_cshell_t containing radial integrals, see below.
224 : !!
225 : !! INPUTS
226 : !! l_max= max. value of ang. momentum l+1;
227 : !! Gaunt coeffs up to [(2*l_max-1,m),(l_max,m),(l_max,m)] are computed
228 : !! lc_max=Max Lc+1 for core states used to contrusct \Sigma_x^\core.
229 : !! lm2_size=Number of symmetrix elements in the (l,m) basis set.
230 : !! NB: lm2_size = (l_max**2)*(l_max**2+1)/2.
231 : !! klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm and jlm for each klm=(ilm,jlm)
232 : !! where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.
233 : !!
234 : !! OUTPUT
235 : !! Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object storing :
236 : !!
237 : !! $ F^{ilsl,ilc}_{klm} = sum_{msl,mc} <li mi|lsl msl;lc mc> <lsl msl;lc mc|lj mj> $
238 : !!
239 : !! where klm = runs over the upper triangle of the ((il,im),(jl,jm)) matrix.
240 : !! ilc runs from 1 up to lc_max and |li-lc| <= lsl <= |li+lc|
241 : !!
242 : !! NOTES
243 : !! Selection rules for F
244 : !! 1) mi = mj
245 : !! 2) In the case of closed shells, one sums for all possible mc"s values from -lc up to +lc.
246 : !! In this particular case, one can use the symmetry properties of Clebsch-Gordan
247 : !! coefficients to show that F is non null only if li==lj. In particular, F can be rewritten as:
248 : !!
249 : !! $ F^{ilsl,ilc}_{klm} = \delta{li,lj}\delta{mi,mj} \times
250 : !! \Gaunt^{lsl,0}_{lc,0;li,0} \sqrt{ \dfrac{(2*lc+1) (2*lsl+1)}{4\pi*(2*li+1)} } $
251 : !!
252 : !! SOURCE
253 :
254 21 : subroutine slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
255 :
256 : !Arguments ------------------------------------
257 : !scalars
258 : integer,intent(in) :: l_max,lc_max,lm2_size
259 : !arrays
260 : integer,intent(in) :: klm2lm(6,lm2_size)
261 : type(slatang_cshell_t),intent(out) :: Slatang3l(lm2_size)
262 :
263 : !Local variables-------------------------------
264 : !scalars
265 : integer :: ilm,ilm0,jlm,lgnt_max,ngnt,ilsl,ilc,lc,ilm0c,lsl
266 : integer :: klm_ci,ilm0sl,li,il,jl,im,jm,ig000
267 : integer :: klm,k0lm_i,k0lm_j,k0lm_c,nsggsel,lslat_max,lslat_min
268 : real(dp) :: dum
269 : !arrays
270 1 : integer,allocatable :: gntselect(:,:)
271 1 : real(dp),allocatable :: realgnt(:),tmp_sgg(:)
272 1 : real(dp),allocatable :: g000(:,:,:)
273 :
274 : ! *************************************************************************
275 :
276 : !@slatang_cshell_t
277 :
278 : ! * Calculate $\Gaunt^{lsl,msl}_{lc,mc;li,mi}$
279 1 : lgnt_max = MAX(l_max,lc_max)
280 3 : ABI_MALLOC( realgnt,((2*lgnt_max-1)**2*(lgnt_max)**4))
281 4 : ABI_MALLOC(gntselect,((2*lgnt_max-1)**2, lgnt_max**2*(lgnt_max**2+1)/2))
282 :
283 1 : call realgaunt(lgnt_max,ngnt,gntselect,realgnt)
284 :
285 : ! Below we need $\Gaunt_{lsl,0}_{lc,0;li,0}$
286 6 : ABI_MALLOC(g000,(2*lgnt_max-1,lc_max,l_max))
287 19 : g000 = zero
288 :
289 4 : do ilsl=1,2*lgnt_max-1
290 3 : lsl = ilsl-1
291 3 : ilm0sl = 1+lsl**2+lsl
292 10 : do il=1,l_max
293 6 : li = il-1
294 6 : ilm0 = 1+li**2+li
295 6 : k0lm_i = ilm0 *(ilm0-1)/2
296 21 : do ilc=1,lc_max
297 12 : lc = ilc-1
298 12 : ilm0c = 1+lc**2+lc
299 12 : k0lm_c= ilm0c * (ilm0c-1)/2
300 12 : if (ilm0c > ilm0) then
301 3 : klm_ci = k0lm_c + ilm0
302 : else
303 9 : klm_ci = k0lm_i + ilm0c
304 : end if
305 12 : ig000 = gntselect(ilm0sl,klm_ci) ! Index of $\Gaunt_{lsl,0}_{lc,0;li,0}$
306 18 : if (ig000 > 0) g000(ilsl,ilc,il)=realgnt(ig000)
307 : end do
308 : end do
309 : end do
310 :
311 1 : ABI_FREE(realgnt)
312 1 : ABI_FREE(gntselect)
313 :
314 : ! === Loop over klm channels in packed form ===
315 11 : do klm=1,lm2_size
316 10 : il = klm2lm(1,klm); im = klm2lm(3,klm)
317 10 : jl = klm2lm(2,klm); jm = klm2lm(4,klm)
318 :
319 10 : nsggsel=0
320 10 : lslat_min = 1 !FIXME find better way
321 10 : lslat_max = il+lc_max-1
322 :
323 10 : Slatang3l(klm)%lslat_min = lslat_min
324 10 : Slatang3l(klm)%lslat_max = lslat_max
325 10 : Slatang3l(klm)%lc_max = lc_max
326 :
327 30 : ABI_MALLOC(tmp_sgg,(lslat_max*lc_max))
328 62 : tmp_sgg = zero
329 :
330 : ! === Calculate F^{lsl,lc}_{li,mi;lj,mj} ===
331 : ! * Selection rule: mi = mj and li==lj
332 10 : if (im == jm .and. il==jl) then
333 4 : li = il-1
334 4 : ilm = klm2lm(5,klm); k0lm_i = ilm *(ilm -1)/2
335 4 : jlm = klm2lm(6,klm); k0lm_j = jlm *(jlm -1)/2
336 :
337 16 : ABI_MALLOC(Slatang3l(klm)%sggselect,(lslat_max,lc_max))
338 34 : Slatang3l(klm)%sggselect = 0
339 :
340 15 : do ilsl=lslat_min,lslat_max
341 : !% do ilsl=lslat_min,lslat_max,2
342 11 : lsl = ilsl-1
343 37 : do ilc=1,lc_max
344 22 : lc = ilc-1
345 22 : dum = SQRT( (two*lc+1)*(two*lsl+1) / (four_pi*(two*li+1)) ) * g000(ilsl,ilc,il)
346 33 : if (ABS(dum)>=tol12) then ! * Store results and progressive index if non null.
347 11 : nsggsel = nsggsel + 1
348 11 : tmp_sgg(nsggsel) = dum
349 11 : Slatang3l(klm)%sggselect(ilsl,ilc) = nsggsel
350 : end if
351 : end do !ilc
352 : end do !ilsl
353 : end if ! Selection rule li=lj and mi == mj
354 : !
355 : ! * Finalize the object.
356 10 : Slatang3l(klm)%nsggsel = nsggsel
357 10 : if (nsggsel > 0) then
358 12 : ABI_MALLOC(Slatang3l(klm)%sgg,(nsggsel))
359 19 : Slatang3l(klm)%sgg = tmp_sgg(1:nsggsel)
360 : end if
361 11 : ABI_FREE(tmp_sgg)
362 : end do !klm
363 :
364 1 : ABI_FREE(g000)
365 :
366 : #if 0
367 : ! Debugging code
368 : do klm=1,lm2_size
369 : if (Slatang3l(klm)%nsggsel>0) then
370 : il = klm2lm(1,klm)
371 : jl = klm2lm(2,klm)
372 : im = klm2lm(3,klm)
373 : jm = klm2lm(4,klm)
374 : write(std_out,*)"--for li, mi",il-1,im-il
375 : lslat_min = Slatang3l(klm)%lslat_min
376 : lslat_max = Slatang3l(klm)%lslat_max
377 :
378 : do ilc=1,lc_max
379 : do ilsl=lslat_min,lslat_max
380 : ii = Slatang3l(klm)%sggselect(ilsl,ilc)
381 : if (ii>0) write(std_out,*)" lc, lslat, sgg",ilc-1,ilsl-1,Slatang3l(klm)%sgg(ii)
382 : end do
383 : end do
384 : end if
385 : end do
386 : #endif
387 :
388 1 : end subroutine slatang_cshell_init
389 : !!***
390 :
391 : !----------------------------------------------------------------------
392 :
393 : !!****f* m_paw_slater/slatang_cshell_free
394 : !! NAME
395 : !! slatang_cshell_free
396 : !!
397 : !! FUNCTION
398 : !! Free the dynamic memory allocated in a structure of type slatang_cshell_t
399 : !!
400 : !! SIDE EFFECTS
401 : !! Slatang3l(lm2_size) <type(slatang_cshell_t)> = Object containing radial integrals
402 : !!
403 : !! SOURCE
404 :
405 1 : subroutine slatang_cshell_free(Slatang3l)
406 :
407 : !Arguments ------------------------------------
408 : !scalars
409 : type(slatang_cshell_t),intent(inout) :: Slatang3l(:)
410 :
411 : !Local variables-------------------------------
412 : integer :: ii
413 : ! *********************************************************************
414 :
415 : !@slatang_cshell_t
416 11 : do ii=1,SIZE(Slatang3l)
417 10 : ABI_SFREE(Slatang3l(ii)%sggselect)
418 11 : ABI_SFREE(Slatang3l(ii)%sgg)
419 : end do
420 :
421 1 : end subroutine slatang_cshell_free
422 : !!***
423 :
424 : !----------------------------------------------------------------------
425 :
426 : !!****f* m_paw_slater/slatrad_cshell_free
427 : !! NAME
428 : !! slatrad_cshell_free
429 : !!
430 : !! FUNCTION
431 : !! Free the dynamic memory allocated in a structure of type slatrad_cshell_t
432 : !!
433 : !! SIDE EFFECTS
434 : !! Slatrad3l(ln2_size) <type(slarad3l_type)> = Object containing radial integrals
435 : !!
436 : !! SOURCE
437 :
438 1 : subroutine slatrad_cshell_free(Slatrad3l)
439 :
440 : !Arguments ------------------------------------
441 : !scalars
442 : type(slatrad_cshell_t),intent(inout) :: Slatrad3l(:)
443 :
444 : !Local variables-------------------------------
445 : integer :: ii
446 : ! *********************************************************************
447 :
448 : !@slatrad_cshell_t
449 11 : do ii=1,SIZE(Slatrad3l)
450 10 : ABI_SFREE(Slatrad3l(ii)%rlphic_select)
451 11 : ABI_SFREE(Slatrad3l(ii)%rlphic_int)
452 : end do
453 :
454 1 : end subroutine slatrad_cshell_free
455 : !!***
456 :
457 : !----------------------------------------------------------------------
458 :
459 : !!****f* m_paw_slater/slatrad_cshell_init
460 : !! NAME
461 : !! slatrad_cshell_init
462 : !!
463 : !! FUNCTION
464 : !! Initialize the structure storing the radial part of Slater"s integrals.
465 : !!
466 : !! INPUTS
467 : !! ln2_size=Number of symmetrical (l,n) channels
468 : !! Pawrad<pawrad_type>=paw radial mesh and related data
469 : !! Pawtab<pawtab_type>=paw tabulated starting data
470 : !! Atm<atomorb_type>=Object containing core orbitals.
471 : !! Atmrad<pawrad_type>=paw radial mesh and related data for the atom.
472 : !! kln_mask
473 : !!
474 : !! OUTPUT
475 : !! Slatrad3l<slatrad_cshell_t>=The object completely initialized.
476 : !!
477 : !! SOURCE
478 :
479 21 : subroutine slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad,kln_mask)
480 :
481 : !Arguments ------------------------------------
482 : !scalars
483 : integer,intent(in) :: ln2_size
484 : !arrays
485 : integer,optional,intent(in) :: kln_mask(ln2_size)
486 : type(atomorb_type),intent(in) :: Atm
487 : type(pawrad_type),target,intent(in) :: Atmrad,Pawrad
488 : type(pawtab_type),target,intent(in) :: Pawtab
489 : type(slatrad_cshell_t),intent(out) :: Slatrad3l(ln2_size)
490 :
491 : !Local variables ---------------------------------------
492 : !scalars
493 : integer :: cmesh_size,dmesh_size
494 : integer :: il,iln,ilnc,isl,in,jl,jln,jn,kln,ll,lnc_size
495 : integer :: lslat_max,lslat_min,lc_max,nintg
496 : integer :: lmn_size,lmn2_size,do_spline,ln_size,whichdenser,isppol
497 : real(dp) :: intg,intg1,ybcbeg,ybcend
498 : logical :: hasameq
499 : !arrays
500 1 : integer,allocatable :: kln2ln(:,:)
501 1 : integer, pointer :: indklmn(:,:),indlmn(:,:)
502 1 : real(dp),allocatable :: ff1(:),ff2(:),tmp_integrals(:)
503 1 : real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:)
504 1 : real(dp),allocatable,target :: phi_spl(:,:)
505 1 : real(dp),allocatable :: der(:),ypp(:)
506 1 : real(dp),ABI_CONTIGUOUS pointer :: crad(:),drad(:),phi_in(:)
507 :
508 : ! *************************************************************************
509 :
510 1 : ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
511 1 : if (PRESENT(kln_mask)) then
512 0 : ABI_ERROR("kln_mask is present")
513 : end if
514 :
515 : !@slatrad_cshell_t
516 1 : lmn_size = Pawtab%lmn_size
517 1 : lmn2_size = Pawtab%lmn2_size
518 1 : ln_size = Pawtab%basis_size
519 :
520 1 : lnc_size = Atm%ln_size
521 1 : lc_max = Atm%l_max
522 :
523 1 : call pawrad_isame(Atmrad,Pawrad,hasameq,whichdenser)
524 :
525 1 : do_spline=0
526 1 : if (.not.hasameq) then
527 1 : do_spline=1
528 1 : if (whichdenser/=1) &
529 1 : & ABI_COMMENT("Pawrad is denser than Atmrad!")
530 : else
531 0 : ABI_CHECK(whichdenser==1,"Pawrad is denser than Atmrad!")
532 : end if
533 :
534 1 : dmesh_size = Atmrad%mesh_size
535 1 : cmesh_size = Pawtab%mesh_size
536 :
537 1 : drad => Atmrad%rad(1:dmesh_size)
538 1 : crad => Pawrad%rad(1:cmesh_size)
539 :
540 : ! === Spline valence basis set onto core mesh (natural spline) ===
541 1 : if (do_spline==1) then
542 1 : ABI_COMMENT("Splining in init_slatrad3l")
543 4 : ABI_MALLOC(phi_spl,(dmesh_size,ln_size))
544 3 : ABI_MALLOC(der,(cmesh_size))
545 2 : ABI_MALLOC(ypp,(cmesh_size))
546 :
547 5 : do iln=1,ln_size
548 4 : phi_in => Pawtab%phi(:,iln)
549 2484 : ypp(:) = zero; ybcbeg = zero; ybcend = zero
550 4 : call spline(crad,phi_in,cmesh_size,ybcbeg,ybcend,ypp)
551 5 : call splint(cmesh_size,crad,phi_in,ypp,dmesh_size,drad,phi_spl(:,iln))
552 : end do
553 :
554 1 : ABI_FREE(der)
555 1 : ABI_FREE(ypp)
556 : end if
557 :
558 1 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
559 1 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
560 :
561 3 : ABI_MALLOC(kln2ln,(6,ln2_size))
562 :
563 1 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
564 :
565 3 : ABI_MALLOC(ff1,(dmesh_size))
566 2 : ABI_MALLOC(ff2,(dmesh_size))
567 :
568 : ! * Loop over the upper triangle of the [(in,il), (jn,il)] matrix.
569 1 : ABI_CHECK(Atm%nsppol==1,"nsppol==2 not tested")
570 :
571 2 : do isppol=1,Atm%nsppol
572 12 : do kln=1,ln2_size
573 10 : il = kln2ln(1,kln)
574 10 : jl = kln2ln(2,kln)
575 10 : in = kln2ln(3,kln)
576 10 : jn = kln2ln(4,kln)
577 10 : iln = kln2ln(5,kln)
578 10 : jln = kln2ln(6,kln)
579 :
580 10 : lslat_max = MAX((il+lc_max),(jl+lc_max)) - 1 ! These are indices, not l-values.
581 : !lslat_min = MIN(ABS(il-lc_max),ABS(jl-lc_max)) + 1
582 10 : lslat_min = 1 ! FIXME find better way
583 :
584 10 : Slatrad3l(kln)%lnc_size = lnc_size
585 10 : Slatrad3l(kln)%lslat_min = lslat_min
586 10 : Slatrad3l(kln)%lslat_max = lslat_max
587 :
588 10 : Slatrad3l(kln)%nrlphic_int = 0
589 :
590 40 : ABI_MALLOC(Slatrad3l(kln)%rlphic_select,(lslat_max,lnc_size))
591 121 : Slatrad3l(kln)%rlphic_select(:,:) = 0
592 :
593 : !if (PRESENT(kln_mask)) then !FIXME THIS IS WRONG, move it below in case
594 : ! if (kln_mask(kln)==0) CYCLE
595 : !end if
596 :
597 10 : if (do_spline==1) then
598 10 : ABI_COMMENT("Performing spline of valence phi")
599 10 : phi_i => phi_spl(:,iln)
600 10 : phi_j => phi_spl(:,jln)
601 : else
602 0 : phi_i => Pawtab%phi(:,iln)
603 0 : phi_j => Pawtab%phi(:,jln)
604 : end if
605 :
606 : ! * Loop over (n,l) channels for Atom orbitals
607 30 : ABI_MALLOC(tmp_integrals,(lslat_max*lnc_size))
608 91 : tmp_integrals(:) = zero
609 : nintg=0
610 :
611 40 : do ilnc=1,lnc_size
612 : ! phicore => Atm%phi(:,ilnc,isppol)
613 41520 : ff1 = phi_i * Atm%phi(:,ilnc,isppol)
614 41520 : ff2 = phi_j * Atm%phi(:,ilnc,isppol)
615 121 : do isl=lslat_min,lslat_max ! L coming from Coulomb expansion
616 81 : ll = isl-1
617 81 : call calc_slatradl(ll,dmesh_size,ff2,ff1,Atmrad,intg1)
618 81 : call calc_slatradl(ll,dmesh_size,ff1,ff2,Atmrad,intg)
619 :
620 : !call calc_slatradl(ll,cmesh_size,ff2,ff1,Pawrad,intg1)
621 : !call calc_slatradl(ll,cmesh_size,ff1,ff2,Pawrad,intg)
622 :
623 81 : if (ABS(intg1-intg)>tol6) write(std_out,*)"DEBUG ",ll,il,in,jl,jn,intg1,intg
624 :
625 : ! * Store results
626 192 : if (ABS(intg)>=tol12) then
627 81 : nintg = nintg +1
628 81 : Slatrad3l(kln)%rlphic_select(isl,ilnc) = nintg
629 81 : tmp_integrals(nintg) = intg
630 : end if
631 : end do !ll
632 : end do ! ilnc
633 :
634 : ! Finalize the object
635 10 : Slatrad3l(kln)%nrlphic_int = nintg
636 30 : ABI_MALLOC(Slatrad3l(kln)%rlphic_int,(nintg))
637 101 : if (nintg>0) Slatrad3l(kln)%rlphic_int(:) = tmp_integrals(1:nintg)
638 :
639 11 : ABI_FREE(tmp_integrals)
640 : end do !kln
641 : end do !isppol
642 :
643 1 : ABI_FREE(ff1)
644 1 : ABI_FREE(ff2)
645 1 : ABI_FREE(kln2ln)
646 :
647 1 : if (do_spline==1) then
648 1 : ABI_FREE(phi_spl)
649 : end if
650 :
651 3 : end subroutine slatrad_cshell_init
652 : !!***
653 :
654 : !----------------------------------------------------------------------
655 :
656 : !!****f* m_paw_slater/paw_sigxcore
657 : !! NAME
658 : !! paw_sigxcore
659 : !!
660 : !! FUNCTION
661 : !! Calculate the integrals:
662 : !! \dfrac{4\pi}{2L+1} \int \phi_{\ni\li}(1) orb_{nl}(1) \dfrac{r_<^L}{r_>^{L+1}} orb_{nl}(2) \phi_{\nj\lj} d1d2!!
663 : !! for given (in,il) and (jn,jl) as a function of (nc,lc) and L = |il-jl|, |il-jl|+2, ..., |il+il|
664 : !!
665 : !! INPUTS
666 : !! cplex_dij=1 if dijexc_core is real, 2 if they are complex
667 : !! lmn2_size=Number of (klmn) channels
668 : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
669 : !! Pawtab<pawtab_type>=paw tabulated starting data
670 : !! Atm<atomorb_type>=Structure containing core orbitals
671 : !! Atmrad<pawrad_type>=The radial mesh for core orbitals
672 : !!
673 : !! OUTPUT
674 : !! dijexc_core(cplex_dij*lmn2_size,ndij)
675 : !!
676 : !! SOURCE
677 :
678 1 : subroutine paw_sigxcore(cplex_dij,lmn2_size,ndij,Pawrad,Pawtab,Atm,Atmrad,dijexc_core)
679 :
680 : !Arguments ------------------------------------
681 : !scalars
682 : integer,intent(in) :: lmn2_size,cplex_dij,ndij
683 : !arrays
684 : real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size,ndij)
685 : type(atomorb_type),intent(in) :: Atm
686 : type(pawrad_type),intent(in) :: Atmrad
687 : type(pawrad_type),intent(in) :: Pawrad
688 : type(pawtab_type),target,intent(in) :: Pawtab
689 :
690 : !Local variables ---------------------------------------
691 : !scalars
692 : integer :: ilnc,ilc,lnc_size,l_max
693 : integer :: lm2_size,ln_size,ln2_size,lmn_size
694 : integer :: lm_size,klmn,kln,klm
695 : integer :: lc_max,ilsl,isgg,israd,opt_l,pawprtvol
696 : real(dp) :: tmp,sgg,intgrl
697 : !character(len=500) :: msg
698 : !arrays
699 : integer :: opt_l_index(0,0),pack2ij(0)
700 1 : integer,allocatable :: kln2ln(:,:),klm2lm(:,:)
701 : integer, pointer :: indklmn(:,:),indlmn(:,:)
702 : type(slatang_cshell_t),allocatable :: Slatang3l(:)
703 1 : type(slatrad_cshell_t),allocatable :: Slatrad3l(:)
704 :
705 : ! *************************************************************************
706 :
707 : ! * Consistency check
708 1 : ABI_CHECK(cplex_dij==1,"cplex_dij must be 1")
709 :
710 1 : ABI_CHECK(ndij==1,"ndij must be 1")
711 :
712 1 : ABI_CHECK(lmn2_size==Pawtab%lmn2_size,"Wrong lmn2_size")
713 :
714 1 : lmn_size = Pawtab%lmn_size
715 1 : ln_size = Pawtab%basis_size
716 1 : ln2_size = Pawtab%ij_size
717 1 : l_max = (Pawtab%l_size-1)/2 +1
718 1 : lm_size = l_max**2
719 1 : lm2_size = lm_size*(lm_size+1)/2
720 :
721 1 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
722 1 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
723 :
724 : ! * Setup of useful tables.
725 3 : ABI_MALLOC(kln2ln,(6,ln2_size))
726 1 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
727 :
728 3 : ABI_MALLOC(klm2lm,(6,lm2_size))
729 1 : call make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)
730 :
731 : ! * Integrate angular part.
732 1 : lnc_size = Atm%ln_size
733 1 : lc_max = Atm%l_max
734 :
735 13 : ABI_MALLOC(Slatang3l,(lm2_size))
736 1 : call slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
737 :
738 1 : ABI_FREE(klm2lm)
739 :
740 : ! * Integrate radial part.
741 13 : ABI_MALLOC(Slatrad3l,(ln2_size))
742 :
743 1 : call slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad)
744 :
745 : ! === Calculate matrix elements of Sigma_x^core ===
746 : ! * $<\phi_i|\Sigma_x^\core|\phi_j>$
747 38 : dijexc_core = zero
748 37 : do klmn=1,lmn2_size
749 36 : klm = Pawtab%indklmn(1,klmn)
750 36 : kln = Pawtab%indklmn(2,klmn)
751 : !
752 : ! * Summing over (lc,nc) and lslat
753 36 : tmp = zero
754 36 : if (Slatang3l(klm)%nsggsel >0) then
755 48 : do ilnc=1,Atm%ln_size
756 36 : ilc = 1+Atm%indln(1,ilnc)
757 147 : do ilsl=1,Slatang3l(klm)%lslat_max !FIXME check this
758 : !do ilsl=Slatang3l(klm)%lslat_min,Slatang3l(klm)%lslat_max
759 99 : isgg = Slatang3l(klm)%sggselect(ilsl,ilc)
760 99 : israd = Slatrad3l(kln)%rlphic_select(ilsl,ilnc)
761 135 : if (isgg>0 .and. israd>0) then
762 45 : sgg = Slatang3l(klm)%sgg(isgg)
763 45 : intgrl = Slatrad3l(kln)%rlphic_int(israd)
764 45 : tmp = tmp + intgrl*sgg
765 : end if
766 : end do
767 : end do
768 : end if
769 :
770 37 : dijexc_core(klmn,1) = -tmp ! Store results.
771 : end do
772 :
773 : ! * Print values
774 1 : call wrtout(std_out," ************** Dij Fock_core ************ ",'COLL')
775 1 : opt_l=-1; pawprtvol=-1
776 1 : call pawio_print_ij(std_out,dijexc_core(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)
777 :
778 : ! * Free memory.
779 1 : ABI_FREE(kln2ln)
780 1 : call slatang_cshell_free(Slatang3l)
781 11 : ABI_FREE(Slatang3l)
782 1 : call slatrad_cshell_free(Slatrad3l)
783 11 : ABI_FREE(Slatrad3l)
784 :
785 1 : end subroutine paw_sigxcore
786 : !!***
787 :
788 : !----------------------------------------------------------------------
789 :
790 : !!****f* m_paw_slater/paw_mkdijexc_core
791 : !! NAME
792 : !! paw_mkdijexc_core
793 : !!
794 : !! FUNCTION
795 : !! Driver routine to calculate the onsite matrix element of the Fock operator between two
796 : !! all-electron partial waves.
797 : !!
798 : !! INPUTS
799 : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
800 : !! cplex_dij=1 if dijexc_core is real, 2 if they are complex
801 : !! lmn2_size_max=Max Number of (klmn) channels over type of atoms.
802 : !! Cryst<crystal_t>=Structure describing the crystal structure and its symmmetries.
803 : !! Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
804 : !! Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
805 : !! pawprtvol=Flags governing the verbosity of the output.
806 : !! filpsp(ntypat)=names of the files containing the all-electron core WF
807 : !!
808 : !! OUTPUT
809 : !! dijexc_core(cplex_dij*lmn2_size_max,ndij,ntypat)= On-site matrix elements $ \<\phi_i|Sigma_x^\core|\phi_j\>
810 : !! for each type of atom.
811 : !!
812 : !! SOURCE
813 :
814 1 : subroutine paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,pawprtvol,filpsp)
815 :
816 : !Arguments ------------------------------------
817 : !scalars
818 : integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max
819 : type(crystal_t),intent(in) :: Cryst
820 : !arrays
821 : real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) !TODO use ragged arrays pawij?
822 : character(len=fnlen) :: filpsp(Cryst%ntypat)
823 : type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
824 : type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
825 :
826 : !Local variables ---------------------------------------
827 : !scalars
828 : integer :: itypat,lmn_size,lmn2_size,ln_size,isppol
829 : real(dp) :: rcut
830 : character(len=500) :: header
831 : !arrays
832 1 : integer,allocatable :: phi_indln(:,:)
833 1 : real(dp),ABI_CONTIGUOUS pointer :: phi(:,:)
834 1 : real(dp),allocatable :: overlap(:,:)
835 1 : type(atomorb_type),allocatable :: Atm(:)
836 1 : type(Pawrad_type),allocatable :: Radatm(:)
837 :
838 : ! *************************************************************************
839 :
840 4 : ABI_MALLOC(Atm,(Cryst%ntypat))
841 4 : ABI_MALLOC(Radatm,(Cryst%ntypat))
842 :
843 1 : ABI_CHECK(ndij==1 ,"spinor+HF not available")
844 1 : ABI_CHECK(cplex_dij==1,"spinor+HF not available")
845 2 : ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
846 :
847 : !allocate(dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) !TODO use ragged arrays pawij?
848 39 : dijexc_core=zero
849 :
850 2 : do itypat=1,Cryst%ntypat
851 :
852 : ! Read core orbitals for this atom type.
853 1 : rcut=Pawtab(itypat)%rpaw
854 1 : call pawpsp_init_core(Atm(itypat),psp_filename=trim(filpsp(itypat)),rcut_in=rcut)
855 1 : call pawrad_copy(Atm(itypat)%radmesh,Radatm(itypat))
856 :
857 1 : write(header,'(a,i4,a)')" === Atom type = ",itypat," === "
858 1 : call print_atomorb(Atm(itypat),header,unit=std_out,prtvol=pawprtvol)
859 : !
860 : ! * Calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for this atom type.
861 1 : lmn_size = Pawtab(itypat)%lmn_size
862 1 : lmn2_size = Pawtab(itypat)%lmn2_size
863 :
864 : call paw_sigxcore(cplex_dij,lmn2_size,ndij,&
865 1 : & Pawrad(itypat),Pawtab(itypat),Atm(itypat),Radatm(itypat),dijexc_core(1:lmn2_size,:,itypat))
866 :
867 1 : ln_size = Pawtab(itypat)%basis_size
868 1 : phi => Pawtab(itypat)%phi
869 :
870 3 : ABI_MALLOC(phi_indln,(2,ln_size))
871 1 : call make_indln(lmn_size,ln_size,Pawtab(itypat)%indlmn(:,:),phi_indln)
872 :
873 4 : ABI_MALLOC(overlap,(Atm(itypat)%ln_size,ln_size))
874 1 : isppol=1 ! hardcoded
875 1 : call get_overlap(Atm(itypat),Radatm(itypat),Pawrad(itypat),isppol,ln_size,phi,phi_indln,overlap)
876 :
877 1 : ABI_FREE(phi_indln)
878 2 : ABI_FREE(overlap)
879 : end do ! ntypat
880 :
881 : ! Free memory
882 1 : call pawrad_free(Radatm)
883 2 : do itypat=1,Cryst%ntypat
884 2 : call destroy_atomorb(Atm(itypat))
885 : end do
886 :
887 2 : ABI_FREE(Atm)
888 2 : ABI_FREE(Radatm)
889 :
890 1 : end subroutine paw_mkdijexc_core
891 : !!****
892 :
893 : !----------------------------------------------------------------------
894 :
895 : !!****f* m_paw_slater/slatrad_free_0D
896 : !! NAME
897 : !! slatrad_free_0D
898 : !!
899 : !! FUNCTION
900 : !! Free the dynamic memory allocated in a structure of type slatrad_t
901 : !!
902 : !! SOURCE
903 :
904 0 : subroutine slatrad_free_0D(Slatrad)
905 :
906 : !Arguments ------------------------------------
907 : !scalars
908 : type(slatrad_t),intent(inout) :: Slatrad
909 :
910 : ! *********************************************************************
911 :
912 : !@slatrad_t
913 0 : ABI_SFREE(Slatrad%intgl_select)
914 0 : ABI_SFREE(Slatrad%intgl)
915 :
916 0 : end subroutine slatrad_free_0D
917 : !!***
918 :
919 : !----------------------------------------------------------------------
920 :
921 : !!****f* m_paw_slater/slatrad_free_1D
922 : !! NAME
923 : !! slatrad_free_1D
924 : !!
925 : !! FUNCTION
926 : !! Free the dynamic memory allocated in a structure of type slatrad_t
927 : !!
928 : !! SOURCE
929 :
930 0 : subroutine slatrad_free_1D(Slatrad)
931 :
932 : !Arguments ------------------------------------
933 : !scalars
934 : type(slatrad_t),intent(inout) :: Slatrad(:)
935 :
936 : !Local variables-------------------------------
937 : integer :: ii
938 : ! *********************************************************************
939 :
940 0 : do ii=1,SIZE(Slatrad)
941 0 : call slatrad_free_0D(Slatrad(ii))
942 : end do
943 :
944 0 : end subroutine slatrad_free_1D
945 : !!***
946 :
947 : !----------------------------------------------------------------------
948 :
949 : !!****f* m_paw_slater/slatrad_init
950 : !! NAME
951 : !! slatrad_init
952 : !!
953 : !! FUNCTION
954 : !! Initialize the structure storing the radial part of Slater"s integrals.
955 : !!
956 : !! INPUTS
957 : !! which_intg= Option defining what kind of integrals have to be calculated:
958 : !! -- 1 for Slater integral of AE partial waves only.
959 : !! A = \frac{4\pi}{2L+1} \int u_i(1) u_j(1) \frac{r_<^L}{r_^{L+1}} u_k(2)u_l(2) d1d2
960 : !1
961 : !! -- 2 for Slater integral of (AE-PS) partial waves
962 : !! B = \frac{4\pi}{2L+1} \int u_i(1) u_j(1) \frac{r_<^L}{r_^{L+1}} u_k(2) u_l(2) d1d2 -
963 : !! \frac{4\pi}{2L+1} \int tu_i(1) tu_j(1) \frac{r_<^L}{r_^{L+1}} tu_k(2) tu_l(2) d1d2
964 : !!
965 : !! -- 3 for Slater integral of (AE-PS-compensation charges)
966 : !! C = A -
967 : !! \frac{4\pi}{2L+1} \int [tu_i(1) tu_j(1) + qhat^L_\ij r_1^2 g^L(1) ] \frac{r_<^L}{r_^{L+1}}*
968 : !! [tu_k(2) tu_l(2) + qhat^L_\kl r_2^2 g^L(2) ] d1d2
969 : !!
970 : !! where u = \phi/r; tu = \tphi/r; and qhat^L_\ij are related to q^\LM_\ij via
971 : !! q^\LM_\ij = \Gaunt_\ij^\LM qhat^L_\ij => qhat^L\ij = \int (u_i*u_j - tu_i*tu_j) r^L
972 : !!
973 : !! ln2_size=Number of symmetrical (l,n) channels for this atom type type.
974 : !! Pawrad<pawrad_type>=paw radial mesh and related data
975 : !! Pawtab<pawtab_type>=paw tabulated starting data
976 : !!
977 : !! OUTPUT
978 : !! Slatrad4<slatrad_t>=The object completely initialized.
979 : !!
980 : !! NOTES
981 : !! Slater integrals S_ij are invariant under exchange of the indices,
982 : !! but the results reported by calc_slatradl are not due to numerical roundoff errors (err < 10^-9).
983 : !! However this does not cause any problem since only the upper triangle of the S_ij matrix
984 : !! is stored and used in the other routines.
985 : !!
986 : !! SOURCE
987 :
988 0 : subroutine slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad,Pawtab)
989 :
990 : !Arguments ------------------------------------
991 : !scalars
992 : integer,intent(in) :: ln2_size,which_intg
993 : !arrays
994 : type(pawrad_type),target,intent(in) :: Pawrad
995 : type(pawtab_type),target,intent(in) :: Pawtab
996 : type(slatrad_t),intent(out) :: Slatrad4(ln2_size*(ln2_size+1)/2)
997 :
998 : !Local variables ---------------------------------------
999 : !scalars
1000 : integer :: mesh_size,il,iln,isl,in,jl,jln,jn,sln1,sln2,l_slat
1001 : integer :: kn,kl,ln,ll,kln,lln,lslat_max,lslat_min,nintgl
1002 : integer :: lmn_size,lmn2_size,ln_size,slt_idx
1003 : real(dp) :: ae_intg,ps_intg,pshat_intg,intg,tqij_L,tqkl_L !intg1
1004 : character(len=500) :: msg
1005 : !arrays
1006 0 : integer,allocatable :: kln2ln(:,:)
1007 0 : integer, pointer :: indklmn(:,:),indlmn(:,:)
1008 0 : real(dp),allocatable :: uiuj(:),ukul(:),tuituj(:),tuktul(:),tuituj_tqgl(:),tuktul_tqgl(:)
1009 0 : real(dp),allocatable :: tmp_integrals(:),ff(:)
1010 0 : real(dp),ABI_CONTIGUOUS pointer :: phi_i(:),phi_j(:),phi_k(:),phi_l(:)
1011 0 : real(dp),ABI_CONTIGUOUS pointer :: tphi_i(:),tphi_j(:),tphi_k(:),tphi_l(:)
1012 0 : real(dp),ABI_CONTIGUOUS pointer :: shapefunc(:),rad(:)
1013 :
1014 : ! *************************************************************************
1015 :
1016 : DBG_ENTER("COLL")
1017 :
1018 0 : ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
1019 :
1020 0 : if ( ALL(which_intg /= (/1,2,3/)) ) then
1021 0 : write(msg,'(a,i0)')"Wrong value for which_intg: ",which_intg
1022 0 : ABI_ERROR(msg)
1023 : end if
1024 :
1025 : !@slatrad_t
1026 0 : lmn_size = Pawtab%lmn_size
1027 0 : lmn2_size = Pawtab%lmn2_size
1028 0 : ln_size = Pawtab%basis_size
1029 0 : mesh_size = Pawtab%mesh_size
1030 : !
1031 : ! Useful table for looping.
1032 0 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
1033 0 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
1034 :
1035 0 : ABI_MALLOC(kln2ln,(6,ln2_size))
1036 0 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
1037 :
1038 0 : ABI_MALLOC(uiuj,(mesh_size))
1039 0 : ABI_MALLOC(ukul,(mesh_size))
1040 0 : ABI_MALLOC(ff,(mesh_size))
1041 0 : ABI_MALLOC(tuituj,(mesh_size))
1042 0 : ABI_MALLOC(tuktul,(mesh_size))
1043 0 : ABI_MALLOC(tuituj_tqgl,(mesh_size))
1044 0 : ABI_MALLOC(tuktul_tqgl,(mesh_size))
1045 0 : rad => Pawrad%rad
1046 : !
1047 : ! * Loop over (k,l) channels in packed form.
1048 0 : do sln2=1,ln2_size
1049 0 : kl = kln2ln(1,sln2)
1050 0 : ll = kln2ln(2,sln2)
1051 0 : kn = kln2ln(3,sln2)
1052 0 : ln = kln2ln(4,sln2)
1053 0 : kln = kln2ln(5,sln2)
1054 0 : lln = kln2ln(6,sln2)
1055 : !write(std_out,*)"sln2, kln, lln",sln2,kln,lln
1056 :
1057 0 : phi_k => Pawtab%phi (:,kln)
1058 0 : tphi_k => Pawtab%tphi(:,kln)
1059 :
1060 0 : phi_l => Pawtab%phi (:,lln)
1061 0 : tphi_l => Pawtab%tphi(:,lln)
1062 : !
1063 : ! * Loop over (i,j) channels in packed form AND only for the upper triangle of sln2, sln1
1064 0 : do sln1=1,sln2
1065 0 : il = kln2ln(1,sln1)
1066 0 : jl = kln2ln(2,sln1)
1067 0 : in = kln2ln(3,sln1)
1068 0 : jn = kln2ln(4,sln1)
1069 0 : iln = kln2ln(5,sln1)
1070 0 : jln = kln2ln(6,sln1)
1071 : !write(std_out,*)"sln1, iln, jln",sln1,iln,jln
1072 :
1073 0 : slt_idx = sln1 + sln2*(sln2-1)/2 ! index for packed storage.
1074 :
1075 0 : phi_i => Pawtab%phi (:,iln)
1076 0 : tphi_i => Pawtab%tphi(:,iln)
1077 :
1078 0 : phi_j => Pawtab%phi (:,jln)
1079 0 : tphi_j => Pawtab%tphi(:,jln)
1080 :
1081 0 : lslat_min = MAX(ABS(il-jl),ABS(kl-ll)) + 1 ! We use indices not l-values.
1082 0 : lslat_max = MIN((il+jl),(kl+ll)) - 1
1083 :
1084 : !lslat_min = MIN(ABS(il-jl),ABS(kl-ll)) + 1
1085 : !lslat_max = MAX((il+jl),(kl+ll)) - 1
1086 :
1087 0 : Slatrad4(slt_idx)%lslat_min = lslat_min
1088 0 : Slatrad4(slt_idx)%lslat_max = lslat_max
1089 :
1090 0 : Slatrad4(slt_idx)%iln = iln
1091 0 : Slatrad4(slt_idx)%jln = jln
1092 0 : Slatrad4(slt_idx)%kln = kln
1093 0 : Slatrad4(slt_idx)%lln = lln
1094 :
1095 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl_select,(lslat_min:lslat_max))
1096 0 : Slatrad4(slt_idx)%intgl_select(:) = 0
1097 0 : Slatrad4(slt_idx)%nintgl = 0
1098 :
1099 0 : if (lslat_min > lslat_max) then
1100 : ! e.g. (1 2) (1 1). Due to angular selection rules, this integral do not contribue
1101 : !write(std_out,*)"lslat_min, lslat_max",lslat_min,lslat_max
1102 : !write(std_out,*)"il,jl,kl,ll",il,jl,kl,ll
1103 : !ABI_ERROR("")
1104 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl,(0))
1105 0 : CYCLE
1106 : end if
1107 :
1108 0 : uiuj = phi_i * phi_j ! The AE part.
1109 0 : ukul = phi_k * phi_l
1110 0 : tuituj = tphi_i * tphi_j ! The pseudized part.
1111 0 : tuktul = tphi_k * tphi_l
1112 : !
1113 : ! Calculate L-depedent integrals where L come from the expansion the Coulomb interaction.
1114 0 : ABI_MALLOC(tmp_integrals,(MAX(lslat_min,lslat_max)))
1115 0 : tmp_integrals=zero
1116 : nintgl=0
1117 :
1118 0 : do isl=lslat_min,lslat_max
1119 : !do isl=lslat_min,lslat_max,2 ! TODO Here I can reduce the number of iterations using a step of 2.
1120 0 : l_slat = isl-1
1121 0 : call calc_slatradl(l_slat,mesh_size,uiuj,ukul,Pawrad,ae_intg)
1122 0 : intg = ae_intg
1123 :
1124 : #if 0
1125 : call calc_slatradl(l_slat,mesh_size,ukul,uiuj,Pawrad,intg1)
1126 : if (ABS(intg1-ae_intg)>tol12) then
1127 : write(msg,'(a,es16.8)')"s_ij and s_ij differ by ",intg1-ae_intg
1128 : ABI_WARNING(msg)
1129 : end if
1130 : #endif
1131 0 : if (which_intg == 2) then ! Subtract the pseudo part.
1132 0 : call calc_slatradl(l_slat,mesh_size,tuituj,tuktul,Pawrad,ps_intg)
1133 0 : intg = intg - ps_intg
1134 :
1135 0 : else if (which_intg == 3) then ! Subtract (pseudo + compensation charges)
1136 : !
1137 : ! Evaluate tqij_L and tqkl_L (without M-dependent part).
1138 0 : ff(1)=zero
1139 0 : ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln1)-pawtab%tphitphj(2:mesh_size,sln1))*rad(2:mesh_size)**l_slat
1140 0 : if (l_slat==0.and.kl==1.and.ll==1) then
1141 0 : call pawrad_deducer0(ff,mesh_size,pawrad)
1142 : end if
1143 0 : call simp_gen(tqij_L,ff,pawrad)
1144 :
1145 0 : ff(1)=zero
1146 0 : ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln2)-pawtab%tphitphj(2:mesh_size,sln2))*rad(2:mesh_size)**l_slat
1147 0 : if (l_slat==0.and.il==1.and.jl==1) then
1148 0 : call pawrad_deducer0(ff,mesh_size,pawrad)
1149 : end if
1150 0 : call simp_gen(tqkl_L,ff,pawrad)
1151 :
1152 0 : shapefunc => Pawtab%shapefunc(:,isl) ! Recheck this part, in particular the convention
1153 0 : tuituj_tqgl = tuituj + tqij_L * shapefunc * rad**2
1154 0 : tuktul_tqgl = tuktul + tqkl_L * shapefunc * rad**2
1155 :
1156 0 : call calc_slatradl(l_slat,mesh_size,tuituj_tqgl,tuktul_tqgl,Pawrad,pshat_intg)
1157 0 : intg = intg - pshat_intg
1158 : end if
1159 : !
1160 : ! * Store results
1161 0 : if (ABS(intg)>=tol12) then
1162 0 : nintgl = nintgl +1
1163 0 : Slatrad4(slt_idx)%intgl_select(isl) = nintgl
1164 0 : tmp_integrals(nintgl) = intg
1165 : end if
1166 : end do !isl
1167 : !
1168 : ! Finalize the object.
1169 0 : Slatrad4(slt_idx)%nintgl = nintgl
1170 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl,(nintgl))
1171 0 : if (nintgl>0) Slatrad4(slt_idx)%intgl(:) = tmp_integrals(1:nintgl)
1172 0 : ABI_FREE(tmp_integrals)
1173 : end do !sln1
1174 : end do !sln2
1175 : !
1176 : ! Free memory
1177 0 : ABI_FREE(kln2ln)
1178 0 : ABI_FREE(uiuj)
1179 0 : ABI_FREE(ukul)
1180 0 : ABI_FREE(tuituj)
1181 0 : ABI_FREE(tuktul)
1182 0 : ABI_FREE(ff)
1183 0 : ABI_FREE(tuituj_tqgl)
1184 0 : ABI_FREE(tuktul_tqgl)
1185 :
1186 : DBG_EXIT("COLL")
1187 :
1188 0 : end subroutine slatrad_init
1189 : !!***
1190 :
1191 : !----------------------------------------------------------------------
1192 :
1193 : !!****f* m_paw_slater/paw_dijhf
1194 : !! NAME
1195 : !! paw_dihf
1196 : !!
1197 : !! FUNCTION
1198 : !! This routine calculates the onsite D_{ij} strengths of the exchange part of the self energy.
1199 : !!
1200 : !! INPUTS
1201 : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
1202 : !! cplex_dij=1 if sigx_dij is real, 2 if they are complex
1203 : !! qphase=2 if dij contains a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
1204 : !! lmn2_size_max=Max Number of (klmn) channels over type of atoms.
1205 : !! my_natom=number of atoms treated by current process
1206 : !! ntypat=number of atom types
1207 : !! Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
1208 : !! Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
1209 : !! Pawang<type(pawang_type)>=paw angular mesh and related data
1210 : !! pawprtvol=Flags governing the verbosity of the output.
1211 : !!
1212 : !! OUTPUT
1213 : !! sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom)=
1214 : !! For each atom, the Pseudopotential strengths of the on-site operator Sigma_x
1215 : !!
1216 : !! NOTES
1217 : !! The on-site contribution to the matrix elements of the exchange part of the self-energy is given by:
1218 : !! <\tpsi_a| [\sum_{ij} |tprj_i\> D_{ij} \<tprj_j|] |\tpsi_b\>.
1219 : !!
1220 : !! When compensation charges are used one obtains:
1221 : !!
1222 : !! D_{ij} = - sum_{kl} \rho_lk [ \Phi_{ikjl} - \Phihat_{ijkl} =
1223 : !!
1224 : !! = - sum_{kl} \rho_lk \sum_{LM} \Gaunt_{ik}^{LM} \Gaunt_{jl}^{LM} [S_{ikjl}^L - tS{ikjl}^L}]
1225 : !!
1226 : !! where S and tS are Slater-like integrals given by
1227 : !!
1228 : !! 1) S_{ijkl}^L = dfrac{4\pi}{2L+1} \iint u_i(1)u_j(1) u_k(2) u_l(2) \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
1229 : !! 1) tS_{ijkl}^L = dfrac{4\pi}{2L+1} \iint [u_i(1)u_j(1)+ tq_{ij}^L g^L(1)]
1230 : !! [u_k(2)u_l(2)+ tq_{kl}^L g^L(2)] \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
1231 : !!
1232 : !! tq_{ij}^L is defined in terms of q_{ij}^L via: q_{ij]^{LM} = tq_{ij}^L \Gaunt_{ij}^{LM}
1233 : !!
1234 : !! SOURCE
1235 :
1236 0 : subroutine paw_dijhf(ndij,cplex_dij,qphase,lmn2_size_max,my_natom,ntypat,Pawtab,Pawrad,Pawang,Pawrhoij,&
1237 0 : & sigx_dij,pawprtvol)
1238 :
1239 : !Arguments ------------------------------------
1240 : !scalars
1241 : integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max,my_natom,ntypat,qphase
1242 : type(pawang_type),intent(in) :: Pawang
1243 : !arrays
1244 : real(dp),target,intent(out) :: sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom) !TODO use ragged arrays pawij?
1245 : type(pawtab_type),intent(in) :: Pawtab(ntypat)
1246 : type(pawrad_type),intent(in) :: Pawrad(ntypat)
1247 : type(pawrhoij_type),intent(in) :: Pawrhoij(my_natom)
1248 :
1249 : !Local variables ---------------------------------------
1250 : !scalars
1251 : integer,parameter :: cplex=1 ! FIXME preliminary implementation
1252 : integer :: cplex_rhoij,iatom,iq,iq0_dij,iq0_rhoij,itypat,lmn_size,lmn2_size,ispden,nspden,ln2_size
1253 : integer :: lm2_size !,isppol ln_size,
1254 : integer :: irhoij,jrhoij
1255 : integer :: rho_lmn !,rho_klm,rho_kln,rho_lmin,rho_lmax,rho_iln,rho_jln
1256 : integer :: klmn
1257 : integer :: i_lmn,j_lmn,k_lmn,l_lmn
1258 : integer :: which_intg,l_max,opt_l
1259 : real(dp) :: ro,slt_ikjl,slt_iljk
1260 : !character(len=500) :: msg
1261 : !arrays
1262 : integer :: opt_l_index(0,0),pack2ij(0)
1263 0 : real(dp), ABI_CONTIGUOUS pointer :: sigx_atm(:,:)
1264 0 : type(slatrad_t),allocatable :: Slatrad4(:)
1265 :
1266 : ! *************************************************************************
1267 :
1268 : DBG_ENTER("COLL")
1269 :
1270 0 : ABI_CHECK(ndij/=4,"ndij=4 not coded")
1271 0 : ABI_CHECK(cplex_dij==1,"cplex_dij/=2 not coded")
1272 0 : ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
1273 :
1274 0 : if (my_natom>0) then
1275 0 : if (pawrhoij(1)%qphase<cplex) then
1276 0 : ABI_BUG('Must have pawrhoij()%qphase >= cplex !')
1277 : end if
1278 : end if
1279 :
1280 0 : sigx_dij=zero
1281 :
1282 0 : do iatom=1,my_natom
1283 0 : itypat =Pawrhoij(iatom)%itypat
1284 0 : lmn_size =Pawtab(itypat)%lmn_size
1285 0 : lmn2_size=Pawtab(itypat)%lmn2_size
1286 0 : l_max =(Pawtab(itypat)%l_size+1)/2
1287 0 : lm2_size = (l_max**2)*(l_max**2+1)/2
1288 0 : cplex_rhoij=Pawrhoij(iatom)%cplex_rhoij
1289 : !write(std_out,*)"in atom ",iatom,"lm2_size=",lm2_size
1290 :
1291 : ! Calculate Slater integral for this atom type.
1292 : ! TODO obviously these tables should be stored in Pawtab!
1293 0 : ln2_size = Pawtab(itypat)%ij_size
1294 0 : ABI_MALLOC(Slatrad4,(ln2_size*(ln2_size+1)/2))
1295 0 : which_intg=3
1296 0 : call slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad(itypat),Pawtab(itypat))
1297 :
1298 0 : sigx_atm => sigx_dij(:,:,iatom)
1299 :
1300 : ! Loop over phase exp(iqr) phase real/imaginary part, if any
1301 0 : do iq=1,qphase
1302 : !First loop: we store the real part in dij(1 -> lmn2_size)
1303 : !2nd loop: we store the imaginary part in dij(lmn2_size+1 -> 2*lmn2_size)
1304 0 : iq0_dij=merge(0,cplex_dij*lmn2_size,iq==1)
1305 0 : iq0_rhoij=merge(0,cplex_rhoij*lmn2_size,iq==1)
1306 :
1307 : ! * Loop over spin components.
1308 0 : nspden=ndij
1309 0 : do ispden=1,ndij
1310 : !
1311 : ! ============================================================
1312 : ! ==== Summing over the non-zero lk channels of \rho_{lk} ====
1313 : ! ============================================================
1314 0 : jrhoij=1+iq0_rhoij
1315 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
1316 0 : rho_lmn=pawrhoij(iatom)%rhoijselect(irhoij)
1317 :
1318 : ! check wheter rho_lmin is consistent with the Indexing used in slatrad
1319 : !rho_klm =pawtab(itypat)%indklmn(1,rho_lmn)
1320 : !rho_kln =pawtab(itypat)%indklmn(2,rho_lmn)
1321 : !rho_lmin=pawtab(itypat)%indklmn(3,rho_lmn)
1322 : !rho_lmax=pawtab(itypat)%indklmn(4,rho_lmn)
1323 :
1324 : ! Retrieve rhoij for this ispden.
1325 0 : if (nspden/=2) then
1326 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1327 : else
1328 0 : ABI_ERROR("Recheck this part")
1329 0 : if (ispden==1) then
1330 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1) + pawrhoij(iatom)%rhoijp(jrhoij,2)
1331 0 : else if (ispden==2) then
1332 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
1333 : end if
1334 : end if
1335 : !
1336 : ! Avoid double-counting the diagonal of rho.
1337 0 : ro=ro*pawtab(itypat)%dltij(rho_lmn)*half
1338 :
1339 0 : call klmn2ijlmn(rho_lmn,lmn_size,k_lmn,l_lmn)
1340 :
1341 : ! Loop over the upper triangle of the D_{ij) matrix and accumulate:
1342 : ! sum_\lk rho_\kl [ \Phi_{ikjl} + \Phi_{iljk} - \Phihat_{ikjl} - \Phihat_{iljk} ]
1343 0 : do klmn=1,lmn2_size
1344 : ! Calculate the indices in the Slatrad4 structure.
1345 0 : call klmn2ijlmn(klmn,lmn_size,i_lmn,j_lmn)
1346 :
1347 : ! My formula
1348 0 : slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,j_lmn,l_lmn)
1349 0 : slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,j_lmn,k_lmn)
1350 :
1351 : !slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,l_lmn,j_lmn)
1352 : !slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,k_lmn,j_lmn)
1353 : !slt_iljk = slt_ikjl
1354 :
1355 0 : sigx_atm(klmn+iq0_dij,ispden) = sigx_atm(klmn,ispden) + ro * (slt_ikjl + slt_iljk)
1356 : end do ! klmn
1357 :
1358 0 : jrhoij=jrhoij+cplex_rhoij
1359 : end do ! irhoij
1360 : end do ! iq
1361 : end do ! ispden
1362 :
1363 0 : if (ABS(pawprtvol)>=1) then ! * Print values
1364 0 : call wrtout(std_out," ************** Dij Fock ************ ",'COLL')
1365 0 : opt_l=-1
1366 0 : call pawio_print_ij(std_out,sigx_atm(:,1),lmn2_size,cplex_dij,lmn_size,opt_l,opt_l_index,0,pawprtvol,pack2ij,-one,1)
1367 : end if
1368 :
1369 0 : call slatrad_free(Slatrad4)
1370 0 : ABI_FREE(Slatrad4)
1371 : end do ! iatom
1372 :
1373 : ! Factor half cancels in the derivation wrt rho_ij.
1374 0 : sigx_dij = - sigx_dij
1375 :
1376 : DBG_EXIT("COLL")
1377 :
1378 0 : end subroutine paw_dijhf
1379 : !!****
1380 :
1381 : !----------------------------------------------------------------------
1382 :
1383 : !!****f* m_paw_slater/summ_2gaunt
1384 : !! NAME
1385 : !! summ_2gaunt
1386 : !!
1387 : !! FUNCTION
1388 : !! Helper function returning \sum_M G_{ij}^{LM} G_{kl}^{LM}
1389 : !!
1390 : !! INPUTS
1391 : !! ij_lm=index of (i_lm,j_lm) element in packed form.
1392 : !! kl_lm=index of (k_lm,l_lm) element in packed form.
1393 : !! ll_idx=Index for L (thus L+1).
1394 : !! Pawang<type(pawang_type)>=paw angular mesh and related data
1395 : !!
1396 : !! SOURCE
1397 :
1398 0 : function summ_2gaunt(Pawang,ij_lm,kl_lm,ll_idx)
1399 :
1400 : !Arguments ------------------------------------
1401 : !scalars
1402 : integer,intent(in) :: ij_lm,kl_lm,ll_idx
1403 : real(dp) :: summ_2gaunt
1404 : type(pawang_type),intent(in) :: Pawang
1405 : !arrays
1406 :
1407 : !Local variables-------------------------------
1408 : !scalars
1409 : integer :: ignt1,ignt2,idx_LM,max_klm,mm,ii,ll
1410 : character(len=500) :: msg
1411 :
1412 : !************************************************************************
1413 :
1414 : ! FIXME: size of gntselect depends on pawxcdev!
1415 : ! Consistency check on input arguments.
1416 0 : max_klm = pawang%l_max**2*(pawang%l_max**2+1)/2
1417 : if (ij_lm>max_klm.or.kl_lm>max_klm.or.ij_lm<1.or.kl_lm<1.or.&
1418 0 : & ll_idx>pawang%l_size_max.or.ll_idx<1) then
1419 0 : write(msg,'(a,3i0)')"Wrong indices, check pawxcdev ",ij_lm,kl_lm,ll_idx
1420 0 : ABI_ERROR(msg)
1421 : end if
1422 :
1423 0 : ll = ll_idx-1
1424 0 : summ_2gaunt=zero; ii=0
1425 0 : do mm=-ll,ll
1426 0 : idx_LM = 1 + ll**2 + ll + mm
1427 0 : ignt1 = Pawang%gntselect(idx_LM,ij_lm)
1428 0 : ignt2 = Pawang%gntselect(idx_LM,kl_lm)
1429 0 : if (ignt1>0 .and. ignt2>0) then
1430 0 : summ_2gaunt = summ_2gaunt + Pawang%realgnt(ignt1)*Pawang%realgnt(ignt2)
1431 0 : ii=ii+1
1432 0 : write(std_out,'(a,4(i2,1x),f8.5,i2)')"ll, mm, ij_lm, kl_lm: ",ll,mm,ij_lm,kl_lm,summ_2gaunt,ii
1433 0 : if (ii/=1) ABI_WARNING("ii>1")
1434 : end if
1435 : end do
1436 :
1437 0 : end function summ_2gaunt
1438 : !!***
1439 :
1440 : !----------------------------------------------------------------------
1441 :
1442 : !!****f* m_paw_slater/slat_intg
1443 : !! NAME
1444 : !! slat_intg
1445 : !!
1446 : !! FUNCTION
1447 : !! Helper function returning the slater integral
1448 : !! \int_\Omega \phi_i(1)\phi_j(1) \dfrac{1}{|1-2|} \phi_k(2)\phi_l(2) d1d2
1449 : !!
1450 : !! INPUTS
1451 : !!
1452 : !! SOURCE
1453 :
1454 0 : function slat_intg(Slatrad4,Pawtab,Pawang,i_lmn,j_lmn,k_lmn,l_lmn)
1455 :
1456 : !Arguments ------------------------------------
1457 : !scalars
1458 : integer,intent(in) :: i_lmn,j_lmn,k_lmn,l_lmn
1459 : real(dp) :: slat_intg
1460 : type(pawtab_type),intent(in) :: Pawtab
1461 : type(pawang_type),intent(in) :: Pawang
1462 : !arrays
1463 : type(slatrad_t),intent(in) :: Slatrad4(:)
1464 :
1465 : !Local variables-------------------------------
1466 : !scalars
1467 : integer :: ij_lmn,kl_lmn,kl_ln,ij_lm,kl_lm,ilsum,ij_ln
1468 : integer :: isel,slt_idx
1469 : integer :: iln,jln,kln,lln,ii
1470 : real(dp) :: sltL_ijkl,angintL_ijkl
1471 : !character(len=500) :: msg
1472 :
1473 : !************************************************************************
1474 :
1475 : ! The lmn packed indices for (ij) and (kl).
1476 0 : if (j_lmn>=i_lmn) then
1477 0 : ij_lmn = i_lmn + j_lmn*(j_lmn-1)/2
1478 : else
1479 0 : ij_lmn = j_lmn + i_lmn*(i_lmn-1)/2
1480 : end if
1481 :
1482 0 : if (l_lmn>=k_lmn) then
1483 0 : kl_lmn = k_lmn + l_lmn*(l_lmn-1)/2
1484 : else
1485 0 : kl_lmn = l_lmn + k_lmn*(k_lmn-1)/2
1486 : end if
1487 : !
1488 : ! The lm indices for (ij) and (kl) in packed storage.
1489 0 : ij_lm = pawtab%indklmn(1,ij_lmn)
1490 0 : ij_ln = pawtab%indklmn(2,ij_lmn)
1491 :
1492 0 : kl_lm = pawtab%indklmn(1,kl_lmn)
1493 0 : kl_ln = pawtab%indklmn(2,kl_lmn)
1494 : !
1495 : ! The index of (ijkl) in the Slatrad4 database.
1496 0 : if (kl_ln>=ij_ln) then
1497 0 : slt_idx = ij_ln +kl_ln*(kl_ln-1)/2
1498 : else
1499 0 : slt_idx = kl_ln +ij_ln*(ij_ln-1)/2
1500 : end if
1501 :
1502 : !BEGIN DEBUG
1503 0 : iln = Slatrad4(slt_idx)%iln
1504 0 : jln = Slatrad4(slt_idx)%jln
1505 0 : kln = Slatrad4(slt_idx)%kln
1506 0 : lln = Slatrad4(slt_idx)%lln
1507 :
1508 0 : ii = kln + lln*(lln-1)/2
1509 0 : if (slt_idx /= (iln + jln*(jln-1)/2 + ii*(ii-1)/2 )) then
1510 0 : write(std_out,*)"slt_idx, iln, jln, kln, lln",slt_idx, iln, jln, kln, lln
1511 0 : ABI_BUG("Check indices")
1512 : end if
1513 : !END DEBUG
1514 : !
1515 : ! Calculate the integral by summing over ilsum.
1516 0 : slat_intg=zero
1517 0 : if (Slatrad4(slt_idx)%nintgl>0) then
1518 0 : do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
1519 : !% do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
1520 0 : isel = Slatrad4(slt_idx)%intgl_select(ilsum)
1521 0 : if (isel/=0) then
1522 0 : sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
1523 0 : angintL_ijkl = summ_2gaunt(Pawang,ij_lm,kl_lm,ilsum)
1524 0 : slat_intg = slat_intg + sltL_ijkl * angintL_ijkl
1525 : end if
1526 : end do
1527 : end if
1528 :
1529 0 : end function slat_intg
1530 : !!***
1531 :
1532 : !----------------------------------------------------------------------
1533 :
1534 0 : END MODULE m_paw_slater
1535 : !!***
|
