  SUBROUTINE egrad_test(nat,e,geom,grad) 
! This program reads cartesian coordinates from the file Ext.in
! and print the energy (in au) and gradient (in au/Ang) in the file Ext.out
! format D25.15

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)

       INTEGER(KINT), INTENT(IN) :: Nat
       REAL(KREAL), INTENT(OUT) :: E,grad(Nat*3)
       REAL(KREAL), INTENT(IN) :: Geom(Nat,3)

!     Bohr --> Angstr
      real(KREAL), parameter :: BOHR   = 0.52917726D+00
!
! Parameters to define the surface
      INTEGER(KINT), DIMENSION(6), PARAMETER :: IECOEF = (/-1,9,-45,45,-9,1/)
      INTEGER(KINT), DIMENSION(6), PARAMETER :: ISCOEF = (/-3,-2,-1,1,2,3/)
      REAL(KREAL), PARAMETER :: hh=0.001d0
      
! Variables
       INTEGER(KINT) :: i,j,iat,jat
       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)
       REAL(KREAL) :: xp,yp

       CHARACTER(132) :: Line
       LOGICAL :: Debug
  
      real(KREAL), parameter :: zero = 0.0_KREAL,   one = 1.0_KREAL
      REAL(KREAL), external :: Exyz 


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

  debug=valid('egrad_test')

 if (debug) WRITE(*,*) '================ Entering Egrad_test ==================='
 if (debug) THEN
    WRITE(*,*) "Cartesian Geometry"
    DO I=1,Nat
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Geom(i,1:3)
    END DO
 END IF
      
      ALLOCATE(XyZTmp(Nat,3),GradTmp(Nat,3))

      e=Exyz(nat,Geom)

! We now calculate the gradients using numerical derivatives
      Grad=0.d0
      GradTmp=0.d0
      do iat=1,3
         do jat=1,nat
            xyztmp=geom
            do i=1,6
               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh
               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*Exyz(nat,xyztmp)
            end do
         end do
      end do
      gradTmp=gradTmp/(60.*hh)
      do iat=1,nat
       grad(3*iat-2:3*iat)=gradTmp(iat,1:3)
      end do

 if (debug) THEN
    WRITE(*,*) "Cartesian gradient GradTmp"
    DO I=1,Nat
       WRITE(*,'(1X,I5,3(1X,F12.6))') I,Gradtmp(i,1:3)
    END DO
    WRITE(*,*) "Cartesian gradient Grad"
    WRITE(*,'(50(1X,F12.6))') Grad

 END IF


      deallocate(xyztmp,gradTmp)


 if (debug) WRITE(*,*) '================ Exiting Egrad_test ==================='

! ======================================================================
      end 



      function Exyz(natoms,Xyz)

        use Path_module, only : order

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


      INTEGER(KINT) ,INTENT(IN) :: natoms
      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
      REAL(KREAL) :: Exyz

      INTEGER(KINT) :: i
      REAL(KREAL) :: dCH,dNH,dCN

      real(KREAL) :: autoA,autoeV
      parameter (autoA=0.52917715d0)
      parameter (autoeV=27.21183d0)


      real(KREAL) :: Q(3)
      real(KREAL) :: TT,gR,pR

      real(KREAL) :: mB,mC,mH
      real(KREAL) :: T,RN,RC,r1,r2,r3

      real(KREAL) :: Z1,Z12,V1
      real(KREAL) :: Z2,Z22,V2
      real(KREAL) :: Z3,Z32,V3

      real(KREAL) :: S1,S2,S3,S12,S22,S32,POLY
      real(KREAL) :: E1,E2,E3,HY1,HY2,HY3,SWITCH,V123

      LOGICAL :: Debug
      LOGICAL, SAVE :: First=.TRUE.

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

  debug=valid('exyz')


      
      dCN=0.
      dCH=0.
      dNH=0.
      do i=1,3
         dCN=dCN+(xyz(Order(3),i)-xyz(order(2),i))**2
         dCH=dCH+(xyz(order(1),i)-xyz(order(2),i))**2
         dNH=dNH+(xyz(order(1),i)-xyz(order(3),i))**2
      end do
      R1=sqrt(dCH)
      R2=sqrt(dCN)
      R3=sqrt(dNH)

      IF (debug.AND.First) THEN
         WRITE(*,*) "Exyz in Test, dCH,dCN,dNH:",R1,R2,R3
      END IF

