TABLE OF CONTENTS
ABINIT/dfpt_eltfrhar [ Functions ]
NAME
dfpt_eltfrhar
FUNCTION
Compute the frozen-wavefunction hartree enegy contribution to the elastic tensor
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DRH, DCA, XG, GM, AR) 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.
INPUTS
rprimd(3,3)=dimensional primitive translation vectors (bohr) gsqcut =Fourier cutoff on G^2 for "large sphere" of radius double that of the basis sphere--appropriate for charge density rho(G), Hartree potential, and pseudopotentials mpi_enreg=informations 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 rhog(2,nfft)=total electron density in G space
OUTPUT
eltfrhar(6,6)=non-symmetrized kinetic energy contribution to the elastic tensor
NOTES
*based largely on hartre.f
PARENTS
respfn
CHILDREN
metric,ptabs_fourdp,timab,xmpi_sum
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 49 subroutine dfpt_eltfrhar(eltfrhar,rprimd,gsqcut,mpi_enreg,nfft,ngfft,rhog) 50 51 use defs_basis 52 use defs_abitypes 53 use m_profiling_abi 54 use m_errors 55 use m_xmpi 56 57 use m_geometry, only : metric 58 use m_mpinfo, only : ptabs_fourdp 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'dfpt_eltfrhar' 64 use interfaces_18_timing 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------------ 70 !scalars 71 integer,intent(in) :: nfft 72 real(dp),intent(in) :: gsqcut 73 type(MPI_type),intent(in) :: mpi_enreg 74 !arrays 75 integer,intent(in) :: ngfft(18) 76 real(dp),intent(in) :: rhog(2,nfft),rprimd(3,3) 77 real(dp),intent(out) :: eltfrhar(6,6) 78 79 !Local variables------------------------------- 80 !scalars 81 integer,parameter :: im=2,re=1 82 integer :: i1,i2,i23,i3,id2,id3,ierr,ig,ig2,ig3,ii,ii1,ing,istr1,istr2,jj 83 integer :: ka,kb,kd,kg,me_fft,n1,n2,n3,nproc_fft 84 real(dp),parameter :: tolfix=1.000000001_dp 85 real(dp) :: cutoff,d2eacc,d2etot,d2gs,deacc01,deacc10,dgs01,dgs10,eacc,fact,gs 86 real(dp) :: term,ucvol 87 !arrays 88 integer,save :: idx(12)=(/1,1,2,2,3,3,3,2,3,1,2,1/) 89 integer :: id(3) 90 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 91 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 92 real(dp) :: d2gm(3,3),dgm01(3,3),dgm10(3,3),gmet(3,3),gprimd(3,3),gqr(3) 93 real(dp) :: rmet(3,3),tsec(2) 94 real(dp),allocatable :: gq(:,:) 95 96 ! ************************************************************************* 97 98 !Compute gmet, gprimd and ucvol from rprimd 99 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 100 101 eltfrhar(:,:)=0.0_dp 102 103 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 104 me_fft=ngfft(11) 105 nproc_fft=ngfft(10) 106 107 !Get the distrib associated with this fft_grid 108 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 109 110 !Initialize a few quantities 111 fact=0.5_dp*ucvol/pi 112 cutoff=gsqcut*tolfix 113 114 !In order to speed the routine, precompute the components of g+q 115 !Also check if the booked space was large enough... 116 ABI_ALLOCATE(gq,(3,max(n1,n2,n3))) 117 do ii=1,3 118 id(ii)=ngfft(ii)/2+2 119 do ing=1,ngfft(ii) 120 ig=ing-(ing/id(ii))*ngfft(ii)-1 121 gq(ii,ing)=ig 122 end do 123 end do 124 125 !Loop over 2nd strain index 126 do istr2=1,6 127 ! Loop over 1st strain index, upper triangle only 128 do istr1=1,istr2 129 130 ka=idx(2*istr1-1);kb=idx(2*istr1);kg=idx(2*istr2-1);kd=idx(2*istr2) 131 132 do ii = 1,3 133 dgm01(:,ii)=-(gprimd(ka,:)*gprimd(kb,ii)+gprimd(kb,:)*gprimd(ka,ii)) 134 dgm10(:,ii)=-(gprimd(kg,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kg,ii)) 135 end do 136 137 d2gm(:,:)=0._dp 138 do ii = 1,3 139 if(ka==kg) d2gm(:,ii)=d2gm(:,ii)& 140 & +gprimd(kb,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(kb,ii) 141 if(ka==kd) d2gm(:,ii)=d2gm(:,ii)& 142 & +gprimd(kb,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(kb,ii) 143 if(kb==kg) d2gm(:,ii)=d2gm(:,ii)& 144 & +gprimd(ka,:)*gprimd(kd,ii)+gprimd(kd,:)*gprimd(ka,ii) 145 if(kb==kd) d2gm(:,ii)=d2gm(:,ii)& 146 & +gprimd(ka,:)*gprimd(kg,ii)+gprimd(kg,:)*gprimd(ka,ii) 147 end do 148 d2gm(:,:)=0.5_dp*d2gm(:,:) 149 150 ! initialize energy accumulator 151 eacc=0._dp 152 deacc01=0._dp 153 deacc10=0._dp 154 d2eacc=0._dp 155 156 id2=n2/2+2 157 id3=n3/2+2 158 ! Triple loop on each dimension 159 do i3=1,n3 160 ig3=i3-(i3/id3)*n3-1 161 gqr(3)=gq(3,i3) 162 do i2=1,n2 163 if (fftn2_distrib(i2)==me_fft) then 164 gqr(2)=gq(2,i2) 165 ig2=i2-(i2/id2)*n2-1 166 i23=n1*(ffti2_local(i2)-1 +(n2/nproc_fft)*(i3-1)) 167 ! Do the test that eliminates the Gamma point outside 168 ! of the inner loop 169 ii1=1 170 if(i23==0 .and. ig2==0 .and. ig3==0)then 171 ii1=2 172 end if 173 174 ! Final inner loop on the first dimension 175 ! (note the lower limit) 176 do i1=ii1,n1 177 gqr(1)=gq(1,i1) 178 gs=(gmet(1,1)*gqr(1)*gqr(1)+gmet(2,2)*gqr(2)*gqr(2)+& 179 & gmet(3,3)*gqr(3)*gqr(3)+2._dp*& 180 & (gmet(1,2)*gqr(1)*gqr(2) + gmet(1,3)*gqr(1)*gqr(3)+& 181 & gmet(2,3)*gqr(2)*gqr(3)) ) 182 ii=i1+i23 183 if(gs<=cutoff)then 184 dgs01=(dgm01(1,1)*gqr(1)*gqr(1)+dgm01(2,2)*gqr(2)*gqr(2)+& 185 & dgm01(3,3)*gqr(3)*gqr(3)+2._dp*& 186 & (dgm01(1,2)*gqr(1)*gqr(2) + dgm01(1,3)*gqr(1)*gqr(3)+& 187 & dgm01(2,3)*gqr(2)*gqr(3)) ) 188 dgs10=(dgm10(1,1)*gqr(1)*gqr(1)+dgm10(2,2)*gqr(2)*gqr(2)+& 189 & dgm10(3,3)*gqr(3)*gqr(3)+2._dp*& 190 & (dgm10(1,2)*gqr(1)*gqr(2) + dgm10(1,3)*gqr(1)*gqr(3)+& 191 & dgm10(2,3)*gqr(2)*gqr(3)) ) 192 d2gs =(d2gm(1,1)*gqr(1)*gqr(1)+d2gm(2,2)*gqr(2)*gqr(2)+& 193 & d2gm(3,3)*gqr(3)*gqr(3)+2._dp*& 194 & (d2gm(1,2)*gqr(1)*gqr(2) + d2gm(1,3)*gqr(1)*gqr(3)+& 195 & d2gm(2,3)*gqr(2)*gqr(3)) ) 196 197 term=(rhog(re,ii)**2+rhog(im,ii)**2)/gs 198 eacc=eacc+term 199 deacc01=deacc01+dgs01*term/gs 200 deacc10=deacc10+dgs10*term/gs 201 d2eacc=d2eacc+(-d2gs+2._dp*dgs01*dgs10/gs)*term/gs 202 end if 203 204 ! End loop on i1 205 end do 206 end if 207 ! End loop on i2 208 end do 209 ! End loop on i3 210 end do 211 212 ! Add contributions taking account diagonal strain terms (from ucvol 213 ! derivatives) 214 d2etot=d2eacc 215 if(istr1<=3) d2etot=d2etot+deacc10 216 if(istr2<=3) d2etot=d2etot+deacc01 217 if(istr1<=3 .and. istr2<=3) d2etot=d2etot+eacc 218 219 eltfrhar(istr1,istr2)=fact*d2etot 220 221 ! End loop on istr1 222 end do 223 ! End loop in istr2 224 end do 225 226 ABI_DEALLOCATE(gq) 227 228 !Init mpi_comm 229 call timab(48,1,tsec) 230 call xmpi_sum(eltfrhar,mpi_enreg%comm_fft,ierr) 231 call timab(48,2,tsec) 232 233 !Fill in lower triangle 234 do jj=2,6 235 do ii=1,jj-1 236 eltfrhar(jj,ii)=eltfrhar(ii,jj) 237 end do 238 end do 239 end subroutine dfpt_eltfrhar