TABLE OF CONTENTS
ABINIT/jellium [ Functions ]
NAME
jellium
FUNCTION
Optionally compute option=1 : local ionic (external) potential due to jellium background option=2 : contribution to the initial density taking into account the jellium slab
COPYRIGHT
Copyright (C) 2007-2018 ABINIT group (SC) 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/Infos/contributors.txt .
INPUTS
gmet(3,3)=metric tensor for G vecs (in bohr**-2) gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential) mpi_enreg=information about MPI parallelization nfft=(effective) number of FFT grid points (for this processor) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nspden=number of spin-density components option= rprimd(3,3)=dimensional primitive translations in real space (bohr) slabwsrad=Wigner-Seitz radius of jellium background slabzstart,slabzend=edges of jellium slab
OUTPUT
(if option==1) vjell(nfft)=external potential due to jellium background (if option==1) rhog(2,nfft), rhor(nfft,nspden)=density of positive charge in reciprocal, real space (only used in setvtr!)
SIDE EFFECTS
(if option==2) rhog(2,nfft), rhor(nfft,nspden)=reciprocal, real space updated initial electronic density
PARENTS
extraprho,gstate,setvtr
CHILDREN
fourdp
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 subroutine jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,nspden,& 56 & option,paral_kgb,slabwsrad,rhog,rhor,rprimd,vjell,slabzstart,slabzend) 57 58 use defs_basis 59 use defs_abitypes 60 use m_errors 61 use m_profiling_abi 62 63 !This section has been created automatically by the script Abilint (TD). 64 !Do not modify the following lines by hand. 65 #undef ABI_FUNC 66 #define ABI_FUNC 'jellium' 67 use interfaces_53_ffts 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 integer,intent(in) :: nfft,nspden,option,paral_kgb 75 real(dp),intent(in) :: gsqcut,slabwsrad,slabzend,slabzstart 76 type(MPI_type),intent(in) :: mpi_enreg 77 !arrays 78 integer,intent(in) :: ngfft(18) 79 real(dp),intent(in) :: gmet(3,3),rprimd(3,3) 80 real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,min(option,nspden)) 81 real(dp),intent(out) :: vjell(nfft) 82 83 !Local variables------------------------------- 84 !scalars 85 integer,parameter :: im=2,re=1 86 integer :: i1,i2,i3,id1,id2,id3,ig1,ig2,ig3,ii,ispden,n1,n2,n3,nfftot 87 real(dp),parameter :: tolfix=1.000000001_dp 88 real(dp) :: areaxy,cgz1,cgz2,cutoff,deltac,deltas,gsquar,gz,rhoave,sfi,sfr 89 real(dp) :: sgz1,sgz2,zcellength 90 character(len=500) :: message 91 !arrays 92 real(dp),allocatable :: rhjg(:,:),rhjr(:),vjelg(:,:) 93 94 ! ********************************************************************* 95 96 !Enforce that nspden<=2 97 if(nspden>2) then 98 MSG_ERROR('Jellium possible only with nspden <= 2.') 99 end if 100 101 !Make sure option is acceptable 102 if (option/=1 .and. option/=2) then 103 write(message, '(a,i0,3a)' )& 104 & 'option=',option,' is not allowed.',ch10,& 105 & 'Must be 1 or 2.' 106 MSG_BUG(message) 107 end if 108 109 zcellength=rprimd(3,3) 110 areaxy=abs(rprimd(1,1)*rprimd(2,2)-rprimd(1,2)*rprimd(2,1)) 111 rhoave=-half*three/(four_pi*slabwsrad**3) 112 113 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 114 id1=n1/2+2 115 id2=n2/2+2 116 id3=n3/2+2 117 nfftot=n1*n2*n3 118 cutoff=gsqcut*tolfix 119 120 ABI_ALLOCATE(rhjg,(2,nfft)) 121 ABI_ALLOCATE(rhjr,(nfft)) 122 rhjg(:,:)=zero 123 if(option==1) then 124 ABI_ALLOCATE(vjelg,(2,nfft)) 125 vjelg(:,:)=zero 126 end if 127 128 !Produce the potential due to the jellium background 129 ii=0 130 do i3=1,n3 131 ig3=i3-(i3/id3)*n3-1 132 do i2=1,n2 133 ig2=i2-(i2/id2)*n2-1 134 do i1=1,n1 135 136 ig1=i1-(i1/id1)*n1-1 137 ii=ii+1 138 gsquar=gsq_jel(ig1,ig2,ig3) 139 140 ! Skip G**2 outside cutoff and use \delta_{G_\|,0}: 141 if (gsquar<=cutoff.and.ig1==0.and.ig2==0) then 142 143 ! N o t e t h a t gz=two_pi*sqrt(gsq(0,0,ig3)) 144 gz=dble(ig3)*two_pi/zcellength 145 146 ! G_z == 0 147 if (ig3==0) then 148 sfr=two*rhoave*(slabzend-slabzstart) 149 sfi=zero 150 ! G_z /= 0 151 else ! of ig3==0 152 sgz2=sin(gz*slabzend) ; sgz1=sin(gz*slabzstart) 153 cgz2=cos(gz*slabzend) ; cgz1=cos(gz*slabzstart) 154 deltas=sgz2-sgz1 155 deltac=cgz2-cgz1 156 sfr=two*rhoave*deltas/gz 157 sfi=two*rhoave*deltac/gz 158 if(option==1) then 159 ! Assemble vjell_G 160 vjelg(re,ii)=four_pi*sfr/gz**2 161 vjelg(im,ii)=four_pi*sfi/gz**2 162 end if 163 end if ! of ig3==0 164 ! Assemble \rho_G 165 rhjg(re,ii)=sfr 166 rhjg(im,ii)=sfi 167 168 end if ! of gsquar ... 169 170 ! End loop on i1 171 end do 172 ! End loop on i2 173 end do 174 ! End loop on i3 175 end do 176 177 rhjg(:,:)=rhjg(:,:)/zcellength 178 if(option==1) vjelg(:,:)=vjelg(:,:)/zcellength 179 180 call fourdp(1,rhjg,rhjr,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 181 if(option==1) then 182 call fourdp(1,vjelg,vjell,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 183 rhog(:,:)=rhjg(:,:) 184 rhor(:,1)=rhjr(:) 185 else 186 ! Update the initial electronic density adding -rhjr 187 rhog(:,:)=rhog(:,:)-rhjg(:,:) 188 do ispden=1,nspden 189 rhor(:,ispden)=rhor(:,ispden)-rhjr(:)/dble(ispden) 190 end do 191 end if 192 193 ABI_DEALLOCATE(rhjg) 194 ABI_DEALLOCATE(rhjr) 195 if(option==1) then 196 ABI_DEALLOCATE(vjelg) 197 end if 198 199 !DEBUG 200 !write(std_out,*)' jellium : exit ' 201 !stop 202 !ENDDEBUG 203 204 contains 205 206 function gsq_jel(i1,i2,i3) 207 208 209 !This section has been created automatically by the script Abilint (TD). 210 !Do not modify the following lines by hand. 211 #undef ABI_FUNC 212 #define ABI_FUNC 'gsq_jel' 213 !End of the abilint section 214 215 real(dp) :: gsq_jel 216 integer,intent(in) :: i1,i2,i3 217 gsq_jel=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 218 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 219 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 220 end function gsq_jel 221 222 end subroutine jellium