Chimie Théorique » Opt'n Path

SUBROUTINE Opt_geom(IGeom,Nat,NCoord,Coord,Geom,E,Flag_Opt_Geom,Step_method)

! This subroutine optimizes a geometry.

! Initially, it was mainly here for  debugging purposes..

!so  It has been designed to be almost independent of the rest of the code.

! It is now an (almost) officiel option.

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

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

  use Path_module, only : maxcyc,sthresh,gthresh,Hinv,Smax,Nom,Atome,OrderInv,indzmat,&

       UMat,ScanCoord,Coordinate,NPrim,&

       BTransInv,BTransInv_local,Xprimitive_t,ScanCoord,&

       Coordinate,Pi,UMat_local,NCart, DynMaxStep,FPrintGeom,GeomFile, AtName,Renum,Order, &

       FormAna,NbVar, PrintGeomFactor,AnaGeom

  use Io_module

  IMPLICIT NONE

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

  CHARACTER(32), INTENT(IN) :: Coord

  CHARACTER(32), INTENT(IN) :: Step_method

  REAL(KREAL), INTENT(INOUT) :: Geom(NCoord)

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

  LOGICAL, INTENT(IN) :: Flag_Opt_Geom

  INTERFACE

     function valid(string) result (isValid)

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

       logical                  :: isValid

     END function VALID

     subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,FOptGeom,GeomOld,GeomCart_old)

       ! This routines calculates the energy E and the gradient Grad of

       ! a molecule with Geometry Geom (may be in internal coordinates),

       ! using for now, either Gaussian or Ext, more general later.

       use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &

            prog,NCart,XyzGeomF,IntCoordF,BTransInv, XyzGeomI, &

            GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &

            , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &

            , Order,OrderInv, XPrimitiveF

       ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)

       ! allocated in Path.f90

       use Io_module

       ! Energy (calculated if F300K=.F., else estimated)

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

       ! NCoord: Number of the degrees of freedom

       ! IGeom: index of the geometry.

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

       ! Geometry at which gradient is calculated (cf Factual also):

       REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)

       ! Gradient calculated at Geom geometry:

       REAL(KREAL), INTENT (OUT) :: Grad(NCoord)

       ! Cartesian geometry corresponding to (Internal Geometry) Geom:

       REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)

