c   FORTRAN test program for the distributed subroutine bkmp2.f
c
c   bkmp2 H3 surface test programme and surface evaluation routines
c   to compile, link and run it:
c      >  f77 -c bkmp2.f
c      >  f77 bkmp2test.f bkmp2.o -o bkmp2test
c      >  bkmp2test
c
      program bkmp2test
c----------------------------------------------------------------------c
c put the bkmp2 through various simple checks.
c see 'subroutine bkmp2' for more details about the surface.
c surfid = 950621.0 (a version number based on date of fitting)
c----------------------------------------------------------------------c
      implicit double precision (a-h,o-z)
      dimension r(3),dvdr(3)
      write(6,*) ' SIMPLE CHECKS OF BKMP2 H3 SURFACE'
      call surfid(h3id)
c-sngl;      checkid=950621.
      checkid=950621.d0
      write(6,9506) h3id,checkid
 9506 format('   surface id number: ',1x,f7.0,' [',1x,f7.0,']')
      write(6,*) ' '
      write(6,*) ' Note: computed values should match the values',
     .           ' in square brackets '
      ideriv = 0
c
c     compute saddle energy and force constants
      call kcalc
c
c     compare analytical and numerical derivatives:
      call compare
c
      eps=1.d-12
c
c minimum H2 potential
      write(6,*) ' '
      write(6,*) ' minimum in H2 potential'
      r0 = 1.40125d0
      do 200 i = 1,9
c-sngl;         r(1)= r0 + dfloat(i-1)*0.00005
         r(1)= r0 + dfloat(i-1)*0.00005d0
         r(2)= 99.d0
         r(3)= 99.d0
         call bkmp2(r,vmin,dvdr,ideriv)
         if(dabs(r(1)-1.40145d0).lt.eps) then
            write(6,202) r(1),vmin
            vminsave=vmin
         else
            write(6,201) r(1),vmin
         endif
 200  continue
 202  format(1x,f9.6,5x,f14.11, ' <----')
 201  format(1x,f9.6,5x,f14.11)
c
c collinear saddle point
      write(6,*) ' '
      write(6,*) ' collinear saddle point, with derivatives'
      r0 = 1.7566d0
      ideriv = 1
      do 300 i = 1,9
c-sngl;         r(1)= r0 + dfloat(i-1)*0.0001
         r(1)= r0 + dfloat(i-1)*0.0001d0
         r(2)= r(1)
         r(3)= r(1)+r(2)
         call bkmp2(r,vsad,dvdr,ideriv)
         if(dabs(r(1)-1.7570d0).lt.eps) then
            write(6,401) r(1),vsad,(dvdr(j),j=1,3)
         else
            write(6,400) r(1),vsad,(dvdr(j),j=1,3)
         endif
 401  format(5(1x,f12.9), ' <----')
 400  format(5(1x,f12.9))
 300  continue
c minimum conical intersection = 120 degree saddle point
      r0 = 1.9707d0
      write(6,*) ' '
      write(6,*) ' minimum energy conical intersection',
     $      ' = 120 deg saddle point'
      do 100 i = 1,9
c-sngl;         r(1)= r0 + dfloat(i-1)*0.0001
         r(1)= r0 + dfloat(i-1)*0.0001d0
         r(2)= r(1)
         r(3)= r(1)
         call bkmp2(r,vcon,dvdr,ideriv)
c-sngl;         hei = (vcon-vminsave)*627.51
         hei = (vcon-vminsave)*627.51d0
         if(dabs(r(1)-1.9711d0).lt.eps) then
            write(6,403) r(1),vcon,hei,(dvdr(j),j=1,3)
         else
            write(6,402) r(1),vcon,hei,(dvdr(j),j=1,3)
         endif
 100  continue
 403  format(6(1x,f12.9), ' <----')
 402  format(6(1x,f12.9))
      stop
      end
c
      subroutine kcalc
c----------------------------------------------------------------------c
C calculate the force constants near the saddle
C point of the H3 potential energy surface.  
C july10/89
c reference values: ks=0.1067 ka=-0.058 kb=0.021 (lsth paper)
C
C to calculate ks use the following points:
C  (a) r = 1.757 1.757 3.514 (saddle-point)
C  (b) r = 1.767 1.767 3.534 (for example)
C to calculate ka use the following two points:
C  (a) r = 1.757 1.757 3.514 (saddle-point)
C  (b) r = 1.747 1.767 3.514 (asymmetric point near saddle)
C to calculate kb, use the following two points:
C  (a) the saddle point
C  (b) r = 1.757 1.757 3.5096084 (bent, theta=0.1 radians)
c----------------------------------------------------------------------c
      implicit double precision (a-h,o-z)
c      parameter(rs=1.757d0, drs=0.01d0)
      parameter(rs=1.757d0, drs=0.005d0)
      dimension rsad(3),rk(3),dvdr(3)
      real*8 ks,ka,kb,ksref,karef,kbref,kscheck,kacheck,kbcheck
