!================================================================
!
! This subroutine create a path using starting geometries
! it is 'original' in the sense that the path is generated
! by interpolating internal coordinates instead of cartesian.
! Comes from the Extrapol.f subroutine of Cart suite.
! Rewritten in F90 to be more flexible
!
!================================================================

SUBROUTINE PathCreate(IOpt)

  use Io_module
  use Path_module, only : Nat, NGeomI, NCoord, NGeomF, IGeomRef, NMaxL, IReparam, IReparamT, &
       Coord, Frozen, Cart, NCart, NFroz, XyzGeomI, atome, r_cov, fact, &
       IndZmat, Renum, Order, OrderInv, IntCoordI, IntCoordF,Pi, BMat_BakerT, Nom, &
       Hess,IntFroz, ISpline, IntTangent, XyzGeomF, NPrim,Xprimitive_t, OptGeom,  &
       UMatF, UMat_local, XyzTangent, Linear, Align, FrozAtoms
  ! BMat_BakerT (3*Nat,NCoord), XyzGeomI(NGeomI,3,Nat), IGeomRef=-1 (default value)
  ! IntCoordI(NGeomI,NCoord)

  IMPLICIT NONE

  REAL(KREAL), PARAMETER :: Inf=1e32
  INTEGER(KINT), INTENT(IN) :: IOpt

  INTEGER(KINT) ::  at,long
  REAL(KREAL), ALLOCATABLE :: val_zmat(:,:), val_zmatTmp(:,:) ! (3,Nat)
  REAL(KREAL), ALLOCATABLE :: X0(:), Y0(:), Z0(:) ! Nat
  REAL(KREAL), ALLOCATABLE ::  XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
  !  REAL(KREAL), ALLOCATABLE ::  XyzTmp3(:,:) ! (3,Nat)
  REAL(KREAL), ALLOCATABLE ::  IntCoordTmp(:) ! (NCoord)
  INTEGER(KINT), ALLOCATABLE :: IndZmat0(:,:) !Nat,5
  REAL(KREAL), ALLOCATABLE :: Coef(:,:) ! (NGeomI,NCoord)
  ! NPrim=Number of primitive coordinates.
  REAL(KREAL), ALLOCATABLE :: XGeom(:) ! NGeomI
  REAL(KREAL), ALLOCATABLE :: XGeomF(:) ! NGeomF
  REAL(KREAL), ALLOCATABLE :: Hess_int(:,:) ! Hess_int (3*Nat,3*Nat), used for Baker case
  INTEGER(KINT) :: NGeomS

  REAL(KREAL) :: Rmsd,MRot(3,3), Delta, normDeltaX_int
  REAL(KREAL), PARAMETER ::  LimAngle=145.d0,TwoPi=360.d0
  REAL(KREAL)  ::s,ds,dist,a_val,d
  INTEGER(KINT) :: I,J,K,NbCycle,igeom,nb,ig,iat
  INTEGER(KINT) :: Itmp, Jat, L
  INTEGER(KINT) :: Idx,IBeg
  CHARACTER(LCHARS) :: Title
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
  LOGICAL ::  debug, print
  LOGICAL, SAVE :: First=.TRUE.
  ! Frozen contains the indices of frozen atoms
  LOGICAL, ALLOCATABLE :: FCart(:) ! Nat
  INTEGER, ALLOCATABLE :: AtCart(:) !Nat
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:),LiaisonsIni(:,:) ! (Nat,0:NMaxL)
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !nat
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !nat,nat
  INTEGER(KINT) :: NbFrag,NFragRef


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

     FUNCTION test_zmat(na,ind_zmat)

       USE Path_module

       logical :: test_zmat
       integer(KINT) :: na
       integer(KINT) :: ind_zmat(Na,5)
     END FUNCTION TEST_ZMAT
  END INTERFACE



  debug=valid("pathcreate")
  print=valid("printgeom")
  if (debug) Call header("Entering PathCreate")
  if (debug.AND.First) WRITE(*,*) "=     First call          ="

  if (debug) WRITE(*,*) "Iopt,IReparam=",Iopt,IReparam

  ALLOCATE(Coef(NGeomI,NCoord),val_zmat(Nat,3),val_zmatTMP(Nat,3))
  ALLOCATE(X0(Nat),Y0(Nat),Z0(Nat),Indzmat0(Nat,5))
  ALLOCATE(X(Nat),Y(Nat),Z(Nat),Xgeom(NGeomI),XgeomF(NGeomF))
  ! ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),XyzTmp3(3,Nat))
  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3))
  ALLOCATE(IntCoordTmp(NCoord))

  Do I=1,NGeomI
     XGeom(I)=FLoat(I)-1.d0
  END DO

  ! First iteration of the optimization:
  IF (First) THEN ! matches at L484
     if (debug) WRITE(*,*) "first iteration in PathCreate.90, L93"

     ! This is the first call, so we might have to generate a Zmat
     First=.FALSE.

     IF ( ((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED") &
          .OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)) THEN

        ! User has not defined the Reference geometry, we have to choose it ourselves!
        ! This is done by counting the number of fragments of each geometry. The 
        ! reference geometry is the one with the least fragments. When there are 
        ! frozen or cart atoms, they are discarded to count the fragments.

        ALLOCATE(Liaisons(Nat,0:NMaxL))
        ALLOCATE(Fragment(nat),NbAtFrag(nat),FragAt(nat,nat),FCart(Nat))
        FCart=.FALSE.
        ! FCart(Nat) was under IF below earlier but is being accessed after this loop.
        IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN
           ALLOCATE(AtCart(Nat),LiaisonsIni(Nat,0:NMaxL))
           FCart=.FALSE.
           NCart=0
           Idx=1
           DO WHILE (Frozen(Idx).NE.0)
              IF (Frozen(Idx).LT.0) THEN
                 DO I=Frozen(Idx-1),abs(Frozen(Idx))
                    IF (.NOT.FCart(I)) THEN
                       NCart=NCart+1
                       AtCart(NCart)=I
                       FCart(I)=.TRUE.
                    END IF
                 END DO
              ELSE IF (.NOT.FCart(Frozen(Idx))) THEN
                 FCart(Frozen(Idx))=.TRUE.
                 NCart=NCart+1
                 AtCart(NCart)=Frozen(Idx)
              END IF
              Idx=Idx+1
           END DO ! matches DO WHILE (Frozen(Idx).NE.0)
           NFroz=NCart
           Idx=1
           DO WHILE (Cart(Idx).NE.0)
              IF (Cart(Idx).LT.0) THEN
                 DO I=Cart(Idx-1),abs(Cart(Idx))
                    IF (.NOT.FCart(I)) THEN
                       NCart=NCart+1
                       AtCart(NCart)=I
                       FCart(I)=.TRUE.
                    END IF
                 END DO
              ELSEIF (.NOT.FCart(Cart(Idx))) THEN
                 FCart(Cart(Idx))=.TRUE.
                 NCart=NCart+1
                 AtCart(NCart)=Cart(Idx)
              END IF
              Idx=Idx+1
           END DO !matches DO WHILE (Cart(Idx).NE.0)
           IF (debug) THEN
              WRITE(*,'(1X,A,I4,A,I4,A)') "DBG PathCreate - GeomRef :: Found ", &
                   NFroz," frozen atoms, and a total of ",NCart, &
                   " atoms described in cartesian coordinates"
              WRITE(*,*) "AtCart:",AtCart(1:NCart)
           END IF
        END IF ! IF ((Frozen(1).NE.0).OR.(Cart(1).NE.0)) THEN

        NFragRef=2*Nat
        DO IGeom=1,NGeomI
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)
           ! we calculate the connectivities
           Call CalcCnct(nat,atome,x,y,z,LIAISONS,r_cov,fact)
           ! if needed, we get rid of links between cart or frozen atoms
           IF (NCart.GE.2) THEN
              LiaisonsIni=Liaisons
              DO I=1,NCart
                 Iat=AtCart(I)
                 if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
                      (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
                 DO J=1,LiaisonsIni(Iat,0)
                    Jat=LiaisonsIni(Iat,J)
                    IF (FCart(Jat)) THEN
                       Call Annul(Liaisons,Iat,Jat)
                       Call Annul(Liaisons,Jat,Iat)
                       Call Annul(LiaisonsIni,Jat,Iat)
                       Liaisons(Iat,0)=Liaisons(Iat,0)-1
                       Liaisons(Jat,0)=Liaisons(Jat,0)-1
                       LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
                    END IF
                 END DO
              END DO
           END IF ! matches IF (NCart.GE.2) THEN
           ! we now calculate the number of fragments.
           Call Decomp_frag(nat,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt)
           IF (NbFrag.LT.NFragRef) THEN
              NFragRef=NbFrag
              ! The reference geometry, IGeomRef, is determined based on the least number
              ! of fragments (here if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD
              !.EQ."MIXED")).AND.(IGeomRef.LE.0))), otherwise it still has the value of 
              ! IGeomRef=-1.
              IGeomRef=IGeom
           END IF
        END DO ! matches DO IGeom=1,NGeomI
        DEALLOCATE(FCart)
     END IF ! matches IF (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD.EQ."MIXED")
     ! .OR.(COORD.EQ."BAKER")).AND.(IGeomRef.LE.0)).

     if (debug) WRITE(*,*) "DBG PathCreate : IGeomRef= ",IGeomRef

     ! we now compute the internal coordinates for this geometry !
     ! The reference geometry, IGeomRef, is determined based on the least number
     ! of fragments if (((COORD.EQ."ZMAT").OR.(COORD.EQ."HYBRID").OR.(COORD
     !.EQ."MIXED")).AND.(IGeomRef.LE.0)), otherwise it has the value of 
     ! IGeomRef=-1.
     x(1:Nat) = XyzGeomI(IGeomRef,1,1:Nat) ! XyzGeomI(NGeomI,3,Nat)
     y(1:Nat) = XyzGeomI(IGeomRef,2,1:Nat)
     z(1:Nat) = XyzGeomI(IGeomRef,3,1:Nat)

     IndZmat=0
     IndZmat0=0

     SELECT CASE(COORD)
     CASE ('BAKER')
        ! atome(Nat), r_cov(0:Max_Z)
        ! Frozen(Nat) contains the indices of frozen atoms
        Call Calc_baker(atome,r_cov,fact,frozen,IGeomRef)
        if (debug) THEN
           WRITE(*,*) "UMat_local after Calc_baker"
           DO I=1,3*Nat-6
              WRITE(*,'(50(1X,F12.8))') UMat_local(:,I)
           END DO
        END IF
        DO IGeom=1,NGeomF
           UMatF(IGeom,:,:)=UMat_Local(:,:)
        END DO

        ALLOCATE(Xprimitive_t(NPrim))
        ! x, y, z corresponds the reference geometry.
        DO I=1,Nat
           Order(I)=I
           OrderInv(I)=I
        END DO
        x0=x
        y0=y
        z0=z
     CASE ('ZMAT','HYBRID')
        ! IN case we are using Hybrid or zmat coordinate system, we have to 
        ! generate the internal coordinates with a Z-Matrix
        ! IN case we have frozen atoms, we generate a constrained Zmat
        IF (Frozen(1).NE.0) THEN 
           Renum=.True.
           call Calc_zmat_const_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact,frozen)
        ELSE
           !no zmat... we create our own
           call Calc_zmat_frag(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact)
           if (valid('old_zmat')) THEN 
              call Calc_zmat(Nat,atome,IndZmat0,val_zmat,x,y,z,r_cov,fact)
              WRITE(*,*) "DBG PathCreate Calc_Zmat.... STOP"
              STOP
           END IF
        END IF

        if (debug) WRITE(*,*) 'Back from Calc_zmat'
        IF (.NOT.test_zmat(nat,IndZmat0)) THEN
           WRITE(IOOUT,*) "I cannot generate a valid Zmat... help"
           STOP
        END IF

        ! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z
        ! However, we need it with the new order so that we can use the zmat 
        ! coordinates to generate cartesian coordinates.
        ! So, we have to reorder it...
        DO I=1,Nat
           Order(IndZmat0(I,1))=I
           OrderInv(I)=IndZmat0(I,1)
           ! I let X0 being the Cartesian coordinates of the reference geometry. Before it was
           ! those of the First geometry (which was also the reference one). This might change!
           X0(i)=X(IndZmat0(i,1))
           Y0(i)=Y(IndZmat0(i,1))
           Z0(i)=Z(IndZmat0(i,1))
        END DO
        IndZmat(1,1)=Order(IndZmat0(1,1))
        IndZmat(2,1)=Order(IndZmat0(2,1))
        IndZmat(2,2)=Order(IndZmat0(2,2))
        IndZmat(3,1)=Order(IndZmat0(3,1))
        IndZmat(3,2)=Order(IndZmat0(3,2))
        IndZmat(3,3)=Order(IndZmat0(3,3))
        DO I=4,Nat
           DO J=1,4
              IndZmat(I,J)=Order(IndZmat0(I,J))
           END DO
        end do
        IF (DEBUG) THEN
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat"
           DO I=1,Nat
              WRITE(IOOUT,*) (IndZmat(I,J),J=1,5)
           END DO
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat0"
           DO I=1,Nat
              WRITE(IOOUT,*) (IndZmat0(I,J),J=1,5)
           END DO

        END IF

        ! It is the first call, we have to create all internal coordinates
        if (debug) WrITE(*,*) "DBG PathCreate L302: First Call, we generate internal coords"
        DO igeom=1,NgeomI
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)

           if (debug) WRITE(*,*) "DBG PathCreate L308: we generate internal coords for geom",IGeom

           Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp)
           IntCoordI(IGeom,1)=val_zmatTmp(2,1)
           IntCoordI(IGeom,2)=val_zmatTmp(3,1)
           IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180.
           Idx=4
           DO i=4,nat
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
              Idx=Idx+1
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
              Idx=Idx+1

              ! Some tests to check that the dihedral angles are similar... this is
              ! important if they are close to Pi.
              Delta=0.

              if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
                 if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
                    Delta=-1.0d0*TwoPi
                 ELSE
                    Delta=TwoPi
                 END IF
                 if (debug) THEN
                    WRITE(*,*) Nom(atome(i)),        &
                         abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
                         val_zmatTMP(i,3),val_zmat(i,3),Delta, &
                         val_zmatTMP(i,3)-Delta
                 END IF
              END IF
              IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
              Idx=Idx+1
           END DO
        END DO  ! End Loop on Igeom

     CASE ('MIXED')
        ! IN case we are using mixed coordinate system, we have to 
        ! we generate the internal coordinates with a Z-Matrix
        Renum=.TRUE.
        ! 
        call Calc_mixed_frag(Nat,atome,IndZmat0,val_zmat,x,y,z, &
             r_cov,fact,frozen,cart,ncart)

        if (debug) WRITE(*,*) 'Back from Calc_Mixed'

        ! How to test this partial zmat ? it must be good as it was automatically generated...
        !        IF (.NOT.test_zmat(nat,IndZmat0)) THEN
        !           WRITE(IOOUT,*) "I cannot generate a valid Zmat... help"
        !           STOP
        !        END IF

        ! The IndZmat0 generated by Calc_zmat refers to the order of the input x,y,z
        ! However, we need it with the new order so that we can use the zmat 
        ! coordinates to generate cartesian coordinates.
        ! So, we have to reorder it...
        DO I=1,Nat
           Order(IndZmat0(I,1))=I
           OrderInv(I)=IndZmat0(I,1)
           ! I let X0 being the Cartesian coordinates of the reference geometry. Before it was
           ! those of the First geometry (which was also the reference one). This might change!
           X0(i)=X(IndZmat0(i,1))
           Y0(i)=Y(IndZmat0(i,1))
           Z0(i)=Z(IndZmat0(i,1))
        END DO
        IndZmat=0
        IndZmat(1:NCart,:)=-1

        DO I=1,Nat
           IndZmat(I,1)=Order(IndZmat0(I,1))
        END DO

        SELECT CASE (NCart)
        CASE (1)
           IndZmat(2,2)=Order(IndZmat0(2,2))
           IndZmat(3,2)=Order(IndZmat0(3,2))
           IndZmat(3,3)=Order(IndZmat0(3,3))
           IdX=4
        CASE (2)
           IndZmat(3,2)=Order(IndZmat0(3,2))
           IndZmat(3,3)=Order(IndZmat0(3,3))
           Idx=4
        CASE DEFAULT 
           Idx=NCart+1
        END SELECT

        DO I=Idx,Nat
           DO J=1,4
              IndZmat(I,J)=Order(IndZmat0(I,J))
           END DO
        end do
        IF (DEBUG) THEN
           WRITE(IOOUT,*) "DBG PathCreate : IndZmat - Mixed"
           DO I=1,Nat
              WRITE(IOOUT,*) (IndZmat(I,J),J=1,5)
           END DO
        END IF

        ! It is the first call, we have to create all internal coordinates
        !            Idx=1
        !            DO I=1,NCart
        !               IntCoordI(1,Idx)=val_zmat(I,1)
        !               IntCoordI(1,Idx+1)=val_zmat(I,2)
        !               IntCoordI(1,Idx+2)=val_zmat(I,3)
        !               Idx=Idx+3
        !            END DO
        !            SELECT CASE (NCART)
        !            CASE (1)
        !               IntCoordI(1,Idx)=val_zmat(2,1)
        !               IntCoordI(1,Idx+1)=val_zmat(3,1)           
        !               IntCoordI(1,Idx+2)=val_zmat(3,2)*Pi/180.
        !               Idx=Idx+3
        !               IBeg=4
        !            CASE (2)
        !               IntCoordI(1,Idx)=val_zmat(3,1)           
        !               IntCoordI(1,Idx+1)=val_zmat(3,2)*Pi/180.
        !               Idx=Idx+2
        !               IBeg=4
        !            CASE DEFAULT
        !               IBeg=NCart+1
        !            END SELECT
        !            DO i=iBeg,Nat
        !               IntCoordI(1,Idx)=val_zmat(i,1)
        !               Idx=Idx+1
        !               DO j=2,3
        !                  IntCoordI(1,Idx)=val_zmat(i,j)*Pi/180.
        !                  Idx=Idx+1
        !               END DO
        !            END DO

        !            If (debug) WRITE(*,*) "Idx, NCoord",Idx-1,NCoord

        DO Igeom=1,NgeomI
           x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
           y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
           z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)

           Call ConvXyzMixed(Nat,Ncart,x,y,z,IndZmat0,val_zmatTmp(1,1))
           Idx=1
           DO I=1,NCart
              IntCoordI(IGeom,Idx)=val_zmatTmp(I,1)
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(I,2)
              IntCoordI(IGeom,Idx+2)=val_zmatTmp(I,3)
              Idx=Idx+3
           END DO

           SELECT CASE (NCART)
           CASE (1)
              IntCoordI(IGeom,Idx)=val_zmatTmp(2,1)
              Idx=Idx+1
              IntCoordI(IGeom,Idx)=val_zmatTmp(3,1)           
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180.
              Idx=Idx+2
              IBeg=4
           CASE (2)
              IntCoordI(IGeom,Idx)=val_zmatTmp(3,1)           
              IntCoordI(IGeom,Idx+1)=val_zmatTmp(3,2)*Pi/180.
              Idx=Idx+2
              IBeg=4
           CASE DEFAULT
              Ibeg=Ncart+1
           END SELECT
           DO i=IBeg,Nat
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
              Idx=Idx+1
              IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
              Idx=Idx+1

              ! Some tests to check that the dihedral angles are similar... this is
              ! important if they are close to Pi.
              Delta=0.

              if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
                 if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
                    Delta=-1.0d0*TwoPi
                 ELSE
                    Delta=TwoPi
                 END IF
                 if (debug) THEN
                    WRITE(*,*) Nom(atome(i)), &
                         abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
                         val_zmatTMP(i,3),val_zmat(i,3),Delta, &
                         val_zmatTMP(i,3)-Delta
                 END IF
              END IF
              IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
              Idx=Idx+1
           END DO

        END DO
     CASE ('CART') 
        DO I=1,Nat
           Order(I)=I
           OrderInv(I)=I
        END DO
        X0=X
        Y0=y
        Z0=z
     END SELECT

  ELSE ! First iteration of the optimization. Match at L616

     IGeom=1
     x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
     y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
     z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)

     SELECT CASE (COORD)
     CASE ('HYBRID')
        IndZmat0=IndZmat
        Call ConvXyzZmat(Nat,x,y,z,IndZmat,val_zmat(1,1))
     CASE ('ZMAT')
        IndZmat0=IndZmat
        val_zmat=0.d0
        val_zmat(2,1)=IntCoordI(1,1)
        val_zmat(3,1)=IntCoordI(1,2)
        val_zmat(3,2)=IntCoordI(1,3)*180./Pi

        Idx=4
        DO iat=4,Nat
           val_zmat(iat,1)=IntCoordI(1,idx)
           Idx=Idx+1
           do j=2,3
              val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi
              Idx=Idx+1
           END DO
        END DO

     CASE ('MIXED')
        IndZmat0=IndZmat
        val_zmat=0.d0
        Idx=1
        DO I=1,NCart
           DO J=1,3
              val_zmat(i,j)= IntCoordI(1,Idx)
              Idx=Idx+1
           END DO
        END DO
        SELECT CASE (NCart)
        CASE(1)
           val_zmat(2,1)=IntCoordI(1,Idx)
           val_zmat(3,1)=IntCoordI(1,Idx+1)
           val_zmat(3,2)=IntCoordI(1,Idx+2)*180./Pi
           Idx=Idx+3
           IBeg=4
        CASE(2)
           val_zmat(3,1)=IntCoordI(1,Idx)
           val_zmat(3,2)=IntCoordI(1,Idx)*180./Pi
           Idx=Idx+2
           IBeg=4
        CASE DEFAULT
           IBeg=NCart+1
        END SELECT

        DO Iat=IBeg,Nat
           val_zmat(iat,1)=IntCoordI(1,idx)
           Idx=Idx+1
           do j=2,3
              val_zmat(iat,J)=IntCoordI(1,idx)*180./Pi
              Idx=Idx+1
           END DO
        END DO
     CASE ('BAKER')
        ! Nothing to be done here.
     CASE ('CART')
        ! Nothing to be done here.
     CASE DEFAULT
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L554.STOP"
        STOP
     END SELECT
     X0=x
     y0=y
     z0=z
  END IF  ! First

  ! Now that we have a zmat, we will generate all the IntCoodI corresponding...
  ! First one
  IF (COORD.EQ.'HYBRID') THEN     ! Matches at L680
     DO IGeom=1,NGeomI
        x(1:Nat) = XyzGeomI(IGeom,1,1:Nat)
        y(1:Nat) = XyzGeomI(IGeom,2,1:Nat)
        z(1:Nat) = XyzGeomI(IGeom,3,1:Nat)

        Call ConvXyzZmat(Nat,x,y,z,IndZmat0,val_zmatTmp(1,1))
        IntCoordI(IGeom,1)=val_zmatTmp(2,1)
        IntCoordI(IGeom,2)=val_zmatTmp(3,1)
        IntCoordI(IGeom,3)=val_zmatTmp(3,2)*Pi/180.
        Idx=4
        DO i=4,nat
           IntCoordI(IGeom,Idx)=val_zmatTmp(i,1)
           Idx=Idx+1
           IntCoordI(IGeom,Idx)=val_zmatTmp(i,2)*Pi/180.
           Idx=Idx+1

           ! Some tests to check that the dihedral angles are similar... this is
           ! important if they are close to Pi.
           Delta=0.

           if (abs(val_zmatTMP(i,3)-val_zmat(i,3))>=180.) THEN
              if ((val_zmatTMP(i,3)-val_zmat(i,3))<=-180.) THEN
                 Delta=-1.0d0*TwoPi
              ELSE
                 Delta=TwoPi
              END IF
              if (debug) THEN
                 WRITE(*,*) Nom(atome(i)), &
                      abs(val_zmatTMP(i,3)-val_zmat(i,3)), &
                      val_zmatTMP(i,3),val_zmat(i,3),Delta, &
                      val_zmatTMP(i,3)-Delta
              END IF
           END IF
           IntCoordI(IGeom,Idx)=(val_zmatTmp(i,3)-Delta)/180.*Pi
           Idx=Idx+1
        END DO
     END DO
  END IF ! Matches Coord=Hybrid

  ! PFL 24 Nov 2008 ->
  ! This should not be necessary...
  !   ! if we have frozen atoms, we make sure that they do not move between geometries
  !   IF (IntFroz.NE.0) THEN
  !      if (debug) WRITE(*,*) 'DBG PathCreate L641: Number of frozen coordinates=', IntFroz
  !      DO IGeom=2,NGeomI ! For IOpt .GT. 0, NGeomI is equal to NGeomF
  !         IntCoordI(IGeom,1:IntFroz)=IntCoordI(1,1:IntFroz)
  !      END DO
  !   END IF
  ! -> PFL 24 Nov 2008

  Idx=min(9,NCoord)

  IF (DEBUG.AND.(COORD.NE.'CART')) THEN
     WRITE(*,*) "DBG PathCreate. L685 IntCoordI(IGeom,:)"
     DO I=1,NGeomI
        WRITE(*,'("Geom:",I5,9(1X,F10.4))') I,IntCoordI(I,1:Idx)
     END DO
     ! We write it also in terms of Zmat
     IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'HYBRID')) THEN
        DO I=1,NGeomI
           WRITE(*,*) 'Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I)
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(1,1)))
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(2,1))), &
                IndZmat(2,2),IntCoordI(I,1)
           WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(3,1))), &
                IndZmat(3,2),IntCoordI(I,2),IndZmat(3,3),IntCoordI(I,3)*180./Pi
           Idx=4
           DO J=4,Nat
              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(IndZmat0(J,1))), &
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3),          &
                   IntCoordI(I,Idx+1)*180./Pi, &
                   IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi
              Idx=Idx+3
           END DO

           XyzTmp2=0.
           ! We convert it into Cartesian coordinates
           XyzTmp2(2,1)=IntCoordI(I,1)
           d=IntCoordI(I,2)
           a_val=IntCoordI(I,3)
           if (Nat.GE.3) THEN
              if (IndZmat(3,2).EQ.1)  THEN
                 XyzTmp2(3,1)=XYzTmp2(1,1)+d*cos(a_val)
              ELSE
                 XyzTmp2(3,1)=XyzTmp2(2,1)-d*cos(a_val)
              ENDIF
              XyzTmp2(3,2)=d*sin(a_val)
           ENDIF

           Idx=4
           DO iat=4,Nat
              val_zmatTMP(iat,1)=IntCoordI(I,idx)
              Idx=Idx+1
              DO j=2,3
                 val_zmatTMP(iat,J)=IntCoordI(I,idx)*180./Pi
                 Idx=Idx+1
              END DO
           END DO
           DO iat=4,Nat
              call ConvertZmat_cart(iat,IndZmat,val_zmatTMP, &
                   XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
           END DO
           DO Iat=1,nat
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),&
                   XyzTmp2(iat,1:3)
           END DO

        END DO
     END IF
     IF (COORD.EQ.'MIXED') THEN
        DO I=1,NGeomI
           WRITE(*,*) 'Cart+Zmat for geom ',I,'/',NGeomI,'Xgeom=',Xgeom(I)
           Idx=1
           XyzTmp2=0.
           DO J=1,NCART
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(1,1))), &
                   IntCoordI(I,Idx:Idx+2)
              XyzTmp2(J,1:3)=IntCoordI(I,Idx:Idx+2)
              Idx=Idx+3
           END DO
           SELECT CASE (NCart)
           CASE (1)
              J=2
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
                   IndZmat(J,2),IntCoordI(I,Idx)
              Idx=Idx+1
              J=3
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
                   IntCoordI(I,Idx+1)*180./Pi
              Idx=Idx+2
              Ibeg=4
           CASE (2)
              J=3
              WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
                   IntCoordI(I,Idx+1)*180./Pi
              Idx=Idx+2
              Ibeg=4
           CASE DEFAULT
              IBeg=NCart+1
           END SELECT
           DO J=IBeg,Nat
              IF (J.EQ.2) &
                   WRITE(*,'(1X,A5,3(I5,F11.6))') Nom(Atome(IndZmat0(J,1))), &
                   IndZmat(J,2),IntCoordI(I,Idx),IndZmat(J,3), &
                   IntCoordI(I,Idx+1)*180./Pi, &
                   IndZmat(J,4),IntCoordI(I,Idx+2)*180./Pi
              Idx=Idx+3
           END DO

           ! We convert it into Cartesian coordinates
           IntCoordTmp=IntCoordI(I,1:NCoord)
           Call Mixed2cart(Nat,IndZmat,IntCoordTmp,XyzTmp2)

           DO Iat=1,nat
              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(IndZmat0(Iat,1))),&
                   XyzTmp2(iat,1:3)
           END DO

        END DO
     END IF
  END IF ! matches IF (DEBUG.AND.(COORD.NE.'CART')) THEN

  if (debug) WRITE(*,*) "Before interpolation, PathCreate.f90, L740, IOpt=",IOpt, &
       "ISpline=", ISpline 
  ! Now comes the Interpolation:
  IF ((NGeomI>2).AND.(IOpt.GE.ISpline)) THEN
     Linear=.FALSE.
     ! We have at least three geometries so we can apply cubic spline
     SELECT CASE (COORD)
     CASE ('ZMAT','HYBRID','MIXED')
        DO I=1,NCoord
           ! We compute the spline coefficients
           call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf, Coef(1,I))
           if (debug) THEN
              WRITE(*,*) 'Coef Spline for coord',I
              WRITE(*,*) Coef(1:NGeomI,I)
           END IF
        END DO
     CASE ('BAKER')
        ! We compute the spline coefficients
        ! xGeom(NGeomI), IntCoordI(NGeomI,3*Nat-6) declared in Path_module.f90
        ! can also be used for the Baker's internal coordinates.
        ! we need to get Baker's internal coordinates from cartesian coordinates.
        ! as we have subroutine ConvXyzmat(...) in Z matrix case.
        ! Coef(NGeomI,NCoord) or (NGeomI,N3at) has INTENT(OUT)
        DO I=1,NCoord ! NCoord=3*Nat-6-symmetry_eliminated.
           call spline(xGeom,IntCoordI(1,I),NGeomI,Inf,Inf,Coef(1,I))
           ! IF (debug) THEN
           !   WRITE(*,*) 'Coef Spline for coord',I
           !  WRITE(*,*) Coef(1:NGeomI,I)
           !    END IF
        END DO
     CASE ('CART')

