1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8
9      module m_rhofft
10
11      use precision
12
13      integer, parameter, public  :: FORWARD = -1
14      integer, parameter, public  :: BACKWARD = +1
15
16      public :: rhofft
17
18      CONTAINS
19
20      subroutine rhofft( CELL, N1, N2, N3, Mesh, nspin, RHO, rhog,
21     $                   direction)
22
23C *********************************************************************
24C WRITTEN BY J.M.SOLER, JUNE 1995.
25C **************** INPUT **********************************************
26C REAL*8  CELL(3,3)     : UNIT CELL VECTORS
27C INTEGER N1,N2,N3      : NUMBER OF LOCAL MESH DIVISIONS IN EACH CELL VECTOR
28C INTEGER Mesh(3)       : Number of global mesh divisions
29C REAL*4  RHO(N1,N2,N3) : DENSITIY AT MESH POINTS
30C *********************************************************************
31
32C     Modules
33      use precision,   only : dp, grid_p
34      use sys,         only : die
35      use m_fft,       only : fft     ! 3-D fast Fourier transform
36      use cellsubs,    only : reclat  ! Finds reciprocal lattice vectors
37      use cellsubs,    only : volcel  ! Finds unit cell volume
38
39      implicit          none
40
41C     Input/output variables
42      integer               :: N1, N2, N3, Mesh(3), nspin
43      integer, intent(in)   :: direction
44      real(grid_p)          :: RHO(N1*N2*N3,nspin)
45      real(dp)              :: CELL(3,3)
46      real(grid_p)          :: rhog(2,N1*N2*N3,nspin)
47
48C     Local variables
49      integer               :: ispin, ng
50      real(dp)              :: VOLUME, b(3,3)
51
52#ifdef DEBUG
53      call write_debug( '    PRE rhofft' )
54#endif
55
56C     Start time counter
57      call timer( 'rhofft', 1 )
58
59C     Find unit cell volume
60      VOLUME = VOLCEL( CELL )
61
62C     Find reciprocal lattice vectors
63      call reclat(CELL, B, 1 )
64
65      NG = mesh(1) * mesh(2) * mesh(3)
66
67      if (direction == FORWARD) then
68         rhog(1,:,:) = RHO(:,:)
69         rhog(2,:,:) = 0.0_grid_p
70
71         do ispin = 1, nspin
72            call fft( rhog(1,1,ispin), Mesh, -1 )
73         enddo
74         ! Units: electrons per unit cell
75         rhog(:,:,:) = rhog(:,:,:) * volume / dble(ng)
76      else if (direction == BACKWARD) then
77         do ispin = 1, nspin
78            call fft( rhog(1,1,ispin), Mesh, +1 )
79         enddo
80         rho(:,:) = rhog(1,:,:) * ng / volume
81         !  print *, "sum(rho)*dvol", sum(rho)*volume/dble(ng)
82      else
83         call die("wrong fft direction")
84      endif
85
86C     Stop time counter
87      call timer( 'rhofft', 2 )
88
89#ifdef DEBUG
90      call write_debug( '    POS rhofft' )
91#endif
92
93      end subroutine rhofft
94
95      end module m_rhofft
96