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 196 : 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 32 : do itypat=1,Psps%ntypat
199 20 : 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 9 : dim1 = Pawtab(itypat)%l_size**2
207 9 : dim2 = Pawtab(itypat)%lmn2_size
208 : case default
209 20 : ABI_BUG(sjoin("Wrong method:", itoa(method)))
210 : end select
211 :
212 20 : Paw_pwff(itypat)%dim1 = dim1
213 20 : Paw_pwff(itypat)%dim2 = dim2
214 260 : Paw_pwff(itypat)%gmet = gmet
215 :
216 : ! Setup of the q-mesh for spline. It can be type-dependent.
217 20 : nq = nq_spl(itypat)
218 20 : dq = qmax(itypat)/(one*(nq-1))
219 : !write(std_out,*)"nq,dq",nq,dq
220 :
221 20 : Paw_pwff(itypat)%nq_spl = nq
222 60 : ABI_MALLOC(Paw_pwff(itypat)%qgrid_spl,(nq))
223 60060 : Paw_pwff(itypat)%qgrid_spl = arth(zero,dq,nq)
224 : !
225 : ! === Calculate form factors depending on method ===
226 100 : 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 32 : Paw_pwff(itypat)%qgrid_spl,Pawrad,Pawtab,Paw_pwff(itypat)%pwff_spl)
230 : end do ! itypat
231 :
232 12 : 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 12 : 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 32 : do ii=1,SIZE(Paw_pwff)
258 20 : ABI_SFREE(Paw_pwff(ii)%qgrid_spl)
259 32 : ABI_SFREE(Paw_pwff(ii)%pwff_spl)
260 : end do
261 :
262 12 : 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 4402 : 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 1434 : type(MPI_type) :: MPI_enreg_seq
315 : !arrays
316 1434 : integer,allocatable :: npwarr(:),dummy_nband(:)
317 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
318 1434 : real(dp),allocatable :: my_qtmp(:,:), ylm_q(:,:),ylmgr_q(:,:,:)
319 :
320 : ! *********************************************************************
321 :
322 : ! ===============================================
323 : ! === Get real spherical harmonics in G space ===
324 : ! ===============================================
325 1434 : 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 1434 : my_mqmem=1; two_lmaxp1=2*Psps%mpsang-1; optder=0
329 :
330 5736 : ABI_MALLOC(ylm_q ,(npw*my_mqmem,two_lmaxp1**2))
331 7170 : ABI_MALLOC(ylmgr_q,(npw*my_mqmem,3+6*(optder/2),two_lmaxp1**2))
332 :
333 1434 : my_nqpt=1
334 1434 : ABI_MALLOC(my_qtmp,(3,my_nqpt))
335 5736 : my_qtmp(:,1)=qpt_in(:)
336 :
337 : ! dummy_nband and dummy_nsppol are not used in sequential mode.
338 1434 : dummy_nsppol=1
339 1434 : ABI_MALLOC(dummy_nband,(my_nqpt*dummy_nsppol))
340 2868 : dummy_nband=0
341 1434 : ABI_MALLOC(npwarr,(my_nqpt))
342 2868 : npwarr(:)=npw
343 :
344 : ! Fake MPI_type for sequential part.
345 1434 : 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 1434 : npwarr,dummy_nsppol,optder,rprimd,ylm_q,ylmgr_q)
349 :
350 1434 : call destroy_mpi_enreg(MPI_enreg_seq)
351 :
352 1434 : ABI_FREE(my_qtmp)
353 1434 : ABI_FREE(dummy_nband)
354 1434 : ABI_FREE(npwarr)
355 1434 : ABI_FREE(ylmgr_q)
356 :
357 : ! Construct the Pwij structure
358 2918 : do itypat=1,Psps%ntypat
359 :
360 1484 : Pwij(itypat)%istpw = 1
361 1484 : Pwij(itypat)%npw = npw
362 :
363 1484 : lmn_size = Pawtab(itypat)%lmn_size
364 1484 : lmn2_size = lmn_size*(lmn_size+1)/2
365 1484 : Pwij(itypat)%lmn_size = lmn_size
366 1484 : Pwij(itypat)%lmn2_size = lmn2_size
367 :
368 5936 : Pwij(itypat)%qpt(:) = qpt_in(:)
369 :
370 : ! Prepare the call to paw_mkrhox
371 1484 : method = Paw_pwff(itypat)%method
372 1484 : dim1 = Paw_pwff(itypat)%dim1
373 1484 : dim2 = Paw_pwff(itypat)%dim2
374 1484 : nq_spl = Paw_pwff(itypat)%nq_spl
375 19292 : gmet = Paw_pwff(itypat)%gmet
376 :
377 5936 : 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 2918 : gmet,qpt_in,npw,gvec,ylm_q,Psps,Pawtab,Pwij(itypat)%mqpgij)
382 : end do ! itypat
383 :
384 1434 : ABI_FREE(ylm_q)
385 :
386 1434 : 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 1434 : 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 2918 : do ii=1,SIZE(Pwij)
410 2918 : ABI_SFREE(Pwij(ii)%mqpgij)
411 : end do
412 :
413 1434 : 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 20 : 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 20 : type(Pawrad_type) :: Tmpmesh
512 : !arrays
513 20 : real(dp),allocatable :: rrshape_l(:),shape_l(:),ff(:),gg(:),rr(:),rrdphi_ij(:)
514 20 : real(dp),allocatable :: dphi_ij(:),tmp_spl(:,:,:,:),tmp_jgl(:,:,:)
515 :
516 : !*************************************************************************
517 :
518 : DBG_ENTER("COLL")
519 :
520 57054114 : pwff_spl=zero
521 :
522 31 : 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 9 : write(msg,'(a,i3)')' paw_mkrhox_spl: Using Shishkin-Kresse expression for atom type ',itypat
654 9 : call wrtout(std_out, msg)
655 9 : l_size = Pawtab(itypat)%l_size
656 9 : nlmn = Pawtab(itypat)%lmn_size
657 :
658 : !allocate(tmp_jgl(nq_spl,2,0:2*(Psps%mpsang-1)))
659 36 : ABI_MALLOC(tmp_jgl,(nq_spl,2,0:l_size-1))
660 :
661 : ! Is mesh beginning with r=0 ?
662 9 : 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 9 : call pawrad_copy(Pawrad(itypat),Tmpmesh)
668 9 : meshsz=Pawtab(itypat)%mesh_size ; mmax=meshsz
669 27 : ABI_MALLOC(ff,(meshsz))
670 6759 : ABI_CALLOC(gg,(meshsz))
671 18 : ABI_MALLOC(rr,(meshsz))
672 18 : ABI_MALLOC(shape_l,(meshsz))
673 18 : ABI_MALLOC(rrshape_l,(meshsz))
674 6750 : rr(1:meshsz)=Tmpmesh%rad(1:meshsz)
675 :
676 198174 : tmp_jgl(:,:,:)=zero
677 9 : rrshape_l(1)=zero
678 9 : 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 42 : do ll=0,l_size-1
684 28611 : shape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**2
685 28611 : 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 33 : tmp_jgl(1,1,ll)=zero ; yp1=zero
689 33 : if (ll==0) then
690 9 : call simp_gen(tmp_jgl(1,1,ll),shape_l,Tmpmesh)
691 : end if
692 33 : if (ll==1) then
693 9 : call simp_gen(yp1,rrshape_l,Tmpmesh) !rr comes from d/dq
694 9 : yp1=yp1*two_pi*third
695 : end if
696 : !
697 : ! === Compute f_l(0<q<qmax) ===
698 33 : if (nq_spl>2) then
699 99000 : do iq=2,nq_spl-1
700 98967 : arg=two_pi*qgrid_spl(iq)
701 85903356 : do ir=1,mmax
702 85804389 : qr=arg*rr(ir)
703 85804389 : call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
704 85903356 : ff(ir)=bes*shape_l(ir)
705 : end do
706 99000 : 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 33 : if (nq_spl>1) then
712 33 : argn=two_pi*qgrid_spl(nq_spl)
713 28644 : do ir=1,mmax
714 28611 : qr=argn*rr(ir)
715 28611 : call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
716 28611 : ff(ir)=bes * shape_l(ir)
717 28644 : gg(ir)=besp*rrshape_l(ir)
718 : end do
719 33 : call simp_gen(tmp_jgl(nq_spl,1,ll),ff,Tmpmesh)
720 28644 : gg(:)=two_pi*gg(:) !two_pi comes from 2\pi|q|
721 33 : 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 42 : 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 93 : do jlmn=1,nlmn
734 561 : do ilmn=1,jlmn
735 468 : klmn=ilmn+(jlmn-1)*jlmn/2
736 : !do ll=0,2*(Psps%mpsang-1)
737 2502 : do ll=0,l_size-1
738 10998 : do mm=-ll,ll
739 8580 : qlm=1+ll**2+ll+mm
740 51524850 : 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 9 : call pawrad_free(Tmpmesh)
747 9 : ABI_FREE(shape_l)
748 9 : ABI_FREE(rrshape_l)
749 9 : ABI_FREE(ff)
750 9 : ABI_FREE(gg)
751 9 : ABI_FREE(rr)
752 9 : ABI_FREE(tmp_jgl)
753 :
754 : CASE DEFAULT
755 20 : ABI_BUG(sjoin('Called with wrong value for method:', itoa(method)))
756 : END SELECT
757 :
758 : DBG_EXIT("COLL")
759 :
760 20 : 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 1484 : subroutine paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,qgrid_spl,pwff_spl,&
795 1484 : 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 1484 : integer,allocatable :: gntselect(:,:)
818 1484 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
819 : real(dp) :: mi_l(2,0:3),qpg(3)
820 1484 : 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 1484 : mpsang = Psps%mpsang
826 1484 : indlmn => Pawtab(itypat)%indlmn
827 :
828 : ! === Pre-calculate (-i)^l ===
829 1484 : mi_l(1,0)=one ; mi_l(2,0)=zero
830 1484 : mi_l(1,1)=zero ; mi_l(2,1)=-one
831 1484 : mi_l(1,2)=-one ; mi_l(2,2)=zero
832 1484 : 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 4452 : ABI_MALLOC(qpg_norm,(npw))
837 98564 : do ig=1,npw
838 388320 : qpg = qpt + gvec(:,ig)
839 1554764 : 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 4555016 : 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 4452 : ABI_MALLOC(wk_ffnl,(nq_spl,2))
852 2968 : ABI_MALLOC(newfun,(npw))
853 2968 : ABI_MALLOC(derfun,(npw))
854 :
855 2931 : 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 163159 : paw_rhox=zero
930 : !
931 : ! === Loop on (jl,jm,jn) channels for this atom ===
932 : !itypat=Cryst%typat(iatm)
933 354 : do jlmn=1,Pawtab(itypat)%lmn_size
934 317 : jl =indlmn(1,jlmn)
935 317 : jm =indlmn(2,jlmn)
936 317 : jlm=indlmn(4,jlmn)
937 317 : jln=indlmn(5,jlmn)
938 :
939 317 : k0lmn=jlmn*(jlmn-1)/2
940 317 : k0lm =jlm *(jlm -1)/2
941 317 : k0ln =jln *(jln -1)/2
942 : !
943 : ! === Loop on (il,im,in) channels; klmn is index for packed form ===
944 1953 : do ilmn=1,jlmn
945 1599 : il =indlmn(1,ilmn)
946 1599 : im =indlmn(2,ilmn)
947 1599 : ilm=indlmn(4,ilmn)
948 1599 : iln=indlmn(5,ilmn)
949 :
950 1599 : klmn=k0lmn+ilmn
951 : ! klm =k0lm +ilm !SF klm incorrect, see AA, but not used in SK
952 1599 : kln =k0ln +iln
953 : !
954 : ! === Summing over allowed (l,m), taking into account Gaunt selection rules ===
955 4607 : do ll_G=ABS(jl-il),jl+il,2
956 2691 : ipow=MOD(ll_G,4)
957 13521 : 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 9231 : ilm_G=1+ll_G**2+ll_G+mm_G
960 9231 : ider=0
961 55432155 : wk_ffnl(:,:)=pwff_spl(:,:,ilm_G-1,klmn) ! Note klmn and ilm_G-1
962 9231 : call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw)
963 : !
964 : ! === Evaluate matrix elements for each plane wave ===
965 337179 : do ig=1,npw
966 325257 : dummy = newfun(ig) * ylm_q(ig,ilm_G)
967 : paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) &
968 325257 : + 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 334488 : + 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 163159 : paw_rhox=four_pi*paw_rhox
980 :
981 : CASE DEFAULT
982 1484 : ABI_BUG(sjoin('Wrong value for method:', itoa(method)))
983 : END SELECT
984 :
985 1484 : ABI_FREE(wk_ffnl)
986 1484 : ABI_FREE(newfun)
987 1484 : ABI_FREE(derfun)
988 1484 : ABI_FREE(qpg_norm)
989 :
990 2968 : 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 294679 : 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 589358 : do iab=1,dim_rtwg
1046 294679 : isp1 = spinor_idxs(1,iab)
1047 294679 : isp2 = spinor_idxs(2,iab)
1048 294679 : spad = npw*(iab-1)
1049 :
1050 13226075 : do ig=1,npw
1051 : tmp(:)=zero
1052 39389559 : do iat=1,cryst%natom
1053 26752842 : itypat = cryst%typat(iat)
1054 26752842 : nlmn = Pwij(itypat)%lmn_size
1055 107011368 : x0(:) = cryst%xred(:,iat)
1056 :
1057 : ! Structure factor e^{-i(q+G)*xred}
1058 107011368 : qpg(:)= Pwij(itypat)%qpt(:) + gvec(:,ig)
1059 107011368 : arg=-two_pi*DOT_PRODUCT(qpg(:),x0)
1060 26752842 : ph3d(1)=COS(arg)
1061 26752842 : ph3d(2)=SIN(arg)
1062 :
1063 : ! Loop over [(jl,jm,jn), (il,im,in)] channels in packed storage mode.
1064 256219008 : do jlmn=1,nlmn
1065 216829449 : k0lmn=jlmn*(jlmn-1)/2
1066 1239012234 : 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 995429943 : +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 995429943 : -Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn)
1076 :
1077 995429943 : 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 995429943 : -Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(2)
1082 :
1083 : im_pw = Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(2) &
1084 995429943 : +Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(1)
1085 :
1086 995429943 : tmp(1)=tmp(1)+ fij * (re_pw*re_psp - im_pw*im_psp)
1087 1212259392 : 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 12931396 : 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 294679 : 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 : !!***
|
