Chimie Théorique » OpenPath

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

  ! 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,XyzGeomI, renum, OrderInv

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

              ! allocated in Path.f90

  use Io_module

  IMPLICIT NONE

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

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

  CHARACTER(LCHARS) :: LINE

  LOGICAL :: debug

  LOGICAL, SAVE :: first=.true.

  REAL(KREAL), SAVE :: Eold=1.e6

  REAL(KREAL), ALLOCATABLE :: GradTmp(:)

  REAL(KREAL), ALLOCATABLE :: GradCart(:)

  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)

  REAL(KREAL) :: d, a_val, Pi

  INTEGER(KINT) :: iat,jat,kat,i,j,IBeg,Idx

  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11

  CHARACTER(LCHARS) :: ListName, CH32SVAR1

  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand

  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.

  LOGICAL   :: FTmp

  INTEGER(kint) :: NbLists,NStep

  INTEGER(KINT), PARAMETER :: NbExtName=4

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

  INTERFACE

     function valid(string) result (isValid)

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

       logical                  :: isValid

     END function VALID

  END INTERFACE

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

  INTEGER(KINT) :: OptGeom

  Namelist/path/OptGeom

  Pi=dacos(-1.0d0)

  debug=valid('EGRAD')

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

  Print *, 'EGrad.f90, L73, IOpt=', IOpt

  Grad=0.

  E=0.

  ALLOCATE(GradCart(3*Nat))

  ALLOCATE(x(Nat),y(Nat),z(Nat))

  SELECT CASE (COORD)

  CASE ('ZMAT')

     ! In order to avoid problem with numbering and linear angles/diehedral, we convert the

     ! zmat into cartesian coordinates

     ! A remplacer par Int2Cart :-)

     Call Int2Cart(nat,indzmat,Geom,GeomCart)

  CASE ('BAKER')

     WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop"

  CASE ('MIXED')

!  write(*,*) "Coucou 4-mixed"

     Call Mixed2Cart(Nat,indzmat,geom,GeomCart)

  CASE ('CART')

!  write(*,*) "Coucou 4-cart"

     if (debug) WRITE(*,*) "Coord=CART... in egrad"

     GeomCart=reshape(Geom,(/Nat,3/))

     DO Iat=1,Nat

        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat-2:3*Iat)

     END DO

  CASE ('HYBRID')

!  write(*,*) "Coucou 4-hybrid"

        GeomCart=reshape(Geom,(/Nat,3/))

  CASE DEFAULT

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

   STOP

  END SELECT

  SELECT CASE (Prog)

     CASE ('GAUSSIAN')

! we call the Gaussian routine.

! it will return the energy and the gradient in cartesian coordinates.

! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)

        Call egrad_gaussian(E,GeomCart,GradCart)

     CASE ('MOPAC')

! we call the Mopac routine.

! it will return the energy and the gradient in cartesian coordinates.

! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)

        Call egrad_mopac(E,GeomCart,GradCart)

     CASE ('VASP')

        Call egrad_vasp(E,Geomcart,GradCart)

     CASE ('TURBOMOLE')

        Call egrad_turbomole(E,Geomcart,GradCart)

     CASE ('EXT')

        Call egrad_ext(E,Geomcart,GradCart)

     CASE ('TEST')

        Call egrad_test(Nat,E,Geomcart,GradCart)

     CASE ('CHAMFRE')

        Call egrad_chamfre(Nat,E,Geomcart,GradCart)

  END SELECT

  if (debug) THEN

     WRITE(*,*) 'DBG EGRAD, GradCart read'

     DO I=1,Nat

        Iat=I

        if (renum) Iat=orderInv(I)

        WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I-2:3*I)

     END DO

  END IF

!  If (PROG=="VASP") STOP

  ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid

  !  but we have to convert it into Zmat if COORD=Zmat

  SELECT CASE (COORD)

     CASE ("ZMAT")

        ALLOCATE(GradTmp(3*Nat))

        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"

        CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)

        if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int"

        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)

        if (debug) WRITE(*,*) "DBG EGRAD Storing Grad"

     Grad(1)=GradTmp(4)

     Grad(2)=GradTmp(7)

! We might have troubles whith angles that are not in the [0:pi] range because,

! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...

! so that the derivative is not good, and a multiplicative factor should be added...

!

! This in fact should be taken care of in B mat calculation...

!

     IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN

        Grad(3)=-1.0d0*GradTmp(8)

     ELSE

        Grad(3)=GradTmp(8)

     END IF

     Idx=4

     DO I=4,Nat

        Grad(Idx)=GradTmp(3*i-2)

        Grad(Idx+1)=GradTmp(3*i-1)

        Grad(Idx+2)=GradTmp(3*i)

        Idx=Idx+3

     END DO

     DEALLOCATE(GradTmp)

   CASE ('BAKER')

        WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop"

     CASE ('MIXED')

        ALLOCATE(GradTmp(3*Nat))

        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"

        CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)

        if (Debug) THEN

           WRITE(*,*) "DzDc"

           DO I=1,Nat

              DO J=1,3

                 WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:)

              END DO

           END DO

        END IF

        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)

        DO I=1,Nat

           WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)

        END DO

        SELECT CASE (NCART)

           CASE (1)

              Grad(1:3)=GradTmp(1:3)

              Grad(4)=GradTmp(4)

              Grad(5)=GradTmp(7)

              Grad(6)=GradTmp(8)

              Idx=7

              IBeg=4

           CASE(2)

              Grad(1:3)=GradTmp(1:3)

              Grad(4:6)=GradTmp(4:6)

              Grad(7)=GradTmp(7)

              Grad(8)=GradTmp(8)

              Idx=9

              IBeg=4

            CASE DEFAULT

               Idx=1

               IBeg=1

            END SELECT

     DO I=IBeg,Nat

        Grad(Idx)=GradTmp(3*i-2)

        Grad(Idx+1)=GradTmp(3*i-1)

        Grad(Idx+2)=GradTmp(3*i)

        Idx=Idx+3

     END DO

     DEALLOCATE(GradTmp)

  CASE ("CART","HYBRID")

     Grad=GradCart

  CASE DEFAULT

     WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP"

     STOP

  END SELECT

  DEALLOCATE(GradCart)

  DEALLOCATE(x,y,z)

  IF (Debug) WRITE(*,*) 'DBG EGRAD grad',grad(1:NCoord)

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

  !999 CONTINUE

  !if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'

  !STOP

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

end subroutine Egrad