Line data Source code
1 : !!****m* ABINIT/m_jellium
2 : !! NAME
3 : !! m_jellium
4 : !!
5 : !! FUNCTION
6 : !! Routines related to jellium
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2007-2026 ABINIT group (SC)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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
34 : !!***
35 :
36 : public :: jellium
37 : !!***
38 :
39 : contains
40 : !!***
41 :
42 : !!****f* m_jellium/jellium
43 : !! NAME
44 : !! jellium
45 : !!
46 : !! FUNCTION
47 : !! Optionally compute
48 : !! option=1 : local ionic (external) potential due to jellium background
49 : !! option=2 : contribution to the initial density taking into account
50 : !! the jellium slab
51 : !!
52 : !! INPUTS
53 : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
54 : !! gsqcut=cutoff on (k+G)^2 (bohr^-2) (sphere for density and potential)
55 : !! mpi_enreg=information about MPI parallelization
56 : !! nfft=(effective) number of FFT grid points (for this processor)
57 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
58 : !! nspden=number of spin-density components
59 : !! option=
60 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
61 : !! slabwsrad=Wigner-Seitz radius of jellium background
62 : !! slabzstart,slabzend=edges of jellium slab
63 : !!
64 : !! OUTPUT
65 : !! (if option==1) vjell(nfft)=external potential due to jellium background
66 : !! (if option==1) rhog(2,nfft), rhor(nfft,nspden)=density of positive charge
67 : !! in reciprocal, real space (only used in setvtr!)
68 : !!
69 : !! SIDE EFFECTS
70 : !! (if option==2) rhog(2,nfft), rhor(nfft,nspden)=reciprocal, real space
71 : !! updated initial electronic density
72 : !!
73 : !! SOURCE
74 :
75 12 : subroutine jellium(gmet,gsqcut,mpi_enreg,nfft,ngfft,nspden,&
76 12 : & 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 12 : real(dp),allocatable :: rhjg(:,:),rhjr(:),vjelg(:,:)
99 :
100 : ! *********************************************************************
101 :
102 : !Enforce that nspden<=2
103 12 : if(nspden>2) then
104 0 : ABI_ERROR('Jellium possible only with nspden <= 2.')
105 : end if
106 :
107 : !Make sure option is acceptable
108 12 : if (option/=1 .and. option/=2) then
109 : write(message, '(a,i0,3a)' )&
110 0 : & 'option=',option,' is not allowed.',ch10,&
111 0 : & 'Must be 1 or 2.'
112 0 : ABI_BUG(message)
113 : end if
114 :
115 12 : zcellength=rprimd(3,3)
116 12 : areaxy=abs(rprimd(1,1)*rprimd(2,2)-rprimd(1,2)*rprimd(2,1))
117 12 : rhoave=-half*three/(four_pi*slabwsrad**3)
118 :
119 12 : n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3)
120 12 : id1=n1/2+2
121 12 : id2=n2/2+2
122 12 : id3=n3/2+2
123 12 : nfftot=n1*n2*n3
124 12 : cutoff=gsqcut*tolfix
125 :
126 36 : ABI_MALLOC(rhjg,(2,nfft))
127 36 : ABI_MALLOC(rhjr,(nfft))
128 629292 : rhjg(:,:)=zero
129 12 : if(option==1) then
130 12 : ABI_MALLOC(vjelg,(2,nfft))
131 314646 : vjelg(:,:)=zero
132 : end if
133 :
134 : !Produce the potential due to the jellium background
135 12 : ii=0
136 540 : do i3=1,n3
137 528 : ig3=i3-(i3/id3)*n3-1
138 10572 : do i2=1,n2
139 10032 : ig2=i2-(i2/id2)*n2-1
140 220320 : do i1=1,n1
141 :
142 209760 : ig1=i1-(i1/id1)*n1-1
143 209760 : ii=ii+1
144 209760 : gsquar=gsq_jel(ig1,ig2,ig3)
145 :
146 : ! Skip G**2 outside cutoff and use \delta_{G_\|,0}:
147 219792 : 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 500 : gz=dble(ig3)*two_pi/zcellength
151 :
152 : ! G_z == 0
153 500 : if (ig3==0) then
154 12 : sfr=two*rhoave*(slabzend-slabzstart)
155 12 : sfi=zero
156 : ! G_z /= 0
157 : else ! of ig3==0
158 488 : sgz2=sin(gz*slabzend) ; sgz1=sin(gz*slabzstart)
159 488 : cgz2=cos(gz*slabzend) ; cgz1=cos(gz*slabzstart)
160 488 : deltas=sgz2-sgz1
161 488 : deltac=cgz2-cgz1
162 488 : sfr=two*rhoave*deltas/gz
163 488 : sfi=two*rhoave*deltac/gz
164 488 : if(option==1) then
165 : ! Assemble vjell_G
166 244 : vjelg(re,ii)=four_pi*sfr/gz**2
167 244 : vjelg(im,ii)=four_pi*sfi/gz**2
168 : end if
169 : end if ! of ig3==0
170 : ! Assemble \rho_G
171 500 : rhjg(re,ii)=sfr
172 500 : 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 629292 : rhjg(:,:)=rhjg(:,:)/zcellength
184 314652 : if(option==1) vjelg(:,:)=vjelg(:,:)/zcellength
185 :
186 12 : call fourdp(1,rhjg,rhjr,1,mpi_enreg,nfft,1,ngfft,0)
187 12 : if(option==1) then
188 6 : call fourdp(1,vjelg,vjell,1,mpi_enreg,nfft,1,ngfft,0)
189 314646 : rhog(:,:)=rhjg(:,:)
190 104886 : rhor(:,1)=rhjr(:)
191 : else
192 : ! Update the initial electronic density adding -rhjr
193 314646 : rhog(:,:)=rhog(:,:)-rhjg(:,:)
194 12 : do ispden=1,nspden
195 104892 : rhor(:,ispden)=rhor(:,ispden)-rhjr(:)/dble(ispden)
196 : end do
197 : end if
198 :
199 12 : ABI_FREE(rhjg)
200 12 : ABI_FREE(rhjr)
201 18 : if(option==1) then
202 6 : ABI_FREE(vjelg)
203 : end if
204 :
205 : !DEBUG
206 : !write(std_out,*)' jellium : exit '
207 : !stop
208 : !ENDDEBUG
209 :
210 : contains
211 :
212 209760 : 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 209760 : & dble(2*i2*i3)*gmet(2,3)+dble(2*i3*i1)*gmet(3,1)
219 209760 : end function gsq_jel
220 :
221 : end subroutine jellium
222 : !!***
223 :
224 : end module m_jellium
225 : !!***
|
