TABLE OF CONTENTS
ABINIT/calc_smeared_density [ Functions ]
NAME
calc_smeared_density
FUNCTION
Calculate the smeared charge density on the an FFT grid in real space.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (GMR, VO, LR, RWG, MG, RShaltaf) 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
rhor - the density on the real space FFT grid kappa_strategy - integer indicating strategy for the smearing parameter 1 - > fixed kappa 2 - > take kappa as n-dependent nfftf - contains all needed information about the 3D FFT grid ngfftf(18) - contains all needed information about the fine 3D FFT grid npw - number of plane waves mpi_enreg - Datatype with parallellisation info paral_kgb - band parallelism variable kappa_in - optional variable setting a fixed value of kappa
OUTPUT
rhotilder - the smeared density on the real space FFT grid
NOTES
This routine transforms a given density rhor(R) -> rhog(G) and calculates
smeared density. This is given by the formula:
\tilde{n}(\mathbf{G}) = n(\mathbf{G})\frac{{\kappa}^2}{{\kappa}^2+|\mathbf{G}|^2}
which corresponds to a convolution of the density with a Yukawa function
in real space:
\tilde{n}(\mathbf{r})=\int \mathrm{d}r'\,\frac{{\kappa}^2}{4\pi} \\
\frac{e^{-{\kappa}|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r'}|} \\
n(\mathbf{r}')
When the convolution has been performed in reciprocal space, the function
is transformed back to real space, giving the smeared density in real space.
PARENTS
CHILDREN
fourdp
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 subroutine calc_smeared_density(rhor,kappa_strategy,rhotilder,nfftf,ngfftf,npw,& 61 & gvec,gprimd,ucvol,mpi_enreg,paral_kgb,kappa_in) 62 63 use defs_basis 64 use defs_abitypes 65 use m_errors 66 use m_profiling_abi 67 68 use m_fft_mesh, only : g2ifft 69 70 !This section has been created automatically by the script Abilint (TD). 71 !Do not modify the following lines by hand. 72 #undef ABI_FUNC 73 #define ABI_FUNC 'calc_smeared_density' 74 use interfaces_53_ffts 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 !scalars 81 integer, intent(in) :: kappa_strategy,nfftf,npw,paral_kgb 82 !integer, optional, intent(in) :: prtvol 83 real(dp), intent(in) :: ucvol 84 real(dp), intent(in), optional :: kappa_in 85 !arrays 86 integer,intent(in) :: ngfftf(18),gvec(3,npw) 87 real(dp), intent(in) :: gprimd(3,3) 88 real(dp), intent(inout) :: rhor(nfftf) 89 real(dp), intent(out) :: rhotilder(nfftf) 90 ! types 91 type(MPI_type),intent(in) :: mpi_enreg 92 93 !Local variables ------------------------------ 94 !scalars 95 integer :: ig,igp,ii,gmgp_idx,ierr 96 real(dp), parameter :: two_pi_sq = two_pi*two_pi 97 real(dp) :: inv_kappasq,abs_gmgp_sq,yukawa_factor,yukawa_denom 98 !arrays 99 real(dp) :: gmet(3,3) 100 integer :: gmgp(3) 101 real(dp), allocatable :: rhog_dp(:,:),rhotildeg_dp(:,:) 102 character(len=500) :: msg 103 104 !************************************************************************* 105 106 DBG_ENTER("COLL") 107 108 ABI_ALLOCATE(rhog_dp,(2,nfftf)) 109 ABI_ALLOCATE(rhotildeg_dp,(2,nfftf)) 110 rhog_dp=zero; rhotildeg_dp=zero 111 112 !write(std_out,*) ' calc. sm. den., npw=',npw 113 !write(std_out,*) ' calc. sm. den., gprimd=',gprimd 114 !write(std_out,*) ' calc. sm. den., ngfftf(1:3)=',ngfftf(1:3) 115 !write(std_out,*) ' calc. sm. den., nfftf=',nfftf 116 117 if (present(kappa_in)) write(std_out,*) ' calc. sm. den., kappa_in=',kappa_in 118 119 call fourdp(1,rhog_dp,rhor,-1,mpi_enreg,nfftf,ngfftf,paral_kgb,0) 120 121 !Compute reciprocal space metrics 122 do ii=1,3 123 gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+& 124 & gprimd(2,ii)*gprimd(2,:)+& 125 & gprimd(3,ii)*gprimd(3,:) 126 end do 127 128 !Check which strategy for determining kappa is used 129 select case(kappa_strategy) 130 131 case(1) ! kappa is a constant, kappa_in is provided in the routine call 132 133 if (present(kappa_in)) then 134 inv_kappasq = 1.0_dp/(kappa_in*kappa_in) 135 else 136 MSG_ERROR(' Argument kappa_in must be set with kappa_strategy=1') 137 end if 138 139 ! Perform the multiplication 140 ierr=ierr+1 141 do ig=1,npw 142 do igp=1,npw 143 ! Calculate |G-G'| 144 gmgp(:) = gvec(:,ig)-gvec(:,igp) 145 ! Make sure the vector is not outside the FFT box 146 if (ANY (gmgp(:)>ngfftf(1:3)/2 .or. gmgp(:)<-(ngfftf(1:3)-1)/2) ) then 147 ierr=ierr+1 148 cycle 149 end if 150 abs_gmgp_sq = two_pi_sq*dot_product(gmgp,MATMUL(gmet,gmgp)) 151 ! Calculate Yukawa factor kappa^2/(kappa^2+|G-G'|^2) 152 ! = 1/(1+(|G-G'|/kappa)^2) 153 yukawa_denom = one + abs_gmgp_sq*inv_kappasq 154 yukawa_factor = one/yukawa_denom 155 ! Calculate smeared density in G-space i.e. n(G-G')*kappa^2/(kappa^2+|G-G'|^2) 156 gmgp_idx = g2ifft(gmgp,ngfftf) 157 ! rhotildeg_dp(:,igfft(ig,igp)) = rhog_dp(:,igfft(ig,igp))*yukawa_factor 158 rhotildeg_dp(:,gmgp_idx) = rhog_dp(:,gmgp_idx)*yukawa_factor 159 end do 160 end do 161 162 if (ierr/=0) then 163 write(msg,'(a,i4,3a)')& 164 & 'Found ',ierr,' G1-G2 vectors falling outside the FFT box. ',ch10,& 165 & 'Enlarge the FFT mesh to get rid of this problem. ' 166 MSG_WARNING(msg) 167 end if 168 169 case(2) ! kappa is a function of the density 170 MSG_ERROR('kappa_strategy=2 not coded yet') 171 172 case default 173 MSG_ERROR('error in kappa_strategy') 174 end select 175 176 !Transform the smeared density back to real space 177 178 call fourdp(1,rhotildeg_dp,rhotilder,1,mpi_enreg,nfftf,ngfftf,paral_kgb,0) 179 180 !Debug 181 abs_gmgp_sq = SUM(rhotilder(:))*ucvol/nfftf 182 write(msg,'(a,ES29.10E3)') ' Integral of rhotilde(r):',abs_gmgp_sq 183 MSG_WARNING(msg) 184 abs_gmgp_sq = SUM(ABS(rhor(:)))*ucvol/nfftf 185 write(msg,'(a,ES29.10E3)') ' Integral of rho(r):',abs_gmgp_sq 186 MSG_WARNING(msg) 187 abs_gmgp_sq = SUM(rhotilder(:)-rhor(:))*ucvol/nfftf 188 write(msg,'(a,ES29.10E3)') ' Integral of difference rhotilde(r)-rho(r):',abs_gmgp_sq 189 MSG_WARNING(msg) 190 abs_gmgp_sq = SUM(ABS(rhotilder(:)-rhor(:)))*ucvol/nfftf 191 write(msg,'(a,ES29.10E3)') ' Integral of abs. difference |rhotilde(r)-rho(r)|:', abs_gmgp_sq 192 MSG_WARNING(msg) 193 194 ABI_DEALLOCATE(rhog_dp) 195 ABI_DEALLOCATE(rhotildeg_dp) 196 197 DBG_EXIT("COLL") 198 199 end subroutine calc_smeared_density