   program check_pes

    use potential

    implicit none

    double precision, parameter :: angstrom = 1.889726124994d0
    double precision, parameter :: bohr     = 0.5291772108d0
    double precision, parameter :: wavenm   = 219474.6313705d0
    double precision, parameter :: pi       = 3.1415926535897932384626433832795d0

    integer, parameter :: natoms=5, npoints=100

    double precision   :: xyz(3,natoms,npoints), energies(npoints)
    double precision   :: scr(3,natoms), e, r1, r2, phi, g(3*natoms), h(3*natoms*3*natoms)
    double precision   :: grad(3*natoms), hess(3*natoms*(3*natoms+1)/2)

    integer :: ipt, i, j, k, l
    character(len=80) :: line

    call pesdata_init()

    open(1,file='checkdata')
    read(1,*) line
    do ipt = 1, npoints
     read(1,*) line
     read(1,*) energies(ipt)
     read(1,*) xyz(1:3,1:natoms,ipt)
    enddo
    read(1,*) line
    read(1,*) energies(1)
    read(1,*) xyz(1:3,1:natoms,1)
    read(1,*) grad
    read(1,*) hess
    close(1)

    open(2,file='scatter')
    write(2,*) '   reference data              pes                  difference'
    do ipt = 1,npoints
      call pes(xyz(1,1,ipt),0,e,g,h)
      write(2,*) energies(ipt)*wavenm,e*wavenm,(e-energies(ipt))*wavenm
    enddo

    write(2,*)

    ipt = 1
    call pes(xyz(1,1,ipt),2,e,g,h)
    write(2,*) '   reference grad            pes grad                difference'
    do i = 1,3*natoms
     write(2,'(3F24.16)') grad(i),g(i),(g(i)-grad(i))/grad(i)
    enddo
    write(2,*) '   reference hess            pes hess                difference'
    do i = 1,3*natoms*(3*natoms+1)/2
     write(2,'(3F24.16)') hess(i),h(i),(h(i)-hess(i))/hess(i)
    enddo

    close(2)

    call pesdata_final()

   end program check_pes
