2003-01-20 05:58:50 +08:00
|
|
|
!
|
|
|
|
! Copyright (C) 2001 PWSCF group
|
|
|
|
! This file is distributed under the terms of the
|
|
|
|
! GNU General Public License. See the file `License'
|
|
|
|
! in the root directory of the present distribution,
|
|
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
|
|
!
|
|
|
|
!
|
|
|
|
!----------------------------------------------------------------------
|
|
|
|
|
2003-02-08 00:04:36 +08:00
|
|
|
subroutine stres_har (sigmahar)
|
2003-01-20 05:58:50 +08:00
|
|
|
!----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
#include "machine.h"
|
2003-02-08 00:04:36 +08:00
|
|
|
use pwcom
|
2003-10-03 22:01:11 +08:00
|
|
|
USE wavefunctions, ONLY : psic
|
2003-01-20 05:58:50 +08:00
|
|
|
implicit none
|
|
|
|
!
|
|
|
|
real(kind=DP) :: sigmahar (3, 3), shart, g2
|
|
|
|
real(kind=DP), parameter :: eps = 1.d-8
|
2003-02-08 00:04:36 +08:00
|
|
|
integer :: is, ig, igl0, l, m
|
2003-01-20 05:58:50 +08:00
|
|
|
|
|
|
|
sigmahar(:,:) = 0.d0
|
|
|
|
psic (:) = (0.d0, 0.d0)
|
2003-02-08 00:04:36 +08:00
|
|
|
do is = 1, nspin
|
|
|
|
call DAXPY (nrxx, 1.d0, rho (1, is), 1, psic, 2)
|
2003-01-20 05:58:50 +08:00
|
|
|
enddo
|
|
|
|
|
2003-02-08 00:04:36 +08:00
|
|
|
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
|
2003-01-20 05:58:50 +08:00
|
|
|
! psic contains now the charge density in G space
|
|
|
|
! the G=0 component is not computed
|
2003-02-08 00:04:36 +08:00
|
|
|
do ig = gstart, ngm
|
|
|
|
g2 = gg (ig) * tpiba2
|
|
|
|
shart = psic (nl (ig) ) * conjg (psic (nl (ig) ) ) / g2
|
|
|
|
do l = 1, 3
|
|
|
|
do m = 1, l
|
2003-01-20 05:58:50 +08:00
|
|
|
sigmahar (l, m) = sigmahar (l, m) + shart * tpiba2 * 2 * &
|
|
|
|
g (l, ig) * g (m, ig) / g2
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
enddo
|
2003-02-21 22:57:00 +08:00
|
|
|
#ifdef __PARA
|
2003-02-08 00:04:36 +08:00
|
|
|
call reduce (9, sigmahar)
|
2003-01-20 05:58:50 +08:00
|
|
|
#endif
|
|
|
|
if (gamma_only) then
|
|
|
|
sigmahar(:,:) = fpi * e2 * sigmahar(:,:)
|
|
|
|
else
|
|
|
|
sigmahar(:,:) = 0.5 * fpi * e2 * sigmahar(:,:)
|
|
|
|
end if
|
2003-02-08 00:04:36 +08:00
|
|
|
do l = 1, 3
|
|
|
|
sigmahar (l, l) = sigmahar (l, l) - ehart / omega
|
2003-01-20 05:58:50 +08:00
|
|
|
enddo
|
2003-02-08 00:04:36 +08:00
|
|
|
do l = 1, 3
|
|
|
|
do m = 1, l - 1
|
|
|
|
sigmahar (m, l) = sigmahar (l, m)
|
2003-01-20 05:58:50 +08:00
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
|
|
|
|
sigmahar(:,:) = -sigmahar(:,:)
|
2003-02-08 00:04:36 +08:00
|
|
|
|
|
|
|
return
|
2003-01-20 05:58:50 +08:00
|
|
|
end subroutine stres_har
|
|
|
|
|