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 [ + + - + : 133 : 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 : 3 : integer,allocatable :: gntselect(:,:)
271 : 3 : real(dp),allocatable :: realgnt(:),tmp_sgg(:)
272 : 3 : real(dp),allocatable :: g000(:,:,:)
273 : :
274 : : ! *************************************************************************
275 : :
276 : : !@slatang_cshell_t
277 : :
278 : : ! * Calculate $\Gaunt^{lsl,msl}_{lc,mc;li,mi}$
279 : 3 : lgnt_max = MAX(l_max,lc_max)
280 [ - + + - : 9 : ABI_MALLOC( realgnt,((2*lgnt_max-1)**2*(lgnt_max)**4))
+ - + - -
+ - + ]
281 [ - + + - : 12 : ABI_MALLOC(gntselect,((2*lgnt_max-1)**2, lgnt_max**2*(lgnt_max**2+1)/2))
+ - - + +
- + - + -
- + - + ]
282 : :
283 : 3 : call realgaunt(lgnt_max,ngnt,gntselect,realgnt)
284 : :
285 : : ! Below we need $\Gaunt_{lsl,0}_{lc,0;li,0}$
286 [ + - - + : 18 : ABI_MALLOC(g000,(2*lgnt_max-1,lc_max,l_max))
+ - - + +
- - + + -
+ - + - -
+ - + ]
287 [ + + + + : 70 : g000 = zero
+ + ]
288 : :
289 [ + + ]: 14 : do ilsl=1,2*lgnt_max-1
290 : 11 : lsl = ilsl-1
291 : 11 : ilm0sl = 1+lsl**2+lsl
292 [ + + ]: 41 : do il=1,l_max
293 : 27 : li = il-1
294 : 27 : ilm0 = 1+li**2+li
295 : 27 : k0lm_i = ilm0 *(ilm0-1)/2
296 [ + + ]: 86 : do ilc=1,lc_max
297 : 48 : lc = ilc-1
298 : 48 : ilm0c = 1+lc**2+lc
299 : 48 : k0lm_c= ilm0c * (ilm0c-1)/2
300 [ + + ]: 48 : if (ilm0c > ilm0) then
301 : 8 : klm_ci = k0lm_c + ilm0
302 : : else
303 : 40 : klm_ci = k0lm_i + ilm0c
304 : : end if
305 : 48 : ig000 = gntselect(ilm0sl,klm_ci) ! Index of $\Gaunt_{lsl,0}_{lc,0;li,0}$
306 [ + + ]: 75 : if (ig000 > 0) g000(ilsl,ilc,il)=realgnt(ig000)
307 : : end do
308 : : end do
309 : : end do
310 : :
311 [ - + ]: 3 : ABI_FREE(realgnt)
312 [ - + ]: 3 : ABI_FREE(gntselect)
313 : :
314 : : ! === Loop over klm channels in packed form ===
315 [ + + ]: 68 : do klm=1,lm2_size
316 : 65 : il = klm2lm(1,klm); im = klm2lm(3,klm)
317 : 65 : jl = klm2lm(2,klm); jm = klm2lm(4,klm)
318 : :
319 : 65 : nsggsel=0
320 : 65 : lslat_min = 1 !FIXME find better way
321 : 65 : lslat_max = il+lc_max-1
322 : :
323 : 65 : Slatang3l(klm)%lslat_min = lslat_min
324 : 65 : Slatang3l(klm)%lslat_max = lslat_max
325 : 65 : Slatang3l(klm)%lc_max = lc_max
326 : :
327 [ - + + - : 195 : ABI_MALLOC(tmp_sgg,(lslat_max*lc_max))
+ - + - -
+ - + ]
328 [ + + ]: 415 : tmp_sgg = zero
329 : :
330 : : ! === Calculate F^{lsl,lc}_{li,mi;lj,mj} ===
331 : : ! * Selection rule: mi = mj and li==lj
332 [ + + ]: 65 : if (im == jm .and. il==jl) then
333 : 17 : li = il-1
334 : 17 : ilm = klm2lm(5,klm); k0lm_i = ilm *(ilm -1)/2
335 : 17 : jlm = klm2lm(6,klm); k0lm_j = jlm *(jlm -1)/2
336 : :
337 [ - + + - : 68 : ABI_MALLOC(Slatang3l(klm)%sggselect,(lslat_max,lc_max))
- + + - +
- + - - +
- + - + ]
338 [ + + + + ]: 138 : Slatang3l(klm)%sggselect = 0
339 : :
340 [ + + ]: 66 : do ilsl=lslat_min,lslat_max
341 : : !% do ilsl=lslat_min,lslat_max,2
342 : 49 : lsl = ilsl-1
343 [ + + ]: 157 : do ilc=1,lc_max
344 : 91 : lc = ilc-1
345 : 91 : dum = SQRT( (two*lc+1)*(two*lsl+1) / (four_pi*(two*li+1)) ) * g000(ilsl,ilc,il)
346 [ + + ]: 140 : if (ABS(dum)>=tol12) then ! * Store results and progressive index if non null.
347 : 41 : nsggsel = nsggsel + 1
348 : 41 : tmp_sgg(nsggsel) = dum
349 : 41 : 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 : 65 : Slatang3l(klm)%nsggsel = nsggsel
357 [ + + ]: 65 : if (nsggsel > 0) then
358 [ - + + - : 51 : ABI_MALLOC(Slatang3l(klm)%sgg,(nsggsel))
+ - - + -
+ - + ]
359 [ - + + + ]: 75 : Slatang3l(klm)%sgg = tmp_sgg(1:nsggsel)
360 : : end if
361 [ - + ]: 68 : ABI_FREE(tmp_sgg)
362 : : end do !klm
363 : :
364 [ - + ]: 3 : 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 : 3 : 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 [ - + ]: 3 : 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 [ + + ]: 68 : do ii=1,SIZE(Slatang3l)
417 [ + + ]: 65 : ABI_SFREE(Slatang3l(ii)%sggselect)
418 [ + + ]: 68 : ABI_SFREE(Slatang3l(ii)%sgg)
419 : : end do
420 : :
421 : 3 : 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 [ - + ]: 3 : 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 [ + + ]: 44 : do ii=1,SIZE(Slatrad3l)
450 [ + - ]: 41 : ABI_SFREE(Slatrad3l(ii)%rlphic_select)
451 [ + - ]: 44 : ABI_SFREE(Slatrad3l(ii)%rlphic_int)
452 : : end do
453 : :
454 : 3 : 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 [ + + - + : 85 : 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 : 3 : integer,allocatable :: kln2ln(:,:)
501 : 3 : integer, pointer :: indklmn(:,:),indlmn(:,:)
502 : 3 : real(dp),allocatable :: ff1(:),ff2(:),tmp_integrals(:)
503 : 3 : real(dp),contiguous, pointer :: phi_i(:),phi_j(:)
504 : 3 : real(dp),allocatable,target :: phi_spl(:,:)
505 : 3 : real(dp),allocatable :: der(:),ypp(:)
506 : 3 : real(dp),contiguous, pointer :: crad(:),drad(:),phi_in(:)
507 : :
508 : : ! *************************************************************************
509 : :
510 [ - + ]: 3 : ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
511 [ - + ]: 3 : if (PRESENT(kln_mask)) then
512 : 0 : ABI_ERROR("kln_mask is present")
513 : : end if
514 : :
515 : : !@slatrad_cshell_t
516 : 3 : lmn_size = Pawtab%lmn_size
517 : 3 : lmn2_size = Pawtab%lmn2_size
518 : 3 : ln_size = Pawtab%basis_size
519 : :
520 : 3 : lnc_size = Atm%ln_size
521 : 3 : lc_max = Atm%l_max
522 : :
523 : 3 : call pawrad_isame(Atmrad,Pawrad,hasameq,whichdenser)
524 : :
525 : 3 : do_spline=0
526 [ + + ]: 3 : 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 : : ! Same mesh equation: the core mesh is a prefix of the PAW mesh (it is usually cut at rpaw, see
532 : : ! paw_mkdijexc_core), so the valence partial waves can be used directly on its first points.
533 : : ! Only the partial waves have to cover the core mesh, not the full PAW mesh.
534 [ + - ]: 2 : ABI_CHECK(Atmrad%mesh_size <= Pawtab%partialwave_mesh_size, "Core mesh is longer than the partial-wave mesh!")
535 : : end if
536 : :
537 : 3 : dmesh_size = Atmrad%mesh_size
538 : 3 : cmesh_size = Pawtab%mesh_size
539 : :
540 : 3 : drad => Atmrad%rad(1:dmesh_size)
541 : 3 : crad => Pawrad%rad(1:cmesh_size)
542 : :
543 : : ! === Spline valence basis set onto core mesh (natural spline) ===
544 [ + + ]: 3 : if (do_spline==1) then
545 : 1 : ABI_COMMENT("Splining in init_slatrad3l")
546 [ - + + - : 4 : ABI_MALLOC(phi_spl,(dmesh_size,ln_size))
- + + - +
- + - - +
- + ]
547 [ - + + - : 3 : ABI_MALLOC(der,(cmesh_size))
+ - + - -
+ - + ]
548 [ - + + - : 2 : ABI_MALLOC(ypp,(cmesh_size))
- + + - ]
549 : :
550 [ + + ]: 5 : do iln=1,ln_size
551 : 4 : phi_in => Pawtab%phi(:,iln)
552 [ + + ]: 2484 : ypp(:) = zero; ybcbeg = zero; ybcend = zero
553 : 4 : call spline(crad,phi_in,cmesh_size,ybcbeg,ybcend,ypp)
554 : 5 : call splint(cmesh_size,crad,phi_in,ypp,dmesh_size,drad,phi_spl(:,iln))
555 : : end do
556 : :
557 [ - + ]: 1 : ABI_FREE(der)
558 [ - + ]: 1 : ABI_FREE(ypp)
559 : : end if
560 : :
561 : 3 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
562 : 3 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
563 : :
564 [ - + + - : 9 : ABI_MALLOC(kln2ln,(6,ln2_size))
+ - + - -
+ - + ]
565 : :
566 [ - + - + ]: 3 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
567 : :
568 [ - + + - : 9 : ABI_MALLOC(ff1,(dmesh_size))
+ - + - -
+ - + ]
569 [ - + + - : 6 : ABI_MALLOC(ff2,(dmesh_size))
- + - + ]
570 : :
571 : : ! * Loop over the upper triangle of the [(in,il), (jn,il)] matrix.
572 [ - + ]: 3 : ABI_CHECK(Atm%nsppol==1,"nsppol==2 not tested")
573 : :
574 [ + + ]: 6 : do isppol=1,Atm%nsppol
575 [ + + ]: 47 : do kln=1,ln2_size
576 : 41 : il = kln2ln(1,kln)
577 : 41 : jl = kln2ln(2,kln)
578 : 41 : in = kln2ln(3,kln)
579 : 41 : jn = kln2ln(4,kln)
580 : 41 : iln = kln2ln(5,kln)
581 : 41 : jln = kln2ln(6,kln)
582 : :
583 : 41 : lslat_max = MAX((il+lc_max),(jl+lc_max)) - 1 ! These are indices, not l-values.
584 : : !lslat_min = MIN(ABS(il-lc_max),ABS(jl-lc_max)) + 1
585 : 41 : lslat_min = 1 ! FIXME find better way
586 : :
587 : 41 : Slatrad3l(kln)%lnc_size = lnc_size
588 : 41 : Slatrad3l(kln)%lslat_min = lslat_min
589 : 41 : Slatrad3l(kln)%lslat_max = lslat_max
590 : :
591 : 41 : Slatrad3l(kln)%nrlphic_int = 0
592 : :
593 [ - + + - : 164 : ABI_MALLOC(Slatrad3l(kln)%rlphic_select,(lslat_max,lnc_size))
- + + - +
- + - - +
- + - + ]
594 [ + + + + ]: 639 : Slatrad3l(kln)%rlphic_select(:,:) = 0
595 : :
596 : : !if (PRESENT(kln_mask)) then !FIXME THIS IS WRONG, move it below in case
597 : : ! if (kln_mask(kln)==0) CYCLE
598 : : !end if
599 : :
600 [ + + ]: 41 : if (do_spline==1) then
601 : 10 : ABI_COMMENT("Performing spline of valence phi")
602 : 10 : phi_i => phi_spl(:,iln)
603 : 10 : phi_j => phi_spl(:,jln)
604 : : else
605 : 31 : phi_i => Pawtab%phi(1:dmesh_size,iln)
606 : 31 : phi_j => Pawtab%phi(1:dmesh_size,jln)
607 : : end if
608 : :
609 : : ! * Loop over (n,l) channels for Atom orbitals
610 [ - + + - : 123 : ABI_MALLOC(tmp_integrals,(lslat_max*lnc_size))
+ - + - -
+ - + ]
611 [ + + ]: 494 : tmp_integrals(:) = zero
612 : : nintg=0
613 : :
614 [ + + ]: 186 : do ilnc=1,lnc_size
615 : : ! phicore => Atm%phi(:,ilnc,isppol)
616 [ + - - + : 215310 : ff1 = phi_i * Atm%phi(:,ilnc,isppol)
- - + + ]
617 [ + - - + : 215310 : ff2 = phi_j * Atm%phi(:,ilnc,isppol)
- - + + ]
618 [ + + ]: 639 : do isl=lslat_min,lslat_max ! L coming from Coulomb expansion
619 : 453 : ll = isl-1
620 : 453 : call calc_slatradl(ll,dmesh_size,ff2,ff1,Atmrad,intg1)
621 : 453 : call calc_slatradl(ll,dmesh_size,ff1,ff2,Atmrad,intg)
622 : :
623 : : !call calc_slatradl(ll,cmesh_size,ff2,ff1,Pawrad,intg1)
624 : : !call calc_slatradl(ll,cmesh_size,ff1,ff2,Pawrad,intg)
625 : :
626 [ - + ]: 453 : if (ABS(intg1-intg)>tol6) write(std_out,*)"DEBUG ",ll,il,in,jl,jn,intg1,intg
627 : :
628 : : ! * Store results
629 [ + - ]: 1051 : if (ABS(intg)>=tol12) then
630 : 453 : nintg = nintg +1
631 : 453 : Slatrad3l(kln)%rlphic_select(isl,ilnc) = nintg
632 : 453 : tmp_integrals(nintg) = intg
633 : : end if
634 : : end do !ll
635 : : end do ! ilnc
636 : :
637 : : ! Finalize the object
638 : 41 : Slatrad3l(kln)%nrlphic_int = nintg
639 [ - + + - : 123 : ABI_MALLOC(Slatrad3l(kln)%rlphic_int,(nintg))
+ - + - -
+ - + -
+ ]
640 [ + - + + ]: 535 : if (nintg>0) Slatrad3l(kln)%rlphic_int(:) = tmp_integrals(1:nintg)
641 : :
642 [ - + ]: 44 : ABI_FREE(tmp_integrals)
643 : : end do !kln
644 : : end do !isppol
645 : :
646 [ - + ]: 3 : ABI_FREE(ff1)
647 [ - + ]: 3 : ABI_FREE(ff2)
648 [ - + ]: 3 : ABI_FREE(kln2ln)
649 : :
650 [ + + ]: 3 : if (do_spline==1) then
651 [ - + ]: 1 : ABI_FREE(phi_spl)
652 : : end if
653 : :
654 [ - + ]: 9 : end subroutine slatrad_cshell_init
655 : : !!***
656 : :
657 : : !----------------------------------------------------------------------
658 : :
659 : : !!****f* m_paw_slater/paw_sigxcore
660 : : !! NAME
661 : : !! paw_sigxcore
662 : : !!
663 : : !! FUNCTION
664 : : !! Calculate the integrals:
665 : : !! \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!!
666 : : !! for given (in,il) and (jn,jl) as a function of (nc,lc) and L = |il-jl|, |il-jl|+2, ..., |il+il|
667 : : !!
668 : : !! INPUTS
669 : : !! cplex_dij=1 if dijexc_core is real, 2 if they are complex
670 : : !! lmn2_size=Number of (klmn) channels
671 : : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
672 : : !! Pawtab<pawtab_type>=paw tabulated starting data
673 : : !! Atm<atomorb_type>=Structure containing core orbitals
674 : : !! Atmrad<pawrad_type>=The radial mesh for core orbitals
675 : : !!
676 : : !! OUTPUT
677 : : !! dijexc_core(cplex_dij*lmn2_size,ndij)
678 : : !!
679 : : !! SOURCE
680 : :
681 : 3 : subroutine paw_sigxcore(cplex_dij,lmn2_size,ndij,Pawrad,Pawtab,Atm,Atmrad,dijexc_core)
682 : :
683 : : !Arguments ------------------------------------
684 : : !scalars
685 : : integer,intent(in) :: lmn2_size,cplex_dij,ndij
686 : : !arrays
687 : : real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size,ndij)
688 : : type(atomorb_type),intent(in) :: Atm
689 : : type(pawrad_type),intent(in) :: Atmrad
690 : : type(pawrad_type),intent(in) :: Pawrad
691 : : type(pawtab_type),target,intent(in) :: Pawtab
692 : :
693 : : !Local variables ---------------------------------------
694 : : !scalars
695 : : integer :: ilnc,ilc,lnc_size,l_max
696 : : integer :: lm2_size,ln_size,ln2_size,lmn_size
697 : : integer :: lm_size,klmn,kln,klm
698 : : integer :: lc_max,ilsl,isgg,israd,opt_l,pawprtvol
699 : : real(dp) :: tmp,sgg,intgrl
700 : : !character(len=500) :: msg
701 : : !arrays
702 : : integer :: opt_l_index(0,0),pack2ij(0)
703 : 3 : integer,allocatable :: kln2ln(:,:),klm2lm(:,:)
704 : : integer, pointer :: indklmn(:,:),indlmn(:,:)
705 : : type(slatang_cshell_t),allocatable :: Slatang3l(:)
706 : 3 : type(slatrad_cshell_t),allocatable :: Slatrad3l(:)
707 : :
708 : : ! *************************************************************************
709 : :
710 : : ! * Consistency check
711 [ - + ]: 3 : ABI_CHECK(cplex_dij==1,"cplex_dij must be 1")
712 : :
713 [ - + ]: 3 : ABI_CHECK(ndij==1,"ndij must be 1")
714 : :
715 [ - + ]: 3 : ABI_CHECK(lmn2_size==Pawtab%lmn2_size,"Wrong lmn2_size")
716 : :
717 : 3 : lmn_size = Pawtab%lmn_size
718 : 3 : ln_size = Pawtab%basis_size
719 : 3 : ln2_size = Pawtab%ij_size
720 : 3 : l_max = (Pawtab%l_size-1)/2 +1
721 : 3 : lm_size = l_max**2
722 : 3 : lm2_size = lm_size*(lm_size+1)/2
723 : :
724 : 3 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
725 : 3 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
726 : :
727 : : ! * Setup of useful tables.
728 [ - + + - : 9 : ABI_MALLOC(kln2ln,(6,ln2_size))
+ - + - -
+ - + ]
729 [ - + - + ]: 3 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
730 : :
731 [ - + + - : 9 : ABI_MALLOC(klm2lm,(6,lm2_size))
+ - + - -
+ - + ]
732 [ - + - + ]: 3 : call make_klm2lm(lmn_size,lmn2_size,lm2_size,indlmn,indklmn,klm2lm)
733 : :
734 : : ! * Integrate angular part.
735 : 3 : lnc_size = Atm%ln_size
736 : 3 : lc_max = Atm%l_max
737 : :
738 [ - + + - : 74 : ABI_MALLOC(Slatang3l,(lm2_size))
+ - + - -
+ - + +
+ ]
739 : 3 : call slatang_cshell_init(Slatang3l,l_max,lm2_size,lc_max,klm2lm)
740 : :
741 [ - + ]: 3 : ABI_FREE(klm2lm)
742 : :
743 : : ! * Integrate radial part.
744 [ - + + - : 50 : ABI_MALLOC(Slatrad3l,(ln2_size))
+ - + - -
+ - + +
+ ]
745 : :
746 : 3 : call slatrad_cshell_init(Slatrad3l,ln2_size,Pawrad,Pawtab,Atm,Atmrad)
747 : :
748 : : ! === Calculate matrix elements of Sigma_x^core ===
749 : : ! * $<\phi_i|\Sigma_x^\core|\phi_j>$
750 [ + + + + ]: 249 : dijexc_core = zero
751 [ + + ]: 246 : do klmn=1,lmn2_size
752 : 243 : klm = Pawtab%indklmn(1,klmn)
753 : 243 : kln = Pawtab%indklmn(2,klmn)
754 : : !
755 : : ! * Summing over (lc,nc) and lslat
756 : 243 : tmp = zero
757 [ + + ]: 243 : if (Slatang3l(klm)%nsggsel >0) then
758 [ + + ]: 234 : do ilnc=1,Atm%ln_size
759 : 183 : ilc = 1+Atm%indln(1,ilnc)
760 [ + + ]: 819 : do ilsl=1,Slatang3l(klm)%lslat_max !FIXME check this
761 : : !do ilsl=Slatang3l(klm)%lslat_min,Slatang3l(klm)%lslat_max
762 : 585 : isgg = Slatang3l(klm)%sggselect(ilsl,ilc)
763 : 585 : israd = Slatrad3l(kln)%rlphic_select(ilsl,ilnc)
764 [ + + ]: 768 : if (isgg>0 .and. israd>0) then
765 : 240 : sgg = Slatang3l(klm)%sgg(isgg)
766 : 240 : intgrl = Slatrad3l(kln)%rlphic_int(israd)
767 : 240 : tmp = tmp + intgrl*sgg
768 : : end if
769 : : end do
770 : : end do
771 : : end if
772 : :
773 : 246 : dijexc_core(klmn,1) = -tmp ! Store results.
774 : : end do
775 : :
776 : : ! * Print values
777 : 3 : call wrtout(std_out," ************** Dij Fock_core ************ ",'COLL')
778 : 3 : opt_l=-1; pawprtvol=-1
779 : 3 : 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)
780 : :
781 : : ! * Free memory.
782 [ - + ]: 3 : ABI_FREE(kln2ln)
783 : 3 : call slatang_cshell_free(Slatang3l)
784 [ + - + + : 68 : ABI_FREE(Slatang3l)
- + - + -
+ ]
785 : 3 : call slatrad_cshell_free(Slatrad3l)
786 [ + - + + : 44 : ABI_FREE(Slatrad3l)
- + - + -
+ ]
787 : :
788 : 3 : end subroutine paw_sigxcore
789 : : !!***
790 : :
791 : : !----------------------------------------------------------------------
792 : :
793 : : !!****f* m_paw_slater/paw_mkdijexc_core
794 : : !! NAME
795 : : !! paw_mkdijexc_core
796 : : !!
797 : : !! FUNCTION
798 : : !! Driver routine to calculate the onsite matrix element of the Fock operator between two
799 : : !! all-electron partial waves.
800 : : !!
801 : : !! INPUTS
802 : : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
803 : : !! cplex_dij=1 if dijexc_core is real, 2 if they are complex
804 : : !! lmn2_size_max=Max Number of (klmn) channels over type of atoms.
805 : : !! Cryst<crystal_t>=Structure describing the crystal structure and its symmmetries.
806 : : !! Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
807 : : !! Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
808 : : !! pawprtvol=Flags governing the verbosity of the output.
809 : : !! filpsp(ntypat)=names of the files containing the all-electron core WF
810 : : !!
811 : : !! OUTPUT
812 : : !! dijexc_core(cplex_dij*lmn2_size_max,ndij,ntypat)= On-site matrix elements $ \<\phi_i|Sigma_x^\core|\phi_j\>
813 : : !! for each type of atom.
814 : : !!
815 : : !! SOURCE
816 : :
817 : 2 : subroutine paw_mkdijexc_core(ndij,cplex_dij,lmn2_size_max,Cryst,Pawtab,Pawrad,dijexc_core,pawprtvol,filpsp)
818 : :
819 : : !Arguments ------------------------------------
820 : : !scalars
821 : : integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max
822 : : type(crystal_t),intent(in) :: Cryst
823 : : !arrays
824 : : real(dp),intent(out) :: dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat) !TODO use ragged arrays pawij?
825 : : character(len=fnlen) :: filpsp(Cryst%ntypat)
826 : : type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat)
827 : : type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat)
828 : :
829 : : !Local variables ---------------------------------------
830 : : !scalars
831 : : integer :: itypat,lmn_size,lmn2_size,ln_size,isppol
832 : : real(dp) :: rcut
833 : : character(len=500) :: header
834 : : !arrays
835 : 2 : integer,allocatable :: phi_indln(:,:)
836 : 2 : real(dp),contiguous, pointer :: phi(:,:)
837 : 2 : real(dp),allocatable :: overlap(:,:)
838 : 2 : type(atomorb_type),allocatable :: Atm(:)
839 : 2 : type(Pawrad_type),allocatable :: Radatm(:)
840 : :
841 : : ! *************************************************************************
842 : :
843 [ + - + - : 9 : ABI_MALLOC(Atm,(Cryst%ntypat))
+ - - + -
+ + + ]
844 [ - + + - : 9 : ABI_MALLOC(Radatm,(Cryst%ntypat))
+ - + - -
+ - + +
+ ]
845 : :
846 [ - + ]: 2 : ABI_CHECK(ndij==1 ,"spinor+HF not available")
847 [ - + ]: 2 : ABI_CHECK(cplex_dij==1,"spinor+HF not available")
848 [ + + - + ]: 5 : ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
849 : :
850 : : !allocate(dijexc_core(cplex_dij*lmn2_size_max,ndij,Cryst%ntypat)) !TODO use ragged arrays pawij?
851 [ + + + + : 386 : dijexc_core=zero
+ + ]
852 : :
853 [ + + ]: 5 : do itypat=1,Cryst%ntypat
854 : :
855 : : ! Read core orbitals for this atom type.
856 : 3 : rcut=Pawtab(itypat)%rpaw
857 [ + - ]: 3 : call pawpsp_init_core(Atm(itypat),psp_filename=trim(filpsp(itypat)),rcut_in=rcut)
858 : 3 : call pawrad_copy(Atm(itypat)%radmesh,Radatm(itypat))
859 : :
860 : 3 : write(header,'(a,i4,a)')" === Atom type = ",itypat," === "
861 : 3 : call print_atomorb(Atm(itypat),header,unit=std_out,prtvol=pawprtvol)
862 : : !
863 : : ! * Calculate $ \<\phi_i|Sigma_x^\core|\phi_j\> $ for this atom type.
864 : 3 : lmn_size = Pawtab(itypat)%lmn_size
865 : 3 : lmn2_size = Pawtab(itypat)%lmn2_size
866 : :
867 : : call paw_sigxcore(cplex_dij,lmn2_size,ndij,&
868 [ + + + - : 40 : & Pawrad(itypat),Pawtab(itypat),Atm(itypat),Radatm(itypat),dijexc_core(1:lmn2_size,:,itypat))
+ + + + +
+ ]
869 : :
870 : 3 : ln_size = Pawtab(itypat)%basis_size
871 : 3 : phi => Pawtab(itypat)%phi
872 : :
873 [ - + + - : 9 : ABI_MALLOC(phi_indln,(2,ln_size))
+ - + - -
+ - + ]
874 : 3 : call make_indln(lmn_size,ln_size,Pawtab(itypat)%indlmn(:,:),phi_indln)
875 : :
876 [ - + + - : 12 : ABI_MALLOC(overlap,(Atm(itypat)%ln_size,ln_size))
- + + - +
- + - - +
- + ]
877 : 3 : isppol=1 ! hardcoded
878 : 3 : call get_overlap(Atm(itypat),Radatm(itypat),Pawrad(itypat),isppol,ln_size,phi,phi_indln,overlap)
879 : :
880 [ - + ]: 3 : ABI_FREE(phi_indln)
881 [ - + ]: 5 : ABI_FREE(overlap)
882 : : end do ! ntypat
883 : :
884 : : ! Free memory
885 : 2 : call pawrad_free(Radatm)
886 [ + + ]: 5 : do itypat=1,Cryst%ntypat
887 : 5 : call destroy_atomorb(Atm(itypat))
888 : : end do
889 : :
890 [ + - + + : 5 : ABI_FREE(Atm)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + ]
891 [ + - + + : 5 : ABI_FREE(Radatm)
- + - + -
+ - + ]
892 : :
893 [ - + ]: 2 : end subroutine paw_mkdijexc_core
894 : : !!****
895 : :
896 : : !----------------------------------------------------------------------
897 : :
898 : : !!****f* m_paw_slater/slatrad_free_0D
899 : : !! NAME
900 : : !! slatrad_free_0D
901 : : !!
902 : : !! FUNCTION
903 : : !! Free the dynamic memory allocated in a structure of type slatrad_t
904 : : !!
905 : : !! SOURCE
906 : :
907 : 0 : subroutine slatrad_free_0D(Slatrad)
908 : :
909 : : !Arguments ------------------------------------
910 : : !scalars
911 : : type(slatrad_t),intent(inout) :: Slatrad
912 : :
913 : : ! *********************************************************************
914 : :
915 : : !@slatrad_t
916 [ # # ]: 0 : ABI_SFREE(Slatrad%intgl_select)
917 [ # # ]: 0 : ABI_SFREE(Slatrad%intgl)
918 : :
919 : 0 : end subroutine slatrad_free_0D
920 : : !!***
921 : :
922 : : !----------------------------------------------------------------------
923 : :
924 : : !!****f* m_paw_slater/slatrad_free_1D
925 : : !! NAME
926 : : !! slatrad_free_1D
927 : : !!
928 : : !! FUNCTION
929 : : !! Free the dynamic memory allocated in a structure of type slatrad_t
930 : : !!
931 : : !! SOURCE
932 : :
933 [ # # ]: 0 : subroutine slatrad_free_1D(Slatrad)
934 : :
935 : : !Arguments ------------------------------------
936 : : !scalars
937 : : type(slatrad_t),intent(inout) :: Slatrad(:)
938 : :
939 : : !Local variables-------------------------------
940 : : integer :: ii
941 : : ! *********************************************************************
942 : :
943 [ # # ]: 0 : do ii=1,SIZE(Slatrad)
944 : 0 : call slatrad_free_0D(Slatrad(ii))
945 : : end do
946 : :
947 : 0 : end subroutine slatrad_free_1D
948 : : !!***
949 : :
950 : : !----------------------------------------------------------------------
951 : :
952 : : !!****f* m_paw_slater/slatrad_init
953 : : !! NAME
954 : : !! slatrad_init
955 : : !!
956 : : !! FUNCTION
957 : : !! Initialize the structure storing the radial part of Slater"s integrals.
958 : : !!
959 : : !! INPUTS
960 : : !! which_intg= Option defining what kind of integrals have to be calculated:
961 : : !! -- 1 for Slater integral of AE partial waves only.
962 : : !! 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
963 : : !1
964 : : !! -- 2 for Slater integral of (AE-PS) partial waves
965 : : !! 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 -
966 : : !! \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
967 : : !!
968 : : !! -- 3 for Slater integral of (AE-PS-compensation charges)
969 : : !! C = A -
970 : : !! \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}}*
971 : : !! [tu_k(2) tu_l(2) + qhat^L_\kl r_2^2 g^L(2) ] d1d2
972 : : !!
973 : : !! where u = \phi/r; tu = \tphi/r; and qhat^L_\ij are related to q^\LM_\ij via
974 : : !! q^\LM_\ij = \Gaunt_\ij^\LM qhat^L_\ij => qhat^L\ij = \int (u_i*u_j - tu_i*tu_j) r^L
975 : : !!
976 : : !! ln2_size=Number of symmetrical (l,n) channels for this atom type type.
977 : : !! Pawrad<pawrad_type>=paw radial mesh and related data
978 : : !! Pawtab<pawtab_type>=paw tabulated starting data
979 : : !!
980 : : !! OUTPUT
981 : : !! Slatrad4<slatrad_t>=The object completely initialized.
982 : : !!
983 : : !! NOTES
984 : : !! Slater integrals S_ij are invariant under exchange of the indices,
985 : : !! but the results reported by calc_slatradl are not due to numerical roundoff errors (err < 10^-9).
986 : : !! However this does not cause any problem since only the upper triangle of the S_ij matrix
987 : : !! is stored and used in the other routines.
988 : : !!
989 : : !! SOURCE
990 : :
991 [ # # # # : 0 : subroutine slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad,Pawtab)
# # # # #
# ]
992 : :
993 : : !Arguments ------------------------------------
994 : : !scalars
995 : : integer,intent(in) :: ln2_size,which_intg
996 : : !arrays
997 : : type(pawrad_type),target,intent(in) :: Pawrad
998 : : type(pawtab_type),target,intent(in) :: Pawtab
999 : : type(slatrad_t),intent(out) :: Slatrad4(ln2_size*(ln2_size+1)/2)
1000 : :
1001 : : !Local variables ---------------------------------------
1002 : : !scalars
1003 : : integer :: mesh_size,il,iln,isl,in,jl,jln,jn,sln1,sln2,l_slat
1004 : : integer :: kn,kl,ln,ll,kln,lln,lslat_max,lslat_min,nintgl
1005 : : integer :: lmn_size,lmn2_size,ln_size,slt_idx
1006 : : real(dp) :: ae_intg,ps_intg,pshat_intg,intg,tqij_L,tqkl_L !intg1
1007 : : character(len=500) :: msg
1008 : : !arrays
1009 : 0 : integer,allocatable :: kln2ln(:,:)
1010 : 0 : integer, pointer :: indklmn(:,:),indlmn(:,:)
1011 : 0 : real(dp),allocatable :: uiuj(:),ukul(:),tuituj(:),tuktul(:),tuituj_tqgl(:),tuktul_tqgl(:)
1012 : 0 : real(dp),allocatable :: tmp_integrals(:),ff(:)
1013 : 0 : real(dp),contiguous, pointer :: phi_i(:),phi_j(:),phi_k(:),phi_l(:)
1014 : 0 : real(dp),contiguous, pointer :: tphi_i(:),tphi_j(:),tphi_k(:),tphi_l(:)
1015 : 0 : real(dp),contiguous, pointer :: shapefunc(:),rad(:)
1016 : :
1017 : : ! *************************************************************************
1018 : :
1019 : : DBG_ENTER("COLL")
1020 : :
1021 [ # # ]: 0 : ABI_CHECK(ln2_size==Pawtab%ij_size,"Wrong ln2_size")
1022 : :
1023 [ # # ]: 0 : if ( ALL(which_intg /= (/1,2,3/)) ) then
1024 : 0 : write(msg,'(a,i0)')"Wrong value for which_intg: ",which_intg
1025 : 0 : ABI_ERROR(msg)
1026 : : end if
1027 : :
1028 : : !@slatrad_t
1029 : 0 : lmn_size = Pawtab%lmn_size
1030 : 0 : lmn2_size = Pawtab%lmn2_size
1031 : 0 : ln_size = Pawtab%basis_size
1032 : 0 : mesh_size = Pawtab%mesh_size
1033 : : !
1034 : : ! Useful table for looping.
1035 : 0 : indlmn => Pawtab%indlmn(1:6,1:lmn_size)
1036 : 0 : indklmn => Pawtab%indklmn(1:8,1:lmn2_size)
1037 : :
1038 [ # # # # : 0 : ABI_MALLOC(kln2ln,(6,ln2_size))
# # # # #
# # # ]
1039 [ # # # # ]: 0 : call make_kln2ln(lmn_size,lmn2_size,ln2_size,indlmn,indklmn,kln2ln)
1040 : :
1041 [ # # # # : 0 : ABI_MALLOC(uiuj,(mesh_size))
# # # # #
# # # ]
1042 [ # # # # : 0 : ABI_MALLOC(ukul,(mesh_size))
# # # # ]
1043 [ # # # # : 0 : ABI_MALLOC(ff,(mesh_size))
# # # # ]
1044 [ # # # # : 0 : ABI_MALLOC(tuituj,(mesh_size))
# # # # ]
1045 [ # # # # : 0 : ABI_MALLOC(tuktul,(mesh_size))
# # # # ]
1046 [ # # # # : 0 : ABI_MALLOC(tuituj_tqgl,(mesh_size))
# # # # ]
1047 [ # # # # : 0 : ABI_MALLOC(tuktul_tqgl,(mesh_size))
# # # # ]
1048 : 0 : rad => Pawrad%rad
1049 : : !
1050 : : ! * Loop over (k,l) channels in packed form.
1051 [ # # ]: 0 : do sln2=1,ln2_size
1052 : 0 : kl = kln2ln(1,sln2)
1053 : 0 : ll = kln2ln(2,sln2)
1054 : 0 : kn = kln2ln(3,sln2)
1055 : 0 : ln = kln2ln(4,sln2)
1056 : 0 : kln = kln2ln(5,sln2)
1057 : 0 : lln = kln2ln(6,sln2)
1058 : : !write(std_out,*)"sln2, kln, lln",sln2,kln,lln
1059 : :
1060 : 0 : phi_k => Pawtab%phi (:,kln)
1061 : 0 : tphi_k => Pawtab%tphi(:,kln)
1062 : :
1063 : 0 : phi_l => Pawtab%phi (:,lln)
1064 : 0 : tphi_l => Pawtab%tphi(:,lln)
1065 : : !
1066 : : ! * Loop over (i,j) channels in packed form AND only for the upper triangle of sln2, sln1
1067 [ # # ]: 0 : do sln1=1,sln2
1068 : 0 : il = kln2ln(1,sln1)
1069 : 0 : jl = kln2ln(2,sln1)
1070 : 0 : in = kln2ln(3,sln1)
1071 : 0 : jn = kln2ln(4,sln1)
1072 : 0 : iln = kln2ln(5,sln1)
1073 : 0 : jln = kln2ln(6,sln1)
1074 : : !write(std_out,*)"sln1, iln, jln",sln1,iln,jln
1075 : :
1076 : 0 : slt_idx = sln1 + sln2*(sln2-1)/2 ! index for packed storage.
1077 : :
1078 : 0 : phi_i => Pawtab%phi (:,iln)
1079 : 0 : tphi_i => Pawtab%tphi(:,iln)
1080 : :
1081 : 0 : phi_j => Pawtab%phi (:,jln)
1082 : 0 : tphi_j => Pawtab%tphi(:,jln)
1083 : :
1084 : 0 : lslat_min = MAX(ABS(il-jl),ABS(kl-ll)) + 1 ! We use indices not l-values.
1085 : 0 : lslat_max = MIN((il+jl),(kl+ll)) - 1
1086 : :
1087 : : !lslat_min = MIN(ABS(il-jl),ABS(kl-ll)) + 1
1088 : : !lslat_max = MAX((il+jl),(kl+ll)) - 1
1089 : :
1090 : 0 : Slatrad4(slt_idx)%lslat_min = lslat_min
1091 : 0 : Slatrad4(slt_idx)%lslat_max = lslat_max
1092 : :
1093 : 0 : Slatrad4(slt_idx)%iln = iln
1094 : 0 : Slatrad4(slt_idx)%jln = jln
1095 : 0 : Slatrad4(slt_idx)%kln = kln
1096 : 0 : Slatrad4(slt_idx)%lln = lln
1097 : :
1098 [ # # # # : 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl_select,(lslat_min:lslat_max))
# # # # #
# # # #
# ]
1099 [ # # ]: 0 : Slatrad4(slt_idx)%intgl_select(:) = 0
1100 : 0 : Slatrad4(slt_idx)%nintgl = 0
1101 : :
1102 [ # # ]: 0 : if (lslat_min > lslat_max) then
1103 : : ! e.g. (1 2) (1 1). Due to angular selection rules, this integral do not contribue
1104 : : !write(std_out,*)"lslat_min, lslat_max",lslat_min,lslat_max
1105 : : !write(std_out,*)"il,jl,kl,ll",il,jl,kl,ll
1106 : : !ABI_ERROR("")
1107 [ # # # # ]: 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl,(0))
1108 : 0 : CYCLE
1109 : : end if
1110 : :
1111 [ # # # # : 0 : uiuj = phi_i * phi_j ! The AE part.
# # # # ]
1112 [ # # # # : 0 : ukul = phi_k * phi_l
# # # # ]
1113 [ # # # # : 0 : tuituj = tphi_i * tphi_j ! The pseudized part.
# # # # ]
1114 [ # # # # : 0 : tuktul = tphi_k * tphi_l
# # # # ]
1115 : : !
1116 : : ! Calculate L-depedent integrals where L come from the expansion the Coulomb interaction.
1117 [ # # # # : 0 : ABI_MALLOC(tmp_integrals,(MAX(lslat_min,lslat_max)))
# # # # #
# ]
1118 [ # # ]: 0 : tmp_integrals=zero
1119 : : nintgl=0
1120 : :
1121 [ # # ]: 0 : do isl=lslat_min,lslat_max
1122 : : !do isl=lslat_min,lslat_max,2 ! TODO Here I can reduce the number of iterations using a step of 2.
1123 : 0 : l_slat = isl-1
1124 : 0 : call calc_slatradl(l_slat,mesh_size,uiuj,ukul,Pawrad,ae_intg)
1125 : 0 : intg = ae_intg
1126 : :
1127 : : #if 0
1128 : : call calc_slatradl(l_slat,mesh_size,ukul,uiuj,Pawrad,intg1)
1129 : : if (ABS(intg1-ae_intg)>tol12) then
1130 : : write(msg,'(a,es16.8)')"s_ij and s_ij differ by ",intg1-ae_intg
1131 : : ABI_WARNING(msg)
1132 : : end if
1133 : : #endif
1134 [ # # ]: 0 : if (which_intg == 2) then ! Subtract the pseudo part.
1135 : 0 : call calc_slatradl(l_slat,mesh_size,tuituj,tuktul,Pawrad,ps_intg)
1136 : 0 : intg = intg - ps_intg
1137 : :
1138 [ # # ]: 0 : else if (which_intg == 3) then ! Subtract (pseudo + compensation charges)
1139 : : !
1140 : : ! Evaluate tqij_L and tqkl_L (without M-dependent part).
1141 : 0 : ff(1)=zero
1142 [ # # ]: 0 : ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln1)-pawtab%tphitphj(2:mesh_size,sln1))*rad(2:mesh_size)**l_slat
1143 [ # # # # ]: 0 : if (l_slat==0.and.kl==1.and.ll==1) then
1144 : 0 : call pawrad_deducer0(ff,mesh_size,pawrad)
1145 : : end if
1146 : 0 : call simp_gen(tqij_L,ff,pawrad)
1147 : :
1148 : 0 : ff(1)=zero
1149 [ # # ]: 0 : ff(2:mesh_size)=(pawtab%phiphj(2:mesh_size,sln2)-pawtab%tphitphj(2:mesh_size,sln2))*rad(2:mesh_size)**l_slat
1150 [ # # # # ]: 0 : if (l_slat==0.and.il==1.and.jl==1) then
1151 : 0 : call pawrad_deducer0(ff,mesh_size,pawrad)
1152 : : end if
1153 : 0 : call simp_gen(tqkl_L,ff,pawrad)
1154 : :
1155 : 0 : shapefunc => Pawtab%shapefunc(:,isl) ! Recheck this part, in particular the convention
1156 [ # # # # : 0 : tuituj_tqgl = tuituj + tqij_L * shapefunc * rad**2
# # # # ]
1157 [ # # # # : 0 : tuktul_tqgl = tuktul + tqkl_L * shapefunc * rad**2
# # # # ]
1158 : :
1159 : 0 : call calc_slatradl(l_slat,mesh_size,tuituj_tqgl,tuktul_tqgl,Pawrad,pshat_intg)
1160 : 0 : intg = intg - pshat_intg
1161 : : end if
1162 : : !
1163 : : ! * Store results
1164 [ # # ]: 0 : if (ABS(intg)>=tol12) then
1165 : 0 : nintgl = nintgl +1
1166 : 0 : Slatrad4(slt_idx)%intgl_select(isl) = nintgl
1167 : 0 : tmp_integrals(nintgl) = intg
1168 : : end if
1169 : : end do !isl
1170 : : !
1171 : : ! Finalize the object.
1172 : 0 : Slatrad4(slt_idx)%nintgl = nintgl
1173 [ # # # # : 0 : ABI_MALLOC(Slatrad4(slt_idx)%intgl,(nintgl))
# # # # #
# # # #
# ]
1174 [ # # # # ]: 0 : if (nintgl>0) Slatrad4(slt_idx)%intgl(:) = tmp_integrals(1:nintgl)
1175 [ # # ]: 0 : ABI_FREE(tmp_integrals)
1176 : : end do !sln1
1177 : : end do !sln2
1178 : : !
1179 : : ! Free memory
1180 [ # # ]: 0 : ABI_FREE(kln2ln)
1181 [ # # ]: 0 : ABI_FREE(uiuj)
1182 [ # # ]: 0 : ABI_FREE(ukul)
1183 [ # # ]: 0 : ABI_FREE(tuituj)
1184 [ # # ]: 0 : ABI_FREE(tuktul)
1185 [ # # ]: 0 : ABI_FREE(ff)
1186 [ # # ]: 0 : ABI_FREE(tuituj_tqgl)
1187 [ # # ]: 0 : ABI_FREE(tuktul_tqgl)
1188 : :
1189 : : DBG_EXIT("COLL")
1190 : :
1191 : 0 : end subroutine slatrad_init
1192 : : !!***
1193 : :
1194 : : !----------------------------------------------------------------------
1195 : :
1196 : : !!****f* m_paw_slater/paw_dijhf
1197 : : !! NAME
1198 : : !! paw_dihf
1199 : : !!
1200 : : !! FUNCTION
1201 : : !! This routine calculates the onsite D_{ij} strengths of the exchange part of the self energy.
1202 : : !!
1203 : : !! INPUTS
1204 : : !! ndij=Usually ndij=nspden, except for spin-orbit (where ndij=nspinor**2)
1205 : : !! cplex_dij=1 if sigx_dij is real, 2 if they are complex
1206 : : !! qphase=2 if dij contains a exp(-i.q.r) phase (as in the q<>0 RF case), 1 if not
1207 : : !! lmn2_size_max=Max Number of (klmn) channels over type of atoms.
1208 : : !! my_natom=number of atoms treated by current process
1209 : : !! ntypat=number of atom types
1210 : : !! Pawtab(ntypat)<pawtab_type>=paw tabulated starting data
1211 : : !! Pawrad(ntypat)<pawrad_type>=paw radial mesh and related data
1212 : : !! Pawang<type(pawang_type)>=paw angular mesh and related data
1213 : : !! pawprtvol=Flags governing the verbosity of the output.
1214 : : !!
1215 : : !! OUTPUT
1216 : : !! sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom)=
1217 : : !! For each atom, the Pseudopotential strengths of the on-site operator Sigma_x
1218 : : !!
1219 : : !! NOTES
1220 : : !! The on-site contribution to the matrix elements of the exchange part of the self-energy is given by:
1221 : : !! <\tpsi_a| [\sum_{ij} |tprj_i\> D_{ij} \<tprj_j|] |\tpsi_b\>.
1222 : : !!
1223 : : !! When compensation charges are used one obtains:
1224 : : !!
1225 : : !! D_{ij} = - sum_{kl} \rho_lk [ \Phi_{ikjl} - \Phihat_{ijkl} =
1226 : : !!
1227 : : !! = - sum_{kl} \rho_lk \sum_{LM} \Gaunt_{ik}^{LM} \Gaunt_{jl}^{LM} [S_{ikjl}^L - tS{ikjl}^L}]
1228 : : !!
1229 : : !! where S and tS are Slater-like integrals given by
1230 : : !!
1231 : : !! 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.
1232 : : !! 1) tS_{ijkl}^L = dfrac{4\pi}{2L+1} \iint [u_i(1)u_j(1)+ tq_{ij}^L g^L(1)]
1233 : : !! [u_k(2)u_l(2)+ tq_{kl}^L g^L(2)] \dfrac{r_<^L/}{r_>^{L+1}} d1d2.
1234 : : !!
1235 : : !! tq_{ij}^L is defined in terms of q_{ij}^L via: q_{ij]^{LM} = tq_{ij}^L \Gaunt_{ij}^{LM}
1236 : : !!
1237 : : !! SOURCE
1238 : :
1239 : 0 : subroutine paw_dijhf(ndij,cplex_dij,qphase,lmn2_size_max,my_natom,ntypat,Pawtab,Pawrad,Pawang,Pawrhoij,&
1240 : 0 : & sigx_dij,pawprtvol)
1241 : :
1242 : : !Arguments ------------------------------------
1243 : : !scalars
1244 : : integer,intent(in) :: pawprtvol,ndij,cplex_dij,lmn2_size_max,my_natom,ntypat,qphase
1245 : : type(pawang_type),intent(in) :: Pawang
1246 : : !arrays
1247 : : real(dp),target,intent(out) :: sigx_dij(cplex_dij*lmn2_size_max,ndij,my_natom) !TODO use ragged arrays pawij?
1248 : : type(pawtab_type),intent(in) :: Pawtab(ntypat)
1249 : : type(pawrad_type),intent(in) :: Pawrad(ntypat)
1250 : : type(pawrhoij_type),intent(in) :: Pawrhoij(my_natom)
1251 : :
1252 : : !Local variables ---------------------------------------
1253 : : !scalars
1254 : : integer,parameter :: cplex=1 ! FIXME preliminary implementation
1255 : : integer :: cplex_rhoij,iatom,iq,iq0_dij,iq0_rhoij,itypat,lmn_size,lmn2_size,ispden,nspden,ln2_size
1256 : : integer :: lm2_size !,isppol ln_size,
1257 : : integer :: irhoij,jrhoij
1258 : : integer :: rho_lmn !,rho_klm,rho_kln,rho_lmin,rho_lmax,rho_iln,rho_jln
1259 : : integer :: klmn
1260 : : integer :: i_lmn,j_lmn,k_lmn,l_lmn
1261 : : integer :: which_intg,l_max,opt_l
1262 : : real(dp) :: ro,slt_ikjl,slt_iljk
1263 : : !character(len=500) :: msg
1264 : : !arrays
1265 : : integer :: opt_l_index(0,0),pack2ij(0)
1266 : 0 : real(dp), contiguous, pointer :: sigx_atm(:,:)
1267 : 0 : type(slatrad_t),allocatable :: Slatrad4(:)
1268 : :
1269 : : ! *************************************************************************
1270 : :
1271 : : DBG_ENTER("COLL")
1272 : :
1273 : 0 : ABI_CHECK(ndij/=4,"ndij=4 not coded")
1274 [ # # ]: 0 : ABI_CHECK(cplex_dij==1,"cplex_dij/=2 not coded")
1275 [ # # # # ]: 0 : ABI_CHECK(lmn2_size_max==MAXVAL(Pawtab(:)%lmn2_size),"Wrong lmn2_size_max")
1276 : :
1277 [ # # ]: 0 : if (my_natom>0) then
1278 [ # # ]: 0 : if (pawrhoij(1)%qphase<cplex) then
1279 : 0 : ABI_BUG('Must have pawrhoij()%qphase >= cplex !')
1280 : : end if
1281 : : end if
1282 : :
1283 [ # # # # : 0 : sigx_dij=zero
# # ]
1284 : :
1285 [ # # ]: 0 : do iatom=1,my_natom
1286 : 0 : itypat =Pawrhoij(iatom)%itypat
1287 : 0 : lmn_size =Pawtab(itypat)%lmn_size
1288 : 0 : lmn2_size=Pawtab(itypat)%lmn2_size
1289 : 0 : l_max =(Pawtab(itypat)%l_size+1)/2
1290 : 0 : lm2_size = (l_max**2)*(l_max**2+1)/2
1291 : 0 : cplex_rhoij=Pawrhoij(iatom)%cplex_rhoij
1292 : : !write(std_out,*)"in atom ",iatom,"lm2_size=",lm2_size
1293 : :
1294 : : ! Calculate Slater integral for this atom type.
1295 : : ! TODO obviously these tables should be stored in Pawtab!
1296 : 0 : ln2_size = Pawtab(itypat)%ij_size
1297 [ # # # # : 0 : ABI_MALLOC(Slatrad4,(ln2_size*(ln2_size+1)/2))
# # # # #
# # # # #
# # # # ]
1298 : 0 : which_intg=3
1299 : 0 : call slatrad_init(Slatrad4,which_intg,ln2_size,Pawrad(itypat),Pawtab(itypat))
1300 : :
1301 : 0 : sigx_atm => sigx_dij(:,:,iatom)
1302 : :
1303 : : ! Loop over phase exp(iqr) phase real/imaginary part, if any
1304 [ # # ]: 0 : do iq=1,qphase
1305 : : !First loop: we store the real part in dij(1 -> lmn2_size)
1306 : : !2nd loop: we store the imaginary part in dij(lmn2_size+1 -> 2*lmn2_size)
1307 [ # # ]: 0 : iq0_dij=merge(0,cplex_dij*lmn2_size,iq==1)
1308 [ # # ]: 0 : iq0_rhoij=merge(0,cplex_rhoij*lmn2_size,iq==1)
1309 : :
1310 : : ! * Loop over spin components.
1311 : 0 : nspden=ndij
1312 [ # # ]: 0 : do ispden=1,ndij
1313 : : !
1314 : : ! ============================================================
1315 : : ! ==== Summing over the non-zero lk channels of \rho_{lk} ====
1316 : : ! ============================================================
1317 : 0 : jrhoij=1+iq0_rhoij
1318 [ # # ]: 0 : do irhoij=1,pawrhoij(iatom)%nrhoijsel
1319 : 0 : rho_lmn=pawrhoij(iatom)%rhoijselect(irhoij)
1320 : :
1321 : : ! check wheter rho_lmin is consistent with the Indexing used in slatrad
1322 : : !rho_klm =pawtab(itypat)%indklmn(1,rho_lmn)
1323 : : !rho_kln =pawtab(itypat)%indklmn(2,rho_lmn)
1324 : : !rho_lmin=pawtab(itypat)%indklmn(3,rho_lmn)
1325 : : !rho_lmax=pawtab(itypat)%indklmn(4,rho_lmn)
1326 : :
1327 : : ! Retrieve rhoij for this ispden.
1328 [ # # ]: 0 : if (nspden/=2) then
1329 : 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden)
1330 : : else
1331 : 0 : ABI_ERROR("Recheck this part")
1332 [ # # ]: 0 : if (ispden==1) then
1333 : 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1) + pawrhoij(iatom)%rhoijp(jrhoij,2)
1334 [ # # ]: 0 : else if (ispden==2) then
1335 : 0 : ro=pawrhoij(iatom)%rhoijp(jrhoij,1)
1336 : : end if
1337 : : end if
1338 : : !
1339 : : ! Avoid double-counting the diagonal of rho.
1340 : 0 : ro=ro*pawtab(itypat)%dltij(rho_lmn)*half
1341 : :
1342 : 0 : call klmn2ijlmn(rho_lmn,lmn_size,k_lmn,l_lmn)
1343 : :
1344 : : ! Loop over the upper triangle of the D_{ij) matrix and accumulate:
1345 : : ! sum_\lk rho_\kl [ \Phi_{ikjl} + \Phi_{iljk} - \Phihat_{ikjl} - \Phihat_{iljk} ]
1346 [ # # ]: 0 : do klmn=1,lmn2_size
1347 : : ! Calculate the indices in the Slatrad4 structure.
1348 : 0 : call klmn2ijlmn(klmn,lmn_size,i_lmn,j_lmn)
1349 : :
1350 : : ! My formula
1351 : 0 : slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,j_lmn,l_lmn)
1352 : 0 : slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,j_lmn,k_lmn)
1353 : :
1354 : : !slt_ikjl = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,k_lmn,l_lmn,j_lmn)
1355 : : !slt_iljk = slat_intg(Slatrad4,Pawtab(itypat),Pawang,i_lmn,l_lmn,k_lmn,j_lmn)
1356 : : !slt_iljk = slt_ikjl
1357 : :
1358 : 0 : sigx_atm(klmn+iq0_dij,ispden) = sigx_atm(klmn,ispden) + ro * (slt_ikjl + slt_iljk)
1359 : : end do ! klmn
1360 : :
1361 : 0 : jrhoij=jrhoij+cplex_rhoij
1362 : : end do ! irhoij
1363 : : end do ! iq
1364 : : end do ! ispden
1365 : :
1366 [ # # ]: 0 : if (ABS(pawprtvol)>=1) then ! * Print values
1367 : 0 : call wrtout(std_out," ************** Dij Fock ************ ",'COLL')
1368 : 0 : opt_l=-1
1369 : 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)
1370 : : end if
1371 : :
1372 : 0 : call slatrad_free(Slatrad4)
1373 [ # # # # : 0 : ABI_FREE(Slatrad4)
# # # # #
# ]
1374 : : end do ! iatom
1375 : :
1376 : : ! Factor half cancels in the derivation wrt rho_ij.
1377 [ # # # # : 0 : sigx_dij = - sigx_dij
# # ]
1378 : :
1379 : : DBG_EXIT("COLL")
1380 : :
1381 [ # # # # : 0 : end subroutine paw_dijhf
# # # # #
# ]
1382 : : !!****
1383 : :
1384 : : !----------------------------------------------------------------------
1385 : :
1386 : : !!****f* m_paw_slater/summ_2gaunt
1387 : : !! NAME
1388 : : !! summ_2gaunt
1389 : : !!
1390 : : !! FUNCTION
1391 : : !! Helper function returning \sum_M G_{ij}^{LM} G_{kl}^{LM}
1392 : : !!
1393 : : !! INPUTS
1394 : : !! ij_lm=index of (i_lm,j_lm) element in packed form.
1395 : : !! kl_lm=index of (k_lm,l_lm) element in packed form.
1396 : : !! ll_idx=Index for L (thus L+1).
1397 : : !! Pawang<type(pawang_type)>=paw angular mesh and related data
1398 : : !!
1399 : : !! SOURCE
1400 : :
1401 : 0 : function summ_2gaunt(Pawang,ij_lm,kl_lm,ll_idx)
1402 : :
1403 : : !Arguments ------------------------------------
1404 : : !scalars
1405 : : integer,intent(in) :: ij_lm,kl_lm,ll_idx
1406 : : real(dp) :: summ_2gaunt
1407 : : type(pawang_type),intent(in) :: Pawang
1408 : : !arrays
1409 : :
1410 : : !Local variables-------------------------------
1411 : : !scalars
1412 : : integer :: ignt1,ignt2,idx_LM,max_klm,mm,ii,ll
1413 : : character(len=500) :: msg
1414 : :
1415 : : !************************************************************************
1416 : :
1417 : : ! FIXME: size of gntselect depends on pawxcdev!
1418 : : ! Consistency check on input arguments.
1419 : 0 : max_klm = pawang%l_max**2*(pawang%l_max**2+1)/2
1420 : : if (ij_lm>max_klm.or.kl_lm>max_klm.or.ij_lm<1.or.kl_lm<1.or.&
1421 [ # # # # : 0 : & ll_idx>pawang%l_size_max.or.ll_idx<1) then
# # # # #
# # # ]
1422 : 0 : write(msg,'(a,3i0)')"Wrong indices, check pawxcdev ",ij_lm,kl_lm,ll_idx
1423 : 0 : ABI_ERROR(msg)
1424 : : end if
1425 : :
1426 : 0 : ll = ll_idx-1
1427 : 0 : summ_2gaunt=zero; ii=0
1428 [ # # ]: 0 : do mm=-ll,ll
1429 : 0 : idx_LM = 1 + ll**2 + ll + mm
1430 : 0 : ignt1 = Pawang%gntselect(idx_LM,ij_lm)
1431 : 0 : ignt2 = Pawang%gntselect(idx_LM,kl_lm)
1432 [ # # ]: 0 : if (ignt1>0 .and. ignt2>0) then
1433 : 0 : summ_2gaunt = summ_2gaunt + Pawang%realgnt(ignt1)*Pawang%realgnt(ignt2)
1434 : 0 : ii=ii+1
1435 : 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
1436 [ # # ]: 0 : if (ii/=1) ABI_WARNING("ii>1")
1437 : : end if
1438 : : end do
1439 : :
1440 : 0 : end function summ_2gaunt
1441 : : !!***
1442 : :
1443 : : !----------------------------------------------------------------------
1444 : :
1445 : : !!****f* m_paw_slater/slat_intg
1446 : : !! NAME
1447 : : !! slat_intg
1448 : : !!
1449 : : !! FUNCTION
1450 : : !! Helper function returning the slater integral
1451 : : !! \int_\Omega \phi_i(1)\phi_j(1) \dfrac{1}{|1-2|} \phi_k(2)\phi_l(2) d1d2
1452 : : !!
1453 : : !! INPUTS
1454 : : !!
1455 : : !! SOURCE
1456 : :
1457 [ # # ]: 0 : function slat_intg(Slatrad4,Pawtab,Pawang,i_lmn,j_lmn,k_lmn,l_lmn)
1458 : :
1459 : : !Arguments ------------------------------------
1460 : : !scalars
1461 : : integer,intent(in) :: i_lmn,j_lmn,k_lmn,l_lmn
1462 : : real(dp) :: slat_intg
1463 : : type(pawtab_type),intent(in) :: Pawtab
1464 : : type(pawang_type),intent(in) :: Pawang
1465 : : !arrays
1466 : : type(slatrad_t),intent(in) :: Slatrad4(:)
1467 : :
1468 : : !Local variables-------------------------------
1469 : : !scalars
1470 : : integer :: ij_lmn,kl_lmn,kl_ln,ij_lm,kl_lm,ilsum,ij_ln
1471 : : integer :: isel,slt_idx
1472 : : integer :: iln,jln,kln,lln,ii
1473 : : real(dp) :: sltL_ijkl,angintL_ijkl
1474 : : !character(len=500) :: msg
1475 : :
1476 : : !************************************************************************
1477 : :
1478 : : ! The lmn packed indices for (ij) and (kl).
1479 [ # # ]: 0 : if (j_lmn>=i_lmn) then
1480 : 0 : ij_lmn = i_lmn + j_lmn*(j_lmn-1)/2
1481 : : else
1482 : 0 : ij_lmn = j_lmn + i_lmn*(i_lmn-1)/2
1483 : : end if
1484 : :
1485 [ # # ]: 0 : if (l_lmn>=k_lmn) then
1486 : 0 : kl_lmn = k_lmn + l_lmn*(l_lmn-1)/2
1487 : : else
1488 : 0 : kl_lmn = l_lmn + k_lmn*(k_lmn-1)/2
1489 : : end if
1490 : : !
1491 : : ! The lm indices for (ij) and (kl) in packed storage.
1492 : 0 : ij_lm = pawtab%indklmn(1,ij_lmn)
1493 : 0 : ij_ln = pawtab%indklmn(2,ij_lmn)
1494 : :
1495 : 0 : kl_lm = pawtab%indklmn(1,kl_lmn)
1496 : 0 : kl_ln = pawtab%indklmn(2,kl_lmn)
1497 : : !
1498 : : ! The index of (ijkl) in the Slatrad4 database.
1499 [ # # ]: 0 : if (kl_ln>=ij_ln) then
1500 : 0 : slt_idx = ij_ln +kl_ln*(kl_ln-1)/2
1501 : : else
1502 : 0 : slt_idx = kl_ln +ij_ln*(ij_ln-1)/2
1503 : : end if
1504 : :
1505 : : !BEGIN DEBUG
1506 : 0 : iln = Slatrad4(slt_idx)%iln
1507 : 0 : jln = Slatrad4(slt_idx)%jln
1508 : 0 : kln = Slatrad4(slt_idx)%kln
1509 : 0 : lln = Slatrad4(slt_idx)%lln
1510 : :
1511 : 0 : ii = kln + lln*(lln-1)/2
1512 [ # # ]: 0 : if (slt_idx /= (iln + jln*(jln-1)/2 + ii*(ii-1)/2 )) then
1513 : 0 : write(std_out,*)"slt_idx, iln, jln, kln, lln",slt_idx, iln, jln, kln, lln
1514 : 0 : ABI_BUG("Check indices")
1515 : : end if
1516 : : !END DEBUG
1517 : : !
1518 : : ! Calculate the integral by summing over ilsum.
1519 : 0 : slat_intg=zero
1520 [ # # ]: 0 : if (Slatrad4(slt_idx)%nintgl>0) then
1521 [ # # ]: 0 : do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max
1522 : : !% do ilsum=Slatrad4(slt_idx)%lslat_min,Slatrad4(slt_idx)%lslat_max,2
1523 : 0 : isel = Slatrad4(slt_idx)%intgl_select(ilsum)
1524 [ # # ]: 0 : if (isel/=0) then
1525 : 0 : sltL_ijkl = Slatrad4(slt_idx)%intgl(isel)
1526 : 0 : angintL_ijkl = summ_2gaunt(Pawang,ij_lm,kl_lm,ilsum)
1527 : 0 : slat_intg = slat_intg + sltL_ijkl * angintL_ijkl
1528 : : end if
1529 : : end do
1530 : : end if
1531 : :
1532 : 0 : end function slat_intg
1533 : : !!***
1534 : :
1535 : : !----------------------------------------------------------------------
1536 : :
1537 [ # # # # : 0 : END MODULE m_paw_slater
# # # # #
# # # #
# ]
1538 : : !!***
|
