Chimie Théorique » OpenPath

SUBROUTINE Extrapol_int(iopt,s,dist,x0,y0,z0,xgeom,Coef,XgeomF)

  ! This subroutine constructs the path, and if dist<>Infinity, it samples

  ! the path to obtain geometries.

  ! Basically, you call it twice: i) dist=infinity, it will calculate the length of the path

  ! ii) dist finite, it will give you the images you want along the path.

  !

  ! For now, it gives equidistant geometries

  !

  ! Default value of FRot=.TRUE.

  ! Default value of NMaxPtPath = 1000

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

  ! XyzGeomF(:,:,:) ! (NGeomF,3,Nat)

  ! Default value of FAlign=.TRUE.

  use Path_module, only : NMaxPtPath, IntCoordI, Pi, IndZmat, XyzGeomF, &

       IntCoordF, IntTangent, Renum, Nom, Order, MassAt, SGeom, XyzGeomI, &

       Atome, Nat, NGeomI, NCoord, NGeomF, OrderInv,NFroz,FrozAtoms

  ! IndZmat(Nat,5)

  use Io_module

  IMPLICIT NONE

  REAL(KREAL), INTENT(OUT) :: s

  ! X0(Nat),Y0(Nat),Z0(Nat): reference geometry.

  REAL(KREAL), INTENT(IN) :: dist,X0(Nat),Y0(Nat),Z0(Nat)

  REAL(KREAL), INTENT(IN) :: Xgeom(NGeomI),Coef(NGeomI,NCoord)

  ! Iopt: Number of the cycles for the optimization

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

  REAL(KREAL), INTENT(OUT) :: XGeomF(NGeomF)

  INTEGER(KINT) :: IdxGeom, I, J, K, Idx

  REAL(KREAL) :: Rmsd,MRot(3,3), ds, u, v

  REAL(KREAL) :: a_val, d

  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:), XyzTmp2(:,:) ! (Nat,3)

  REAL(KREAL), ALLOCATABLE :: ValZmat(:,:),DerInt(:) ! (Nat,3)

  LOGICAL ::  debug, print,printspline

  LOGICAL, EXTERNAL :: valid

  INTEGER(KINT) :: NSpline

  CHARACTER(LCHARS) :: FileSpline,TmpChar

  !	LOGICAL :: FAlign=.FALSE., FRot=.FALSE.

  ! We will calculate the length of the path, in MW coordinates...

  ! this is done is a stupid way: we interpolate the zmatrix values,

  ! convert them into cartesian, weight the cartesian

  ! and calculate the evolution of the distance !

  ! We have to follow the same procedure for every geometry

  ! so even for the first one, we have to convert it from zmat

  ! to cartesian !

  debug=valid("pathcreate")

  print=valid("printgeom")

  printspline=(valid("printspline").AND.(dist<=1e30))

  if (debug) Call Header("Entering Extrapol_int")

