Chimie Théorique » Opt'n Path

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

!

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

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

  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,  Nom, &

       ISpline, IntTangent,  NPrim,Xprimitive_t, OptGeom,  &

       UMatF, UMat_local, XyzTangent, Linear, Align, FrozAtoms,AtName

  ! 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

  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

  INTEGER(KINT) :: NGeomS

  REAL(KREAL) :: Rmsd, MRot(3, 3), Delta

  REAL(KREAL), PARAMETER ::  TwoPi=360.d0

  REAL(KREAL)  ::s, dist, a_val, d

  INTEGER(KINT) :: I, J, igeom, iat

  INTEGER(KINT) :: Itmp, Jat

  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

    SUBROUTINE die(routine, msg, file, line, unit)

      Use VarTypes

      Use io_module

      implicit none

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

    SUBROUTINE Warning(routine, msg, file, line, unit)

      Use VarTypes

      Use io_module

      implicit none

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE Warning

  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

     if (Print) THEN

        WRITE(*,*) "PathCreate - L121 - Initial geometries "

        DO J=1,Nat

           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,J)

        END DO

     END IF

  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 (debug) THEN

              WRITE(*,*) 'Debug PathCreat Line190'

              WRITE(*,*) 'Igeom,NbFrag, NbFragRef=',IGeom,NbFrag,NFragRef

           END IF

           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 ((COORD.EQ."CART").AND.(IGeomREF.LE.0)) THEN

        IGeomRef=1

        CALL Warning('PathCreate L209','IGeomRef<=0',UNIT=IOOUT)

     END IF

     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

     if (Print) THEN

        WRITE(*,*) "PathCreate - L631 - geometries "

        DO I=1,NGeomI

           DO J=1,Nat

              If (renum) THEN

                 Iat=Order(J)

                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)

              ELSE

                 Iat=OrderInv(J)

                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)

              END IF

           END DO

        END DO

     END IF

  ! 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

           WRITE(*,*) "PathCreate L786: NCart,IBeg",NCart,IBeg

           DO J=IBeg,Nat

              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

     if (Print) THEN

        WRITE(*,*) "PathCreate - L811 - geometries "

        DO I=1,NGeomI

           DO J=1,Nat

              If (renum) THEN

                 Iat=Order(J)

                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(J)),XyzGeomI(I,1:3,Iat)

              ELSE

                 Iat=OrderInv(J)

                 WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),XyzGeomI(I,1:3,J)

              END IF

           END DO

        END DO

     END IF

  ! 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