      subroutine CalcHess (nfree,geom, hess,Igeom,Iopt)

! This routine calculate the Hessian (or its inverse) numerically using Egrad.
! It should be general and should work with any combination of Coord and Prog
! however, it has not been tested with Prog different from 'test'

! It comes mostly from the hess.f subroutine of the IRC package.


  use Path_module, only : Nat,Coord,Hinv

  use Io_module
!
      IMPLICIT NONE

!!!!!!!!!!! Parameters
  ! NFree==NCoord: Number of the degrees of freedom
  ! IGeom: index of the geometry.
  ! Iopt: step in the optimization process
  INTEGER(KINT), INTENT (IN) :: NFree, IGeom, IOpt
! Geometry at which Hessian is calculated (or its inverse)
   REAL(KREAL), INTENT(IN) :: Geom(Nfree)
! Hessian in compacted form
      REAL(KREAL), INTENT(OUT) ::  Hess(Nfree,Nfree)
!
! Parameters used to calculate the Hessian
! numerically with increasing accuray
      INTEGER(KINT), DIMENSION(6,0:3), PARAMETER :: IECOEF  &
      = RESHAPE(SOURCE=(/           &
          -1,1,0,0,0,0,             &
          -1,1,0,0,0,0,             &
          1,-8,8,-1,0,0,            & 
          -1,9,-45,45,-9,1          &
       /), SHAPE=(/6,4/))

      INTEGER(KINT), DIMENSION(6,0:3), PARAMETER :: ISCOEF   &
      = RESHAPE(SOURCE=(/    &
         0,1,0,0,0,0,         &
          -1,1,0,0,0,0,       & 
          -2,-1,1,2,0,0,      &
           -3,-2,-1,1,2,3     &
         /),SHAPE=(/6,4/))
      INTEGER(KINT), DIMENSION(0:3), PARAMETER :: NbDiv = (/1,2,12,60/)
      INTEGER(KINT), DIMENSION(0:3), PARAMETER :: NbCoef = (/2,2,4,6/)


      INTEGER(KINT), PARAMETER :: HessAcc=2
      REAL(KREAL), PARAMETER :: StepHess=0.01d0
! ======================================================================
      logical           :: debug

      REAL(KREAL) :: Ener
      INTEGER(KINT) :: I, J, IAT, JAT, IJ0, IJ1, IConst

       REAL(KREAL), ALLOCATABLE :: Temp(:),Temp2(:,:) 


       REAL(KREAL), ALLOCATABLE :: Value(:), ValueOld(:), GradTMP(:)
       REAL(KREAL), ALLOCATABLE :: RTmp1(:,:), RTmp3(:,:)  ! Nat,3
       REAL(KREAL), ALLOCATABLE :: RTmp2(:)  ! NFree
! ======================================================================

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID
  END INTERFACE

! ======================================================================

      hess=0.d0

      debug=valid('CALCHESS').OR.valid('hess')

      if (debug) THEN
         WRITE(*,*) 'Entering CalcHess'
         WRITE(*,*) 'StepHess=',StepHess
         WRITE(*,*) "HessAcc=",HessAcc
      END IF

      allocate(Value(Nfree), ValueOld(Nfree), GradTmp(Nfree))
      ALLOCATE(RTmp1(Nat,3),RTmp3(Nat,3),RTmp2(NFree))
      allocate(temp2(Nfree,Nfree))

      RTmp1=0.d0
      RTmp3=0.d0
      RTmp2=Geom

!     Now, we calculate the second derivatives for the active variables               
      Value=Geom
      ValueOld=Geom

      Do IConst=1,Nfree
            DO I=1,NbCoef(HessAcc)
               Value(IConst)=valueOld(IConst)+StepHess*ISCOEF(I,HessAcc)
!     We 'propagate this new value in the Zmatrix'
               call egrad(ener,value,gradtmp,NFree,IGeom,Iopt,RTmp1,RTmp2,RTmp3,.FALSE.)
               Hess(1:Nfree,IConst)= Hess(1:Nfree,IConst)+   &
                    gradtmp(1:Nfree)*IECOEF(I,HessAcc)
            end do
         END DO
        deallocate(Value, ValueOld, GradTmp)
         
        HESS=Hess/(NbDiv(HessAcc)*StepHess)

         if (debug) then
            WRITE(*,*) 'Non symmetrized hessian'
            DO iat=1,Nfree
               WRITE(*,'(15(1X,F12.5))') Hess(iat,1:Nfree)
            END DO
         END if

         DO i=1,Nfree
           Temp2(i,i)=Hess(i,i)
           DO j=i+1,Nfree
             Temp2(i,j)=0.5*(Hess(i,j)+Hess(j,i))
             Temp2(j,i)=Temp2(i,j)
           END DO
         END DO

         Hess=Hess-Temp2

         if (debug) then
            WRITE(*,*) 'Symmetrized hessian'
            DO iat=1,Nfree
               WRITE(*,'(15(1X,F12.5))') Temp2(iat,1:Nfree)
            END DO

            WRITE(*,*) 'Difference between non-syma and Symmetrized hessian'
            DO iat=1,Nfree
               WRITE(*,'(15(1X,F12.5))') Hess(iat,1:Nfree)
            END DO
         end if

! We calculate inverse of the Hessian only if HInv=T

         Hess=Temp2
      if (Hinv) then
        call Geninv(Nfree,Hess,Temp2,NFree)
        Hess=Temp2
      end if

      deallocate(RTmp1,RTmp2,RTmp3,Temp2)

! ======================================================================
      end subroutine CalcHess
