Chimie Théorique » Opt'n Path

!C  HEAT is never used, not even in call of Space(...)

!C  Geom = input parameter vector (Geometry).

!C  Grad = input gradient vector.

!C  Geom_new = New Geometry.

    SUBROUTINE Step_diis_all(NGeomF,IGeom,Step,Geom,Grad,HP,HEAT,Hess,NCoord,allocation_flag,Tangent)

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

! This routine was adapted from the public domain mopac6 diis.f

!  source file (c) 2009, Stewart Computational Chemistry.

!  <http://www.openmopac.net/Downloads/Downloads.html>

!

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

!  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, only : IoOut

    USE Path_module, only : Vfree

    IMPLICIT NONE

    INTEGER, parameter :: KINT = kind(1)

    INTEGER, parameter :: KREAL = kind(1.0d0)

    !      INCLUDE 'SIZES'

    INTEGER(KINT) :: NGeomF,IGeom

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

    REAL(KREAL) :: Geom(NCoord), Grad(NCoord)

    REAL(KREAL) :: Hess(NCoord*NCoord),Step(NCoord)

    REAL(KREAL) :: HEAT,HP

    LOGICAL :: allocation_flag

    REAL(KREAL), INTENT(INOUT) :: Tangent(Ncoord)

!************************************************************************

!*                                                                      *

!*     DIIS PERFORMS DIRECT INVERSION IN THE ITERATIVE SUBSPACE         *

!*                                                                      *

!*     THIS INVOLVES SOLVING FOR C IN Geom(NEW) = Geom' - HG'      *

!*                                                                      *

!*  WHERE Geom' = SUM(C(I)Geom(I), THE C COEFFICIENTES COMING FROM      *

!*                                                                      *

!*                   | B   1 | . | C | = | 0 |                          *

!*                   | 1   0 |   |-L |   | 1 |                          *

!*                                                                      *

!*  WHERE B(I,J) =GRAD(I)H(T)HGRAD(J)  GRAD(I) = GRADIENT ON CYCLE I    *

!*                              Hess    = INVERSE HESSIAN        *

!*                                                                      *

!*                          REFERENCE                                   *

!*                                                                      *

!*  P. CSASZAR, P. PULAY, J. MOL. STRUCT. (THEOCHEM), 114, 31 (1984)    *

!*                                                                      *

!************************************************************************

!************************************************************************

!*                                                                      *

!*     GEOMETRY OPTIMIZATION USING THE METHOD OF DIRECT INVERSION IN    *

!*     THE ITERATIVE SUBSPACE (GDIIS), COMBINED WITH THE BFGS OPTIMIZER *

!*     (A VARIABLE METRIC METHOD)                                       *

!*                                                                      *

!*     WRITTEN BY PETER L. CUMMINS, UNIVERSITY OF SYDNEY, AUSTRALIA     *

!*                                                                      *

!*                              REFERENCE                               *

!*                                                                      *

!*      "COMPUTATIONAL STRATEGIES FOR THE OPTIMIZATION OF EQUILIBRIUM   *

!*     GEOMETRIES AND TRANSITION-STATE STRUCTURES AT THE SEMIEMPIRICAL  *

!*     LEVEL", PETER L. CUMMINS, JILL E. GREADY, J. COMP. CHEM., 10,    *

!*     939-950 (1989).                                                  *

!*                                                                      *

!*     MODIFIED BY JJPS TO CONFORM TO EXISTING MOPAC CONVENTIONS        *

!*                                                                      *

