SUBROUTINE Extrapol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)

  ! This subroutine constructs the path, and if dist<>Infinity, it samples
  ! the path to obtain geometries.
  ! Basically, you call it twice: i) dist=infinity, it will calculate the length of the path
  ! ii) dist finite, it will give you the images you want along the path.
  !
  ! For now, it gives equidistant geometries
  !

  ! Default value of FRot=.TRUE.
  ! Default value of NMaxPtPath = 1000
  ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)
  ! XyzGeomF(:,:,:) ! (NGeomF,3,Nat)
  ! Default value of FAlign=.TRUE.
  use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, &
       IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, &
       Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms
  ! IndZmat(Nat,5)

  use Io_module


  IMPLICIT NONE


  REAL(KREAL), INTENT(OUT) :: s
  ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.
  REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)
  REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)
  ! Iopt: Number of the cycles for the optimization
  INTEGER(KINT), INTENT(IN) :: Iopt
  REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF)

  INTEGER(KINT) :: IdxGeom, I, J, K, Idx
  REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v
  REAL(KREAL) :: a_val, d

  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)
  REAL(KREAL), ALLOCATABLE :: ValZmat(:,:),DerInt(:) ! (Nat,3)


  LOGICAL ::  debug, print,printspline
  LOGICAL, EXTERNAL :: valid

  INTEGER(KINT) :: NSpline
  CHARACTER(LCHARS) :: FileSpline,TmpChar


  !	LOGICAL :: FAlign=.FALSE., FRot=.FALSE.

  ! We will calculate the length of the path, in MW coordinates...
  ! this is done is a stupid way: we interpolate the zmatrix values,
  ! convert them into cartesian, weight the cartesian
  ! and calculate the evolution of the distance ! 
  ! We have to follow the same procedure for every geometry
  ! so even for the first one, we have to convert it from zmat
  ! to cartesian !

  debug=valid("pathcreate")
  print=valid("printgeom")
  printspline=(valid("printspline").AND.(dist<=1e30))

  if (debug) Call Header("Entering Extrapol_int")

