TABLE OF CONTENTS
ABINIT/calc_ubare [ Functions ]
NAME
calc_ubare
FUNCTION
Calculate the bare interaction on atomic orbitals
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon) 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
itypatcor = value of itypat for correlated species lpawu = angular momentum for correlated species pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: pawang %lmax=Maximum value of angular momentum l+1 %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data: %mesh_size=Dimension of radial mesh %rad(mesh_size)=The coordinates of all the points of the radial mesh %radfact(mesh_size)=Factor used to compute radial integrals
OUTPUT
NOTES
PARENTS
pawpuxinit
CHILDREN
poisson,simp_gen,wrtout
SOURCE
41 #if defined HAVE_CONFIG_H 42 #include "config.h" 43 #endif 44 45 #include "abi_common.h" 46 47 subroutine calc_ubare(itypatcor,lpawu,pawang,pawrad,pawtab,rmax) 48 49 use defs_basis 50 use m_profiling_abi 51 use m_xmpi 52 use m_errors 53 54 use m_pawang, only : pawang_type, pawang_init, pawang_free 55 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen, poisson, pawrad_ifromr 56 use m_pawtab, only : pawtab_type 57 58 !This section has been created automatically by the script Abilint (TD). 59 !Do not modify the following lines by hand. 60 #undef ABI_FUNC 61 #define ABI_FUNC 'calc_ubare' 62 use interfaces_14_hidewrite 63 !End of the abilint section 64 65 implicit none 66 !Arguments ------------------------------------ 67 integer, intent(in) :: itypatcor,lpawu 68 type(pawang_type),intent(in) :: pawang 69 type(pawrad_type),intent(in) :: pawrad 70 type(pawtab_type),target,intent(in) :: pawtab 71 real(dp), optional, intent(in) :: rmax 72 73 !Local variables ------------------------------ 74 !scalars 75 integer :: ilmn,ilmn1,iln,iln1,isel,isel1,itypat,jlmn,jlmn1,jln,jln1 76 integer :: klm,klm1,klmn,klmn1,ll,lm0 77 integer :: lmin,lmax,lmn2_size,mesh_size,meshsz,mm 78 real(dp) :: norm,r_for_intg,rg,rg1,ubare,uint,uint_tmp 79 character(len=800) :: message 80 !arrays 81 real(dp),allocatable :: ff(:),gg(:),phiphj(:),phiphj1(:) 82 83 !************************************************************************ 84 85 itypat=itypatcor 86 write(message,'(11a,f12.4,2a,i7,2a,f12.4,2a,i7,2a,f12.4)') & 87 & ch10," =======================================================================",ch10, & 88 & " == Calculation of diagonal bare Coulomb interaction on ATOMIC orbitals ",ch10, & 89 & " (it is assumed that the wavefunction for the first reference ",ch10, & 90 & " energy in PAW atomic data is an atomic eigenvalue)",ch10,ch10, & 91 & " Max value of the radius in atomic data file =", pawrad%rmax ,ch10, & 92 & " Max value of the mesh in atomic data file =", pawrad%mesh_size,ch10, & 93 & " PAW radius is =", pawtab%rpaw,ch10, & 94 & " PAW value of the mesh for integration is =", pawrad%int_meshsz,ch10, & 95 & " Integral of atomic wavefunction until rpaw =", pawtab%ph0phiint(1) 96 if(.not.present(rmax)) then 97 call wrtout(ab_out,message,'COLL') 98 call wrtout(std_out,message,'COLL') 99 end if 100 101 mesh_size=pawrad%mesh_size 102 103 ! Definition of the mesh used for integration. 104 if(present(rmax)) then 105 if(rmax>pawrad%rmax) then 106 write(message, '(a)' ) 'calc_ubare: the radius cannot be larger than the maximum radius of the mesh' 107 MSG_ERROR(message) 108 end if 109 meshsz=pawrad_ifromr(pawrad,rmax)+5 110 r_for_intg=rmax 111 else 112 meshsz=pawtab%partialwave_mesh_size 113 r_for_intg=pawrad%rad(meshsz) ! (we could use r_for_intg=-1) 114 end if 115 116 lmn2_size=pawtab%lmn2_size 117 ABI_ALLOCATE(ff,(mesh_size)) 118 ABI_ALLOCATE(gg,(mesh_size)) 119 ABI_ALLOCATE(phiphj,(mesh_size)) 120 ABI_ALLOCATE(phiphj1,(mesh_size)) 121 do klmn=1,lmn2_size 122 ilmn=pawtab%indklmn(7,klmn);jlmn=pawtab%indklmn(8,klmn) 123 ! Select lpawu and first projectors il=jl=lpawu and first proj only 124 if (( pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. & 125 & (-pawtab%indklmn(3,klmn)+pawtab%indklmn(4,klmn)==2*lpawu).and. & 126 & (pawtab%indlmn(3,ilmn)==1).and.(pawtab%indlmn(3,jlmn)==1) ) then 127 klm=pawtab%indklmn(1,klmn);iln=pawtab%indlmn(5,ilmn);jln=pawtab%indlmn(5,jlmn) 128 lmin=pawtab%indklmn(3,klmn);lmax=pawtab%indklmn(4,klmn) 129 phiphj(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln) 130 !write(6,*) "A",klmn,pawtab%klmntomn(1,klmn),pawtab%klmntomn(2,klmn),& 131 !&pawtab%indklmn(7,klmn),pawtab%indklmn(8,klmn),pawtab%klmntomn(3,klmn),pawtab%klmntomn(4,klmn) 132 do ll=lmin,lmin,2 133 lm0=ll*ll+ll+1 134 ff(1:meshsz)=phiphj(1:meshsz) 135 call simp_gen(norm,ff,pawrad,r_for_intg=r_for_intg) 136 call poisson(ff,ll,pawrad,gg) 137 do klmn1=klmn,lmn2_size 138 ilmn1=pawtab%indklmn(7,klmn);jlmn1=pawtab%indklmn(8,klmn) 139 ! Select lpawu and first projectors il=jl=lpawu and first proj only 140 if (( pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. & 141 & (-pawtab%indklmn(3,klmn1)+pawtab%indklmn(4,klmn1)==2*lpawu).and. & 142 & (pawtab%indlmn(3,ilmn1)==1).and.(pawtab%indlmn(3,jlmn1)==1) ) then 143 !write(6,*) "A1",klmn1,pawtab%klmntomn(1,klmn1),pawtab%klmntomn(2,klmn1),& 144 !&pawtab%indklmn(7,klmn1),pawtab%indklmn(8,klmn1),pawtab%klmntomn(3,klmn1),pawtab%klmntomn(4,klmn1) 145 klm1=pawtab%indklmn(1,klmn1);iln1=pawtab%indlmn(5,ilmn1);jln1=pawtab%indlmn(5,jlmn1) 146 phiphj1(1:meshsz)=pawtab%phi(1:meshsz,iln1)*pawtab%phi(1:meshsz,jln1) 147 uint_tmp=zero 148 if ((ll==lmin)) then 149 ff(1)=zero 150 ff(2:meshsz)=phiphj1(2:meshsz)*gg(2:meshsz)*two/pawrad%rad(2:meshsz) 151 call simp_gen(uint_tmp,ff,pawrad,r_for_intg=r_for_intg) 152 end if 153 uint=zero 154 do mm=-ll,ll 155 isel =pawang%gntselect(lm0+mm,klm) 156 isel1=pawang%gntselect(lm0+mm,klm1) 157 if (isel>0.and.isel1>0) then 158 rg =pawang%realgnt(isel) 159 rg1=pawang%realgnt(isel1) 160 uint=uint+uint_tmp*rg*rg1*two_pi 161 end if 162 end do 163 if((pawtab%indklmn(5,klmn)==pawtab%indklmn(6,klmn)).and.& 164 & (pawtab%indklmn(5,klmn1)==pawtab%indklmn(6,klmn1)).and.& 165 & (pawtab%indklmn(5,klmn)==pawtab%indklmn(5,klmn1))) then 166 ubare=uint*Ha_eV 167 end if 168 end if 169 end do 170 end do 171 end if 172 end do 173 ABI_DEALLOCATE(gg) 174 ABI_DEALLOCATE(ff) 175 ABI_DEALLOCATE(phiphj) 176 ABI_DEALLOCATE(phiphj1) 177 178 write(message,'(a,3(a,f12.4,a),2a,f12.4,a)') ch10," For an atomic wfn truncated at rmax =",r_for_intg,ch10,& 179 & " The norm of the wfn is =",norm,ch10,& 180 & " The bare interaction (no renormalization) =",ubare," eV",ch10,& 181 & " The bare interaction (for a renorm. wfn ) =",ubare/norm/norm," eV" 182 call wrtout(ab_out,message,'COLL') 183 call wrtout(std_out,message,'COLL') 184 if(r_for_intg < 10_dp .and. .not.present(rmax)) then 185 write(message,'(a,f6.2,4a)') ' ( WARNING: The radial mesh in the atomic data file is cut at',r_for_intg,ch10,& 186 & ' Use XML atomic data files to compute the bare Coulomb interaction',ch10,& 187 & ' on a true normalized atomic wavefunction )' 188 call wrtout(ab_out,message,'COLL') 189 call wrtout(std_out,message,'COLL') 190 end if 191 if(present(rmax)) then 192 write(message,'(2a)') " =======================================================================",ch10 193 call wrtout(ab_out,message,'COLL') 194 call wrtout(std_out,message,'COLL') 195 end if 196 197 end subroutine calc_ubare