!!! Optional, just for geometry optimization with Baker coordinates

       REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)

       REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)

       ! FOptGeom is a flag indicating if we are doing a geom optimization

       ! it can be omitted so that we use a local flag: Flag_Opt_Geom

       LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom

       ! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F

       LOGICAL  :: Flag_Opt_Geom

     END subroutine Egrad

  END INTERFACE

  LOGICAL  :: debug

  LOGICAL :: FiniS,FiniG,Fini

  LOGICAL, SAVE :: FRST=.TRUE.

  ! Variables

  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG

  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6

  REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat-6

  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)

  REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)

  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6

  REAL(KREAL), ALLOCATABLE :: Hess_local(:,:)

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

  REAL(KREAL), ALLOCATABLE :: GeomCart_old(:,:) !(Nat,3)

  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)

  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))

  REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6

  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)

  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis

  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep,NormGrad

  REAL(KREAL) :: Eold, Eini

  ! Values contains the values for the geometry analizes

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

  CHARACTER(LCHARS) :: Line

  debug=valid('OptGeom')

  E=0.

  Eold=0.

  Eini=0.

  MaxStep=SMax

  if (debug) Call Header("Entering Geom Optimization")

  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord))

  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))

  ALLOCATE(GeomRef(NCoord))

  ALLOCATE(Hess_local(NCoord,NCoord))

  ALLOCATE(GeomCart(Nat,3),GeomCart_old(Nat,3))

  ALLOCATE(NewGeom(NCoord))

  ALLOCATE(NewGradTmp(NCoord))

  ALLOCATE(Hess_local_inv(NCoord,NCoord))

  if (NbVar>0) THEN

     ALLOCATE(Values(NbVar))

  END IF

  FormAna(5:8)=" I5 "

  IOpt=0

  ZeroVec=0.d0

  ! Initialize the Hessian

  Hess_local=0.

  SELECT CASE (COORD)

 CASE ('ZMAT')

     ! We use the fact that we know that approximate good hessian values

     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles

     Hess_local(1,1)=0.5d0

     Hess_local(2,2)=0.5d0

     Hess_local(3,3)=0.1d0

     DO J=1,Nat-3

        Hess_local(3*J+1,3*J+1)=0.5d0

        Hess_local(3*J+2,3*J+2)=0.1d0

        Hess_local(3*J+3,3*J+3)=0.02d0

     END DO

     IF (HInv) THEN

        DO I=1,NCoord

           Hess_local(I,I)=1.d0/Hess_local(I,I)

        END DO

     END IF

     IF (Step_method .EQ. "RFO") Then

        Hess_local=0.

        DO I=1,NCoord

           Hess_local(I,I)=1.d0

        END DO

     END IF

  CASE ('BAKER')

     ! UMat(NPrim,3*Nat-6)

     BTransInv_local = BTransInv

     UMat_local = UMat

     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,NCoord))

     Hprim=0.d0

     ScanCoord=>Coordinate

     I=0

     DO WHILE (Associated(ScanCoord%next))

        I=I+1

        SELECT CASE (ScanCoord%Type)

        CASE ('BOND')

           Hprim(I,I) = 0.5d0

        CASE ('ANGLE')

           Hprim(I,I) = 0.2d0

        CASE ('DIHEDRAL')

           Hprim(I,I) = 0.1d0

        END SELECT

        ScanCoord => ScanCoord%next

     END DO

     ! Hprim U:

     HprimU=0.d0

     DO I=1, NCoord

        DO J=1,NPrim

           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)

        END DO

     END DO

     ! Hess = U^T Hprim U:

     Hess_local=0.d0

     DO I=1, NCoord

        DO J=1,NPrim

           ! UMat^T is needed below.

           Hess_local(:,I) = Hess_local(:,I) + UMat(J,:)*HprimU(J,I)

        END DO

     END DO

     !Print *, 'Hprim='

     DO I=1,NPrim

        !   WRITE(*,'(10(1X,F6.3))') Hprim(I,:)

     END DO

     !Print *, 'UMat='

     DO I=1,NPrim

        !  WRITE(*,'(3(1X,F6.3))') UMat(I,1:3)

     END DO

     !Print *, 'HprimU='

     DO I=1,NPrim

        !  WRITE(*,'(3(1X,F6.3))') HprimU(I,1:3)

     END DO

     !Print *, 'Hess_local='

     DO I=1,NCoord

        !  WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)

     END DO

     DEALLOCATE(Hprim,HprimU)

     IF (HInv) THEN

        Call GenInv(NCoord,Hess_local,Hess_local_inv,NCoord)

        Hess_local = Hess_local_inv

     END IF

     !Print *, 'Hess_local after inversion='

!     DO I=1,NCoord

        !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)

