TABLE OF CONTENTS
ABINIT/paw_mknewh0 [ Functions ]
NAME
paw_mknewh0
FUNCTION
Calculates the new bare PAW Hamiltonian in the case of quasi-particle self-consistent GW calculations.
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
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 nsppol=1 for unpolarized, 2 for spin-polarized nspden=number of spin-density components nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid pawspnorb=flag: 1 if spin-orbit coupling is activated pawprtvol=control print volume and debugging output for PAW Cryst<crystal_t>=Info on unit cell and its symmetries Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data Paw_an(natom) <type(paw_an_type)>=paw arrays given on angular mesh Pawang<type(pawang_type)>=paw angular mesh and related data Pawfgrtab(natom) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid vxc(nfftf,nspden)=exchange-correlation potential vxc_val(nfftf,nspden)=valence only exchange-correlation potential vtrial(nfftf,nspden)=potential (Hartree+XC+loc)
SIDE EFFECTS
Paw_ij(natom*usepaw)<Paw_ij_type)>=paw arrays given on (i,j) channels At output: new value for Paw_ij()%dij
PARENTS
calc_vhxc_me
CHILDREN
free_my_atmtab,get_my_atmtab,pawgylm,symdij,symdij_all,wrtout
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 subroutine paw_mknewh0(my_natom,nsppol,nspden,nfftf,pawspnorb,pawprtvol,Cryst,& 51 & Pawtab,Paw_an,Paw_ij,Pawang,Pawfgrtab,vxc,vxc_val,vtrial,& 52 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 53 54 use defs_basis 55 use m_profiling_abi 56 use m_errors 57 use m_xmpi, only : xmpi_comm_self 58 59 use m_crystal, only : crystal_t 60 use m_pawang, only : pawang_type 61 use m_pawtab, only : pawtab_type 62 use m_paw_an, only : paw_an_type 63 use m_paw_ij, only : paw_ij_type 64 use m_pawfgrtab, only : pawfgrtab_type 65 use m_pawdij, only : symdij, symdij_all 66 use m_paw_finegrid, only : pawgylm 67 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'paw_mknewh0' 73 use interfaces_14_hidewrite 74 !End of the abilint section 75 76 implicit none 77 78 !Arguments ------------------------------------ 79 !scalars 80 integer,intent(in) :: my_natom,nsppol,nspden,nfftf,pawprtvol,pawspnorb 81 integer,optional,intent(in) :: comm_atom 82 !arrays 83 integer,optional,target,intent(in) :: mpi_atmtab(:) 84 real(dp),intent(in) :: vxc(nfftf,nspden),vxc_val(nfftf,nspden),vtrial(nfftf,nspden) 85 type(crystal_t),intent(in) :: Cryst 86 type(Pawang_type),intent(in) :: Pawang 87 type(Pawtab_type),target,intent(in) :: Pawtab(Cryst%ntypat) 88 type(Paw_an_type),intent(in) :: Paw_an(my_natom) 89 type(Paw_ij_type),intent(inout) :: Paw_ij(my_natom) 90 type(Pawfgrtab_type),intent(inout) :: Pawfgrtab(my_natom) 91 92 !Local variables------------------------------- 93 !scalars 94 integer,parameter :: ipert0=0 95 integer :: iat,iat_tot,idij,cplex,ndij,option_dij 96 integer :: itypat,lmn_size,j0lmn,jlmn,ilmn,klmn,klmn1,klm 97 integer :: lmin,lmax,mm,isel,lm_size,lmn2_size,my_comm_atom,cplex_dij 98 integer :: ils,ilslm,ic,lm0 99 integer :: nsploop,is2fft 100 real(dp) :: gylm,qijl 101 logical :: ltest,my_atmtab_allocated,paral_atom 102 character(len=500) :: msg 103 !arrays 104 integer, pointer :: indklmn_(:,:) 105 integer,pointer :: my_atmtab(:) 106 real(dp) :: rdum(1),rdum2(1) 107 real(dp),allocatable :: prod_hloc(:,:),prodhxc_core(:,:) 108 real(dp),allocatable :: dijhl_hat(:,:),dijhmxc_val(:,:) 109 110 ! ************************************************************************* 111 112 DBG_ENTER("COLL") 113 114 call wrtout(std_out,'Assembling PAW strengths for the bare Hamiltonian','COLL') 115 116 !== Set up parallelism over atoms === 117 paral_atom=(present(comm_atom).and.(my_natom/=Cryst%natom)) 118 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 119 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 120 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,Cryst%natom,my_natom_ref=my_natom) 121 122 if (my_natom>0) then 123 124 ! === Test if required pointers in paw_ij are allocated === 125 ltest = (allocated(Paw_ij(1)%dijxc).and.allocated(Paw_ij(1)%dijxc_val) ) 126 ABI_CHECK(ltest,'dijxc or dijxc_val not calculated') 127 128 ltest=(allocated(Paw_ij(1)%dijhat)) !.and.Paw_ij(1)%has_dijhat==2) 129 ABI_CHECK(ltest,'dijhat not calculated') 130 131 ltest=(allocated(Paw_ij(1)%dijhartree)) !.and.Paw_ij(1)%has_dijhartree==2) 132 ABI_CHECK(ltest,'dijhartree not calculated') 133 134 if (ANY(Pawtab(:)%usepawu>0)) then 135 do iat=1,my_natom 136 iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat) 137 itypat=Cryst%typat(iat_tot) 138 if (Pawtab(itypat)%usepawu>0) then 139 ltest=(allocated(Paw_ij(iat)%dijU) ) !.and.Paw_ij(iat)%has_dijU==2) 140 write(msg,'(a,i3,a)')" For atom no. ",iat," %dijU(iat) has not been calculated." 141 ABI_CHECK(ltest,msg) 142 end if 143 end do 144 end if 145 146 if (pawspnorb>0) then 147 do iat=1,my_natom 148 ltest=(allocated(Paw_ij(iat)%dijso) ) !.and.Paw_ij(iat)%has_dijso==2) 149 write(msg,'(a,i3,a)')" For atom no. ",iat," %dijso(iat) has not been calculated." 150 ABI_CHECK(ltest,msg) 151 end do 152 end if 153 end if ! my_natom>0 154 155 !== Construct the new PAW H0 Hamiltonian === 156 do iat=1,my_natom 157 iat_tot=iat;if (paral_atom) iat_tot=my_atmtab(iat) 158 159 itypat = Cryst%typat(iat_tot) 160 lmn_size = Pawtab(itypat)%lmn_size 161 lmn2_size = Pawtab(itypat)%lmn2_size 162 lm_size = Paw_an(iat)%lm_size 163 cplex = Paw_ij(iat)%cplex 164 cplex_dij = Paw_ij(iat)%cplex_dij 165 ndij = Paw_ij(iat)%ndij 166 167 ABI_CHECK(cplex==1,'cplex/=1 not implemented') 168 ABI_CHECK(cplex_dij==1,'cplex_dij/=1 not implemented') 169 ! 170 ! Eventually compute g_l(r).Y_lm(r) factors for the current atom (if not already done) 171 if (Pawfgrtab(iat)%gylm_allocated==0) then 172 if (allocated(Pawfgrtab(iat)%gylm)) then 173 ABI_DEALLOCATE(Pawfgrtab(iat)%gylm) 174 end if 175 ABI_ALLOCATE(Pawfgrtab(iat)%gylm,(Pawfgrtab(iat)%nfgd,lm_size)) 176 Pawfgrtab(iat)%gylm_allocated=2 177 178 call pawgylm(Pawfgrtab(iat)%gylm,rdum,rdum2,lm_size,& 179 & Pawfgrtab(iat)%nfgd,1,0,0,Pawtab(itypat),Pawfgrtab(iat)%rfgd) 180 end if 181 182 ! === Calculate LM contribution to dijhmxc_val for this atom === 183 ! * Dijxc contains also the Hat term on the FFT mesh while Dijxc_val does not 184 ! contain neither the hat term nor the LM sum of onsite terms (they should cancel each other) 185 ! FIXME change paw_dij, otherwise I miss tnc in vxc 186 ! * prodhxc_core is used to assemble $\int g_l Ylm (vtrial - vxc_val[tn+nhat] dr$ on the FFT mesh === 187 ! * The following quantities do not depend on ij 188 ABI_ALLOCATE(prod_hloc ,(lm_size,ndij)) 189 ABI_ALLOCATE(prodhxc_core,(lm_size,ndij)) 190 prod_hloc =zero 191 prodhxc_core=zero 192 do idij=1,ndij 193 do ilslm=1,lm_size 194 do ic=1,Pawfgrtab(iat)%nfgd 195 is2fft=Pawfgrtab(iat)%ifftsph(ic) 196 gylm=Pawfgrtab(iat)%gylm(ic,ilslm) 197 prod_hloc (ilslm,idij)=prod_hloc (ilslm,idij) + (vtrial(is2fft,idij)-vxc(is2fft,idij))*gylm 198 ! prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vxc_val(is2fft,idij))*gylm 199 prodhxc_core(ilslm,idij)=prodhxc_core(ilslm,idij) + (vtrial(is2fft,idij)-vxc_val(is2fft,idij))*gylm 200 end do 201 end do 202 end do !idij 203 204 ! === Assembly the "Hat" contribution for this atom ==== 205 ABI_ALLOCATE(dijhl_hat ,(cplex_dij*lmn2_size,ndij)) 206 ABI_ALLOCATE(dijhmxc_val,(cplex_dij*lmn2_size,ndij)) 207 dijhl_hat =zero 208 dijhmxc_val=zero 209 indklmn_ => Pawtab(itypat)%indklmn(1:6,1:lmn2_size) 210 211 do idij=1,ndij 212 do klmn=1,lmn2_size 213 klm =indklmn_(1,klmn) 214 lmin=indklmn_(3,klmn) 215 lmax=indklmn_(4,klmn) 216 217 ! === $\sum_lm q_ij^l prod* for each idij$ === 218 do ils=lmin,lmax,2 219 lm0=ils**2+ils+1 220 do mm=-ils,ils 221 ilslm=lm0+mm 222 isel=Pawang%gntselect(lm0+mm,klm) 223 if (isel>0) then 224 qijl=Pawtab(itypat)%qijl(ilslm,klmn) 225 dijhl_hat (klmn,idij)=dijhl_hat (klmn,idij) + prod_hloc (ilslm,idij)*qijl 226 dijhmxc_val(klmn,idij)=dijhmxc_val(klmn,idij) +prodhxc_core(ilslm,idij)*qijl 227 end if 228 end do 229 end do 230 end do 231 end do 232 233 ABI_DEALLOCATE(prod_hloc) 234 ABI_DEALLOCATE(prodhxc_core) 235 236 ! * Normalization factor due to integration on the FFT mesh 237 dijhl_hat = dijhl_hat *Cryst%ucvol/DBLE(nfftf) 238 dijhmxc_val= dijhmxc_val*Cryst%ucvol/DBLE(nfftf) 239 240 ! === Now assembly the bare Hamiltonian === 241 ! * Loop over density components overwriting %dij 242 nsploop=nsppol; if (Paw_ij(iat)%ndij==4) nsploop=4 243 244 do idij=1,nsploop 245 klmn1=1 246 247 do jlmn=1,lmn_size 248 j0lmn=jlmn*(jlmn-1)/2 249 do ilmn=1,jlmn 250 klmn=j0lmn+ilmn 251 252 ! The following gives back the input dij. 253 ! since dijxc contains the hat term done on the FFT mesh 254 if (.FALSE.) then 255 Paw_ij(iat)%dij(klmn,idij) = & 256 & Pawtab(itypat)%dij0 (klmn) & 257 & +Paw_ij(iat)%dijhartree(klmn) & 258 & +Paw_ij(iat)%dijxc (klmn,idij) & 259 & +dijhl_hat (klmn,idij) 260 261 else 262 ! === Make nonlocal part of h0 removing the valence contribution === 263 ! Remeber that XC contains already the Hat contribution 264 Paw_ij(iat)%dij(klmn,idij) = & 265 & Pawtab(itypat)%dij0 (klmn) & 266 & +Paw_ij(iat)%dijhartree(klmn) & 267 & +Paw_ij(iat)%dijxc (klmn,idij) & ! 2 lines to get the d1-dt1 XC core contribution + XC hat (core+val) 268 & -Paw_ij(iat)%dijxc_val (klmn,idij) & ! I suppose that the "hat" term on the FFT mesh in included in both. 269 & +dijhmxc_val(klmn,idij) ! Local + Hartree - XC val contribution to the "hat" term. 270 271 ! Add the U contribution to the 272 ! if (.FALSE. .and. Pawtab(itypat)%usepawu>0) then 273 if (.TRUE. .and. Pawtab(itypat)%usepawu>0) then 274 Paw_ij(iat)%dij(klmn,idij) = Paw_ij(iat)%dij(klmn,idij) + Paw_ij(iat)%dijU(klmn,idij) 275 end if 276 end if 277 ! TODO dijso, dijU, vpawx? 278 ! Just to be consistent, update some values. 279 !$Paw_ij(iat)%dijhat(klmn,idij)=Paw_ij(iat)%dijhat(klmn,idij)-dijhmxc_val(klmn,idij) 280 281 end do !ilmn 282 end do !jlmn 283 end do !idij 284 285 ! this is to be consistent? 286 ! deallocate(Paw_ij(iat)%dijvxc_val) 287 ABI_DEALLOCATE(dijhl_hat) 288 ABI_DEALLOCATE(dijhmxc_val) 289 end do !iat 290 291 !=== Symmetrize total Dij === 292 option_dij=0 ! For total Dij. 293 #if 0 294 if (paral_atom) then 295 call symdij(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 296 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,& 297 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 298 else 299 call symdij(Cryst%gprimd,,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,option_dij,& 300 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 301 end if 302 #else 303 if (paral_atom) then 304 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 305 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec,& 306 & comm_atom=my_comm_atom,mpi_atmtab=my_atmtab) 307 else 308 call symdij_all(Cryst%gprimd,Cryst%indsym,ipert0,my_natom,Cryst%natom,Cryst%nsym,Cryst%ntypat,& 309 & Paw_ij,Pawang,pawprtvol,Pawtab,Cryst%rprimd,Cryst%symafm,Cryst%symrec) 310 end if 311 #endif 312 313 !Destroy atom table used for parallelism 314 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 315 316 DBG_EXIT("COLL") 317 318 end subroutine paw_mknewh0