!#######################################################################
! In this routine a list of surfaces is available
! 1. it is checked that modes in calculation are in correct order.
! 2. A number is assigned identifying the surface
! 3. Any default cut-offs are also assigned 
! 4. Any options are assigned
!
! ifile=201
!
! Routines on this modul
! defsrf   : defines the surface by reading and interpreting the label.
! vpoint   : evaluates the potential at one point by calling the 
!            appropriate surface routine. 
! surfinfo : prints information on the surfaces and their coordinates
!            to log file. 
!
!  For user specific surfaces use modul usersrf.F
!
!#######################################################################

      module funcsrf

use decimal, only: long, dop
use logdat
use global
use channels
use yzclmod
use NH3pot, potc_nh3=>pot, potc_nh3_init=>prepot
use h2copot, pot_h2co=>Getpot
use maxv, only: maxhop
use keyunitsmod
use dvrdatmod
use griddatmod
use operdef
use runopermod
use datenmod

use coortransmod, only:coortransj,coortrans,coortransb
use splinepot2dmod
use splinepot3dmod

use funcshep

      implicit none
      private
      public :: defsrf,vpoint,surfinfo,read1dinfo
      logical,save::NH3init=.false.
      integer,dimension(7),save::NH3map
      contains

!#######################################################################

      subroutine defsrf(hopilab,label,ifile,pesopts,zpesopts,hopipar, &
                 hoprpar,maxhpar,pesmin,pesmax,lcutpes,ndof,ierr)


      implicit none

      integer(long)                                     :: nopt,ilbl,jlbl,np, &
                                                           icoord,imod,jj,nk, &
                                                           nofpars,length,ll,i1,i2
      integer(long), intent(in)                         :: maxhpar,ndof
      integer(long), intent(out)                        :: hopilab, ifile,ierr
      integer(long), dimension(*), intent(out)          :: hopipar
      integer(long), dimension(2*maxhpar,*), intent(in) :: zpesopts

      real(dop)                                         :: tfac,tfac1,tfac2, &
                                                           width
      real(dop), intent(inout)                          :: pesmin,pesmax
      real(dop), dimension(*), intent(out)              :: hoprpar
      character(len=*), intent(in)                      :: label
      character(len=*), dimension(*), intent(in)        :: pesopts
      character(len=c5)                                 :: filename
      logical(kind=4), intent(out)                      :: lcutpes
      logical(kind=4)                                   :: ldoferr,ascii
      logical(kind=4)                                   :: NH3isCart,NH3isNM
      integer::rrd

      integer(long), save   :: ich=170

      hopilab=0
      ifile=0
      ierr=0
      ldoferr=.false.
      nofpars=0

!-----------------------------------------------------------------------
! Determine label length
!-----------------------------------------------------------------------
      ilbl=index(label,'{')-1
      jlbl=index(label,' ')-1
      if(ilbl.le.0) ilbl=32
      if(jlbl.le.0) jlbl=32
      ll = min(ilbl,jlbl)

!-----------------------------------------------------------------------
! NOCl S1 surface from Uwe Manthe
!-----------------------------------------------------------------------
      if (label .eq. 'nocl1srf') then
         hopilab=1
         lcutpes=.true.
         if (pesmin .lt. -99999.0_dop) pesmin=-1.0_dop/ev
         if (pesmax .gt.  99999.0_dop) pesmax=5.0_dop/ev

!-----------------------------------------------------------------------
! LSTH surface
! imod: 1 collinear
!     : 2 3D
!     : 3 2D 
! icoord : 1 binding coordinates
!        : 2 Jacobian coordinates
!        : 3 kind of hyperspherical coordinates (K(T) calculations).
!        : for 3 see equation (14) in F. Matzkies and Ywe Manthe
!        : J. Chem. Phys. 106, 2646 (1997) 
!-----------------------------------------------------------------------
!        RATE CONSTANT CALCULATIONS OF H + H2
!-----------------------------------------------------------------------
!        : 3 kind of hyperspherical coordinates.
!        : for 3 see equation (14) in F. Matzkies and Uwe Manthe
!        : J. Chem. Phys. 106, 2646 (1997) 
!        : 4 potential in two dimensions for fixed reaction coordinate.
!        : 5 Cartesian Jacobian coordinates (3 of)
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'lsth') then
!
! set default options
!
         if (ndof .eq. 2 ) then
            imod=1
         else
            imod=2
         endif
         icoord=2
         tfac=0.5_dop
!
! read options
!
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   10       continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'tfac') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) tfac
               if (ierr .ne. 0) then
                  write(6,'(i3,a)') np, pesopts(np)(1:72)
                  write(6,'(3i5,a)') ilbl,jlbl,pesopts(np)(ilbl:jlbl)
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read argument to LSTH tfac keyword'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'hyper3D') then
               icoord=3
            else if (pesopts(np)(ilbl:jlbl) .eq. 'hyper2D') then
               icoord=4
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'cartjacobian') then
               icoord=5
            else if (pesopts(np)(ilbl:jlbl) .eq. 'colinear') then
               imod=1
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 10
         endif
!
! now assign surface
!
         if (imod .eq. 1 .and. icoord .eq. 2) then
            hopilab=2   !  LSTH PES : collinear model in Jacobian Coord.  
         elseif (imod .eq. 1 .and. icoord .eq. 1) then
            hopilab=3   !  LSTH PES : collinear model in Binding Coord.
         elseif (imod .eq. 2 .and. icoord .eq. 2) then
            hopilab=4   !  LSTH PES : 3D model in Jacobian Coordinates.
         elseif (imod .eq. 2 .and. icoord .eq. 1) then
            hopilab=5   !  LSTH PES : 3D model in Binding Coordinates.
         elseif (imod .eq. 2 .and. icoord .eq. 3) then
            hopilab=6  !  LSTH PES : 3D model in Hyperspherical Coord.
         elseif (imod .eq. 2 .and. icoord .eq. 4) then
            hopilab=7  !  LSTH PES : 2D model in Hyperspherical Coord.
         elseif (imod .eq. 2 .and. icoord .eq. 5) then
            hopilab=32 !  LSTH PES : 3D model in 3D Cartesian representation of Jacobis (x1,z1,z3)
         endif

         lcutpes=.true.
         if (pesmax .gt.  99999.0_dop) pesmax=0.0_dop

         hoprpar(1)=tfac

!-----------------------------------------------------------------------
! Default surface for elsewhere created vpot files
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'none') then
         hopilab=8

!-----------------------------------------------------------------------
! YZCL1,YZCL2 Surfaces for H3O from Yang,Zhang,Collins and Lee
!             (JCP, 115 174 (2001)).
!
! isurf  : 1 - YZCL1 Surface
!          2 - YZCL2 Surface
! icoord : 1 - Binding coordinates
!          2 - Jacobian coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll-1) .eq. 'yzcl') then

           read(label(5:5),*) isurf
           if (isurf .ne. 1 .and. isurf .ne. 2) then
              routine='defsrf (opfuncs/funcsrf.F)'
              message='Surface '//label(1:5)//' not found'
              ierr=1
              go to 999
           endif
! - Set default options
           icoord = 2
           tfac1 = 0.940735    !#  15.9994/(1.00794+15.9994)
           tfac2 = 0.5_dop

! - Read options
           ilbl=index(label,'{')+1
           jlbl=index(label,'}')-1
           if (ilbl .gt. 1) then
              read(label(ilbl:jlbl),*) np
              nopt=1
              ilbl=1
              jlbl=zpesopts(1,np)
   40         continue
              if (pesopts(np)(ilbl:jlbl) .eq. 'tfac1') then
                 ilbl=zpesopts(nopt+1,np)+1
                 jlbl=zpesopts(nopt+2,np)
                 nopt=nopt+2
                 read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) tfac1
                 if (ierr .ne. 0) then
                    routine='defsrf (opfuncs/funcsrf.F)'
                    message='Cannot read argument to YZCL tfac1 keyword'
                    go to 999
                 endif
              else if (pesopts(np)(ilbl:jlbl) .eq. 'tfac2') then
                 ilbl=zpesopts(nopt+1,np)+1
                 jlbl=zpesopts(nopt+2,np)
                 nopt=nopt+2
                 read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) tfac2
                 if (ierr .ne. 0) then
                    routine='defsrf (opfuncs/funcsrf.F)'
                    message='Cannot read argument to YZCL tfac2 keyword'
                    go to 999
                 endif
              else if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
                 icoord=1
              else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
                 icoord=2
              endif
              nopt=nopt+1
              ilbl=jlbl+1
              jlbl=zpesopts(nopt,np)
              if (jlbl .gt. 0) go to 40
           endif

! - Assign surface

           if (icoord.eq.1) then
              hopilab=9 ! YZCL PES in Binding Coordinates
           else if (icoord.eq.2) then
              hopilab=10  !  YZCL PES in Jacobian Coordinates
           endif
           if (tfac1 .gt. 1.0_dop .or. tfac1 .le. 0.0_dop) then
              routine ='defsrf (opfuncs/funcsrf.F)'
              message='tfac1 in YZCL must be in interval (0,1) .'
              ierr=1
              go to 999
           endif
           if (tfac2 .gt. 1.0_dop .or. tfac2 .le. 0.0_dop) then
              routine='defsrf (opfuncs/funcsrf.F)'
              message='tfac2 in YZCL must be in interval (0,1) .'
              ierr=1
              go to 999
           endif

! -- Read data for the surfaces --
           call yzclinput

           hoprpar(1)=tfac1
           hoprpar(2)=tfac2
           hopipar(1)=isurf

!-----------------------------------------------------------------------

      else if (label(1:ll) .eq. 'cpp') then

! - Set default options
         nk=2
         jj=1

! - Read options
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   71       continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'jtot') then
               ilbl=zpesopts(2,np)+1
               jlbl=zpesopts(3,np)
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) jj
               if (ierr .ne. 0) then
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read jtot argument for cpp keyword'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'dim') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) nk
               if (ierr .ne. 0) then
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read dim argument for cpp keyword'
                  go to 999
               endif
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 71
         endif
         if (jj.lt.1) then
            routine='defsrf (opfuncs/funcsrf.F)'
            message='Jtot must be positive.'
            call errormsg
         endif

!  -- Assign surface
         hopilab=11
         hopipar(1)=jj
         hopipar(2)=nk

!-----------------------------------------------------------------------

      else if (label(1:ll) .eq. 'cmm') then

! - Set default options
         nk=2
         jj=1

! - Read options
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   72       continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'jtot') then
               ilbl=zpesopts(2,np)+1
               jlbl=zpesopts(3,np)
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) jj
               if (ierr .ne. 0) then
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read jtot argument for cmm keyword'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'dim') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) nk
               if (ierr .ne. 0) then
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read dim argument for cmm keyword'
                  go to 999
               endif
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 72
         endif
         if (jj.lt.1) then
            routine='defsrf (opfuncs/funcsrf.F)'
            message='Jtot must be positive.'
            call errormsg
         endif

!  -- Assign surface
         hopilab=12
         hopipar(1)=jj
         hopipar(2)=nk

!-----------------------------------------------------------------------

      else if (label(1:ll) .eq. 'readsrf'  .or. &
               label(1:ll) .eq. 'read1d'  ) then
!
! read options
!
         ascii=.false.
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   12       continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'ascii') then
               ascii=.true.
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binary') then
               ascii=.false.
            else
               filename=pesopts(np)(ilbl:jlbl)
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 12
         else
            routine='defsrf (opfuncs/funcsrf.F)'
            message='No argument (filename) given for '//label(1:ll)
            call errormsg
         endif
!
         ich = ich + 1  ! ich = channel for readsrf/read1d files 171...380
         length=index(filename,' ')-1
         call abspath(filename,length)
         if(ascii) then
            hopilab=13
            open(ich,file=filename,&
              form='formatted',status='old',err=200)
         else
            hopilab=14
            open(ich,file=filename,&
              form='unformatted',status='old',err=200)
         endif

!  -- Assign surface
         hopipar(1)=ich

!-----------------------------------------------------------------------
! Gauss 1D surface.
!-----------------------------------------------------------------------
      elseif (label(1:ll) .eq. 'gauss1d') then
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         width=1.0_dop
         hopipar(1)=0
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
  19        continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'width') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2              
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) width
               if (ierr .ne.  0) then
                  routine='defsrf (opfuncs/funcsrf.F)'
                  message='Cannot read argument "width" of gauss1D'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'periodic') then
               hopipar(1)=1
            else
               message='Unknown parameter for gauss1d: ' &
                        //pesopts(np)(ilbl:jlbl)
               routine='defsrf (opfuncs/funcsrf.F)'
               call errormsg
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 19
         endif
         hoprpar(1)=width
         hopilab=15

!-----------------------------------------------------------------------
! Gauss 2D surface.
!-----------------------------------------------------------------------
      elseif (label(1:ll) .eq. 'gauss2d') then
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         width=1.0_dop
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
            if (pesopts(np)(ilbl:jlbl) .eq. 'width') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) width
               if (ierr .ne.  0) then
                   routine='defsrf (opfuncs/funcsrf.F)'
                   message='Cannot read argument "width" of gauss2D'
                   go to 999
               endif
            else
               message='Unknown parameter for gauss2d: ' &
                        //pesopts(np)(ilbl:jlbl)
               routine='defsrf (opfuncs/funcsrf.F)'
               call errormsg
            endif
         endif
         hoprpar(1)=width
         hopilab=16