!     END DO

     IF (Step_method .EQ. "RFO") Then

        Hess_local=0.

        DO I=1,NCoord

           Hess_local(I,I)=0.5d0

        END DO

     END IF

  CASE ('MIXED','CART','HYBRID')

     DO J=1,NCoord

        Hess_local(J,J)=1.

     END DO

  CASE DEFAULT

     WRITE (*,*) "Coord=",TRIM(COORD)," not recognized, opt_Geom.f90, L164. Stop"

     STOP

  END SELECT

  ! Go to optimization

  GeomOld=0.d0 ! Internal coordinates

  GeomCart_old=0.d0 ! Cartesian coordinates

  IF (FPrintGeom) THEN

     OPEN(IOGeom,File=GeomFile)

  END IF

  Fini=.FALSE.

  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))

     IOpt=IOpt+1

     Write(Line,'(1X,A,I5)') "Iteration ",Iopt

     Call Header(TRIM(Line))

     WRITE(IoOut,*) "Current Geometry"

     Line=""

     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)

     if (COORD/="CART") THEN

        WRITE(IoOut,*) "Current Geometry in Cartesian coordinates"

        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)

     END IF

     WRITE(IoOut,*) "Computing energy and gradient"

     ! Calculate the  energy and gradient

     IF (debug) THEN

        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt

        WRITE(*,*) "Energy and Coordinates:"

        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))

        WRITE(*,*) "Geom Old:"

        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,13))

        WRITE(*,*) "Grad Old:"

        WRITE(*,'(F12.6,12(1X,F8.3))') GradOld(1:min(NCoord,13))

     END IF

     !Call EGrad(E,Geom,GradTmp,NCoord)

     IF (COORD.EQ.'BAKER') THEN

        ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).

        ! GeomCart has INTENT(OUT)

        ! GeomRef is modified in Egrag_baker, so we should not use GeomOld variable

        GeomRef=GeomOld

        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &

             GeomRef,GeomCart_old)

        GeomCart_old=GeomCart

     ELSE

        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)

     END IF

!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!

! We have a pb with the format of GradTmp(Nat,3) whereas it is 3,Nat in Egrad

!

! This is a path for CART coordinates !!!

     IF (COORD=="CART") THEN

        Step=reshape(reshape(GradTmp,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))

        GradTmp=Step

     END IF