! We want 100 points along the interpolating path

  NSpline=int(NMaxPtPath/100)

  if (printspline) THEN

     WRITE(TmpChar,'(I5)') Iopt

     FileSpline=Trim(adjustL(PathName)) // '_spline.' // AdjustL(TRIM(TmpChar))

     OPEN(IOTMP,FILE=FileSpline)

  END IF

  ALLOCATE(XyzTmp(Nat,3),XyzTmp2(Nat,3),valzmat(Nat,3),DerInt(NCoord))

  IdxGeom=1

  XYZTmp2(1,1)=0.

  XYZTmp2(1,2)=0.

  XYZTmp2(1,3)=0.

  XYZTmp2(2,2)=0.

  XYZTmp2(2,3)=0.

  XYZTmp2(3,3)=0.

  valzmat=0.

  valzmat(2,1)=IntCoordI(1,1)

  valzmat(3,1)=IntCoordI(1,2)

  valzmat(3,2)=IntCoordI(1,3)*180./Pi

  Idx=4

  DO I=4,Nat

     ValZmat(I,1)=IntCoordI(1,Idx)

     Idx=Idx+1

     DO J=2,3

        ValZmat(I,J)=IntCoordI(1,Idx)*180./Pi

        Idx=Idx+1

     END DO

  END DO

  XyzTmp2(2,1)=valzmat(2,1)

  d=valzmat(3,1)

  a_val=valzmat(3,2)/180.*Pi

  !              write(*,*) "aval,pi",a_val,valzmat(3,2),pi

  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

  !              i=1

  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)

  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)

  !              i=2

  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)

  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)

  !              i=3

  !                WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)

  !                WRITE(*,*) 'TOTOZma:',i, (valzmat(I,J),J=1,3)

  DO i=4,Nat

     call ConvertZmat_cart(i,IndZmat,valzmat,                &

          XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))

     !                  WRITE(*,*) 'TOTOZma:',i,IndZmat(I,1),            &

     !                        (IndZmat(I,J+1),valzmat(I,J),J=1,3)

     !                   WRITE(*,*) 'TOTOCart:',i, (XYzTMP2(I,J),J=1,3)

  END DO

  ! We align this geometry with the original one

  ! PFL 17/July/2006: only if we have more than 4 atoms.

  IF (Nat.GE.4) THEN

     ! The rotation matrix MRot has INTENT(OUT) in subroutine rotation_matrix(...),

     ! which is called in the CalcRmsd(...).

     ! PFL 24 Nov 2008 ->

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

     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

  XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))

  IntCoordF(IdxGeom,:)=IntCoordI(1,:)

  ! We calculate the first derivatives

  u=0.d0

  DO Idx=1,NCoord

     call splintder(u,v,DerInt(iDX),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))

  END DO

  IntTangent(IdxGeom,:)=DerInt

  if (print.AND.(Dist.LE.1e20)) THEN

     WRITE(IOOUT,'(1X,I5)') Nat

     WRITE(IOOUT,*) "# Cartesian Coordinates for geom",IdxGeom

     DO i=1,Nat

        If (Renum) THEN

           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),      &

                (XyzTmp2(Order(I),J),J=1,3)

        ELSE

           WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),      &

                (XyzTmp2(I,J),J=1,3)

        END IF

     END DO

  END IF

  ! We initialize the first geometry

  ! PFL 26 Nov 2008 ->

  ! Now, I copy XyzTmp2 into XyzTmp because XyzTmp2 has been aligned

  ! with the reference geometry

  XyzTmp=XyzTmp2

  !   XyzTmp(1,1)=0.

  !   XyzTmp(1,2)=0.

  !   XyzTmp(1,3)=0.

  !   XyzTmp(2,2)=0.

  !   XyzTmp(2,3)=0.

  !   XyzTmp(3,3)=0.

  !   XyzTmp(2,1)=valzmat(2,1)

  !   d=valzmat(3,1)

  !   a_val=valzmat(3,2)/180.*Pi

  !   if (Nat.GE.3) THEN

  !      if (IndZmat(3,2).EQ.1)  THEN

  !         XyzTmp(3,1)=XYzTmp(1,1)+d*cos(a_val)

  !      ELSE

  !         XyzTmp(3,1)=XyzTmp(2,1)-d*cos(a_val)

  !      ENDIF

  !      XyzTmp(3,2)=d*sin(a_val)

  !   ENDIF

  !   DO i=4,Nat

  !      call ConvertZmat_cart(i,IndZmat,valzmat,        &

  !           XYzTMP(1,1), XYzTMP(1,2),XYzTMP(1,3))

  !   END DO

  ! <- PFL 26 Nov 2008

  s=0.

  if (printspline) THEN

     SELECT CASE (Nat)

     CASE(2)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)

     CASE (3)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)

     CASE(4:)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)

     END SELECT

  END IF

  valzmat=0.d0

  DO K=1,NMaxPtPath

     u=real(K)/NMaxPtPath*(NGeomI-1.)

     XYZTmp2(1,1)=0.

     XYZTmp2(1,2)=0.

     XYZTmp2(1,3)=0.

     XYZTmp2(2,2)=0.

     XYZTmp2(2,3)=0.

     XYZTmp2(3,3)=0.

     ! We generate the interpolated Zmatrix

     call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))

     valzmat(2,1)=v

     call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))

     valzmat(3,1)=v

     call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))

     valzmat(3,2)=v*180./Pi

     Idx=4

     DO I=4,Nat

        call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))

        valzmat(I,1)=v

        Idx=Idx+1

        DO J=2,3

           call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))

           valzmat(I,J)=v*180./pi

           Idx=Idx+1

        END DO

     END DO ! matches DO I=4,Nat

     ! We convert it into Cartesian coordinates

     XyzTmp2(2,1)=valzmat(2,1)

     d=valzmat(3,1)

     a_val=valzmat(3,2)/180.*Pi

     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

     DO I=4,Nat

        call ConvertZmat_cart(I,IndZmat,valzmat,       &

             XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))

     END DO

     IF (Nat.GE.4) THEN

        ! PFL 24 Nov 2008 ->

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

        if (NFroz.GT.0) THEN

           Call AlignPartial(Nat,x0,y0,z0,                     &

                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &

                FrozAtoms,MRot)

        ELSE

           ! PFL 26 Nov 2008!

           ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here

           ! we align on the previous geom... what is the best ? Is there a difference ?

           Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &

                xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &

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

        END IF

        ! <- PFL 24 Nov 2008

     END IF

     ds=0.

     DO I=1,Nat

        DO J=1,3

           ds=ds+MassAt(I)*(XYZTMp2(I,J)-XYZTmp(I,J))**2

           XYZTmp(I,J)=XyzTMP2(I,J)

        END DO

     END DO

     s=s+sqrt(ds)

     if ((printspline).AND.(MOD(K,NSpline).EQ.0)) THEN

        SELECT CASE (Nat)

        CASE(2)

           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)

        CASE (3)

           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)

        CASE(4:)

           WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)

        END SELECT

     END IF

     !         if (debug) WRITE(*,*) "Debug u,s,dist",u,s,dist

     if (s>=dist) THEN

        if (debug) THEN

           WRITE(*,*) "DBG Zmat",s

           WRITE(*,'(1X,I5)') IndZmat(1,1)

           WRITE(*,'(1X,I5,I5,1X,F10.4)') IndZmat(2,1),IndZmat(2,2),valzmat(2,1)

           WRITE(*,'(1X,I5,2(I5,1X,F10.4))') IndZmat(3,1),IndZmat(3,2),valzmat(3,1)  &

                ,IndZmat(3,3),valzmat(3,2)

           DO I=4,Nat

              WRITE(*,'(1X,I5,3(I5,1X,F10.4))') IndZmat(I,1),IndZmat(I,2),valzmat(I,1)  &

                   ,IndZmat(I,3),valzmat(I,2),IndZmat(I,4),valzmat(I,3)

           END DO

        END IF ! matches if (debug) THEN

        s=s-dist

        IdxGeom=IdxGeom+1

        SGeom(IdxGeom)=s+IdxGeom*dist

        XgeomF(IdxGeom)=u

        XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))

        IntCoordF(IdxGeom,1)=valzmat(2,1)

        IntCoordF(IdxGeom,2)=valzmat(3,1)

        IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi

        Idx=4

        DO I=4,Nat

           IntCoordF(IdxGeom,Idx)=valzmat(I,1)

           IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi

           IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi

           Idx=Idx+3

        END DO

        IntTangent(IdxGeom,:)=DerInt

        if (print) THEN

           WRITE(IOOUT,'(1X,I5)') Nat

           WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K

           ! PFL 17/July/2006: only if we have more than 4 atoms.

           IF (Nat.GE.4) THEN

              ! PFL 24 Nov 2008 ->

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

              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

           DO I=1,Nat

              If (Renum) THEN

                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),   &

                      (XyzTmp2(Order(I),J),J=1,3)

              ELSE

                 WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),    &

                      (XyzTmp2(I,J),J=1,3)

              END IF

           END DO

        END IF ! matches if (print) THEN

     END IF ! matches if (s>=dist) THEN

  END DO ! matches DO K=1,NMaxPtPath, Is it correct??

  if (printspline) THEN

     SELECT CASE (Nat)

     CASE(2)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1)

     CASE (3)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2)

     CASE(4:)

        WRITE(IOTMP,'(15(1X,F10.5))') s,valzmat(2,1),valzmat(3,1:2),valzmat(4:Nat,1:3)

     END SELECT

  END IF

  if (s>=0.9*dist) THEN

     if (debug) WRITE(*,*) "DBG Extrapol_int L383: Adding last geom"

     write(*,*) "Extrapol_int u,xgeom(NGeomI),s,dist,s-dist",u,xgeom(NGeomI),s,dist,s-dist

