TABLE OF CONTENTS


ABINIT/pawmkaewf [ Functions ]

[ Top ] [ Functions ]

NAME

 pawmkaewf

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 .

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

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