TABLE OF CONTENTS


ABINIT/m_paw_pwaves_lmn [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_pwaves_lmn

FUNCTION

  This module contains the definition of the paw_pwaves_lmn structured datatype,
  as well as related functions and methods.
  paw_pwaves_lmn is used to store the 3D values of the all-electron and of
  the pseudized part of the PAW partial waves on the set of FFT points falling
  inside the spheres around each atom.

COPYRIGHT

 Copyright (C) 2008-2024 ABINIT group (MG,MT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

NOTES

 * Routines tagged with "@type_name" are strongly connected to the definition of the data type.
   Strongly connected means that the proper functioning of the implementation relies on the
   assumption that the tagged procedure is consistent with the type declaration.
   Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
   that all the strongly connected routines are changed accordingly to accommodate the modification of the data type.
   Typical examples of strongly connected routines are creation, destruction or reset methods.

SOURCE

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

m_paw_pwaves_lmn/paw_pwaves_lmn_free [ Functions ]

[ Top ] [ m_paw_pwaves_lmn ] [ Functions ]

NAME

 paw_pwaves_lmn_free

FUNCTION

INPUTS

OUTPUT

SOURCE

429 subroutine paw_pwaves_lmn_free(Paw_onsite)
430 
431  implicit none
432 
433 !Arguments ------------------------------------
434 !arrays
435  type(paw_pwaves_lmn_t),intent(inout) :: Paw_onsite(:)
436 
437 !Local variables-------------------------------
438 !scalars
439  integer :: iatom
440 
441 ! *************************************************************************
442 
443 !@paw_pwaves_lmn_t
444 
445  do iatom=LBOUND(Paw_onsite,DIM=1),UBOUND(Paw_onsite,DIM=1)
446    if (allocated(Paw_onsite(iatom)%phi)) then
447      ABI_FREE(Paw_onsite(iatom)%phi)
448    end if
449    if (allocated(Paw_onsite(iatom)%tphi)) then
450      ABI_FREE(Paw_onsite(iatom)%tphi)
451    end if
452    if (allocated(Paw_onsite(iatom)%r0shift)) then
453      ABI_FREE(Paw_onsite(iatom)%r0shift)
454    end if
455    if (allocated(Paw_onsite(iatom)%phi_gr )) then
456      ABI_FREE(Paw_onsite(iatom)%phi_gr)
457    end if
458    if (allocated(Paw_onsite(iatom)%tphi_gr)) then
459      ABI_FREE(Paw_onsite(iatom)%tphi_gr)
460    end if
461  end do
462 
463 end subroutine paw_pwaves_lmn_free

m_paw_pwaves_lmn/paw_pwaves_lmn_init [ Functions ]

[ Top ] [ m_paw_pwaves_lmn ] [ Functions ]

NAME

 paw_pwaves_lmn_init

FUNCTION

INPUTS

  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms

OUTPUT

SOURCE

126 subroutine paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,ntypat,rprimd,xcart,Pawtab, &
127 &                              Pawrad,local_pawfgrtab,optgrad,&
128 &                              mpi_atmtab,comm_atom) ! optional arguments (parallelism)
129 
130  implicit none
131 
132 !Arguments ------------------------------------
133 !scalars
134  integer,intent(in) :: my_natom,natom,ntypat
135  integer,optional,intent(in) :: optgrad,comm_atom
136 !arrays
137  integer,optional,target,intent(in) :: mpi_atmtab(:)
138  real(dp),intent(in) :: xcart(3,natom),rprimd(3,3)
139  type(pawtab_type),target,intent(in) :: Pawtab(ntypat)
140  type(pawrad_type),intent(in) :: Pawrad(ntypat)
141  type(pawfgrtab_type),intent(in) :: local_pawfgrtab(my_natom)
142  type(paw_pwaves_lmn_t),intent(out) :: Paw_onsite(my_natom)
143 
144 !Local variables-------------------------------
145 !scalars
146  integer :: itypat,ln_size,lmn_size,mesh_size,inl,iatom,iatom1,my_comm_atom,my_optgrad
147  integer :: nfgd,ifgd,ipsang,option_ylmr,normchoice,ii,jlmn,jl,jm,jlm,jln
148  logical :: my_atmtab_allocated,paral_atom
149  real(dp) :: phj,rR,tphj,ybcbeg,ybcend
150 !arrays
151  integer, allocatable :: iperm(:)
152  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
153  integer,pointer :: my_atmtab(:)
154  real(dp) :: yvals(4),red(3),shift(3)
155  real(dp),allocatable :: ff(:),nrm(:),nrm_sort(:),phigrd(:,:),tphigrd(:,:),ylm_tmp(:,:),ylm(:,:),ylm_gr(:,:,:)
156  real(dp),allocatable :: rsph_red(:,:),rsph_cart(:,:),phigrd_gr(:,:),tphigrd_gr(:,:),gg(:)
157  type(paw_pwaves_lmn_t),allocatable :: Paw_lmn_spline(:)
158 
159 ! *************************************************************************
160 
161 !@paw_pwaves_lmn_t
162 
163  if (my_natom==0) return
164  my_optgrad = -1; if (present(optgrad)) my_optgrad = optgrad
165 
166  ABI_CHECK(all(local_pawfgrtab(:)%rfgd_allocated==1),"R vectors not allocated in pawfgrtab")
167 
168  ! Set up parallelism over atoms
169  paral_atom=(present(comm_atom).and.(my_natom/=natom))
170  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
171  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
172  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
173 
174  ! Prepare the spline. Calculate 2nd derivatives of partial waves for each atom type.
175  ABI_MALLOC(Paw_lmn_spline,(ntypat))
176 
177  do itypat=1,ntypat
178    ln_size  =Pawtab(itypat)%basis_size
179    mesh_size=Pawtab(itypat)%mesh_size
180 
181    ABI_MALLOC(Paw_lmn_spline(itypat)%phi ,(mesh_size,ln_size))
182    ABI_MALLOC(Paw_lmn_spline(itypat)%tphi,(mesh_size,ln_size))
183 
184    do inl=1,ln_size ! Calculate 2nd derivatives of %phi and %tphi for each ln component.
185      ybcbeg=zero; ybcend=zero
186      call paw_spline(Pawrad(itypat)%rad,Pawtab(itypat)%phi(:,inl), mesh_size,&
187 &                    ybcbeg,ybcend,Paw_lmn_spline(itypat)%phi(:,inl))
188 
189      ybcbeg=zero; ybcend=zero
190      call paw_spline(Pawrad(itypat)%rad,Pawtab(itypat)%tphi(:,inl),mesh_size,&
191 &                    ybcbeg,ybcend,Paw_lmn_spline(itypat)%tphi(:,inl))
192    end do
193  end do
194  !
195  ! === spline for each atom ===
196  ! * FFT points within PAW sphere depend on the atom site.
197  do iatom1=1,my_natom
198    iatom=iatom1;if (paral_atom) iatom=my_atmtab(iatom1)
199 
200    itypat   = local_pawfgrtab(iatom1)%itypat
201    ln_size  = Pawtab(itypat)%basis_size
202    lmn_size = Pawtab(itypat)%lmn_size
203    mesh_size= Pawrad(itypat)%mesh_size
204    nfgd     = local_pawfgrtab(iatom1)%nfgd ! no. of points in the fine grid for this PAW sphere
205    indlmn   => Pawtab(itypat)%indlmn
206 
207    ! The points in the PAW sphere might belong to a different unit cell, in this case one has to
208    ! reconstruct the contribution to the AE psi_k given by the lattice-symmetric atom in another sphere.
209    ! Here I wrap rr back into the first unit cell keeping trace of the lattice vector L
210    ! needed so that rr = rr_first_ucell + L
211    ! The contribution to the AE u(r) in the first unit cell has to be multiplied by e^{-ikL}.
212 
213    Paw_onsite(iatom1)%nfgd     = nfgd
214    Paw_onsite(iatom1)%lmn_size = lmn_size
215    !
216    ABI_MALLOC(Paw_onsite(iatom1)%r0shift,(3,nfgd))
217 
218    ABI_MALLOC(rsph_red,(3,nfgd))
219    ABI_MALLOC(rsph_cart,(3,nfgd))
220    do ifgd=1,nfgd
221      rsph_cart(:,ifgd) = local_pawfgrtab(iatom1)%rfgd(:,ifgd) + xcart(:,iatom)
222    end do
223    call xcart2xred(nfgd,rprimd,rsph_cart,rsph_red) ! go to reduced coordinates.
224    ABI_FREE(rsph_cart)
225 
226    do ifgd=1,nfgd
227      call wrap2_zero_one(rsph_red(1,ifgd),red(1),shift(1)) ! rr = r_cell + shift
228      call wrap2_zero_one(rsph_red(2,ifgd),red(2),shift(2))
229      call wrap2_zero_one(rsph_red(3,ifgd),red(3),shift(3))
230      Paw_onsite(iatom1)%r0shift(:,ifgd) = NINT(shift)
231      !if (ANY( ABS(shift) > tol12)) then
232        !ABI_WARNING("rmR_red is outside the first unit cell.")
233        !write(ab_out,*)rsph_red(:,ifgd),shift
234      !end if
235    end do
236    ABI_FREE(rsph_red)
237    !
238    ! * Obtain |r-R| on fine grid, note that rfgd is given in Cartesian coordinates.
239    ABI_MALLOC(nrm,(nfgd))
240    do ifgd=1,nfgd
241      nrm(ifgd) = sqrt(dot_product(local_pawfgrtab(iatom1)%rfgd(:,ifgd),local_pawfgrtab(iatom1)%rfgd(:,ifgd)))
242    end do
243    !
244    ! * Compute Ylm for each r-R vector.
245    ipsang = 1 + (Pawtab(itypat)%l_size-1)/2 ! recall l_size=2*l_max-1 where l_max is shifted by 1.
246    ABI_MALLOC(ylm_tmp,(ipsang**2,nfgd))
247    normchoice = 1 ! Use computed norms of input vectors.
248    if (my_optgrad==1) then
249      option_ylmr=2 ! Compute Ylm(r-R) and its gradient
250      ABI_MALLOC(ylm_gr,(3,ipsang**2,nfgd))
251    else
252      option_ylmr= 1 ! To compute Ylm(r-R).
253      ABI_MALLOC(ylm_gr,(3,3,0))
254    end if
255    call initylmr(ipsang,normchoice,nfgd,nrm,option_ylmr,local_pawfgrtab(iatom1)%rfgd,ylm_tmp,ylm_gr)
256    !
257    !  Exchange dimensions for better memory access.
258    ABI_MALLOC(ylm,(nfgd,ipsang**2))
259    do ii=1,ipsang**2
260      ylm(:,ii) = ylm_tmp(ii,:)
261    end do
262    ABI_FREE(ylm_tmp)
263    !
264    ! In order to do spline fits, the |r-R| data must be sorted
265    ! Here we sort the nrm points, keeping track of which goes where
266    ABI_MALLOC(nrm_sort,(nfgd))
267    nrm_sort = nrm
268 
269    ABI_MALLOC(iperm,(nfgd))
270    do ifgd=1,nfgd
271      iperm(ifgd)=ifgd
272    end do
273 
274    call sort_dp(nfgd,nrm_sort,iperm,tol8)
275 
276    ! Now make spline fits of phi and tphi  onto the fine grid around the atom
277    ABI_MALLOC(phigrd,(nfgd,ln_size))
278    ABI_MALLOC(tphigrd,(nfgd,ln_size))
279    ABI_MALLOC(ff,(nfgd))
280    if (my_optgrad==1) then
281      ABI_MALLOC(phigrd_gr,(nfgd,ln_size))
282      ABI_MALLOC(tphigrd_gr,(nfgd,ln_size))
283      ABI_MALLOC(gg,(mesh_size))
284    end if
285 
286    do inl=1,ln_size
287      !
288      ! * splint phi onto points and reorder indices.
289      call paw_splint(mesh_size,Pawrad(itypat)%rad,Pawtab(itypat)%phi(:,inl),&
290 &                    Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
291      do ifgd=1,nfgd
292        ii=iperm(ifgd)
293        phigrd(ii,inl) = ff(ifgd)
294      end do
295      !
296      ! * compute d phi/dr, interpolate onto points and reorder indices.
297      if (my_optgrad==1) then
298        ybcbeg=zero; ybcend=zero
299        call nderiv_gen(gg,Pawtab(itypat)%phi(:,inl),Pawrad(itypat))
300        call paw_spline(Pawrad(itypat)%rad,gg,mesh_size,ybcbeg,ybcend,&
301 &                      Paw_lmn_spline(itypat)%phi(:,inl))
302        call paw_splint(mesh_size,Pawrad(itypat)%rad,Paw_lmn_spline(itypat)%phi(:,inl),&
303 &                      Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
304        do ifgd=1,nfgd
305          ii=iperm(ifgd)
306          phigrd_gr(ii,inl)  = ff(ifgd)
307        end do
308      end if
309      !
310      ! * compute d tphi/dr, interpolate onto points and reorder indices.
311      call paw_splint(mesh_size,Pawrad(itypat)%rad,Pawtab(itypat)%tphi(:,inl),&
312 &                    Paw_lmn_spline(itypat)%tphi(:,inl),nfgd,nrm_sort,ff)
313      do ifgd=1,nfgd
314        ii=iperm(ifgd)
315        tphigrd(ii,inl) = ff(ifgd)
316      end do
317      if (my_optgrad==1) then
318        ybcbeg=zero; ybcend=zero
319        call nderiv_gen(gg,Pawtab(itypat)%tphi(:,inl),Pawrad(itypat))
320        call paw_spline(Pawrad(itypat)%rad,gg,mesh_size,ybcbeg,ybcend,&
321 &                      Paw_lmn_spline(itypat)%tphi(:,inl))
322        call paw_splint(mesh_size,Pawrad(itypat)%rad,Paw_lmn_spline(itypat)%tphi(:,inl),&
323 &                      Paw_lmn_spline(itypat)%phi(:,inl),nfgd,nrm_sort,ff)
324        do ifgd=1,nfgd
325          ii=iperm(ifgd)
326          tphigrd_gr(ii,inl)  = ff(ifgd)
327        end do
328      end if
329    end do !inl
330 
331    ABI_FREE(ff)
332    if (my_optgrad==1) then
333      ABI_FREE(gg)
334    end if
335    !
336    ! === Calculate AE and PS partial waves inside the sphere ===
337    ! * recall that <r|phi>=u(r)*Slm(r^)/r, hence avoid division by zero except for s-waves.
338    ABI_MALLOC(Paw_onsite(iatom1)%phi ,(nfgd,lmn_size))
339    ABI_MALLOC(Paw_onsite(iatom1)%tphi,(nfgd,lmn_size))
340 
341    if (my_optgrad==1) then
342      ABI_MALLOC(Paw_onsite(iatom1)%phi_gr ,(3,nfgd,lmn_size))
343      ABI_MALLOC(Paw_onsite(iatom1)%tphi_gr,(3,nfgd,lmn_size))
344    end if
345 
346    do jlmn=1,lmn_size
347      jl  = indlmn(1,jlmn)
348      jm  = indlmn(2,jlmn)
349      jlm = indlmn(4,jlmn)
350      jln = indlmn(5,jlmn)
351 
352      do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere
353        !if (nrm(ifgd)>tol16) then
354        if (nrm(ifgd)>tol10) then ! tol10 to be consistent with initylmr.
355          rR  = nrm(ifgd) ! value of |r-R|
356          !write(ab_out,*) 'rR:',rR,' phigrd:',phigrd(ifgd,jln),' tphigrd:',tphigrd(ifgd,jln),' ylm:',ylm(ifgd,jlm)
357          phj = phigrd (ifgd,jln)*ylm(ifgd,jlm)/rR
358          tphj= tphigrd(ifgd,jln)*ylm(ifgd,jlm)/rR
359          Paw_onsite(iatom1)%phi (ifgd,jlmn) =  phj
360          Paw_onsite(iatom1)%tphi(ifgd,jlmn) = tphj
361 
362          if (my_optgrad==1) then
363            Paw_onsite(iatom1)%phi_gr (1:3,ifgd,jlmn) = phigrd (ifgd,jln)*ylm_gr(1:3,jlm,ifgd) &
364 &            + phigrd_gr(ifgd,jln)*local_pawfgrtab(iatom1)%rfgd(1:3,ifgd)*ylm(ifgd,jlm)
365            Paw_onsite(iatom1)%tphi_gr (1:3,ifgd,jlmn) = tphigrd (ifgd,jln)*ylm_gr(1:3,jlm,ifgd) &
366 &            + tphigrd_gr(ifgd,jln)*local_pawfgrtab(iatom1)%rfgd(1:3,ifgd)*ylm(ifgd,jlm)
367          end if
368 
369        else ! Extrapolate if the point is at the origin
370          yvals(1) = zero
371          if (jl==0) then
372            yvals(2:4) = Pawtab(itypat)%phi(2:4,jln)/Pawrad(itypat)%rad(2:4)
373            call pawrad_deducer0(yvals,4,Pawrad(itypat))
374          end if
375          Paw_onsite(iatom1)%phi(ifgd,jlmn) = yvals(1) * ylm(ifgd,jlm)
376          yvals(1) = zero
377          if (jl==0) then
378            yvals(2:4) = Pawtab(itypat)%tphi(2:4,jln)/Pawrad(itypat)%rad(2:4)
379            call pawrad_deducer0(yvals,4,pawrad(itypat))
380          end if
381          Paw_onsite(iatom1)%tphi(ifgd,jlmn) = yvals(1) * ylm(ifgd,jlm)
382          ! The gradient is expected to go to zero at the origin
383          if (my_optgrad==1) then
384            Paw_onsite(iatom1)%phi_gr (1:3,ifgd,jlmn) = zero
385            Paw_onsite(iatom1)%tphi_gr (1:3,ifgd,jlmn) = zero
386          end if
387        end if
388 
389      end do !nfgd
390    end do !jlmn
391 
392    ABI_FREE(nrm)
393    ABI_FREE(nrm_sort)
394    ABI_FREE(iperm)
395    ABI_FREE(phigrd)
396    ABI_FREE(tphigrd)
397    ABI_FREE(ylm)
398    ABI_FREE(ylm_gr)
399    if (my_optgrad==1) then
400      ABI_FREE(phigrd_gr)
401      ABI_FREE(tphigrd_gr)
402    end if
403  end do !iatom
404  !
405  !* Free 2nd derivates used for spline.
406  call paw_pwaves_lmn_free(Paw_lmn_spline)
407  ABI_FREE(Paw_lmn_spline)
408 
409  ! Destroy atom table used for parallelism
410  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
411 
412 end subroutine paw_pwaves_lmn_init

m_paw_pwaves_lmn/paw_pwaves_lmn_t [ Types ]

[ Top ] [ m_paw_pwaves_lmn ] [ Types ]

NAME

 paw_pwaves_lmn_t

FUNCTION

  Datatype used to store the 3D values of the all-electron and of the pseudized part of the
  PAW partial waves on the set of FFT points falling inside the spheres around each atom.
  The data is mainly used for plotting the true PAW wavefunctions in real space.

SOURCE

 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