c-sngl;      data ksref/0.1067/, karef/-0.058/, kbref/0.021/
c-sngl;      data kscheck/0.106798/, kacheck/-0.057904/, kbcheck/0.021013/
      data ksref/0.1067d0/, karef/-0.058d0/, kbref/0.021d0/
      data kscheck/0.106798d0/, kacheck/-0.057904d0/,
     .     kbcheck/0.021013d0/
      ideriv = 0
      rsad(1) = rs
      rsad(2) = rs
      rsad(3) = rs + rs  
      call bkmp2(rsad,esaddle,dvdr,ideriv)
c-sngl;      check = -0.159176000
      check = -0.159176000d0
      write(6,*) ' '
      write(6,*) ' saddle point energy (hartrees)'
      write(6,6002) esaddle,check
 6002 format(5x,f12.9,'  [',f12.9,']')
c
      write(6,*) ' '
      write(6,*) ' nominal minimum of H2 well at 1.401 (hartrees)'
c-sngl;      rk(1) = 1.401
c-sngl;      rk(2) = 99.0
c-sngl;      rk(3) = 99.0
      rk(1) = 1.401d0
      rk(2) = 99.0d0
      rk(3) = 99.0d0
      call bkmp2(rk,vmin,dvdr,ideriv)
c-sngl;      check = -0.174495730
      check = -0.174495730d0
      write(6,6001) vmin,check
 6001 format(5x,f12.9,'  [',f12.9,']')
c
      write(6,*) ' '
      write(6,*) ' classical barrier height',
     $      ' using this minumum (kcal/mol)'
c-sngl;      barrier = (esaddle-vmin)*627.51
c-sngl;      check = 9.613283501
      barrier = (esaddle-vmin)*627.51d0
      check = 9.613283501d0
      write(6,6003) barrier,check
 6003 format(5x,f12.9,'  [',f12.9,']')
c
cd    write(6,*)'ksymmetric value:'
      rk(1) = rs + drs
      rk(2) = rs + drs
      rk(3) = rk(1) + rk(2)
      call bkmp2(rk,ee,dvdr,ideriv)
      ds = ee-esaddle
c-sngl;      xs = rk(1)+rk(2)-2.0*rs
      xs = rk(1)+rk(2)-2.0d0*rs
      xs = xs*dsqrt(3.d0)/2.d0
      ks = 2.d0*ds/xs**2
cd    write(6,*) 'kasymmetric value:'
      rk(1) = rs + drs
      rk(2) = rs - drs
      rk(3) = rk(1) + rk(2)
      call bkmp2(rk,ee,dvdr,ideriv)
      da = ee-esaddle
      xa = 0.5d0*(rk(1)-rk(2))
      ka = 2.d0*da/xa**2
c
cd    write(6,*) 'bending force constant: '   
      rk(1) = rs
      rk(2) = rs
      rk(3) = 3.5096084d0
      call bkmp2(rk,ee,dvdr,ideriv)
      db = ee-esaddle
c-sngl;      xb = 0.5d0*rs*0.1
      xb = 0.5d0*rs*0.1d0
      kb = 2.d0*db/xb**2
      write(6,*) ' '
      write(6,*) ' force constants'
      write(6,6000)
     .              ks,ka,kb,
     .              kscheck,kacheck,kbcheck,
     .              ksref,karef,kbref
 6000 format(
     .  5x,'ksymmetric    kasymmetric    kbending ',/,   
     .  5x,   f9.6, 5x, f9.6, 5x, f9.6,/,
     .'    [',f9.6, 5x, f9.6, 5x, f9.6,']',/,
     .  5x,   f9.6, 5x, f9.6, 5x, f9.6,'   from liu and siegbahn')
      return
      end
c
      subroutine compare
c----------------------------------------------------------------------c
c compare analytical and numerical derivatives for a specified
c geometry containing some compact distances
c----------------------------------------------------------------------c
      implicit double precision (a-h,o-z)
      dimension r(3),dvdr(3),dvdrnum(3),dummy(3),dvdrref(3)
      parameter(dr = 1.d-07)
      data r/ 0.9d0, 1.1d0, 1.5d0 /
c-sngl;      data dvdrref/ -0.524470263, -0.271446000, -0.108650526 /
      data dvdrref/ -0.524470263d0, -0.271446000d0, -0.108650526d0 /
      ideriv = 1
      call bkmp2(r,v0,dvdr,ideriv)
      ideriv = 0
      r(1) = r(1) + dr
      call bkmp2(r,v1,dummy,ideriv)
      r(1) = r(1) - dr
      r(2) = r(2) + dr
      call bkmp2(r,v2,dummy,ideriv)
      r(2) = r(2) - dr
      r(3) = r(3) + dr
      call bkmp2(r,v3,dummy,ideriv)
      r(3) = r(3) - dr
      dvdrnum(1) = (v1-v0)/dr
      dvdrnum(2) = (v2-v0)/dr
      dvdrnum(3) = (v3-v0)/dr
c
      write(6,*) ' '
      write(6,*) ' numerical check of analytical derivatives'
      write(6,6000) r,dvdrnum,dvdr,dvdrref
 6000 format(5x,3(1x,f12.6),'   r/bohr    ',/,
     .       5x,3(1x,f12.9),'   numerical (dr = 1.d-7)',/,
     .       5x,3(1x,f12.9),'   analytical',/,
     .   4x,'[',3(1x,f12.9),']  reference values (bkmp2 analytical)')
      return
      end 
c
