TABLE OF CONTENTS
ABINIT/m_jellium [ Modules ]
NAME
m_jellium
FUNCTION
Routines related to jellium
COPYRIGHT
Copyright (C) 2007-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_jellium 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 use defs_abitypes, only : MPI_type 29 use m_fft, only : fourdp 30 31 implicit none 32 33 private
m_jellium/jellium [ Functions ]
[ Top ] [ m_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
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
SOURCE
75 subroutine jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,nspden,& 76 & option,slabwsrad,rhog,rhor,rprimd,vjell,slabzstart,slabzend) 77 78 !Arguments ------------------------------------ 79 !scalars 80 integer,intent(in) :: nfft,nspden,option 81 real(dp),intent(in) :: gsqcut,slabwsrad,slabzend,slabzstart 82 type(MPI_type),intent(in) :: mpi_enreg 83 !arrays 84 integer,intent(in) :: ngfft(18) 85 real(dp),intent(in) :: gmet(3,3),rprimd(3,3) 86 real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,min(option,nspden)) 87 real(dp),intent(out) :: vjell(nfft) 88 89 !Local variables------------------------------- 90 !scalars 91 integer,parameter :: im=2,re=1 92 integer :: i1,i2,i3,id1,id2,id3,ig1,ig2,ig3,ii,ispden,n1,n2,n3,nfftot 93 real(dp),parameter :: tolfix=1.000000001_dp 94 real(dp) :: areaxy,cgz1,cgz2,cutoff,deltac,deltas,gsquar,gz,rhoave,sfi,sfr 95 real(dp) :: sgz1,sgz2,zcellength 96 character(len=500) :: message 97 !arrays 98 real(dp),allocatable :: rhjg(:,:),rhjr(:),vjelg(:,:) 99 100 ! ********************************************************************* 101 102 !Enforce that nspden<=2 103 if(nspden>2) then 104 ABI_ERROR('Jellium possible only with nspden <= 2.') 105 end if 106 107 !Make sure option is acceptable 108 if (option/=1 .and. option/=2) then 109 write(message, '(a,i0,3a)' )& 110 & 'option=',option,' is not allowed.',ch10,& 111 & 'Must be 1 or 2.' 112 ABI_BUG(message) 113 end if 114 115 zcellength=rprimd(3,3) 116 areaxy=abs(rprimd(1,1)*rprimd(2,2)-rprimd(1,2)*rprimd(2,1)) 117 rhoave=-half*three/(four_pi*slabwsrad**3) 118 119 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 120 id1=n1/2+2 121 id2=n2/2+2 122 id3=n3/2+2 123 nfftot=n1*n2*n3 124 cutoff=gsqcut*tolfix 125 126 ABI_MALLOC(rhjg,(2,nfft)) 127 ABI_MALLOC(rhjr,(nfft)) 128 rhjg(:,:)=zero 129 if(option==1) then 130 ABI_MALLOC(vjelg,(2,nfft)) 131 vjelg(:,:)=zero 132 end if 133 134 !Produce the potential due to the jellium background 135 ii=0 136 do i3=1,n3 137 ig3=i3-(i3/id3)*n3-1 138 do i2=1,n2 139 ig2=i2-(i2/id2)*n2-1 140 do i1=1,n1 141 142 ig1=i1-(i1/id1)*n1-1 143 ii=ii+1 144 gsquar=gsq_jel(ig1,ig2,ig3) 145 146 ! Skip G**2 outside cutoff and use \delta_{G_\|,0}: 147 if (gsquar<=cutoff.and.ig1==0.and.ig2==0) then 148 149 ! N o t e t h a t gz=two_pi*sqrt(gsq(0,0,ig3)) 150 gz=dble(ig3)*two_pi/zcellength 151 152 ! G_z == 0 153 if (ig3==0) then 154 sfr=two*rhoave*(slabzend-slabzstart) 155 sfi=zero 156 ! G_z /= 0 157 else ! of ig3==0 158 sgz2=sin(gz*slabzend) ; sgz1=sin(gz*slabzstart) 159 cgz2=cos(gz*slabzend) ; cgz1=cos(gz*slabzstart) 160 deltas=sgz2-sgz1 161 deltac=cgz2-cgz1 162 sfr=two*rhoave*deltas/gz 163 sfi=two*rhoave*deltac/gz 164 if(option==1) then 165 ! Assemble vjell_G 166 vjelg(re,ii)=four_pi*sfr/gz**2 167 vjelg(im,ii)=four_pi*sfi/gz**2 168 end if 169 end if ! of ig3==0 170 ! Assemble \rho_G 171 rhjg(re,ii)=sfr 172 rhjg(im,ii)=sfi 173 174 end if ! of gsquar ... 175 176 ! End loop on i1 177 end do 178 ! End loop on i2 179 end do 180 ! End loop on i3 181 end do 182 183 rhjg(:,:)=rhjg(:,:)/zcellength 184 if(option==1) vjelg(:,:)=vjelg(:,:)/zcellength 185 186 call fourdp(1,rhjg,rhjr,1,mpi_enreg,nfft,1,ngfft,0) 187 if(option==1) then 188 call fourdp(1,vjelg,vjell,1,mpi_enreg,nfft,1,ngfft,0) 189 rhog(:,:)=rhjg(:,:) 190 rhor(:,1)=rhjr(:) 191 else 192 ! Update the initial electronic density adding -rhjr 193 rhog(:,:)=rhog(:,:)-rhjg(:,:) 194 do ispden=1,nspden 195 rhor(:,ispden)=rhor(:,ispden)-rhjr(:)/dble(ispden) 196 end do 197 end if 198 199 ABI_FREE(rhjg) 200 ABI_FREE(rhjr) 201 if(option==1) then 202 ABI_FREE(vjelg) 203 end if 204 205 !DEBUG 206 !write(std_out,*)' jellium : exit ' 207 !stop 208 !ENDDEBUG 209 210 contains 211 212 function gsq_jel(i1,i2,i3) 213 214 real(dp) :: gsq_jel 215 integer,intent(in) :: i1,i2,i3 216 gsq_jel=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+& 217 & dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+& 218 & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1) 219 end function gsq_jel 220 221 end subroutine jellium