TABLE OF CONTENTS


ABINIT/strhar [ Functions ]

[ Top ] [ 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