TABLE OF CONTENTS
ABINIT/strhar [ Functions ]
NAME
strhar
FUNCTION
Compute Hartree energy contribution to stress tensor (Cartesian coordinates).
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 .
INPUTS
ehart=Hartree energy (hartree) gsqcut=cutoff value on $G^2$ for (large) sphere inside fft box. $gsqcut=(boxcut^2)*ecut/(2._dp*(\pi^2))$ 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)=Fourier transform of charge density (bohr^-3) rhog(2,nfft)= optional argument: Fourier transform of a second charge density (bohr^-3) rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
harstr(6)=components of Hartree part of stress tensor (Cartesian coordinates, symmetric tensor) in hartree/bohr^3 Definition of symmetric tensor storage: store 6 unique components in the order 11, 22, 33, 32, 31, 21 (suggested by Xavier Gonze).
PARENTS
stress
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 strhar(ehart,gsqcut,harstr,mpi_enreg,nfft,ngfft,rhog,rprimd,& 50 & rhog2) ! optional argument 51 52 use defs_basis 53 use defs_abitypes 54 use m_errors 55 use m_profiling_abi 56 use m_xmpi 57 58 use m_geometry, only : metric 59 use m_mpinfo, only : ptabs_fourdp 60 61 !This section has been created automatically by the script Abilint (TD). 62 !Do not modify the following lines by hand. 63 #undef ABI_FUNC 64 #define ABI_FUNC 'strhar' 65 use interfaces_18_timing 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 integer,intent(in) :: nfft 73 real(dp),intent(in) :: ehart,gsqcut 74 type(MPI_type),intent(in) :: mpi_enreg 75 !arrays 76 integer,intent(in) :: ngfft(18) 77 real(dp),intent(in) :: rprimd(3,3),rhog(2,nfft) 78 real(dp),intent(in),optional :: rhog2(2,nfft) 79 real(dp),intent(out) :: harstr(6) 80 81 !Local variables------------------------------- 82 !scalars 83 integer,parameter :: im=2,re=1 84 integer :: i1,i2,i3,id1,id2,id3,ierr,ig1,ig2,ig3,ii,irho2,me_fft,n1,n2,n3,nproc_fft 85 real(dp) :: cutoff,gsquar,rhogsq,tolfix=1.000000001_dp,ucvol 86 !arrays 87 real(dp) :: gcart(3),gmet(3,3),gprimd(3,3),rmet(3,3),tsec(2) 88 integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:) 89 integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:) 90 91 ! ************************************************************************* 92 93 call timab(568,1,tsec) 94 95 harstr(:)=zero 96 !ehtest=0.0_dp (used for testing) 97 98 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 99 100 irho2=0;if (present(rhog2)) irho2=1 101 102 !Conduct looping over all fft grid points to find G vecs inside gsqcut 103 !Include G**2 on surface of cutoff sphere as well as inside: 104 cutoff=gsqcut*tolfix 105 n1=ngfft(1) ; n2=ngfft(2) ; n3=ngfft(3) 106 me_fft=ngfft(11) 107 nproc_fft=ngfft(10) 108 id1=n1/2+2 109 id2=n2/2+2 110 id3=n3/2+2 111 ii=0 112 113 ! Get the distrib associated with this fft_grid 114 call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local) 115 116 do i3=1,n3 117 ig3=i3-(i3/id3)*n3-1 118 do i2=1,n2 119 ig2=i2-(i2/id2)*n2-1 120 if (fftn2_distrib(i2)==me_fft) then 121 do i1=1,n1 122 ig1=i1-(i1/id1)*n1-1 123 ! ii=ii+1 124 ii=i1+n1*(ffti2_local(i2)-1+(n2/nproc_fft)*(i3-1)) 125 ! ** GET RID OF THIS IF STATEMENT LATER for speed if needed 126 ! Avoid G=0: 127 ! if (ii>1) then 128 if (ig1==0 .and. ig2==0 .and. ig3==0) cycle 129 ! Compute cartesian components of G 130 gcart(1)=gprimd(1,1)*dble(ig1)+gprimd(1,2)*dble(ig2)+gprimd(1,3)*dble(ig3) 131 gcart(2)=gprimd(2,1)*dble(ig1)+gprimd(2,2)*dble(ig2)+gprimd(2,3)*dble(ig3) 132 gcart(3)=gprimd(3,1)*dble(ig1)+gprimd(3,2)*dble(ig2)+gprimd(3,3)*dble(ig3) 133 ! Compute |G|^2 134 gsquar=gcart(1)**2+gcart(2)**2+gcart(3)**2 135 136 ! Keep only G**2 inside larger cutoff (not sure this is needed): 137 if (gsquar<=cutoff) then 138 ! take |rho(G)|^2 for complex rhog 139 if (irho2==0) then 140 rhogsq=rhog(re,ii)**2+rhog(im,ii)**2 141 else 142 rhogsq=rhog(re,ii)*rhog2(re,ii)+rhog(im,ii)*rhog2(im,ii) 143 end if 144 harstr(1)=harstr(1)+(rhogsq/gsquar**2)*gcart(1)*gcart(1) 145 harstr(2)=harstr(2)+(rhogsq/gsquar**2)*gcart(2)*gcart(2) 146 harstr(3)=harstr(3)+(rhogsq/gsquar**2)*gcart(3)*gcart(3) 147 harstr(4)=harstr(4)+(rhogsq/gsquar**2)*gcart(3)*gcart(2) 148 harstr(5)=harstr(5)+(rhogsq/gsquar**2)*gcart(3)*gcart(1) 149 harstr(6)=harstr(6)+(rhogsq/gsquar**2)*gcart(2)*gcart(1) 150 end if 151 ! end if 152 end do 153 end if 154 end do 155 end do 156 157 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 158 #ifdef FC_IBM 159 write(std_out,*)' strhar : before mpi_comm, harstr=',harstr 160 #endif 161 162 !Init mpi_comm 163 if(mpi_enreg%nproc_fft>1)then 164 call timab(48,1,tsec) 165 call xmpi_sum(harstr,mpi_enreg%comm_fft ,ierr) 166 call timab(48,2,tsec) 167 end if 168 169 #ifdef FC_IBM 170 !DO not remove : seems needed to avoid problem with pathscale compiler, in parallel 171 write(std_out,*)' strhar : after mpi_comm, harstr=',harstr 172 write(std_out,*)' strhar : ehart,ucvol=',ehart,ucvol 173 #endif 174 175 !Normalize and add term -ehart/ucvol on diagonal 176 harstr(1)=harstr(1)/pi-ehart/ucvol 177 harstr(2)=harstr(2)/pi-ehart/ucvol 178 harstr(3)=harstr(3)/pi-ehart/ucvol 179 harstr(4)=harstr(4)/pi 180 harstr(5)=harstr(5)/pi 181 harstr(6)=harstr(6)/pi 182 183 call timab(568,2,tsec) 184 185 end subroutine strhar