!! PFL 2011 June 22 ->
!! The alignment procedure from Extrapol_cart has been moved
!! here, before computation of the spline coefficients

!  XyzTmp2=Reshape(XyzGeomI(1,:,:),(/Nat,3/),ORDER=(/2,1/))

  if (debug) THEN
     WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Initial geometries"
     DO IGeom=1,NGeomI
        WRITE(*,*) 'XyzGeomI, IGeom=',IGeom
        DO I=1,Nat
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
                (XyzGeomI(IGeom,J,I),J=1,3)
        END DO
     END DO
  END IF


! In order to mesure only the relevant distance between two points
! we align all geometries on the original one

   DO IGeom=1,NGeomI

      XyzTmp2=Reshape(XyzGeomI(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/))
  ! We align this geometry with the original one
  ! PFL 17/July/2006: only if we have more than 4 atoms.
!  IF (Nat.GT.4) THEN
! PFL 24 Nov 2008 ->
! If we have frozen atoms we align only those ones.
! PFL 8 Feb 2011 ->
! I add a flag to see if we should align or not.
! For small systems, it might be better to let the user align himself
      IF (Align) THEN
         if (NFroz.GT.0) THEN
            Call AlignPartial(Nat,x0,y0,z0,                     &
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
                 FrozAtoms,MRot)
         ELSE
            Call  CalcRmsd(Nat, x0,y0,z0,                              &
                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
                 MRot,rmsd,.TRUE.,.TRUE.)
         END IF
