Minor cleanup; integration routine prepared for dealing with even

number of grid point (still commented out). I think we should figure 
out which integration routine is the best and stick to it: there are
two simpson-style routines that yield slightly different results


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@9416 c92efa57-630b-4861-b058-cf58834340f0

flib/simpsn.f90

@@ -1,16 +1,16 @@
 !
-! Copyright (C) 2001 PWSCF group
+! Copyright (C) 2001-2012 Quantum ESPRESSO group
 ! This file is distributed under the terms of the
 ! GNU General Public License. See the file `License'
 ! in the root directory of the present distribution,
 ! or http://www.gnu.org/copyleft/gpl.txt .
 !
 !-----------------------------------------------------------------------
-subroutine simpson (mesh, func, rab, asum)
+SUBROUTINE simpson(mesh, func, rab, asum)
   !-----------------------------------------------------------------------
   !
   !     simpson's rule integration. On input:
-  !       mesh = mhe number of grid points (should be odd)
+  !       mesh = the number of grid points (should be odd)
   !       func(i)= function to be integrated
   !       rab(i) = r(i) * dr(i)/di * di
   !     For the logarithmic grid not including r=0 :
@@ -20,37 +20,41 @@ subroutine simpson (mesh, func, rab, asum)
   !     Output in asum = \sum_i c_i f(i)*rab(i) = \int_0^\infty f(r) dr
   !     where c_i are alternativaly 2/3, 4/3 except c_1 = c_mesh = 1/3
   !
-  use kinds, ONLY: DP
+  USE kinds, ONLY: DP
-  implicit none
+  IMPLICIT NONE
-  integer, intent(in) :: mesh
+  INTEGER, INTENT(in) :: mesh
-  real(DP), intent(in) :: rab (mesh), func (mesh)
+  real(DP), INTENT(in) :: rab (mesh), func (mesh)
-  real(DP), intent(out):: asum
+  real(DP), INTENT(out):: asum
   !
   real(DP) :: f1, f2, f3, r12
-  integer :: i
+  INTEGER :: i
   !
-  !     routine assumes that mesh is an odd number so run check
-  !     if ( mesh+1 - ( (mesh+1) / 2 ) * 2 .ne. 1 ) then
-  !       write(*,*) '***error in subroutine radlg'
-  !       write(*,*) 'routine assumes mesh is odd but mesh =',mesh+1
-  !       stop
-  !     endif
   asum = 0.0d0
-  r12 = 1.0d0 / 12.0d0
+  r12 = 1.0d0 / 3.0d0
   f3 = func (1) * rab (1) * r12
 
-  do i = 2, mesh - 1, 2
+  DO i = 2, mesh - 1, 2
      f1 = f3
      f2 = func (i) * rab (i) * r12
      f3 = func (i + 1) * rab (i + 1) * r12
-     asum = asum + 4.0d0 * f1 + 16.0d0 * f2 + 4.0d0 * f3
+     asum = asum + f1 + 4.0d0 * f2 + f3
-  enddo
+  ENDDO
+  !
+  ! if mesh is not odd, use open formula instead:
+  ! ... 2/3*f(n-5) + 4/3*f(n-4) + 13/12*f(n-3) + 0*f(n-2) + 27/12*f(n-1)
+  !!! Under testing
+  !
+  !IF ( MOD(mesh,2) == 0 ) THEN
+  !   print *, 'mesh even: correction:', f1*5.d0/4.d0-4.d0*f2+23.d0*f3/4.d0, &
+  !                                      func(mesh)*rab(mesh), asum
+  !   asum = asum + f1*5.d0/4.d0 - 4.d0*f2 + 23.d0*f3/4.d0
+  !END IF
 
-  return
+  RETURN
-end subroutine simpson
+END SUBROUTINE simpson
 
 !=-----------------------------------------------------------------------
-subroutine simpson_cp90( mesh, func, rab, asum )
+SUBROUTINE simpson_cp90( mesh, func, rab, asum )
   !-----------------------------------------------------------------------
   !
   !    This routine computes the integral of a function defined on a
@@ -65,16 +69,16 @@ subroutine simpson_cp90( mesh, func, rab, asum )
   !
   !    last revised 12 May 1995 by Andrea Dal Corso
   !
-  use kinds, ONLY: DP
+  USE kinds, ONLY: DP
-  implicit none
+  IMPLICIT NONE
-  integer, intent(in) :: mesh
+  INTEGER, INTENT(in) :: mesh
-  real(DP), intent(in) :: rab (mesh), func (mesh)
+  real(DP), INTENT(in) :: rab (mesh), func (mesh)
-  real(DP), intent(out):: asum
+  real(DP), INTENT(out):: asum
   !
   real(DP) :: c(4)
-  integer ::i
+  INTEGER ::i
   !
-  if ( mesh < 8 ) call errore ('simpson_cp90','few mesh points',8)
+  IF ( mesh < 8 ) CALL errore ('simpson_cp90','few mesh points',8)
 
   c(1) = 109.0d0 / 48.d0
   c(2) = -5.d0 / 48.d0
@@ -85,12 +89,12 @@ subroutine simpson_cp90( mesh, func, rab, asum )
        + ( func(2)*rab(2) + func(mesh-1)*rab(mesh-1) )*c(2) &
        + ( func(3)*rab(3) + func(mesh-2)*rab(mesh-2) )*c(3) &
        + ( func(4)*rab(4) + func(mesh-3)*rab(mesh-3) )*c(4)
-  do i=5,mesh-4
+  DO i=5,mesh-4
      asum = asum + func(i)*rab(i)
-  end do
+  ENDDO
 
-  return
+  RETURN
-end subroutine simpson_cp90
+END SUBROUTINE simpson_cp90
 !
 !-----------------------------------------------------------------------
 SUBROUTINE herman_skillman_int(mesh,func,rab,asum)
@@ -98,12 +102,12 @@ SUBROUTINE herman_skillman_int(mesh,func,rab,asum)
   !     simpson rule integration for herman skillman mesh (obsolescent)
   !     Input as in "simpson". BEWARE: "func" is overwritten!!!
   !
-  use kinds, ONLY: DP
+  USE kinds, ONLY: DP
   IMPLICIT NONE
-  integer, intent(in) :: mesh
+  INTEGER, INTENT(in) :: mesh
-  real(DP), intent(in) :: rab (mesh)
+  real(DP), INTENT(in) :: rab (mesh)
-  real(DP), intent(inout) :: func (mesh)
+  real(DP), INTENT(inout) :: func (mesh)
-  real(DP), intent(out):: asum
+  real(DP), INTENT(out):: asum
   !
   INTEGER :: i, j, k, i1, nblock
   REAL(DP) :: a1, a2e, a2o, a2es
@@ -125,10 +129,10 @@ SUBROUTINE herman_skillman_int(mesh,func,rab,asum)
         func(i1)=asum+a1*rab(i1)
         a1=a1-a2es+8.0d0*a2o+5.0d0*a2e
         func(i)=asum+a1*rab(i)
-     END DO
+     ENDDO
      asum=func(i)
      a1=0.0d0
-  END DO
+  ENDDO
   !
   RETURN
 END SUBROUTINE herman_skillman_int