Chimie Théorique » OpenPath

!This programs reads a set of geometries, and print the tangent

! to the path they form

SUBROUTINE CalcTangent(XgeomF,NGeomS,XGeom,Coef)

  use Path_module, only: IntCoordI, Nat, Coord, XyzGeomI, Atome, Order, &

                         OrderInv, IndZmat, NMaxPtPath, IntTangent, XyzTangent, &

             NGeomF, NCoord

  use Io_module, ONly : IoTmp,PathName, Kint,KReal

  use VarTypes

  IMPLICIT NONE

  ! NGeomS: number of geometries.

  INTEGER(KINT), INTENT(IN) :: NGeomS

  REAL(KREAL), INTENT(IN) :: XGeomF(NGeomF),Xgeom(NGeomS),Coef(NGeomS,NCoord)

  REAL(KREAL), ALLOCATABLE :: Path(:,:) ! NGeomS,NCoord

  REAL(KREAL), ALLOCATABLE :: PathTan(:,:) ! NGeomF,NCoord

  REAL(KREAL), ALLOCATABLE :: DbgTgt(:),GeomTmp(:) ! NCoord

  INTEGER(KINT) :: i, j, n3at, Idx, iat,igeom,k

  INTEGER(KINT), SAVE :: icall=-1

  INTEGER(KINT) :: NMaxPt

  REAL(KREAL), PARAMETER :: Inf=1e32, h=0.001d0

  CHARACTER(120) :: Line,File

  CHARACTER(LCHARS) :: TITLE

  LOGICAL :: Debug=.False.

  REAL(KREAL) :: utmp,der

  LOGICAL, EXTERNAL :: valid

  debug=valid('calctangent')

  if (debug) WRITE(*,*) "======================================== Entering CalcTangen ============================"

  n3at=3*nat

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

 icall=icall+2

 write(Line,'(I5)') icall

 File=TRIM(adjustl(PathName)) // '_dbgtgt_' // TRIM(adjustl(Line)) // '.dat'

 OPEN(IOTMP,File=File)

  ! We construct path, depending on COORD:

  ! Path(NGeomS,Ncoord), NGeoms is equal to NGeomI

  ! IntCoordI(:,:) ! (NGeomI,3*Nat-6)

 SELECT CASE(Coord)

    CASE ('ZMAT','MIXED')

     Path(:,:)=IntCoordI(:,:)

  CASE ('CART','HYBRID')

     Path=reshape(XyzGeomI,(/NGeomS,NCoord/))

  CASE ('BAKER')

   Path(:,:)=IntCoordI(:,:)

  CASE DEFAULT

     WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in CalcTangent. STOP"

     STOP

  END SELECT

  if (debug) THEN

     WRITE(*,*) 'DBG CalcTangent. Geometries from Path'

     DO I=1,NGeomS

        GeomTmp=Path(I,1:NCoord)

        WRITE(Title,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS

        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)

!         IF (COORD.EQ.'ZMAT') THEN

!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomS

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &

!                 Nom(Atome(OrderInv(indzmat(2,2)))),Path(I,1)

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &

!                 Nom(Atome(OrderInv(indzmat(3,2)))),Path(I,2), &

!                 Nom(Atome(OrderInv(indzmat(3,3)))),Path(I,3)

!            Idx=4

!            DO Iat=4,Nat

!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

!                    Nom(Atome(OrderInv(indzmat(iat,2)))),Path(I,Idx), &

!                    Nom(Atome(OrderInv(indzmat(iat,3)))),Path(I,Idx+1), &

!                    Nom(Atome(OrderInv(indzmat(iat,4)))),Path(I,Idx+2)

!               Idx=Idx+3

!            END DO

!         ELSE

!            WRITE(*,*) Nat

!            WRITE(*,"('Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF

!            DO Iat=1,Nat

!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

!                    Path(I,3*Iat-2:3*Iat)

!            END DO

!         END IF

     END DO

  END IF

!Debug info

 IF (debug.AND.(COORD.EQ.'ZMAT')) THEN

     ALLOCATE(DbgTgt(NCoord))

     WRITE(*,*) 'DBG CalcTangent. Path written in ' // TRIM(File)

     WRITE(*,'(A,12(1X,F10.5))') "DBG CalcTangent, XgeomF:",XGeomF

     NMaxPt=min(200,NMaxPtPath)-1

    DO K=0,NMaxPt

       utmp=K*Xgeom(NGeomS)/NMaxPt

       DO J=1,NCoord

       Call Splint(utmp,DbgTgt(J),NGeomS,Xgeom,Path(1,J),Coef(1,J))

       END DO

       WRITE(IOTMP,'(12(1X,F10.5))') utmp,DbgTgt

    END DO

    DEALLOCATE(DbgTgt)

 END IF

!We construct the derivatives

DO I=1,NGeomF

  DO J=1,NCoord

        Call SplinDer(XgeomF(I),PathTan(I,J),NGeomS,Xgeom,Path(1,J),Coef(1,J))

  END DO

END DO

 SELECT CASE(Coord)

    CASE ('ZMAT','MIXED')

     IntTangent=PathTan

  CASE ('CART','HYBRID')

     XyzTangent=Pathtan

  END SELECT

  if (debug) THEN

     WRITE(*,*) 'DBG CalcTangent. Tangents'

     DO I=1,NGeomF

        GeomTmp=PathTan(I,1:NCoord)

        WRITE(Title,"('Tangents for geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF

        Call PrintGeom(Title,Nat,Ncoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)

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

!         IF (COORD.EQ.'ZMAT') THEN

!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &

!                 Nom(Atome(OrderInv(indzmat(2,2)))),IntTangent(I,1)

!            WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &

!                 Nom(Atome(OrderInv(indzmat(3,2)))),IntTangent(I,2), &

!                 Nom(Atome(OrderInv(indzmat(3,3)))),IntTangent(I,3)

!            Idx=4

!            DO Iat=4,Nat

!               WRITE(*,'(1X,A5,3(1X,A4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

!                    Nom(Atome(OrderInv(indzmat(iat,2)))),IntTangent(I,Idx), &

!                    Nom(Atome(OrderInv(indzmat(iat,3)))),IntTangent(I,Idx+1), &

!                    Nom(Atome(OrderInv(indzmat(iat,4)))),IntTangent(I,Idx+2)

!               Idx=Idx+3

!            END DO

!         ELSE

!            WRITE(*,*) Nat

!            WRITE(*,"('Tangent for Geometry ',I3,'/',I3,' for iteration ',I3)") I,NgeomF

!            DO Iat=1,Nat

!               WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

!                    XyzTangent(I,3*Iat-2:3*Iat)

!            END DO

!         END IF

     END DO

  END IF

DeALLOCATE(Path) 

 CLOSE(IOTMP)

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

END SUBROUTINE CalcTangent