Line data Source code
1 : !!****m* ABINIT/m_paw_pwaves_lmn
2 : !! NAME
3 : !! m_paw_pwaves_lmn
4 : !!
5 : !! FUNCTION
6 : !! This module contains the definition of the paw_pwaves_lmn structured datatype,
7 : !! as well as related functions and methods.
8 : !! paw_pwaves_lmn is used to store the 3D values of the all-electron and of
9 : !! the pseudized part of the PAW partial waves on the set of FFT points falling
10 : !! inside the spheres around each atom.
11 : !!
12 : !! COPYRIGHT
13 : !! Copyright (C) 2008-2026 ABINIT group (MG,MT)
14 : !! This file is distributed under the terms of the
15 : !! GNU General Public License, see ~abinit/COPYING
16 : !! or http://www.gnu.org/copyleft/gpl.txt .
17 : !!
18 : !! NOTES
19 : !! * Routines tagged with "@type_name" are strongly connected to the definition of the data type.
20 : !! Strongly connected means that the proper functioning of the implementation relies on the
21 : !! assumption that the tagged procedure is consistent with the type declaration.
22 : !! Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
23 : !! that all the strongly connected routines are changed accordingly to accommodate the modification of the data type.
24 : !! Typical examples of strongly connected routines are creation, destruction or reset methods.
25 : !!
26 : !! SOURCE
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : MODULE m_paw_pwaves_lmn
35 :
36 : use defs_basis
37 : use m_errors
38 : use m_abicore
39 : use m_xmpi
40 : use m_sort
41 :
42 : use m_numeric_tools, only : wrap2_zero_one
43 : use m_geometry, only : xcart2xred
44 : use m_pawrad, only : pawrad_type, nderiv_gen, pawrad_deducer0
45 : use m_pawtab, only : pawtab_type
46 : use m_pawfgrtab, only : pawfgrtab_type
47 : use m_paw_numeric, only : paw_spline, paw_splint
48 : use m_paw_sphharm, only : initylmr
49 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
50 :
51 : implicit none
52 :
53 : private
54 : !!***
55 :
56 : !----------------------------------------------------------------------
57 :
58 : !!****t* m_paw_pwaves_lmn/paw_pwaves_lmn_t
59 : !! NAME
60 : !! paw_pwaves_lmn_t
61 : !!
62 : !! FUNCTION
63 : !! Datatype used to store the 3D values of the all-electron and of the pseudized part of the
64 : !! PAW partial waves on the set of FFT points falling inside the spheres around each atom.
65 : !! The data is mainly used for plotting the true PAW wavefunctions in real space.
66 : !!
67 : !! SOURCE
68 :
69 : type,public :: paw_pwaves_lmn_t
70 :
71 : integer :: nfgd
72 :
73 : integer :: lmn_size
74 :
75 : !$integer :: ngfft(18)
76 :
77 : integer,allocatable :: r0shift(:,:)
78 : ! r0shift(3,nfgd)
79 :
80 : !real(dp),allocatable:: phk_atm(:,:)
81 : ! phk_atmt(2,nfgd)
82 :
83 : real(dp),allocatable :: phi(:,:)
84 : ! phi (nfgd,lmn_size)
85 : ! \phi_{nlm}(ifgd) for each point of the FFT mesh located inside the PAW sphere (see pawfgrtab_type).
86 :
87 : real(dp),allocatable :: tphi(:,:)
88 : ! tphi (nfgd,lmn_size)
89 : ! \tphi_{nlm}(ifgd) for each point of the FFT mesh located inside the PAW sphere (see pawfgrtab_type).
90 :
91 : real(dp),allocatable :: phi_gr(:,:,:)
92 : ! phi_gr (3,nfgd,lmn_size)
93 : ! gradient, in cartesian coordinates, of \phi_{nlm}(ifgd) for each point of the FFT mesh
94 : ! located inside the PAW sphere (see pawfgrtab_type).
95 :
96 : real(dp),allocatable :: tphi_gr(:,:,:)
97 : ! tphi_gr (3,nfgd,lmn_size)
98 : ! gradient, in cartesian coordinates, of \tphi_{nlm}(ifgd) for each point of the FFT mesh
99 : ! located inside the PAW sphere (see pawfgrtab_type).
100 :
101 : end type paw_pwaves_lmn_t
102 :
103 : public :: paw_pwaves_lmn_init
104 : public :: paw_pwaves_lmn_free
105 : !!***
106 :
107 : CONTAINS !===========================================================
108 : !!***
109 :
110 : !----------------------------------------------------------------------
111 :
112 : !!****f* m_paw_pwaves_lmn/paw_pwaves_lmn_init
113 : !! NAME
114 : !! paw_pwaves_lmn_init
115 : !!
116 : !! FUNCTION
117 : !!
118 : !! INPUTS
119 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
120 : !! comm_atom=--optional-- MPI communicator over atoms
121 : !!
122 : !! OUTPUT
123 : !!
124 : !! SOURCE
125 :
126 3 : subroutine paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,ntypat,rprimd,xcart,Pawtab, &
127 2 : & Pawrad,local_pawfgrtab,optgrad,&
128 1 : & mpi_atmtab,comm_atom) ! optional arguments (parallelism)
129 :
130 : !Arguments ------------------------------------
131 : !scalars
132 : integer,intent(in) :: my_natom,natom,ntypat
133 : integer,optional,intent(in) :: optgrad,comm_atom
134 : !arrays
135 : integer,optional,target,intent(in) :: mpi_atmtab(:)
136 : real(dp),intent(in) :: xcart(3,natom),rprimd(3,3)
137 : type(pawtab_type),target,intent(in) :: Pawtab(ntypat)
138 : type(pawrad_type),intent(in) :: Pawrad(ntypat)
139 : type(pawfgrtab_type),intent(in) :: local_pawfgrtab(my_natom)
140 : type(paw_pwaves_lmn_t),intent(out) :: Paw_onsite(my_natom)
141 :
142 : !Local variables-------------------------------
143 : !scalars
144 : integer :: itypat,ln_size,lmn_size,mesh_size,inl,iatom,iatom1,my_comm_atom,my_optgrad
145 : integer :: nfgd,ifgd,ipsang,option_ylmr,normchoice,ii,jlmn,jl,jm,jlm,jln
146 : logical :: my_atmtab_allocated,paral_atom
147 : real(dp) :: phj,rR,tphj,ybcbeg,ybcend
148 : !arrays
149 1 : integer, allocatable :: iperm(:)
150 1 : integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
151 1 : integer,pointer :: my_atmtab(:)
152 : real(dp) :: yvals(4),red(3),shift(3)
153 1 : real(dp),allocatable :: ff(:),nrm(:),nrm_sort(:),phigrd(:,:),tphigrd(:,:),ylm_tmp(:,:),ylm(:,:),ylm_gr(:,:,:)
154 1 : real(dp),allocatable :: rsph_red(:,:),rsph_cart(:,:),phigrd_gr(:,:),tphigrd_gr(:,:),gg(:)
155 1 : type(paw_pwaves_lmn_t),allocatable :: Paw_lmn_spline(:)
156 :
157 : ! *************************************************************************
158 :
159 : !@paw_pwaves_lmn_t
160 :
161 1 : if (my_natom==0) return
162 1 : my_optgrad = -1; if (present(optgrad)) my_optgrad = optgrad
163 :
164 2 : ABI_CHECK(all(local_pawfgrtab(:)%rfgd_allocated==1),"R vectors not allocated in pawfgrtab")
165 :
166 : ! Set up parallelism over atoms
167 1 : paral_atom=(present(comm_atom).and.(my_natom/=natom))
168 1 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
169 1 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
170 1 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
171 :
172 : ! Prepare the spline. Calculate 2nd derivatives of partial waves for each atom type.
173 4 : ABI_MALLOC(Paw_lmn_spline,(ntypat))
174 :
175 2 : do itypat=1,ntypat
176 1 : ln_size =Pawtab(itypat)%basis_size
177 1 : mesh_size=Pawtab(itypat)%mesh_size
178 :
179 4 : ABI_MALLOC(Paw_lmn_spline(itypat)%phi ,(mesh_size,ln_size))
180 3 : ABI_MALLOC(Paw_lmn_spline(itypat)%tphi,(mesh_size,ln_size))
181 :
182 7 : do inl=1,ln_size ! Calculate 2nd derivatives of %phi and %tphi for each ln component.
183 5 : ybcbeg=zero; ybcend=zero
184 : call paw_spline(Pawrad(itypat)%rad,Pawtab(itypat)%phi(:,inl), mesh_size,&
185 5 : & ybcbeg,ybcend,Paw_lmn_spline(itypat)%phi(:,inl))
186 :
187 : ybcbeg=zero; ybcend=zero
188 : call paw_spline(Pawrad(itypat)%rad,Pawtab(itypat)%tphi(:,inl),mesh_size,&
189 6 : & ybcbeg,ybcend,Paw_lmn_spline(itypat)%tphi(:,inl))
190 : end do
191 : end do
192 : !
193 : ! === spline for each atom ===
194 : ! * FFT points within PAW sphere depend on the atom site.
195 2 : do iatom1=1,my_natom
196 1 : iatom=iatom1;if (paral_atom) iatom=my_atmtab(iatom1)
197 :
198 1 : itypat = local_pawfgrtab(iatom1)%itypat
199 1 : ln_size = Pawtab(itypat)%basis_size
200 1 : lmn_size = Pawtab(itypat)%lmn_size
201 1 : mesh_size= Pawrad(itypat)%mesh_size
202 1 : nfgd = local_pawfgrtab(iatom1)%nfgd ! no. of points in the fine grid for this PAW sphere
203 1 : indlmn => Pawtab(itypat)%indlmn
204 :
205 : ! The points in the PAW sphere might belong to a different unit cell, in this case one has to
206 : ! reconstruct the contribution to the AE psi_k given by the lattice-symmetric atom in another sphere.
207 : ! Here I wrap rr back into the first unit cell keeping trace of the lattice vector L
208 : ! needed so that rr = rr_first_ucell + L
209 : ! The contribution to the AE u(r) in the first unit cell has to be multiplied by e^{-ikL}.
210 :
211 1 : Paw_onsite(iatom1)%nfgd = nfgd
212 1 : Paw_onsite(iatom1)%lmn_size = lmn_size
213 : !
214 3 : ABI_MALLOC(Paw_onsite(iatom1)%r0shift,(3,nfgd))
215 :
216 3 : ABI_MALLOC(rsph_red,(3,nfgd))
217 2 : ABI_MALLOC(rsph_cart,(3,nfgd))
218 1430 : do ifgd=1,nfgd
219 5717 : rsph_cart(:,ifgd) = local_pawfgrtab(iatom1)%rfgd(:,ifgd) + xcart(:,iatom)
220 : end do
221 1 : call xcart2xred(nfgd,rprimd,rsph_cart,rsph_red) ! go to reduced coordinates.
222 1 : ABI_FREE(rsph_cart)
223 :
224 1430 : do ifgd=1,nfgd
225 1429 : call wrap2_zero_one(rsph_red(1,ifgd),red(1),shift(1)) ! rr = r_cell + shift
226 1429 : call wrap2_zero_one(rsph_red(2,ifgd),red(2),shift(2))
227 1429 : call wrap2_zero_one(rsph_red(3,ifgd),red(3),shift(3))
228 5717 : Paw_onsite(iatom1)%r0shift(:,ifgd) = NINT(shift)
229 : !if (ANY( ABS(shift) > tol12)) then
230 : !ABI_WARNING("rmR_red is outside the first unit cell.")
231 : !write(ab_out,*)rsph_red(:,ifgd),shift
232 : !end if
233 : end do
234 1 : ABI_FREE(rsph_red)
235 : !
236 : ! * Obtain |r-R| on fine grid, note that rfgd is given in Cartesian coordinates.
237 3 : ABI_MALLOC(nrm,(nfgd))
238 1430 : do ifgd=1,nfgd
239 5717 : nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom1)%rfgd(:,ifgd),local_pawfgrtab(iatom1)%rfgd(:,ifgd)))
240 : end do
241 : !
242 : ! * Compute Ylm for each r-R vector.
243 1 : ipsang = 1 + (Pawtab(itypat)%l_size-1)/2 ! recall l_size=2*l_max-1 where l_max is shifted by 1.
244 4 : ABI_MALLOC(ylm_tmp,(ipsang**2,nfgd))
245 1 : normchoice = 1 ! Use computed norms of input vectors.
246 1 : if (my_optgrad==1) then
247 0 : option_ylmr=2 ! Compute Ylm(r-R) and its gradient
248 0 : ABI_MALLOC(ylm_gr,(3,ipsang**2,nfgd))
249 : else
250 1 : option_ylmr= 1 ! To compute Ylm(r-R).
251 1 : ABI_MALLOC(ylm_gr,(3,3,0))
252 : end if
253 1 : call initylmr(ipsang,normchoice,nfgd,nrm,option_ylmr,local_pawfgrtab(iatom1)%rfgd,ylm_tmp,ylm_gr)
254 : !
255 : ! Exchange dimensions for better memory access.
256 3 : ABI_MALLOC(ylm,(nfgd,ipsang**2))
257 10 : do ii=1,ipsang**2
258 12871 : ylm(:,ii) = ylm_tmp(ii,:)
259 : end do
260 1 : ABI_FREE(ylm_tmp)
261 : !
262 : ! In order to do spline fits, the |r-R| data must be sorted
263 : ! Here we sort the nrm points, keeping track of which goes where
264 2 : ABI_MALLOC(nrm_sort,(nfgd))
265 1431 : nrm_sort = nrm
266 :
267 3 : ABI_MALLOC(iperm,(nfgd))
268 1430 : do ifgd=1,nfgd
269 1430 : iperm(ifgd)=ifgd
270 : end do
271 :
272 1 : call sort_dp(nfgd,nrm_sort,iperm,tol8)
273 :
274 : ! Now make spline fits of phi and tphi onto the fine grid around the atom
275 4 : ABI_MALLOC(phigrd,(nfgd,ln_size))
276 3 : ABI_MALLOC(tphigrd,(nfgd,ln_size))
277 2 : ABI_MALLOC(ff,(nfgd))
278 1 : if (my_optgrad==1) then
279 0 : ABI_MALLOC(phigrd_gr,(nfgd,ln_size))
280 0 : ABI_MALLOC(tphigrd_gr,(nfgd,ln_size))
281 0 : ABI_MALLOC(gg,(mesh_size))
282 : end if
283 :
284 6 : do inl=1,ln_size
285 : !
286 : ! * splint phi onto points and reorder indices.
287 : call paw_splint(mesh_size,Pawrad(itypat)%rad,Pawtab(itypat)%phi(:,inl),&
288 5 : & Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
289 7150 : do ifgd=1,nfgd
290 7145 : ii=iperm(ifgd)
291 7150 : phigrd(ii,inl) = ff(ifgd)
292 : end do
293 : !
294 : ! * compute d phi/dr, interpolate onto points and reorder indices.
295 5 : if (my_optgrad==1) then
296 0 : ybcbeg=zero; ybcend=zero
297 0 : call nderiv_gen(gg,Pawtab(itypat)%phi(:,inl),Pawrad(itypat))
298 : call paw_spline(Pawrad(itypat)%rad,gg,mesh_size,ybcbeg,ybcend,&
299 0 : & Paw_lmn_spline(itypat)%phi(:,inl))
300 : call paw_splint(mesh_size,Pawrad(itypat)%rad,Paw_lmn_spline(itypat)%phi(:,inl),&
301 0 : & Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
302 0 : do ifgd=1,nfgd
303 0 : ii=iperm(ifgd)
304 0 : phigrd_gr(ii,inl) = ff(ifgd)
305 : end do
306 : end if
307 : !
308 : ! * compute d tphi/dr, interpolate onto points and reorder indices.
309 : call paw_splint(mesh_size,Pawrad(itypat)%rad,Pawtab(itypat)%tphi(:,inl),&
310 5 : & Paw_lmn_spline(itypat)%tphi(:,inl),nfgd,nrm_sort,ff)
311 7150 : do ifgd=1,nfgd
312 7145 : ii=iperm(ifgd)
313 7150 : tphigrd(ii,inl) = ff(ifgd)
314 : end do
315 6 : if (my_optgrad==1) then
316 0 : ybcbeg=zero; ybcend=zero
317 0 : call nderiv_gen(gg,Pawtab(itypat)%tphi(:,inl),Pawrad(itypat))
318 : call paw_spline(Pawrad(itypat)%rad,gg,mesh_size,ybcbeg,ybcend,&
319 0 : & Paw_lmn_spline(itypat)%tphi(:,inl))
320 : call paw_splint(mesh_size,Pawrad(itypat)%rad,Paw_lmn_spline(itypat)%tphi(:,inl),&
321 0 : & Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
322 0 : do ifgd=1,nfgd
323 0 : ii=iperm(ifgd)
324 0 : tphigrd_gr(ii,inl) = ff(ifgd)
325 : end do
326 : end if
327 : end do !inl
328 :
329 1 : ABI_FREE(ff)
330 1 : if (my_optgrad==1) then
331 0 : ABI_FREE(gg)
332 : end if
333 : !
334 : ! === Calculate AE and PS partial waves inside the sphere ===
335 : ! * recall that <r|phi>=u(r)*Slm(r^)/r, hence avoid division by zero except for s-waves.
336 4 : ABI_MALLOC(Paw_onsite(iatom1)%phi ,(nfgd,lmn_size))
337 3 : ABI_MALLOC(Paw_onsite(iatom1)%tphi,(nfgd,lmn_size))
338 :
339 1 : if (my_optgrad==1) then
340 0 : ABI_MALLOC(Paw_onsite(iatom1)%phi_gr ,(3,nfgd,lmn_size))
341 0 : ABI_MALLOC(Paw_onsite(iatom1)%tphi_gr,(3,nfgd,lmn_size))
342 : end if
343 :
344 14 : do jlmn=1,lmn_size
345 13 : jl = indlmn(1,jlmn)
346 13 : jm = indlmn(2,jlmn)
347 13 : jlm = indlmn(4,jlmn)
348 13 : jln = indlmn(5,jlmn)
349 :
350 18591 : do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere
351 : !if (nrm(ifgd)>tol16) then
352 18590 : if (nrm(ifgd)>tol10) then ! tol10 to be consistent with initylmr.
353 18564 : rR = nrm(ifgd) ! value of |r-R|
354 : !write(ab_out,*) 'rR:',rR,' phigrd:',phigrd(ifgd,jln),' tphigrd:',tphigrd(ifgd,jln),' ylm:',ylm(ifgd,jlm)
355 18564 : phj = phigrd (ifgd,jln)*ylm(ifgd,jlm)/rR
356 18564 : tphj= tphigrd(ifgd,jln)*ylm(ifgd,jlm)/rR
357 18564 : Paw_onsite(iatom1)%phi (ifgd,jlmn) = phj
358 18564 : Paw_onsite(iatom1)%tphi(ifgd,jlmn) = tphj
359 :
360 18564 : if (my_optgrad==1) then
361 : Paw_onsite(iatom1)%phi_gr (1:3,ifgd,jlmn) = phigrd (ifgd,jln)*ylm_gr(1:3,jlm,ifgd) &
362 0 : & + phigrd_gr(ifgd,jln)*local_pawfgrtab(iatom1)%rfgd(1:3,ifgd)*ylm(ifgd,jlm)
363 : Paw_onsite(iatom1)%tphi_gr (1:3,ifgd,jlmn) = tphigrd (ifgd,jln)*ylm_gr(1:3,jlm,ifgd) &
364 0 : & + tphigrd_gr(ifgd,jln)*local_pawfgrtab(iatom1)%rfgd(1:3,ifgd)*ylm(ifgd,jlm)
365 : end if
366 :
367 : else ! Extrapolate if the point is at the origin
368 13 : yvals(1) = zero
369 13 : if (jl==0) then
370 8 : yvals(2:4) = Pawtab(itypat)%phi(2:4,jln)/Pawrad(itypat)%rad(2:4)
371 2 : call pawrad_deducer0(yvals,4,Pawrad(itypat))
372 : end if
373 13 : Paw_onsite(iatom1)%phi(ifgd,jlmn) = yvals(1) * ylm(ifgd,jlm)
374 13 : yvals(1) = zero
375 13 : if (jl==0) then
376 8 : yvals(2:4) = Pawtab(itypat)%tphi(2:4,jln)/Pawrad(itypat)%rad(2:4)
377 2 : call pawrad_deducer0(yvals,4,pawrad(itypat))
378 : end if
379 13 : Paw_onsite(iatom1)%tphi(ifgd,jlmn) = yvals(1) * ylm(ifgd,jlm)
380 : ! The gradient is expected to go to zero at the origin
381 13 : if (my_optgrad==1) then
382 0 : Paw_onsite(iatom1)%phi_gr (1:3,ifgd,jlmn) = zero
383 0 : Paw_onsite(iatom1)%tphi_gr (1:3,ifgd,jlmn) = zero
384 : end if
385 : end if
386 :
387 : end do !nfgd
388 : end do !jlmn
389 :
390 1 : ABI_FREE(nrm)
391 1 : ABI_FREE(nrm_sort)
392 1 : ABI_FREE(iperm)
393 1 : ABI_FREE(phigrd)
394 1 : ABI_FREE(tphigrd)
395 1 : ABI_FREE(ylm)
396 1 : ABI_FREE(ylm_gr)
397 2 : if (my_optgrad==1) then
398 0 : ABI_FREE(phigrd_gr)
399 0 : ABI_FREE(tphigrd_gr)
400 : end if
401 : end do !iatom
402 : !
403 : !* Free 2nd derivates used for spline.
404 1 : call paw_pwaves_lmn_free(Paw_lmn_spline)
405 2 : ABI_FREE(Paw_lmn_spline)
406 :
407 : ! Destroy atom table used for parallelism
408 1 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
409 :
410 2 : end subroutine paw_pwaves_lmn_init
411 : !!***
412 :
413 : !----------------------------------------------------------------------
414 :
415 : !!****f* m_paw_pwaves_lmn/paw_pwaves_lmn_free
416 : !! NAME
417 : !! paw_pwaves_lmn_free
418 : !!
419 : !! FUNCTION
420 : !!
421 : !! INPUTS
422 : !!
423 : !! OUTPUT
424 : !!
425 : !! SOURCE
426 :
427 2 : subroutine paw_pwaves_lmn_free(Paw_onsite)
428 :
429 : !Arguments ------------------------------------
430 : !arrays
431 : type(paw_pwaves_lmn_t),intent(inout) :: Paw_onsite(:)
432 :
433 : !Local variables-------------------------------
434 : !scalars
435 : integer :: iatom
436 : ! *************************************************************************
437 :
438 : !@paw_pwaves_lmn_t
439 :
440 6 : do iatom=LBOUND(Paw_onsite,DIM=1),UBOUND(Paw_onsite,DIM=1)
441 2 : ABI_SFREE(Paw_onsite(iatom)%phi)
442 2 : ABI_SFREE(Paw_onsite(iatom)%tphi)
443 2 : ABI_SFREE(Paw_onsite(iatom)%r0shift)
444 2 : ABI_SFREE(Paw_onsite(iatom)%phi_gr)
445 4 : ABI_SFREE(Paw_onsite(iatom)%tphi_gr)
446 : end do
447 :
448 2 : end subroutine paw_pwaves_lmn_free
449 : !!***
450 :
451 : !----------------------------------------------------------------------
452 :
453 0 : END MODULE m_paw_pwaves_lmn
454 : !!***
|
