1!
2! Copyright (C) 2001-2009 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 setlocq (xq, mesh, msh, rab, r, vloc_at, zp, tpiba2, ngm, &
11     g, omega, vloc)
12  !----------------------------------------------------------------------
13  !
14  !    This routine computes the Fourier transform of the local
15  !    part of the pseudopotential in the q+G vectors.
16  !
17  !    The local pseudopotential of the US case is always in
18  !    numerical form, expressed in Ry units.
19  !
20  USE kinds, only  : DP
21  USE constants, ONLY : e2, fpi, pi
22  USE Coul_cut_2D, ONLY: do_cutoff_2D
23  !
24  implicit none
25  !
26  !    first the dummy variables
27  !
28  integer :: ngm, mesh, msh
29  ! input: the number of G vectors
30  ! input: the dimensions of the mesh
31  ! input: mesh points for radial integration
32
33  real(DP) :: xq (3), zp, rab (mesh), r (mesh), vloc_at(mesh), tpiba2,&
34       omega, g (3, ngm), vloc (ngm)
35  ! input: the q point
36  ! input: valence pseudocharge
37  ! input: the derivative of mesh points
38  ! input: the mesh points
39  ! input: the pseudo on the radial
40  ! input: 2 pi / alat
41  ! input: the volume of the unit cell
42  ! input: the g vectors coordinates
43  ! output: the fourier transform of the potential
44  !
45  !    and the local variables
46  !
47  real(DP), parameter :: eps = 1.d-8
48  real(DP) :: vlcp, vloc0, fac, g2a, aux (mesh), &
49       aux1 (mesh), gx
50  ! auxiliary variables
51  ! gx = modulus of g vectors
52  real(DP), external :: qe_erf
53  ! the erf function
54  integer :: ig, ir
55  ! counters
56  !
57  ! Pseudopotentials in numerical form (Vnl(lloc) contain the local part)
58  ! in order to perform the Fourier transform, a term erf(r)/r is
59  ! subtracted in real space and added again in G space
60  !
61  ! first the G=0 term
62  !
63  !
64  IF (do_cutoff_2D) THEN
65     do ir = 1, msh
66        aux (ir) = r (ir) * (r (ir) * vloc_at (ir) + zp * e2    &
67                   * qe_erf (r (ir) ) )
68     enddo
69  ELSE
70      do ir = 1, msh
71         aux (ir) = r (ir) * (r (ir) * vloc_at (ir) + zp * e2)
72      enddo
73  ENDIF
74  !
75  call simpson (msh, aux, rab, vloc0)
76  !
77  !   here the G<>0 terms, we first compute the part of the integrand func
78  !   indipendent of |G| in real space
79  !
80  do ir = 1, msh
81     aux1 (ir) = r (ir) * vloc_at (ir) + zp * e2 * qe_erf (r (ir) )
82  enddo
83  fac = zp * e2 / tpiba2
84  !
85  !    and here we perform the integral, after multiplying for the |G|
86  !    dependent  part
87  !
88  do ig = 1, ngm
89     g2a = (xq (1) + g (1, ig) ) **2 + (xq (2) + g (2, ig) ) **2 + &
90          (xq (3) + g (3, ig) ) **2
91     if (g2a < eps) then
92        vloc (ig) = vloc0
93     else
94        gx = sqrt (g2a * tpiba2)
95        do ir = 1, msh
96           aux (ir) = aux1 (ir) * sin (gx * r (ir) ) / gx
97        enddo
98        call simpson (msh, aux, rab, vlcp)
99        !
100        !     here we add the analytic fourier transform of the erf function
101        !
102        !  if 2D cutoff calculation, do not re-add the FT of erf function
103        IF (.not. do_cutoff_2D) vlcp = vlcp - fac * exp ( - g2a * tpiba2 * 0.25d0) / g2a
104        vloc (ig) = vlcp
105     endif
106  enddo
107
108  vloc(:) = vloc(:) * fpi / omega
109
110  return
111end subroutine setlocq
112
113!----------------------------------------------------------------------
114subroutine setlocq_coul (xq, zp, tpiba2, ngm, g, omega, vloc)
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, ONLY: DP
121 USE constants, ONLY : fpi, e2, eps8
122 implicit none
123 !
124 integer, intent(in) :: ngm
125 real(DP) :: xq (3), zp, tpiba2, omega, g(3,ngm)
126 real(DP), intent (out) :: vloc(ngm)
127 !
128 real(DP) :: g2a
129 integer :: ig
130
131 do ig = 1, ngm
132  g2a = (xq (1) + g (1, ig) ) **2 + (xq (2) + g (2, ig) ) **2 + &
133        (xq (3) + g (3, ig) ) **2
134  if (g2a < eps8) then
135       vloc (ig) = 0.d0
136  else
137       vloc (ig) = - fpi * zp *e2 / omega / tpiba2 / g2a
138  endif
139 enddo
140
141end subroutine setlocq_coul
142