! We want 100 points along the interpolating path
  NSpline=int(NMaxPtPath/100)
  if (printspline) THEN
     WRITE(TmpChar,'(I5)') Iopt
     FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))
     OPEN(IOTMP,FILE=FileSpline)

  END IF

  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),valzmat(Nat,3),DerInt(NCoord))
  IdxGeom=1

  XYZTmp2(1,1)=0.
  XYZTmp2(1,2)=0.
  XYZTmp2(1,3)=0.
  XYZTmp2(2,2)=0.
  XYZTmp2(2,3)=0.
  XYZTmp2(3,3)=0.

  valzmat=0.
  valzmat(2,1)=IntCoordI(1,1)
  valzmat(3,1)=IntCoordI(1,2)
  valzmat(3,2)=IntCoordI(1,3)*180./Pi
  Idx=4
  DO I=4,Nat
     ValZmat(I,1)=IntCoordI(1,Idx)
     Idx=Idx+1
     DO J=2,3
        ValZmat(I,J)=IntCoordI(1,Idx)*180./Pi
        Idx=Idx+1
     END DO
  END DO

  XyzTmp2(2,1)=valzmat(2,1)
  d=valzmat(3,1)
  a_val=valzmat(3,2)/180.*Pi
  !              write(*,*) "aval,pi",a_val,valzmat(3,2),pi
  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
  !              i=1
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)
  !              i=2
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)
  !              i=3
  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)

  DO i=4,Nat
     call ConvertZmat_cart(i,IndZmat,valzmat,                &
          XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
     !                  WRITE(*,*) 'TOTOZma:',i,IndZmat(I,1),            &
     !                        (IndZmat(I,J+1),valzmat(I,J),J=1,3)
     !                   WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)
  END DO

  ! We align this geometry with the original one
  ! PFL 17/July/2006: only if we have more than 4 atoms.
  IF (Nat.GE.4) THEN
     ! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...),
     ! which is called in the CalcRmsd(...).
     ! PFL 24 Nov 2008 ->
     ! If we have frozen atoms we align only those ones.
     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

  XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
  IntCoordF(IdxGeom,:)=IntCoordI(1,:)

  ! We calculate the first derivatives
  u=0.d0
  DO Idx=1,NCoord
     call splintder(u,v,DerInt(iDX),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
  END DO
  IntTangent(IdxGeom,:)=DerInt

  if (print.AND.(Dist.LE.1e20)) THEN
     WRITE(IOOUT,'(1X,I5)') Nat
     WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom
     DO i=1,Nat
        If (Renum) THEN
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &
                (XyzTmp2(Order(I),J),J=1,3)
        ELSE
           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),      &
                (XyzTmp2(I,J),J=1,3)
        END IF
     END DO
  END IF

  ! We initialize the first geometry
  ! PFL 26 Nov 2008 ->
  ! Now, I copy XyzTmp2 into XyzTmp because XyzTmp2 has been aligned
  ! with the reference geometry

  XyzTmp=XyzTmp2

  !   XyzTmp(1,1)=0.
  !   XyzTmp(1,2)=0.
  !   XyzTmp(1,3)=0.
  !   XyzTmp(2,2)=0.
  !   XyzTmp(2,3)=0.     
  !   XyzTmp(3,3)=0.

  !   XyzTmp(2,1)=valzmat(2,1)
  !   d=valzmat(3,1)
  !   a_val=valzmat(3,2)/180.*Pi

  !   if (Nat.GE.3) THEN
  !      if (IndZmat(3,2).EQ.1)  THEN
  !         XyzTmp(3,1)=XYzTmp(1,1)+d*cos(a_val)
  !      ELSE
  !         XyzTmp(3,1)=XyzTmp(2,1)-d*cos(a_val)
  !      ENDIF
  !      XyzTmp(3,2)=d*sin(a_val)
  !   ENDIF

  !   DO i=4,Nat
  !      call ConvertZmat_cart(i,IndZmat,valzmat,        &
  !           XYzTMP(1,1), XYzTMP(1,2),XYzTMP(1,3))
  !   END DO

  ! <- PFL 26 Nov 2008

  s=0.
  if (printspline) THEN
     SELECT CASE (Nat)
     CASE(2)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
     CASE (3)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
     CASE(4:)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
     END SELECT
  END IF
  valzmat=0.d0

  DO K=1,NMaxPtPath
     u=real(K)/NMaxPtPath*(NGeomI-1.)

     XYZTmp2(1,1)=0.
     XYZTmp2(1,2)=0.
     XYZTmp2(1,3)=0.
     XYZTmp2(2,2)=0.
     XYZTmp2(2,3)=0.
     XYZTmp2(3,3)=0.

     ! We generate the interpolated Zmatrix
     call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))
     valzmat(2,1)=v
     call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))
     valzmat(3,1)=v
     call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))
     valzmat(3,2)=v*180./Pi
     Idx=4
     DO I=4,Nat
        call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
        valzmat(I,1)=v
        Idx=Idx+1
        DO J=2,3
           call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
           valzmat(I,J)=v*180./pi
           Idx=Idx+1
        END DO
     END DO ! matches DO I=4,Nat

     ! We convert it into Cartesian coordinates
     XyzTmp2(2,1)=valzmat(2,1)
     d=valzmat(3,1)
     a_val=valzmat(3,2)/180.*Pi
     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

     DO I=4,Nat
        call ConvertZmat_cart(I,IndZmat,valzmat,       &
             XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
     END DO

     IF (Nat.GE.4) THEN
        ! PFL 24 Nov 2008 ->
        ! If we have frozen atoms we align only those ones.
        if (NFroz.GT.0) THEN
           Call AlignPartial(Nat,x0,y0,z0,                     &
                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
                FrozAtoms,MRot)
        ELSE
           ! PFL 26 Nov 2008!
           ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here
           ! we align on the previous geom... what is the best ? Is there a difference ?
           Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &
                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
                MRot,rmsd,.TRUE.,.TRUE.)
        END IF
        ! <- PFL 24 Nov 2008
     END IF

     ds=0.
     DO I=1,Nat
        DO J=1,3
           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2
           XYZTmp(I,J)=XyzTMP2(I,J)
        END DO
     END DO


     s=s+sqrt(ds)

     if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN
        SELECT CASE (Nat)
        CASE(2)
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
        CASE (3)
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
        CASE(4:)
           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
        END SELECT
     END IF

     !         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist
     if (s>=dist) THEN
        if (debug) THEN
           WRITE(*,*) "DBG Zmat",s
           WRITE(*,'(1X,I5)') IndZmat(1,1)
           WRITE(*,'(1X,I5,I5,1X,F10.4)') IndZmat(2,1),IndZmat(2,2),valzmat(2,1)
           WRITE(*,'(1X,I5,2(I5,1X,F10.4))') IndZmat(3,1),IndZmat(3,2),valzmat(3,1)  &
                ,IndZmat(3,3),valzmat(3,2)
           DO I=4,Nat
              WRITE(*,'(1X,I5,3(I5,1X,F10.4))') IndZmat(I,1),IndZmat(I,2),valzmat(I,1)  &
                   ,IndZmat(I,3),valzmat(I,2),IndZmat(I,4),valzmat(I,3)
           END DO
        END IF ! matches if (debug) THEN

        s=s-dist
        IdxGeom=IdxGeom+1
        SGeom(IdxGeom)=s+IdxGeom*dist
        XgeomF(IdxGeom)=u
        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
        IntCoordF(IdxGeom,1)=valzmat(2,1)
        IntCoordF(IdxGeom,2)=valzmat(3,1)
        IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi
        Idx=4
        DO I=4,Nat
           IntCoordF(IdxGeom,Idx)=valzmat(I,1)
           IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi
           IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi
           Idx=Idx+3
        END DO
        IntTangent(IdxGeom,:)=DerInt

        if (print) THEN
           WRITE(IOOUT,'(1X,I5)') Nat
           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
           ! PFL 17/July/2006: only if we have more than 4 atoms.
           IF (Nat.GE.4) THEN
              ! PFL 24 Nov 2008 ->
              ! If we have frozen atoms we align only those ones.
              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

           DO I=1,Nat
              If (Renum) THEN
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),   &
                      (XyzTmp2(Order(I),J),J=1,3)
              ELSE
                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),    &
                      (XyzTmp2(I,J),J=1,3)
              END IF
           END DO

        END IF ! matches if (print) THEN
     END IF ! matches if (s>=dist) THEN
  END DO ! matches DO K=1,NMaxPtPath, Is it correct??


  if (printspline) THEN
     SELECT CASE (Nat)
     CASE(2)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)
     CASE (3)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)
     CASE(4:)
        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)
     END SELECT
  END IF

  if (s>=0.9*dist) THEN
     if (debug) WRITE(*,*) "DBG Extrapol_int L383: Adding last geom"
     write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist
