  SUBROUTINE egrad_test_2D(nat,e,geom,grad) 
! This program computes the energy and the gradient
! to cheat PATH we use a 3 atoms system

!----------------------------------------------------------------------
!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 
!  Centre National de la Recherche Scientifique,
!  Université Claude Bernard Lyon 1. All rights reserved.
!
!  This work is registered with the Agency for the Protection of Programs 
!  as IDDN.FR.001.100009.000.S.P.2014.000.30625
!
!  Authors: P. Fleurat-Lessard, P. Dayal
!  Contact: optnpath@gmail.com
!
! This file is part of "Opt'n Path".
!
!  "Opt'n Path" is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Affero General Public License as
!  published by the Free Software Foundation, either version 3 of the License,
!  or (at your option) any later version.
!
!  "Opt'n Path" is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU Affero General Public License for more details.
!
!  You should have received a copy of the GNU Affero General Public License
!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.
!
! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,
! for commercial licensing opportunities.
!----------------------------------------------------------------------

   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.01d0
      INTEGER(KINT), PARAMETER :: IOTMP=99
! Variables
       INTEGER(KINT) :: i, iat, jat,j
       REAL(KREAL), ALLOCATABLE :: Xyztmp(:,:),GradTmp(:,:)

       LOGICAL :: Debug
       LOGICAL, save :: First
  
      REAL(KREAL) :: x,y



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


      function E2D(natoms,Xyz)

        use Path_module, only : order
        use Io_module, only : IOTMP

   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) :: E2D
    END function E2D

     
  END INTERFACE

  debug=valid('egrad_test_2D')

 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))

      if (first) THEN
         OPEN(IOTMP,File="tt")
         DO I=1, 501
            DO j=1,501
               x=(i-1.)/100.
               y=(j-1.)/100.
               xyztmp(2,1)=x
               xyztmp(3,1)=y
               write(IOTMP,*) x,y,E2D(nat,xyztmp)
            END DO
! we skip a blank line
            WRITe(IOTMP,*) 
         END DO
         CLOSE(IOTMP)
         First=.FALSE.
      END IF

      e=E2D(nat,Geom)

! We now calculate the gradients using numerical derivatives
      Grad=0.d0
      GradTmp=0.d0
      do iat=1,1
         do jat=2,nat
            xyztmp=geom
            do i=1,6
               xyztmp(jat,iat)=geom(jat,iat)+ISCoef(i)*hh
               gradTmp(jat,iat)=gradTmp(jat,iat)+IECoef(i)*E2D(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_2D ==================='

! ======================================================================
   end SUBROUTINE egrad_test_2D



      function E2D(natoms,Xyz)

        use Io_module, only : IOTMP

   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) :: E,x,y, E2D
      REAL(KREAL), save :: Data(501,501)

      INTEGER(KINT) :: i,j

      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('e2D')

 if (First) THEN
! we read the data
   open(IOTMP,File="CCSDT.subtr.reg")
   DO I=1, 501
    DO j=1,501
     READ(IOTMP,*) x,y,E
     Data(int(x*100+0.1)+1,int(y*100+0.1)+1)=E
     write(*,*) x,y,int(x*100+0.1)+1, int(y*100+0.1)+1,E
    END DO
! we skip a blank line
     READ(IOTMP,*) 
  END DO
  CLOSE(IOTMP)
  First=.FALSE.
 END IF

    x=Xyz(2,1)
    y=Xyz(3,1)
    E2D=Data(int(x*100+0.1)+1,int(y*100+0.1)+1)

      return
    end function E2D
