Branch data Line data Source code
1 : : !!****m* ABINIT/m_pawpwij
2 : : !! NAME
3 : : !! m_pawpwij
4 : : !!
5 : : !! FUNCTION
6 : : !! This module defines methods to calculate the onsite contribution of a plane wave in the PAW method:
7 : : !!
8 : : !! - pawpwff_t: Form factors used to calculate the onsite contributions of a plane wave.
9 : : !! - pawpwij_t: Onsite matrix elements of a plane wave for a given atom type.
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2008-2026 ABINIT group (MG,GKA)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !!
17 : : !! SOURCE
18 : :
19 : : #if defined HAVE_CONFIG_H
20 : : #include "config.h"
21 : : #endif
22 : :
23 : : #include "abi_common.h"
24 : :
25 : : MODULE m_pawpwij
26 : :
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_fft
31 : :
32 : : use m_fstrings, only : sjoin, itoa
33 : : use defs_datatypes, only : pseudopotential_type
34 : : use defs_abitypes, only : MPI_type
35 : : use m_numeric_tools, only : arth
36 : : use m_geometry, only : metric
37 : : use m_crystal, only : crystal_t
38 : : use m_paw_numeric, only : paw_jbessel_4spline, paw_spline
39 : : use m_splines, only : splfit
40 : : use m_pawang, only : pawang_type
41 : : use m_paw_sphharm, only : realgaunt
42 : : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free, pawrad_copy, simp_gen
43 : : use m_pawtab, only : pawtab_type
44 : : use m_pawcprj, only : pawcprj_type
45 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
46 : : use m_initylmg, only : initylmg
47 : :
48 : : implicit none
49 : :
50 : : private
51 : : !!***
52 : :
53 : : !!****t* m_pawpwij/pawpwff_t
54 : : !! NAME
55 : : !! pawpwff_t
56 : : !!
57 : : !! FUNCTION
58 : : !! PAW form factors used to evaluate $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
59 : : !!
60 : : !! SOURCE
61 : :
62 : : type,public :: pawpwff_t
63 : :
64 : : integer :: method = -1
65 : : ! 1 For Arnaud-Alouani"s exact expression.
66 : : ! 2 For Shishkin-Kresse"s approximated expression.
67 : :
68 : : integer :: dim1 = -1, dim2 = -1
69 : : ! Dimensions of pwff_spl, depending on method.
70 : :
71 : : integer :: nq_spl = -1
72 : : ! Number of points in the reciprocal space grid on which
73 : : ! the radial integrals are evaluated.
74 : :
75 : : real(dp) :: gmet(3,3) = zero
76 : : ! Reciprocal space metric tensor in Bohr**-2
77 : :
78 : : real(dp),allocatable :: qgrid_spl(:)
79 : : ! qgrid_spl(nq_spl)
80 : : ! The coordinates of the points of the radial grid for the integrals used in the spline.
81 : :
82 : : real(dp),allocatable :: pwff_spl(:,:,:,:)
83 : : ! pwff_spl(nq_spl,2,0:dim1,dim2)
84 : : ! The different integrals on the radial |q| grid, for a given atom type.
85 : :
86 : : end type pawpwff_t
87 : :
88 : : public :: pawpwff_init ! Initialize form factors for spline.
89 : : public :: pawpwff_free ! Deallocate dynamic memory.
90 : : !!***
91 : :
92 : : !----------------------------------------------------------------------
93 : :
94 : : !!****t* m_pawpwij/pawpwij_t
95 : : !! NAME
96 : : !! pawpwij_t
97 : : !!
98 : : !! FUNCTION
99 : : !! For PAW, object storing $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
100 : : !! for a given q-point, for a particular TYPE of atom. Therefore the phase factor
101 : : !! e^{i(q+G).R_at} has to be considered to have the onsite contribution of a particular atom.
102 : : !!
103 : : !! SOURCE
104 : :
105 : : type,public :: pawpwij_t
106 : :
107 : : integer :: istpw
108 : : ! Storage mode (similar to istwfk), not used at present
109 : :
110 : : integer :: npw
111 : : ! The number of plane waves
112 : :
113 : : integer :: lmn_size
114 : : ! Number of (l,m,n) elements for the paw basis
115 : :
116 : : integer :: lmn2_size
117 : : ! lmn2_size=lmn_size*(lmn_size+1)/2
118 : : ! where lmn_size is the number of (l,m,n) elements for the paw basis
119 : :
120 : : real(dp) :: qpt(3)
121 : : ! The q-point in e^{-i(q+G)}.r}
122 : :
123 : : real(dp),allocatable :: mqpgij(:,:,:)
124 : : ! pwij(2,npw,lmn2_size)
125 : : ! $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
126 : :
127 : : end type pawpwij_t
128 : :
129 : : public :: pawpwij_init ! Calculate onsite matrix elements of a set of plane waves.
130 : : public :: pawpwij_free ! Deallocate dynamic memory in the structure.
131 : : public :: paw_rho_tw_g ! Calculate the PAW contribution to the oscillator matrix element.
132 : : public :: paw_cross_rho_tw_g ! Calculate the PAW cross term contribution to the oscillator matrix element.
133 : : !!***
134 : :
135 : : interface pawpwij_free
136 : : module procedure pawpwij_free_d1
137 : : module procedure pawpwij_free_d2
138 : : end interface pawpwij_free
139 : :
140 : : !----------------------------------------------------------------------
141 : :
142 : : integer,parameter :: PWIJ_ARNAUD = 1 ! Arnaud-Alouani exact expression. PRB 62. 4464 [[cite:Arnaud2000]]
143 : : integer,parameter :: PWIJ_SHISHKIN = 2 ! Shishkin-Kresse approximated expression. PRB 74. 035101 [[cite:Shishkin2006]]
144 : :
145 : : CONTAINS !========================================================================================
146 : : !!***
147 : :
148 : : !!****f* m_pawpwij/pawpwff_init
149 : : !! NAME
150 : : !! pawpwff_init
151 : : !!
152 : : !! FUNCTION
153 : : !! Initialize the structure containing integrals used to evaluate the onsite
154 : : !! matrix elements of a plane wave by means of a spline fit technique.
155 : : !!
156 : : !! INPUTS
157 : : !! method=1 for Arnaud-Alouani, 2 for Shishkin-Kresse.
158 : : !! nq_spl(%ntypat)=Number of points in the mesh used for the spline.
159 : : !! qmax(%ntypat)=Max |q| for the mesh
160 : : !! gmet(3,3)=reciprocal space metric tensor in bohr**-2.
161 : : !! Pawrad(%ntypat) <type(pawrad_type)>=paw radial mesh and related data....
162 : : !! Pawtab(%ntypat) <type(pawtab_type)>=paw tabulated starting data.
163 : : !! %lsize=1+maximum value of l leading to non zero Gaunt coeffs.
164 : : !! %ij_size=Number of (i,j) elements for the symetric paw basis
165 : : !! %lmn2_size=lmn_size*(lmn_size+1)/2
166 : : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
167 : : !! %ntypat=Number of type of atoms.
168 : : !!
169 : : !! OUTPUT
170 : : !! Paw_pwff(%ntypat) <pawpwff_t>=Object storing the form factors
171 : : !! for the spline used in pawpwij_init.
172 : : !! SOURCE
173 : :
174 [ + + - + : 247 : subroutine pawpwff_init(Paw_pwff,method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
- + + + +
+ + + +
- ]
175 : :
176 : : !Arguments ------------------------------------
177 : : !scalars
178 : : integer,intent(in) :: method
179 : : type(Pseudopotential_type),intent(in) :: Psps
180 : : !arrays
181 : : integer,intent(in) :: nq_spl(Psps%ntypat)
182 : : real(dp),intent(in) :: gmet(3,3)
183 : : real(dp),intent(in) :: qmax(Psps%ntypat)
184 : : type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat)
185 : : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
186 : : type(pawpwff_t),intent(out) :: Paw_pwff(Psps%ntypat)
187 : :
188 : : !Local variables-------------------------------
189 : : !scalars
190 : : integer :: dim1,dim2,itypat,nq
191 : : real(dp) :: dq
192 : :
193 : : !************************************************************************
194 : :
195 : : !@pawpwff_t
196 : :
197 : : ! === Evaluate form factors for the radial part of phi.phj-tphi.tphj ===
198 [ + + ]: 41 : do itypat=1,Psps%ntypat
199 : 26 : Paw_pwff(itypat)%method = method
200 : :
201 : 11 : select case (method)
202 : : case (PWIJ_ARNAUD)
203 : 11 : dim1 = Pawtab(itypat)%l_size-1
204 : 11 : dim2 = Pawtab(itypat)%ij_size
205 : : case (PWIJ_SHISHKIN)
206 : 15 : dim1 = Pawtab(itypat)%l_size**2
207 : 15 : dim2 = Pawtab(itypat)%lmn2_size
208 : : case default
209 [ + + - ]: 26 : ABI_BUG(sjoin("Wrong method:", itoa(method)))
210 : : end select
211 : :
212 : 26 : Paw_pwff(itypat)%dim1 = dim1
213 : 26 : Paw_pwff(itypat)%dim2 = dim2
214 [ + + + + ]: 338 : Paw_pwff(itypat)%gmet = gmet
215 : :
216 : : ! Setup of the q-mesh for spline. It can be type-dependent.
217 : 26 : nq = nq_spl(itypat)
218 : 26 : dq = qmax(itypat)/(one*(nq-1))
219 : : !write(std_out,*)"nq,dq",nq,dq
220 : :
221 : 26 : Paw_pwff(itypat)%nq_spl = nq
222 [ - + + - : 78 : ABI_MALLOC(Paw_pwff(itypat)%qgrid_spl,(nq))
+ - + - -
+ - + -
+ ]
223 [ + - + - : 78078 : Paw_pwff(itypat)%qgrid_spl = arth(zero,dq,nq)
- + - - +
+ ]
224 : : !
225 : : ! === Calculate form factors depending on method ===
226 [ - + + - : 130 : ABI_MALLOC(Paw_pwff(itypat)%pwff_spl,(nq,2,0:dim1,dim2))
- + + - -
+ + - + -
+ - - + -
+ - + ]
227 : :
228 : : call paw_mkrhox_spl(itypat,Psps%ntypat,method,dim1,dim2,nq, &
229 : 41 : Paw_pwff(itypat)%qgrid_spl,Pawrad,Pawtab,Paw_pwff(itypat)%pwff_spl)
230 : : end do ! itypat
231 : :
232 : 15 : end subroutine pawpwff_init
233 : : !!***
234 : :
235 : : !----------------------------------------------------------------------
236 : :
237 : : !!****f* m_pawpwij/pawpwff_free
238 : : !! NAME
239 : : !! pawpwff_free
240 : : !!
241 : : !! FUNCTION
242 : : !! Free memory
243 : : !!
244 : : !! SOURCE
245 : :
246 [ - + ]: 15 : subroutine pawpwff_free(Paw_pwff)
247 : :
248 : : !Arguments ------------------------------------
249 : : !scalars
250 : : type(pawpwff_t),intent(inout) :: Paw_pwff(:)
251 : :
252 : : !Local variables-------------------------------
253 : : integer :: ii
254 : :
255 : : !************************************************************************
256 : :
257 [ + + ]: 41 : do ii=1,SIZE(Paw_pwff)
258 [ + - ]: 26 : ABI_SFREE(Paw_pwff(ii)%qgrid_spl)
259 [ + - ]: 41 : ABI_SFREE(Paw_pwff(ii)%pwff_spl)
260 : : end do
261 : :
262 : 15 : end subroutine pawpwff_free
263 : : !!***
264 : :
265 : : !----------------------------------------------------------------------
266 : :
267 : : !!****f* m_pawrhox/pawpwij_init
268 : : !! NAME
269 : : !! pawpwij_init
270 : : !!
271 : : !! FUNCTION
272 : : !! Calculate the onsite matrix elements
273 : : !!
274 : : !! $ <phj|e^{-i(q+G)}|phi> - <tphj|e^{-i(q+G)}|tphi> $
275 : : !!
276 : : !! for a given q and a set of g in gvec for a given TYPE of atom.
277 : : !! Phase factors arising from atom positions are therefore not included.
278 : : !!
279 : : !! INPUTS
280 : : !! npw=Number of plane waves
281 : : !! Psps <type(pseudopotential_type)>=variables related to pseudopotentials
282 : : ! gvec(3,npw)=Plane wave reduced components.
283 : : !! qpt_in(3)=The reduced components of the q-point.
284 : : !! rprim(3,3)=dimensionless real space primitive translations
285 : : !! Pawtab(%ntypat) <type(pawtab_type)>=paw tabulated starting data
286 : : !! Paw_pwff(%ntypat) <pawpwff_t>=Object storing the form factors for the spline used in pawpwij_init.
287 : : !! Psps<type(pseudopotential_type)>=variables related to pseudopotentials
288 : : !!
289 : : !! OUTPUT
290 : : !! Pwij(%ntypat)<pawpwij_t>=Structure containing the onsite matrix elements of e^{-i(q+G).r}.
291 : : !! Completely initialized in output.
292 : : !!
293 : : !! SOURCE
294 : :
295 [ + + - + : 4477 : subroutine pawpwij_init(Pwij, npw, qpt_in, gvec, rprimd, Psps, Pawtab, Paw_pwff)
+ + + - ]
296 : :
297 : : !Arguments ------------------------------------
298 : : !scalars
299 : : integer,intent(in) :: npw
300 : : type(Pseudopotential_type),intent(in) :: Psps
301 : : !arrays
302 : : integer,intent(in) :: gvec(3,npw)
303 : : real(dp),intent(in) :: qpt_in(3),rprimd(3,3)
304 : : type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
305 : : type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat)
306 : : type(pawpwij_t),intent(out) :: Pwij(Psps%ntypat)
307 : :
308 : : !Local variables-------------------------------
309 : : !scalars
310 : : integer,parameter :: unkg0=0,unylm0=0
311 : : integer :: dim1,dim2,method, my_mqmem,my_nqpt,optder,two_lmaxp1,itypat
312 : : integer :: dummy_nsppol,lmn_size,lmn2_size,nq_spl,ierr
313 : : real(dp) :: ucvol
314 : 1449 : type(MPI_type) :: MPI_enreg_seq
315 : : !arrays
316 : 1449 : integer,allocatable :: npwarr(:),dummy_nband(:)
317 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
318 : 1449 : real(dp),allocatable :: my_qtmp(:,:), ylm_q(:,:),ylmgr_q(:,:,:)
319 : :
320 : : ! *********************************************************************
321 : :
322 : : ! ===============================================
323 : : ! === Get real spherical harmonics in G space ===
324 : : ! ===============================================
325 : 1449 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
326 : :
327 : : ! Set up of REAL Ylm(q+G) up to 2*l_max for this q-point.
328 : 1449 : my_mqmem=1; two_lmaxp1=2*Psps%mpsang-1; optder=0
329 : :
330 [ - + + - : 5796 : ABI_MALLOC(ylm_q ,(npw*my_mqmem,two_lmaxp1**2))
- + + - +
- + - - +
- + ]
331 [ - + + - : 7245 : ABI_MALLOC(ylmgr_q,(npw*my_mqmem,3+6*(optder/2),two_lmaxp1**2))
+ - - + +
- + - + -
- + - + ]
332 : :
333 : 1449 : my_nqpt=1
334 [ + - ]: 1449 : ABI_MALLOC(my_qtmp,(3,my_nqpt))
335 [ + + ]: 5796 : my_qtmp(:,1)=qpt_in(:)
336 : :
337 : : ! dummy_nband and dummy_nsppol are not used in sequential mode.
338 : 1449 : dummy_nsppol=1
339 [ + - ]: 1449 : ABI_MALLOC(dummy_nband,(my_nqpt*dummy_nsppol))
340 [ + + ]: 2898 : dummy_nband=0
341 [ + - ]: 1449 : ABI_MALLOC(npwarr,(my_nqpt))
342 [ + + ]: 2898 : npwarr(:)=npw
343 : :
344 : : ! Fake MPI_type for sequential part.
345 : 1449 : call initmpi_seq(MPI_enreg_seq)
346 : :
347 : : call initylmg(gprimd,gvec,my_qtmp,my_mqmem,MPI_enreg_seq,two_lmaxp1,npw,dummy_nband,my_nqpt,&
348 : 1449 : npwarr,dummy_nsppol,optder,rprimd,ylm_q,ylmgr_q)
349 : :
350 : 1449 : call destroy_mpi_enreg(MPI_enreg_seq)
351 : :
352 [ - + ]: 1449 : ABI_FREE(my_qtmp)
353 [ - + ]: 1449 : ABI_FREE(dummy_nband)
354 [ - + ]: 1449 : ABI_FREE(npwarr)
355 [ - + ]: 1449 : ABI_FREE(ylmgr_q)
356 : :
357 : : ! Construct the Pwij structure
358 [ + + ]: 2963 : do itypat=1,Psps%ntypat
359 : :
360 : 1514 : Pwij(itypat)%istpw = 1
361 : 1514 : Pwij(itypat)%npw = npw
362 : :
363 : 1514 : lmn_size = Pawtab(itypat)%lmn_size
364 : 1514 : lmn2_size = lmn_size*(lmn_size+1)/2
365 : 1514 : Pwij(itypat)%lmn_size = lmn_size
366 : 1514 : Pwij(itypat)%lmn2_size = lmn2_size
367 : :
368 [ + + ]: 6056 : Pwij(itypat)%qpt(:) = qpt_in(:)
369 : :
370 : : ! Prepare the call to paw_mkrhox
371 : 1514 : method = Paw_pwff(itypat)%method
372 : 1514 : dim1 = Paw_pwff(itypat)%dim1
373 : 1514 : dim2 = Paw_pwff(itypat)%dim2
374 : 1514 : nq_spl = Paw_pwff(itypat)%nq_spl
375 [ + + + + ]: 19682 : gmet = Paw_pwff(itypat)%gmet
376 : :
377 [ - + + - : 6056 : ABI_MALLOC_OR_DIE(Pwij(itypat)%mqpgij,(2,npw,lmn2_size), ierr)
- + + - +
- + - + -
+ - - + +
- - + ]
378 : :
379 : : ! Evaluate oscillator matrix elements mqpgij
380 : : call paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,Paw_pwff(itypat)%qgrid_spl,Paw_pwff(itypat)%pwff_spl,&
381 : 2963 : gmet,qpt_in,npw,gvec,ylm_q,Psps,Pawtab,Pwij(itypat)%mqpgij)
382 : : end do ! itypat
383 : :
384 [ - + ]: 1449 : ABI_FREE(ylm_q)
385 : :
386 [ - + - + : 1449 : end subroutine pawpwij_init
- + - + -
+ - + - +
- + - + -
+ - + ]
387 : : !!***
388 : :
389 : : !----------------------------------------------------------------------
390 : :
391 : : !!****f* m_pawpwij/pawpwij_free_d1
392 : : !! NAME
393 : : !! pawpwij_free_d1
394 : : !!
395 : : !! FUNCTION
396 : : !! Free all memory allocated in a structure of type pawpwij_t
397 : : !!
398 : : !! SOURCE
399 : :
400 [ - + ]: 1449 : subroutine pawpwij_free_d1(Pwij)
401 : :
402 : : !Arguments ------------------------------------
403 : : type(pawpwij_t),intent(inout) :: Pwij(:)
404 : :
405 : : !Local variables-------------------------------
406 : : integer :: ii
407 : : !************************************************************************
408 : :
409 [ + + ]: 2963 : do ii=1,SIZE(Pwij)
410 [ + - ]: 2963 : ABI_SFREE(Pwij(ii)%mqpgij)
411 : : end do
412 : :
413 : 1449 : end subroutine pawpwij_free_d1
414 : : !!***
415 : :
416 : : !----------------------------------------------------------------------
417 : :
418 : : !!****f* m_pawpwij/pawpwij_free_d2
419 : : !! NAME
420 : : !! pawpwij_free_d2
421 : : !!
422 : : !! FUNCTION
423 : : !! Free all memory allocated in a structure of type pawpwij_t
424 : : !!
425 : : !! SOURCE
426 : :
427 [ # # ]: 0 : subroutine pawpwij_free_d2(Pwij)
428 : :
429 : : !Arguments ------------------------------------
430 : : !scalars
431 : : type(pawpwij_t),intent(inout) :: Pwij(:,:)
432 : :
433 : : !Local variables-------------------------------
434 : : integer :: jj
435 : :
436 : : !************************************************************************
437 : :
438 [ # # ]: 0 : do jj=1,SIZE(Pwij,DIM=2)
439 : 0 : call pawpwij_free_d1(Pwij(:,jj))
440 : : end do
441 : :
442 : 0 : end subroutine pawpwij_free_d2
443 : : !!***
444 : :
445 : : !----------------------------------------------------------------------
446 : :
447 : : !!****f* m_pawpwij/paw_mkrhox_spl
448 : : !! NAME
449 : : !! paw_mkrhox_spl
450 : : !!
451 : : !! FUNCTION
452 : : !! Evaluate PAW form factor ff^{aL}_{ij}(q) for each angular momentum L,
453 : : !! each type of atom, a, and each [(in,il),(jn,jl)] channel. These quantities
454 : : !! are used in paw_mkrhox to evaluate $<phi|e^{-i(q+G)}|phj>-<tphi|e^{-i(q+G)}|tphj>$
455 : : !! for an arbitrary q+G vector.
456 : : !!
457 : : !! INPUTS
458 : : !! dim1
459 : : !! = 2*(Pawtab(itypat)%l_size-1) if method=1
460 : : !! = MAXVAL(Pawtab(:)%l_size)**2 if method=2
461 : : !! dim2
462 : : !! = Pawtab(itypat)%ij_size if method=1
463 : : !! = MAXVAL(Pawtab(:)%lmn2_size) if method=2
464 : : !! method=integer flag defining the approach used:
465 : : !! 1 --> Expression based on the expansion on a plane wave in terms of Bessel functions
466 : : !! and spherical harmonics (Arnaud-Alouani's methos, see PRB 62, 4464 [[cite:Arnaud2000]]
467 : : !! 2 --> Approximate expression with correct description of the multipoles. Eq. 9 in PRB 74, 035101 [[cite:Shishkin2006]]
468 : : !! nq_spl=number of grid points in the q-mesh
469 : : !! qgrid_spl(nq_spl)=values where form factors are returned
470 : : !! ntypat=number of type of atoms
471 : : !! Pawrad<type(Pawrad_type)>=datatype containing radial grid information
472 : : !! Pawtab(ntypat)<type(pawtab_type)>=PAW tabulated starting data
473 : : !!
474 : : !! OUTPUT
475 : : !! pwff_spl(nq_spl,2,0:dim1,dim1_rhox2,ntypat)
476 : : !! form factors ff^{aL}_{ij}(q) and second derivative in packed storage mode
477 : : !! === if method=1 ===
478 : : !! $ff_^{aL}_{ij}(q) =
479 : : !! \int_0^{r_a} j_L(2\pi qr) [phi_{n_i,l_i}.phi_{n_j l_j}(r) - tphi_{n_i l_i}.tph_{n_j l_j}(r)] dr$
480 : : !! === if method=2 ===
481 : : !! $ff_^{aL}_{ij}(q) = q_ij^{LM} \int_0^{r_a} j_L(2\pi qr) g_L(r) r^2 dr$
482 : : !!
483 : : !! NOTES
484 : : !! * $j_L(2\pi q)$ is a spherical Bessel function
485 : : !! * Output matrix elements are stored in packed storage mode
486 : : !! * Inspired by pawpsp_nl
487 : : !!
488 : : !! TODO
489 : : !! One might save CPU time taking into account Gaunt selection rules!
490 : : !!
491 : : !! SOURCE
492 : :
493 : 26 : subroutine paw_mkrhox_spl(itypat,ntypat,method,dim1,dim2,nq_spl,qgrid_spl,Pawrad,Pawtab,pwff_spl)
494 : :
495 : : !Arguments ------------------------------------
496 : : !scalars
497 : : integer,intent(in) :: itypat,ntypat,method,dim2,dim1,nq_spl
498 : : !arrays
499 : : real(dp),intent(in) :: qgrid_spl(nq_spl)
500 : : real(dp),intent(out) :: pwff_spl(nq_spl,2,0:dim1,dim2)
501 : : type(Pawrad_type),intent(in) :: Pawrad(ntypat)
502 : : type(Pawtab_type),intent(in) :: Pawtab(ntypat)
503 : :
504 : : !Local variables-------------------------------
505 : : !scalars
506 : : integer :: mm,nlmn,jlmn,ilmn,klmn,ij_size,l_size !,ider
507 : : integer :: iq,ir,ll,meshsz,mmax,iln,jln,nln,k0ln,kln,qlm
508 : : real(dp),parameter :: EPS=tol14**4,TOLJ=0.001_dp
509 : : real(dp) :: arg,argn,bes,besp,qr,yp1,ypn
510 : : character(len=500) :: msg
511 : 26 : type(Pawrad_type) :: Tmpmesh
512 : : !arrays
513 : 26 : real(dp),allocatable :: rrshape_l(:),shape_l(:),ff(:),gg(:),rr(:),rrdphi_ij(:)
514 : 26 : real(dp),allocatable :: dphi_ij(:),tmp_spl(:,:,:,:),tmp_jgl(:,:,:)
515 : :
516 : : !*************************************************************************
517 : :
518 : : DBG_ENTER("COLL")
519 : :
520 [ + + + + : 143634831 : pwff_spl=zero
+ + + + ]
521 : :
522 : 37 : SELECT CASE (method)
523 : :
524 : : CASE (PWIJ_ARNAUD)
525 : : ! === Arnaud-Alouani exact expression PRB 62. 4464 [[cite:Arnaud2000]] ===
526 : : ! * $ff_^{aL}_{ij}(q) =
527 : : ! \int_0^{r_a} j_L(2\pi qr) [phi_{n_i,l_i}.phi_{n_j l_j}(r)-tphi_{n_i l_i}.tph_{n_j l_j}(r)]dr$
528 : : ! * It does not descrive correctly the multipoles of the AE charge density if low cutoff on G
529 : 11 : write(msg,'(a,i3)')' paw_mkrhox_spl: Using Arnaud-Alouani expression for atom type: ',itypat
530 : 11 : call wrtout(std_out,msg)
531 : :
532 : 11 : nln = Pawtab(itypat)%basis_size
533 : 11 : ij_size = Pawtab(itypat)%ij_size
534 : 11 : l_size = Pawtab(itypat)%l_size
535 : :
536 [ - + + - : 55 : ABI_MALLOC(tmp_spl,(nq_spl,2,0:l_size-1,ij_size))
- + + - -
+ + - + -
+ - - + -
+ ]
537 : :
538 : : ! Is mesh beginning with r=0 ?
539 [ - + ]: 11 : if (abs(Pawrad(itypat)%rad(1)) > tol10) then
540 : 0 : ABI_ERROR("Radial mesh starts with r/=0")
541 : : end if
542 : : !
543 : : ! === Initialize temporary arrays and variables ===
544 : 11 : meshsz = Pawtab(itypat)%mesh_size ; mmax=meshsz
545 : :
546 [ - + + - : 33 : ABI_MALLOC(dphi_ij,(meshsz))
+ - + - -
+ - + ]
547 [ - + + - : 22 : ABI_MALLOC(rrdphi_ij,(meshsz))
- + - + ]
548 [ - + + - : 22 : ABI_MALLOC(ff,(meshsz))
- + - + ]
549 [ - + + - : 8003 : ABI_CALLOC(gg,(meshsz))
- + - + +
+ ]
550 [ - + + - : 22 : ABI_MALLOC(rr,(meshsz))
- + - + ]
551 [ + + ]: 7992 : rr(1:meshsz) = Pawrad(itypat)%rad(1:meshsz)
552 : : !
553 : : ! === Loop on (jln,iln) channels for this type. Packed form ===
554 [ + + + + : 2720397 : tmp_spl=zero
+ + + + ]
555 : :
556 [ + + ]: 57 : do jln=1,nln
557 : 46 : k0ln=jln*(jln-1)/2
558 [ + + ]: 178 : do iln=1,jln
559 : 121 : kln=k0ln+iln
560 : :
561 [ + + ]: 99165 : dphi_ij(1:meshsz) = Pawtab(itypat)%phiphj(1:meshsz,kln)-Pawtab(itypat)%tphitphj(1:meshsz,kln)
562 [ + + ]: 99165 : rrdphi_ij(1:meshsz) = rr(1:meshsz)*dphi_ij(1:meshsz) ! will be used for first derivative
563 : :
564 : : ir=meshsz
565 [ + + ]: 703 : do while (ABS(dphi_ij(ir))<EPS)
566 : 582 : ir=ir-1
567 : : end do
568 : 121 : mmax=MIN(ir+1,meshsz)
569 : : ! ir is equal to meshsz if no point are below EPS
570 [ + + ]: 121 : if (mmax/=Pawrad(itypat)%int_meshsz) then ! mmax=meshsz
571 : : call pawrad_init(Tmpmesh,mesh_size=Pawtab(itypat)%mesh_size,mesh_type=Pawrad(itypat)%mesh_type, &
572 : 3 : & rstep=Pawrad(itypat)%rstep,lstep=Pawrad(itypat)%lstep,r_for_intg=rr(mmax))
573 : : else
574 : 118 : call pawrad_copy(Pawrad(itypat),Tmpmesh)
575 : : end if
576 : : !
577 : : ! === Loop on l for Bessel function. Note the starting point ===
578 : : ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm!
579 : : ! and only on lmax for this atom
580 [ + + ]: 574 : do ll=0,l_size-1
581 : : !
582 : : ! === Compute f_l(q=0) only if l=0, and first derivative fp_l(q=0) (nonzero only if ll==1) ===
583 : 453 : tmp_spl(1,1,ll,kln)=zero; yp1=zero
584 [ + + ]: 453 : if (ll==0) then
585 : 121 : call simp_gen(tmp_spl(1,1,ll,kln),dphi_ij,Tmpmesh)
586 : : end if
587 [ + + ]: 453 : if (ll==1) then
588 : 121 : call simp_gen(yp1,rrdphi_ij,Tmpmesh)
589 : 121 : yp1=yp1*two_pi*third
590 : : end if
591 : : !
592 : : ! === Compute f_l(0<q<qmax) ===
593 [ + - ]: 453 : if (nq_spl>2) then
594 [ + + ]: 1359000 : do iq=2,nq_spl-1
595 : 1358547 : arg=two_pi*qgrid_spl(iq)
596 [ + + ]: 1264258440 : do ir=1,mmax
597 : 1262899893 : qr=arg*rr(ir)
598 : 1262899893 : call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
599 : 1264258440 : ff(ir)=bes*dphi_ij(ir)
600 : : end do
601 : 1359000 : call simp_gen(tmp_spl(iq,1,ll,kln),ff,Tmpmesh)
602 : : end do
603 : : end if
604 : : !
605 : : ! === Compute f_l(q=qmax) and first derivative ===
606 [ + - ]: 453 : if (nq_spl>1) then
607 : 453 : argn=two_pi*qgrid_spl(nq_spl)
608 [ + + ]: 421560 : do ir=1,mmax
609 : 421107 : qr=argn*rr(ir)
610 : 421107 : call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
611 : 421107 : ff(ir)=bes * dphi_ij(ir)
612 : 421560 : gg(ir)=besp*rrdphi_ij(ir)
613 : : end do
614 : 453 : call simp_gen(tmp_spl(nq_spl,1,ll,kln),ff,Tmpmesh)
615 [ + + ]: 423405 : gg(:)=two_pi*gg(:) ! two_pi comes from 2\pi|q| in the Bessel function
616 : 453 : call simp_gen(ypn,gg,Tmpmesh)
617 : : else
618 : 0 : ypn=yp1
619 : : end if
620 : : !
621 : : ! === Compute second derivative of ff^{al}_{ij)(q) ===
622 : : !yp1=zero; ypn=zero
623 : 574 : call paw_spline(qgrid_spl,tmp_spl(:,1,ll,kln),nq_spl,yp1,ypn,tmp_spl(:,2,ll,kln))
624 : : end do !ll
625 : :
626 : 167 : call pawrad_free(Tmpmesh)
627 : :
628 : : end do !iln
629 : : end do !jln
630 : : !
631 : : ! === Save values for this atom type, each ll and kln channel ===
632 [ + + + + : 2720397 : pwff_spl = tmp_spl
+ + + + ]
633 : :
634 [ - + ]: 11 : ABI_FREE(dphi_ij)
635 [ - + ]: 11 : ABI_FREE(rrdphi_ij)
636 [ - + ]: 11 : ABI_FREE(ff)
637 [ - + ]: 11 : ABI_FREE(gg)
638 [ - + ]: 11 : ABI_FREE(rr)
639 [ - + ]: 11 : ABI_FREE(tmp_spl)
640 : :
641 : : !if (.FALSE.) then ! write form factors on file for plotting purpose.
642 : : ! ll=0
643 : : ! do iq=1,nq_spl
644 : : ! write(777+itypat,'(50(es16.8))')qgrid_spl(iq),((pwff_spl(iq,ider,ll,iln),ider=1,2),iln=1,dim2)
645 : : ! end do
646 : : !end if
647 : :
648 : : CASE (PWIJ_SHISHKIN)
649 : : ! ==== Shishkin-Kresse approximated expression ====
650 : : ! $ff_^{aL}_{ij}(q) = q_ij^{LM} \int_0^{r_a} j_L(2\pi qr) g_L(r) r^2 dr$
651 : : ! * Better description of multipoles of AE charge,
652 : : ! * Better results for energy degeneracies in GW band structure
653 : 15 : write(msg,'(a,i3)')' paw_mkrhox_spl: Using Shishkin-Kresse expression for atom type ',itypat
654 : 15 : call wrtout(std_out, msg)
655 : 15 : l_size = Pawtab(itypat)%l_size
656 : 15 : nlmn = Pawtab(itypat)%lmn_size
657 : :
658 : : !allocate(tmp_jgl(nq_spl,2,0:2*(Psps%mpsang-1)))
659 [ - + + - : 60 : ABI_MALLOC(tmp_jgl,(nq_spl,2,0:l_size-1))
- + + - +
- + - - +
- + ]
660 : :
661 : : ! Is mesh beginning with r=0 ?
662 [ - + ]: 15 : if (ABS(Pawrad(itypat)%rad(1))>tol10) then
663 : 0 : ABI_ERROR("Radial mesh starts with r/=0")
664 : : end if
665 : :
666 : : ! === Initialize temporary arrays and variables ===
667 : 15 : call pawrad_copy(Pawrad(itypat),Tmpmesh)
668 : 15 : meshsz=Pawtab(itypat)%mesh_size ; mmax=meshsz
669 [ - + + - : 45 : ABI_MALLOC(ff,(meshsz))
+ - + - -
+ - + ]
670 [ - + + - : 15285 : ABI_CALLOC(gg,(meshsz))
- + - + +
+ ]
671 [ - + + - : 30 : ABI_MALLOC(rr,(meshsz))
- + - + ]
672 [ - + + - : 30 : ABI_MALLOC(shape_l,(meshsz))
- + - + ]
673 [ - + + - : 30 : ABI_MALLOC(rrshape_l,(meshsz))
- + - + ]
674 [ + + ]: 15270 : rr(1:meshsz)=Tmpmesh%rad(1:meshsz)
675 : :
676 [ + + + + : 342300 : tmp_jgl(:,:,:)=zero
+ + ]
677 : 15 : rrshape_l(1)=zero
678 : 15 : shape_l(1)=zero
679 : : !
680 : : ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm!
681 : : ! and only lmax for this atom
682 : : !do ll=0,2*(Psps%mpsang-1)
683 [ + + ]: 72 : do ll=0,l_size-1
684 [ + + ]: 63315 : shape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**2
685 [ + + ]: 63315 : rrshape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**3
686 : : !
687 : : ! === Compute f_l(q=0) and first derivative fp_l(q=0) (only if ll==1) ===
688 : 57 : tmp_jgl(1,1,ll)=zero ; yp1=zero
689 [ + + ]: 57 : if (ll==0) then
690 : 15 : call simp_gen(tmp_jgl(1,1,ll),shape_l,Tmpmesh)
691 : : end if
692 [ + + ]: 57 : if (ll==1) then
693 : 15 : call simp_gen(yp1,rrshape_l,Tmpmesh) !rr comes from d/dq
694 : 15 : yp1=yp1*two_pi*third
695 : : end if
696 : : !
697 : : ! === Compute f_l(0<q<qmax) ===
698 [ + - ]: 57 : if (nq_spl>2) then
699 [ + + ]: 171000 : do iq=2,nq_spl-1
700 : 170943 : arg=two_pi*qgrid_spl(iq)
701 [ + + ]: 190052628 : do ir=1,mmax
702 : 189881685 : qr=arg*rr(ir)
703 : 189881685 : call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
704 : 190052628 : ff(ir)=bes*shape_l(ir)
705 : : end do
706 : 171000 : call simp_gen(tmp_jgl(iq,1,ll),ff,Tmpmesh)
707 : : end do
708 : : end if
709 : : !
710 : : ! === Compute f_l(q=qmax) and first derivative ===
711 [ + - ]: 57 : if (nq_spl>1) then
712 : 57 : argn=two_pi*qgrid_spl(nq_spl)
713 [ + + ]: 63372 : do ir=1,mmax
714 : 63315 : qr=argn*rr(ir)
715 : 63315 : call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
716 : 63315 : ff(ir)=bes * shape_l(ir)
717 : 63372 : gg(ir)=besp*rrshape_l(ir)
718 : : end do
719 : 57 : call simp_gen(tmp_jgl(nq_spl,1,ll),ff,Tmpmesh)
720 [ + + ]: 63372 : gg(:)=two_pi*gg(:) !two_pi comes from 2\pi|q|
721 : 57 : call simp_gen(ypn,gg,Tmpmesh)
722 : : else
723 : 0 : ypn=yp1
724 : : end if
725 : : !
726 : : ! === Compute second derivative of ff_^{al}_{ij)(q) ===
727 : 72 : call paw_spline(qgrid_spl,tmp_jgl(:,1,ll),nq_spl,yp1,ypn,tmp_jgl(:,2,ll))
728 : : end do !ll
729 : : !
730 : : ! === Save values for this type, each ll and ilmn,jlm channels ===
731 : : ! * Here we assembly q_{ij}^{lm} \int_0^{r_a} j_l(2\pi(q+G)r) g_l(r)r^2 dr
732 : : ! * Some of the contributions from qijl are zero due to Gaunt selection rules
733 [ + + ]: 177 : do jlmn=1,nlmn
734 [ + + ]: 1266 : do ilmn=1,jlmn
735 : 1089 : klmn=ilmn+(jlmn-1)*jlmn/2
736 : : !do ll=0,2*(Psps%mpsang-1)
737 [ + + ]: 6090 : do ll=0,l_size-1
738 [ + + ]: 28305 : do mm=-ll,ll
739 : 22377 : qlm=1+ll**2+ll+mm
740 [ + + + + ]: 134378724 : pwff_spl(:,:,qlm-1,klmn)=tmp_jgl(:,:,ll)*Pawtab(itypat)%qijl(qlm,klmn)
741 : : end do
742 : : end do
743 : : end do
744 : : end do
745 : :
746 : 15 : call pawrad_free(Tmpmesh)
747 : 15 : ABI_FREE(shape_l)
748 : 15 : ABI_FREE(rrshape_l)
749 [ - + ]: 15 : ABI_FREE(ff)
750 [ - + ]: 15 : ABI_FREE(gg)
751 [ - + ]: 15 : ABI_FREE(rr)
752 [ - + ]: 15 : ABI_FREE(tmp_jgl)
753 : :
754 : : CASE DEFAULT
755 [ + + - ]: 26 : ABI_BUG(sjoin('Called with wrong value for method:', itoa(method)))
756 : : END SELECT
757 : :
758 : : DBG_EXIT("COLL")
759 : :
760 [ - + - + : 26 : end subroutine paw_mkrhox_spl
- + - + -
+ - + - +
- + - + ]
761 : : !!***
762 : :
763 : : !----------------------------------------------------------------------
764 : :
765 : : !!****f* m_pawpwij/paw_mkrhox
766 : : !! NAME
767 : : !! paw_mkrhox
768 : : !!
769 : : !! FUNCTION
770 : : !! Evaluate $<phj|e^{-i(q+G)}|phi>-<tphj|e^{-i(q+G)}|tphi>$
771 : : !! for a fixed q-point and npw G vectors. Matrix elements are stored in packed storage mode.
772 : : !!
773 : : !! INPUTS
774 : : !! gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
775 : : !! gvec(3,npw)=G vectors in reduced coordinates
776 : : !! npw=numper of G vectors
777 : : !! Psps<pseudopotential_type>:
778 : : !! %lnmax=Max. number of (l,n) components over all type of PAW datasets
779 : : !! nq_spl=Number of points in the reciprocal space grid on which the radial functions pwff_spl are specified
780 : : !! qgrid_spl(nq_spl)=values at which form factors have been evaluated
781 : : !! %mpsang=1+maximum angular momentum
782 : : !! qpt(3)= q-point in reduced coordinates
783 : : !! ylm_q(npw,(2*Psps%mpsang-1)**2)=real spherical harmonics Ylm(q+G) for q-point qpt up to l=2*l_max
784 : : !! pwff_spl(nq_spl,2,0:2*(Psps%mpsang-1),Psps%lnmax*(Psps%lnmax+1)/2))
785 : : !! Pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
786 : : !! %indklmn(8,lmn2_size)=array giving klm, kln, abs(il-jl) and (il+jl), ilm and jlm, ilmn and jlmn for each klmn=(ilmn,jlmn)
787 : : !!
788 : : !! OUTPUT
789 : : !! paw_rhox(2,npw,lmn2_size): $<phj|e^{-i(q+G).r}|phi>-<tphj|e^{-i(q+G).r}|tphi>$ in packed form for
790 : : !! a type itypat (phase factor arising from atom position is not included)
791 : : !!
792 : : !! SOURCE
793 : :
794 : 1514 : subroutine paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,qgrid_spl,pwff_spl,&
795 : 1514 : gmet,qpt,npw,gvec,ylm_q,Psps,Pawtab,paw_rhox)
796 : :
797 : : !Arguments ------------------------------------
798 : : !scalars
799 : : integer,intent(in) :: itypat,dim1,dim2,method,npw,nq_spl,lmn2_size
800 : : type(Pseudopotential_type),intent(in) :: Psps
801 : : !arrays
802 : : integer,intent(in) :: gvec(3,npw)
803 : : real(dp),intent(in) :: gmet(3,3)
804 : : real(dp),intent(in) :: pwff_spl(nq_spl,2,0:dim1,dim2)
805 : : real(dp),intent(in) :: qpt(3),ylm_q(npw,(2*Psps%mpsang-1)**2)
806 : : real(dp),intent(out) :: paw_rhox(2,npw,lmn2_size)
807 : : real(dp),intent(in) :: qgrid_spl(nq_spl)
808 : : type(Pawtab_type),target,intent(in) :: Pawtab(Psps%ntypat)
809 : :
810 : : !Local variables-------------------------------
811 : : !scalars
812 : : integer :: ider,ig,ignt,il,ilm,ilm_G,ilmn,iln,im,ipow,jl,jlm
813 : : integer :: jlmn,jln,jm,k0lm,k0lmn,k0ln,klm,klmn,kln,ll_G,mm_G,mpsang,ngnt
814 : : real(dp) :: rgnt,dummy
815 : : character(len=500) :: msg
816 : : !arrays
817 : 1514 : integer,allocatable :: gntselect(:,:)
818 : 1514 : integer,contiguous, pointer :: indlmn(:,:)
819 : : real(dp) :: mi_l(2,0:3),qpg(3)
820 : 1514 : real(dp),allocatable :: derfun(:),newfun(:),qpg_norm(:),realgnt(:),wk_ffnl(:,:)
821 : :
822 : : ! *************************************************************************
823 : :
824 : : !write(std_out,*)itypat,dim1,dim2,method,npw,nq_spl,lmn2_size,Psps%mpsang
825 : 1514 : mpsang = Psps%mpsang
826 : 1514 : indlmn => Pawtab(itypat)%indlmn
827 : :
828 : : ! === Pre-calculate (-i)^l ===
829 : 1514 : mi_l(1,0)=one ; mi_l(2,0)=zero
830 : 1514 : mi_l(1,1)=zero ; mi_l(2,1)=-one
831 : 1514 : mi_l(1,2)=-one ; mi_l(2,2)=zero
832 : 1514 : mi_l(1,3)=zero ; mi_l(2,3)=one
833 : :
834 : : ! === Calculate |q+G| ===
835 : : ! * 2\pi is not included to be consistent with the spline.
836 [ + - + - : 4542 : ABI_MALLOC(qpg_norm,(npw))
+ - - + -
+ ]
837 [ + + ]: 101972 : do ig=1,npw
838 [ + + ]: 401832 : qpg = qpt + gvec(:,ig)
839 [ + + + + : 1608842 : qpg_norm(ig)=SQRT(DOT_PRODUCT(qpg,MATMUL(gmet,qpg)))
+ + ]
840 : : end do
841 : :
842 : : ! Check q-grid as %qgrid_spl must be large enoung.
843 [ + - - + : 4648514 : if (MAXVAL(qpg_norm)>MAXVAL(qgrid_spl)) then
- - + + +
- - + - -
+ + - + ]
844 : : write(msg,'(3a,f8.4,a,f8.4,2a)')&
845 : 0 : ' Function values are being requested outside range of data. ',ch10,&
846 [ # # # # : 0 : ' Max qpg_norm = ',MAXVAL(qpg_norm),' Max qgrid_spl = ',MAXVAL(qgrid_spl),ch10,&
# # # # #
# # # # #
# # ]
847 : 0 : ' Increase ecut(wfn), check qrid_ff and gsqcut '
848 : 0 : ABI_ERROR(msg)
849 : : end if
850 : :
851 [ - + + - : 4542 : ABI_MALLOC(wk_ffnl,(nq_spl,2))
+ - + - -
+ - + ]
852 [ - + + - : 3028 : ABI_MALLOC(newfun,(npw))
- + - + ]
853 [ - + + - : 3028 : ABI_MALLOC(derfun,(npw))
- + - + ]
854 : :
855 : 2961 : SELECT CASE (method)
856 : : CASE (PWIJ_ARNAUD)
857 : : ! === Arnaud-Alouani exact expression ===
858 : : ! * It does not describe the multipoles of the AE charge density
859 : : ! * $ 4\pi \sum_{LM} (-i)^l Y_M^L(q+G) G_{\li\mi\lj\mj}^{\LM} ff^{aL}_{ij}(|q+G|) $
860 : : ! where f has been calculated in paw_mkrhox_spl
861 : : !
862 : : ! === Re-evaluate Gaunt coefficients, just to be on the safe side ===
863 : : ! * Note that gntselect is in packed form, thanks to invariance under permutation.
864 : : ! * Could use Pawang% but size of gntselect depends on pawxcdev!
865 : :
866 [ - + + - : 4341 : ABI_MALLOC( realgnt,((2*mpsang-1)**2*(mpsang)**4))
+ - + - -
+ - + ]
867 [ - + + - : 5788 : ABI_MALLOC(gntselect,((2*mpsang-1)**2, mpsang**2*(mpsang**2+1)/2))
+ - - + +
- + - + -
- + - + ]
868 : 1447 : call realgaunt(mpsang,ngnt,gntselect,realgnt)
869 [ + + + + : 10960825 : paw_rhox=zero
+ + ]
870 : :
871 : : ! Loop on (jl,jm,jn) channels for this atom
872 [ + + ]: 12996 : do jlmn=1,Pawtab(itypat)%lmn_size
873 : 11549 : jl =indlmn(1,jlmn)
874 : 11549 : jm =indlmn(2,jlmn)
875 : 11549 : jlm=indlmn(4,jlmn)
876 : 11549 : jln=indlmn(5,jlmn)
877 : :
878 : 11549 : k0lmn=jlmn*(jlmn-1)/2
879 : 11549 : k0lm =jlm *(jlm -1)/2
880 : 11549 : k0ln =jln *(jln -1)/2
881 : :
882 : : ! === Loop on (il,im,in) channels; klmn is index for packed form ===
883 [ + + ]: 65139 : do ilmn=1,jlmn
884 : 52143 : il =indlmn(1,ilmn)
885 : 52143 : im =indlmn(2,ilmn)
886 : 52143 : ilm=indlmn(4,ilmn)
887 : 52143 : iln=indlmn(5,ilmn)
888 : :
889 : 52143 : klmn=k0lmn+ilmn
890 : 52143 : klm = Pawtab(itypat)%indklmn(1,klmn) !SF klm =k0lm +ilm is incorrect,
891 : : ! because ilm can be larger than jlm if more than 1 proj per angular momentum
892 : 52143 : kln =k0ln +iln ! This problem does not occur for kln, because always jln>=iln
893 : :
894 : : ! === Summing over allowed (L,M), taking into account Gaunt selection rules ===
895 [ - + + + ]: 146507 : do ll_G=ABS(jl-il),jl+il,2
896 : 82815 : ipow=MOD(ll_G,4)
897 : 82815 : ider=0
898 [ + + + + ]: 497304075 : wk_ffnl(:,:)=pwff_spl(:,:,ll_G,kln)
899 : 82815 : call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw)
900 [ + + ]: 377481 : do mm_G=-ll_G,ll_G
901 : 242523 : ilm_G=1+ll_G**2+ll_G+mm_G
902 : 242523 : ignt=gntselect(ilm_G,klm)
903 [ + + ]: 242523 : if (ignt==0) CYCLE
904 : 74314 : rgnt=realgnt(ignt)
905 : :
906 : : ! === Evaluate matrix elements for each plane wave ===
907 [ + + ]: 5452051 : do ig=1,npw
908 : 5294922 : dummy = newfun(ig) * ylm_q(ig,ilm_G) * rgnt
909 : 5294922 : paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) + ( dummy * mi_l(1,ipow) )
910 : 5537445 : paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) + ( dummy * mi_l(2,ipow) )
911 : : end do
912 : : end do !mm_G
913 : : end do !ll_G
914 : :
915 : : end do !ilmn
916 : : end do !jlmn
917 : :
918 : : ! Multiply by 4\pi arising from the expansion of the plane wave
919 [ + + + + : 10960825 : paw_rhox = four_pi*paw_rhox
+ + ]
920 [ - + ]: 1447 : ABI_FREE(realgnt)
921 [ - + ]: 2894 : ABI_FREE(gntselect)
922 : :
923 : : CASE (PWIJ_SHISHKIN)
924 : : ! === Shishkin-Kresse approximated expression ====
925 : : ! * Better description of multipoles of AE charge,
926 : : ! * Better results for energy degeneracies in GW band structure
927 : : ! * $4\pi \sum_{LM} q_ij^{LM} Y_M^L(q+G) f^{aL}_{ij}(q+G)$ where f has been calculated in paw_mkrhox_spl
928 : : !
929 [ + + + + : 1215163 : paw_rhox=zero
+ + ]
930 : : !
931 : : ! === Loop on (jl,jm,jn) channels for this atom ===
932 : : !itypat=Cryst%typat(iatm)
933 [ + + ]: 774 : do jlmn=1,Pawtab(itypat)%lmn_size
934 : 707 : jl =indlmn(1,jlmn)
935 : 707 : jm =indlmn(2,jlmn)
936 : 707 : jlm=indlmn(4,jlmn)
937 : 707 : jln=indlmn(5,jlmn)
938 : :
939 : 707 : k0lmn=jlmn*(jlmn-1)/2
940 : 707 : k0lm =jlm *(jlm -1)/2
941 : 707 : k0ln =jln *(jln -1)/2
942 : : !
943 : : ! === Loop on (il,im,in) channels; klmn is index for packed form ===
944 [ + + ]: 5478 : do ilmn=1,jlmn
945 : 4704 : il =indlmn(1,ilmn)
946 : 4704 : im =indlmn(2,ilmn)
947 : 4704 : ilm=indlmn(4,ilmn)
948 : 4704 : iln=indlmn(5,ilmn)
949 : :
950 : 4704 : klmn=k0lmn+ilmn
951 : : ! klm =k0lm +ilm !SF klm incorrect, see AA, but not used in SK
952 : 4704 : kln =k0ln +iln
953 : : !
954 : : ! === Summing over allowed (l,m), taking into account Gaunt selection rules ===
955 [ - + + + ]: 14387 : do ll_G=ABS(jl-il),jl+il,2
956 : 8976 : ipow=MOD(ll_G,4)
957 [ + + ]: 50736 : do mm_G=-ll_G,ll_G
958 : : ! here I can move splfit before the loop over mm_G but I have to change paw_rhox_spl
959 : 37056 : ilm_G=1+ll_G**2+ll_G+mm_G
960 : 37056 : ider=0
961 [ + + + + ]: 222521280 : wk_ffnl(:,:)=pwff_spl(:,:,ilm_G-1,klmn) ! Note klmn and ilm_G-1
962 : 37056 : call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw)
963 : : !
964 : : ! === Evaluate matrix elements for each plane wave ===
965 [ + + ]: 3504384 : do ig=1,npw
966 : 3458352 : dummy = newfun(ig) * ylm_q(ig,ilm_G)
967 : : paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) &
968 : 3458352 : + dummy * mi_l(1,ipow) !(ph3d(1,ig)*mi_l(1,ipow)-ph3d(2,ig)*mi_l(2,ipow))
969 : : paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) &
970 : 3495408 : + dummy * mi_l(2,ipow) !(ph3d(1,ig)*mi_l(2,ipow)+ph3d(2,ig)*mi_l(1,ipow))
971 : : end do
972 : : end do !mm_G
973 : : end do !ll_G
974 : :
975 : : end do !ilmn
976 : : end do !jlmn
977 : :
978 : : ! Multiply by 4\pi arising from the expansion of the plane wave
979 [ + + + + : 1215163 : paw_rhox=four_pi*paw_rhox
+ + ]
980 : :
981 : : CASE DEFAULT
982 [ + + - ]: 1514 : ABI_BUG(sjoin('Wrong value for method:', itoa(method)))
983 : : END SELECT
984 : :
985 [ - + ]: 1514 : ABI_FREE(wk_ffnl)
986 [ - + ]: 1514 : ABI_FREE(newfun)
987 [ - + ]: 1514 : ABI_FREE(derfun)
988 [ - + ]: 1514 : ABI_FREE(qpg_norm)
989 : :
990 [ - + - + : 3028 : end subroutine paw_mkrhox
- + ]
991 : : !!***
992 : :
993 : : !----------------------------------------------------------------------
994 : :
995 : : !!****f* m_pawpwij/paw_rho_tw_g
996 : : !! NAME
997 : : !! paw_rho_tw_g
998 : : !!
999 : : !! FUNCTION
1000 : : !! Evaluates the PAW onsite contribution to the oscillator strengths:
1001 : : !!
1002 : : !! sum_{i,j} <\tpsi_{k-q,b1}|\cprj_i> <\cprj_j|\tpsi_{k,b2}>*
1003 : : !! \[ <\phi_i|e^{-i(q+G).r}|\phi_j> - <\tilde\phi_i|e^{-i(q+G).r}|\tilde\phi_j> \].
1004 : : !!
1005 : : !! INPUTS
1006 : : !! dim_rtwg=Define the size of the array rhotwg
1007 : : !! === for nspinor==1 ===
1008 : : !! dim_rtwg=1
1009 : : !! === for nspinor==2 ===
1010 : : !! dim_rtwg=2 if only <up|up>, <dwn|dwn> matrix elements are required
1011 : : !! dim_rtwg=4 for <up|up>, <dwn|dwn>, <up|dwn> and <dwn|up>.
1012 : : !! nspinor=number of spinorial components.
1013 : : !! npw=number of plane waves for oscillator matrix elements
1014 : : !! Cprj_kmqb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
1015 : : !! projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
1016 : : !! wavefunctions (k-q,b1,s) and (k,b2,s), respectively.
1017 : : !!
1018 : : !! SIDE EFFECTS
1019 : : !! rhotwg(npw*dim_rtwg)=Updated oscillator strengths with the on-site PAW contributions added.
1020 : : !!
1021 : : !! SOURCE
1022 : :
1023 : 295726 : pure subroutine paw_rho_tw_g(cryst, pwij, npw, dim_rtwg, nspinor, gvec, Cprj_kmqb1, Cprj_kb2, rhotwg)
1024 : :
1025 : : !Arguments ------------------------------------
1026 : : !scalars
1027 : : type(crystal_t),intent(in) :: cryst
1028 : : integer,intent(in) :: npw,nspinor,dim_rtwg
1029 : : !arrays
1030 : : integer,intent(in) :: gvec(3,npw)
1031 : : type(pawcprj_type),intent(in) :: Cprj_kmqb1(cryst%natom,nspinor), Cprj_kb2(cryst%natom,nspinor)
1032 : : type(pawpwij_t),intent(in) :: Pwij(cryst%ntypat)
1033 : : complex(gwp),intent(inout) :: rhotwg(npw*dim_rtwg)
1034 : :
1035 : : !Local variables-------------------------------
1036 : : !scalars
1037 : : integer :: ig,iat,nlmn,ilmn,jlmn,k0lmn,klmn,iab,isp1,isp2,spad,itypat
1038 : : real(dp) :: fij,re_psp,im_psp,re_pw,im_pw,arg
1039 : : !arrays
1040 : : integer,parameter :: spinor_idxs(2,4) = RESHAPE([1,1,2,2,1,2,2,1], [2, 4])
1041 : : real(dp) :: tmp(2),qpg(3),x0(3),ph3d(2)
1042 : : ! *************************************************************************
1043 : :
1044 : : ! Loop over the four spinorial combinations
1045 [ + + ]: 591452 : do iab=1,dim_rtwg
1046 : 295726 : isp1 = spinor_idxs(1,iab)
1047 : 295726 : isp2 = spinor_idxs(2,iab)
1048 : 295726 : spad = npw*(iab-1)
1049 : :
1050 [ + + ]: 13296578 : do ig=1,npw
1051 : : tmp(:)=zero
1052 [ + + ]: 39594786 : do iat=1,cryst%natom
1053 : 26889660 : itypat = cryst%typat(iat)
1054 : 26889660 : nlmn = Pwij(itypat)%lmn_size
1055 [ + + ]: 107558640 : x0(:) = cryst%xred(:,iat)
1056 : :
1057 : : ! Structure factor e^{-i(q+G)*xred}
1058 [ + + ]: 107558640 : qpg(:)= Pwij(itypat)%qpt(:) + gvec(:,ig)
1059 [ + + ]: 107558640 : arg=-two_pi*DOT_PRODUCT(qpg(:),x0)
1060 : 26889660 : ph3d(1)=COS(arg)
1061 : 26889660 : ph3d(2)=SIN(arg)
1062 : :
1063 : : ! Loop over [(jl,jm,jn), (il,im,in)] channels in packed storage mode.
1064 [ + + ]: 258202869 : do jlmn=1,nlmn
1065 : 218608083 : k0lmn=jlmn*(jlmn-1)/2
1066 [ + + ]: 1255088349 : do ilmn=1,jlmn
1067 : : re_psp = Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) &
1068 : : +Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) &
1069 : : +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn) &
1070 : 1009590606 : +Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn)
1071 : :
1072 : : im_psp = Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) &
1073 : : -Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) &
1074 : : +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn) &
1075 : 1009590606 : -Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn)
1076 : :
1077 [ + + ]: 1009590606 : klmn=k0lmn+ilmn; fij=one; if (jlmn==ilmn) fij=half
1078 : :
1079 : : ! Multiply by the phase due to the atom position.
1080 : : re_pw = Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(1) &
1081 : 1009590606 : -Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(2)
1082 : :
1083 : : im_pw = Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(2) &
1084 : 1009590606 : +Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(1)
1085 : :
1086 : 1009590606 : tmp(1)=tmp(1)+ fij * (re_pw*re_psp - im_pw*im_psp)
1087 : 1228198689 : tmp(2)=tmp(2)+ fij * (re_pw*im_psp + im_pw*re_psp)
1088 : : end do !ilmn
1089 : : end do !jlmn
1090 : : end do !iat
1091 : :
1092 : : !if(ig==1) write(std_out,*) " TOTAL PW osc str = ",rhotwg(ig+spad)
1093 : : !if(ig==1) write(std_out,*) " TOTAL PAW osc str = ",tmp(1),tmp(2)
1094 : : ! Update input data using the appropriate index.
1095 : 13000852 : rhotwg(ig+spad) = rhotwg(ig+spad) + CMPLX(tmp(1),tmp(2),kind=gwp)
1096 : : !if(ig==1) write(std_out,*) " TOTAL PW+PAW osc str = ",rhotwg(ig+spad)
1097 : : end do !ig
1098 : :
1099 : : end do ! dim_rtwg
1100 : :
1101 : 295726 : end subroutine paw_rho_tw_g
1102 : : !!***
1103 : :
1104 : : !----------------------------------------------------------------------
1105 : :
1106 : : !!****f* ABINIT/paw_cross_rho_tw_g
1107 : : !! NAME
1108 : : !! paw_cross_rho_tw_g
1109 : : !!
1110 : : !! FUNCTION
1111 : : !! Compute the cross term between the PAW onsite part and plane-wave part in rho_tw
1112 : : !!
1113 : : !! INPUTS
1114 : : !!
1115 : : !! OUTPUT
1116 : : !!
1117 : : !! SOURCE
1118 : :
1119 : 0 : subroutine paw_cross_rho_tw_g(nspinor,npwvec,nr,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
1120 : 0 : ur_ae1,ur_ae_onsite1,ur_ps_onsite1,i1,ktabr1,ktabp1,spinrot1,&
1121 : 0 : ur_ae2,ur_ae_onsite2,ur_ps_onsite2,i2,ktabr2,ktabp2,spinrot2,&
1122 [ # # ]: 0 : dim_rtwg,rhotwg)
1123 : :
1124 : : !Arguments ------------------------------------
1125 : : !scalars
1126 : : integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft
1127 : : complex(dp),intent(in) :: ktabp1,ktabp2
1128 : : !arrays
1129 : : integer,intent(in) :: gbound(:,:)
1130 : : integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18)
1131 : : integer,intent(in) :: ktabr1(nr),ktabr2(nr)
1132 : : real(dp),intent(in) :: spinrot1(4),spinrot2(4)
1133 : : complex(gwp),intent(in) :: ur_ae1(nr),ur_ae2(nr)
1134 : : complex(gwp),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr)
1135 : : complex(gwp),intent(in) :: ur_ps_onsite1(nr),ur_ps_onsite2(nr)
1136 : : complex(gwp),intent(inout) :: rhotwg(npwvec*dim_rtwg)
1137 : :
1138 : : !Local variables-------------------------------
1139 : : !scalars
1140 : : integer,parameter :: ndat1 = 1, fftcache0 = 0, gpu_option_0 = 0
1141 : : integer :: ig,igfft,nx,ny,nz,ldx,ldy,ldz,mgfft,isprot1,isprot2
1142 : : type(fftbox_plan3_t) :: plan
1143 : : !arrays
1144 : 0 : complex(dp),allocatable :: usk(:),uu(:),rho(:)
1145 : : ! *************************************************************************
1146 : :
1147 : 0 : SELECT CASE (nspinor)
1148 : :
1149 : : CASE (1) ! Collinear case.
1150 : :
1151 [ # # # # : 0 : ABI_MALLOC(uu,(nr))
# # # # #
# # # ]
1152 [ # # # # : 0 : ABI_MALLOC(usk,(nr))
# # # # ]
1153 [ # # # # : 0 : ABI_MALLOC(rho,(nr))
# # # # ]
1154 : :
1155 [ # # # # : 0 : uu = ur_ae1(ktabr1)*ktabp1 - ur_ae_onsite1(ktabr1)*ktabp1; if (i1==1) uu = CONJG(uu)
# # # # ]
1156 [ # # # # : 0 : usk = ur_ae_onsite2(ktabr2)*ktabp2 - ur_ps_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk)
# # # # ]
1157 [ # # # # ]: 0 : rho = uu * usk
1158 : :
1159 [ # # # # : 0 : uu = ur_ae_onsite1(ktabr1)*ktabp1 - ur_ps_onsite1(ktabr1)*ktabp1; if (i1==1) uu = CONJG(uu)
# # # # ]
1160 [ # # # # : 0 : usk = ur_ae2(ktabr2)*ktabp2 - ur_ae_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk)
# # # # ]
1161 [ # # # # ]: 0 : rho = rho + uu * usk
1162 : :
1163 : 0 : SELECT CASE (map2sphere)
1164 : :
1165 : : CASE (0)
1166 : : ! Need results on the full FFT box thus cannot use zero-padded FFT.
1167 : 0 : call plan%init(ndat1, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
1168 : 0 : call plan%execute(rho, -1, ndat1)
1169 : 0 : call plan%free()
1170 : :
1171 [ # # ]: 0 : rhotwg=rhotwg + rho
1172 : :
1173 : : CASE (1)
1174 : : ! Need results on the G-sphere. Call zero-padded FFT routines if required.
1175 [ # # ]: 0 : if (use_padfft==1) then
1176 [ # # ]: 0 : nx =ngfft(1); ny =ngfft(2); nz =ngfft(3); mgfft = MAXVAL(ngfft(1:3))
1177 : 0 : ldx=nx ; ldy=ny ; ldz=nz
1178 [ # # # # : 0 : call fftpad(rho,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat1,mgfft,-1,gbound)
# # # # #
# ]
1179 : : else
1180 : 0 : call plan%init(ndat1, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0)
1181 : 0 : call plan%execute(rho, -1, ndat1)
1182 : 0 : call plan%free()
1183 : : end if
1184 : :
1185 : : ! Have to map FFT to G-sphere.
1186 [ # # ]: 0 : do ig=1,npwvec
1187 : 0 : igfft=igfftg0(ig)
1188 : : ! If G-G0 belong to the FFT mesh.
1189 [ # # ]: 0 : if (igfft/=0) rhotwg(ig)=rhotwg(ig)+rho(igfft)
1190 : : end do
1191 : :
1192 : : CASE DEFAULT
1193 [ # # # ]: 0 : ABI_BUG(sjoin("Wrong map2sphere:", itoa(map2sphere)))
1194 : : END SELECT
1195 : :
1196 : 0 : RETURN
1197 : :
1198 : : CASE (2)
1199 : : ! Spinorial case.
1200 : 0 : isprot1=spinrot1(1); isprot2=spinrot2(1) ! This is to bypass abirule
1201 : 0 : ABI_ERROR("Spinorial case not implemented yet")
1202 : :
1203 : 0 : SELECT CASE (map2sphere)
1204 : :
1205 : : CASE (0) ! Need results on the full FFT box thus cannot use zero-padded FFT.
1206 : : CASE (1) ! Need results on the G-sphere. Call zero-padded FFT routines if required.
1207 : : CASE DEFAULT
1208 [ # # ]: 0 : ABI_BUG("Wrong map2sphere")
1209 : : END SELECT
1210 : 0 : RETURN
1211 : :
1212 : : CASE DEFAULT
1213 : 0 : ABI_BUG(sjoin('Wrong nspinor:', itoa(nspinor)))
1214 : : END SELECT
1215 : :
1216 [ # # # # : 0 : end subroutine paw_cross_rho_tw_g
# # # # #
# # # # ]
1217 : : !!***
1218 : :
1219 : : !----------------------------------------------------------------------
1220 : :
1221 [ # # # # : 0 : END MODULE m_pawpwij
# # # # #
# # # #
# ]
1222 : : !!***
|