!!!!!!!!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!!!!!!!!

     WRITE(IoOut,'(1X," Energy E=",F15.6,A)') E,TRIM(UnitProg)

     WRITE(IoOut,*) "Gradient for current geometry"

     Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)

     If (Iopt==1) EIni=E

     IF (debug) THEN

        WRITE(*,*) 'L198, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt

        WRITE(*,*) "Energy and Coordinates:"

        WRITE(*,'(F12.6,12(1X,F8.3))') E,Geom(1:min(NCoord,12))

        WRITE(*,*) "Geom Old:"

        WRITE(*,'(F12.6,12(1X,F8.3))') GEomOld(1:min(NCoord,12))

        WRITE(*,*) "Grad:"

        WRITE(*,'(F12.6,12(1X,F8.3))') GradTmp(1:min(NCoord,12))

     END IF

     IF (FPrintGeom) THEN

        WRITE(IoGeom,*) Nat

        WRITE(IoGeom,"(' Geom. ',I5,' Energy= ',F15.6)") IOpt, E

        DO I=1,Nat

           If (renum) THEN

              Iat=Order(I)

              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(I)),GeomCart(Iat,:),GradTmp(3*Iat-2:3*Iat)

           ELSE

              Iat=OrderInv(I)

              WRITE(IoGeom,'(1X,A10,3(1X,F15.8),5X,3(1X,F15.8))') Trim(AtName(Iat)),GeomCart(I,:),GradTmp(3*I-2:3*I)

           END IF

        END DO

     END IF

     ! do we have to analyze geometries ?

     If (AnaGeom) THEN

        If (NbVar>0) THEN

           Call AnalyzeGeom(GeomCart,Values)

           WRITE(IoDat,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E

           if (debug) WRITE(*,FormAna) Iopt,Values*PrintGeomFactor,(E-Eini)*ConvE,E

        ELSE

           WRITE(IoDat,FormAna) Iopt,(E-Eini)*ConvE,E

           if (debug) WRITE(*,FormAna) Iopt,(E-Eini)*ConvE,E

        END IF

     END IF

     IF (IOpt.GE.2) THEN

        ! We have enough geometries and gradient to update the hessian (or its

        ! inverse)

        GradTmp2=GradTmp-GradOld

        GeomTmp2=Geom-GeomOld

        if (debug) Write(*,*) "Call Hupdate_all, Geom"

        IF (debug) THEN

           WRITE(*,*) "dx before calling",SHAPE(GeomTmp2)

           WRITE(*,'(12(1X,F8.4))') GeomTmp2(1:NCoord)

           WRITE(*,*) "dgrad before calling",SHAPE(GradTmp2)

           WRITE(*,'(12(1X,F8.4))') GradTmp2(1:NCoord)

        END IF

        Call Hupdate_all(NCoord,GeomTmp2,GradTmp2,Hess_local)

     END IF ! matches IF (IOpt.GE.2) THEN

     GradOld=GradTmp

     GeomOld=Geom

     ! GradTmp becomes Step in Step_RFO_all.

     SELECT CASE (Step_method)

     CASE ('RFO')

        Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec)

        GradTmp=Step

     CASE ('GDIIS')

        Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)

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

        Geom=0.d0

        DO I=1, NCoord

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

           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)

        END DO

        Geom(:) = NewGeom(:) - Geom(:)

        ! GradTmp now becomes "step" below:

        GradTmp = Geom - GeomOld

     CASE ('GDIIST')

        Call Step_GDIIS_Simple_Err(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,NCoord,FRST)

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

        Geom=0.d0

        DO I=1, NCoord

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

           Geom(:) = Geom(:) + Hess_local(:,I)*NewGradTmp(I)

        END DO

        Geom(:) = NewGeom(:) - Geom(:)

        ! GradTmp now becomes "step" below:

        GradTmp = Geom - GeomOld

     CASE ('GEDIIS')

        Call Step_GEDIIS(NewGeom,Geom,GradTmp,E,Hess_local,NCoord,FRST)

        ! GradTmp is actually "step" below:

        GradTmp = NewGeom - GeomOld

        !Call Step_GEDIIS_All(3,1,GradTmp,Geom,GradTmp,E,Hess_local,NCoord,FRST,ZeroVec)

     CASE DEFAULT

        WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, opt_Geom.f90, L225. Stop"

        STOP

     END SELECT

     Fini=.TRUE.

     !   If (IOpt.LT.5) GradTmp=(IOpt*GradTmp+(5-IOpt)*0.01*Grad(IGeom,:))/5.d0

     NormStep=sqrt(DOT_Product(GradTmp,GradTmp))

     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh

     FactStep=min(1.d0,MaxStep/NormStep)

     FiniS=((NormStep*FactStep).LE.SThresh)

     NormGrad=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.

     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormGrad (Gradient)=', NormGrad, 'Gradient Threshold=', GThresh

     FiniG=(NormGrad.LE.GThresh)

     Fini=FiniS.AND.FiniG

     if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep

     GradTmp=GradTmp*FactStep ! GradTmp is step  here, from Step_RFO_all.

     WRITE(IoOut,*) " Converged ?"

     WRITE(IoOut,'(1X,A18,1X,A15,1X,A13,5X,A3)') " Criterion","Current Value","Threshold","Ok?"

     IF (FiniS) THEN

        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"YES"

     ELSE

        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"NO"

     END IF

     IF (FiniG) THEN

        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"YES"

     ELSE

        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"NO"

     END IF

     if (DynMaxStep.AND.(IOpt>1)) THEN

        If (E<EOld) Then

           MaxStep=min(1.1*MaxStep,2.*SMax)

        ELSE

           MaxStep=max(0.8*MaxStep,SMax/2.)

        END IF

        if (debug) WRITE(*,*) "Dynamically updated  Maximum step in Opt_Geom:",MaxStep,MaxStep/SMax

     END IF

     Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)

