TABLE OF CONTENTS
ABINIT/pawnhatfr [ Functions ]
NAME
pawnhatfr
FUNCTION
PAW: Compute frozen part of charge compensation density nhat nhatfr(r)=Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)] Depends on q wave vector but not on first-order wave-function.
COPYRIGHT
Copyright (C) 2011-2018 ABINIT group (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.
INPUTS
ider=0: computes frozen part of compensation density 1: computes frozen part of compensation density and cartesian gradients idir=direction of atomic displacement (in case of phonons perturb.) ipert=nindex of perturbation mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor natom=total number of atoms in cell nspden=number of spin-density components ntypat=number of types of atoms 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)>= Ground-State paw rhoij occupancies and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data rprimd(3,3)=dimensional primitive translations for real space
OUTPUT
pawfgrtab(iatom)%nhatfr(nfgd,nspden) frozen part of charge compensation density (inside PAW spheres) =Sum_ij,lm[rhoij_ij.q_ij^l.(g_l(r).Y_lm(r))^(1)] === If ider==1 pawfgrtab(iatom)%nhatfrgr(3,nfgd,nspden) gradients of frozen part of charge compensation density (inside PAW spheres) =Sum_ij,lm[rhoij_ij.q_ij^l . d/dr((g_l(r).Y_lm(r))^(1))]
PARENTS
dfpt_nstpaw,dfpt_scfcv,pawmknhat
CHILDREN
free_my_atmtab,get_my_atmtab,pawgylm
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine pawnhatfr(ider,idir,ipert,my_natom,natom,nspden,ntypat,& 60 & pawang,pawfgrtab,pawrhoij,pawtab,rprimd, & 61 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 62 63 use defs_basis 64 use m_profiling_abi 65 use m_errors 66 use m_xmpi, only : xmpi_comm_self 67 68 use m_pawang, only : pawang_type 69 use m_pawtab, only : pawtab_type 70 use m_pawfgrtab, only : pawfgrtab_type 71 use m_pawrhoij, only : pawrhoij_type 72 use m_paw_finegrid, only : pawgylm 73 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 74 75 !This section has been created automatically by the script Abilint (TD). 76 !Do not modify the following lines by hand. 77 #undef ABI_FUNC 78 #define ABI_FUNC 'pawnhatfr' 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: ider,idir,ipert,my_natom,natom,nspden,ntypat 86 integer,optional,intent(in) :: comm_atom 87 type(pawang_type),intent(in) :: pawang 88 !arrays 89 integer,optional,target,intent(in) :: mpi_atmtab(:) 90 real(dp),intent(in) :: rprimd(3,3) 91 type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom) 92 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 93 type(pawtab_type),intent(in) :: pawtab(ntypat) 94 95 !Local variables------------------------------- 96 !scalars 97 integer :: iatom,iatom_tot,ic,ils,ilslm,irhoij,isel,ispden,istr,itypat,jrhoij 98 integer :: klm,klmn,lm_size,lmn2_size,lm0,lmax,lmin,mua,mub,mm,mu,my_comm_atom,nfgd,nu,optgr0,optgr1,optgr2 99 logical :: my_atmtab_allocated,my_pert,paral_atom 100 real(dp) :: contrib,ro 101 !arrays 102 integer,parameter :: voigt(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/)) 103 integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/) 104 integer,pointer :: my_atmtab(:) 105 real(dp),allocatable :: nhatfr_tmp(:,:),nhatfrgr_tmp(:,:,:) 106 107 ! ************************************************************************* 108 109 DBG_ENTER("COLL") 110 111 !Only relevant for atomic displacement and strain perturbation 112 if (ipert>natom.and.ipert/=natom+3.and.ipert/=natom+4) return 113 114 !Set up parallelism over atoms 115 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 116 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 117 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 118 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 119 120 !Compatibility tests 121 if (my_natom>0) then 122 if ((pawfgrtab(1)%gylm_allocated==0.or.pawfgrtab(1)%gylmgr_allocated==0).and. & 123 & pawfgrtab(1)%rfgd_allocated==0) then 124 MSG_BUG(' pawfgrtab()%rfgd array must be allocated !') 125 end if 126 end if 127 128 my_pert = (ipert<=natom).or.ipert==natom+3.or.ipert==natom+4 129 130 !Get correct index of strain pertubation 131 if (ipert==natom+3) istr = idir 132 if (ipert==natom+4) istr = idir + 3 133 134 !Loops over atoms 135 do iatom=1,my_natom 136 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 137 138 ! Eventually allocate frozen nhat points 139 if (my_pert) then 140 if (pawfgrtab(iatom)%nhatfr_allocated==0) then 141 if (allocated(pawfgrtab(iatom)%nhatfr)) then 142 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfr) 143 end if 144 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfr,(pawfgrtab(iatom)%nfgd,nspden)) 145 pawfgrtab(iatom)%nhatfr_allocated=1 146 end if 147 if (ider==1.and.pawfgrtab(iatom)%nhatfrgr_allocated==0) then 148 if (allocated(pawfgrtab(iatom)%nhatfrgr)) then 149 ABI_DEALLOCATE(pawfgrtab(iatom)%nhatfrgr) 150 end if 151 ABI_ALLOCATE(pawfgrtab(iatom)%nhatfrgr,(3,pawfgrtab(iatom)%nfgd,nspden)) 152 pawfgrtab(iatom)%nhatfrgr_allocated=1 153 end if 154 end if 155 156 ! Select if frozen part of nhat exists for the current perturbation 157 if ((.not.my_pert).or.(pawfgrtab(iatom)%nhatfr_allocated==0)) cycle 158 159 ! Some atom-dependent quantities 160 itypat=pawfgrtab(iatom)%itypat 161 lm_size=pawfgrtab(iatom)%l_size**2 162 lmn2_size=pawtab(itypat)%lmn2_size 163 164 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 165 nfgd=pawfgrtab(iatom)%nfgd 166 if ((pawfgrtab(iatom)%gylmgr_allocated==0).or. & 167 & (pawfgrtab(iatom)%gylmgr2_allocated==0.and.ider==1)) then 168 optgr0=0;optgr1=0;optgr2=0 169 if(ipert==natom+3.or.ipert==natom+4)then 170 if (pawfgrtab(iatom)%gylm_allocated==0) then 171 if (allocated(pawfgrtab(iatom)%gylm)) then 172 ABI_DEALLOCATE(pawfgrtab(iatom)%gylm) 173 end if 174 ABI_ALLOCATE(pawfgrtab(iatom)%gylm,(nfgd,lm_size)) 175 pawfgrtab(iatom)%gylm_allocated=2;optgr0=1 176 end if 177 end if 178 if (pawfgrtab(iatom)%gylmgr_allocated==0) then 179 if (allocated(pawfgrtab(iatom)%gylmgr)) then 180 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr) 181 end if 182 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(3,nfgd,lm_size)) 183 pawfgrtab(iatom)%gylmgr_allocated=2;optgr1=1 184 end if 185 if (ider==1.and.pawfgrtab(iatom)%gylmgr2_allocated==0) then 186 if (allocated(pawfgrtab(iatom)%gylmgr2)) then 187 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2) 188 end if 189 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(6,nfgd,lm_size)) 190 pawfgrtab(iatom)%gylmgr2_allocated=2;optgr2=1 191 end if 192 if (optgr0+optgr1+optgr2>0) then 193 call pawgylm(pawfgrtab(iatom)%gylm,pawfgrtab(iatom)%gylmgr,pawfgrtab(iatom)%gylmgr2,& 194 & lm_size,nfgd,optgr0,optgr1,optgr2,pawtab(itypat),pawfgrtab(iatom)%rfgd) 195 end if 196 end if 197 198 199 ! ============ Phonons ==================================== 200 if (ipert<=natom) then 201 202 ! Loop over spin components 203 do ispden=1,nspden 204 205 ABI_ALLOCATE(nhatfr_tmp,(3,nfgd)) 206 nhatfr_tmp=zero 207 if (ider==1) then 208 ABI_ALLOCATE(nhatfrgr_tmp,(3,nfgd,3)) 209 nhatfrgr_tmp=zero 210 end if 211 212 jrhoij=1 213 do irhoij=1,pawrhoij(iatom)%nrhoijsel 214 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 215 klm =pawtab(itypat)%indklmn(1,klmn) 216 lmin=pawtab(itypat)%indklmn(3,klmn) 217 lmax=pawtab(itypat)%indklmn(4,klmn) 218 if (nspden/=2) then 219 ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 220 else 221 if (ispden==1) then 222 ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2) 223 else if (ispden==2) then 224 ro=pawrhoij(iatom)%rhoijp(jrhoij,1) 225 end if 226 end if 227 ro=pawtab(itypat)%dltij(klmn)*ro 228 do ils=lmin,lmax,2 229 lm0=ils**2+ils+1 230 do mm=-ils,ils 231 ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm) 232 if (isel>0) then 233 do ic=1,nfgd 234 do mu=1,3 235 contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn)& 236 & *pawfgrtab(iatom)%gylmgr(mu,ic,ilslm) 237 nhatfr_tmp(mu,ic)=nhatfr_tmp(mu,ic)+contrib 238 end do 239 end do 240 if (ider==1) then 241 do ic=1,nfgd 242 do nu=1,3 243 do mu=1,3 244 contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) & 245 & *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm) 246 nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib 247 end do 248 end do 249 end do 250 end if 251 end if 252 end do 253 end do 254 jrhoij=jrhoij+pawrhoij(iatom)%cplex 255 end do 256 257 ! Convert from cartesian to reduced coordinates 258 do ic=1,nfgd 259 pawfgrtab(iatom)%nhatfr(ic,ispden)= & 260 & rprimd(1,idir)*nhatfr_tmp(1,ic) & 261 & +rprimd(2,idir)*nhatfr_tmp(2,ic) & 262 & +rprimd(3,idir)*nhatfr_tmp(3,ic) 263 end do 264 if (ider==1) then 265 do nu=1,3 266 do ic=1,nfgd 267 pawfgrtab(iatom)%nhatfrgr(nu,ic,ispden)= & 268 & rprimd(1,idir)*nhatfrgr_tmp(1,ic,nu) & 269 & +rprimd(2,idir)*nhatfrgr_tmp(2,ic,nu) & 270 & +rprimd(3,idir)*nhatfrgr_tmp(3,ic,nu) 271 end do 272 end do 273 end if 274 ABI_DEALLOCATE(nhatfr_tmp) 275 if (ider==1) then 276 ABI_DEALLOCATE(nhatfrgr_tmp) 277 end if 278 ! End loop over spin components 279 end do ! ispden 280 281 282 ! ============ Elastic tensor =============================== 283 else if (ipert==natom+3.or.ipert==natom+4) then 284 ! Loop over spin components 285 pawfgrtab(iatom)%nhatfr(:,:) = zero 286 do ispden=1,nspden 287 jrhoij=1 288 do irhoij=1,pawrhoij(iatom)%nrhoijsel 289 klmn=pawrhoij(iatom)%rhoijselect(irhoij) 290 klm =pawtab(itypat)%indklmn(1,klmn) 291 lmin=pawtab(itypat)%indklmn(3,klmn) 292 lmax=pawtab(itypat)%indklmn(4,klmn) 293 if (nspden/=2) then 294 ro=pawrhoij(iatom)%rhoijp(jrhoij,ispden) 295 else 296 if (ispden==1) then 297 ro=pawrhoij(iatom)%rhoijp(jrhoij,1)+pawrhoij(iatom)%rhoijp(jrhoij,2) 298 else if (ispden==2) then 299 ro=pawrhoij(iatom)%rhoijp(jrhoij,1) 300 end if 301 end if 302 ro=pawtab(itypat)%dltij(klmn)*ro 303 do ils=lmin,lmax,2 304 lm0=ils**2+ils+1 305 do mm=-ils,ils 306 ilslm=lm0+mm;isel=pawang%gntselect(lm0+mm,klm) 307 if (isel>0) then 308 ! Sum{[Q_ij_q^LM^(1)]} 309 do ic=1,nfgd 310 mua=alpha(istr);mub=beta(istr) 311 pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+& 312 & ro*pawtab(itypat)%qijl(ilslm,klmn)*half*(& 313 & pawfgrtab(iatom)%gylmgr(mua,ic,ilslm)*pawfgrtab(iatom)%rfgd(mub,ic)& 314 & +pawfgrtab(iatom)%gylmgr(mub,ic,ilslm)*pawfgrtab(iatom)%rfgd(mua,ic)) 315 end do 316 ! Add volume contribution 317 if(istr<=3)then 318 do ic=1,nfgd 319 pawfgrtab(iatom)%nhatfr(ic,ispden) = pawfgrtab(iatom)%nhatfr(ic,ispden)+& 320 & ro*pawtab(itypat)%qijl(ilslm,klmn)*pawfgrtab(iatom)%gylm(ic,ilslm) 321 end do 322 end if 323 if (ider==1) then 324 MSG_ERROR("nhatgr not implemented for strain perturbationxs") 325 ! do ic=1,nfgd 326 ! do nu=1,6 327 ! do mu=1,6 328 ! contrib=-ro*pawtab(itypat)%qijl(ilslm,klmn) & 329 !& *pawfgrtab(iatom)%gylmgr2(voigt(mu,nu),ic,ilslm) 330 ! nhatfrgr_tmp(mu,ic,nu)=nhatfrgr_tmp(mu,ic,nu)+contrib 331 ! end do 332 ! end do 333 ! end do 334 end if 335 end if 336 end do 337 end do 338 jrhoij=jrhoij+pawrhoij(iatom)%cplex 339 end do 340 end do ! ispden 341 end if 342 343 ! Eventually free temporary space for g_l(r).Y_lm(r) gradients and exp(-i.q.r) 344 if (pawfgrtab(iatom)%gylmgr_allocated==2) then 345 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr) 346 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr,(0,0,0)) 347 pawfgrtab(iatom)%gylmgr_allocated=0 348 end if 349 if (pawfgrtab(iatom)%gylmgr2_allocated==2) then 350 ABI_DEALLOCATE(pawfgrtab(iatom)%gylmgr2) 351 ABI_ALLOCATE(pawfgrtab(iatom)%gylmgr2,(0,0,0)) 352 pawfgrtab(iatom)%gylmgr2_allocated=0 353 end if 354 355 ! End loop on atoms 356 end do 357 358 359 !Destroy atom table used for parallelism 360 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 361 362 DBG_EXIT("COLL") 363 364 end subroutine pawnhatfr