Chimie Théorique » Opt'n Path

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, renum, &

       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

  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)

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

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

  LOGICAL :: debug

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

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

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

  REAL(KREAL) ::  Pi

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

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

  INTEGER(KINT), PARAMETER :: NbExtName=4

  REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

!!!!!!!! PFL Debug

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

  LOGICAL ::  FALOBBT,FALOBPrimt,FAloBBTInv

  LOGICAL :: DebugPFL

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

  INTERFACE

     function valid(string) result (isValid)

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

       logical                  :: isValid

     END function VALID

  END INTERFACE

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

  Pi=dacos(-1.0d0)

  if (present(FOptGeom)) THEN

     Flag_Opt_Geom=FOptGeom

  ELSE

     Flag_Opt_Geom=.FALSE.

  END IF

  debug=valid('EGRAD')

  debugPFL=valid('bakerPFL')

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

  if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom

  Grad=0.

  E=0.

  ALLOCATE(GradCart(3*Nat))

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

  ALLOCATE(XPrimRef(NPrim),XPrim(NPrim))

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

     XPrimRef=XPrimitiveF(IGeom,:)

     IF (Flag_Opt_Geom) THEN

        IF (IOpt .EQ. 1) THEN

           GeomCart(:,1) = XyzGeomI(IGeom,1,:)

           GeomCart(:,2) = XyzGeomI(IGeom,2,:)

           GeomCart(:,3) = XyzGeomI(IGeom,3,:)

        ELSE

           if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90'

           if (.NOT.present(GeomOld)) THEN

              WRITE(*,*) "ERROR: GeomOld not passed as an argument"

              STOP

           END IF

           WRITE(*,*) 'GeomOld='

           WRITE(*,'(12(1X,F6.3))') GeomOld

           WRITE(*,*) 'Geom='

           WRITE(*,'(12(1X,F6.3))') Geom

           if (.NOT.present(GeomCart_Old)) THEN

              WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument"

              STOP

           END IF

           WRITE(*,*) 'GeomCart_old='

           WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3)

           Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), &

                GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef)

           GeomCart(:,1) = x(:)

           GeomCart(:,2) = y(:)

           GeomCart(:,3) = z(:)

        END IF

     ELSE

        IF (IOpt .EQ. 0) THEN

           GeomCart(:,1) = XyzGeomF(IGeom,1,:)

           GeomCart(:,2) = XyzGeomF(IGeom,2,:)

           GeomCart(:,3) = XyzGeomF(IGeom,3,:)

        ELSE

           ! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat)

           ! Geom has to be converted into x,y,z

           ! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp.

           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()

              BTransInv_local(I,:) = BTransInvF(IGeom,I,:)

              UMat_local(:,I) = UMatF(IGeom,:,I)

           END DO

           if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90'

           WRITE(*,*) 'GeomOld='

           WRITE(*,'(12(1X,F6.3))') GeomOld_all(IGeom,:)

           WRITE(*,*) 'Geom='

           WRITE(*,'(12(1X,F6.3))') Geom

           WRITE(*,*) 'DBG Egrad L165 GeomCart_old='

           DO I=1,Nat

              WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I)

           END DO

           Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), &

          XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef)

           if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90'

           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()

              BTransInvF(IGeom,I,:) = BTransInv_local(I,:)

           END DO

           ! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine.

           GeomCart(:,1) = x(:)

           GeomCart(:,2) = y(:)

           GeomCart(:,3) = z(:)

        END IF ! matches IF (IOpt .EQ. 0) THEN

     END IF

  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

  !WRITE(*,*) Nat

  !WRITE(*,*) 'GeomCart in Egrad_baker'

  !DO I=1,Nat

  !  WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3)

  !END DO

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

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

     ! but we have to convert it into internal coordinates if COORD=Baker

!!!!!!!!!!!!!!!!!!!!

     !

     !   PFL Debug

     !

