TABLE OF CONTENTS
ABINIT/pawmknhat [ Functions ]
NAME
pawmknhat
FUNCTION
PAW only: Compute compensation charge density (and derivatives) on the fine FFT grid Can also compute first-order compensation charge density (RF calculations)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (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
cplex: if 1, real space 1-order functions on FFT grid are REAL, if 2, COMPLEX distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism gprimd(3,3)=dimensional primitive translations for reciprocal space ider= 0: nhat(r) is computed 1: cartesian derivatives of nhat(r) are computed 2: nhat(r) and derivatives are computed idir=direction of atomic displacement (in case of atomic displ. perturb.) ipert=index of perturbation; must be 0 for ground-state calculations 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 nhatgrdim= -PAW only- 0 if pawgrnhat 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 pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data (1st-order occupancies if ipert>0) pawrhoij0(my_natom) <type(pawrhoij_type)>= GS paw rhoij occupancies and related data (used only if ipert>0) set equat to pawrhoij for GS calculations pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data qphon(3)=wavevector of the phonon (RF only) rprimd(3,3)=dimensional primitive translations for real space ucvol=volume of the unit cell xred(3,natom)= reduced atomic coordinates
OUTPUT
=== if ider=0 or 2
compch_fft=compensation charge inside spheres computed over fine fft grid
pawnhat(nfft,ispden)=nhat on fine rectangular grid
=== if ider=1 or 2
pawgrnhat(nfft,ispden,3)=derivatives of nhat on fine rectangular grid (and derivatives)
PARENTS
bethe_salpeter,dfpt_scfcv,energy,nres2vres,odamix,paw_qpscgw,pawmkrho respfn,scfcv,screening,setup_positron,sigma
CHILDREN
destroy_distribfft,fourdp,free_my_atmtab,get_my_atmtab init_distribfft_seq,initmpi_seq,mean_fftr,pawexpiqr,pawgylm,pawnhatfr set_mpi_enreg_fft,timab,unset_mpi_enreg_fft,xmpi_sum,zerosym
SOURCE
69 #if defined HAVE_CONFIG_H 70 #include "config.h" 71 #endif 72 73 #include "abi_common.h" 74 75 subroutine pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,& 76 & my_natom,natom,nfft,ngfft,nhatgrdim,nspden,ntypat,pawang,pawfgrtab,& 77 & pawgrnhat,pawnhat,pawrhoij,pawrhoij0,pawtab,qphon,rprimd,ucvol,usewvl,xred,& 78 & mpi_atmtab,comm_atom,comm_fft,mpi_comm_wvl,me_g0,paral_kgb,distribfft) ! optional arguments 79 80 81 use defs_basis 82 use m_profiling_abi 83 use m_errors 84 use m_xmpi 85 use m_cgtools 86 87 use defs_abitypes, only : mpi_type 88 use m_pawang, only : pawang_type 89 use m_pawtab, only : pawtab_type 90 use m_pawfgrtab, only : pawfgrtab_type 91 use m_pawrhoij, only : pawrhoij_type 92 use m_paw_finegrid, only : pawgylm, pawexpiqr 93 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 94 use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft 95 use m_fft, only : zerosym 96 use m_distribfft, only : distribfft_type, init_distribfft_seq, destroy_distribfft 97 98 !This section has been created automatically by the script Abilint (TD). 99 !Do not modify the following lines by hand. 100 #undef ABI_FUNC 101 #define ABI_FUNC 'pawmknhat' 102 use interfaces_18_timing 103 use interfaces_51_manage_mpi 104 use interfaces_53_ffts 105 use interfaces_65_paw, except_this_one => pawmknhat 106 !End of the abilint section 107 108 implicit none 109 110 !Arguments --------------------------------------------- 111 !scalars 112 integer,intent(in) :: cplex,ider,idir,ipert,izero,my_natom,natom,nfft 113 integer,intent(in) :: usewvl 114 integer,intent(in) :: nhatgrdim,nspden,ntypat 115 integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,mpi_comm_wvl,paral_kgb 116 real(dp),intent(in) :: ucvol 117 real(dp),intent(out) :: compch_fft 118 type(distribfft_type),optional,intent(in),target :: distribfft 119 type(pawang_type),intent(in) :: pawang 120 !arrays 121 integer,intent(in) :: ngfft(18) 122 integer,optional,target,intent(in) :: mpi_atmtab(:) 123 real(dp),intent(in) :: gprimd(3,3),qphon(3),rprimd(3,3),xred(3,natom) 124 real(dp),intent(out) :: pawgrnhat(cplex*nfft,nspden,3*nhatgrdim) 125 real(dp),intent(inout) :: pawnhat(cplex*nfft,nspden) !vz_i 126 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 127 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom),pawrhoij0(my_natom) 128 type(pawtab_type),intent(in) :: pawtab(ntypat) 129 130 !Local variables --------------------------------------- 131 !scalars 132 integer :: dplex,iatom,iatom_tot,ic,ierr,ii,ils,ilslm,irhoij,ispden,itypat 133 integer :: jc,jrhoij,kc,klm,klmn,lmax,lmin,lm_size,mfgd,mm,mpi_comm_sphgrid 134 integer :: my_comm_atom,my_comm_fft,nfgd,nfftot,option,optgr0,optgr1,optgr2,paral_kgb_fft 135 logical :: compute_grad,compute_nhat,my_atmtab_allocated,need_frozen,paral_atom,qeq0 136 logical :: compute_phonons,has_phase 137 type(distribfft_type),pointer :: my_distribfft 138 type(mpi_type) :: mpi_enreg_fft 139 !arrays 140 integer,pointer :: my_atmtab(:) 141 real(dp) :: ro(cplex),ro_ql(cplex),tmp_compch_fft(nspden),tsec(2) 142 real(dp),allocatable :: pawgrnhat_atm(:,:),pawnhat_atm(:),work(:,:) 143 144 ! ************************************************************************* 145 146 DBG_ENTER("COLL") 147 148 compute_nhat=(ider==0.or.ider==2) 149 compute_grad=(ider==1.or.ider==2) 150 compute_phonons=(ipert>0.and.ipert<=natom) 151 152 !Compatibility tests 153 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-15) 154 if (present(comm_fft)) then 155 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 156 MSG_BUG('Need paral_kgb and me_g0 with comm_fft !') 157 end if 158 end if 159 if(ider>0.and.nhatgrdim==0) then 160 MSG_BUG(' Gradients of nhat required but not allocated !') 161 end if 162 if (my_natom>0) then 163 if(nspden>1.and.nspden/=pawrhoij(1)%nspden) then 164 MSG_BUG(' Wrong values for nspden and pawrhoij%nspden !') 165 end if 166 if(nspden>1.and.nspden/=pawfgrtab(1)%nspden) then 167 MSG_BUG(' Wrong values for nspden and pawfgrtab%nspden !') 168 end if 169 if(pawrhoij(1)%cplex<cplex) then 170 MSG_BUG(' Must have pawrhoij()%cplex >= cplex !') 171 end if 172 if (compute_phonons.and.(.not.qeq0)) then 173 if (pawfgrtab(1)%rfgd_allocated==0) then 174 MSG_BUG(' pawfgrtab()%rfgd array must be allocated !') 175 end if 176 if (compute_grad.and.(.not.compute_nhat)) then 177 MSG_BUG(' When q<>0, nhat gradients need nhat !') 178 end if 179 end if 180 end if 181 182 !nhat1 does not have to be computed for ddk or d2dk 183 if (ipert==natom+1.or.ipert==natom+10) return 184 185 !Set up parallelism over atoms 186 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 187 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 188 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 189 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,& 190 & my_natom_ref=my_natom) 191 192 !Initialisations 193 if ((.not.compute_nhat).and.(.not.compute_grad)) return 194 mfgd=zero;if (my_natom>0) mfgd=maxval(pawfgrtab(1:my_natom)%nfgd) 195 dplex=cplex-1 196 if (compute_nhat) then 197 ABI_ALLOCATE(pawnhat_atm,(cplex*mfgd)) 198 pawnhat=zero 199 end if 200 if (compute_grad) then 201 ABI_ALLOCATE(pawgrnhat_atm,(cplex*mfgd,3)) 202 pawgrnhat=zero 203 end if 204 205 !mpi communicators for spherical grid: 206 mpi_comm_sphgrid=xmpi_comm_self !no communicators passed 207 if(present(comm_fft) .and. usewvl==0) mpi_comm_sphgrid=comm_fft 208 if(present(mpi_comm_wvl) .and. usewvl==1) mpi_comm_sphgrid=mpi_comm_wvl 209 210 !------------------------------------------------------------------------ 211 !----- Loop over atoms 212 !------------------------------------------------------------------------ 213 214 do iatom=1,my_natom 215 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 216 217 itypat=pawrhoij(iatom)%itypat 218 lm_size=pawfgrtab(iatom)%l_size**2 219 need_frozen=((compute_nhat).and.(ipert==iatom_tot.or.ipert==natom+3.or.ipert==natom+4)) 220 nfgd=pawfgrtab(iatom)%nfgd 221 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 222 if (((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)).or.& 223 & ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)).or.& 224 & ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0))) then 225 optgr0=0;optgr1=0;optgr2=0 226 if ((compute_nhat).and.(pawfgrtab(iatom)%gylm_allocated==0)) then 227 if (allocated(pawfgrtab(iatom)%gylm)) then 228 ABI_DEALLOCATE(pawfgrtab(iatom)%gylm) 229 end if 230 ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,pawfgrtab(iatom)%l_size**2)) 231 pawfgrtab(iatom)%gylm_allocated=2;optgr0=1 232 end if 233 if ((compute_grad).and.(pawfgrtab(iatom)%gylmgr_allocated==0)) then 234 if (allocated(pawfgrtab(iatom)%gylmgr)) then 235 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr) 236 end if 237 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,pawfgrtab(iatom)%l_size**2)) 238 pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1 239 end if 240 if ((compute_grad.and.need_frozen).and.(pawfgrtab(iatom)%gylmgr2_allocated==0)) then 241 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 242 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2) 243 end if 244 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,pawfgrtab(iatom)%l_size**2)) 245 pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1 246 end if 247 if (optgr0+optgr1+optgr2>0) then 248 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,& 249 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd) 250 end if 251 end if 252 253 254 ! Eventually compute exp(-i.q.r) factors for the current atom (if not already done) 255 if (compute_phonons.and.(.not.qeq0).and.pawfgrtab(iatom)%expiqr_allocated==0) then 256 if (allocated(pawfgrtab(iatom)%expiqr)) then 257 ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr) 258 end if 259 ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(2,nfgd)) 260 call pawexpiqr(pawfgrtab(iatom)%expiqr,gprimd,nfgd,qphon,& 261 & pawfgrtab(iatom)%rfgd,xred(:,iatom_tot)) 262 pawfgrtab(iatom)%expiqr_allocated=2 263 end if 264 has_phase=(compute_phonons.and.pawfgrtab(iatom)%expiqr_allocated/=0) 265 266 ! Eventually compute frozen part of nhat for the current atom (if not already done) 267 if ((need_frozen).and.((pawfgrtab(iatom)%nhatfr_allocated==0).or.& 268 & (compute_grad.and.pawfgrtab(iatom)%nhatfrgr_allocated==0))) then 269 if (allocated(pawfgrtab(iatom)%nhatfr)) then 270 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr) 271 end if 272 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(nfgd,pawfgrtab(iatom)%nspden)) 273 option=0;pawfgrtab(iatom)%nhatfr_allocated=2 274 if (compute_grad) then 275 option=1 276 if (allocated(pawfgrtab(iatom)%nhatfrgr)) then 277 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr) 278 end if 279 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(3,nfgd,pawfgrtab(iatom)%nspden)) 280 pawfgrtab(iatom)%nhatfrgr_allocated=2 281 end if 282 call pawnhatfr(option,idir,ipert,1,natom,nspden,ntypat,pawang,pawfgrtab(iatom),& 283 & pawrhoij0(iatom),pawtab,rprimd) 284 end if 285 286 ! ------------------------------------------------------------------------ 287 ! ----- Loop over density components 288 ! ------------------------------------------------------------------------ 289 290 do ispden=1,nspden 291 292 if (compute_nhat) pawnhat_atm(1:cplex*nfgd)=zero 293 if (compute_grad) pawgrnhat_atm(1:cplex*nfgd,1:3)=zero 294 295 ! ------------------------------------------------------------------------ 296 ! ----- Loop over ij channels (basis components) 297 ! ------------------------------------------------------------------------ 298 jrhoij=1 299 do irhoij=1,pawrhoij(iatom)%nrhoijsel 300 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 301 klm =pawtab(itypat)%indklmn(1,klmn) 302 lmin=pawtab(itypat)%indklmn(3,klmn) 303 lmax=pawtab(itypat)%indklmn(4,klmn) 304 305 ! Retrieve rhoij 306 if (nspden/=2) then 307 ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden) 308 else 309 if (ispden==1) then 310 ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)& 311 & +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2) 312 else if (ispden==2) then 313 ro(1:cplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1) 314 end if 315 end if 316 ro(1:cplex)=pawtab(itypat)%dltij(klmn)*ro(1:cplex) 317 if (compute_nhat) then 318 if (cplex==1) then 319 do ils=lmin,lmax,2 320 do mm=-ils,ils 321 ilslm=ils*ils+ils+mm+1 322 if (pawang%gntselect(ilslm,klm)>0) then 323 ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn) 324 do ic=1,nfgd 325 pawnhat_atm(ic)=pawnhat_atm(ic)+ro_ql(1)*pawfgrtab(iatom)%gylm(ic,ilslm) 326 end do 327 end if 328 end do 329 end do 330 else 331 do ils=lmin,lmax,2 332 do mm=-ils,ils 333 ilslm=ils*ils+ils+mm+1 334 if (pawang%gntselect(ilslm,klm)>0) then 335 ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn) 336 do ic=1,nfgd 337 jc=2*ic-1 338 pawnhat_atm(jc:jc+1)=pawnhat_atm(jc:jc+1)+ro_ql(1:2)*pawfgrtab(iatom)%gylm(ic,ilslm) 339 end do 340 end if 341 end do 342 end do 343 end if 344 end if 345 346 if (compute_grad) then 347 if (cplex==1) then 348 do ils=lmin,lmax,2 349 do mm=-ils,ils 350 ilslm=ils*ils+ils+mm+1 351 if (pawang%gntselect(ilslm,klm)>0) then 352 ro_ql(1)=ro(1)*pawtab(itypat)%qijl(ilslm,klmn) 353 do ic=1,nfgd 354 pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 355 pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 356 pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+ro_ql(1)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 357 end do 358 end if 359 end do 360 end do 361 else 362 do ils=lmin,lmax,2 363 do mm=-ils,ils 364 ilslm=ils*ils+ils+mm+1 365 if (pawang%gntselect(ilslm,klm)>0) then 366 ro_ql(1:2)=ro(1:2)*pawtab(itypat)%qijl(ilslm,klmn) 367 do ic=1,nfgd 368 jc=2*ic-1 369 pawgrnhat_atm(jc:jc+1,1)=pawgrnhat_atm(jc:jc+1,1) & 370 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(1,ic,ilslm) 371 pawgrnhat_atm(jc:jc+1,2)=pawgrnhat_atm(jc:jc+1,2) & 372 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(2,ic,ilslm) 373 pawgrnhat_atm(jc:jc+1,3)=pawgrnhat_atm(jc:jc+1,3) & 374 & +ro_ql(1:2)*pawfgrtab(iatom)%gylmgr(3,ic,ilslm) 375 end do 376 end if 377 end do 378 end do 379 end if 380 end if 381 382 ! ------------------------------------------------------------------------ 383 ! ----- End loop over ij channels 384 ! ------------------------------------------------------------------------ 385 jrhoij=jrhoij+pawrhoij(iatom)%cplex 386 end do 387 388 ! If RF calculation, add frozen part of 1st-order compensation density 389 if (need_frozen) then 390 if (cplex==1) then 391 do ic=1,nfgd 392 pawnhat_atm(ic)=pawnhat_atm(ic)+pawfgrtab(iatom)%nhatfr(ic,ispden) 393 end do 394 else 395 do ic=1,nfgd 396 jc=2*ic-1 397 pawnhat_atm(jc)=pawnhat_atm(jc)+pawfgrtab(iatom)%nhatfr(ic,ispden) 398 end do 399 end if 400 if (compute_grad) then 401 if (cplex==1) then 402 do ic=1,nfgd 403 pawgrnhat_atm(ic,1)=pawgrnhat_atm(ic,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden) 404 pawgrnhat_atm(ic,2)=pawgrnhat_atm(ic,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden) 405 pawgrnhat_atm(ic,3)=pawgrnhat_atm(ic,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden) 406 end do 407 else 408 do ic=1,nfgd 409 jc=2*ic-1 410 pawgrnhat_atm(jc,1)=pawgrnhat_atm(jc,1)+pawfgrtab(iatom)%nhatfrgr(1,ic,ispden) 411 pawgrnhat_atm(jc,2)=pawgrnhat_atm(jc,2)+pawfgrtab(iatom)%nhatfrgr(2,ic,ispden) 412 pawgrnhat_atm(jc,3)=pawgrnhat_atm(jc,3)+pawfgrtab(iatom)%nhatfrgr(3,ic,ispden) 413 end do 414 end if 415 end if 416 end if 417 418 ! If needed, multiply eventually by exp(-i.q.r) phase 419 if (has_phase) then 420 if (cplex==1) then 421 do ic=1,nfgd 422 pawnhat_atm(ic)=pawnhat_atm(ic)*pawfgrtab(iatom)%expiqr(1,ic) 423 end do 424 else 425 do ic=1,nfgd 426 jc=2*ic-1 427 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 428 ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic) 429 ro(1:2)=pawnhat_atm(jc:jc+1) 430 pawnhat_atm(jc )=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 431 pawnhat_atm(jc+1)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 432 end do 433 end if 434 if (compute_grad) then 435 if (cplex==1) then 436 do ic=1,nfgd 437 pawgrnhat_atm(ic,1:3)=pawgrnhat_atm(ic,1:3)*pawfgrtab(iatom)%expiqr(1,ic) 438 end do 439 else 440 do ic=1,nfgd 441 jc=2*ic-1 442 ! dn^hat(r)/dr_i * exp(-i.q.r) 443 ro_ql(1)= pawfgrtab(iatom)%expiqr(1,ic) 444 ro_ql(2)=-pawfgrtab(iatom)%expiqr(2,ic) 445 do ii=1,3 446 ro(1:2)=pawgrnhat_atm(jc:jc+1,ii) 447 pawgrnhat_atm(jc ,ii)=ro(1)*ro_ql(1)-ro(2)*ro_ql(2) 448 pawgrnhat_atm(jc+1,ii)=ro(2)*ro_ql(1)+ro(1)*ro_ql(2) 449 end do 450 ! -i.q_i * [n^hat(r).exp(-i.q.r)] 451 ro(1:2)=pawnhat_atm(jc:jc+1) 452 do ii=1,3 453 pawgrnhat_atm(jc ,ii)=pawgrnhat_atm(kc ,ii)+qphon(ii)*ro(2) 454 pawgrnhat_atm(jc+1,ii)=pawgrnhat_atm(kc+1,ii)-qphon(ii)*ro(1) 455 end do 456 end do 457 end if 458 end if 459 end if 460 461 ! Add the contribution of the atom to the compensation charge 462 if (compute_nhat) then 463 if (cplex==1) then 464 do ic=1,nfgd 465 kc=pawfgrtab(iatom)%ifftsph(ic) 466 pawnhat(kc,ispden)=pawnhat(kc,ispden)+pawnhat_atm(ic) 467 end do 468 else 469 do ic=1,nfgd 470 jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1 471 pawnhat(kc:kc+1,ispden)=pawnhat(kc:kc+1,ispden)+pawnhat_atm(jc:jc+1) 472 end do 473 end if 474 end if 475 if (compute_grad) then 476 if (cplex==1) then 477 do ic=1,nfgd 478 kc=pawfgrtab(iatom)%ifftsph(ic) 479 pawgrnhat(kc,ispden,1:3)=pawgrnhat(kc,ispden,1:3)+pawgrnhat_atm(ic,1:3) 480 end do 481 else 482 do ic=1,nfgd 483 jc=2*ic-1;kc=2*pawfgrtab(iatom)%ifftsph(ic)-1 484 do ii=1,3 485 pawgrnhat(kc:kc+1,ispden,ii)=pawgrnhat(kc:kc+1,ispden,ii)+pawgrnhat_atm(jc:jc+1,ii) 486 end do 487 end do 488 end if 489 end if 490 ! ------------------------------------------------------------------------ 491 ! ----- End loop over density components 492 ! ------------------------------------------------------------------------ 493 end do 494 495 if (pawfgrtab(iatom)%gylm_allocated==2) then 496 ABI_DEALLOCATE(pawfgrtab(iatom)%gylm) 497 ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(0,0)) 498 pawfgrtab(iatom)%gylm_allocated=0 499 end if 500 if (pawfgrtab(iatom)%gylmgr_allocated==2) then 501 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr) 502 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0)) 503 pawfgrtab(iatom)%gylmgr_allocated=0 504 end if 505 if (pawfgrtab(iatom)%gylmgr2_allocated==2) then 506 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2) 507 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(0,0,0)) 508 pawfgrtab(iatom)%gylmgr2_allocated=0 509 end if 510 if (pawfgrtab(iatom)%nhatfr_allocated==2) then 511 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr) 512 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(0,0)) 513 pawfgrtab(iatom)%nhatfr_allocated=0 514 end if 515 if (pawfgrtab(iatom)%nhatfrgr_allocated==2) then 516 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr) 517 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(0,0,0)) 518 pawfgrtab(iatom)%nhatfrgr_allocated=0 519 end if 520 if (pawfgrtab(iatom)%expiqr_allocated==2) then 521 ABI_DEALLOCATE(pawfgrtab(iatom)%expiqr) 522 ABI_ALLOCATE(pawfgrtab(iatom)%expiqr,(0,0)) 523 pawfgrtab(iatom)%expiqr_allocated=0 524 end if 525 526 ! ------------------------------------------------------------------------ 527 ! ----- End loop over atoms 528 ! ------------------------------------------------------------------------ 529 end do 530 531 !----- Free some memory 532 if (compute_nhat) then 533 ABI_DEALLOCATE(pawnhat_atm) 534 end if 535 if (compute_grad) then 536 ABI_DEALLOCATE(pawgrnhat_atm) 537 end if 538 539 !----- Reduction in case of parallelism 540 if (paral_atom) then 541 call timab(48,1,tsec) 542 if (compute_nhat) then 543 call xmpi_sum(pawnhat,my_comm_atom,ierr) 544 end if 545 if (compute_grad) then 546 call xmpi_sum(pawgrnhat,my_comm_atom,ierr) 547 end if 548 call timab(48,2,tsec) 549 end if 550 551 !----- Avoid unbalanced g-components numerical errors 552 if (izero==1.and.compute_nhat.and.usewvl==0) then 553 ! Create fake mpi_enreg to wrap fourdp 554 if (present(distribfft)) then 555 my_distribfft => distribfft 556 else 557 ABI_DATATYPE_ALLOCATE(my_distribfft,) 558 call init_distribfft_seq(my_distribfft,'f',ngfft(2),ngfft(3),'fourdp') 559 end if 560 call initmpi_seq(mpi_enreg_fft) 561 ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft) 562 if (present(comm_fft)) then 563 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 564 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 565 else 566 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 567 mpi_enreg_fft%distribfft => my_distribfft 568 end if 569 ! do FFT 570 ABI_ALLOCATE(work,(2,nfft)) 571 do ispden=1,min(2,nspden) 572 call fourdp(cplex,work,pawnhat(:,ispden),-1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 573 call zerosym(work,2,ngfft(1),ngfft(2),ngfft(3),comm_fft=my_comm_fft,distribfft=my_distribfft) 574 call fourdp(cplex,work,pawnhat(:,ispden),+1,mpi_enreg_fft,nfft,ngfft,paral_kgb_fft,0) 575 end do 576 ABI_DEALLOCATE(work) 577 ! Destroy fake mpi_enreg 578 call unset_mpi_enreg_fft(mpi_enreg_fft) 579 if (.not.present(distribfft)) then 580 call destroy_distribfft(my_distribfft) 581 ABI_DATATYPE_DEALLOCATE(my_distribfft) 582 end if 583 end if 584 585 !----- Computation of compensation charge over real space grid 586 if (compute_nhat.and.ipert==0) then 587 nfftot=PRODUCT(ngfft(1:3)) 588 call mean_fftr(pawnhat,tmp_compch_fft,nfft,nfftot,1,mpi_comm_sphgrid) 589 compch_fft = tmp_compch_fft(1) 590 compch_fft=compch_fft*ucvol 591 end if 592 593 !Destroy atom table used for parallelism 594 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 595 596 DBG_EXIT("COLL") 597 598 end subroutine pawmknhat