TABLE OF CONTENTS
ABINIT/make_fc_paw [ Functions ]
NAME
make_fc_paw
FUNCTION
Compute the Fermi-contact term due to the PAW cores
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (JZ,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 .
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 natom=number of atoms in cell. nspden=number of spin ensity component ntypat=number of atom types pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
fc(nspden,natom)=the Fermi-contact interaction at each site due to PAW for each spin density
NOTES
The Fermi contact interaction is the electron density evaluated exactly at the nuclear site.
For a nuclear site at R, we are thus computing the expectation value of $\delta^3(R)$, the
the three-dimensional delta function at vector position $R$. In terms of the radial variable only
the delta function is $\delta(r)/4\pi r^2$. Because this observable is
absolutely confined within the PAW radius, only the response due to the AE PAW functions is
needed, the pseudo wavefunctions and pseudo PAW functions cancel each other out. We then
must compute the integral of $u_i/r times u_j/r \delta(R)d^3r$, for the $l=0$ angular momentum
states only. This is simplified with the use of L'H\^{o}spital's theorem to take the limit
as $r\rightarrow 0$, yielding $u_i'(r) u_j'(r)$. To compute the derivatives we just fit the
first 5 points of the $u$ functions to a line through the origin, using the least squares
procedure resulting from $\chi = sum_i (y_i - m*x_i)^2$ . This is more stable than
computing the derivative of the whole function and extrapolating it to zero.
See Zwanziger, J. Phys. Conden. Matt. 21, 15024-15036 (2009).
PARENTS
calc_fc
CHILDREN
free_my_atmtab,get_my_atmtab,xmpi_sum
SOURCE
52 #if defined HAVE_CONFIG_H 53 #include "config.h" 54 #endif 55 56 #include "abi_common.h" 57 58 subroutine make_fc_paw(fc,my_natom,natom,nspden,ntypat,pawrhoij,pawrad,pawtab,& 59 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 60 61 62 use defs_basis 63 use m_profiling_abi 64 use m_errors 65 use m_xmpi, only : xmpi_comm_self 66 67 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 68 use m_pawrad, only : pawrad_type 69 use m_pawtab, only : pawtab_type 70 use m_pawrhoij, only : pawrhoij_type 71 use m_xmpi, only : xmpi_sum 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'make_fc_paw' 77 !End of the abilint section 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 integer,intent(in) :: my_natom,natom,nspden,ntypat 84 integer,optional,intent(in) :: comm_atom 85 !arrays 86 integer,optional,target,intent(in) :: mpi_atmtab(:) 87 real(dp),intent(out) :: fc(nspden,natom) 88 type(pawrad_type),intent(in) :: pawrad(ntypat) 89 type(pawrhoij_type),intent(in) :: pawrhoij(my_natom) 90 type(pawtab_type),target,intent(in) :: pawtab(ntypat) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: iatom,iatom_tot,ierr,irhoij,islope,ispden,itypat 95 integer :: ilmn,il,iln,ilm,im,jl,jlm,jlmn,jln,jm,j0lmn 96 integer :: klmn,kln,mesh_size,my_comm_atom,nslope 97 logical :: my_atmtab_allocated,paral_atom 98 real(dp) :: mi,mj,xi,xxsum,xysumi,xysumj,yi,yj 99 !arrays 100 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 101 integer,pointer :: my_atmtab(:) 102 103 ! ************************************************************************ 104 105 DBG_ENTER("COLL") 106 107 !Set up parallelism over atoms 108 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 109 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 110 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 111 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 112 113 !Initialization 114 fc(:,:)=zero 115 116 !number of points to use in computing initial slopes of radial functions 117 nslope = 5 118 119 !loop over atoms in cell 120 do iatom = 1, my_natom 121 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 122 itypat = pawrhoij(iatom)%itypat 123 mesh_size=pawtab(itypat)%mesh_size 124 indlmn => pawtab(itypat)%indlmn 125 126 ! loop over spin components 127 do ispden=1,nspden 128 129 ! loop over basis elements for this atom 130 ! ---- 131 do jlmn=1,pawtab(itypat)%lmn_size 132 jl=indlmn(1,jlmn) 133 jm=indlmn(2,jlmn) 134 jlm=indlmn(4,jlmn) 135 jln=indlmn(5,jlmn) 136 j0lmn=jlmn*(jlmn-1)/2 137 do ilmn=1,jlmn 138 il=indlmn(1,ilmn) 139 im=indlmn(2,ilmn) 140 iln=indlmn(5,ilmn) 141 ilm=indlmn(4,ilmn) 142 klmn=j0lmn+ilmn 143 kln = pawtab(itypat)%indklmn(2,klmn) ! need this for mesh selection below 144 145 if (jl==0 .and. il==0) then ! select only s-states 146 147 ! Loop over non-zero elements of rhoij 148 do irhoij=1,pawrhoij(iatom)%nrhoijsel 149 if (klmn==pawrhoij(iatom)%rhoijselect(irhoij)) then ! rho_ij /= 0 for this klmn 150 xxsum = 0 ! these three variables will be used to compute the slopes 151 xysumi = 0 152 xysumj = 0 153 do islope=1, nslope 154 xi=0 155 if(pawrad(itypat)%mesh_type == 1) xi = (islope - 1)*pawrad(itypat)%rstep 156 if(pawrad(itypat)%mesh_type == 2) xi = pawrad(itypat)%rstep * & 157 & (exp(pawrad(itypat)%lstep * (islope - 1)) - 1) 158 if(pawrad(itypat)%mesh_type == 3) then 159 if (islope == 1) then 160 xi = 0 161 else 162 xi = pawrad(itypat)%rstep * exp(pawrad(itypat)%lstep*(islope-1)) 163 end if 164 end if 165 if(pawrad(itypat)%mesh_type == 4) xi = & 166 & -pawrad(itypat)%rstep*log(1.0-(islope-1)/pawrad(itypat)%mesh_size) 167 yi = pawtab(itypat)%phi(islope,iln) ! function value for u_i 168 yj = pawtab(itypat)%phi(islope,jln) ! function value for u_j 169 xxsum = xxsum + xi*xi 170 xysumi = xysumi + xi*yi 171 xysumj = xysumj + xi*yj 172 end do 173 ! the slopes of the radial functions are obtained by minimizing 174 ! chi = sum(y_i - m*x_i)^2 (in other words, a linear least squares 175 ! fit constrained to go through the origin) 176 ! the result is m = sum(y_i*x_i)/sum(x_i*x_i) 177 mi = xysumi/xxsum 178 mj = xysumj/xxsum 179 ! accumulate the rho_ij contribution to the fermi contact for this spin density: 180 if (pawrhoij(iatom)%cplex == 1) then 181 fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+& 182 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(irhoij,ispden)*mi*mj/four_pi 183 else 184 fc(ispden,iatom_tot)=fc(ispden,iatom_tot)+& 185 & pawtab(itypat)%dltij(klmn)*pawrhoij(iatom)%rhoijp(2*irhoij-1,ispden)*mi*mj/four_pi 186 end if 187 end if ! end selection on klmn for nonzero rho_ij 188 end do ! end loop over nonzero rho_ij 189 end if ! end l=l'=0 selection 190 end do ! end loop over ilmn 191 end do ! end loop over jlmn 192 end do ! end loop over spin densities 193 end do ! Loop on atoms 194 195 !Reduction in case of parallelisation over atoms 196 if (paral_atom) then 197 call xmpi_sum(fc,my_comm_atom,ierr) 198 end if 199 200 !Destroy atom table used for parallelism 201 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 202 203 DBG_EXIT("COLL") 204 205 end subroutine make_fc_paw