!     u=xgeom(NGeomI)
     s=s-dist


!      ! We generate the interpolated Zmatrix
!      call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))
!      valzmat(2,1)=v
!      call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))
!      valzmat(3,1)=v
!      call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))
!      valzmat(3,2)=v*180./Pi
!      Idx=4
!      DO I=4,Nat
!         call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
!         valzmat(I,1)=v
!         Idx=Idx+1
!         DO J=2,3
!            call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))
!            valzmat(I,J)=v*180./pi
!            Idx=Idx+1
!         END DO
!      END DO ! matches DO I=4,Nat

!      ! We convert it into Cartesian coordinates
!      XyzTmp2(2,1)=valzmat(2,1)
!      d=valzmat(3,1)
!      a_val=valzmat(3,2)/180.*Pi
!      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

!      DO I=4,Nat
!         call ConvertZmat_cart(I,IndZmat,valzmat,       &
!              XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))
!      END DO

!      IF (Nat.GE.4) THEN
!         ! PFL 24 Nov 2008 ->
!         ! If we have frozen atoms we align only those ones.
!         if (NFroz.GT.0) THEN
!            Call AlignPartial(Nat,x0,y0,z0,                     &
!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
!                 FrozAtoms,MRot)
!         ELSE
!            ! PFL 26 Nov 2008!
!            ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here
!            ! we align on the previous geom... what is the best ? Is there a difference ?
!            Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &
!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &
!                 MRot,rmsd,.TRUE.,.TRUE.)
!         END IF
!         ! <- PFL 24 Nov 2008
!      END IF

     IdxGeom=IdxGeom+1



     IF (IdxGeom.GT.NGeomF) THEN
        WRITE(IOOUT,*) "!!! ERROR in Extrapol_int !!!!"
        WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"
        WRITE(*,*) "** PathCreate ***"
        WRITE(*,*) "Distribution of points along the path is wrong."
        WRITE(*,*) "Increase value of NMaxPtPath in the input file"
        WRITE(*,*) "Present value is:", NMaxPtPath
        STOP
     END IF

     SGeom(IdxGeom)=s+IdxGeom*dist
     XgeomF(IdxGeom)=min(u,NGeomI-1.d0)
     XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))
     !  XyzGeomF(IdxGeom,:,:)=XyzTmp2(:,:)
     IntCoordF(IdxGeom,1)=valzmat(2,1)
     IntCoordF(IdxGeom,2)=valzmat(3,1)
     IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi
     Idx=4
     DO I=4,Nat
        IntCoordF(IdxGeom,Idx)=valzmat(I,1)
        IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi
        IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi
        Idx=Idx+3
     END DO
     IntTangent(IdxGeom,:)=DerInt     

     if (print) THEN
        WRITE(IOOUT,'(1X,I5)') Nat
        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K
        ! PFL 17/July/2006: only if we have more than 4 atoms.
        IF (Nat.GE.4) THEN
           ! PFL 24 Nov 2008 ->
           ! If we have frozen atoms we align only those ones.
           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

        DO I=1,Nat
           If (Renum) THEN
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &
                   (XyzTmp2(Order(I),J),J=1,3)
           ELSE
              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),     &
                   (XyzTmp2(I,J),J=1,3)
           END IF
        END DO ! matches DO I=1,Nat
     END IF ! matches if (print) THEN
  END IF ! matches if (s>=0.9*dist) THEN

  if (debug) WRITE(*,*) 's final =',s

  DEALLOCATE(XyzTmp,XyzTmp2,valzmat)

  if (printspline) CLOSE(IOTMP)

  if (debug) Call Header("Extrapol_int Over")


END SUBROUTINE EXTRAPOL_INT
