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!----------------------------------------------------------------------------
9SUBROUTINE deriv_drhoc( ngl, gl, omega, tpiba2, mesh, r, rab, rhoc, drhocg )
10  !--------------------------------------------------------------------------
11  !! Calculates the Fourier transform of \(d\text{Rho}_c/dG\).
12  !
13  USE kinds
14  USE constants,  ONLY : pi, fpi
15  !
16  IMPLICIT NONE
17  !
18  INTEGER :: ngl
19  !! input: the number of g shell
20  INTEGER :: mesh
21  !! input: the number of radial mesh points
22  REAL(DP), INTENT(IN) :: gl(ngl)
23  !! input: the number of G shells
24  REAL(DP), INTENT(IN) :: r(mesh)
25  !! input: the radial mesh
26  REAL(DP), INTENT(IN) :: rab(mesh)
27  !! input: the derivative of the radial mesh
28  REAL(DP), INTENT(IN) :: rhoc(mesh)
29  !! input: the radial core charge
30  REAL(DP), INTENT(IN) :: omega
31  !! input: the volume of the unit cell
32  REAL(DP), INTENT(IN) :: tpiba2
33  !! input: 2 times pi / alat
34  REAL(DP), INTENT(OUT) :: drhocg(ngl)
35  !! Fourier transform of d Rho_c/dG
36  !
37  ! ... local variables
38  !
39  REAL(DP) :: gx, rhocg1
40  ! the modulus of g for a given shell
41  ! the fourier transform
42  REAL(DP), ALLOCATABLE :: aux(:)
43  ! auxiliary memory for integration
44  INTEGER :: ir, igl, igl0
45  ! counter on radial mesh points
46  ! counter on g shells
47  ! lower limit for loop on ngl
48  !
49  ! G=0 term
50  !
51  IF (gl(1) < 1.0d-8) THEN
52     drhocg(1) = 0.0d0
53     igl0 = 2
54  ELSE
55     igl0 = 1
56  ENDIF
57  !
58  ! G <> 0 term
59  !
60!$omp parallel private(aux, gx, rhocg1)
61  !
62  ALLOCATE( aux(mesh) )
63!$omp do
64  DO igl = igl0, ngl
65     gx = SQRT( gl(igl) * tpiba2 )
66     DO ir = 1, mesh
67        aux(ir) = r(ir)*rhoc(ir)*( r(ir) * COS(gx*r(ir)) /       &
68                                      gx - SIN(gx*r(ir)) / gx**2 )
69     ENDDO
70     CALL simpson( mesh, aux, rab, rhocg1 )
71     drhocg(igl) = fpi / omega * rhocg1
72  ENDDO
73!$omp end do nowait
74  DEALLOCATE( aux )
75  !
76!$omp end parallel
77  !
78  RETURN
79  !
80END SUBROUTINE deriv_drhoc
81
82