C***********************************************************************
C
C                      Subroutine  CMDVR
C
C
C   Calculate Hamiltonian matrix in Colbert/Miller DVR basis.
C            (see J. Chem. Phys. 96, 1982  (1992) )
C   Kinetic energy elements are that of Eqs. A6 of Appendix A.
C
C   Input:
C     a     :  Number of basis functions and grid points
C              (Note that 'a' corresponds to 'N-1' of the paper)
C     vpot  :  One-dimensional array of potential values in au
C     length:  Length of grid in au (is the term (b-a) in the
C               Appendix A of the paper)
C     mass  :  Mass in au as used by the kinetic energy term
C
C   Output: 
C     ham: Hamiltonian matrix elements in Colbert-Miller DVR basis
C
C                 (hbar is taken as 1)
C
C***********************************************************************


      subroutine cmdvr(ham,a,vpot,length,mass)

      implicit none

      integer i,j,a,n
      real*8  ham(a,a),vpot(a),length,mass,pi,fac1,fac2

      Pi   = 4d0*atan(1d0)
      fac1 = Pi**2/(4d0*mass*length**2)
      n    = a+1
      fac2 = pi/(2d0*(a+1))
      do i=1,a
         do j=1,i-1
            ham(j,i) = fac1*(-1d0)**(j-i)
     +           *(1d0/sin(fac2*(j-i))**2-1d0/sin(fac2*(j+i))**2)
            ham(i,j) = ham(j,i)
         enddo
         ham(i,i) = fac1*((2d0*n**2+1d0)/3d0-1d0/sin(Pi*i/(1d0*n))**2)
         ham(i,i) = ham(i,i)+vpot(i)
      enddo

      return
      end
