Branch data 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 : : !!***
|