!************************************************************************

      ! MRESET = maximum number of iterations.

      INTEGER(KINT), PARAMETER :: MRESET=15, M2=(MRESET+1)*(MRESET+1) !M2 = 256

      REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet(:,:),GradSet(:,:),ERR(:,:) ! MRESET*NCoord

      REAL(KREAL), SAVE :: ESET(MRESET)

    REAL(KREAL), ALLOCATABLE, SAVE :: DXTMP(:,:),GSAVE(:,:) !NCoord, why DXTMP has SAVE attribute??

      REAL(KREAL) :: B(M2), BS(M2)

      LOGICAL DEBUG, PRINT

      INTEGER(KINT), ALLOCATABLE, SAVE :: MSET(:)

    LOGICAL, ALLOCATABLE, SAVE :: FRST(:)

      INTEGER(KINT) :: NDIIS, MPLUS, INV, ITERA, MM, NFree, I, J, K

      INTEGER(KINT) :: JJ, KJ, JNV, II, IONE, IJ, INK,ITmp, Isch, Idx

      REAL(KREAL) :: XMax, XNorm, S, DET, THRES, Norm

    REAL(KREAL), PARAMETER :: eps=1e-12

    REAL(KREAL), PARAMETER :: crit=1e-8

    REAL(KREAL), ALLOCATABLE :: Tanf(:) ! NCoord

    REAL(KREAL), ALLOCATABLE :: HFree(:) ! NFree*NFree

    REAL(KREAL), ALLOCATABLE :: Htmp(:,:) ! NCoord,NFree

    REAL(KREAL), ALLOCATABLE :: Grad_free(:),Grad_new_free_inter(:),Step_free(:) ! NFree

    REAL(KREAL), ALLOCATABLE :: Geom_free(:),Geom_new_free_inter(:),Geom_new_free(:) ! NFree

    REAL(KREAL), ALLOCATABLE, SAVE :: GeomSet_free(:,:),GradSet_free(:,:)

      DEBUG=.TRUE.

      PRINT=.TRUE.

      IF (PRINT)  WRITE(*,'(/,''      BEGIN GDIIS   '')')

      ! Initialization

      IF (allocation_flag) THEN ! allocation_flag = .TRUE. at the begining and effective for all geometries in path.

      ! FRST(IGeom) will be set to False in Space, so no need to modify it here

         IF (ALLOCATED(GeomSet)) THEN

            IF (PRINT)  WRITE(*,'(/,''    In FRST, GDIIS Dealloc  '')')

            DEALLOCATE(GeomSet,GradSet,ERR,DXTMP,GSave,GeomSet_free,GradSet_free)

            RETURN

         ELSE

        ! these allocated arrays need to be properly deallocated.

            IF (PRINT)  WRITE(*,'(/,''     In FRST,  GDIIS Alloc  '')')

            ALLOCATE(GeomSet(NGeomF,MRESET*NCoord),GradSet(NGeomF,MRESET*NCoord),ERR(NGeomF,MRESET*NCoord))

      ALLOCATE(GeomSet_free(NGeomF,MRESET*NCoord),GradSet_free(NGeomF,MRESET*NCoord))

            ALLOCATE(DXTMP(NGeomF,NCoord),GSAVE(NGeomF,NCoord),MSET(NGeomF),FRST(NGeomF))

      DO I=1,NGeomF

         FRST(I) = .TRUE.

         GeomSet(I,:) = 0.d0

         GradSet(I,:) = 0.d0

         ERR(I,:) = 0.d0

         GeomSet_free(I,:) = 0.d0

         GradSet_free(I,:) = 0.d0

         DXTMP(I,:)=0.d0

         GSAVE(I,:)=0.d0

      END DO

      MSET(:)=0

         END IF

     allocation_flag = .FALSE.

      END IF

      ! Addded from here:

    Call FreeMv(NCoord,Vfree) ! VFree(Ncoord,Ncoord)

      ! we orthogonalize Vfree to the tangent vector of this geom only if Tangent/=0.d0

      Norm=sqrt(dot_product(Tangent,Tangent))

      IF (Norm.GT.eps) THEN

         ALLOCATE(Tanf(NCoord))

         ! We normalize Tangent

         Tangent=Tangent/Norm

         ! We convert Tangent into Vfree only displacements. This is useless for now (2007.Apr.23)

     ! as Vfree=Id matrix but it will be usefull as soon as we introduce constraints.

         DO I=1,NCoord

            Tanf(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Tangent)

         END DO

         Tangent=0.d0

         DO I=1,NCoord

            Tangent=Tangent+Tanf(I)*Vfree(:,I)

         END DO

         ! first we subtract Tangent from vfree

         DO I=1,NCoord

            Norm=dot_product(reshape(vfree(:,I),(/NCoord/)),Tangent)

            Vfree(:,I)=Vfree(:,I)-Norm*Tangent

     END DO

         Idx=0

         ! Schmidt orthogonalization of the Vfree vectors

         DO I=1,NCoord

            ! We subtract the first vectors, we do it twice as the Schmidt procedure is not numerically stable.

            DO Isch=1,2

               DO J=1,Idx

                  Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,J),(/NCoord/)))

                  Vfree(:,I)=Vfree(:,I)-Norm*Vfree(:,J)

               END DO

            END DO

            Norm=dot_product(reshape(Vfree(:,I),(/NCoord/)),reshape(Vfree(:,I),(/NCoord/)))

            IF (Norm.GE.crit) THEN

               Idx=Idx+1

               Vfree(:,Idx)=Vfree(:,I)/sqrt(Norm)

            END IF

         END DO

         Print *, 'Idx=', Idx

         IF (Idx/= NCoord-1) THEN

            WRITE(*,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom

            WRITE(IoOut,*) "Pb in orthogonalizing Vfree to tangent for geom",IGeom

            STOP

         END IF

         DEALLOCATE(Tanf)

         NFree=Idx

      ELSE ! Tangent =0, matches IF (Norm.GT.eps) THEN

         if (debug) WRITE(*,*) "Tangent=0, using full displacement"

         NFree=NCoord

      END IF !IF (Norm.GT.eps) THEN

      if (debug) WRITE(*,*) 'DBG Step_DIIS_All, IGeom, NFree=', IGeom, NFree

      ! We now calculate the new step

      ! we project the hessian onto the free vectors

      ALLOCATE(HFree(NFree*NFree),Htmp(NCoord,NFree),Grad_free(NFree),Grad_new_free_inter(NFree))

    ALLOCATE(Geom_free(NFree),Step_free(NFree),Geom_new_free_inter(NFree),Geom_new_free(NFree))

      DO J=1,NFree

        DO I=1,NCoord

           Htmp(I,J)=0.d0

           DO K=1,NCoord

              Htmp(I,J)=Htmp(I,J)+Hess(((I-1)*NCoord)+K)*Vfree(K,J)

           END DO

        END DO

      END DO

      DO J=1,NFree

        DO I=1,NFree

       HFree(I+((J-1)*NFree))=0.d0

           DO K=1,NCoord

        HFree(I+((J-1)*NFree))=HFree(I+((J-1)*NFree))+Vfree(K,I)*Htmp(K,J)

           END DO

        END DO

      END DO

      DO I=1,NFree

         Grad_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Grad)

     Geom_free(I)=dot_product(reshape(Vfree(:,I),(/NCoord/)),Geom)

      END DO

      !Added Ends here.***********************************************