!-----------------------------------------------------------------------
! Coulomb 1D surface. 1/SQRT( (a*x1-b*x2 +c)**2 + d)
!-----------------------------------------------------------------------
      elseif (label(1:ll) .eq. 'coulomb1d') then
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         hoprpar(1)=1.0_dop
         hoprpar(2)=1.0_dop
         hoprpar(3)=0.0_dop
         hoprpar(4)=1.0_dop
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
 86      continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'a') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hoprpar(1)
               if (ierr .ne.  0) then
                   routine='defsrf (opfuncs/funcsrf.F)'
                   message='Cannot read argument "a" of coulomb1d'
                   go to 999
               endif
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'b') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hoprpar(2)
               if (ierr .ne.  0) then
                   routine='defsrf (opfuncs/funcsrf.F)'
                   message='Cannot read argument "b" of coulomb1d'
                   go to 999
               endif
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'c') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hoprpar(3)
               if (ierr .ne.  0) then
                   routine='defsrf (opfuncs/funcsrf.F)'
                   message='Cannot read argument "c" of coulomb1d'
                   go to 999
               endif
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'd') then
               ilbl=zpesopts(nopt+1,np)+1
               jlbl=zpesopts(nopt+2,np)
               nopt=nopt+2
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hoprpar(4)
               if (ierr .ne.  0) then
                   routine='defsrf (opfuncs/funcsrf.F)'
                   message='Cannot read argument "d" of coulomb1d'
                   go to 999
               endif
            else
               message='Unknown parameter for coulomb1d: ' &
                        //pesopts(np)(ilbl:jlbl)
               routine='defsrf (opfuncs/funcsrf.F)'
               call errormsg
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl.gt.0) go to 86
         endif
         hopilab=17

!
! Uwe Manthe NOCl
!
      else if (label(1:ll) .eq. 'nocl0um'  .or. label(1:ll) .eq. 'nocl1um') then
         hopilab = 18
!
! -- Set default options --
!
         icoord=1
!
! -- Read options --
!
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
            do while (jlbl > 0)
               if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
                  icoord=1
               else if (pesopts(np)(ilbl:jlbl) .eq. 'jac-cartesian') then
                  icoord=2
               else if (pesopts(np)(ilbl:jlbl) .eq. '3D-cartesian') then
                  icoord=3
               else if (pesopts(np)(ilbl:jlbl) .eq. '5D-cartesian') then
                  icoord=4
               else if (pesopts(np)(ilbl:jlbl) .eq. '9D-cartesian') then
                  icoord=5
               endif
               nopt=nopt+1
               ilbl=jlbl+1
               jlbl=zpesopts(nopt,np)
            end do
         endif

         if (label(1:ll) .eq. 'nocl1um') then
            hopipar(2)=1
         else
            hopipar(2)=0
         endif
         hopipar(1)=icoord
         lcutpes=.true.
         if (pesmin .lt. -99999.0_dop) pesmin=-1.0_dop/ev
         if (pesmax .gt.  99999.0_dop) pesmax=5.0_dop/ev

! ! -- Read data for the surface --
!          call innocl1sch
!
! -- Factor tfac  = 1 + mCl / (mN + mO) --
! -- Factor tfac1 = mO / (mN + mO)      --
! -- Factor tfac2 = mN / (mN + mO)      --
!
         tfac=2.165706579
         tfac1=0.5331995695
         tfac2=0.4668004305
         hoprpar(1)=tfac
         hoprpar(2)=tfac1
         hoprpar(3)=tfac2


!-----------------------------------------------------------------------
! Zero potential
!-----------------------------------------------------------------------
      else if (label .eq. 'zero_pot') then

         hopilab = 19   ! zero_pot

!-----------------------------------------------------------------------
! NH3 Surface from POTLIB site by Truhlar Group
! Z H Li, R Valero and D G Truhlar, Theor Chem Acc 118, 9-24(2007)
! 'Cartesian' or '12D' options for Cartesians
! 'normal' or 'NM' options for planar TS normal modes
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'nh3' .or. label(1:ll) .eq. 'nd3') then

         if (label(1:ll) .eq. 'nh3') then
            hoprpar(1) = 8
         else if (label(1:ll) .eq. 'nd3') then
            hoprpar(1) = 9
         endif

!  -  Set default options
         hopipar(1) = 4
         NH3isCart=.false.
         NH3isNM=.false.
         NH3map=(/1,2,3,4,5,6,6/)

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            rrd=0
            jlbl=zpesopts(1,np)
   97       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'U11') then
               hopipar(1) = 1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U22') then
               hopipar(1) = 2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U12') then
               hopipar(1) = 3
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V1') then
               hopipar(1) = 4
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V2') then
               hopipar(1) = 5
            else if (pesopts(np)(ilbl:jlbl) .eq. 'shepard') then
               hopipar(1) = 6
            endif
            if ((pesopts(np)(ilbl:jlbl)=='Cartesian') &
                .or.(pesopts(np)(ilbl:jlbl)=='cartesian') &
                .or.(pesopts(np)(ilbl:jlbl)=='12D') &
                .or.(pesopts(np)(ilbl:jlbl)=='12d')) NH3isCart=.true.
            if ((pesopts(np)(ilbl:jlbl)=='Normal') &
                .or.(pesopts(np)(ilbl:jlbl)=='normal') &
                .or.(pesopts(np)(ilbl:jlbl)=='NM') &
                .or.(pesopts(np)(ilbl:jlbl)=='nm')) NH3isNM=.true.
            if (rrd>0) then
              if (pesopts(np)(ilbl:jlbl)=='v1') then
                NH3map(rrd)=1
              elseif (pesopts(np)(ilbl:jlbl)=='v2') then
                NH3map(rrd)=2
              elseif (pesopts(np)(ilbl:jlbl)=='v3') then
                NH3map(rrd)=3
              elseif (pesopts(np)(ilbl:jlbl)=='v4') then
                NH3map(rrd)=4
              elseif (pesopts(np)(ilbl:jlbl)=='v5') then
                NH3map(rrd)=5
              elseif (pesopts(np)(ilbl:jlbl)=='v6') then
                NH3map(rrd)=6
              endif
              rrd=rrd+1
              if (rrd>NH3map(7)) rrd=-1
            endif
!Note that potfit feeds in blanks after pesopts but quantics doesn't.
            if ((len_trim(pesopts(np)(ilbl:jlbl))==2).and.(rrd==0)) then
              jj=jlbl
              if (pesopts(np)(jj:jj)==' ') jj=ilbl+1
              if ((pesopts(np)(jj:jj)=='d').or.(pesopts(np)(jlbl:jlbl)=='D')) then
                rrd=1
                read(pesopts(np)(ilbl:ilbl),*)NH3map(7)
              endif
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 97
         endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a)') '                                    '

           if (hopipar(1).eq.1) then
              write(ilog,'(a)') 'Diabatic energy 1 '
           else if (hopipar(1).eq.2) then
              write(ilog,'(a)') 'Diabatic energy 2 '
           else if (hopipar(1).eq.3) then
              write(ilog,'(a)') 'Coupling energy between two Diabats'
           else if (hopipar(1).eq.4) then
              write(ilog,'(a)') 'Adiabatic energy 1'
           else if (hopipar(1).eq.5) then
              write(ilog,'(a)') 'Adiabatic energy 2'
           else if (hopipar(1).eq.6) then
              write(ilog,'(a)') 'Shepard interpolation from POT'
              call defshep(i1,'shepard',i2,pesopts,zpesopts,ierr)
           endif
           if (NH3isCart) write(ilog,*)'12D Cartesian coordinates'
           if (NH3isNM) write(ilog,*)'Inversion TS normal mode coordinates'

         if (NH3isCart) then
           hopilab=25   ! Truhlar NH3  (Cartesian)
         elseif (NH3isNM) then
           hopilab=33   ! Truhlar NH3  (normal modes)
           open(unit=69,file='nh3map')
           write(69,*)NH3map
           close(69)
         else
           hopilab=22   ! Truhlar NH3
         endif
!-----------------------------------------------------------------------
! NH3 Surface from POTLIB site by Truhlar Group
! Z H Li, R Valero and D G Truhlar, Theor CHem Acc 118, 9-24(2007)
! in Valence coordinate
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'nh3valence') then

!  -  Set default options
         hopipar(1) = 1

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   98       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'U11') then
               hopipar(1) = 1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U22') then
               hopipar(1) = 2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U12') then
               hopipar(1) = 3
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V1') then
               hopipar(1) = 4
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V2') then
               hopipar(1) = 5
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 98
         endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a)') '                                    '

           if (hopipar(1).eq.1) then
              write(ilog,'(a)') 'Diabatic energy 1 '
           else if (hopipar(1).eq.2) then
              write(ilog,'(a)') 'Diabatic energy 2 '
           else if (hopipar(1).eq.3) then
              write(ilog,'(a)') 'Coupling energy between two Diabats'
           else if (hopipar(1).eq.4) then
              write(ilog,'(a)') 'Adiabatic energy 1'
           else if (hopipar(1).eq.5) then
              write(ilog,'(a)') 'Adiabatic energy 2'
           endif

         hopilab=23     !  Truhlar NH3 (valence)

!-----------------------------------------------------------------------
! NH3 Surface from POTLIB site by Truhlar Group
! Z H Li, R Valero and D G Truhlar, Theor CHem Acc 118, 9-24(2007)
! PES in Jacobi Cartesian Coordinates with COM at N
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'nh3jcart') then

!  -  Set default options
         hopipar(1) = 1

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
   99       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'U11') then
               hopipar(1) = 1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U22') then
               hopipar(1) = 2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'U12') then
               hopipar(1) = 3
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V1') then
               hopipar(1) = 4
            else if (pesopts(np)(ilbl:jlbl) .eq. 'V2') then
               hopipar(1) = 5
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 99
         endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a)') '                                    '

           if (hopipar(1).eq.1) then
              write(ilog,'(a)') 'Diabatic energy 1 '
           else if (hopipar(1).eq.2) then
              write(ilog,'(a)') 'Diabatic energy 2 '
           else if (hopipar(1).eq.3) then
              write(ilog,'(a)') 'Coupling energy between two Diabats'
           else if (hopipar(1).eq.4) then
              write(ilog,'(a)') 'Adiabatic energy 1'
           else if (hopipar(1).eq.5) then
              write(ilog,'(a)') 'Adiabatic energy 2'
           endif

         hopilab=24     ! Truhlar NH3 (Cartesian Jacobi Coordinates)

!-----------------------------------------------------------------------
! NH3 energy surfaces
!-----------------------------------------------------------------------
      else if (label(1:10) .eq. 'nh3_gwr_e_') then
!
! get which function
!C
         if(len(trim(label)).le.12)then
!(1:12) .eq. 'nh3_gwr_e_s0'.or.
!     &        label(1:12) .eq. 'nh3_gwr_e_s1') then
            call defnh3_gwr(hopipar(1),label(1:12))
         else
            call defnh3_gwr(hopipar(1),label(1:13))
         endif
         hopilab = 26

!-----------------------------------------------------------------------
! NH3 dipole and polarisability surfaces.
!-----------------------------------------------------------------------
      elseif (label(1:9).eq.'nh3_gwr_d') then

         if(len(trim(label)).le.13)then
            call defnh3_gwr(hopipar(1),label(1:13))
         else
            call defnh3_gwr(hopipar(1),label(1:14))
         endif
!         call defnh3_gwr(hopipar(1),label(1:13))
         hopilab = 27   ! GWR NH3 dipole

      elseif (label(1:9).eq.'nh3_gwr_p') then

         if(len(trim(label)).le.14)then
            call defnh3_gwr(hopipar(1),label(1:14))
         else
            call defnh3_gwr(hopipar(1),label(1:15))
         endif
!         call defnh3_gwr(hopipar(1),label(1:14))
         hopilab = 28   ! GWR NH3 polarisability

!-----------------------------------------------------------------------
! NOF surfaces Emma Chapman
!-----------------------------------------------------------------------
!!     else if (label(1:6) .eq. 'nof_ec') then
      else if (label(1:10) .eq. 'nof_ec_old') then

!  -  Set default options

         icoord = 2
         tfac  = 0.5_dop

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)

  94       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
               icoord=1
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 94
         endif
!
! get which function
!
         call defnof_ec(hopipar(1),label(1:ll))

         hopilab=29  ! NOF
         hopipar(2)=icoord

!-----------------------------------------------------------------------
! HCN Surface by Bowman
!-----------------------------------------------------------------------
      else if (label(1:ll)  .eq. 'hcn') then
         hopilab = 30  ! default
! read options
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            nopt=1
            ilbl=1
            jlbl=zpesopts(1,np)
            if (pesopts(np)(ilbl:jlbl) .eq. 'theta') then
               hopilab = 30
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'cos') then
               hopilab = 31
            else
               message= &
               'Only "theta" or "cos" are allowed as arguments to hcn'
               go to 999
            endif
         endif

!-----------------------------------------------------------------------
! Temporary hack to do nD SN2 backside attack modified Shepard
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'sn2') then
         hopilab=34

!  -  Set default options
         hopipar(1:13)=(/12,1,2,3,4,5,6,7,8,9,10,11,12/)

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            ilbl=1
            jlbl=zpesopts(1,np)
            read(pesopts(np)(ilbl:jlbl),*)hopipar(1)
            do i1=2,hopipar(1)+1
              read(pesopts(np)(zpesopts(i1-1,np)+1:zpesopts(i1,np)),*)hopipar(i1)
            enddo
         endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a)') '                                    '

           write(ilog,'(a)') 'Shepard interpolation from POT'
           write(ilog,'(i2,a)')hopipar(1),'D TS normal mode coordinates'
           call defshep(i1,'shepard',i2,pesopts,zpesopts,ierr)
           doforce=.false.

!-----------------------------------------------------------------------
! H2CO by Bowman J. Phys. Chem. A 2004, 108, 8980-8986
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h2co') then
         hopilab=35

!  -  Set default options
         hopipar(1:7)=(/6,1,2,3,4,5,6/)

!  -  Read option
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            ilbl=1
            jlbl=zpesopts(1,np)
            read(pesopts(np)(ilbl:jlbl),*)hopipar(1)
            do i1=2,hopipar(1)+1
              read(pesopts(np)(zpesopts(i1-1,np)+1:zpesopts(i1,np)),*)hopipar(i1)
            enddo
         endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a)') '                                    '

           write(ilog,'(i2,a)')hopipar(1),'D normal mode coordinates'


