!This programs reads a set of geometries, and print the tangent
! to the path they form

SUBROUTINE CalcTangent(XgeomF,NGeomS,XGeom,Coef)

  use Path_module, only: IntCoordI, Nat, Coord, XyzGeomI, Atome, Order, &
                         OrderInv, IndZmat, NMaxPtPath, IntTangent, XyzTangent, &
						 NGeomF, NCoord
  use Io_module, ONly : IoTmp,PathName, Kint,KReal
  use VarTypes

  IMPLICIT NONE
  ! NGeomS: number of geometries.
  INTEGER(KINT), INTENT(IN) :: NGeomS
  REAL(KREAL), INTENT(IN) :: XGeomF(NGeomF),Xgeom(NGeomS),Coef(NGeomS,NCoord)
  
  REAL(KREAL), ALLOCATABLE :: Path(:,:) ! NGeomS,NCoord
  REAL(KREAL), ALLOCATABLE :: PathTan(:,:) ! NGeomF,NCoord
  REAL(KREAL), ALLOCATABLE :: DbgTgt(:),GeomTmp(:) ! NCoord

  INTEGER(KINT) :: i, j, n3at, Idx, iat,igeom,k
  INTEGER(KINT), SAVE :: icall=-1
  INTEGER(KINT) :: NMaxPt
  REAL(KREAL), PARAMETER :: Inf=1e32, h=0.001d0
  CHARACTER(120) :: Line,File
  CHARACTER(LCHARS) :: TITLE
  LOGICAL :: Debug=.False.
  REAL(KREAL) :: utmp,der

  LOGICAL, EXTERNAL :: valid

  debug=valid('calctangent')

  if (debug) WRITE(*,*) "======================================== Entering CalcTangen ============================"
  n3at=3*nat

  ALLOCATE(Path(NGeomS,Ncoord),PathTan(NGeomF,NCoord),GeomTMP(NCoord))

 icall=icall+2
 write(Line,'(I5)') icall
 File=TRIM(adjustl(PathName)) // '_dbgtgt_' // TRIM(adjustl(Line)) // '.dat'
 OPEN(IOTMP,File=File)

  ! We construct path, depending on COORD:
  ! Path(NGeomS,Ncoord), NGeoms is equal to NGeomI
  ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)
 SELECT CASE(Coord)
    CASE ('ZMAT','MIXED')
     Path(:,:)=IntCoordI(:,:)
  CASE ('CART','HYBRID')
     Path=reshape(XyzGeomI,(/NGeomS,NCoord/))
  CASE ('BAKER')
	 Path(:,:)=IntCoordI(:,:)
  CASE DEFAULT
     WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in CalcTangent. STOP"
     STOP
  END SELECT


  if (debug) THEN
     WRITE(*,*) 'DBG CalcTangent. Geometries from Path'
     DO I=1,NGeomS
        GeomTmp=Path(I,1:NCoord)
        WRITE(Title,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS
        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
!         IF (COORD.EQ.'ZMAT') THEN
!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
!                 Nom(Atome(OrderInv(indzmat(2,2)))),Path(I,1)
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
!                 Nom(Atome(OrderInv(indzmat(3,2)))),Path(I,2), &
!                 Nom(Atome(OrderInv(indzmat(3,3)))),Path(I,3)
!            Idx=4
!            DO Iat=4,Nat
!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
!                    Nom(Atome(OrderInv(indzmat(iat,2)))),Path(I,Idx), &
!                    Nom(Atome(OrderInv(indzmat(iat,3)))),Path(I,Idx+1), &
!                    Nom(Atome(OrderInv(indzmat(iat,4)))),Path(I,Idx+2)
!               Idx=Idx+3
!            END DO
!         ELSE
!            WRITE(*,*) Nat
!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
!            DO Iat=1,Nat
!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
!                    Path(I,3*Iat-2:3*Iat)
!            END DO
!         END IF
     END DO
  END IF

!Debug info
 IF (debug.AND.(COORD.EQ.'ZMAT')) THEN
     ALLOCATE(DbgTgt(NCoord))
     WRITE(*,*) 'DBG CalcTangent. Path written in ' // TRIM(File)
     WRITE(*,'(A,12(1X,F10.5))') "DBG CalcTangent, XgeomF:",XGeomF

     NMaxPt=min(200,NMaxPtPath)-1
    DO K=0,NMaxPt
       utmp=K*Xgeom(NGeomS)/NMaxPt
       DO J=1,NCoord
       Call Splint(utmp,DbgTgt(J),NGeomS,Xgeom,Path(1,J),Coef(1,J))
       END DO
       WRITE(IOTMP,'(12(1X,F10.5))') utmp,DbgTgt
    END DO
    DEALLOCATE(DbgTgt)
 END IF

!We construct the derivatives
DO I=1,NGeomF
  DO J=1,NCoord
        Call SplinDer(XgeomF(I),PathTan(I,J),NGeomS,Xgeom,Path(1,J),Coef(1,J))
  END DO
END DO

 SELECT CASE(Coord)
    CASE ('ZMAT','MIXED')
     IntTangent=PathTan
  CASE ('CART','HYBRID')
     XyzTangent=Pathtan
  END SELECT


  if (debug) THEN
     WRITE(*,*) 'DBG CalcTangent. Tangents'
     DO I=1,NGeomF
        GeomTmp=PathTan(I,1:NCoord)
        WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)

        WRITe(*,'(A,12(1X,F10.6))') "Test:",GeomTmp(1:NCoord)

!         IF (COORD.EQ.'ZMAT') THEN
!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
!                 Nom(Atome(OrderInv(indzmat(2,2)))),IntTangent(I,1)
!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
!                 Nom(Atome(OrderInv(indzmat(3,2)))),IntTangent(I,2), &
!                 Nom(Atome(OrderInv(indzmat(3,3)))),IntTangent(I,3)
!            Idx=4
!            DO Iat=4,Nat
!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
!                    Nom(Atome(OrderInv(indzmat(iat,2)))),IntTangent(I,Idx), &
!                    Nom(Atome(OrderInv(indzmat(iat,3)))),IntTangent(I,Idx+1), &
!                    Nom(Atome(OrderInv(indzmat(iat,4)))),IntTangent(I,Idx+2)
!               Idx=Idx+3
!            END DO
!         ELSE
!            WRITE(*,*) Nat
!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
!            DO Iat=1,Nat
!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
!                    XyzTangent(I,3*Iat-2:3*Iat)
!            END DO
!         END IF
     END DO
  END IF


DeALLOCATE(Path) 

 CLOSE(IOTMP)
  if (debug) WRITE(*,*) "======================================== Ending CalcTangent ============================"

END SUBROUTINE CalcTangent