!C  SPACE SIMPLY LOADS THE CURRENT VALUES OF Geom AND GRAD INTO

!C  THE ARRAYS GeomSet AND GradSet

!C  HEAT is never used, not even in Space_all(...)

      CALL Space_all(NGeomF,IGeom,MRESET,MSET,Geom,Grad,HEAT,NCoord,GeomSet,GradSet,ESET,FRST)

      IF (PRINT)  WRITE(*,'(/,''       GDIIS after Space  '')')

    DO J=1,MSet(IGeom)

         DO K=1,NFree

            GradSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&

               GradSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))

        GeomSet_free(IGeom,((J-1)*NFree)+K)=dot_product(reshape(Vfree(:,K),(/NCoord/)),&

                GeomSet(IGeom,((J-1)*NCoord)+1:((J-1)*NCoord)+NCoord))

         END DO

    END DO

!C

!C     INITIALIZE SOME VARIABLES AND CONSTANTS

!C

      NDIIS = MSET(IGeom)

      MPLUS = MSET(IGeom) + 1

      MM = MPLUS * MPLUS

!C

!C     COMPUTE THE APPROXIMATE ERROR VECTORS

!C

      INV=-NFree

      DO 30 I=1,MSET(IGeom)

         INV = INV + NFree

         DO 30 J=1,NFree

            S = 0.D0

            KJ=(J*(J-1))/2

            DO 10 K=1,J

               KJ = KJ+1

   10       S = S - HFree(KJ) * GradSet_free(IGeom,INV+K)

            DO 20 K=J+1,NFree

               KJ = (K*(K-1))/2+J

   20       S = S - HFree(KJ) * GradSet_free(IGeom,INV+K)

   30 ERR(IGeom,INV+J) = S

!C

!C     CONSTRUCT THE GDIIS MATRIX

