Chimie Théorique » OpenPath

      SUBROUTINE ConvGrad_Cart2Int(n,BMat,indzmat,GradCart,GradInt)

        use Io_module

        use Path_module, only : atome

        IMPLICIT NONE

!

        INTEGER(KINT), INTENT(IN) ::  n, indzmat(n,5)

        REAL(KREAL), INTENT(IN)  :: BMat(3,n,3,n), GradCart(3,n)

        REAL(KREAL), INTENT(INOUT)  :: GradInt(3,n)

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

!

! This routines converts a gradient in cartesian coordinates into

! the gradient in internal coordinates (either Zmat, Mixed or Baker)

!

! It uses a LAPACK subroutine (dgelsd) which purpose is:

! *  DGELSD computes the minimum-norm solution to a real linear least

! *  squares problem:

! *      minimize 2-norm(| b - A*x |)

! *  using the singular value decomposition (SVD) of A. A is an M-by-N

! *  matrix which may be rank-deficient.

! Variables:

! input  : n, number of atoms

!          BMat: the dZ/dx matrix

!          GradCart: Gradient in cartesian coordinates

!

! output : GradInt: Gradient in internal coordintaes

!

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

      logical    debug, init

      REAL(KREAL), parameter :: zero = 0d0, tol = 1d-8,  eps = 1d-12

      REAL(KREAL) , ALLOCATABLE :: AEigVec(:,:), AEigVal(:),ATA(:)

      INTEGER(KINT) :: n3, iat, k

      INTEGER(KINT) :: i

!

      save  debug, init

      data  init /.false./

!     .. Local Arrays  for LAPACK

      INTEGER(KINT) :: Info, LWork, Rank

      REAL(KREAL), ALLOCATABLE :: S(:), Work(:), GradCTmp(:,:)

!     IWORK dimension should be at least 3*MIN(M,N)*NLVL + 11*MIN(M,N),

!     where NLVL = MAX( 0, INT( LOG_2( MIN(M,N)/(SMLSIZ+1) ) )+1 )

!     and SMLSIZ = 25 

      INTEGER(KINT) :: LIWork, NLVL, SMLSIZ=12 ! We take smaller value just to be sure

      INTEGER(KINT), ALLOCATABLE :: IWork(:) ! IWork

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

!

  INTERFACE

     function valid(string) result (isValid)

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

       logical                  :: isValid

     END function VALID

     SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, &

                        WORK, LWORK, IWORK, INFO )

! *

! *  -- LAPACK driver routine (version 3.2.2) --

! *  -- LAPACK is a software package provided by Univ. of Tennessee,    --

! *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--

! *     June 2010

! *

! *     .. Scalar Arguments ..

      INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK

      DOUBLE PRECISION   RCOND

!*     ..

!*     .. Array Arguments ..

      INTEGER            IWORK( * )

      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), S( * ), WORK( * )

    END SUBROUTINE DGELSD

  END INTERFACE

!

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

      if (.not.init) then

        init = .true.

        debug = valid ('GEOSTEP').OR.valid('ConvGrad')

      endif

!

      Allocate(GradCTmp(3,n))

      if (debug) WRITE(*,*) "================ Entering ConvGrad_Cart2Int ============="

!

      n3 = 3 * n

!

      NLVL= MAX( 0, INT( LOG( n3*1.d0/(SMLSIZ+1))/LOG(2.d0) ) +1 )

      LIWork=Max(1,3*n3*NLVL + 11*n3+1)

      Allocate(IWork(LIWork),WOrk(1))

      LWork=-1

      GradInt=0.d0

        call DGELSD( n3, n3, 1, BMat, n3, GradCart, n3, GradInt, Tol, RANK, &

                        WORK, LWORK, IWORK, INFO )

        LWork=Work(1)

        if (debug) WRITE(*,*) "DEBUG CONVGRAD_CART2INT LWORK,N3=",LWORK,N3

        DEALLOCATE(Work)

        ALLOCATE(Work(LWork))

        call DGELSD( n3, n3, 1, BMat, n3, GradCart, n3, GradInt, Tol, RANK, &

                        WORK, LWORK, IWORK, INFO )

        GradInt=GradCart

      if (debug) then

         WRITE(*,*) "DEBUG CONVGRAD_CART2INT -- GradInt after DGELSD"

         WRITE(*,*) "INFO=",INFO

         DO I=1,N

            WRITE(*,'(1X,I5,3(1X,F13.6))') I, GradInt(1:3,I)

         END DO

      END IF

! In case we used a Zmat, we might have dummy atoms.

! We are extra cautious and put their gradient components to 0

      do iat = 1, n

        do k = 1, 3

          if (atome(indzmat(iat,1)) == 0) then

            GradInt(k,iat) = zero

          endif

        end do

      end do

      where (abs(GradInt)<=eps) GradInt=zero

      if (debug) write (*,*) "DBG ConvGrad_Cart2Int GradInt=",GradInt

      DEALLOCATE(IWork,Work,GradCTmp)

!

!

      if (debug) WRITE(*,*) "================  ConvGrad_Cart2Int Over ============="

      return

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

    end SUBROUTINE ConvGrad_Cart2Int