!-----------------------------------------------------------------------
! ADD NEW FUNCTION HERE
! !!! next free hopilab is 36 !!! 
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! Look for spline fit data hopilab=20,21 for 3D,2D  
! hopilab = -1 returned for 1D and this assigned at end to 34, ifile=303 
! i.e. to func1D             
!-----------------------------------------------------------------------
      else
         call defspline(label(1:ll),hopilab,hopipar,hoprpar)

!-----------------------------------------------------------------------
! end of if loop
!-----------------------------------------------------------------------
      endif

!-----------------------------------------------------------------------
! if label is assigned, assign ifile now.
!-----------------------------------------------------------------------
      if (hopilab .gt. 0) then
         ifile=201
      else if (hopilab .eq. -1) then
! assign 1D splines to func1d
         hopilab = 17
         ifile = 303
      endif


      return

  999 continue
      routine='defsrf (opfuncs/funcsrf.F)'
      call errormsg
      return

  200 continue
      message='Cannot open file '//filename(1:length)
      routine='defsrf (opfuncs/funcsrf.F)'
      call errormsg
      return

      end subroutine defsrf

!#######################################################################
!
! search spline fits
!
!#######################################################################

      subroutine defspline(hoplab,hopilab,hopipar,hoprpar)

use global
use runopermod

      implicit none

      integer(long)                :: ilbl,jlbl,n,nspl
      integer(long), intent(out)   :: hopilab
      integer(long), dimension(*), intent(out) :: hopipar
      real(dop), dimension(*), intent(out)     :: hoprpar
      character(len=*), intent(in) :: hoplab
      character(len=c2)            :: label1
      integer(long), external      :: slen
      logical(kind=4)              :: lexist
      character(len=c5)            :: filename,fileopts
      character                    :: cc

!-----------------------------------------------------------------------
! No spline fits for direct dynamics calculations
!-----------------------------------------------------------------------
      if (ldd) return

      hopilab = 0
      hopipar(1)=0

!-----------------------------------------------------------------------
! search list of spline fits. Format of name is (i5,a) XXXXXname.ndspl
! where XXXXX defines the directory and nd the dimensionality
!
! see if label is in the list of spline files 
!-----------------------------------------------------------------------
      n=1
      nspl=0
      do while (n .le. nsplnames)
         ilbl=slen(splnames(n))
         if (hoplab .eq. splnames(n)(6:ilbl-6)) then
            read(splnames(n)(1:5),'(i5)') nspl
            exit
         endif
         n=n+1
      enddo

