TABLE OF CONTENTS
ABINIT/hartrestr [ Functions ]
NAME
hartrestr
FUNCTION
To be called for strain perturbation only Compute the inhomogenous terms generated by the strain derivative of Hartree potential due to the ground state charge rho(G) FFT of (rho(G)/pi)*[d(1/G**2)/d(strain) - delta(diagonal strain)*(1/G**2)]
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, 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
*based largely on hartre.f *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
gsqcut=cutoff value on G**2 for sphere inside fft box. idir=direction of the current perturbation ipert=type of the perturbation mpi_enreg=informations about MPI parallelization natom=number of atoms in cell. nfft=number of fft grid points (gsqcut=(boxcut**2)*ecut/(2._dp*(Pi**2)) ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft rhog(2,nfft)=array for Fourier transform of GS electron density rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
vhartr1(nfft)=Inhomogeneous term in strain-perturbation-induced Hartree potential in real space,
PARENTS
dfpt_nselt,dfpt_nstpaw,dfpt_rhotov
CHILDREN
fourdp,metric,ptabs_fourdp
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 62 subroutine hartrestr(gsqcut,idir,ipert,mpi_enreg,natom,nfft,ngfft,& 63 & paral_kgb,rhog,rprimd,vhartr1) 64 65 use defs_basis 66 use defs_abitypes 67 use m_errors 68 use m_profiling_abi 69 70 use m_mpinfo, only : ptabs_fourdp 71 use m_geometry, only : metric 72 73 !This section has been created automatically by the script Abilint (TD). 74 !Do not modify the following lines by hand. 75 #undef ABI_FUNC 76 #define ABI_FUNC 'hartrestr' 77 use interfaces_53_ffts 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: idir,ipert,natom,nfft,paral_kgb 85 real(dp),intent(in) :: gsqcut 86 type(MPI_type),intent(in) :: mpi_enreg 87 !arrays 88 integer,intent(in) :: ngfft(18) 89 real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3) 90 real(dp),intent(out) :: vhartr1(nfft) 91 92 !Local variables------------------------------- 93 !scalars 94 integer,parameter :: im=2,re=1 95 integer :: i1,i2,i23,i3,id2,id3,ig,ig2,ig3,ii,ii1,ing,istr,ka,kb,n1,n2,n3 96 real(dp),parameter :: tolfix=1.000000001_dp 97 real(dp) :: cutoff,ddends,den,dgsds,gqg2p3,gqgm12,gqgm13,gqgm23,gs,gs2,gs3 98 real(dp) :: term,ucvol 99 character(len=500) :: message 100 !arrays 101 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 102 integer :: id(3) 103 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 104 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 105 real(dp) :: dgmetds(3,3),gmet(3,3),gprimd(3,3),gqr(3),rmet(3,3) 106 real(dp),allocatable :: gq(:,:),work1(:,:) 107 108 ! ************************************************************************* 109 110 if( .not. (ipert==natom+3 .or. ipert==natom+4))then 111 write(message, '(a,i0,a,a)' )& 112 & 'From the calling routine, ipert=',ipert,ch10,& 113 & 'so this routine for the strain perturbation should not be called.' 114 MSG_BUG(message) 115 end if 116 117 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 118 119 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 120 121 !Get the distrib associated with this fft_grid 122 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 123 124 !Initialize a few quantities 125 cutoff=gsqcut*tolfix 126 127 istr=idir + 3*(ipert-natom-3) 128 129 if(istr<1 .or. istr>6)then 130 write(message, '(a,i10,a,a,a)' )& 131 & 'Input dir gives istr=',istr,' not allowed.',ch10,& 132 & 'Possible values are 1,2,3,4,5,6 only.' 133 MSG_BUG(message) 134 end if 135 136 ka=idx(2*istr-1);kb=idx(2*istr) 137 do ii = 1,3 138 dgmetds(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 139 end do 140 !For historical reasons: 141 dgmetds(:,:)=0.5_dp*dgmetds(:,:) 142 143 !In order to speed the routine, precompute the components of g+q 144 !Also check if the booked space was large enough... 145 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 146 do ii=1,3 147 id(ii)=ngfft(ii)/2+2 148 do ing=1,ngfft(ii) 149 ig=ing-(ing/id(ii))*ngfft(ii)-1 150 gq(ii,ing)=ig 151 end do 152 end do 153 154 ABI_ALLOCATE(work1,(2,nfft)) 155 id2=n2/2+2 156 id3=n3/2+2 157 !Triple loop on each dimension 158 do i3=1,n3 159 ig3=i3-(i3/id3)*n3-1 160 ! Precompute some products that do not depend on i2 and i1 161 gqr(3)=gq(3,i3) 162 gs3=gq(3,i3)*gq(3,i3)*gmet(3,3) 163 gqgm23=gq(3,i3)*gmet(2,3)*2 164 gqgm13=gq(3,i3)*gmet(1,3)*2 165 166 do i2=1,n2 167 if (fftn2_distrib(i2)==mpi_enreg%me_fft) then 168 gqr(2)=gq(2,i2) 169 gs2=gs3+ gq(2,i2)*(gq(2,i2)*gmet(2,2)+gqgm23) 170 gqgm12=gq(2,i2)*gmet(1,2)*2 171 gqg2p3=gqgm13+gqgm12 172 ig2=i2-(i2/id2)*n2-1 173 ! i23=n1*((i2-1)+n2*(i3-1)) 174 i23=n1*((ffti2_local(i2)-1)+(n2/mpi_enreg%nproc_fft)*(i3-1)) 175 ! Do the test that eliminates the Gamma point outside 176 ! of the inner loop 177 ii1=1 178 if(i23==0 .and. ig2==0 .and. ig3==0)then 179 ii1=2 180 work1(re,1+i23)=0.0_dp 181 work1(im,1+i23)=0.0_dp 182 end if 183 184 ! Final inner loop on the first dimension 185 ! (note the lower limit) 186 do i1=ii1,n1 187 gs=gs2+ gq(1,i1)*(gq(1,i1)*gmet(1,1)+gqg2p3) 188 ii=i1+i23 189 if(gs<=cutoff)then 190 den=piinv/gs 191 gqr(1)=gq(1,i1) 192 dgsds=& 193 & (gqr(1)*(dgmetds(1,1)*gqr(1)+dgmetds(1,2)*gqr(2)+dgmetds(1,3)*gqr(3))+ & 194 & gqr(2)*(dgmetds(2,1)*gqr(1)+dgmetds(2,2)*gqr(2)+dgmetds(2,3)*gqr(3))+ & 195 & gqr(3)*(dgmetds(3,1)*gqr(1)+dgmetds(3,2)*gqr(2)+dgmetds(3,3)*gqr(3)) ) 196 ddends=-piinv*dgsds/gs**2 197 if(istr<=3)then 198 term=2.0_dp*ddends-den 199 else 200 term=2.0_dp*ddends 201 end if 202 work1(re,ii)=rhog(re,ii)*term 203 work1(im,ii)=rhog(im,ii)*term 204 else 205 work1(re,ii)=0.0_dp 206 work1(im,ii)=0.0_dp 207 end if 208 209 end do ! End loop on i1 210 end if 211 end do ! End loop on i2 212 end do ! End loop on i3 213 214 ABI_DEALLOCATE(gq) 215 216 !Fourier Transform Vhartree. 217 !Vh in reciprocal space was stored in work1 218 call fourdp(1,work1,vhartr1,1,mpi_enreg,nfft,ngfft,paral_kgb,0) 219 220 ABI_DEALLOCATE(work1) 221 222 end subroutine hartrestr