!C

      DO 40 I=1,MM

   40 B(I) = 1.D0

      JJ=0

      INV=-NFree

      DO 50 I=1,MSET(IGeom)

         INV=INV+NFree

         JNV=-NFree

         DO 50 J=1,MSET(IGeom)

            JNV=JNV+NFree

            JJ = JJ + 1

            B(JJ)=0.D0

            DO 50 K=1,NFree

      !Print *, 'B(',JJ,')=', B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)

   50 B(JJ) = B(JJ) + ERR(IGeom,INV+K) * ERR(IGeom,JNV+K)

     ! The following shifting is required to correct indices of B_ij elements in the GDIIS matrix.

   ! The correction is needed because the last coloumn of the matrix contains all 1 and one zero.

      DO 60 I=MSET(IGeom)-1,1,-1

         DO 60 J=MSET(IGeom),1,-1

   60 B(I*MSET(IGeom)+J+I) = B(I*MSET(IGeom)+J)

      ! For the last row and last column of GEDIIS matrix:

      DO 70 I=1,MPLUS

         B(MPLUS*I) = 1.D0

   70 B(MPLUS*MSET(IGeom)+I) = 1.D0

      B(MM) = 0.D0

!C

!C     ELIMINATE ERROR VECTORS WITH THE LARGEST NORM

!C

   80 CONTINUE

      DO 90 I=1,MM

   90 BS(I) = B(I)

      IF (NDIIS .EQ. MSET(IGeom)) GO TO 140

      DO 130 II=1,MSET(IGeom)-NDIIS

         XMAX = -1.D10

         ITERA = 0

         DO 110 I=1,MSET(IGeom)

            XNORM = 0.D0

            INV = (I-1) * MPLUS

            DO 100 J=1,MSET(IGeom)

  100       XNORM = XNORM + ABS(B(INV + J))

            IF (XMAX.LT.XNORM .AND. XNORM.NE.1.0D0) THEN

               XMAX = XNORM

               ITERA = I

               IONE = INV + I

            ENDIF

  110    CONTINUE

         DO 120 I=1,MPLUS

            INV = (I-1) * MPLUS

            DO 120 J=1,MPLUS

               JNV = (J-1) * MPLUS

               IF (J.EQ.ITERA) B(INV + J) = 0.D0

               B(JNV + I) = B(INV + J)

         !Print *,'B(JNV + I)=',B(JNV + I)

  120    CONTINUE

         B(IONE) = 1.0D0

  130 CONTINUE

  140 CONTINUE

!C

!C     OUTPUT THE GDIIS MATRIX

!C

      IF (DEBUG) THEN

         WRITE(*,'(/5X,'' GDIIS MATRIX'')')

         ITmp=min(12,MPLUS)

         DO IJ=1,MPLUS

            WRITE(*,'(12(F12.4,1X))') B((IJ-1)*MPLUS+1:(IJ-1)*MPLUS+ITmp)

         END DO

      ENDIF

!C

!C     SCALE DIIS MATRIX BEFORE INVERSION

!C

      DO 160 I=1,MPLUS

         II = MPLUS * (I-1) + I

     !Print *, 'B(',II,')=', B(II)

     !Print *, 'GSave(',IGeom,',',I,')=', 1.D0 / DSQRT(1.D-20+DABS(B(II)))

  160 GSAVE(IGeom,I) = 1.D0 / DSQRT(1.D-20+DABS(B(II)))

      GSAVE(IGeom,MPLUS) = 1.D0

      !Print *, 'GSave(',IGeom,',',MPlus,')=1.D0'

      DO 170 I=1,MPLUS

         DO 170 J=1,MPLUS

            IJ = MPLUS * (I-1) + J

  170 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)

!C

!C     OUTPUT SCALED GDIIS MATRIX

!C

      IF (DEBUG) THEN

         WRITE(*,'(/5X,'' GDIIS MATRIX (SCALED)'')')

         ITmp=min(12,MPLUS)

         DO IJ=1,MPLUS

            WRITE(*,'(12(F12.4,1X))') B((IJ-1)*MPLUS+1:(IJ-1)*MPLUS+ITmp)

         END DO

      ENDIF

!C

!C     INVERT THE GDIIS MATRIX B

!C

      CALL MINV(B,MPLUS,DET) ! matrix inversion.

      DO 190 I=1,MPLUS

         DO 190 J=1,MPLUS

            IJ = MPLUS * (I-1) + J

      !Print *, 'B(',IJ,')=', B(IJ)

      !Print *, 'GSAVE(',IGeom,',',I,')=', GSAVE(IGeom,I)

      !Print *, 'GSAVE(',IGeom,',',J,')=', GSAVE(IGeom,J)

      !Print *, 'B(',IJ,')=', B(IJ) * GSAVE(I) * GSAVE(J)

  190 B(IJ) = B(IJ) * GSAVE(IGeom,I) * GSAVE(IGeom,J)

