1!
2! Copyright (C) 2001-2007 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!----------------------------------------------------------------------
10SUBROUTINE stres_har( sigmahar )
11  !--------------------------------------------------------------------
12  !! Calculates the Hartree contribution to the stress
13  !
14  USE kinds,           ONLY: DP
15  USE constants,       ONLY: e2, fpi
16  USE cell_base,       ONLY: omega, tpiba2
17  USE ener,            ONLY: ehart
18  USE fft_base,        ONLY: dfftp
19  USE fft_interfaces,  ONLY: fwfft
20  USE gvect,           ONLY: ngm, gstart, g, gg
21  USE scf,             ONLY: rho
22  USE control_flags,   ONLY: gamma_only
23  USE wavefunctions,   ONLY: psic
24  USE mp_bands,        ONLY: intra_bgrp_comm
25  USE mp,              ONLY: mp_sum
26  USE Coul_cut_2D,     ONLY: do_cutoff_2D, cutoff_stres_sigmahar
27  !
28  IMPLICIT NONE
29  !
30  REAL(DP) :: sigmahar(3,3)
31  !! Hartree term of the stress tensor
32  !
33  ! ... local variables
34  !
35  REAL(DP) :: shart, g2
36  REAL(DP), PARAMETER :: eps = 1.d-8
37  INTEGER :: ig, l, m
38  !
39  sigmahar(:,:) = 0.0_DP
40  psic(:) = CMPLX( rho%of_r(:,1), KIND=DP )
41  !
42  CALL fwfft('Rho', psic, dfftp)
43  ! psic contains now the charge density in G space
44  ! the  G=0 component is not computed
45  IF (do_cutoff_2D) THEN
46     CALL cutoff_stres_sigmahar( psic, sigmahar )
47  ELSE
48     DO ig = gstart, ngm
49        g2 = gg(ig) * tpiba2
50        shart = psic(dfftp%nl(ig)) * CONJG(psic(dfftp%nl(ig))) / g2
51        DO l = 1, 3
52           DO m = 1, l
53              sigmahar(l,m) = sigmahar(l,m) + shart * tpiba2 * 2 * &
54                              g(l,ig) * g(m,ig) / g2
55           ENDDO
56        ENDDO
57     ENDDO
58  ENDIF
59  !
60  CALL mp_sum( sigmahar, intra_bgrp_comm )
61  !
62  IF (gamma_only) THEN
63     sigmahar(:,:) = fpi * e2 * sigmahar(:,:)
64  ELSE
65     sigmahar(:,:) = fpi * e2 * sigmahar(:,:) * 0.5_DP
66  ENDIF
67  !
68  DO l = 1, 3
69     sigmahar(l,l) = sigmahar(l,l) - ehart / omega
70  ENDDO
71  !
72  DO l = 1, 3
73     DO m = 1, l-1
74        sigmahar(m,l) = sigmahar(l,m)
75     ENDDO
76  ENDDO
77  !
78  sigmahar(:,:) = -sigmahar(:,:)
79  !
80  !
81  RETURN
82  !
83END SUBROUTINE stres_har
84
85