!      R1=R1*autoA
!      R2=R2*autoA
!      R3=R3*autoA

!.....CH
      Z1=R1-1.0823
      Z12=Z1*Z1
      V1=-2.8521*(1.0+5.5297*Z1+8.7166*Z12+5.3082*Z1*Z12)*dexp(-5.5297*Z1)
!.....CN
      Z2=R2-1.1718
      Z22=Z2*Z2 
      V2=-7.9282*(1.0+5.2448*Z2+7.3416*Z22+4.9785*Z2*Z22)*dexp(-5.2448*Z2)
!.....NH
      Z3=R3-1.0370
      Z32=Z3*Z3
      V3=-3.9938*(1.0+3.0704*Z3)*dexp(-3.0704*Z3)

!     write(6,*) v1,v2,v3

!.....THREE BODY TERMS
      Z1=R1-1.9607
      Z2=R2-2.2794
      Z3=R3-1.8687
      S1=0.4436*Z1+0.6091*Z2+0.6575*Z3
      S2=-.8941*Z1+0.2498*Z2+0.3718*Z3
      S3=0.0622*Z1-0.7527*Z2+0.6554*Z3

      S12=S1*S1
      S22=S2*S2
      S32=S3*S3
      POLY=-3.0578*(1.0+1.9076*S1-0.5008*S2-0.0149*S3+0.6695*S12  &
          -1.3535*S22-1.0501*S32+0.2698*S1*S2-1.1120*S1*S3        &
          +1.9310*S2*S3-0.0877*S1*S12+0.0044*S2*S22+0.0700        &
          *S3*S32+0.0898*S12*S2-1.0186*S1*S22-0.0911*S12          &
          *S3+0.0017*S1*S32+0.4567*S22*S3-0.8840*S2*S32           &
          +0.3333*S1*S2*S3-0.0367*S12*S12+0.4821*S22*S22          &
          +0.2564*S32*S32-0.0017*S12*S1*S2-0.2278*S12*S22         &
          -0.1287*S1*S2*S22+0.1759*S1*S12*S3-0.0399*S12*S32       &
          -0.1447*S1*S3*S32-0.3147*S2*S22*S3+0.1233*S22*S32       &
          +0.3161*S2*S3*S32+0.0919*S12*S2*S3-0.0954*S1*S22*S3     &
          +0.1778*S1*S2*S32-0.1892*S22*S22*S2)


      E1=dexp(3.9742*Z1/2.)
      E2=dexp(4.3688*Z2/2.)
      E3=dexp(1.5176*Z3/2.)
      HY1=(E1-1.0/E1)/(E1+1.0/E1)
      HY2=(E2-1.0/E2)/(E2+1.0/E2)
      HY3=(E3-1.0/E3)/(E3+1.0/E3)
      SWITCH=(1.0-HY1)*(1.0-HY2)*(1.0-HY3)
      V123=SWITCH*POLY

!     write(6,*) v123,V1+V2+V3+V123

      Exyz=(V1+V2+V3+V123)/autoeV

      First=.FALSE.

      return
      end function Exyz

      function Exyz2(natoms,Xyz)

   IMPLICIT NONE
     integer, parameter :: KINT = kind(1)
     integer, parameter :: KREAL = kind(1.0d0)


      INTEGER(KINT) ,INTENT(IN) :: natoms
      REAL(KREAL) ,INTENT(IN) :: Xyz(natoms,3)
      REAL(KREAL) :: Exyz

     INTEGER(KINT) :: i,j,iat,jat
     REAL(KREAL) :: Exyz2

      real(KREAL) :: autoA,autoeV
      parameter (autoA=0.52917715d0)
      parameter (autoeV=27.21183d0)


      real(KREAL) :: Q(3)
      real(KREAL) :: TT,gR,pR

      real(KREAL) :: mB,mC,mH
      real(KREAL) :: T,RN,RC,r1,r2,r3

      real(KREAL) :: Z1,Z12,V1
      real(KREAL) :: Z2,Z22,V2
      real(KREAL) :: Z3,Z32,V3

      real(KREAL) :: S1,S2,S3,S12,S22,S32,POLY
      real(KREAL) :: E1,E2,E3,HY1,HY2,HY3,SWITCH,V123

      Exyz=0.
      do jat=1,natoms
         do iat=1,3
         Exyz=Exyz+(3*iat-3+jat)*xyz(jat,iat)**2
         end do
      end do

      Exyz2=Exyz

      return
      end function Exyz2
