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 dvloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, gl, &
11     omega, dvloc)
12  !----------------------------------------------------------------------
13  !
14  ! dvloc = D Vloc (g^2) / D g^2 = (1/2g) * D Vloc(g) / D g
15  !
16  USE kinds
17  USE constants , ONLY : pi, fpi, e2, eps8
18  USE Coul_cut_2D, ONLY: do_cutoff_2D
19  USE esm, ONLY : do_comp_esm, esm_bc
20  implicit none
21  !
22  !    first the dummy variables
23  !
24  integer, intent(in) :: ngl, mesh, msh
25  ! the number of shell of G vectors
26  ! max number of mesh points
27  ! number of mesh points for radial integration
28
29  real(DP), intent(in) :: zp, rab (mesh), r (mesh), vloc_at (mesh), &
30                          tpiba2, omega, gl (ngl)
31  ! valence pseudocharge
32  ! the derivative of the radial grid
33  ! the radial grid
34  ! the pseudo on the radial grid
35  ! 2 pi / alat
36  ! the volume of the unit cell
37  ! the moduli of g vectors for each s
38  !
39  real(DP), intent(out) ::  dvloc (ngl)
40  ! the fourier transform dVloc/dG
41  !
42  real(DP) :: vlcp, g2a, gx
43  real(DP), allocatable ::  aux (:), aux1 (:)
44  real(DP), external ::  qe_erf
45
46  integer :: i, igl, igl0
47  ! counter on erf functions or gaussians
48  ! counter on g shells vectors
49  ! first shell with g != 0
50
51  ! the  G=0 component is not computed
52  if (gl (1) < eps8) then
53     dvloc (1) = 0.0d0
54     igl0 = 2
55  else
56     igl0 = 1
57  endif
58
59  ! Pseudopotentials in numerical form (Vloc contains the local part)
60  ! In order to perform the Fourier transform, a term erf(r)/r is
61  ! subtracted in real space and added again in G space
62
63  allocate (aux1( mesh))
64  !
65  !   This is the part of the integrand function
66  !   indipendent of |G| in real space
67  !
68  do i = 1, msh
69     aux1 (i) = r (i) * vloc_at (i) + zp * e2 * qe_erf (r (i) )
70  enddo
71  !
72!$omp parallel private(aux, gx, vlcp, g2a)
73  !
74  allocate (aux( mesh))
75  !
76!$omp do
77  do igl = igl0, ngl
78     gx = sqrt (gl (igl) * tpiba2)
79     !
80     !    and here we perform the integral, after multiplying for the |G|
81     !    dependent  part
82     !
83     ! DV(g)/Dg = Integral of r (Dj_0(gr)/Dg) V(r) dr
84     do i = 1, msh
85        aux (i) = aux1 (i) * (r (i) * cos (gx * r (i) ) / gx - sin (gx &
86             * r (i) ) / gx**2)
87     enddo
88     call simpson (msh, aux, rab, vlcp)
89     ! DV(g^2)/Dg^2 = (DV(g)/Dg)/2g
90     vlcp = fpi / omega / 2.0d0 / gx * vlcp
91
92    ! for ESM stress
93     ! In ESM, vloc and dvloc have only short term.
94     IF ( ( .not. do_comp_esm ) .or. ( esm_bc .eq. 'pbc' ) ) THEN
95        ! subtract the long-range term
96        IF (.not.do_cutoff_2D) then ! 2D cutoff: do not re-add LR part here (re-added later in stres_loc)
97         g2a = gl (igl) * tpiba2 / 4.d0
98         vlcp = vlcp + fpi / omega * zp * e2 * exp ( - g2a) * (g2a + &
99             1.d0) / (gl (igl) * tpiba2) **2
100        ENDIF
101     END IF
102     dvloc (igl) = vlcp
103  enddo
104!$omp end do nowait
105  deallocate (aux)
106!$omp end parallel
107  !
108  deallocate (aux1)
109
110  return
111end subroutine dvloc_of_g
112!
113!----------------------------------------------------------------------
114subroutine dvloc_coul (zp, tpiba2, ngl, gl, omega, dvloc)
115  !----------------------------------------------------------------------
116  !
117  !    Fourier transform of the Coulomb potential - For all-electron
118  !    calculations, in specific cases only, for testing purposes
119  !
120  USE kinds
121  USE constants , ONLY : fpi, e2, eps8
122  implicit none
123  !
124  integer, intent(in) :: ngl
125  ! the number of shell of G vectors
126  real(DP), intent(in) :: zp, tpiba2, omega, gl (ngl)
127  ! valence pseudocharge
128  ! 2 pi / alat
129  ! the volume of the unit cell
130  ! the moduli of g vectors for each s
131  real(DP), intent(out) :: dvloc (ngl)
132  ! fourier transform: dvloc = D Vloc (g^2) / D g^2 = 4pi e^2/omegai /G^4
133  !
134  integer :: igl0
135  ! first shell with g != 0
136
137  ! the  G=0 component is 0
138  if (gl (1) < eps8) then
139     dvloc (1) = 0.0d0
140     igl0 = 2
141  else
142     igl0 = 1
143  endif
144
145  dvloc  (igl0:ngl) = fpi * zp * e2 / omega / ( tpiba2 * gl (igl0:ngl) ) ** 2
146
147return
148end subroutine dvloc_coul
149
150