Chimie Théorique » Opt'n Path

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

!

! This subroutine analyzes the geometries along a path

! and prints it to FileUnit

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

SUBROUTINE AnaPath(IOpt,FileUnit)

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

!  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,  NGeomF, XyzGeomF, SGeom,MassAt, &

       IReparam,NbVar, FormAna, Energies,Order, &

       Renum, PrintGeomFactor, &

       Align, NFroz,  FrozAtoms, AnaFile,NbVar

  IMPLICIT NONE

  ! Input

  ! Iteration number

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

  ! Unit file to print

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

  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

  LOGICAL :: Debug, CalcDist

  INTEGER(KINT) :: I,J,K,Iat

  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:),GeomCartPrec(:,:) ! Nat,3

  REAL(KREAL), ALLOCATABLE :: xRef(:),yRef(:),zRef(:) ! Nat

  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar

  REAL(KREAL) ::  ds

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

  CHARACTER(SCHARS) :: Plot="plot"

  debug=valid("AnaPath")

  if (debug) Call header("Entering AnaPath")

  !  if (debug) WRITE(*,*) "DBG AnaPaht PrintGeomFactor:",PrintGeomFactor

  ALLOCATE(GeomCart(Nat,3))

  If (NbVar>0) THEN

     ALLOCATE(Values(NbVar))

  END IF

  CalcDist=.FALSE.

  WRITE(FileUnit,'("# Iteration ",I5)') Iopt

  DO J=1,Nat

     Iat=J

     IF (renum) Iat=Order(J)

     GeomCart(J,:)=XyzGeomF(1,:,Iat)

  END DO

  ! Do we know the curvilinear values ?

  IF (MOD(IOpt,IReparam)/=0) THEN

     CalcDist=.TRUE.

     SGeom=0.

     ALLOCATE(xRef(Nat),yRef(Nat),zref(Nat),GeomCartPrec(Nat,3))

     xRef=GeomCart(:,1)

     yRef=GeomCart(:,2)

     zRef=GeomCart(:,3)

     GeomCartPrec=GeomCart

  END IF

  if (debug) THEN

     WRITE(*,*) "AnaPath  - GeomCart - I=",1

     DO K=1,Nat

        WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)

     END DO

     WRITE(*,*) "AnaPath  - XYzGeomF(I,:,:) - I=",1

     DO K=1,Nat

        WRITE(*,'(1X,I5,3(1X,F15.8))') K,XyzGeomF(1,:,K)

     END DO

  END IF

  If (NbVar>0) THEN

     Call AnalyzeGeom(GeomCart,Values)

     WRITE(FileUnit,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)

     if (debug) WRITE(*,FormAna) SGeom(1),Values*PrintGeomFactor,0.d0,Energies(1)

  ELSE

     WRITE(FileUnit,FormAna) SGeom(1),0.d0,Energies(1)

     if (debug) WRITE(*,FormAna) SGeom(1),0.d0,Energies(1)

  END IF

  DO I=2,NGeomF

     DO J=1,Nat

        Iat=J

        IF (renum) Iat=Order(J)

        GeomCart(J,:)=XyzGeomF(I,:,Iat)

     END DO

     If (CalcDist) THEN

        ! PFL 2013 Feb 26

        ! For now, I do _NOT_ align the geometries

        if (debug) THEN

           WRITE(*,*) "AnaPath  - GeomCart avant align - I=",I

           DO K=1,Nat

              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)

           END DO

        END IF

        IF ((Nat.GE.4).OR.Align) THEN

           if (debug) WRITE(*,*) "DBG AnaPath: Aglin 2 geoms"

           ! If we have frozen atoms we align only those ones.

           if (NFroz.GT.0) THEN

              Call AlignPartial(Nat,xRef,yRef,zRef,                 &

                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),       &

                   FrozAtoms,MRot)

           ELSE 

              Call  CalcRmsd(Nat, xRef,yRef,zRef,               &

                   GeomCart(1,1),GeomCart(1,2),GeomCart(1,3),   &

                   MRot,rmsd,.TRUE.,.TRUE.)

           END IF

           ! <- PFL 24 Nov 2008

        END IF

        if (debug) WRITE(*,*) "Mass:",MassAt(1:Nat)

        ds=0.

        DO J=1,Nat

           Iat=J

           !        IF (renum) Iat=Order(J)

           if (debug) WRITE(*,*) "Mass(J):",J,MassAt(Iat)

           do K=1,3

              ds=ds+MassAt(Iat)*(GeomCartPrec(J,K)-GeomCart(J,K))**2

           END DO

           if (debug) WRITE(*,*) "DBG DS:",J,ds

        END DO

        ds=sqrt(ds) 

        if (debug) THEN

           WRITE(*,*) "AnaPath  - GeomCartPrec - I=",I

           DO K=1,Nat

              WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCartPrec(K,1:3)

           END DO

           WRITE(*,*) "DBG Anapaht, ds=",ds

        END IF

        SGeom(I)=SGeom(I-1)+ds

        GeomCartPrec=GeomCart

        xRef=GeomCart(:,1)

        yRef=GeomCart(:,2)

        zRef=GeomCart(:,3)

     END IF

     if (debug) THEN

        WRITE(*,*) "AnaPath  - GeomCart - I=",I

        DO K=1,Nat

           WRITE(*,'(1X,I5,3(1X,F15.8))') K,GeomCart(K,1:3)

        END DO

     END IF

     If (NbVar>0) THEN

        Call AnalyzeGeom(GeomCart,Values)

        WRITE(FileUnit,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)

        IF (Debug) WRITE(*,FormAna) SGeom(i),Values*PrintGeomFactor,(Energies(i)-Energies(1))*ConvE,Energies(i)

     ELSE

        WRITE(FileUnit,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)

        If (debug) WRITE(*,FormAna) SGeom(i),(Energies(i)-Energies(1))*ConvE,Energies(i)

     END IF

  END DO

  WRITE(FileUnit,*) ""

  WRITE(FileUnit,*) ""

  DeAllocate(GeomCart)

  If (NbVar>0) DeAllocate(Values)

  If (CalcDist) Deallocate(GeomCartPrec,xRef,yRef,zRef)

  Write(IoGplot,'(A,I5,A)') TRIM(Plot) // ' "' // TRIM(AnaFile) // '"  i ',Iopt," u xcol:Ecol w lp "

  Plot="replot"

  if (debug) Call header("AnaPath over")

END SUBROUTINE AnaPath