!     u=xgeom(NGeomI)

     s=s-dist

!      ! We generate the interpolated Zmatrix

!      call splintder(u,v,DerInt(1),NGeomI,xgeom(1),IntCoordI(1,1),Coef(1,1))

!      valzmat(2,1)=v

!      call splintder(u,v,DerInt(2),NGeomI,xgeom(1),IntCoordI(1,2),Coef(1,2))

!      valzmat(3,1)=v

!      call splintder(u,v,DerInt(3),NGeomI,xgeom(1),IntCoordI(1,3),Coef(1,3))

!      valzmat(3,2)=v*180./Pi

!      Idx=4

!      DO I=4,Nat

!         call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))

!         valzmat(I,1)=v

!         Idx=Idx+1

!         DO J=2,3

!            call splintder(u,v,DerInt(Idx),NGeomI,xgeom(1),IntCoordI(1,Idx),Coef(1,Idx))

!            valzmat(I,J)=v*180./pi

!            Idx=Idx+1

!         END DO

!      END DO ! matches DO I=4,Nat

!      ! We convert it into Cartesian coordinates

!      XyzTmp2(2,1)=valzmat(2,1)

!      d=valzmat(3,1)

!      a_val=valzmat(3,2)/180.*Pi

!      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

!      DO I=4,Nat

