Chimie Théorique » OpenPath

!================================================================

!

! 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

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

        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