Chimie Théorique » Opt'n Path

      subroutine CalcHess (nfree,geom, hess,Igeom,Iopt)

! This routine calculate the Hessian (or its inverse) numerically using Egrad.

! It should be general and should work with any combination of Coord and Prog

! however, it has not been tested with Prog different from 'test'

! It comes mostly from the hess.f subroutine of the IRC package.

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

!  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 : Nat,Hinv

  use Io_module

!

      IMPLICIT NONE

!!!!!!!!!!! Parameters

  ! NFree==NCoord: Number of the degrees of freedom

  ! IGeom: index of the geometry.

  ! Iopt: step in the optimization process

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

! Geometry at which Hessian is calculated (or its inverse)

   REAL(KREAL), INTENT(IN) :: Geom(Nfree)

! Hessian in compacted form

      REAL(KREAL), INTENT(OUT) ::  Hess(Nfree,Nfree)

!

! Parameters used to calculate the Hessian

! numerically with increasing accuray

      INTEGER(KINT), DIMENSION(6,0:3), PARAMETER :: IECOEF  &

      = RESHAPE(SOURCE=(/           &

          -1,1,0,0,0,0,             &

          -1,1,0,0,0,0,             &

          1,-8,8,-1,0,0,            & 

          -1,9,-45,45,-9,1          &

       /), SHAPE=(/6,4/))

      INTEGER(KINT), DIMENSION(6,0:3), PARAMETER :: ISCOEF   &

      = RESHAPE(SOURCE=(/    &

         0,1,0,0,0,0,         &

          -1,1,0,0,0,0,       & 

          -2,-1,1,2,0,0,      &

           -3,-2,-1,1,2,3     &

         /),SHAPE=(/6,4/))

      INTEGER(KINT), DIMENSION(0:3), PARAMETER :: NbDiv = (/1,2,12,60/)

      INTEGER(KINT), DIMENSION(0:3), PARAMETER :: NbCoef = (/2,2,4,6/)

      INTEGER(KINT), PARAMETER :: HessAcc=2

      REAL(KREAL), PARAMETER :: StepHess=0.01d0

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

      logical           :: debug

      REAL(KREAL) :: Ener

      INTEGER(KINT) :: I, J, IAT, IConst

       REAL(KREAL), ALLOCATABLE :: Temp2(:, :) 

       REAL(KREAL), ALLOCATABLE :: Value(:), ValueOld(:), GradTMP(:)

       REAL(KREAL), ALLOCATABLE :: RTmp1(:,:), RTmp3(:,:)  ! Nat,3

       REAL(KREAL), ALLOCATABLE :: RTmp2(:)  ! NFree

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

  INTERFACE

     function valid(string) result (isValid)

       CHARACTER(*), intent(in) :: string

       logical                  :: isValid

     END function VALID

  END INTERFACE

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

      hess=0.d0

      debug=valid('CALCHESS').OR.valid('hess')

      if (debug) THEN

         WRITE(*,*) 'Entering CalcHess'

         WRITE(*,*) 'StepHess=',StepHess

         WRITE(*,*) "HessAcc=",HessAcc

      END IF

      allocate(Value(Nfree), ValueOld(Nfree), GradTmp(Nfree))

      ALLOCATE(RTmp1(Nat,3),RTmp3(Nat,3),RTmp2(NFree))

      allocate(temp2(Nfree,Nfree))

      RTmp1=0.d0

      RTmp3=0.d0

      RTmp2=Geom

!     Now, we calculate the second derivatives for the active variables               

      Value=Geom

      ValueOld=Geom

      Do IConst=1,Nfree

            DO I=1,NbCoef(HessAcc)

               Value(IConst)=valueOld(IConst)+StepHess*ISCOEF(I,HessAcc)

!     We 'propagate this new value in the Zmatrix'

               call egrad(ener,value,gradtmp,NFree,IGeom,Iopt,RTmp1,RTmp2,RTmp3,.FALSE.)

               Hess(1:Nfree,IConst)= Hess(1:Nfree,IConst)+   &

                    gradtmp(1:Nfree)*IECOEF(I,HessAcc)

            end do

         END DO

        deallocate(Value, ValueOld, GradTmp)

        HESS=Hess/(NbDiv(HessAcc)*StepHess)

         if (debug) then

            WRITE(*,*) 'Non symmetrized hessian'

            DO iat=1,Nfree

               WRITE(*,'(15(1X,F12.5))') Hess(iat,1:Nfree)

            END DO

         END if

         DO i=1,Nfree

           Temp2(i,i)=Hess(i,i)

           DO j=i+1,Nfree

             Temp2(i,j)=0.5*(Hess(i,j)+Hess(j,i))

             Temp2(j,i)=Temp2(i,j)

           END DO

         END DO

         Hess=Hess-Temp2

         if (debug) then

            WRITE(*,*) 'Symmetrized hessian'

            DO iat=1,Nfree

               WRITE(*,'(15(1X,F12.5))') Temp2(iat,1:Nfree)

            END DO

            WRITE(*,*) 'Difference between non-syma and Symmetrized hessian'

            DO iat=1,Nfree

               WRITE(*,'(15(1X,F12.5))') Hess(iat,1:Nfree)

            END DO

         end if

! We calculate inverse of the Hessian only if HInv=T

         Hess=Temp2

      if (Hinv) then

        call Geninv(Nfree,Hess,Temp2,NFree)

        Hess=Temp2

      end if

      deallocate(RTmp1,RTmp2,RTmp3,Temp2)

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

      end subroutine CalcHess