!-----------------------------------------------------------------------
! Check for a spline in the namedir if splnames has not been set up
! properly.  This is a hack for potfit.
!-----------------------------------------------------------------------
      if (nsplpaths==-1) then
        cc=''
        inquire(file=name(1:laenge)//'/'//hoplab//'.1dspl',exist=lexist)
        if (lexist) cc='1'
        inquire(file=name(1:laenge)//'/'//hoplab//'.2dspl',exist=lexist)
        if (lexist) cc='2'
        inquire(file=name(1:laenge)//'/'//hoplab//'.3dspl',exist=lexist)
        if (lexist) cc='3'
        if (.not.(cc=='')) then
          splpaths(1)=name(1:laenge)
          splnames(1)='XXXXX'//hoplab//'.'//cc//'dspl'
          splunits='au'
          n=1
          nspl=1
        endif
      endif

!-----------------------------------------------------------------------
! return if not found
! if it is, hopipar(1) set to index of splpath (path where file is stored)
! hopipar(2) set to index of splname 
!-----------------------------------------------------------------------
      if (nspl .eq. 0) return
      hopipar(1) = nspl
      hopipar(2) = n

!-----------------------------------------------------------------------
! test file exists
!-----------------------------------------------------------------------
      lexist=.false.
      ilbl=slen(splnames(n))
      jlbl=slen(splpaths(nspl))
      filename=splpaths(nspl)(1:jlbl)//'/'//splnames(n)(6:ilbl)
      inquire(file=filename,exist=lexist)

      if (.not. lexist) then
         routine = 'Defsplines (funcsrf1)'
         write(message,'(2a)') 'Spline file not found : ',filename(1:ilbl+jlbl)
         call errormsg
      endif

!-----------------------------------------------------------------------
! define function depending on dimensionality
!-----------------------------------------------------------------------
      label1 = splnames(n)(ilbl-4:ilbl-3)
      call keylowc(label1)
      if (label1 .eq. '3d') then
         hopilab = 20
      else if (label1 .eq. '2d') then
         hopilab = 21
      else if (label1 .eq. '1d') then
         hopilab = -1
      endif

!-----------------------------------------------------------------------
! get options / parameters
! set defaults
! first look in splpaths for defaults
! then look in file header
!
! hoprpar(1) : vzero (shift)
! hoprpar(2) : vmax (cutoff)
!-----------------------------------------------------------------------
      hoprpar(1) = 0.0
      hoprpar(2) = 9999.0/ev
      lexist=.false.
      fileopts=splpaths(nspl)(1:jlbl)//'/spline.opt '
      inquire(file=fileopts,exist=lexist)
      if (lexist) call getsplopts(fileopts,hoprpar,splunits(1,n))

      call getsplopts(filename,hoprpar,splunits(1,n))

      return
      end subroutine defspline

!#######################################################################
!
! filename is name of file - either splpath/spline.opts or splpath/splname.spl
! hoprpar: parameters for operator
! splunits(n) : units for dofs / energy in splien data
! 
! If any new options are added, these should be documented in the
! HOWTO on spline fits doc/howtos/spline_fits.html
!#######################################################################

      subroutine getsplopts(filename,hoprpar,splunits)

use channels
use global

      real(dop),dimension(*),intent(out)           :: hoprpar
      integer(long)                                :: i,j,ierr,ic,ilbl,jlbl
      integer(long), external                        :: slen
      character(*), intent(in)                     :: filename
      character(len=c1), dimension(4), intent(out) :: splunits
      real(dop)                                    :: vmax,vzero
      character(len=c4)                            :: line
      integer(long),parameter                      :: max_nof_tokens=30
      integer(long)                                :: nof_tokens
      integer(long), dimension(max_nof_tokens,2)   :: aof_tbounds
      character(len=c1),dimension(max_nof_tokens)  :: token

      open(idump,file=filename,status='old',err=1000)

      lines: do
! read line
         read(idump,'(a)',end=999) line
         if (line(1:1) .ne. '#') exit

! break up into tokens (keywords)
         call get_tokens(line,max_nof_tokens,nof_tokens,aof_tbounds)
         do i=1,nof_tokens
            ilbl=aof_tbounds(i,1)
            jlbl=aof_tbounds(i,2)
            token(i)=line(ilbl:jlbl)
            call keylowc(token(i))
         enddo

! search tokens and store options
         i=0
         ic = nof_tokens
         tokens: do while (i .le. ic)
            i=i+1
            if (token(i) .eq. 'units') then
               if (token(i+1) .ne. '=') then
                  routine='Getsplnames (einopermod)'
                  write(message,'(a)') &
                   'No units given in Options file: ',filename
                  call errormsg
               endif
               j = 0
               do
                  i=i+2
                  j=j+1
                  splunits(j) = token(i)(1:c1)
                  if (token(i+1) .ne. ',') exit
               enddo

            else if (token(i) .eq. 'vzero') then
               if (token(i+1) .ne. '=') then
                  routine='Getsplnames (einopermod)'
                  write(message,'(a)') &
                   'No ezero given in Options file: ',filename
                  call errormsg
               endif

               i=i+2

               read(token(i),*) vzero
               if (token(i+1) .eq. ',') then
                  i=i+2
                  call keyunits(vzero,token(i),ierr)
                  if (ierr .eq. 1) then
                     routine='Getsplnames (einopermod)'
                     message = 'no unit after comma for vzero'
                     call errormsg
                  endif
               endif
               hoprpar(1)=vzero

            else if (token(i) .eq. 'vmax') then
               if (token(i+1) .ne. '=') then
                  routine='Getsplnames (einopermod)'
                  write(message,'(a)') &
                   'No ezero given in Options file: ',filename
                  call errormsg
               endif

               i=i+2
               read(token(i),*) vmax
               if (token(i+1) .eq. ',') then
                  i=i+2
                  call keyunits(vmax,token(i),ierr)
                  if (ierr .eq. 1) then
                     routine='Getsplnames (einopermod)'
                     message = 'no unit after comma for vmax'
                     call errormsg
                  endif
               endif
               hoprpar(2)=vmax

            endif
         enddo tokens

      enddo lines

  999 continue
! Close file
      close(idump)

      return

 1000 jlbl=slen(filename)
      routine='Getsplopts (funcsrf1)'
      message='Error opening file: '//filename(1:jlbl)
      call errormsg

      end subroutine getsplopts

!#######################################################################
! This routine calculates the value of a potential at a grid point.
! The point on the MCTDH grid is contained in gpoint. 
!
! To add new surfaces, add the new ifunc number, and then any code
! needed to manipulate the grid point and call the potential surface
! routine. For user specific routines please use uvpoint (and udefsrf)
! on usersrf.F  
!
!#######################################################################

      subroutine vpoint (v,gpoint,ndof1,ifunc,hopipar,hoprpar,label)

use channels
      implicit none

      integer(long)                           :: imod,icoll,jj,ich, &
                                                 i,j,n,icoord
      integer(long), dimension(*), intent(in) :: hopipar
      integer(long), intent(in)               :: ifunc,ndof1
      real(dop)                               :: r1,r2,r3,rd,rv,theta,r,phi, &
                                                 th1,th2,tfac,tfac1,tfac2,&
                                                 rho,alf,&
                                                 xO,yO,zO,xN,yN,zN,xCl,yCl,&
                                                 zCl,rdz,rvx,rvz,u

      real(dop), dimension(12)               :: gU11,gU22,gU12,gV11,gV22,xc
      real(dop), dimension(5)                :: UV
      real(dop)                              :: U11,U22,U12,V11,V22,mh,mn,&
                                                cbe1,cbe2,gamma,a_r2,a_r3,&
                                                b_r2,b_r3,rms,totsum


      real(dop), intent(out)                  :: v
      real(dop), dimension(*), intent(in)     :: hoprpar
      real(dop), dimension(3)                 :: a,b,c,d
      real(dop), dimension(ndof1), intent(in) :: gpoint
      real(dop), dimension(6)                 :: gp
      real(dop), dimension(12)                :: gp12
      real(dop), dimension(18)                :: x18
      character(len=*), intent(in)            :: label
      logical,dimension(12)                   :: m12
      logical,dimension(6)                    :: m6
      integer(long)                           :: natom
      real(dop), dimension(36)                :: xmass
      real(dop), dimension(36,3)              :: xx
      real(dop), parameter                    :: amu2au=1822.888486192d0
      real(dop), parameter                    :: hcm=219474.63d0
      real(dop)                               :: v1,v2,s_h2co
      integer(long), save                     :: iz = 0
      real(dop), save                         :: vold
      character(len=16)                       :: vvv    

     select case (ifunc)

     case default
      routine='vpoint (opfuncs/funcsrf.F)'
      write(message,'(a,i5)') 'Function not found :',ifunc
      call errormsg

     case (1)
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      call nocl1um(v,rd,rv,theta)
      return

     case (2)
      imod=0
      icoll=1
      rd=gpoint(1)
      rv=gpoint(2)

      tfac=hoprpar(1)
      r1 = max(0.4_dop,rd-tfac*rv)
      r2 = max(0.4_dop,rv)
      r3 = r1+r2
      call lsth(r1,r2,r3,v,imod,icoll)
      return

     case (3)
      imod=0
      icoll=1
      r1=gpoint(1)
      r2=gpoint(2)

      r1 = max(0.4_dop,r1)
      r2 = max(0.4_dop,r2)
      r3 = r1+r2
      call lsth(r1,r2,r3,v,imod,icoll)
      return

     case (4)
      imod=0
      icoll=0
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

!     x=rv**2/4_dop+rd**2
!     y=rd*rv*cos(theta)

!     r1 = max(0.4_dop,sqrt(x-y))
!     r2 = max(0.4_dop,rv)
!     r3 = max(0.4_dop,sqrt(x+y))

      tfac=hoprpar(1)
      r1=rd**2+(tfac*rv)**2-2.0_dop*tfac*rd*rv*cos(theta)
      r1=sqrt(r1)
      r2=rv
      r3=rd**2+((1.0_dop-tfac)*rv)**2+2.0_dop*(1.0_dop-tfac)*rd*rv*cos(theta)
      r3=sqrt(r3)

      r1 = max(0.4_dop,r1)
      r2 = max(0.4_dop,r2)
      r3 = max(0.4_dop,r3)

      call lsth(r1,r2,r3,v,imod,icoll)
      return

     case (5)
      imod=0
      icoll=0
      r1=gpoint(1)
      r2=gpoint(2)
      phi=gpoint(3)

      r3 = sqrt(r1**2+r2**2-2*r1*r2*cos(phi))

      r1 = max(0.4_dop,r1)
      r2 = max(0.4_dop,r2)
      r3 = max(0.4_dop,r3)

      call lsth(r1,r2,r3,v,imod,icoll)
      return

     case (6)
      imod=0
      icoll=0
      rho=gpoint(1)
      alf=gpoint(2)
      theta=gpoint(3)
       
      r1=(rho*tan(alf))**2/(1+tan(alf)**2)
      r1=sqrt(r1)
      r2=rho**2/(1.0_dop+tan(alf)**2)
      r2=sqrt(r2)
      r3=r1**2+r2**2-2.0_dop*r1*r2*cos(theta)
      r3=sqrt(r3)
      
      r1=max(0.4_dop,r1)
      r2=max(0.4_dop,r2)
      r3=max(0.4_dop,r3)
      call lsth(r1,r2,r3,v,imod,icoll)                            
      return
      
     case (7)
      imod=0
      icoll=0
      rho=gpoint(1)
      alf=0.8_dop
      theta=gpoint(3)
       
      r1=(rho*tan(alf))**2/(1+tan(alf)**2)
      r1=sqrt(r1)
      r2=rho**2/(1.0_dop+tan(alf)**2)
      r2=sqrt(r2)
      r3=r1**2+r2**2-2.0_dop*r1*r2*cos(theta)
      r3=sqrt(r3)
      
      r1=max(0.4_dop,r1)
      r2=max(0.4_dop,r2)
      r3=max(0.4_dop,r3)
      call lsth(r1,r2,r3,v,imod,icoll)                            
      return

     case (8)
      return

     case (9)
      r1 =gpoint(1)
      r2 =gpoint(2)
      r3 =gpoint(3)
      th1=gpoint(4)
      th2=gpoint(5)
      phi=gpoint(6)
 
      call coortransb(r1,r2,r3,th1,th2,phi,b,a,c,d)
      call yzcl(a,b,c,d,v)
      return
 
     case (10)
      R  =gpoint(1)
      r1 =gpoint(2)
      r2 =gpoint(3)
      th1=gpoint(4)
      th2=gpoint(5)
      phi=gpoint(6)
 
      tfac1=hoprpar(1)
      tfac2=hoprpar(2)
 
      call coortransj(R,r1,r2,th1,th2,phi,tfac1,tfac2,b,a,c,d)
      call yzcl(a,b,c,d,v)
      return

     case (11)!.....cpp
      r=gpoint(1)
      do j=2,hopipar(2)
         r=r+gpoint(j)
      enddo
      r=r*(r+1.d0)
      jj=hopipar(1)
      if (r .lt. dble(jj*(jj+1))) then
         v=sqrt(dble(jj*(jj+1))-r)
      else
         v=0.d0
      endif
      return

     case (12)!.....cmm
      r=gpoint(1)
      do j=2,hopipar(2)
         r=r+gpoint(j)
      enddo
      r=r*(r-1.d0)
      jj=hopipar(1)
      if (r .lt. dble(jj*(jj+1))) then
         v=sqrt(dble(jj*(jj+1))-r)
      else
         v=0.d0
      endif
      return

     case (13)!.....readsrf (ascii)
      ich=hopipar(1)
      iz=iz+1
      read(ich,*,err=200,end=210) v
      vold=v
      return

     case (14)!.....readsrf (binary)
      ich=hopipar(1)
      iz=iz+1
      read(ich,err=200,end=210) v
      vold=v
      return

     case (15)!.....Normalized 1D Gaussian.  
      r2 = gpoint(1)-gpoint(2)
      if(hopipar(1) .eq. 1) then
         if( r2 .gt.  pi ) r2 = r2 - 2.0_dop*pi
         if( r2 .lt. -pi ) r2 = r2 + 2.0_dop*pi
      endif
      r2 = 0.5_dop*(r2/hoprpar(1))**2
      v  = exp(-r2)/(SQRT(2.0_dop*pi)*hoprpar(1))
      return

     case (16)!.....Normalized 2D Gaussian.  
      r2 = 0.5_dop*((gpoint(1)-gpoint(3))**2 + (gpoint(2)-gpoint(4))**2) &
                /hoprpar(1)**2
      v  = exp(-r2)/(2.0_dop*pi*hoprpar(1)**2)
      return

     case (17)!.....Regularized 1D Coulomb.  
      r2 = hoprpar(1)*gpoint(1)-hoprpar(2)*gpoint(2)+hoprpar(3)
      v  = 1.0_dop/SQRT(r2*r2 + hoprpar(4))
      return

   case (18)
      icoord=hopipar(1)
      if (icoord .eq. 1) then
         rd=gpoint(1)
         rv=gpoint(2)
         theta=gpoint(3)

      else if (icoord .eq. 2) then
         !.....NOCl S0 PES from McDonald in the 3D Cartesian coordinates
         !.....of the Jacobi vectors in the Jacobi-embedded BF frame
         rdz = gpoint(1)
         rvx = gpoint(2)
         rvz = gpoint(3)
         !.....Transformation Cartesian -> Jacobian
         rd = ABS(rdz)
         rv = SQRT(rvx**2+rvz**2)
         theta = ACOS(rvz/rv)

      else if (icoord .eq. 3) then
         !.....NOCl S0 PES from McDonald in the 3D-full-constrained Cartesian
         !.....coordinates of the nuclei in the Jacobi-embedded BF frame
         xO = gpoint(1)
         zO = gpoint(2)
         zCl= gpoint(3)
         tfac = hoprpar(1)
         tfac2 = hoprpar(3)
         !.....Transformation Cartesian -> Jacobian
         rd = tfac*ABS(zCl)
         rv = SQRT(xO**2+(zO+(tfac-1)*zCl)**2)/tfac2
         theta = ACOS((zO+(tfac-1)*zCl)/tfac2/rv)

      else if (icoord .eq. 4) then
         !.....NOCl S0 PES from McDonald in the 5D-semi-constrained Cartesian
         !.....coordinates of the nuclei in the 2D(NO)/1D(Cl) frame
         xO = gpoint(1)
         zO = gpoint(2)
         xN = gpoint(3)
         zN = gpoint(4)
         zCl= gpoint(5)
         tfac = hoprpar(1)
         tfac1 = hoprpar(2)
         tfac2 = hoprpar(3)
         !.....Transformation Cartesian -> Jacobian
         rd = SQRT((tfac1*xO+tfac2*xN)**2+(zCl-tfac1*zO-tfac2*zN)**2)
         rv = SQRT((xO-xN)**2+(zO-zN)**2)
         theta = ACOS(((zCl-tfac1*zO-tfac2*zN)*(zO-zN) &
              -(tfac1*xO+tfac2*xN)*(xO-xN))/rd/rv)

      else if (icoord .eq. 5) then
         !.....NOCl S0 PES from McDonald in the 9D free Cartesian
         !.....coordinates of the nuclei in the SF frame
         xO = gpoint(1)
         yO = gpoint(2)
         zO = gpoint(3)
         xN = gpoint(4)
         yN = gpoint(5)
         zN = gpoint(6)
         xCl= gpoint(7)
         yCl= gpoint(8)
         zCl= gpoint(9)
         tfac = hoprpar(1)
         tfac1 = hoprpar(2)
         tfac2 = hoprpar(3)
         !.....Transformation Cartesian -> Jacobian
         rd = SQRT((xCl-tfac1*xO-tfac2*xN)**2+ &
              (yCl-tfac1*yO-tfac2*yN)**2+ &
              (zCl-tfac1*zO-tfac2*zN)**2)
         rv = SQRT((xO-xN)**2+(yO-yN)**2+(zO-zN)**2)
         theta = ACOS(((xCl-tfac1*xO-tfac2*xN)*(xO-xN) &
              +(yCl-tfac1*yO-tfac2*yN)*(yO-yN) &
              +(zCl-tfac1*zO-tfac2*zN)*(zO-zN))/rd/rv)
      else if (icoord .eq. 6) then
         !.. Jacobi coordinates with the third coordinate utheta=cos(theta)
         rd=gpoint(1)
         rv=gpoint(2)
         theta=ACOS(gpoint(3))
      endif

      if (hopipar(2) == 0) then
         call nocl0um(rv,rd,theta,v)
      else
         call nocl1um(v,rd,rv,theta)
      end if

      return

     case (19)! Zero potential
      v = 0.0_dop
      return

     case (20)! 3D Spline Fits
      call splinepot3d(v,gpoint,label,hopipar,hoprpar)
      return

     case (21)! 2D Spline Fits
      call splinepot2d(v,gpoint,label,hopipar,hoprpar)
      return

     case (22)! 6D NH3 in Jacobians
      r1    = gpoint(1)
      r2    = gpoint(2)
      r3    = gpoint(3)
      gamma = gpoint(4)
      phi   = gpoint(5)
      theta = gpoint(6)

! 11     r1    = gpoint(1)
!
      mh = 999.0
      if (hoprpar(1) .eq. 8) then
         mh = 1.007825_dop
!         if (hopipar(1).eq.1) then
!          r2    = 1.96_dop
!          r3    = 1.92_dop
!          gamma = 1.883_dop
!          theta = 1.915_dop
!         else if (hopipar(1).eq.2) then
!          r2 = 2.032_dop
!          r3 = 1.96_dop
!          gamma = 2.09_dop
!          theta = 2.151_dop
!         else if (hopipar(1).eq.3) then
!          r2 = 1.977_dop
!          r3 = 1.929_dop
!          gamma = 2.179_dop
!          theta = 1.981_dop
!         endif
      else if (hoprpar(1) .eq. 9) then
          mh = 2.014_dop
!         if (hopipar(1).eq.1) then
!          r2    = 2.0_dop
!          r3    = 1.92_dop
!          gamma = 1.90_dop
!          theta = 1.968_dop
!         else if (hopipar(1).eq.2) then
!          r2 = 2.098_dop
!          r3 = 1.96_dop
!          gamma = 2.09_dop
!          theta = 2.196_dop
!         else if (hopipar(1).eq.3) then
!          r2 = 2.025_dop
!          r3 = 1.929_dop
!          gamma = 2.175_dop
!          theta = 2.032_dop
!         endif
      endif
!c
         mn = 14.0067_dop

!cccc 3D
!      gamma = gpoint(2)
!      phi   = gpoint(3)
!cccc 2D
!      phi   = gpoint(2)

!  Jacobi -> Cartesian -> potential routine coordinates

      b_r3=r3*(mn/(mn+mh))
      a_r3=r3-b_r3
      a_r2=(r2*mh)/(mh+mn+mh)
      b_r2=r2-a_r2

! ATOM 4(c)=N (x,y,z)
      Xc(1)=-a_r3*COS(theta)-a_r2
      Xc(2)=-a_r3*SIN(theta)
      Xc(3)=0.0_dop
! ATOM 3(a)=H (x,y,z)
      Xc(4)=b_r3*COS(theta)-a_r2
      Xc(5)=b_r3*SIN(theta)
      Xc(6)=0.0_dop
! ATOM 2(d)=H (x,y,z)
      Xc(7)=b_r2
      Xc(8)=0.0_dop
      Xc(9)=0.0_dop
! ATOM 1(b)=H (x,y,z)
      Xc(10)=r1*COS(gamma)
      Xc(11)=(r1*SIN(gamma)*COS(phi))
      Xc(12)=(r1*SIN(gamma)*SIN(phi))

      if (.not.NH3init) then
        call potc_nh3_init
        NH3init=.true.
      endif

      call potc_nh3(Xc,U11,U22,U12,V11,V22,&
                   gU11,gU22,gU12,gV11,gV22)


      if (hopipar(1).eq.1) then
         v = U11
      else if (hopipar(1).eq.2) then
         v = U22
      else if (hopipar(1).eq.3) then
         v = U12
      else if (hopipar(1).eq.4) then
         v = V11
      else if (hopipar(1).eq.5) then
         v = V22
      endif

      return
 
     case (23)! 6D NH3 in Valence Coordinate
      r1    = gpoint(1)
      r2    = gpoint(2)
      r3    = gpoint(3)
      cbe1  = gpoint(4)
      cbe2  = gpoint(5)
      phi   = gpoint(6)

! Valence -> Cartesian -> potential routine coordinates
! ATOM 4(c) = N (x,y,z)
      Xc(1) = 0.0_dop
      Xc(2) = 0.0_dop
      Xc(3) = 0.0_dop
! ATOM 3(d) = H3 (x,y,z)
      Xc(4) = 0.0_dop
      Xc(5) = 0.0_dop
      Xc(6) = r3
! ATOM 2(a) = H2 (x,y,z)
      Xc(7) = r2 * (1.0 - cbe2*cbe2)**0.5
      Xc(8) = 0.0_dop
      Xc(9) = r2 * cbe2
! ATOM 1(b) = H1 (x,y,z)
!      Xc(10) = r1 * cbe1
!      Xc(11) = (r1 * (1.0 - cbe1*cbe1)**0.5 * SIN(phi))
!      Xc(12) = (r1 * (1.0 - cbe1*cbe1)**0.5 * COS(phi))

      Xc(10) = (r1 * (1.0 - cbe1*cbe1)**0.5 * COS(phi))
      Xc(11) = (r1 * (1.0 - cbe1*cbe1)**0.5 * SIN(phi))
      Xc(12) = r1 * cbe1

!      Xc(11) = (r1 * (1.0 - cbe2*cbe2)**0.5 * COS(phi))
!      Xc(12) = (r1 * (1.0 - cbe2*cbe2)**0.5 * SIN(phi))

      if (.not.NH3init) then
        call potc_nh3_init
        NH3init=.true.
      endif

      call potc_nh3(xc,U11,U22,U12,V11,V22,gU11,gU22,gU12,gV11,gV22)

      if (hopipar(1).eq.1) then
         v = U11
      else if (hopipar(1).eq.2) then
         v = U22
      else if (hopipar(1).eq.3) then
         v = U12
      else if (hopipar(1).eq.4) then
         v = V11
      else if (hopipar(1).eq.5) then
         v = V22
      endif

      return
 
     case (24)! 6D NH3 in Jacobi Cartesian Coordinates with COM at N
      xc(4)  = gpoint(1)
      xc(7)  = gpoint(2)
      xc(8)  = gpoint(3)
      xc(10) = gpoint(4)
      xc(11) = gpoint(5)
      xc(12) = gpoint(6)

      xc(1) = 0.0_dop
      xc(2) = 0.0_dop
      xc(3) = 0.0_dop

      xc(5) = 0.0_dop
      xc(6) = 0.0_dop
      xc(9) = 0.0_dop

      mh = 1.007825_dop
      mn = 14.0067_dop
      xc(10) = xc(10) + (mh /(mh + mn))*xc(4) + (mh/(2.0_dop*mh+mn))*xc(7)
      xc(11) = xc(11) + (mh /(2.0_dop*mh + mn))*xc(8)

      xc(7) = xc(7) + (mh /(mh + mn))*xc(4)

      if (.not.NH3init) then
        call potc_nh3_init
        NH3init=.true.
      endif

      call potc_nh3(xc,U11,U22,U12,V11,V22,gU11,gU22,gU12,gV11,gV22)

      if (hopipar(1).eq.1) then
         v = U11
      else if (hopipar(1).eq.2) then
         v = U22
      else if (hopipar(1).eq.3) then
         v = U12
      else if (hopipar(1).eq.4) then
         v = V11
      else if (hopipar(1).eq.5) then
         v = V22
      endif

      return

     case (25)  ! 12D cartesian NH3
      if (.not.NH3init) then
        call potc_nh3_init
        NH3init=.true.
      endif
      call potc_nh3(gpoint,xc(1),xc(2),xc(3),xc(4),xc(5),&
                   gU11,gU22,gU12,gV11,gV22)
      v=xc(hopipar(1))
      return

     case (26)! NH3 energies
      call func_nh3_gwr(v,gpoint,hopipar(1),hopipar(2))
      return

     case (27)! NH3 dipoles.
      call func_nh3_gwr(v,gpoint,hopipar(1),hopipar(2))
      return

     case (28)! NH3 polarisabilities.
      call func_nh3_gwr(v,gpoint,hopipar(1),hopipar(2))
      return

     case (29)! NOF
      call func_nof_ec(v,gpoint,hopipar(1),hopipar(2))
      return

     case(30) ! HCN surface with u=cos(theta)
      r1=gpoint(1)
      r2=gpoint(2)
      u=gpoint(3)
      call hcn(v,r1,r2,u)
      return

     case(31) ! HCN surface with theta
      r1=gpoint(1)
      r2=gpoint(2)
      th1=gpoint(3)
      call hcn1(v,r1,r2,th1)
      return

     case (32) ! LSTH surface with 3D Cartesian rep of Jacobian coords
      imod=0
      icoll=0

      r1=sqrt(gpoint(1)**2+(gpoint(2)+2*gpoint(3))**2)
      r2=sqrt(4*gpoint(1)**2+(2*gpoint(2)+gpoint(3))**2)
      r3=sqrt(gpoint(1)**2+(gpoint(3)-gpoint(1))**2)

      r1 = max(0.4_dop,r1)
      r2 = max(0.4_dop,r2)
      r3 = max(0.4_dop,r3)

      call lsth(r1,r2,r3,v,imod,icoll)
      return

     case (33)  ! NH3, normal modes of planar inversion TS
      if ((hopipar(1)<6).and..not.NH3init) then
        call potc_nh3_init
        open(unit=69,file='nh3map')
        read(69,*)NH3map
        close(69)
        NH3init=.true.
      endif
      gp=0
      gp(NH3map(1:NH3map(7)))=gpoint(1:NH3map(7))
       ! TS geom
      xc=0
      xc(4)=1.88968815_8
      xc(7:8)=1.88968815_8*(/-0.5_8,sqrt(3.0_8)/2.0_8/)
      xc(10:11)=1.88968815_8*(/-0.5_8,-sqrt(3.0_8)/2.0_8/)
       ! add NM contributions
      xc((/3,6,9,12/))=xc((/3,6,9,12/))+gp(1)/100.0_8 &
                *(/-0.26370490245146216_8,1.2214664578977195_8,1.2214664578977195_8,1.2214664578977195_8/)
      xc=xc+gp(2)/100.0_8*(/-0.16273634160599129_8,5.7551903532009507E-002_8,0.0_8, &
                        -0.20619421867846482_8,-0.60607527528225114_8,0.0_8, &
                        0.93976170941952819_8,0.73453879924069843_8,0.0_8, &
                        1.5277895897479574_8,-0.92819518492059760_8,0.0_8/)
      xc=xc+gp(3)/100.0_8*(/-5.7551903572669927E-002_8,-0.16273634183557750_8,0.0_8, &
                        -7.2920834638153686E-002_8,1.7137656061729178_8,0.0_8, &
                        1.2676932392510134_8,0.56780967695923978_8,0.0_8, &
                        -0.39504074359242641_8,-2.0218203176711032E-002_8,0.0_8/)
      xc(4)=xc(4)-gp(4)/100.0_8*1.3468804690168004_8
      xc((/7,10/))=xc((/7,10/))+gp(4)/100.0_8*0.673440236_8
      xc((/8,11/))=xc((/8,11/))+gp(4)/100.0_8*(/-1,1/)*1.166432703_8
      xc=xc+gp(5)/100.0_8*(/-1.4036260641945086E-002_8,0.19886658958432066_8,0.0_8, &
                        0.12708714988220579_8,-4.1699565929883095E-002_8,0.0_8, &
                        0.79559608081280819_8,-1.4146146122646188_8,0.0_8, &
                        -0.72763769217016838_8,-1.3071028068985187_8,0.0_8/)
      xc=xc+gp(6)/100.0_8*(/-0.19886658965071452_8,-1.4036260638368282E-002_8,0.0_8, &
                        1.8005784248158148_8,2.9432092925098120E-003_8,0.0_8, &
                        0.42766337707681007_8,-0.66556572129947045_8,0.0_8, &
                        0.53517518306214962_8,0.85766805034581273_8,0.0_8/)
      if (hopipar(1)/=6) then
        call potc_nh3(xc,UV(1),UV(2),UV(3),UV(4),UV(5),&
                   gU11,gU22,gU12,gV11,gV22)
        v=UV(hopipar(1))
      else
        call vshepard(v,xc,12,rms,totsum)
        if (lsheptraj) write(itmpin,*)gpoint,v,rms,totsum
      endif
      return

     case (34)  ! SN2 (CH3F + Cl-), normal modes of TS
      gp12=0
      gp12(hopipar(2:hopipar(1)+1))=gpoint(1:ndof1-1)
      m12=.false.
      m12(hopipar(2:hopipar(1)+1))=.true.
       ! This is based on MP2/aug-cc-pVTZ
       ! TS geom
      x18=0
      x18((/10,13/))=1.742172144984748_8*(/-1,1/)
      x18(8)=2.011687113763213_8
      x18((/11,14/))=-1.005843556881606_8
      x18(3)=-0.976365651573553_8
      x18(6)=-4.751627963929321_8
      x18(9:15:3)=-1.176536570232136_8
      x18(18)=3.018313517099574_8
       ! add NM contributions
      if (m12( 1)) x18=x18+gp12( 1)*(/ &
             -0.148981287776009E-12_8,-0.109023267277396E-09_8,-0.606211512618424E-02_8, &
             -0.283465958036520E-13_8, 0.279912505677514E-10_8, 0.195475147704614E-02_8, &
              0.181803165547612E-12_8, 0.547482469979973E-03_8, 0.470657885621496E-03_8, &
             -0.474133996154392E-03_8,-0.273741167547902E-03_8, 0.470657320534490E-03_8, &
              0.474133995216815E-03_8,-0.273741168350803E-03_8, 0.470657319204778E-03_8, &
              0.883072432205173E-13_8, 0.183408843519226E-10_8, 0.977590080633163E-03_8/)
      if (m12( 2)) x18=x18+gp12( 2)*(/ &
             -0.520803625986115E-02_8,-0.201831515193173E-06_8, 0.507490860752270E-12_8, &
              0.218060523843415E-02_8, 0.845068514851418E-07_8, 0.141087490790589E-11_8, &
             -0.552474235212917E-02_8,-0.210976280365149E-06_8,-0.355848818020024E-07_8, &
             -0.546418852141769E-02_8, 0.347454432789509E-04_8, 0.795231746917154E-03_8, &
             -0.546419123108857E-02_8,-0.351720889741612E-04_8,-0.795196161505022E-03_8, &
              0.107668062045852E-02_8, 0.417255421616968E-07_8,-0.955957618147215E-12_8/)
      if (m12( 3)) x18=x18+gp12( 3)*(/ &
             -0.201831493876129E-06_8, 0.520803626844216E-02_8,-0.540055565003516E-10_8, &
              0.845068675340529E-07_8,-0.218060521745431E-02_8, 0.147660441948212E-09_8, &
             -0.214112505753786E-06_8, 0.544400647373165E-02_8, 0.918233982638800E-03_8, &
             -0.351711308158806E-04_8, 0.550455914688191E-02_8,-0.459086197042389E-03_8, &
              0.347476211647373E-04_8, 0.550455644968282E-02_8,-0.459147832246071E-03_8, &
              0.417255313530019E-07_8,-0.107668063381177E-02_8,-0.603460244369892E-10_8/)
      if (m12( 4)) x18=x18+gp12( 4)*(/ &
              0.191592437842684E-11_8,-0.163859736550234E-09_8, 0.879935564136868E-03_8, &
             -0.753344260956841E-12_8, 0.340526411683359E-10_8, 0.404994597825256E-02_8, &
              0.439556198051789E-13_8,-0.347687139024910E-04_8, 0.517943654275254E-03_8, &
              0.301104631653463E-04_8, 0.173841194840843E-04_8, 0.517943597401008E-03_8, &
             -0.301104573037856E-04_8, 0.173841228324527E-04_8, 0.517943598177710E-03_8, &
             -0.418386095490316E-12_8, 0.513213052772781E-10_8,-0.254705877270918E-02_8/)
      if (m12( 5)) x18=x18+gp12( 5)*(/ &
              0.826666034385362E-03_8, 0.131578845183824E-07_8, 0.225202640414975E-12_8, &
              0.250957146014550E-03_8, 0.399443933112708E-08_8,-0.585216441160496E-13_8, &
             -0.215152789847197E-02_8,-0.156137155501062E-07_8,-0.297435461181154E-06_8, &
             -0.127359173442750E-02_8, 0.506842433479647E-03_8, 0.161834582738921E-01_8, &
             -0.127360787003285E-02_8,-0.506901608522645E-03_8,-0.161831608382600E-01_8, &
             -0.284604196492217E-03_8,-0.452999133186307E-08_8,-0.504118245143361E-13_8/)
      if (m12( 6)) x18=x18+gp12( 6)*(/ &
              0.131579069435417E-07_8,-0.826666027317947E-03_8, 0.869361914333929E-10_8, &
              0.399443237104438E-08_8,-0.250957139017178E-03_8,-0.162251738069886E-09_8, &
             -0.342455567748841E-07_8, 0.980957396753426E-03_8, 0.186868778340266E-01_8, &
             -0.506892484838898E-03_8, 0.185889308657173E-02_8,-0.934317898290559E-02_8, &
              0.506851941672201E-03_8, 0.185887695112692E-02_8,-0.934369415632275E-02_8, &
             -0.452999625686723E-08_8, 0.284604192238952E-03_8,-0.769898171541305E-10_8/)
      if (m12( 7)) x18=x18+gp12( 7)*(/ &
             -0.665484338249237E-14_8, 0.227995748108841E-10_8,-0.398769525188194E-03_8, &
             -0.159933453796659E-14_8, 0.193720521356177E-10_8,-0.836913247947599E-03_8, &
             -0.439497071817864E-13_8, 0.116679083773005E-02_8, 0.131018263530465E-01_8, &
             -0.101047107822235E-02_8,-0.583395161269302E-03_8, 0.131018299348724E-01_8, &
              0.101047107832165E-02_8,-0.583395161197327E-03_8, 0.131018299351127E-01_8, &
              0.155735640281664E-14_8,-0.331989252800097E-10_8,-0.541276082379388E-03_8/)
      if (m12( 8)) x18=x18+gp12( 8)*(/ &
              0.746861334604103E-06_8, 0.203623761780262E-02_8,-0.101847587604415E-09_8, &
             -0.134830270192914E-07_8,-0.367600278461477E-04_8,-0.122367750291229E-10_8, &
             -0.654815588240497E-05_8, 0.204324614733401E-02_8, 0.265542547783681E-02_8, &
              0.861418598328443E-02_8,-0.128756638416161E-01_8,-0.132855529530379E-02_8, &
             -0.861633591311507E-02_8,-0.128819837377489E-01_8,-0.132686832990545E-02_8, &
              0.171491328623874E-08_8, 0.467553197264594E-05_8,-0.117954745555190E-10_8/)
      if (m12( 9)) x18=x18+gp12( 9)*(/ &
             -0.203623762786098E-02_8, 0.746861400518525E-06_8,-0.343906205377883E-13_8, &
              0.367600317554995E-04_8,-0.134830429877338E-07_8, 0.829429039254012E-14_8, &
              0.178528470835092E-01_8, 0.749431724628296E-06_8, 0.973970188381993E-06_8, &
              0.293393708151045E-02_8,-0.861998470491087E-02_8, 0.229917846865258E-02_8, &
              0.292761718520222E-02_8, 0.861053718731640E-02_8,-0.230015243810940E-02_8, &
             -0.467552828713256E-05_8, 0.171490386141436E-08_8,-0.137886410271226E-13_8/)
      if (m12(10)) x18=x18+gp12(10)*(/ &
              0.103170652270431E-12_8, 0.324972655678385E-09_8,-0.284512344096856E-03_8, &
             -0.177299929688807E-15_8,-0.182635418201617E-11_8,-0.351438225034719E-05_8, &
             -0.486756162183226E-13_8,-0.134081285794992E-01_8, 0.115801204353633E-02_8, &
              0.116117767239552E-01_8, 0.670406245912580E-02_8, 0.115801138063312E-02_8, &
             -0.116117767251360E-01_8, 0.670406245944759E-02_8, 0.115801138077655E-02_8, &
              0.124659791023066E-15_8,-0.501593008004489E-11_8,-0.580616306659551E-06_8/)
      if (m12(11)) x18=x18+gp12(11)*(/ &
              0.172137035607097E-04_8, 0.226528733096790E-02_8, 0.597987627095880E-10_8, &
             -0.109648315014095E-06_8,-0.144294888367725E-04_8,-0.119008944028452E-11_8, &
              0.197619531720735E-05_8,-0.178647928803422E-01_8, 0.170386797100106E-02_8, &
             -0.794961594675688E-02_8,-0.433079068134412E-02_8,-0.863147169801635E-03_8, &
              0.774697472299399E-02_8,-0.421151380302867E-02_8,-0.840721385732835E-03_8, &
             -0.642302888538035E-07_8,-0.845257261889062E-05_8,-0.302749255251370E-11_8/)
      if (m12(12)) x18=x18+gp12(12)*(/ &
             -0.226528732995351E-02_8, 0.172137035394839E-04_8, 0.443044638901506E-12_8, &
              0.144294886568518E-04_8,-0.109648313491268E-06_8, 0.219328541689631E-14_8, &
             -0.260063174534766E-03_8,-0.135752867044052E-03_8, 0.129475311072232E-04_8, &
              0.132739418323233E-01_8, 0.781583761470168E-02_8, 0.146911933562458E-02_8, &
              0.133932187350050E-01_8,-0.788074977566231E-02_8,-0.148206687107667E-02_8, &
              0.845257155930045E-05_8,-0.642302951070966E-07_8,-0.280061603045535E-13_8/)
        call vshepard(v,x18,18,rms,totsum)
        if (lsheptraj) write(itmpin,*)gpoint,v,rms,totsum

      return

     case (35)  ! H2CO, normal modes min(xe)
!    H   0.0000000000000000        1.7935587134812330        3.4022052420726365     
!    O   0.0000000000000000        0.0000000000000000       -1.1149325371329212E-002
!    C   0.0000000000000000        0.0000000000000000        2.2771079783816446     
!    H   0.0000000000000000       -1.7935587134812330        3.4022052420726365
      natom=4
      xmass(1)=1.00794d0*amu2au
      xmass(2)=15.9994d0*amu2au
      xmass(3)=12.0107d0*amu2au
      xmass(4)=xmass(1)
      gp=0
      gp(hopipar(2:hopipar(1)+1))=gpoint(1:ndof1-1)
      m6=.false.
      m6(hopipar(2:hopipar(1)+1))=.true.
      xc=0
      xc((/2,11/))=1.7935587134812330_8*(/1,-1/)
      xc((/3,12/))=3.4022052420726365_8
      xc(6)=-1.1149325371329212E-02_8
      xc(9)=2.2771079783816446_8
       ! add NM contributions
       ! L-1CL=w^2
       ! mwq=LQ Q;normal mode coordinate
       ! mwq=sqrt(mass)*(xc-xe)
       ! (xc-xe)=M-1LQ M-1=1/sqrt(mass)
       ! The below matrix is M-1L
      if (m6( 1)) xc=xc+gp( 1)*(/ &
          -1.3923729400233877E-002_8, 2.5864108925088192E-020_8,-3.7175984008226607E-020_8, &
          -8.6186057845532497E-004_8, 5.4207274736354631E-021_8, 7.4937394448412666E-022_8, &
           3.4850442212372244E-003_8,-1.5448174724765660E-020_8, 1.1680976797028083E-021_8, &
          -1.3923729404232916E-002_8,-1.0233991604494739E-020_8, 9.1368201899463010E-021_8/)
      if (m6( 2)) xc=xc+gp( 2)*(/ &
          -8.0273166498587166E-020_8, 5.1341632492475494E-003_8,-1.3087638347630336E-002_8, &
           1.0756212260911306E-018_8, 1.6108114868904990E-003_8, 3.9504064783415436E-013_8, &
           3.4300354936632149E-019_8,-3.0074744454436568E-003_8,-1.2428276919283222E-013_8, &
           1.8181196739092853E-019_8, 5.1341632500636622E-003_8, 1.3087638347169375E-002_8/)
      if (m6( 3)) xc=xc+gp( 3)*(/ &
           8.4503927301724110E-019_8,-6.5906228660033254E-003_8, 1.2096188138732462E-002_8, &
          -4.3882301611999422E-019_8, 4.2831877726550387E-014_8,-2.8407343530058265E-003_8, &
          -5.0255208221171419E-020_8,-2.6168303654000186E-013_8, 1.7539013659286994E-003_8, &
          -2.9544244533668541E-018_8, 6.5906228668014682E-003_8, 1.2096188139193727E-002_8/)
      if (m6( 4)) xc=xc+gp( 4)*(/ &
          -9.6969075151832615E-019_8, 7.2157738742053262E-003_8,-5.3327697514571441E-003_8, &
           3.0690671248686198E-019_8, 2.4640989540470349E-014_8,-2.8165692060788242E-003_8, &
          -1.7984428104877017E-019_8,-1.9720995166565368E-015_8, 4.6469930211708371E-003_8, &
           5.3063304833048711E-018_8,-7.2157738736256952E-003_8,-5.3327697496006058E-003_8/)
      if (m6( 5)) xc=xc+gp( 5)*(/ &
           1.0686796776890006E-018_8,-1.3219303354129415E-002_8,-8.2380774230118008E-003_8, &
          -7.9427134343421207E-019_8,-5.1890661941718952E-006_8,-6.7726883940278911E-014_8, &
           5.3891956632829536E-020_8, 2.2256447400222485E-003_8, 5.2410951080142914E-013_8, &
          -2.3717272516651444E-018_8,-1.3219303344903521E-002_8, 8.2380774179591776E-003_8/)
      if (m6( 6)) xc=xc+gp( 6)*(/ &
          -1.8006307996201875E-018_8, 1.3290054395131199E-002_8, 8.8939797102369021E-003_8, &
           1.2578034741003223E-018_8,-9.6399297740052256E-015_8, 1.2050493459058443E-004_8, &
          -1.9130423425616706E-019_8, 6.6578707041412296E-013_8,-1.6532926275126883E-003_8, &
           7.9091602444586302E-019_8,-1.3290054403506505E-002_8, 8.8939797162615637E-003_8/)
      n=0
      do i=1,natom
        do j=1,3
          n=n+1
          xx(i,j)=xc(n)
        enddo
      enddo
!      do i=1,natom
!        write(1000,*) i,(xx(i,j),j=1,3)
!      enddo
      call pot_h2co(natom,xmass,xx,v,v1,v2,s_h2co)
      v=v/hcm
!      write(6,*)"pot",v
      return

   end select

 200  write(6,'(a,i5)') ' Error reading readsrf-file from channel:',ich
      write(2,'(a,i5)') ' Error reading readsrf-file from channel:',ich
      write(6,'(a,i5,a)') ' Error at line',iz,' of readsrf-file.'
      write(2,'(a,i5,a)') ' Error at line',iz,' of readsrf-file.'
      backspace ich
      if(ifunc.eq.13) read(ich,'(a)') vvv
      if(ifunc.eq.14) read(ich) vvv
      write(6,'(1x,a,f15.9,2x,a)') &
        'Before last and last lines read:',vold,trim(vvv)
      write(2,'(1x,a,f15.9,2x,a)') &
        'Before last and last lines read:',vold,trim(vvv)
      call stopnow

 210  write(6,'(a,i5)') ' Unexpected EOF for readsrf-file, channel:',ich
      write(2,'(a,i5)') ' Unexpected EOF for readsrf-file, channel:',ich
      write(6,'(a,i5,a)') ' There are only',iz-1,' lines on readsrf-file.'
      write(2,'(a,i5,a)') ' There are only',iz-1,' lines on readsrf-file.'
      write(6,'(a,f15.9)') ' Last energy values read:',vold
      write(2,'(a,f15.9)') ' Last energy values read:',vold
      write(6,'(a,4e12.4)') ' gpoint:',gpoint(1),gpoint(2),gpoint(3),gpoint(4)
      write(2,'(a,4e12.4)') ' gpoint:',gpoint(1),gpoint(2),gpoint(3),gpoint(4)
      call stopnow


      end subroutine vpoint

!##########################################################################

      subroutine splinepot3d(v,gpoint,label,hopipar,hoprpar)

use runopermod

      implicit none

      integer(long)                           :: n,nspl
      integer(long)                           :: icoord,ilbl,jlbl,ll,ierr
      integer(long), dimension(*), intent(in) :: hopipar
      real(dop)                               :: r1,r2,alpha,rv,rd,b,arg,&
                                                 theta,rtmp
      real(dop)                               :: vmax,vzero,tfac
      real(dop), intent(out)                  :: v
      real(dop), dimension(3), intent(in)     :: gpoint
      real(dop), dimension(*), intent(in)     :: hoprpar
      character(len=120)                      :: filename
      character(len=*), intent(in)            :: label

!-----------------------------------------------------------------------
! get splinepath index
!-----------------------------------------------------------------------
      nspl = hopipar(1)
      n = hopipar(2)
      vzero = hoprpar(1)
      vmax  = hoprpar(2)
      icoord = 1
      tfac = 0.5_dop

!-----------------------------------------------------------------------
! coordinate unit trafo
!-----------------------------------------------------------------------
      if (icoord .le. 1) then
         r1 = gpoint(1)
         r2 = gpoint(2)
         alpha = gpoint(3)

      else if (icoord .eq. 2) then
         rd = gpoint(1)
         rv = gpoint(2)
         theta = gpoint(3)
!
! Transformation auf Molekularkoordinaten (r1,r2,alpha), in denen
! die Potentialdaten gegeben sind.

         r2 = rv
         b = tfac*rv
         r1=dsqrt(b*b+rd*rd+2.0_dop*b*rd*dcos(theta))
         arg=(b**2+r1**2-rd**2)/(2.0_dop*b*r1)
         if(arg.gt.1.0_dop) then
            alpha=0.0_dop
         else if(arg.lt.-1.0_dop) then
            alpha=3.141592654_dop
         else
            alpha=dacos(arg)
         end if

         alpha = alpha

      endif

      rtmp=1.0_dop
      call keylunits(rtmp,splunits(1,n),ierr)
      r1=r1/rtmp
      rtmp=1.0_dop
      call keylunits(rtmp,splunits(2,n),ierr)
      r2=r2/rtmp
      rtmp=1.0_dop
      call keylunits(rtmp,splunits(3,n),ierr)
      alpha=alpha/rtmp

!-----------------------------------------------------------------------
! calculate spline fit
!-----------------------------------------------------------------------
      ilbl=index(label,'{')-1
      jlbl=index(label,' ')-1
      if(ilbl.le.0) ilbl=99
      if(jlbl.le.0) jlbl=99
      ll = min(ilbl,jlbl)

      ilbl=index(splpaths(nspl),' ')-1
      filename=splpaths(nspl)(1:ilbl)//'/'//label(1:ll)//'.3dspl'
      call splinepot3d_css(r1,r2,alpha,v,filename)
      call keyeunits(v,splunits(4,n),ierr)

!-----------------------------------------------------------------------
! add energy zero and cut off high values
!-----------------------------------------------------------------------
      v = v+vzero
      v = min(v,vmax)

      return
      end subroutine splinepot3d

!##########################################################################

      subroutine splinepot2d(v,gpoint,label,hopipar,hoprpar)

use runopermod

      implicit none

      integer(long)                           :: n,nspl
      integer(long)                           :: icoord,ilbl,jlbl,ll,ierr
      integer(long), dimension(*), intent(in) :: hopipar
      real(dop)                               :: r1,r2,rv,rd,b,arg,theta,rtmp
      real(dop)                               :: vmax,vzero,tfac
      real(dop), intent(out)                  :: v
      real(dop), dimension(3), intent(in)     :: gpoint
      real(dop), dimension(*), intent(in)     :: hoprpar
      character(len=120)                      :: filename
      character(len=*), intent(in)            :: label

!-----------------------------------------------------------------------
! get splinepath index
!-----------------------------------------------------------------------
      nspl = hopipar(1)
      n = hopipar(2)
      vzero = hoprpar(1)
      vmax  = hoprpar(2)
      icoord = 1
      tfac = 0.5_dop

!-----------------------------------------------------------------------
! coordinate unit trafo
!-----------------------------------------------------------------------
      if (icoord .le. 1) then
         r1 = gpoint(1)
         r2 = gpoint(2)

      else if (icoord .eq. 2) then
         rd = gpoint(1)
         rv = gpoint(2)
         theta = 0.0_dop
!
! Transformation auf Molekularkoordinaten (r1,r2,alpha), in denen
! die Potentialdaten gegeben sind.

         r2 = rv
         b = tfac*rv
         r1=dsqrt(b*b+rd*rd+2.0_dop*b*rd*dcos(theta))
         arg=(b**2+r1**2-rd**2)/(2.0_dop*b*r1)
      endif

      rtmp=1.0_dop
      call keylunits(rtmp,splunits(1,n),ierr)
      r1=r1/rtmp
      rtmp=1.0_dop
      call keylunits(rtmp,splunits(2,n),ierr)
      r2=r2/rtmp

!-----------------------------------------------------------------------
! calculate spline fit
!-----------------------------------------------------------------------
      ilbl=index(label,'{')-1
      jlbl=index(label,' ')-1
      if(ilbl.le.0) ilbl=99
      if(jlbl.le.0) jlbl=99
      ll = min(ilbl,jlbl)

      ilbl=index(splpaths(nspl),' ')-1
      filename=splpaths(nspl)(1:ilbl)//'/'//label(1:ll)//'.2dspl'
      call splinepot2d_css(r1,r2,v,filename)
      call keyeunits(v,splunits(3,n),ierr)

!-----------------------------------------------------------------------
! add energy zero and cut off high values
!-----------------------------------------------------------------------
      v = v+vzero
      v = min(v,vmax)

      return 
      end subroutine splinepot2d

!##########################################################################
!      SURFINFO
! This routine prints information on the surfaces to log file.
! It is called from getpes and htmsetup.
! For user implemented surfaces use the routine usersurfinfo on usersrf.F
!##########################################################################

      subroutine surfinfo (h,hl,done)

use global
use maxv, only: maxdim,maxhop
use dvrdatmod
use griddatmod
use operdef
use runopermod
use channels

use funcqchemmod, only:qcheminfo
use funcqch1mod, only:qchem1info
use funcshep

implicit none

      integer(long)                              :: i,j,muld,mode,f,ilbl,&
                                                    jlbl,np,nn,icoord,ll,n,&
                                                    nspl
      integer(long)                              :: slen
      integer(long),intent(in)                   :: h,hl
      integer(long), dimension(2*maxhop), intent(inout) :: done
!integer(long), dimension(:), allocatable, save :: done
      character(len=16), dimension(maxdim)       :: mlabel 
      character(len=c2)                          :: label
      character(len=80)                          :: filename
 
!      logical(kind=4), save                      :: initial=.true.

      integer(long), save :: nz = 0

!      if (initial) then
!         allocate(done(maxhop))
!         done(:) = 0
!         initial=.false.
!      endif

      if(nz.eq.0) write(ilog,'(//a)') &  
        '********* MULTIDIMENSIONAL POTENTIALS *********'
      do i=1,nz
         if(done(i).eq.hl) goto 520
      enddo
      goto 510
 520  if(hopfile(hl).ne.201) return
      if(hopilab(hl).eq.36.or.hopilab(hl).eq.37) rewind(hopipar(1,hl))
      return

      label=hoplab(hl)
 510  nz=nz+1
      done(nz)=hl
      muld=htmmuld(h)
      mode=htmmode(h)
      if(muld.gt.0) then
         nn=nmulddof(muld) 
         do i = 1, nn
            f = mulddof2(i,muld)
            if( f.eq.0 ) then
                mlabel(i) = 'UNDEFINED'
            else
                mlabel(i) = modelabel(f)
            endif
         enddo
      elseif(mode.gt.0) then
         nn=nspfdof(mode)
         do i = 1, nn
            f = spfdof(i,mode)
            if( f.eq.0 ) then
                mlabel(i) = 'UNDEFINED'
            else
                mlabel(i) = modelabel(f)
            endif
         enddo
      else
         return
      endif
      if(hopilab(hl).eq.36 .or. hopilab(hl).eq.37) then
         if(hoplab(hl)(1:6).eq.'read1d') then
            routine='surfinfo (opfuncs/funcsrf.F)'
            message='read1d for muld? Do you mean readsrf?'
            call warnmsg
         endif 
      endif 

      ilbl=index(hoplab(hl),' ')-1
      if(hopidef(hl).gt.0) then
         np=hopidef(hl)
         jlbl=index(deflabel(np),' ')-1
         write(ilog,'(/5a)') '-------  ',deflabel(np)(1:jlbl), &
                      '  =  ',hoplab(hl)(1:ilbl),'  ---------' 
       else
          write(ilog,'(/3a)')  &
          '-------  ',hoplab(hl)(1:ilbl),'  ---------'
       endif

      if( hopfile(hl) .eq. 201 ) then
hopilab201:select case (hopilab(hl))

        case default 
         routine='surfinfo (opfuncs/funcsrf.F)'
         write(message,'(a,i5)') 'Surface not found :',hopilab(hl)
         call errormsg

        case (1)
         write(ilog,'(a)') 'NOCl S1 PES from U. Manthe '
         write(ilog,'(a)') 'Jacobian Coordinates '
         write(ilog,'(2a)') 'Dissociative Coordinate  :',mlabel(1)
         write(ilog,'(2a)') 'Vibrational Coordinate   :',mlabel(2)
         write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(3)

        case (2)
         write(ilog,'(a)') &
             'LSTH PES : collinear model in Jacobian Coordinates'
         write(ilog,'(2a)')'Dissociative Coordinate : ',mlabel(1)
         write(ilog,'(2a)')'Vibrational Coordinate  : ',mlabel(2)

        case (3)
         write(ilog,'(a)') &
              'LSTH PES : collinear model in Binding Coordinates'
         write(ilog,'(2a)') 'Bond 1     : ',mlabel(1)
         write(ilog,'(2a)') 'Bond 2     : ',mlabel(2)

        case (4)
         write(ilog,'(2a,f9.4)') 'LSTH PES :', &
            '3D model in Jacobian Coordinates. tfac =',hoprpar(1,hl)
         write(ilog,'(2a)')'Dissociative Coordinate : ',mlabel(1)
         write(ilog,'(2a)')'Vibrational Coordinate  : ',mlabel(2)
         write(ilog,'(2a)')'Angular Coordinate      : ',mlabel(3)

        case (5)
         write(ilog,'(a)') &
              'LSTH PES : 3D model in Binding Coordinates'
         write(ilog,'(2a)') 'Bond 1     : ',mlabel(1)
         write(ilog,'(2a)') 'Bond 2     : ',mlabel(2)
         write(ilog,'(2a)') 'Bond Angle : ',mlabel(3)


        case (6)
         write(ilog,'(a)') &
               'LSTH PES : 3D model in Hyperspherical Coordinates'
         write(ilog,'(2a)') 'Hyper-radius : ',mlabel(1)
         write(ilog,'(2a)') 'Hyper-angle  : ',mlabel(2)
         write(ilog,'(2a)') 'bond Angle   : ',mlabel(3)

        case (7)
         write(ilog,'(a)') &
               'LSTH PES : 2D model in Hyperspherical Coordinates'
         write(ilog,'(2a)') 'Hyper-radius : ',mlabel(1)
         write(ilog,'(2a)') 'Hyper-angle at TST'
         write(ilog,'(2a)') 'Bond Angle : ',mlabel(2)


        case (8)
         return


        case (9)
         write(ilog,'(a,i1,2a)') &
             'YZCL',hopipar(1,hl),' PES in Binding Coordinates ', &
             '(r1,r2,r3,th1,th2,phi)'
         write(ilog,'(2a)')'Vibrational Coordinate r1=r(H1-O)  : ', &
                                mlabel(1)
         write(ilog,'(2a)')'Dissociative Coordinate r2=r(O-H2) : ', &
                                mlabel(2)
         write(ilog,'(2a)')'Vibrational Coordinate r3=r(H2-H3) : ', &
                                mlabel(3)
         write(ilog,'(2a)')'Angular Coordinate th1=(H1-O-H2)   : ', &
                                mlabel(4)
         write(ilog,'(2a)')'Angular Coordinate th2=(O-H2-H3)   : ', &
                                mlabel(5)
         write(ilog,'(2a)')'Torsional Angular Coordinate (phi) : ', &
                                mlabel(6)

        case (10)
         write(ilog,'(a,i1,2a)') &
             'YZCL',hopipar(1,hl),' PES in Jacobian Coordinates ', &
             '(R,r1,r2,th1,th2,phi)'
         write(ilog,'(2a)')'Scattering coordinate  (R)         : ', &
                                mlabel(1)
         write(ilog,'(2a)')'First Vibrational Coordinate (r1)  : ', &
                                mlabel(2)
         write(ilog,'(2a)')'Second Vibrational Coordinate (r2) : ', &
                                mlabel(3)
         write(ilog,'(2a)')'First Angular Coordinate  (th1)    : ', &
                                mlabel(4)
         write(ilog,'(2a)')'Second Angular Coordinate  (th2)   : ', &
                                mlabel(5)
         write(ilog,'(2a)')'Torsional Angular Coordinate (phi) : ', &
                                mlabel(6)
         write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1,hl), &
                                     '  tfac2 = ',hoprpar(2,hl)

        case (11)
         write(ilog,'(a,i3)') 'C_+ Operator term for J=',hopipar(1,hl)
         write(ilog,'(2a)') ' First coordinate  : ',mlabel(1)
         write(ilog,'(2a)') ' Second coordinate : ',mlabel(2)

        case (12)
         write(ilog,'(/a,i3)') 'C_- Operator term for J=',hopipar(1,hl)
         write(ilog,'(2a)') ' First  coordinate : ',mlabel(1)
         write(ilog,'(2a)') ' Second coordinate : ',mlabel(2)

        case (13)
         write(ilog,'(a,i4,2x,a)') &
               'readsrf, channel =',hopipar(1,hl),'ASCII'
         do i=1,nn
            write(ilog,'(i3,2a)') i,'. Variable : ',mlabel(i)
            if(muld.gt.0) then
               if(mulddof2(i,muld).ne.mulddof(i,muld)) then
                  write(6,'(/a,30i4)') 'mulddof2:', &
                       (mulddof2(j,muld),j=1,nmulddof(muld))
                  write(6,'(a,30i4)') 'mulddof :', &
                       (mulddof(j,muld),j=1,nmulddof(muld))
                  write(6,'(a,30(2x,a5))') 'Use order :', &
                   (modelabel(mulddof(j,muld)),j=1,nmulddof(muld))
                  write(ilog,'(/a,/30(2x,a5))')'Use order for readsrf:', &
                   (modelabel(mulddof(j,muld)),j=1,nmulddof(muld))
                  routine='surfinfo (opfuncs/funcsrf.F)' 
                  message='readsrf does not allow a re-ordering of '// &
                  'the DOFs. The DOFs must hence be correctly ordered'// &
                  ' (SPF-basis ordering) via a |i&j&k construct.' 
                  call errormsg
               endif
            endif
         enddo
         if(nn.eq.1) then
            routine='surfinfo (opfuncs/funcsrf.F)'
            message='readsrf for 1D? Do you mean read1d?'
            call warnmsg
         endif

        case (14)
         write(ilog,'(a,i4,2x,a)') &
               'readsrf, channel =',hopipar(1,hl),'binary'
         do i=1,nn
            write(ilog,'(i3,2a)') i,'. Variable : ',mlabel(i)
            if(muld.gt.0) then
               if(mulddof2(i,muld).ne.mulddof(i,muld)) then
                  write(6,'(/a,30i4)') 'mulddof2:', &
                       (mulddof2(j,muld),j=1,nmulddof(muld))
                  write(6,'(a,30i4)') 'mulddof :', &
                       (mulddof(j,muld),j=1,nmulddof(muld))
                  write(6,'(a,30(2x,a5))') 'Use order :', &
                   (modelabel(mulddof(j,muld)),j=1,nmulddof(muld))
                  write(ilog,'(/a,/30(2x,a5))')'Use order for readsrf:', &
                   (modelabel(mulddof(j,muld)),j=1,nmulddof(muld))
                  routine='surfinfo (opfuncs/funcsrf.F)' 
                  message='readsrf does not allow a re-ordering of '// &
                  'the DOFs. The DOFs must hence be correctly ordered'// &
                  ' (SPF-basis ordering) via a |i&j&k construct.' 
                  call errormsg
               endif
            endif
         enddo
         if(nn.eq.1) then
            routine='surfinfo (opfuncs/funcsrf.F)'
            message='readsrf for 1D? Do you mean read1d?'
            call warnmsg
         endif

        case (15)
         write(ilog,'(a,f12.6)') 'Normalized 1D Gaussian,  sigma =', &
                                   hoprpar(1,hl)
         write(ilog,'(2a)') 'x1 : ',mlabel(1)
         write(ilog,'(2a)') 'x2 : ',mlabel(2)

        case (16)
         write(ilog,'(a,f12.6)') 'Normalized 2D Gaussian,  sigma =', &
                                   hoprpar(1,hl)
         write(ilog,'(2a)') 'x1 : ',mlabel(1)
         write(ilog,'(2a)') 'y1 : ',mlabel(2)
         write(ilog,'(2a)') 'x2 : ',mlabel(3)
         write(ilog,'(2a)') 'y2 : ',mlabel(4)

        case (17)
         write(ilog,'(a,/4(a,f12.6))')  &
         'Regularized 1D Coulomb, V=1/SQRT((a*x1-b*x2+c)**2 + d)', &
           'a =',hoprpar(1,hl),',  b =',hoprpar(2,hl), &
         '  c =',hoprpar(3,hl),',  d =',hoprpar(4,hl)
         write(ilog,'(2a)') 'x1 : ',mlabel(1)
         write(ilog,'(2a)') 'x2 : ',mlabel(2)

       case (18)
         icoord=hopipar(1,hl)
         write(ilog,'(a)') 'NOCl S0 PES from McDonald '
         if (icoord .eq. 1) then
            write(ilog,'(a)') 'Jacobian Coordinates '
            write(ilog,'(2a)') 'Dissociative Coordinate : ',mlabel(1)
            write(ilog,'(2a)') 'Vibrational Coordinate  : ',mlabel(2)
            write(ilog,'(2a)') 'Angular Coordinate      : ',mlabel(3)
         else if (icoord .eq. 2) then
            write(ilog,'(a)') 'Jacobi-Cartesian Coordinates (BF frame) '
            write(ilog,'(2a)') 'z_rd >= 0 : ',mlabel(1)
            write(ilog,'(2a)') 'x_rv >= 0 : ',mlabel(2)
            write(ilog,'(2a)') 'z_rv      : ',mlabel(3)
         else if (icoord .eq. 3) then
            write(ilog,'(a)') '3D Cartesian Coordinates (BF frame) '
            write(ilog,'(2a)') 'x_O >= 0  : ',mlabel(1)
            write(ilog,'(2a)') 'z_O       : ',mlabel(2)
            write(ilog,'(2a)') 'z_Cl >= 0 : ',mlabel(3)
         else if (icoord .eq. 4) then
            write(ilog,'(a)') '5D Cartesian Coordinates (2D/1D frame) '
            write(ilog,'(2a)') 'x_O >= 0  : ',mlabel(1)
            write(ilog,'(2a)') 'z_O       : ',mlabel(2)
            write(ilog,'(2a)') 'x_N =< 0  : ',mlabel(3)
            write(ilog,'(2a)') 'z_N       : ',mlabel(4)
            write(ilog,'(2a)') 'z_Cl >= 0 : ',mlabel(5)
         else if (icoord .eq. 5) then
            write(ilog,'(a)') '9D Cartesian Coordinates (SF frame) '
            write(ilog,'(2a)') 'x_O  : ',mlabel(1)
            write(ilog,'(2a)') 'y_O  : ',mlabel(2)
            write(ilog,'(2a)') 'z_O  : ',mlabel(3)
            write(ilog,'(2a)') 'x_N  : ',mlabel(4)
            write(ilog,'(2a)') 'y_N  : ',mlabel(5)
            write(ilog,'(2a)') 'z_N  : ',mlabel(6)
            write(ilog,'(2a)') 'x_Cl : ',mlabel(7)
            write(ilog,'(2a)') 'y_Cl : ',mlabel(8)
            write(ilog,'(2a)') 'z_Cl : ',mlabel(9)
         endif

       case (19)
         write(ilog,'(a)') 'Zero potential '

       case (20)
         label=hoplab(hl)
         nspl=hopipar(1,hl)
         n=hopipar(2,hl)
         if (allocated(splpaths)) then
            ilbl=index(splpaths(nspl),' ')-1
            ll=slen(label)
            filename=splpaths(nspl)(1:ilbl)//'/'//label(1:ll)//'.3dspl'
            write(ilog,'(a)') '3D Spline fit to: '
            write(ilog,'(a)') filename
            write(ilog,'(a,4a10)') 'Units : ',(splunits(i,n),i=1,4)
            write(ilog,'(2(a,f12.6,a,5x))') 'Vzero : ',hoprpar(1,hl)*ev,' eV', &
                                            'Vmax  : ',hoprpar(2,hl)*ev,' eV'
         else
            write(ilog,'(a)') '3D Spline fit'
         endif

       case (21)
         label=hoplab(hl)
         n=hopipar(1,hl)
         if (allocated(splpaths)) then
            ilbl=index(splpaths(n),' ')-1
            ll=slen(label)
            filename=splpaths(n)(1:ilbl)//'/'//label(1:ll)//'.2dspl'
            write(ilog,'(a)') '2D Spline fit to: '
            write(ilog,'(a)') filename
            write(ilog,'(a,4a10)') 'Units : ',(splunits(i,n),i=1,3)
            write(ilog,'(2(a,f12.6,a,5x))') 'Vzero : ',hoprpar(1,hl)*ev,' eV', &
                                            'Vmax  : ',hoprpar(2,hl)*ev,' eV'
         else
            write(ilog,'(a)') '2D Spline fit'
         endif

       case (22)
        write(ilog,'(a)') 'NH3 Surface :'
        write(ilog,'(a)') 'in Jacobian Coordinates '
        write(ilog,'(2a)') 'Translational Coordinate :',mlabel(1)
        write(ilog,'(2a)') 'Dissociative Coordinate  :',mlabel(2)
        write(ilog,'(2a)') 'Vibrational Coordinate   :',mlabel(3)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(4)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(5)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(6)

       case (23)
        write(ilog,'(a)') 'NH3 Surface :'
        write(ilog,'(a)') 'in Valence Coordinates '
        write(ilog,'(2a)') 'Translational Coordinate :',mlabel(1)
        write(ilog,'(2a)') 'Dissociative Coordinate  :',mlabel(2)
        write(ilog,'(2a)') 'Vibrational Coordinate   :',mlabel(3)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(4)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(5)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(6)

       case (24)
        write(ilog,'(a)') 'NH3 Surface :'
        write(ilog,'(a)') 'in Jacobian Cartesian Coordinates '
        write(ilog,'(2a)') 'Translational Coordinate :',mlabel(1)
        write(ilog,'(2a)') 'Dissociative Coordinate  :',mlabel(2)
        write(ilog,'(2a)') 'Vibrational Coordinate   :',mlabel(3)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(4)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(5)
        write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(6)

      case (25)
        write(ilog,'(a)') 'NH3 Surface :'
        write(ilog,'(a)') 'in Cartesian Coordinates '
        write(ilog,'(a,3('' '',a))') 'N atom coordinates   :', &
                  trim(mlabel(1)),trim(mlabel(2)),trim(mlabel(3))
        write(ilog,'(a,3('' '',a))') 'H atom coordinates 1 :',&
                  trim(mlabel(4)),trim(mlabel(5)),trim(mlabel(6))
        write(ilog,'(a,3('' '',a))') 'H atom coordinates 2 :',&
                  trim(mlabel(7)),trim(mlabel(8)),trim(mlabel(9))
        write(ilog,'(a,3('' '',a))') 'H atom coordinates 3 :',&
                  trim(mlabel(10)),trim(mlabel(11)),trim(mlabel(12))

     case(26)
      write(ilog,'(a)') 'NH3 2D CASSCF Energy Surfaces'

     case(27)
      write(ilog,'(a)') 'NH3 2D CASSCF Dipole Surfaces'

     case(28)
      write(ilog,'(a)') 'NH3 2D CASSCF Polarisability Surfaces'

     case (29)
      if (hopipar(1,hl) .eq. 1) then
         write(ilog,'(a)') 'NOF S0 PES from E. Chapman '
      else if (hopipar(1,hl) .eq. 2) then
         write(ilog,'(a)') 'NOF S1 PES from E. Chapman '
      endif
      if (hopipar(2,hl) .eq. 1) then
      write(ilog,'(a)') 'Binding Coordinates '
         write(ilog,'(2a)') 'R1             :',mlabel(1)
         write(ilog,'(2a)') 'R2             :',mlabel(2)
         write(ilog,'(2a)') 'Bond Angle     :',mlabel(3)
      else if (hopipar(2,hl) .eq. 2) then
         write(ilog,'(a)') 'Jacobian Coordinates '
         write(ilog,'(2a)') 'Dissociative Coordinate  :',mlabel(1)
         write(ilog,'(2a)') 'Vibrational Coordinate   :',mlabel(2)
         write(ilog,'(2a)') 'Angular Coordinate       :',mlabel(3)
      endif

     case(30)
      write(ilog,'(2a,/a,/a)') &
       'Potential surface of HCN in Jacobian coordinates:', &
       ' R,r,cos(theta)', &
       'J.M. Bowman, B. Gazdy, J.A. Bentley, T.J. Lee, C.E. Dateo', &
       'J. Chem. Phys. 99 (1993) 308'

      write(ilog,'(2a)')'R (in a.u.) -- the h-(cn) distance : ',mlabel(1)
      write(ilog,'(2a)')'r (in a.u.) -- the cn distance : ',mlabel(2)
      write(ilog,'(2a)')'cos(theta) : ',mlabel(3)

     case(31)
      write(ilog,'(a,/a,/a)') &
      'Potential surface of HCN in Jacobian coordinates: R,r,theta', &
       'J.M. Bowman, B. Gazdy, J.A. Bentley, T.J. Lee, C.E. Dateo', &
       'J. Chem. Phys. 99 (1993) 308'

      write(ilog,'(2a)')'R (in a.u.) -- the h-(cn) distance : ',mlabel(1)
      write(ilog,'(2a)')'r (in a.u.) -- the cn distance : ',mlabel(2)
      write(ilog,'(2a)')'theta (in radians) : ',mlabel(3)

        case (32)
          write(ilog,'(2a,f9.4)') 'LSTH PES :', &
            '3D model in Cartesian Jacobian Coordinates.'
          write(ilog,'(2a)')'Atom 1 x Coordinate : ',mlabel(1)
          write(ilog,'(2a)')'Atom 1 z Coordinate : ',mlabel(2)
          write(ilog,'(2a)')'Atom 3 z Coordinate : ',mlabel(3)

       case (33)
        write(ilog,'(a)') 'NH3 Surface :'
        write(ilog,'(i1,a)') NH3map(7),'D Normal Mode Coordinates '
        do i=1,NH3map(7)
          if (NH3map(i)==1) then
            write(ilog,'(2a)') 'Inversion Coordinate          :',mlabel(i)
          elseif (NH3map(i)==2) then
            write(ilog,'(2a)') 'Bending Coordinate            :',mlabel(i)
          elseif (NH3map(i)==3) then
            write(ilog,'(2a)') 'Bending Coordinate            :',mlabel(i)
          elseif (NH3map(i)==4) then
            write(ilog,'(2a)') 'Symmetric Stretch Coordinate  :',mlabel(i)
          elseif (NH3map(i)==5) then
            write(ilog,'(2a)') 'Asymmetric Stretch Coordinate :',mlabel(i)
          elseif (NH3map(i)==6) then
            write(ilog,'(2a)') 'Asymmetric Stretch Coordinate :',mlabel(i)
          endif
        enddo

       case (34)
        write(ilog,'(a)') 'CH3FCl- Surface :'
        write(ilog,'(i2,a)') hopipar(1,hl),'D Normal Mode Coordinates '
        do i=1,hopipar(1,hl)
          select case (hopipar(i+1,hl))
            case (1)
              write(ilog,'(2a)') 'TS RC coordinate            :',mlabel(i)
            case (2)
              write(ilog,'(2a)') 'X-C-Y bend 1                :',mlabel(i)
            case (3)
              write(ilog,'(2a)') 'X-C-Y bend 2                :',mlabel(i)
            case (4)
              write(ilog,'(2a)') 'X-C-Y stretch coordinate    :',mlabel(i)
            case (5)
              write(ilog,'(2a)') 'CH3 rotation 1              :',mlabel(i)
            case (6)
              write(ilog,'(2a)') 'CH3 rotation 2              :',mlabel(i)
            case (7)
              write(ilog,'(2a)') 'CH3 inversion coordinate    :',mlabel(i)
            case (8)
              write(ilog,'(2a)') 'CH3 in-plane bend 1         :',mlabel(i)
            case (9)
              write(ilog,'(2a)') 'CH3 in-plane bend 2         :',mlabel(i)
            case (10)
              write(ilog,'(2a)') 'CH3 sym. stretch coordinate :',mlabel(i)
            case (11)
              write(ilog,'(2a)') 'CH3 asym. stretch 1         :',mlabel(i)
            case (12)
              write(ilog,'(2a)') 'CH3 asym. stretch 2         :',mlabel(i)
          end select
        enddo

       case (35)
        write(ilog,'(a)') 'H2CO Surface by Bowman 2004:'
        write(ilog,'(i2,a)') hopipar(1,hl),'D Normal Mode Coordinates '
        do i=1,hopipar(1,hl)
          select case (hopipar(i+1,hl))
            case (1)
              write(ilog,'(2a)') 'H-C-H out-of-plane wagging :',mlabel(i)
            case (2)
              write(ilog,'(2a)') 'H-C-H in-plane rocking     :',mlabel(i)
            case (3)
              write(ilog,'(2a)') 'H-C-H in-plane scissoring  :',mlabel(i)
            case (4)
              write(ilog,'(2a)') 'C-O stretch                :',mlabel(i)
            case (5)
              write(ilog,'(2a)') 'C-H asymmetric stretch      :',mlabel(i)
            case (6)
              write(ilog,'(2a)') 'C-H symmetric stretch     :',mlabel(i)
          end select
        enddo

      end select hopilab201

      elseif( hopfile(hl) == 202 ) then
         goto(501,501,502,502,502,502,502,502) hopilab(hl)

         routine='surfinfo (opfuncs/funcsrf.F), line 2466'
         write(message,'(a,i5)') 'Surface not found :',hopilab(hl)
         call errormsg

 501     return
 502     write(ilog,'(2a)') ' theta : ',mlabel(1)
         write(ilog,'(2a)') ' phi   : ',mlabel(2)

!hopilab202:select case (hopilab(hl))
!
!         case default
!          routine='surfinfo (opfuncs/funcsrf.F)'
!          write(message,'(a,i5)') 'Surface not found :',hopilab(hl)
!          call errormsg
!
!         case (501)
!          return
!
!         case(502)
!          write(ilog,'(2a)') ' theta : ',mlabel(1)
!          write(ilog,'(2a)') ' phi   : ',mlabel(2)
!
!         end select hopilab202


      elseif( hopfile(hl) .eq. 203 ) then
         call usersurfinfo(hopilab(hl),hopipar(1,hl),hoprpar(1,hl), &
                           mlabel,nn,ilog)

      elseif( hopfile(hl) .eq. 204 ) then
         call srflibinfo(hopilab(hl),hopipar(1,hl),hoprpar(1,hl), &
                           mlabel,nn,ilog)

      elseif( hopfile(hl) .eq. 205 ) then
         call qcheminfo(hl)

      elseif( hopfile(hl) .eq. 206 ) then
         call shepardinfo

      elseif( hopfile(hl) .eq. 207 ) then
         call qchem1info(hl)

      elseif( hopfile(hl) .eq. 208 ) then
         call qcheminfo(hl)

      elseif( hopfile(hl) .eq. 209 ) then
         call qcheminfo(hl)

      endif

      return

      end subroutine surfinfo


!##########################################################################
!      READ1DINFO
! This routine prints information on read1d to log file.
! It is called from htmsetup.
!##########################################################################

      subroutine read1dinfo (h,hl,done)

      implicit none

      integer(long)                          :: i,f,ilbl,jlbl,np
      integer(long),intent(in)               :: h,hl
      integer(long), dimension(maxhop), intent(inout) :: done
      integer(long), save                    :: nz = 0

      done(1) = 0

      do i=1,nz
         if(done(i).eq.hl) return
      enddo

      nz=nz+1
      done(nz)=hl
      f = htmdof(h)
      if(htmmuld(h).gt.0 .or. htmmode(h).gt.0) then
         routine='read1dinfo (opfuncs/funcsrf.F)'
         message='read1d for muld? Do you mean readsrf?'
         call warnmsg
      endif
      ilbl=index(hoplab(hl),' ')-1
      if(hopidef(hl).gt.0) then
         np=hopidef(hl)
         jlbl=index(deflabel(np),' ')-1
         write(ilog,'(/5a)') '-------  ',deflabel(np)(1:jlbl), &
                      '  =  ',hoplab(hl)(1:ilbl),'  ---------'
       else
          write(ilog,'(/3a)') &
          '-------  ',hoplab(hl)(1:ilbl),'  ---------'
       endif

       goto(36,37) hopilab(hl)-35
 36    write(ilog,'(a,i4,2x,a,i3,3x,a)') &
       'read1d, channel =',hopipar(1,hl),'ASCII    f =',f,modelabel(f)
       return

 37    write(ilog,'(a,i4,2x,a,i3,3x,a)') &
       'read1d, channel =',hopipar(1,hl),'binary   f =',f,modelabel(f)
       return
       end subroutine read1dinfo


      end module funcsrf