!         call ConvertZmat_cart(I,IndZmat,valzmat,       &

!              XYzTMP2(1,1), XYzTMP2(1,2),XYzTMP2(1,3))

!      END DO

!      IF (Nat.GE.4) THEN

!         ! PFL 24 Nov 2008 ->

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

!         if (NFroz.GT.0) THEN

!            Call AlignPartial(Nat,x0,y0,z0,                     &

!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &

!                 FrozAtoms,MRot)

!         ELSE

!            ! PFL 26 Nov 2008!

!            ! Note to myself: in Extrapol_cart we always align on x0,y0,z0 but here

!            ! we align on the previous geom... what is the best ? Is there a difference ?

!            Call  CalcRmsd(Nat, XyzTmp(1:Nat,1),XyzTmp(1:Nat,2),XyzTmp(1:Nat,3), &

!                 xyzTmp2(1,1),xyzTmp2(1,2),xyzTMP2(1,3),       &

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

!         END IF

!         ! <- PFL 24 Nov 2008

!      END IF

     IdxGeom=IdxGeom+1

     IF (IdxGeom.GT.NGeomF) THEN

        WRITE(IOOUT,*) "!!! ERROR in Extrapol_int !!!!"

        WRITE(IOOUT,*) "Too many structures. Increase NMaxPath"

        WRITE(*,*) "** PathCreate ***"

        WRITE(*,*) "Distribution of points along the path is wrong."

        WRITE(*,*) "Increase value of NMaxPtPath in the input file"

        WRITE(*,*) "Present value is:", NMaxPtPath

        STOP

     END IF

     SGeom(IdxGeom)=s+IdxGeom*dist

     XgeomF(IdxGeom)=min(u,NGeomI-1.d0)

     XyzGeomF(IdxGeom,:,:)=Reshape(XyzTmp2(:,:),(/3,Nat/),ORDER=(/2,1/))

     !  XyzGeomF(IdxGeom,:,:)=XyzTmp2(:,:)

     IntCoordF(IdxGeom,1)=valzmat(2,1)

     IntCoordF(IdxGeom,2)=valzmat(3,1)

     IntCoordF(IdxGeom,3)=valzmat(3,2)/180.*Pi

     Idx=4

     DO I=4,Nat

        IntCoordF(IdxGeom,Idx)=valzmat(I,1)

        IntCoordF(IdxGeom,Idx+1)=valzmat(I,2)/180.*Pi

        IntCoordF(IdxGeom,Idx+2)=valzmat(I,3)/180.*Pi

        Idx=Idx+3

     END DO

     IntTangent(IdxGeom,:)=DerInt

     if (print) THEN

        WRITE(IOOUT,'(1X,I5)') Nat

        WRITE(IOOUT,*) "# Cartesian coord for Geometry ",IdxGeom,K

        ! PFL 17/July/2006: only if we have more than 4 atoms.

        IF (Nat.GE.4) THEN

           ! PFL 24 Nov 2008 ->

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

           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

        DO I=1,Nat

           If (Renum) THEN

              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(I)),    &

                   (XyzTmp2(Order(I),J),J=1,3)

           ELSE

              WRITE(IOOUT,'(1X,A2,3(1X,F15.6))') Nom(Atome(OrderInv(I))),     &

                   (XyzTmp2(I,J),J=1,3)

           END IF

        END DO ! matches DO I=1,Nat

     END IF ! matches if (print) THEN

  END IF ! matches if (s>=0.9*dist) THEN

  if (debug) WRITE(*,*) 's final =',s

  DEALLOCATE(XyzTmp,XyzTmp2,valzmat)

  if (printspline) CLOSE(IOTMP)

  if (debug) Call Header("Extrapol_int Over")

END SUBROUTINE EXTRAPOL_INT