!C

!C     COMPUTE THE INTERMEDIATE INTERPOLATED PARAMETER AND GRADIENT VECTORS

!C

      !Print *, 'MSET(',IGeom,')=', MSET(IGeom), ' MPLUS=', MPLUS

      DO 200 K=1,NFree

         Geom_new_free_inter(K) = 0.D0

         Grad_new_free_inter(K) = 0.D0

         DO 200 I=1,MSET(IGeom)

            INK = (I-1) * NFree + K

            Geom_new_free_inter(K) = Geom_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GeomSet_free(IGeom,INK)

      !Print *, 'Geom_new_free_inter(',K,')=', Geom_new_free_inter(K)

      !Print *, 'B(MPLUS*MSET(',IGeom,')+',I,')=', B(MPLUS*MSET(IGeom)+I)

      !Print *, 'GeomSet_free(',IGeom,',',INK,')=', GeomSet_free(IGeom,INK)

  200 Grad_new_free_inter(K) = Grad_new_free_inter(K) + B(MPLUS*MSET(IGeom)+I) * GradSet_free(IGeom,INK)

      HP=0.D0

      DO 210 I=1,MSET(IGeom)

  210 HP=HP+B(MPLUS*MSET(IGeom)+I)*ESET(I)

      DO 220 K=1,NFree

  220 DXTMP(IGeom,K) = Geom_free(K) - Geom_new_free_inter(K)

      XNORM = SQRT(DOT_PRODUCT(DXTMP(IGeom,1:NFree),DXTMP(IGeom,1:NFree)))

      IF (PRINT) THEN

         WRITE (6,'(/10X,''DEVIATION IN X '',F10.6, 8X,''DETERMINANT '',G9.3)') XNORM,DET

         WRITE(*,'(10X,''GDIIS COEFFICIENTS'')')

         WRITE(*,'(10X,5F12.5)') (B(MPLUS*MSET(IGeom)+I),I=1,MSET(IGeom))

      ENDIF

!C     THE FOLLOWING TOLERENCES FOR XNORM AND DET ARE SOMEWHAT ARBITRARY!

      THRES = MAX(10.D0**(-NFree), 1.D-25)

      IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN

         IF (PRINT)THEN

            WRITE(*,*) "THE DIIS MATRIX IS ILL CONDITIONED"

            WRITE(*,*) " - PROBABLY, VECTORS ARE LINEARLY DEPENDENT - "

            WRITE(*,*) "THE DIIS STEP WILL BE REPEATED WITH A SMALLER SPACE"

         END IF

         DO 230 K=1,MM

  230    B(K) = BS(K)

         NDIIS = NDIIS - 1

         IF (NDIIS .GT. 0) GO TO 80

         IF (PRINT) WRITE(*,'(10X,''NEWTON-RAPHSON STEP TAKEN'')')

         DO 240 K=1,NFree

            Geom_new_free_inter(K) = Geom_free(K)

  240       Grad_new_free_inter(K) = Grad_free(K)

      ENDIF ! matches IF (XNORM.GT.2.D0 .OR. DABS(DET).LT. THRES) THEN, L378

    ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}

    Geom_new_free=0.d0

    DO I = 1, NFree

       DO J = 1, NFree

          ! If Hinv=.False., then we need to invert Hess

        !Geom_new_free(:) = Geom_new_free(:) + HFree(:,I)*Grad_new_free_inter(I)

       Geom_new_free(J) = Geom_new_free(J) + HFree(I+((J-1)*NFree))*Grad_new_free_inter(I)

     END DO

    END DO

    Geom_new_free(:) = Geom_new_free_inter(:) - Geom_new_free(:)

    Step_free = Geom_new_free - Geom_free

    Step = 0.d0

      DO I=1,NFree

     Step = Step + Step_free(I)*Vfree(:,I)

      END DO

      DEALLOCATE(Hfree,Htmp,Grad_free,Grad_new_free_inter,Step_free,Geom_free)

      DEALLOCATE(Geom_new_free_inter,Geom_new_free)

    IF (PRINT)  WRITE(*,'(/,''       END GDIIS  '',/)')

      END SUBROUTINE Step_diis_all