TABLE OF CONTENTS
ABINIT/pawmknhat_psipsi [ 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