! We add the step to the old geom

     Geom=GeomOld+GradTmp

     EOld=E

     IF (debug) THEN

        WRITE(*,*) "L235, DBG OptGeom, IGeom=",IGeom," IOpt=", IOpt

        SELECT CASE (COORD)

        CASE ('ZMAT','BAKER')

           !WRITE(*,'(12(1X,F8.3))') Geom(1:min(12,NCoord))

        CASE('CART','HYBRID')

           DO Iat=1,Nat

              ! PFL 30 Apr 2008 ->

              ! Deleted ref to orderinv that is created only for coord in zmat, baker, hybrid.

              !                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), &

                   Geom(3*Iat-2:3*Iat)

              ! <- PFL 30 Apr 2008

           END DO

        CASE('MIXED')

           DO Iat=1,NCart

              WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

                   Geom(3*Iat-2:3*Iat)

           END DO

           SELECT CASE (NCart)

           CASE(1)

              if (NAt.GE.2) THEN

                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &

                      IndZmat(2,2),Geom(4)

                 IBEG=5

              END IF

              IF (NAT.GE.3) THEN

                 WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &

                      IndZmat(3,2),Geom(5),IndZmat(3,3),Geom(6)

                 IBeg=7

              END IF

           CASE(2)

              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))')  Nom(Atome(OrderInv(indzmat(3,1)))), &

                   IndZmat(3,2),Geom(7),IndZmat(3,3),Geom(8)

              IBeg=9

           CASE DEFAULT

              IBeg=3*NCart+1

           END SELECT

           DO Iat=max(4,NCart),Nat

              WRITE(*,'(1X,A5,3(1X,I5,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &

                   IndZmat(iat,1),Geom(IBeg),IndZmat(3,2),Geom(IBeg+1), &

                   IndZmat(3,3),Geom(IBeg+2)*180./pi

              IBeg=IBeg+3

           END DO

        CASE DEFAULT

           WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in OptGeom. L294."

           STOP

        END SELECT

     END IF ! matches IF (debug) THEN

  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))

  IF (Fini) THEN

     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles"

  ELSE

     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles"

  END IF

  Call Header(Line)

  WRITE(IoOut,*) "Last Energy E=",E

!  WRITE(*,*) "Initial Geometry:"

!  DO I=1, Nat

!     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)

!  END DO

  WRITE(IoOut,*) "Final Geometry:"

     Line=""

     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)

     if (COORD/="CART") THEN

        WRITE(IoOut,*) "Last Geometry in Cartesian coordinates"

        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)

     END IF

!  DO I=1, Nat

!     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)

     !IF (I .GT. 1) Then

     !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&

     !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&

     !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))

     !END IF

!  END DO

  IF (COORD .EQ. "BAKER") Then

     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).

     ScanCoord=>Coordinate

     DO WHILE (Associated(ScanCoord%next))

        I=I+1

        SELECT CASE (ScanCoord%Type)

        CASE ('BOND')

           WRITE(*,'(1X,I3,":",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,ScanCoord%At2,Xprimitive_t(I)

        CASE ('ANGLE')

           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1, &

                ScanCoord%At2, ScanCoord%At3,Xprimitive_t(I)*180./Pi

        CASE ('DIHEDRAL')

           WRITE(*,'(1X,I3,":",I5," - ",I5," - ",I5," - ",I5," = ",F10.4)') I,ScanCoord%At1,&

                ScanCoord%At2, ScanCoord%At3,ScanCoord%At4,Xprimitive_t(I)*180./Pi

        END SELECT

        ScanCoord => ScanCoord%next

     END DO ! matches DO WHILE (Associated(ScanCoord%next))

  END IF

  DEALLOCATE(GeomCart)

  DEALLOCATE(GradTmp2,GeomTmp2,GradTmp,Step)

  DEALLOCATE(GradOld,GeomOld)

  DEALLOCATE(Hess_local)

  DEALLOCATE(GeomCart_old)

  DEALLOCATE(NewGeom,NewGradTmp)

  DEALLOCATE(Hess_local_inv)

  if (debug) Call Header("Geom Optimization Over")

END SUBROUTINE Opt_geom