Cleaned up the CHARMM gradient code a bit.

gradient.f

      subroutine usualgrd(nn,nbeweg,ndim,izp,lmax,ind,period,doscf,
     *  x,ev,a,b,occ,shift,shift3,shift3A,shift3B,shiftE,shiftE2,
     *  dacc,boxsiz,xinvbox,grad,lldep,lcolspin,spinfac,spinshift)
      use sccdftbsrc, only: mdim, ldim, lscc3rd, nndim

      implicit none

      real*8 deltax
      parameter ( deltax = 0.00001 )
      logical period,doscf,lcolspin,lldep ! MG_QC_UW1207,MG_UW1210
      integer nn,nbeweg,ndim,izp(*),lmax(*),ind(*)
      real*8 x(3,*),ev(*),a(mdim,mdim),b(mdim,mdim),occ(*),grad(3,*)
      real*8 dacc,boxsiz(3,3),xinvbox(3,3),shift(nndim,3) !MG_UW1210
      real*8 shift3(nndim,3),shift3A(nndim,3),shift3b(nndim),dgr3 ! MG+Guanhua_QC_UW1205,MG_UW1210
      real*8 shiftE(*),shiftE2(*) ! MG+Guanhua_QC_UW1206
      real*8 spinfac,spinshift(nndim,3) ! MG_QC_UW1207
      integer m,n,i,j,k,lj,lk,mj,mk
      integer mu,nu,izpj,izpk,indj,indk
      real*8 au(ldim,ldim),bu(ldim,ldim),ocmcc,dgrh,rcdx,xhelp,dgrs,dgr
      real*8 auh(ldim,ldim),buh(ldim,ldim),dgrspin ! MG_QC_UW1207

      real*8 cpe_dqdr(3,ndim)
      real*8 cpe_dsdr

      logical qlink
      common/qlinkl/qlink

      rcdx=1.0/deltax

      ! Reset grad array
      do i=1,nn
        do k=1,3
          grad(k,i) = 0.0d0
        end do
      end do

      dgr=0.0d0
      dgr3=0.0d0

      do n = 1,ndim
        do m = 1,ndim
          b(m,n) = 0.0d0
        end do
      end do

      ! Clear CPE dq/dr tangent matrix
      do n = 1,ndim
        do m = 1,3
          cpe_dqdr(m,n) = 0.0d0
        end do
      end do
!
! for GHO, we have to loop over all orbitals, becuase
! auxiliary orbitals are also occupied ... PJ 7/2004
!
      do i = 1,ndim


        if (.not. qlink) then
          if (occ(i) .lt. dacc) exit
        endif
!
        do m = 1,ndim
          do n = 1,m-1
            ocmcc = occ(i)*a(m,i)*a(n,i)
            b(n,m) = b(n,m) + ocmcc*ev(i)
            b(m,n) = b(m,n) + ocmcc
          end do
        end do

      end do

      do n = 1,ndim
        do m = 1,ndim
          if(abs(b(m,n)) .lt. dacc)then
            b(m,n) = 0.0
          endif
        end do
      end do

      ! Hint: nbeweg = number of atoms
      ! Also: nbeweg = nn
      do j = 1,nbeweg
        indj=ind(j)
        izpj=izp(j)

        ! Hint: nn = number of atoms
        do k = 1,nn
          if(k.ne.j)then
            indk=ind(k)
            izpk=izp(k)

            ! Hint: i = x or y or z
            do i = 1,3
              xhelp = x(i,j)
              x(i,j) = xhelp + deltax
              call slkmatrices(k,j,x,au,bu)
              x(i,j) = xhelp - deltax
              call slkmatrices(k,j,x,auh,buh)
              x(i,j) = xhelp

              do lj = 1,lmax(izpj)
                do mj=1,2*lj-1

                  n = (lj-1)**2 + mj
                  nu = n + indj

                  do lk = 1,lmax(izpk)
                    do mk = 1,2*lk-1

                      m = (lk-1)**2 + mk
                      mu = m + indk
                      dgrh = (au(m,n)-auh(m,n))*rcdx
                      dgrs = -(bu(m,n)-buh(m,n))*rcdx

                      ! DFTB2 term
                      if(doscf)then
                        dgr = -0.5d0*dgrs*(shift(k,1)+shift(j,1)
     &                    - shiftE(k)-shiftE(j)-shiftE2(k)-shiftE2(j))
                      endif

                      ! DFTB3 term
                      if (lscc3rd) then
                        dgr3=-0.5d0*dgrs/3.0d0*(2.0d0*shift3(k,1)
     &                   +shift3A(k,1)+2.0d0*shift3(j,1)+shift3A(j,1))
                      endif

                      cpe_dqdr(i,j) = cpe_dqdr(i,j) + dgrs

                      if(mu.gt.nu)then
                        dgrh = dgrh * b(mu,nu)
                        dgrs = dgrs * b(nu,mu)

                        if(doscf)then
                          dgr = dgr * b(mu,nu)
                        endif

                        if (lscc3rd) then
                          dgr3= dgr3* b(mu,nu)
                        endif

                      else ! not if(mu.gt.nu)
                        dgrh = dgrh * b(nu,mu)
                        dgrs = dgrs * b(mu,nu)
    
                        if(doscf)then
                          dgr = dgr * b(nu,mu)
                        endif

                        if (lscc3rd) then
                          dgr3=dgr3*b(nu,mu)
                        endif

                      endif !if(mu.gt.nu)

                      grad(i,j) = grad(i,j) + dgrh + dgrs

                      if(doscf)then
                        grad(i,j) = grad(i,j) + dgr
                      endif

                      if (lscc3rd) then
                        grad(i,j) = grad(i,j) + dgr3
                      endif

                    end do !nu/atomk
                  end do !lmax/atomk
                end do !mu/atomj
              end do ! lmax/atomj
            end do ! loop i = x/y/z
          endif ! if(k.ne.j)
        end do ! loop atom k
      end do ! loop atom j

      write (*,*) "CPEDQDR"
      write (*,*) cpe_dqdr
      write (*,*) "CEDDQDR"

      end