Branch data 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 : : !!***
|
