TABLE OF CONTENTS


ABINIT/pawmknhat_psipsi [ Functions ]

[ Top ] [ Functions ]

NAME

 pawmknhat_psipsi

FUNCTION

 PAW only:
 Compute on the fine FFT grid the compensation charge density (and derivatives) associated
 to the product of two wavefunctions n_{12}(r) = \Psi_1* \Psi_2. Based on pawmknhat.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (MG, FJ, 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  cprj1(natom,nspinor), cprj2(natom,nspinor) <type(pawcprj_type)>=
   projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
   the \Psi_1 and \Psi_2, respectively.
  distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
  ider= 0: nhat(r) is computed
        1: cartesian derivatives of nhat(r) are computed
        2: nhat(r) and derivatives are computed
        3: nhat(r) and gradients of nhat  wrt atomic coordinates are computed
        Note: ider>0 not compatible with ipert>0
  izero=if 1, unbalanced components of nhat(g) have to be set to zero
  me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
  mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
  comm_atom=--optional-- MPI communicator over atoms
  comm_fft=--optional-- MPI communicator over FFT components
  my_natom=number of atoms treated by current processor
  natom=total number of atoms in cell
  nfft=number of point on the rectangular fft grid
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nhat12_grdim= 0 if grnhat12 array is not used ; 1 otherwise
  ntypat=number of types of atoms in unit cell.
  paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawfgrtab(my_natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  === if ider=0 or 2
    nhat12(2,nfft,nspinor**2)=nhat on fine rectangular grid*exp(iqr)
  === if ider=1 or 2
    grnhat12(nfft,nspinor**2,3)=gradient of (nhat*exp(iqr)) on fine rectangular grid (derivative versus r)
  === if ider=3
    grnhat_12(nfft,nspinor**2,3,natom*(ider/3))=derivatives of nhat on fine rectangular grid versus R*exp(iqr)

PARENTS

      calc_sigx_me,fock_getghc,prep_calc_ucrpa

CHILDREN

      destroy_distribfft,fourdp,free_my_atmtab,get_my_atmtab
      init_distribfft_seq,initmpi_seq,pawexpiqr,pawgylm,set_mpi_enreg_fft
      timab,unset_mpi_enreg_fft,xmpi_sum,zerosym

SOURCE

 62 #if defined HAVE_CONFIG_H
 63 #include "config.h"
 64 #endif
 65 
 66 #include "abi_common.h"
 67 
 68 subroutine pawmknhat_psipsi(cprj1,cprj2,ider,izero,my_natom,natom,nfft,ngfft,nhat12_grdim,&
 69 &          nspinor,ntypat,pawang,pawfgrtab,grnhat12,nhat12,pawtab, &
 70 &          gprimd,grnhat_12,qphon,xred,atindx,mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments
 71 
 72  use defs_basis
 73  use m_profiling_abi
 74  use m_errors
 75  use m_xmpi
 76 
 77  use defs_abitypes,    only : mpi_type
 78  use m_mpinfo,         only : set_mpi_enreg_fft,unset_mpi_enreg_fft
 79  use m_lmn_indices,    only : klmn2ijlmn
 80  use m_pawang,         only : pawang_type
 81  use m_pawtab,         only : pawtab_type
 82  use m_pawfgrtab,      only : pawfgrtab_type
 83  use m_pawcprj,        only : pawcprj_type
 84  use m_paw_finegrid,   only : pawgylm,pawexpiqr
 85  use m_paral_atom,     only : get_my_atmtab, free_my_atmtab
 86  use m_fft,            only : zerosym
 87  use m_distribfft,     only : distribfft_type, init_distribfft_seq, destroy_distribfft
 88 
 89 ! use m_lmn_indices,  only : klmn2ijlmn
 90 
 91 !This section has been created automatically by the script Abilint (TD).
 92 !Do not modify the following lines by hand.
 93 #undef ABI_FUNC
 94 #define ABI_FUNC 'pawmknhat_psipsi'
 95  use interfaces_18_timing
 96  use interfaces_51_manage_mpi
 97  use interfaces_53_ffts
 98 !End of the abilint section
 99 
100  implicit none
101 
102 !Arguments ---------------------------------------------
103 !scalars
104  integer,intent(in) :: ider,izero,my_natom,natom,nfft,nhat12_grdim,ntypat,nspinor
105  integer,optional,intent(in) :: me_g0,comm_fft,paral_kgb
106  integer,optional,intent(in) :: comm_atom
107  type(distribfft_type),optional,intent(in),target :: distribfft
108  type(pawang_type),intent(in) :: pawang
109 !arrays
110  integer,intent(in) :: ngfft(18)
111  integer,optional,intent(in) ::atindx(natom)
112  integer,optional,target,intent(in) :: mpi_atmtab(:)
113  real(dp),optional, intent(in) ::gprimd(3,3),qphon(3),xred(3,natom)
114  real(dp),intent(out) :: grnhat12(2,nfft,nspinor**2,3*nhat12_grdim)
115  real(dp),optional,intent(out) :: grnhat_12(2,nfft,nspinor**2,3,natom*(ider/3))
116  real(dp),intent(out) :: nhat12(2,nfft,nspinor**2)
117  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom)
118  type(pawtab_type),intent(in) :: pawtab(ntypat)
119  type(pawcprj_type),intent(in) :: cprj1(natom,nspinor),cprj2(natom,nspinor)
120 
121 !Local variables ---------------------------------------
122 !scalars
123  integer :: iatm,iatom,iatom_tot,ic,ierr,ils,ilslm,isp1,isp2,isploop,itypat,jc,klm,klmn
124  integer :: lmax,lmin,lm_size,mm,my_comm_atom,my_comm_fft,optgr0,optgr1,paral_kgb_fft
125  integer :: cplex,ilmn,jlmn,lmn_size,lmn2_size
126  real(dp) :: re_p,im_p
127  logical :: compute_grad,compute_grad1,compute_nhat,my_atmtab_allocated,paral_atom,qeq0,compute_phonon,order
128  type(distribfft_type),pointer :: my_distribfft
129  type(mpi_type) :: mpi_enreg_fft
130 !arrays
131  integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
132  integer,pointer :: my_atmtab(:)
133  real(dp) :: rdum(1),cpf(2),cpf_ql(2),tsec(2),ro(2),ro_ql(2),nhat12_atm(2,nfft,nspinor**2)
134  real(dp),allocatable :: work(:,:), qijl(:,:)
135 
136 ! *************************************************************************
137 
138  DBG_ENTER("COLL")
139 
140 !Compatibility tests
141  if (present(comm_fft)) then
142    if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then
143      MSG_BUG('Need paral_kgb and me_g0 with comm_fft !')
144    end if
145    if (present(paral_kgb)) then
146      if (paral_kgb/=0) then
147        MSG_BUG('paral_kgb/=0 not coded!')
148      end if
149    end if
150  end if
151  if (ider>0.and.nhat12_grdim==0) then
152 !   MSG_BUG('Gradients of nhat required but not allocated !')
153  end if
154  if (nspinor==2) then
155    MSG_BUG('nspinor==2 not coded!')
156  end if
157 
158  compute_phonon=.false.;qeq0=.false.
159  if (present(gprimd).and.present(qphon).and.present(xred)) compute_phonon=.true.
160  if (compute_phonon) qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15)
161  if (present(atindx)) order=.true.
162 !Set up parallelism over atoms
163  paral_atom=(present(comm_atom).and.(my_natom/=natom))
164  nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
165  my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
166  call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
167 
168 !Initialisations
169  compute_nhat=(ider==0.or.ider==2.or.ider==3)
170  compute_grad=(ider==1.or.ider==2)
171  compute_grad1=(ider==3)
172  if ((.not.compute_nhat).and.(.not.compute_grad)) return
173 
174  if (compute_nhat) nhat12=zero
175  if (compute_grad) grnhat12=zero
176  if (compute_grad1) grnhat_12=zero
177 
178  if (compute_grad) then
179 !   MSG_BUG('compute_grad not tested!')
180  end if
181 
182 !------------------------------------------------------------------------
183 !----- Loop over atoms
184 !------------------------------------------------------------------------
185  do iatom=1,my_natom
186    iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
187    iatm=iatom_tot
188    if (order) iatm=atindx(iatom_tot)
189    itypat    = pawfgrtab(iatom)%itypat
190    lm_size   = pawfgrtab(iatom)%l_size**2
191    lmn_size  = pawtab(itypat)%lmn_size
192    lmn2_size = pawtab(itypat)%lmn2_size
193    ABI_ALLOCATE(qijl,(lm_size,lmn2_size))
194    qijl=zero
195    qijl=pawtab(itypat)%qijl
196    if (compute_nhat) nhat12_atm=zero
197 !  Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done)
198    if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.&
199 &   (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0))) then
200      optgr0=0; optgr1=0
201      if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then
202        if (allocated(pawfgrtab(iatom)%gylm))  then
203          ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
204        end if
205        ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
206        pawfgrtab(iatom)%gylm_allocated=2;optgr0=1
207      end if
208 
209      if (((compute_grad).or.(compute_grad1)).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then
210        if (allocated(pawfgrtab(iatom)%gylmgr))  then
211          ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
212        end if
213        ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,pawfgrtab(iatom)%nfgd,pawfgrtab(iatom)%l_size**2))
214        pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1
215      end if
216 
217      if (optgr0+optgr1>0) then
218        call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,rdum,&
219 &       lm_size,pawfgrtab(iatom)%nfgd,optgr0,optgr1,0,pawtab(itypat),&
220 &       pawfgrtab(iatom)%rfgd)
221      end if
222 
223    end if
224    if (compute_phonon.and.(.not.qeq0).and.(pawfgrtab(iatom)%expiqr_allocated==0)) then
225      if (allocated(pawfgrtab(iatom)%expiqr))  then
226        ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
227      end if
228      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,pawfgrtab(iatom)%nfgd))
229      call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,pawfgrtab(iatom)%nfgd,qphon,&
230 &     pawfgrtab(iatom)%rfgd,xred(:,iatom_tot))
231      pawfgrtab(iatom)%expiqr_allocated=2
232    end if
233 
234    do isploop=1,nspinor**2    ! Loop over density components of the compensation charge.
235 !    TODO Here we might take advantage of symmetry relations between the four components if nspinor==2
236      isp1=spinor_idxs(1,isploop)
237      isp2=spinor_idxs(2,isploop)
238 
239      do klmn=1,lmn2_size  ! Loop over ij channels of this atom type.
240        klm =pawtab(itypat)%indklmn(1,klmn)
241        lmin=pawtab(itypat)%indklmn(3,klmn)  ! abs(il-jl)
242        lmax=pawtab(itypat)%indklmn(4,klmn)  ! il+jl
243        ilmn=pawtab(itypat)%indklmn(7,klmn)
244        jlmn=pawtab(itypat)%indklmn(8,klmn)
245 !       call klmn2ijlmn(klmn,lmn_size,ilmn,jlmn)  ! This mapping should be stored in pawtab_type
246 
247 !      Retrieve the factor due to the PAW projections.
248        re_p =  cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) &
249 &       +cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) &
250 &       +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn) &
251 &       +cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn)
252 
253        im_p =  cprj1(iatm,isp1)%cp(1,ilmn) * cprj2(iatm,isp2)%cp(2,jlmn) &
254 &       -cprj1(iatm,isp1)%cp(2,ilmn) * cprj2(iatm,isp2)%cp(1,jlmn) &
255 &       +cprj1(iatm,isp1)%cp(1,jlmn) * cprj2(iatm,isp2)%cp(2,ilmn) &
256 &       -cprj1(iatm,isp1)%cp(2,jlmn) * cprj2(iatm,isp2)%cp(1,ilmn)
257 
258        cpf(1)=re_p*pawtab(itypat)%dltij(klmn)*half
259        cpf(2)=im_p*pawtab(itypat)%dltij(klmn)*half
260 
261        if (compute_nhat) then
262          do ils=lmin,lmax,2   ! Sum over (L,M)
263            do mm=-ils,ils
264              ilslm=ils*ils+ils+mm+1
265              if (pawang%gntselect(ilslm,klm)>0) then
266                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
267                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
268                do ic=1,pawfgrtab(iatom)%nfgd
269                  jc=pawfgrtab(iatom)%ifftsph(ic)
270                  nhat12_atm(1,jc,isploop)=nhat12_atm(1,jc,isploop)+cpf_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm)
271                  nhat12_atm(2,jc,isploop)=nhat12_atm(2,jc,isploop)+cpf_ql(2)*pawfgrtab(iatom)%gylm(ic,ilslm)
272                end do
273              end if
274            end do
275          end do
276        end if ! compute_nhat
277 
278        if (compute_grad) then
279          do ils=lmin,lmax,2  ! Sum over (L,M)
280            do mm=-ils,ils
281              ilslm=ils*ils+ils+mm+1
282              if (pawang%gntselect(ilslm,klm)>0) then
283                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
284                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
285                do ic=1,pawfgrtab(iatom)%nfgd
286                  jc=pawfgrtab(iatom)%ifftsph(ic)
287                  grnhat12(1,jc,isploop,1)=grnhat12(1,jc,isploop,1)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
288                  grnhat12(1,jc,isploop,2)=grnhat12(1,jc,isploop,2)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
289                  grnhat12(1,jc,isploop,3)=grnhat12(1,jc,isploop,3)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
290 
291                  grnhat12(2,jc,isploop,1)=grnhat12(2,jc,isploop,1)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
292                  grnhat12(2,jc,isploop,2)=grnhat12(2,jc,isploop,2)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
293                  grnhat12(2,jc,isploop,3)=grnhat12(2,jc,isploop,3)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
294                end do
295              end if
296            end do
297          end do
298        end if ! compute_grad
299        if (compute_grad1) then
300          do ils=lmin,lmax,2  ! Sum over (L,M)
301            do mm=-ils,ils
302              ilslm=ils*ils+ils+mm+1
303              if (pawang%gntselect(ilslm,klm)>0) then
304                cpf_ql(1)=cpf(1)*qijl(ilslm,klmn)
305                cpf_ql(2)=cpf(2)*qijl(ilslm,klmn)
306                do ic=1,pawfgrtab(iatom)%nfgd
307                  jc=pawfgrtab(iatom)%ifftsph(ic)
308                  grnhat_12(1,jc,isploop,1,iatom)=grnhat_12(1,jc,isploop,1,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
309                  grnhat_12(1,jc,isploop,2,iatom)=grnhat_12(1,jc,isploop,2,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
310                  grnhat_12(1,jc,isploop,3,iatom)=grnhat_12(1,jc,isploop,3,iatom)+cpf_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
311 
312                  grnhat_12(2,jc,isploop,1,iatom)=grnhat_12(2,jc,isploop,1,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm)
313                  grnhat_12(2,jc,isploop,2,iatom)=grnhat_12(2,jc,isploop,2,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm)
314                  grnhat_12(2,jc,isploop,3,iatom)=grnhat_12(2,jc,isploop,3,iatom)+cpf_ql(2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm)
315                end do
316              end if
317            end do
318          end do
319        end if ! compute_grad1
320      end do  ! klmn (ij channels)
321 !    If needed, multiply eventually by exp(-i.q.r) phase
322      if (compute_nhat) then
323        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
324          do ic=1,pawfgrtab(iatom)%nfgd
325            jc=pawfgrtab(iatom)%ifftsph(ic)
326            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
327            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
328            ro(1:2)=nhat12_atm(1:2,jc,isploop)
329            nhat12_atm(1,jc,isploop)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
330            nhat12_atm(2,jc,isploop)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
331          end do
332        end if
333      end if
334 
335      if (compute_grad) then
336        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
337          do ic=1,pawfgrtab(iatom)%nfgd
338            jc=pawfgrtab(iatom)%ifftsph(ic)
339            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
340            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
341            ro(1)=grnhat12(1,jc,isploop,1)-qphon(1)*nhat12_atm(2,jc,isploop)
342            ro(2)=grnhat12(2,jc,isploop,1)+qphon(1)*nhat12_atm(1,jc,isploop)
343            grnhat12(1,jc,isploop,1)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
344            grnhat12(2,jc,isploop,1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
345            ro(1)=grnhat12(1,jc,isploop,2)-qphon(2)*nhat12_atm(2,jc,isploop)
346            ro(2)=grnhat12(2,jc,isploop,2)+qphon(2)*nhat12_atm(1,jc,isploop)
347            grnhat12(1,jc,isploop,2)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
348            grnhat12(2,jc,isploop,2)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
349            ro(1)=grnhat12(1,jc,isploop,3)-qphon(3)*nhat12_atm(2,jc,isploop)
350            ro(2)=grnhat12(2,jc,isploop,3)+qphon(3)*nhat12_atm(1,jc,isploop)
351            grnhat12(1,jc,isploop,3)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
352            grnhat12(2,jc,isploop,3)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
353          end do
354        end if
355      end if
356      if (compute_grad1) then
357        if(compute_phonon.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated/=0) then
358          do ic=1,pawfgrtab(iatom)%nfgd
359            jc=pawfgrtab(iatom)%ifftsph(ic)
360            ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic)
361            ro_ql(2)= pawfgrtab(iatom)%expiqr(2,ic)
362            ro(1)=grnhat_12(1,jc,isploop,1,iatom)
363            ro(2)=grnhat_12(2,jc,isploop,1,iatom)
364            grnhat_12(1,jc,isploop,1,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
365            grnhat_12(2,jc,isploop,1,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
366            ro(1)=grnhat_12(1,jc,isploop,2,iatom)
367            ro(2)=grnhat_12(2,jc,isploop,2,iatom)
368            grnhat_12(1,jc,isploop,2,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
369            grnhat_12(2,jc,isploop,2,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
370            ro(1)=grnhat_12(1,jc,isploop,3,iatom)
371            ro(2)=grnhat_12(2,jc,isploop,3,iatom)
372            grnhat_12(1,jc,isploop,3,iatom)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2)
373            grnhat_12(2,jc,isploop,3,iatom)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2)
374          end do
375        end if
376      end if
377    end do ! isploop (density components of the compensation charge)
378 ! accumlate nhat12 for all the atoms
379    if (compute_nhat) nhat12=nhat12+nhat12_atm
380 
381    if (pawfgrtab(iatom)%gylm_allocated==2) then
382      ABI_DEALLOCATE(pawfgrtab(iatom)%gylm)
383      ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(0,0))
384      pawfgrtab(iatom)%gylm_allocated=0
385    end if
386    if (pawfgrtab(iatom)%gylmgr_allocated==2) then
387      ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr)
388      ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0))
389      pawfgrtab(iatom)%gylmgr_allocated=0
390    end if
391    ABI_DEALLOCATE(qijl)
392    if (pawfgrtab(iatom)%expiqr_allocated==2) then
393      ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr)
394      ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(0,0))
395      pawfgrtab(iatom)%expiqr_allocated=0
396    end if
397 
398  end do ! iatom
399 
400  if (compute_grad1) grnhat_12=-grnhat_12
401 
402 !----- Reduction in case of parallelism -----!
403  if (paral_atom)then
404    call timab(48,1,tsec)
405    if (compute_nhat) then
406      call xmpi_sum(nhat12,my_comm_atom,ierr)
407    end if
408    if (compute_grad) then
409      call xmpi_sum(grnhat12,my_comm_atom,ierr)
410    end if
411    if (compute_grad1) then
412      call xmpi_sum(grnhat_12,my_comm_atom,ierr)
413    end if
414    call timab(48,2,tsec)
415  end if
416 
417 !----- Avoid unbalanced g-components numerical errors -----!
418 
419  if (izero==1.and.compute_nhat) then
420 !  Create fake mpi_enreg to wrap fourdp
421    if (present(distribfft)) then
422      my_distribfft => distribfft
423    else
424      ABI_DATATYPE_ALLOCATE(my_distribfft,)
425      call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp')
426    end if
427    call initmpi_seq(mpi_enreg_fft)
428    ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft)
429    if (present(comm_fft)) then
430      call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb)
431      my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb
432    else
433      my_comm_fft=xmpi_comm_self;paral_kgb_fft=0;
434      mpi_enreg_fft%distribfft => my_distribfft
435    end if
436 !  Do FFT
437    ABI_ALLOCATE(work,(2,nfft))
438    cplex=2
439    do isp1=1,MIN(2,nspinor**2)
440      call fourdp(cplex,work,nhat12(:,:,isp1),-1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
441      call zerosym(work,cplex,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft)
442      call fourdp(cplex,work,nhat12(:,:,isp1),+1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0)
443    end do
444    ABI_DEALLOCATE(work)
445 !  Destroy fake mpi_enreg
446    call unset_mpi_enreg_fft(mpi_enreg_fft)
447    if (.not.present(distribfft)) then
448      call destroy_distribfft(my_distribfft)
449      ABI_DATATYPE_DEALLOCATE(my_distribfft)
450    end if
451  end if
452 
453 !Destroy atom table used for parallelism
454  call free_my_atmtab(my_atmtab,my_atmtab_allocated)
455 
456  DBG_EXIT("COLL")
457 
458 end subroutine pawmknhat_psipsi