Chimie Théorique » Opt'n Path

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

! to the path they form

SUBROUTINE CalcTangent(XgeomF,NGeomS,XGeom,Coef)

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

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

  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