TABLE OF CONTENTS


ABINIT/m_jellium [ Modules ]

[ Top ] [ Modules ]

NAME

  m_jellium

FUNCTION

  Routines related to jellium

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_jellium
28 
29  use defs_basis
30  use defs_abitypes
31  use m_errors
32  use m_abicore
33 
34  use m_fft,      only : fourdp
35 
36  implicit none
37 
38  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

PARENTS

      extraprho,gstate,setvtr

CHILDREN

      fourdp

SOURCE

 86 subroutine jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,nspden,&
 87 &  option,paral_kgb,slabwsrad,rhog,rhor,rprimd,vjell,slabzstart,slabzend)
 88 
 89 
 90 !This section has been created automatically by the script Abilint (TD).
 91 !Do not modify the following lines by hand.
 92 #undef ABI_FUNC
 93 #define ABI_FUNC 'jellium'
 94 !End of the abilint section
 95 
 96  implicit none
 97 
 98 !Arguments ------------------------------------
 99 !scalars
100  integer,intent(in) :: nfft,nspden,option,paral_kgb
101  real(dp),intent(in) :: gsqcut,slabwsrad,slabzend,slabzstart
102  type(MPI_type),intent(in) :: mpi_enreg
103 !arrays
104  integer,intent(in) :: ngfft(18)
105  real(dp),intent(in) :: gmet(3,3),rprimd(3,3)
106  real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,min(option,nspden))
107  real(dp),intent(out) :: vjell(nfft)
108 
109 !Local variables-------------------------------
110 !scalars
111  integer,parameter :: im=2,re=1
112  integer :: i1,i2,i3,id1,id2,id3,ig1,ig2,ig3,ii,ispden,n1,n2,n3,nfftot
113  real(dp),parameter :: tolfix=1.000000001_dp
114  real(dp) :: areaxy,cgz1,cgz2,cutoff,deltac,deltas,gsquar,gz,rhoave,sfi,sfr
115  real(dp) :: sgz1,sgz2,zcellength
116  character(len=500) :: message
117 !arrays
118  real(dp),allocatable :: rhjg(:,:),rhjr(:),vjelg(:,:)
119 
120 ! *********************************************************************
121 
122 !Enforce that nspden<=2
123  if(nspden>2) then
124    MSG_ERROR('Jellium possible only with nspden <= 2.')
125  end if
126 
127 !Make sure option is acceptable
128  if (option/=1 .and. option/=2) then
129    write(message, '(a,i0,3a)' )&
130 &   'option=',option,' is not allowed.',ch10,&
131 &   'Must be 1 or 2.'
132    MSG_BUG(message)
133  end if
134 
135  zcellength=rprimd(3,3)
136  areaxy=abs(rprimd(1,1)*rprimd(2,2)-rprimd(1,2)*rprimd(2,1))
137  rhoave=-half*three/(four_pi*slabwsrad**3)
138 
139  n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
140  id1=n1/2+2
141  id2=n2/2+2
142  id3=n3/2+2
143  nfftot=n1*n2*n3
144  cutoff=gsqcut*tolfix
145 
146  ABI_ALLOCATE(rhjg,(2,nfft))
147  ABI_ALLOCATE(rhjr,(nfft))
148  rhjg(:,:)=zero
149  if(option==1) then
150    ABI_ALLOCATE(vjelg,(2,nfft))
151    vjelg(:,:)=zero
152  end if
153 
154 !Produce the potential due to the jellium background
155  ii=0
156  do i3=1,n3
157    ig3=i3-(i3/id3)*n3-1
158    do i2=1,n2
159      ig2=i2-(i2/id2)*n2-1
160      do i1=1,n1
161 
162        ig1=i1-(i1/id1)*n1-1
163        ii=ii+1
164        gsquar=gsq_jel(ig1,ig2,ig3)
165 
166 !      Skip G**2 outside cutoff and use \delta_{G_\|,0}:
167        if (gsquar<=cutoff.and.ig1==0.and.ig2==0) then
168 
169 !        N o t e   t h a t   gz=two_pi*sqrt(gsq(0,0,ig3))
170          gz=dble(ig3)*two_pi/zcellength
171 
172 !        G_z == 0
173          if (ig3==0) then
174            sfr=two*rhoave*(slabzend-slabzstart)
175            sfi=zero
176 !          G_z /= 0
177          else ! of ig3==0
178            sgz2=sin(gz*slabzend) ; sgz1=sin(gz*slabzstart)
179            cgz2=cos(gz*slabzend) ; cgz1=cos(gz*slabzstart)
180            deltas=sgz2-sgz1
181            deltac=cgz2-cgz1
182            sfr=two*rhoave*deltas/gz
183            sfi=two*rhoave*deltac/gz
184            if(option==1) then
185 !            Assemble vjell_G
186              vjelg(re,ii)=four_pi*sfr/gz**2
187              vjelg(im,ii)=four_pi*sfi/gz**2
188            end if
189          end if ! of ig3==0
190 !        Assemble \rho_G
191          rhjg(re,ii)=sfr
192          rhjg(im,ii)=sfi
193 
194        end if ! of gsquar ...
195 
196 !      End loop on i1
197      end do
198 !    End loop on i2
199    end do
200 !  End loop on i3
201  end do
202 
203  rhjg(:,:)=rhjg(:,:)/zcellength
204  if(option==1) vjelg(:,:)=vjelg(:,:)/zcellength
205 
206  call fourdp(1,rhjg,rhjr,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
207  if(option==1) then
208    call fourdp(1,vjelg,vjell,1,mpi_enreg,nfft,ngfft,paral_kgb,0)
209    rhog(:,:)=rhjg(:,:)
210    rhor(:,1)=rhjr(:)
211  else
212 !  Update the initial electronic density adding -rhjr
213    rhog(:,:)=rhog(:,:)-rhjg(:,:)
214    do ispden=1,nspden
215      rhor(:,ispden)=rhor(:,ispden)-rhjr(:)/dble(ispden)
216    end do
217  end if
218 
219  ABI_DEALLOCATE(rhjg)
220  ABI_DEALLOCATE(rhjr)
221  if(option==1) then
222    ABI_DEALLOCATE(vjelg)
223  end if
224 
225 !DEBUG
226 !write(std_out,*)' jellium : exit '
227 !stop
228 !ENDDEBUG
229 
230  contains
231 
232    function gsq_jel(i1,i2,i3)
233 
234 
235 !This section has been created automatically by the script Abilint (TD).
236 !Do not modify the following lines by hand.
237 #undef ABI_FUNC
238 #define ABI_FUNC 'gsq_jel'
239 !End of the abilint section
240 
241    real(dp) :: gsq_jel
242    integer,intent(in) :: i1,i2,i3
243    gsq_jel=dble(i1*i1)*gmet(1,1)+dble(i2*i2)*gmet(2,2)+&
244 &   dble(i3*i3)*gmet(3,3)+dble(2*i1*i2)*gmet(1,2)+&
245 &   dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
246  end function gsq_jel
247 
248 end subroutine jellium