! <- PFL 24 Nov 2008
      END IF
! -> PFL 8 Feb 2011
!  END IF

  XyzGeomI(IGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
 END DO

   if (debug) THEN
     WRITE(*,*) "DBG PathCreate - Calc Spline coeff CART- Aligned geometries"
      DO J=1, NGeomI
         WRITE(IOOUT,*) Nat
         WRITE(IOOUT,*) " Aligned geometry ",J,"/",NGeomI
           DO i=1,Nat
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
                   XyzGeomI(J,1:3,I)
           END DO
        END DO
     END IF

!! <- PFL 2011 June 22

        DO I=1,Nat
           ! We compute the spline coefficients
           call spline(xGeom,XyzGeomI(1,1,I),NGeomI,Inf,Inf, Coef(1,3*I-2))
           call spline(xGeom,XyzGeomI(1,2,I),NGeomI,Inf,Inf, Coef(1,3*I-1))
           call spline(xGeom,XyzGeomI(1,3,I),NGeomI,Inf,Inf, Coef(1,3*I))
           if (debug) THEN
              WRITE(*,*) 'Coef Spline for coord',i
              WRITE(*,*) Coef(1:NGeomI,I)
           END IF
        END DO
     CASE DEFAULT
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate.L783.STOP"
        STOP
     END SELECT
  ELSE ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN
     ! On n'a que deux images. On fera une interpolation lineaire,
     ! qui est un cas particulier avec Coef=0
     if (debug.AND.(NGeomI.EQ.2)) WRITE(*,*) "DBG : Only 2 starting geometries"
     if (debug.AND.(Iopt.LE.ISpline)) WRITE(*,*) "DBG : Not enough path -> linear int"
     Coef=0.d0
     Linear=.TRUE.
  END IF ! matches IF ((NGeomI>2).AND.(Iopt.Ge.ISpline)) THEN
  ! Shall we redistribute the points along the path ?
  if (debug) WRITE(*,*) "PathCreate.f90, L795, IOpt=", IOpt, "IReparam=", IReparam

  IF (MOD(IOpt,IReparam).EQ.0) THEN
     ! We generate the Path, in two steps:
     ! i) we calculate the length of the mass weighted (MW) path
     ! ii) we sample the path to get the geometries

     NGeomS=NGeomI
     SELECT CASE (COORD)
     CASE ('ZMAT','HYBRID')
        ! We call the extrapolation once to get the path length
        dist=Inf
        Call ExtraPol_int(IOpt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
        ! we have the length, we calculate the distance between two points

        if (s.LE.1e-6) THEN
           IF (OptGeom.LE.0) THEN
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
              WRITE(*,*) "ERROR !!! Stop !"
              STOP
           END IF
        ELSE
           ! we have the length, we calculate the distance between two points
           dist=s/Real(NGeomf-1)
           if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist
           ! we call it again to obtain the geometries !
           Call ExtraPol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
        END IF
     CASE ('BAKER')
        ! We generate the Path, in two steps:
        ! i) we calculate the length of the mass weighted (MW) path
        ! ii) we sample the path to get the geometries

        ! We call the extrapolation first time to get the path length
        ! iopt: Number of the cycles for the optimization.
        ! s: a real number and has INTENT(OUT), probably returns the 
        ! length of the path.
        ! dist=Inf=> ExtraPol_baker has to calculate the length of the path.
        ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.
        ! Xgeom: Xgeom(NGeomI) has INTENT(IN)
        ! Coef(NGeomI,NCoord) is obtained from spline interpolation subroutine.
        ! XgeomF: XGeomF(NGeomF) has INTENT(OUT)
        ! Which reference coordinate geometry X0, Y0, Z0 to use???
        dist=Inf
        !if (debug) WRITE(*,*) "Before the call of ExtraPol_baker in redistribution &
        !of the points, in PathCreate.f90, L836"
        !Print *, 'Before ExtraPol_baker, PathCreate.90, L848, IntCoordI='
        !Print *, IntCoordI
        Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
        !Print *, 'After ExtraPol_baker, PathCreate.90, L848, IntCoordI='
        !Print *, IntCoordI
        !if (debug) WRITE(*,*) "After the call of ExtraPol_baker in the redistribution &
        ! of the points, in PathCreate.f90, L843"

        ! we have the lenght, we calculate the distance between two points
        !Print *, 'After first call of ExtraPol_baker in PathCreate.f90, s=', s

        if (s.LE.1e-6) THEN
           IF (OptGeom.LE.0) THEN
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
              WRITE(*,*) "ERROR !!! Stop !"
              STOP
           END IF
        ELSE

           dist=s/Real(NGeomf-1)
           !if (debug) WRITE(*,*) 'Dbg PathCreate ZMAT s, dist:',s,dist

           ! we call it again to obtain the geometries !
           Call ExtraPol_baker(s,dist,x0,y0,z0,xgeom,Coef,XgeomF)
           WRITE(*,*) "At the end of Baker case in the redistribution of the points, in PathCreate.f90, L966"
        END IF

        DO IGeom=1,NGeomF
           WRITE(*,*) "UMatF for IGeom=",IGeom
           DO I=1,3*Nat-6
              WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
           END DO
        END DO

     CASE ('MIXED')
        !           WRITE(*,*) "IOIOIOIOIOIOIOI"
        ! We call the extrapolation once to get the path length
        Dist=Inf
        Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef)

        if (s.LE.1e-6) THEN
           IF (OptGeom.LE.0) THEN
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
              WRITE(*,*) "ERROR !!! Stop !"
              STOP
           END IF
        ELSE

           ! we have the length, we calculate the distance between two points
           dist=s/Real(NGeomf-1)
           if (debug) WRITE(*,*) 'Dbg PathCreate Mixed s, dist:',s,dist

           ! we call it again to obtain the geometries !
           Call ExtraPol_Mixed(s,dist,x0,y0,z0,xgeom,Coef)

           IF (debug) THEN
              DO I=1,NGeomF
                 WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord)
              END DO
           END IF
        END IF
     CASE ('CART')
        ! We call the extrapolation once to get the path length
        Dist=Inf
        Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)

        if (s.LE.1e-6) THEN
           IF (OptGeom.LE.0) THEN
              WRITE(*,*) "Length of the path is 0 and you are not doing Geometry Optimization."
              WRITE(*,*) "ERROR !!! Stop !"
              STOP
           END IF
        ELSE
           ! we have the lenght, we calculate the distance between two points
           dist=s/Real(NGeomf-1)
           if (debug) WRITE(*,*) 'Debug s, dist:',s,dist

           ! we call it again to obtain the geometries !
           Call ExtraPol_cart(iopt,s,dist,x0,y0,z0,xgeom,Coef)
        END IF
     END SELECT
  ELSE
     ! PFL 3 Jan 2008 -> 
     ! Tangents are not recalculated if the points are not reparameterized
     ! PFL 6 Apr 2008 -> 
     ! Unless user has asked to reparameterized them !
     IF (mod(Iopt,iReparamT).EQ.0) THEN
        XGeomF=XGeom
        NGeomS=NGeomI
        Call CalcTangent(XgeomF,NGeomS,xgeom,Coef)
        !
     END IF
     ! <- PFL 3 Jan 2008 / 6 Apr 2008
  END IF ! if (mod(Iop,iReparam).EQ.0)
  if (debug) WRITE(*,*) "PathCreate.f90, L885"

  ! PFL 22.Nov.2007 ->
  ! We do not write tangents for baker coordinates because they are hard to 
  ! interpret for a chemist.

  IF (debug.AND.(COORD.NE."BAKER")) THEN
     ! we write the tangents :-)
     DO I=1,NGeomF
        SELECT CASE(Coord)
        CASE ('ZMAT','MIXED')
           IntCoordTmp=IntTangent(I,:)
           WRITE(*,'(1X,A,12(1X,F10.5))') "IntCoordTmp:",IntCoordTmp(1:NCoord)
           WRITE(*,'(1X,A,12(1X,F10.5))') "IntTangent:",IntTangent(I,1:NCoord)
        CASE ('CART','HYBRID')
           IntCoordTmp=XyzTangent(I,:)
        CASE DEFAULT
           WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in PathCreate. STOP"
           STOP
        END SELECT
        WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF
        Call PrintGeom(Title,Nat,Ncoord,IntCoordTmp,Coord,6,Atome,Order,OrderInv,IndZmat) 
     END DO
  END IF


  IF (DEBUG.AND.(COORD.NE.'CART')) THEN
     Idx=min(12,NCoord)
     WRITE(*,*) "DBG PathCreate. IntCoordF(IGeom,:)"
     DO I=1,NGeomF
        WRITE(*,'("Geom:",I5,12(1X,F10.4))') I,IntCoordF(I,1:Idx)
     END DO
  END IF

  DEALLOCATE(Coef,val_zmat,val_zmatTMP)
  DEALLOCATE(X0,Y0,Z0,Indzmat0)
  DEALLOCATE(X,Y,Z,Xgeom,XgeomF)
  DEALLOCATE(XyzTmp,XyzTmp2)

  if (debug) Call header("PathCreate over")

END SUBROUTINE PathCreate

