TABLE OF CONTENTS
ABINIT/hartre [ Functions ]
NAME
hartre
FUNCTION
Given rho(G), compute Hartree potential (=FFT of rho(G)/pi/(G+q)**2) When cplex=1, assume q=(0 0 0), and vhartr will be REAL When cplex=2, q must be taken into account, and vhartr will be COMPLEX
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR). 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/doc/developers/contributors.txt .
NOTES
*Modified code to avoid if statements inside loops to skip G=0. Replaced if statement on G^2>gsqcut to skip G s outside where rho(G) should be 0. Effect is negligible but gsqcut should be used to be strictly consistent with usage elsewhere in code. *The speed-up is provided by doing a few precomputations outside the inner loop. One variable size array is needed for this (gq).
INPUTS
cplex= if 1, vhartr is REAL, if 2, vhartr is COMPLEX gsqcut=cutoff value on G**2 for sphere inside fft box. (gsqcut=(boxcut**2)*ecut/(2.d0*(Pi**2)) izero=if 1, unbalanced components of Vhartree(g) are set to zero 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 [qpt(3)=reduced coordinates for a wavevector to be combined with the G vectors (needed if cplex==2).] rhog(2,nfft)=electron density in G space rprimd(3,3)=dimensional primitive translations in real space (bohr) divgq0= [optional argument] value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq
OUTPUT
vhartr(cplex*nfft)=Hartree potential in real space, either REAL or COMPLEX
PARENTS
dfpt_rhotov,dfptnl_loop,energy,fock_getghc,m_kxc,nonlinear,nres2vres odamix,prcref,prcref_PMA,respfn,rhotov,setup_positron,setvtr,tddft
CHILDREN
fourdp,metric,ptabs_fourdp,timab,zerosym
SOURCE
51 #if defined HAVE_CONFIG_H 52 #include "config.h" 53 #endif 54 55 #include "abi_common.h" 56 57 subroutine hartre(cplex,gsqcut,izero,mpi_enreg,nfft,ngfft,paral_kgb,rhog,rprimd,vhartr,& 58 & divgq0,qpt) ! Optional argument 59 60 use defs_basis 61 use defs_abitypes 62 use m_errors 63 use m_profiling_abi 64 65 use m_geometry, only : metric 66 use m_mpinfo, only : ptabs_fourdp 67 use m_fft, only : zerosym 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'hartre' 73 use interfaces_18_timing 74 use interfaces_53_ffts 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 !scalars 81 integer,intent(in) :: cplex,izero,nfft,paral_kgb 82 real(dp),intent(in) :: gsqcut 83 type(MPI_type),intent(in) :: mpi_enreg 84 !arrays 85 integer,intent(in) :: ngfft(18) 86 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft) 87 real(dp),intent(in),optional :: divgq0 88 real(dp),intent(in),optional :: qpt(3) 89 real(dp),intent(out) :: vhartr(cplex*nfft) 90 91 !Local variables------------------------------- 92 !scalars 93 integer,parameter :: im=2,re=1 94 integer :: i1,i2,i23,i2_local,i3,id1,id2,id3 95 integer :: ig,ig1min,ig1,ig1max,ig2,ig2min,ig2max,ig3,ig3min,ig3max 96 integer :: ii,ii1,ing,n1,n2,n3,qeq0,qeq05,me_fft,nproc_fft 97 real(dp),parameter :: tolfix=1.000000001e0_dp 98 real(dp) :: cutoff,den,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3,ucvol 99 character(len=500) :: message 100 !arrays 101 integer :: id(3) 102 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 103 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 104 real(dp) :: gmet(3,3),gprimd(3,3),qpt_(3),rmet(3,3),tsec(2) 105 real(dp),allocatable :: gq(:,:),work1(:,:) 106 107 ! ************************************************************************* 108 109 ! Keep track of total time spent in hartre 110 call timab(10,1,tsec) 111 112 ! Check that cplex has an allowed value 113 if(cplex/=1 .and. cplex/=2)then 114 write(message, '(a,i0,a,a)' )& 115 & 'From the calling routine, cplex=',cplex,ch10,& 116 & 'but the only value allowed are 1 and 2.' 117 MSG_BUG(message) 118 end if 119 120 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 121 122 n1=ngfft(1); n2=ngfft(2); n3=ngfft(3) 123 nproc_fft = mpi_enreg%nproc_fft; me_fft = mpi_enreg%me_fft 124 125 ! Get the distrib associated with this fft_grid 126 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 127 128 ! Initialize a few quantities 129 cutoff=gsqcut*tolfix 130 if(present(qpt))then 131 qpt_=qpt 132 else 133 qpt_=zero 134 end if 135 qeq0=0 136 if(qpt_(1)**2+qpt_(2)**2+qpt_(3)**2<1.d-15) qeq0=1 137 qeq05=0 138 if (qeq0==0) then 139 if (abs(abs(qpt_(1))-half)<tol12.or.abs(abs(qpt_(2))-half)<tol12.or. & 140 & abs(abs(qpt_(3))-half)<tol12) qeq05=1 141 end if 142 143 ! If cplex=1 then qpt_ should be 0 0 0 144 if (cplex==1.and. qeq0/=1) then 145 write(message,'(a,3e12.4,a,a)')& 146 & 'cplex=1 but qpt=',qpt_,ch10,& 147 & 'qpt should be 0 0 0.' 148 MSG_BUG(message) 149 end if 150 151 ! If FFT parallelism then qpt should not be 1/2 152 if (nproc_fft>1.and.qeq05==1) then 153 write(message, '(a,3e12.4,a,a)' )& 154 & 'FFT parallelism selected but qpt',qpt_,ch10,& 155 & 'qpt(i) should not be 1/2...' 156 MSG_ERROR(message) 157 end if 158 159 ! In order to speed the routine, precompute the components of g+q 160 ! Also check if the booked space was large enough... 161 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 162 do ii=1,3 163 id(ii)=ngfft(ii)/2+2 164 do ing=1,ngfft(ii) 165 ig=ing-(ing/id(ii))*ngfft(ii)-1 166 gq(ii,ing)=ig+qpt_(ii) 167 end do 168 end do 169 ig1max=-1;ig2max=-1;ig3max=-1 170 ig1min=n1;ig2min=n2;ig3min=n3 171 172 ABI_ALLOCATE(work1,(2,nfft)) 173 id1=n1/2+2;id2=n2/2+2;id3=n3/2+2 174 175 ! Triple loop on each dimension 176 do i3=1,n3 177 ig3=i3-(i3/id3)*n3-1 178 ! Precompute some products that do not depend on i2 and i1 179 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 180 gqgm23=gq(3,i3)*gmet(2,3)*2 181 gqgm13=gq(3,i3)*gmet(1,3)*2 182 183 do i2=1,n2 184 ig2=i2-(i2/id2)*n2-1 185 if (fftn2_distrib(i2) == me_fft) then 186 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 187 gqgm12=gq(2,i2)*gmet(1,2)*2 188 gqg2p3=gqgm13+gqgm12 189 190 i2_local = ffti2_local(i2) 191 i23=n1*(i2_local-1 +(n2/nproc_fft)*(i3-1)) 192 ! Do the test that eliminates the Gamma point outside of the inner loop 193 ii1=1 194 if(i23==0 .and. qeq0==1 .and. ig2==0 .and. ig3==0)then 195 ii1=2 196 work1(re,1+i23)=zero 197 work1(im,1+i23)=zero 198 ! If the value of the integration of the Coulomb singularity 4pi\int_BZ 1/q^2 dq is given, use it 199 if (PRESENT(divgq0)) then 200 work1(re,1+i23)=rhog(re,1+i23)*divgq0*piinv 201 work1(im,1+i23)=rhog(im,1+i23)*divgq0*piinv 202 end if 203 end if 204 205 ! Final inner loop on the first dimension (note the lower limit) 206 do i1=ii1,n1 207 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 208 ii=i1+i23 209 210 if(gs<=cutoff)then 211 ! Identify min/max indexes (to cancel unbalanced contributions later) 212 ! Count (q+g)-vectors with similar norm 213 if ((qeq05==1).and.(izero==1)) then 214 ig1=i1-(i1/id1)*n1-1 215 ig1max=max(ig1max,ig1); ig1min=min(ig1min,ig1) 216 ig2max=max(ig2max,ig2); ig2min=min(ig2min,ig2) 217 ig3max=max(ig3max,ig3); ig3min=min(ig3min,ig3) 218 end if 219 220 den=piinv/gs 221 work1(re,ii)=rhog(re,ii)*den 222 work1(im,ii)=rhog(im,ii)*den 223 else 224 ! gs>cutoff 225 work1(re,ii)=zero 226 work1(im,ii)=zero 227 end if 228 229 end do ! End loop on i1 230 end if 231 end do ! End loop on i2 232 end do ! End loop on i3 233 234 ABI_DEALLOCATE(gq) 235 236 if (izero==1) then 237 ! Set contribution of unbalanced components to zero 238 239 if (qeq0==1) then !q=0 240 call zerosym(work1,2,n1,n2,n3,comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 241 242 else if (qeq05==1) then 243 !q=1/2; this doesn't work in parallel 244 ig1=-1;if (mod(n1,2)==0) ig1=1+n1/2 245 ig2=-1;if (mod(n2,2)==0) ig2=1+n2/2 246 ig3=-1;if (mod(n3,2)==0) ig3=1+n3/2 247 if (abs(abs(qpt_(1))-half)<tol12) then 248 if (abs(ig1min)<abs(ig1max)) ig1=abs(ig1max) 249 if (abs(ig1min)>abs(ig1max)) ig1=n1-abs(ig1min) 250 end if 251 if (abs(abs(qpt_(2))-half)<tol12) then 252 if (abs(ig2min)<abs(ig2max)) ig2=abs(ig2max) 253 if (abs(ig2min)>abs(ig2max)) ig2=n2-abs(ig2min) 254 end if 255 if (abs(abs(qpt_(3))-half)<tol12) then 256 if (abs(ig3min)<abs(ig3max)) ig3=abs(ig3max) 257 if (abs(ig3min)>abs(ig3max)) ig3=n3-abs(ig3min) 258 end if 259 call zerosym(work1,2,n1,n2,n3,ig1=ig1,ig2=ig2,ig3=ig3,& 260 & comm_fft=mpi_enreg%comm_fft,distribfft=mpi_enreg%distribfft) 261 end if 262 end if 263 264 ! Fourier Transform Vhartree. Vh in reciprocal space was stored in work1 265 call fourdp(cplex,work1,vhartr,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 266 267 ABI_DEALLOCATE(work1) 268 269 call timab(10,2,tsec) 270 271 end subroutine hartre