!!!!!!!!!!!!!!!!!!!!

     if (debugPFL) THEN

        WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc"

     if (.not.ALLOCATED(indzmat)) THEN

        ALLOCATE(indzmat(Nat,5))

     IndZmat=0

     Indzmat(1,1)=1

     IndZmat(2,1)=2

     IndZmat(2,2)=1

     IndZmat(3,1)=3

     IndZmat(3,2)=2

     IndZmat(3,3)=1

     IndZmat(4,1)=4

     IndZmat(4,2)=3

     IndZmat(4,3)=2

     IndZmat(4,4)=1

     END IF

     IF (.NOT.ALLOCATED(DzDc)) THEN

        ALLOCATE(DzDc(3,nat,3,Nat))

     END IF

     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))

     DO I=1,Nat

     GeomTmp(1,I)=GeomCart(OrderInv(I),1)

     GeomTmp(2,I)=GeomCart(OrderInv(i),2)

     GeomTmp(3,I)=GeomCart(OrderInv(i),3)

     END DO

     DzDc=0.d0

     CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome)

     END IF ! if (debugPFL)

!!!!!!!!!!!!!!!!!!!!!!!!!

     ! Debugging PFL

!!!!!!!!!!!!!!!!!!!!!!!!

     WRITE(*,*) "BTransInv_local Trans that is used originally"

     DO J=1,3*Nat

        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)

     END DO

     WRITE(*,*) "Calculating actual values using GeomCart"

     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))

     GeomTmp(1,:)=GeomCart(1:Nat,1)

     GeomTmp(2,:)=GeomCart(1:Nat,2)

     GeomTmp(3,:)=GeomCart(1:Nat,3)

     ! Computing B^prim:

     FAloBBT=.NOT.ALLOCATED(BBT)

     FAloBBTInv=.NOT.ALLOCATED(BBT_inv)

     FAloBPrimT=.NOT.ALLOCATED(BprimT)

     if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim))

     if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord))

     if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord))

     BprimT=0.d0

     ScanCoord=>Coordinate

     I=0

     DO WHILE (Associated(ScanCoord%next))

        I=I+1

        SELECT CASE (ScanCoord%Type)

        CASE ('BOND')

           CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I))

        CASE ('ANGLE')

           CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I))

        CASE ('DIHEDRAL')

           CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I))

        END SELECT

        ScanCoord => ScanCoord%next

     END DO

     if (debug) THEN

        WRITE(*,*) "Bprim "

        DO J=1,3*Nat

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

        END DO

     END IF

     BMat_BakerT = 0.d0

     DO I=1,NCoord

        DO J=1,NPrim

           !BprimT is transpose of B^prim.

           !B = UMat^T B^prim, B^T = (B^prim)^T UMat

           BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I)

        END DO

     END DO

     BBT = 0.d0

     DO I=1,NCoord

        DO J=1,3*Nat   ! BBT(:,I) forms BB^T

           BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)

        END DO

     END DO

     BBT_inv = 0.d0

     Call GenInv(NCoord,BBT,BBT_inv,NCoord)

!     ALLOCATE(V(NCoord,NCoord))

!     tol = 1e-10

!     ! BBT is destroyed by GINVSE

!     Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT)

!     DEALLOCATE(V)

!     IF(Fail) Then

!        WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90"

!        STOP

!     END IF

     !Print *, 'BBT_inv='

     DO J=1,NCoord

        ! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J)

        ! Print *, BBT_inv(:,J)

     END DO

     !Print *, 'End of BBT_inv'

     ! Calculation of (B^T)^-1 = (BB^T)^-1B:

     BTransInv_local = 0.d0

     DO I=1,3*Nat

        DO J=1,NCoord

           BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J)

        END DO

     END DO

     Print *, 'BMat_BakerT='

     DO J=1,3*Nat

        WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J)

     END DO

     WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart"

     DO J=1,3*Nat

        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)

     END DO

     if (FAloBPrimT) DEALLOCATE(Bprimt)

     if (FAloBBT) DEALLOCATE(BBt)

     if (FAloBBTInv) DEALLOCATE(BBt_inv)

     Grad=0.d0

     DO I=1, 3*Nat

        ! Here Grad is gradient in Baker coordinates. GradCart(3*Nat)

        Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90

     END DO

     !WRITE(*,*) 'In Egrad.f90, L263, Gradient:'

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

     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