Chimie Théorique » OpenPath

subroutine Egrad_baker(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom,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

  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)

  REAL(KREAL), INTENT (INOUT) :: GeomOld(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)

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

  LOGICAL,INTENT (IN) :: Flag_Opt_Geom

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

  CHARACTER(LCHARS) :: LINE

  LOGICAL :: debug

  LOGICAL, SAVE :: first=.true.

  LOGICAL, ALLOCATABLE :: Done(:)

  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, n3at,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,LastIt,Nat4,NStep, Firstit

  INTEGER(KINT), PARAMETER :: NbExtName=4

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

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

!

!!!!!!!! PFL Debug

!

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

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

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

  REAL(KREAL), ALLOCATABLE :: XyzTmp(:,:) ! xyzatm(3,nat)

  REAL(KREAL) :: tol

  LOGICAL :: FAIL , FALOBBT,FALOBPrimt,FAloBBTInv

  LOGICAL :: DebugPFL

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

  LOGICAL, EXTERNAL :: valid

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

  INTEGER(KINT) :: OptGeom

  Namelist/path/OptGeom

  Pi=dacos(-1.0d0)

  n3at=3*nat

  debug=valid('EGRAD_baker')

  debugPFL=valid('bakerPFL')

  if (debug) WRITE(*,*) '================ Entering Egrad_baker ================='

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

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

  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; L82, Egrad.f90'

           WRITE(*,*) 'GeomOld='

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

           WRITE(*,*) 'Geom='

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

           WRITE(*,*) 'GeomCart_old='

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

           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; L108, Egrad.f90'

           WRITE(*,*) 'GeomOld='

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

           WRITE(*,*) 'Geom='

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

           WRITE(*,*) 'DBG Egrad_baker L129 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 (*,*) "Coord=",TRIM(COORD),": use Egrad.f90, instead of Egrad_baker.f90. Stop"

  CASE ('CART')

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

  CASE ('HYBRID')

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

  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 Iat=1,Nat

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

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

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

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

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

  CASE ("CART","HYBRID")

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

  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_baker grad',grad(1:NCoord)

  if (debug) WRITE(*,*) '================  Egrad_baker Over ===================='

  !999 CONTINUE

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

  !STOP

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

end subroutine Egrad_baker