TABLE OF CONTENTS


ABINIT/m_paw_mkaewf [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_mkaewf

FUNCTION

 Construct complete AE wave functions on the fine FFT grid adding onsite PAW corrections.

COPYRIGHT

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

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_paw_mkaewf
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_abicore
34  use m_xmpi
35  use m_hide_blas
36  use m_splines
37  use m_errors
38  use m_nctk
39  use m_hdr
40 #ifdef HAVE_NETCDF
41  use netcdf
42 #endif
43 
44  use m_io_tools,       only : flush_unit
45  use m_numeric_tools,  only : wrap2_zero_one
46  use m_fftcore,        only : sphereboundary
47  use m_geometry,       only : xcart2xred
48  use m_crystal,        only : crystal_t
49  use m_crystal_io,     only : crystal_ncwrite
50  use m_ebands,         only : ebands_ncwrite
51  use m_pawrad,         only : pawrad_type
52  use m_pawtab,         only : pawtab_type, pawtab_get_lsize
53  use m_pawfgrtab,      only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free, pawfgrtab_print
54  use m_pawcprj,        only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free
55  use m_paw_pwaves_lmn, only : paw_pwaves_lmn_t, paw_pwaves_lmn_init, paw_pwaves_lmn_free
56  use m_paral_atom,     only : get_my_atmtab, free_my_atmtab
57  use m_paw_nhat,       only : nhatgrid
58  use m_mpinfo,         only : proc_distrb_cycle
59  use m_fft,            only : fourwf
60 
61  implicit none
62 
63  private
64 
65  public :: pawmkaewf
66 
67 CONTAINS  !========================================================================================

m_paw_mkaewf/pawmkaewf [ Functions ]

[ Top ] [ m_paw_mkaewf ] [ Functions ]

NAME

 pawmkaewf

FUNCTION

 Construct complete AE wave functions on the fine FFT grid adding onsite PAW corrections.

INPUTS

 crystal<crystal_t>=Crystalline structure
 ebands<ebands_t>=Electronic energies
 dimcprj(natom)=array of dimensions of array cprj (not ordered)
 mband=maximum number of bands
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
 mkmem=number of k points treated by this node.
 mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
 [comm_atom]= MPI communicator over atoms
 mpw=maximum dimensioned size of npw.
 my_natom=number of atoms treated by current processor
 natom=number of atoms in cell
 ntypat=number of types of atoms in the cell
 nkpt=Total number of k-points
 nsppol=1 for unpolarized, 2 for spin-polarized
 unks=unit number for G vectors.
 nband(nkpt*nsppol)=Number of bands for each k-point and spin.
 istwfk(nkpt)=Storage mode at each k-point.
 Pawfgrtab(natom) <type(pawfgrtab_type)> : data about the fine grid around each atom
 Pawrad(ntypat) <type(pawrad_type)> : radial mesh data for each type of atom
 Pawtab(ntypat) <type(pawtab_type)> : PAW functions around each type of atom
 Dtfil <type(datafiles_type)>=variables related to files
 cg(2,mcg)=planewave coefficients of wavefunctions.
 Cprj(natom,nspinor*mband*mkmem*nsppol)=<p_lmn|Cnk> coefficients for each WF |Cnk>
   and each |p_lmn> non-local projector
 npwarr(nkpt)=Number of plane waves at each k-point
 ngfftf(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  Note that ngfftf refers to the fine mesh.
 kg(3,mpw*mkmem)=reduced planewave coordinates
 Hdr<hdr_type>=the header of wf, den and pot files
 kpt(3,nkpt)=reduced coordinates of k points.

OUTPUT

  ierr=Status error
  Main output is written on file (ETSF_IO file format).

NOTES

 In PAW calculations, the pseudized wavefunction us represented
 on a relatively small plane wave basis set and is not normalized
 as it does not include the on-site PAW contributions which is described
 in terms of real spherical harmonics and radial functions.
 For post-processing and proper visualization, it is necessary
 to use the full electronic wave function, which is what this subroutine constructs.
 Specifically, it computes the pseudo part by doing an FFT from G- to r-space
 using the dense mesh defined by pawecutdg. The on-site PAW terms are also
 computed in real space inside each sphere and added to the pseudo part.
 Notice that this formula is expressed on the fine grid, and requires
 interpolating the PAW radial functions onto this grid, as well as calling
 initylmr in order to get the angular functions on the grid points.

PARENTS

      outscfcv

CHILDREN

      flush_unit,fourwf,free_my_atmtab,get_my_atmtab,nhatgrid
      paw_pwaves_lmn_free,paw_pwaves_lmn_init,pawcprj_alloc,pawcprj_free
      pawfgrtab_free,pawfgrtab_init,pawfgrtab_print,pawtab_get_lsize
      sphereboundary,wrap2_zero_one,wrtout,xcart2xred,xmpi_barrier,xmpi_max
      xmpi_sum

SOURCE

140 subroutine pawmkaewf(Dtset,crystal,ebands,my_natom,mpw,mband,mcg,mcprj,nkpt,mkmem,nsppol,nband,&
141 & istwfk,npwarr,kpt,ngfftf,kg,dimcprj,Pawfgrtab,Pawrad,Pawtab,&
142 & Hdr,Dtfil,cg,Cprj,MPI_enreg,ierr,pseudo_norms,set_k,set_band , &
143 & mpi_atmtab,comm_atom) ! Optional arguments
144 
145 
146 !This section has been created automatically by the script Abilint (TD).
147 !Do not modify the following lines by hand.
148 #undef ABI_FUNC
149 #define ABI_FUNC 'pawmkaewf'
150 !End of the abilint section
151 
152  implicit none
153 
154 !Arguments ------------------------------------
155 !scalars
156  integer,intent(in) :: my_natom,mband,mcg,mcprj,mkmem,mpw,nsppol,nkpt
157  integer,intent(in),optional :: comm_atom,set_k,set_band
158  integer,intent(out) :: ierr
159  type(Datafiles_type),intent(in) :: Dtfil
160  type(MPI_type),intent(in) :: MPI_enreg
161  type(hdr_type),intent(inout) :: Hdr
162  type(dataset_type),intent(in) :: Dtset
163  type(crystal_t),intent(in) :: crystal
164  type(ebands_t),intent(in) :: ebands
165 !arrays
166  integer,intent(in) :: nband(nkpt*nsppol),istwfk(nkpt),npwarr(nkpt),dimcprj(crystal%natom)
167  integer,intent(in) :: ngfftf(18),kg(3,mpw*mkmem)
168  integer,optional,target,intent(in) :: mpi_atmtab(:)
169  real(dp),intent(in) :: cg(2,mcg)
170  real(dp),intent(in) :: kpt(3,nkpt)
171  real(dp),optional,intent(out) :: pseudo_norms(nsppol,nkpt,mband)
172  type(pawfgrtab_type),intent(in) :: Pawfgrtab(my_natom)
173  type(pawrad_type),intent(in) :: Pawrad(crystal%ntypat)
174  type(pawtab_type),intent(in) :: Pawtab(crystal%ntypat)
175  type(pawcprj_type),intent(in) :: Cprj(crystal%natom,mcprj)
176 
177 !Local variables-------------------------------
178 !scalars
179  integer,parameter :: tim_fourwf0=0,tim_rwwf0=0,master=0
180  integer :: bdtot_index,iband,icg,mgfftf,paral_kgb
181  integer :: iatom,iatom_tot,ifgd,ifftsph,ifft,itypat,ispinor,ipw,ndat,ii,i1,i2,i3
182  integer :: jl,jm,jlmn,natom
183  integer :: max_nfgd,nfgd,ln_size,lmn_size,my_comm_atom,option
184  integer :: iorder_cprj,comm_cell,me_kpt,ibsp,ibg,isppol,ikpt,nband_k,cplex
185  integer :: n1,n2,n3,n4,n5,n6,ikg,npwout,istwf_k,npw_k
186  integer :: nfftot,nprocs
187  integer :: optcut,optgr0,optgr1,optgr2,optrad,start_band,start_kpt,stop_kpt,stop_band
188  logical :: my_atmtab_allocated,paral_atom
189  real(dp),parameter :: weight1=one
190  real(dp) :: phj,tphj,re_p,im_p,norm,norm_rerr,max_rerr,imur,reur,arg
191  character(len=500) :: msg
192  character(len=nctk_slen) :: shape_str
193 !arrays
194  integer,allocatable :: l_size_atm(:)
195  integer, pointer :: my_atmtab(:)
196  integer,allocatable :: gbound(:,:),kg_k(:,:)
197  real(dp) :: red(3),shift(3),rfft(3),kpoint(3),cp_fact(2)
198  real(dp),allocatable :: r0shift(:,:,:),phk_atm(:,:,:)
199  real(dp),allocatable :: buf_tmp(:,:,:),fofgin(:,:),fofgin_down(:,:),fofgout(:,:)
200  real(dp),allocatable :: denpot(:,:,:),fofr(:,:,:,:),fofr_down(:,:,:,:),phkr(:,:)
201  real(dp),allocatable :: ur_ae(:,:), ur_pw(:,:),ur_ae_onsite(:,:),ur_ps_onsite(:,:)
202  real(dp),allocatable :: ur_mask(:),dummy_1d(:),rsph_red(:,:),rsph_cart(:,:)
203  type(pawcprj_type),allocatable :: Cprj_k(:,:)
204  type(pawfgrtab_type) :: local_pawfgrtab(my_natom)
205  type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:)
206 #ifdef HAVE_NETCDF
207  integer :: fform,ncerr,ncid,ae_ncid,pw_ncid,aeons_ncid,psons_ncid
208  character(len=fnlen) :: fname
209 #endif
210 
211 ! ************************************************************************
212 
213  DBG_ENTER("COLL")
214 
215 !Init parallelism
216  comm_cell = MPI_enreg%comm_cell; nprocs = xmpi_comm_size(comm_cell)
217  me_kpt = MPI_enreg%me_kpt; paral_kgb=mpi_enreg%paral_kgb
218 
219 !Compatibility tests
220  ABI_CHECK(mkmem/=0, "mkmem==0 not supported anymore!")
221  ABI_CHECK(MPI_enreg%paral_kgb == 0, "paral_kgb/=0 not coded")
222  ABI_CHECK(SIZE(dimcprj)>0, "dimcprj should be allocated")
223  ABI_CHECK(mpi_enreg%paral_spinor==0, "parallelisation over spinors not implemented")
224  ABI_CHECK(nprocs==1, "k spin parallelism not yet active")
225  ABI_CHECK(dtset%nspinor==1, "nspinor == 2 is buggy")
226 
227  natom = crystal%natom
228 
229 !Set up parallelism over atoms
230  paral_atom=(present(comm_atom).and.(my_natom/=natom))
231  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
232  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
233  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
234 
235 !If collection of pseudo norms is enabled, make sure the array is initialised
236  if (present(pseudo_norms)) pseudo_norms = zero
237 
238 !use a local copy of pawfgrtab to make sure we use the correction in the paw spheres
239 !the usual pawfgrtab uses r_shape which may not be the same as r_paw
240  if (my_natom>0) then
241    if (paral_atom) then
242      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat,mpi_atmtab=my_atmtab)
243      call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat,&
244 &     mpi_atmtab=my_atmtab,comm_atom=my_comm_atom)
245    else
246      call pawtab_get_lsize(pawtab,l_size_atm,my_natom,Dtset%typat)
247      call pawfgrtab_init(local_pawfgrtab,Pawfgrtab(1)%cplex,l_size_atm,Dtset%nspden,Dtset%typat)
248    end if
249    ABI_FREE(l_size_atm)
250  end if
251  optcut = 1 ! use rpaw to construct local_pawfgrtab
252  optgr0 = 0; optgr1 = 0; optgr2 = 0 ! dont need gY terms locally
253  optrad = 1 ! do store r-R
254 
255  if (paral_atom) then
256    call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,&
257 &   optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred,&
258 &   comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
259  else
260    call nhatgrid(crystal%atindx1,crystal%gmet,my_natom,natom,crystal%nattyp,ngfftf,crystal%ntypat,&
261 &   optcut,optgr0,optgr1,optgr2,optrad,local_pawfgrtab,pawtab,crystal%rprimd,Dtset%typat,crystal%ucvol,Hdr%xred)
262  end if
263 !now local_pawfgrtab is ready to use
264 
265  max_nfgd=MAXVAL(local_pawfgrtab(:)%nfgd) ! MAX no. of points in the fine grid for this PAW sphere
266  ABI_MALLOC(r0shift,(3,max_nfgd,my_natom))
267  ABI_MALLOC(phk_atm,(2,max_nfgd,my_natom))
268 
269  do iatom=1,my_natom
270    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
271 
272    nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
273    ABI_MALLOC(rsph_red,(3,nfgd))
274    ABI_MALLOC(rsph_cart,(3,nfgd))
275    do ifgd=1,nfgd
276      rsph_cart(:,ifgd) = local_pawfgrtab(iatom)%rfgd(:,ifgd) + crystal%xcart(:,iatom_tot)
277    end do
278    call xcart2xred(nfgd,crystal%rprimd,rsph_cart,rsph_red) ! we work in reduced coordinates.
279    do ifgd=1,nfgd
280      call wrap2_zero_one(rsph_red(1,ifgd),red(1),shift(1)) ! num = red + shift
281      call wrap2_zero_one(rsph_red(2,ifgd),red(2),shift(2))
282      call wrap2_zero_one(rsph_red(3,ifgd),red(3),shift(3))
283      r0shift(:,ifgd,iatom) = shift
284      !if (ANY( ABS(shift) > tol12)) then
285      !  MSG_WARNING("rmR_red is outside the first unit cell.")
286      !  write(std_out,*)rsph_red(:,ifgd),shift
287      !end if
288    end do
289    ABI_FREE(rsph_red)
290    ABI_FREE(rsph_cart)
291  end do
292 
293  if (.not.paral_atom .and. my_natom>0) then
294    call pawfgrtab_print(local_pawfgrtab,natom=natom,unit=std_out,&
295 &   prtvol=Dtset%prtvol,mode_paral="COLL")
296  end if
297 
298  ierr=0
299 #ifndef HAVE_NETCDF
300  ierr = -1
301  write(msg,'(3a)')&
302 & "netcdf support must be enabled in order to output AE PAW wavefunction. ",ch10,&
303 & "No output will be produced, use --enable-netcdf at configure-time. "
304  MSG_WARNING(msg)
305  return
306 !These statements are necessary to avoid the compiler complain about unused variables:
307  ii=Dtset%usepaw;ii=Dtfil%unpaw;ii=Hdr%usepaw
308 #endif
309 
310 !FIXME check ordering in cprj and Eventually in external file
311 !why is iorder_cprj not stored in the file for crosschecking purpose?
312 !Here Im assuming cprj are not ordered!
313  iorder_cprj=0
314 
315 !n4,n5,n6 are FFT dimensions, modified to avoid cache trashing
316  n1=ngfftf(1); n2=ngfftf(2); n3=ngfftf(3)
317  n4=ngfftf(4); n5=ngfftf(5); n6=ngfftf(6)
318  nfftot=PRODUCT(ngfftf(1:3))
319  mgfftf=MAXVAL(ngfftf(1:3))
320 
321  ABI_MALLOC(phkr,(2,nfftot))
322  ABI_MALLOC(gbound,(2*mgfftf+8,2))
323 
324 #ifdef HAVE_NETCDF
325 !=== Initialize ETSF_IO files ===
326 ! FIXME: nspinor == 2 is buggy
327 
328  fname = trim(dtfil%filnam_ds(4))//'_PAWAVES.nc'
329  write(msg,'(2a)')' Opening file for AE PAW wave functions: ',trim(fname)
330  call wrtout(std_out,msg,'PERS')
331  call wrtout(ab_out,msg,'PERS')
332 
333  if (xmpi_comm_rank(comm_cell) == master) then
334    NCF_CHECK(nctk_open_create(ncid, fname, xmpi_comm_self))
335 
336    fform=602
337    NCF_CHECK(hdr_ncwrite(hdr, ncid, fform, nc_define=.True.))
338 
339    ! Define wavefunctions in real space on the dense FFT mesh
340    ! Fortran layout:
341    !real_space_wavefunctions: double 8d array with shape:
342    !  [real_or_complex_wavefunctions]
343    !  [number_of_grid_points_vector1][number_of_grid_points_vector2][number_of_grid_points_vector3]
344    !  [number_of_spinor_components]
345    !  [max_number_of_states][number_of_kpoints][number_of_spins]
346 
347    ncerr = nctk_def_dims(ncid, [ &
348    nctkdim_t("real_or_complex_wavefunctions", 2),  &
349    nctkdim_t("number_of_grid_points_vector1", n1), &
350    nctkdim_t("number_of_grid_points_vector2", n2), &
351    nctkdim_t("number_of_grid_points_vector3", n3)  &
352    ], defmode=.True.)
353    NCF_CHECK(ncerr)
354 
355    shape_str = "real_or_complex_wavefunctions, &
356 &   number_of_grid_points_vector1, number_of_grid_points_vector2, number_of_grid_points_vector3, &
357 &   number_of_spinor_components, &
358 &   max_number_of_states, number_of_kpoints, number_of_spins"
359 
360    ! Define wavefunctions in real space.
361    ncerr = nctk_def_arrays(ncid, [&
362    nctkarr_t('ur_ae', "dp", shape_str),&
363    nctkarr_t('ur_pw', "dp", shape_str),&
364    nctkarr_t('ur_ae_onsite', "dp", shape_str),&
365    nctkarr_t('ur_ps_onsite', "dp", shape_str) &
366    ], defmode=.True.)
367    NCF_CHECK(ncerr)
368 
369    ! Complete the geometry information.
370    NCF_CHECK(crystal_ncwrite(crystal, ncid))
371    NCF_CHECK(ebands_ncwrite(ebands, ncid))
372 
373    ! Close the file.
374    NCF_CHECK(nf90_close(ncid))
375  end if
376 
377  call xmpi_barrier(comm_cell)
378 
379  ! Reopen the file in parallel inside comm_cell
380  ! Note that we use individual IO thus there's no need to handle idle processes
381  ! if paral_kgb == 0 and nprocs > nkpt * nsppol
382  NCF_CHECK(nctk_open_modify(ncid, fname, comm_cell))
383  ae_ncid = nctk_idname(ncid, "ur_ae")
384  pw_ncid = nctk_idname(ncid, "ur_pw")
385  aeons_ncid = nctk_idname(ncid, "ur_ae_onsite")
386  psons_ncid = nctk_idname(ncid, "ur_ps_onsite")
387 
388  NCF_CHECK(nctk_set_datamode(ncid))
389 #endif
390 
391 !Init structure storing phi_{nlm} and tphi_(nlm} on the dense FFT points located in the PAW spheres.
392  ABI_DT_MALLOC(Paw_onsite,(natom))
393  if (paral_atom) then
394    call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,&
395    Pawtab,Pawrad,local_pawfgrtab, comm_atom=my_comm_atom,mpi_atmtab=my_atmtab)
396  else
397    call paw_pwaves_lmn_init(Paw_onsite,my_natom,natom,crystal%ntypat,crystal%rprimd,crystal%xcart,&
398    Pawtab,Pawrad,local_pawfgrtab)
399  end if
400 
401  bdtot_index=0; icg=0; ibg=0; norm_rerr=smallest_real
402 
403  ! === Loop over spin ===
404  do isppol=1,nsppol
405    ikg=0; start_kpt=1; stop_kpt=nkpt
406 
407    ! Check if k-point was specified (only serial)
408    if (present(set_k) .and. nprocs==1) then
409      if (set_k/=0) then
410        start_kpt = set_k
411        stop_kpt = set_k
412        !MSG_ERROR("set_k")
413      end if
414    end if
415 
416    ! === Loop over k points ===
417    do ikpt=start_kpt,stop_kpt
418      kpoint  = kpt(:,ikpt)
419      nband_k = nband(ikpt+(isppol-1)*nkpt)
420      npw_k   = npwarr(ikpt)
421      istwf_k = istwfk(ikpt)
422 
423      if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then
424        bdtot_index=bdtot_index+nband_k
425        !MSG_ERROR("cycle in seq!")
426        cycle
427      end if
428 
429      do i3=0,n3-1
430        rfft(3)=DBLE(i3)/n3
431        do i2=0,n2-1
432          rfft(2)=DBLE(i2)/n2
433          do i1=0,n1-1
434            rfft(1)=DBLE(i1)/n1
435            ifft = 1 +i1 +i2*n1 +i3*n1*n2
436            phkr(1,ifft) = COS(two_pi*dot_product(kpoint,rfft))
437            phkr(2,ifft) = SIN(two_pi*dot_product(kpoint,rfft))
438          end do
439        end do
440      end do
441      ! phkr(1,:)=one; phkr(2,:)=zero
442 
443 !    Calculate the phase for the onsite PAW contributions.
444      do iatom=1,my_natom
445        nfgd=local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
446        do ifgd=1,nfgd
447          arg = -two_pi* dot_product(r0shift(:,ifgd,iatom),kpoint)
448          phk_atm(1,ifgd,iatom) = COS(arg)
449          phk_atm(2,ifgd,iatom) = SIN(arg)
450        end do
451      end do
452 
453      ABI_DT_MALLOC(Cprj_k,(natom,dtset%nspinor*nband_k))
454      call pawcprj_alloc(Cprj_k,0,dimcprj)
455 
456 !    Extract cprj for this k-point.
457      ibsp=0
458      do iband=1,nband_k
459        do ispinor=1,dtset%nspinor
460          ibsp=ibsp+1
461          do iatom=1,natom
462            Cprj_k(iatom,ibsp)%cp(:,:)=Cprj(iatom,ibsp+ibg)%cp(:,:)
463          end do
464        end do
465      end do
466 
467      ABI_MALLOC(kg_k,(3,npw_k))
468 
469      ! Extract G-vectors.
470      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
471      call sphereboundary(gbound,istwf_k,kg_k,mgfftf,npw_k)
472 
473      ! If a single band is requested, neuter the loop (only serial)
474      start_band = 1; stop_band = nband_k
475      if (present(set_band).AND.nprocs==1) then
476        if (set_band/=0) then
477          start_band = set_band
478          stop_band = set_band
479          !MSG_ERROR("set_band")
480        end if
481      end if
482 
483      ! Loop over bands.
484      do iband=start_band,stop_band
485 
486        ! Fourier transform on the real fft box of the smooth part.
487        ndat=Dtset%nspinor
488        ABI_MALLOC(fofgin,(2,npw_k*ndat))
489        ABI_MALLOC(fofr,(2,n4,n5,n6*ndat))
490 
491        do ipw=1,npw_k*dtset%nspinor
492          fofgin(:,ipw)=cg(:,ipw+(iband-1)*npw_k*dtset%nspinor+icg)
493        end do
494 
495        ! Complex can be set to 0 with this option(0) of fourwf
496        option=0; cplex=0; npwout=1
497        ABI_MALLOC(denpot,(cplex*n4,n5,n6))
498        ABI_MALLOC(fofgout,(2,npwout*ndat))
499 
500        call fourwf(cplex,denpot,fofgin(:,1:npw_k),fofgout,fofr(:,:,:,1:n6),gbound,gbound,istwf_k,kg_k,kg_k,&
501 &       mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,paral_kgb,tim_fourwf0,weight1,weight1,&
502 &       use_gpu_cuda=Dtset%use_gpu_cuda)
503 
504 !      Here I do not know if fourwf works in the case of spinors,
505 !      It seems that not all fftalg option support ndata! should check!
506 !      Do not forget to declare real(dp)::fofgin_down(:,:) to use the following statements
507        if (Dtset%nspinor==2) then
508          ABI_MALLOC(fofgin_down,(2,npw_k))
509          ABI_MALLOC(fofr_down,(2,n4,n5,n6))
510          fofgin_down(:,:)=fofgin(:,1+npw_k:2*npw_k)
511 !        Complex can be set to 0 with this option(0) of fourwf
512 !        cplex=1; option=1; npwout=1; ndat=1
513 !        NOTE: fofr_down can NOT be replaced by fofr(:,:,:,n6+1:2*n6), or else
514 !        the data in fofr(:,:,:,1:n6) will be the same with fofr(:,:,:,n6+1:2*n6)
515          call fourwf(cplex,denpot,fofgin_down,fofgout,fofr_down,gbound,gbound,istwf_k,kg_k,kg_k,&
516 &         mgfftf,MPI_enreg,1,ngfftf,npw_k,npwout,n4,n5,n6,option,paral_kgb,tim_fourwf0,weight1,weight1)
517          ABI_FREE(fofgin_down)
518        end if
519 
520        ABI_MALLOC(ur_ae,(2,n1*n2*n3*ndat))
521        ABI_MALLOC(ur_ae_onsite,(2,n1*n2*n3))
522        ABI_MALLOC(ur_ps_onsite,(2,n1*n2*n3))
523        ABI_MALLOC(ur_pw,(2,n1*n2*n3*ndat))
524        ABI_MALLOC(ur_mask,(n1*n2*n3))
525 
526        ur_ae=zero;ur_ae_onsite=zero;ur_ps_onsite=zero;ur_pw=zero;ur_mask=zero
527 
528        ! * Add phase e^{ikr} since it is contained in cprj.
529        do i3=1,n3
530          do i2=1,n2
531            do i1=1,n1
532              ii = i1 + n1*(i2-1)+ n1*n2*(i3-1)
533              ur_pw(:,ii)=fofr(:,i1,i2,i3) ! Save pw part separately without the phase.
534              ur_ae(1,ii)= fofr(1,i1,i2,i3) * phkr(1,ii) - fofr(2,i1,i2,i3) * phkr(2,ii)
535              ur_ae(2,ii)= fofr(1,i1,i2,i3) * phkr(2,ii) + fofr(2,i1,i2,i3) * phkr(1,ii)
536              if(Dtset%nspinor==2) then
537                ur_pw(:,ii+n1*n2*n3)=fofr_down(:,i1,i2,i3) ! Save pw part separately without the phase.
538                ur_ae(1,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(1,ii) - fofr_down(2,i1,i2,i3) * phkr(2,ii)
539                ur_ae(2,ii+n1*n2*n3)= fofr_down(1,i1,i2,i3) * phkr(2,ii) + fofr_down(2,i1,i2,i3) * phkr(1,ii)
540              end if
541            end do
542          end do
543        end do
544        ABI_FREE(fofr)
545 
546        if(Dtset%nspinor==2) then
547          ABI_FREE(fofr_down)
548        end if
549 
550        ! === Add onsite term on the augmented FFT mesh ===
551        do iatom=1,my_natom
552          itypat  =local_pawfgrtab(iatom)%itypat
553          lmn_size=Pawtab(itypat)%lmn_size
554          ln_size =Pawtab(itypat)%basis_size   ! no. of nl elements in PAW basis
555          nfgd    =local_pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere
556 
557          ibsp=(iband-1)*dtset%nspinor
558          do ispinor=1,dtset%nspinor
559            ibsp=ibsp+1
560            do jlmn=1,lmn_size
561              jl=Pawtab(itypat)%indlmn(1,jlmn)
562              jm=Pawtab(itypat)%indlmn(2,jlmn)
563              cp_fact(1) = Cprj_k(iatom,ibsp)%cp(1,jlmn) *sqrt(crystal%ucvol) ! Magic factor
564              cp_fact(2) = Cprj_k(iatom,ibsp)%cp(2,jlmn) *sqrt(crystal%ucvol)
565 
566              do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere.
567                ifftsph = local_pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid
568                phj  = Paw_onsite(iatom)% phi(ifgd,jlmn)
569                tphj = Paw_onsite(iatom)%tphi(ifgd,jlmn)
570                ! old code
571                !re_p = cp_fact(1); im_p = cp_fact(2)
572                ! apply the phase
573                re_p = cp_fact(1) * phk_atm(1,ifgd,iatom) - cp_fact(2) * phk_atm(2,ifgd,iatom)
574                im_p = cp_fact(1) * phk_atm(2,ifgd,iatom) + cp_fact(2) * phk_atm(1,ifgd,iatom)
575 
576                ur_ae(1,ifftsph+(ispinor-1)*nfftot) = ur_ae(1,ifftsph+(ispinor-1)*nfftot) + re_p * (phj-tphj)
577                ur_ae(2,ifftsph+(ispinor-1)*nfftot) = ur_ae(2,ifftsph+(ispinor-1)*nfftot) + im_p * (phj-tphj)
578                ur_ae_onsite(1,ifftsph) = ur_ae_onsite(1,ifftsph) + re_p * phj
579                ur_ae_onsite(2,ifftsph) = ur_ae_onsite(2,ifftsph) + im_p * phj
580                ur_ps_onsite(1,ifftsph) = ur_ps_onsite(1,ifftsph) + re_p * tphj
581                ur_ps_onsite(2,ifftsph) = ur_ps_onsite(2,ifftsph) + im_p * tphj
582                ur_mask(ifftsph) = one
583              end do
584 
585            end do !jlmn
586          end do !ispinor
587        end do !iatom
588 
589        if (paral_atom) then
590          ABI_MALLOC(buf_tmp,(2,n1*n2*n3,3))
591          buf_tmp(:,:,1) = ur_ae
592          buf_tmp(:,:,2) = ur_ae_onsite
593          buf_tmp(:,:,3) = ur_ps_onsite
594          call xmpi_sum(buf_tmp,my_comm_atom,ierr)
595          ur_ae = buf_tmp(:,:,1)
596          ur_ae_onsite= buf_tmp(:,:,2)
597          ur_ps_onsite= buf_tmp(:,:,3)
598          ABI_FREE(buf_tmp)
599        end if
600 
601 !      * Remove the phase e^{ikr}, we store u(r).
602        do i3=1,n3
603          do i2=1,n2
604            do i1=1,n1
605              ii = i1 + n1*(i2-1)+ n1*n2*(i3-1)
606              reur=ur_ae(1,ii)
607              imur=ur_ae(2,ii)
608              ur_ae(1,ii)=  reur * phkr(1,ii) + imur * phkr(2,ii)
609              ur_ae(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
610              if(Dtset%nspinor==2) then
611                reur=ur_ae(1,ii+nfftot)    ! Important!
612                imur=ur_ae(2,ii+nfftot)
613                ur_ae(1,ii+nfftot)=  reur * phkr(1,ii) + imur * phkr(2,ii)
614                ur_ae(2,ii+nfftot)= -reur * phkr(2,ii) + imur * phkr(1,ii)
615              end if
616              reur=ur_ae_onsite(1,ii)
617              imur=ur_ae_onsite(2,ii)
618              ur_ae_onsite(1,ii)=  reur * phkr(1,ii) + imur * phkr(2,ii)
619              ur_ae_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
620              reur=ur_ps_onsite(1,ii)
621              imur=ur_ps_onsite(2,ii)
622              ur_ps_onsite(1,ii)=  reur * phkr(1,ii) + imur * phkr(2,ii)
623              ur_ps_onsite(2,ii)= -reur * phkr(2,ii) + imur * phkr(1,ii)
624            end do
625          end do
626        end do
627 
628        norm=zero
629        do ii=1,npw_k*Dtset%nspinor
630          norm=norm+fofgin(1,ii)**2+fofgin(2,ii)**2
631        end do
632        write(std_out,'(a,2i5,f22.16)',advance='no') 'ikpt,iband, norm (G,PSWF)=',ikpt,iband,norm
633        norm=zero
634        do ifft=1,nfftot*Dtset%nspinor
635          norm = norm + ur_ae(1,ifft)**2+ur_ae(2,ifft)**2
636        end do
637        norm=norm/nfftot
638        norm_rerr = MAX((ABS(norm-one))*100,norm_rerr)
639        write(std_out,*)"norm (R,AEWF)= ",norm
640        call flush_unit(std_out)
641 
642 !      MS: Various testing and debugging options
643        if (.TRUE..and.nprocs==1) then
644          if (present(pseudo_norms)) then
645 !          Check the supposedly zero overlap |\tilde{Psi_n}-\tilde{Psi_n^1}|^2
646            ABI_MALLOC(dummy_1d,(n1*n2*n3))
647            dummy_1d = zero
648            norm = zero
649            do ifft = 1, nfftot
650              dummy_1d(ifft) = ((ur_pw(1,ifft)-ur_ps_onsite(1,ifft))**2 &
651 &             +  (ur_pw(2,ifft)-ur_ps_onsite(2,ifft))**2) * ur_mask(ifft)
652              norm = norm + dummy_1d(ifft)
653            end do
654            norm = norm / nfftot
655            pseudo_norms(isppol,ikpt,iband) = norm
656            ABI_FREE(dummy_1d)
657          end if
658 
659        else
660          write(msg,'(5a)')&
661 &         "The option to print PAW all-electron wavefunctions is on, but execution ",ch10,&
662 &         "is in parallel on two or more processors. XcrysDen files with individual con-",ch10,&
663 &         "tributions will not be written. In order to enable this you must run in serial."
664          MSG_WARNING(msg)
665        end if ! Check if serial run
666 
667 #ifdef HAVE_NETCDF
668        ncerr = nf90_put_var(ncid, ae_ncid, ur_ae, &
669 &       start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
670        NCF_CHECK(ncerr)
671 
672        ncerr = nf90_put_var(ncid, pw_ncid, ur_pw, &
673 &       start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
674        NCF_CHECK(ncerr)
675 
676        ncerr = nf90_put_var(ncid, aeons_ncid, ur_ae_onsite, &
677 &       start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
678        NCF_CHECK(ncerr)
679 
680        ncerr = nf90_put_var(ncid, psons_ncid, ur_ps_onsite, &
681 &       start=[1,1,1,1,1,iband,ikpt,isppol], count=[2,n1,n2,n3,1,1,1,1])
682        NCF_CHECK(ncerr)
683 #endif
684 
685        ABI_FREE(ur_ae)
686        ABI_FREE(ur_ae_onsite)
687        ABI_FREE(ur_ps_onsite)
688        ABI_FREE(ur_pw)
689        ABI_FREE(ur_mask)
690 
691        ABI_FREE(fofgin)
692        ABI_FREE(fofgout)
693        ABI_FREE(denpot)
694      end do !nband_k
695 
696      bdtot_index=bdtot_index+nband_k
697 
698      if (mkmem/=0) then
699        ibg=ibg+dtset%nspinor*nband_k
700        icg=icg+npw_k*dtset%nspinor*nband_k
701        ikg=ikg+npw_k
702      end if
703 
704      ABI_FREE(kg_k)
705 
706      call pawcprj_free(Cprj_k)
707      ABI_DT_FREE(Cprj_k)
708 
709    end do !ikpt
710  end do !nsppol
711 
712  ABI_FREE(phkr)
713  ABI_FREE(gbound)
714 
715  ! Free augmentation waves.
716  call paw_pwaves_lmn_free(Paw_onsite)
717  ABI_DT_FREE(Paw_onsite)
718 
719  ! Maximum relative error over CPUs.
720  call xmpi_max(norm_rerr,max_rerr,comm_cell,ierr)
721  write(std_out,*)"max_rerr=",max_rerr
722 
723  if (max_rerr > ten) then
724    write(msg,'(7a)')&
725 &   "Inaccuracy on the normalization of the wave funtions exceeds 10%. ",ch10,&
726 &   "Likely due to the use of a too coarse FFT mesh or unconverged wavefunctions. ",ch10,&
727 &   "Numerical values inside the augmentation regions might be inaccurate. ",ch10,&
728 &   "Action: increase pawecutdg in your input file. "
729    MSG_COMMENT(msg)
730  end if
731 
732  ABI_FREE(r0shift)
733  ABI_FREE(phk_atm)
734  call pawfgrtab_free(local_pawfgrtab)
735 
736 !Destroy atom table used for parallelism
737  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
738 
739  DBG_EXIT("COLL")
740 
741 end subroutine pawmkaewf