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

!----------------------------------------------------------------------
!  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 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
 
      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 :: 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=int(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
