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