!#######################################################################
! 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=204
!
! 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
!
!#######################################################################

      subroutine defsrflib(hopilab,label,ifile,pesopts,zpesopts,hopipar, &
                 hoprpar,maxhpar,pesmin,pesmax,lcutpes,ndof,ierr,ilog, &
                 opdatadef,mpirank)

      implicit none

! Quantics sizes
      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4) 
      
      integer(long)                                     :: nopt,ilbl,jlbl,&
                                                           np,icoord,ll,imod,&
                                                           interact,matdim,&
                                                           nofpars,dipcomponent,i1
      integer(long), intent(in)                         :: maxhpar,ndof,ilog
      integer(long), intent(in)                         :: mpirank
      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,bmkpr0
      real(dop), intent(out)                            :: pesmin,pesmax
      real(dop), dimension(*), intent(out)              :: hoprpar
      character(len=*), intent(in)                      :: label
      character(len=*), dimension(*), intent(in)        :: pesopts
      character(len=10)                                 :: matrow,matcol,raregas
      character(len=80)                                 :: message
      character(len=*), intent(in)                      :: opdatadef
      logical(kind=4)                                   :: ldoferr,rlog,zlog
      logical(kind=4), intent(out)                      :: lcutpes

      real(dop), parameter :: ev = 27.21138386d0

      logical(kind=4), save                             :: lfirst = .true.
      
!-----------------------------------------------------------------------
! Set defaults
!-----------------------------------------------------------------------
      hopilab=0
      ifile=0  
      ierr=0
      ldoferr=.false.

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

!-----------------------------------------------------------------------
! read and set standard options
!-----------------------------------------------------------------------

      icoord = 0
      tfac  = -999.0_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)
           
   10    continue
           
! specify coordinates
         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. 'tfac') then
            nopt=nopt+2
            ilbl=zpesopts(nopt-1,np)+1
            jlbl=zpesopts(nopt,np)
            read(pesopts(np)(ilbl:jlbl),*) tfac
         endif
         nopt=nopt+1
         ilbl=jlbl+1
         jlbl=zpesopts(nopt,np)
         if (jlbl .gt. 0) go to 10
      endif


!-----------------------------------------------------------------------
! Ar-H2O Surface from Richard Doyle(July, 2007) 
!-----------------------------------------------------------------------
      if (label .eq. 'arh2o') then

!  -  Set default options
            
         icoord = 2
!         tfac  = 0.5_dop
         lfirst=.false.
         
!  -  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)
           
   90       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 90
          endif


         if (icoord .eq. 2) then
            hopilab=1  !  Ar-H2O in Jacobian Coordinates
         endif

!-----------------------------------------------------------------------
! Schwenke Surface for H2O from Partridge and Schwenke.
!       (JCP 106, 15 (1996))
!
! icoord : 1 binding coordinates
!        : 2 Jacobian coordinates - first variant:  (OH)H
!        : 3 Jacobian coordinates - second variant: (HH)O
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'SchwenkeH2O') then


!  -  Set default options
          icoord = 2
          tfac   = 1.0_dop/16.9994_dop
          lfirst = .false.
!
! 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)

   91       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian1') then
               icoord=2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian2') then
               icoord=3
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if(jlbl .gt. 0) go to 91
         endif


! .........set default options
         if (icoord.eq.0) icoord=1
         if (tfac.eq.0.0_dop) then
            if (icoord.eq.2) tfac=1.0_dop/16.9994_dop
            if (icoord.eq.3) tfac=0.5_dop
         endif
!
! now assign surface
!
         if (icoord .eq. 1) then
            hopilab=2  !  SCHWENKE PES in Binding Coordinates
         elseif (icoord .eq. 2) then
            hopilab=3  !  SCHWENKE PES in Jacobian Coordinates: H-(OH)
         elseif (icoord .eq. 3) then
            hopilab=4  !  SCHWENKE PES in Jacobian Coordinates: O-(HH)
         endif

         lcutpes=.true.
!         if (pesmax .gt.  99999.0_dop) pesmax=0.0_dop
         hoprpar(1)=tfac

!-----------------------------------------------------------------------
! H3O Surface from ftp://rsc.anu.edu.au/pub/collins/ 
! R P A Bettens, M A Collins, M J T Jordan, D H Zhang JCP 112, 10162(2000)
!
! Only PES3 data points have been used for this
! icoord : 1 - Jacobian coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h3oyzcl') then

! - Set default options
           icoord = 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)
   92       continue
           
            if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jcartcoord') then
               icoord=2
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if(jlbl .gt. 0) go to 92
         endif

! now assign surface
!
         if (icoord .eq. 1) then
            hopilab=5  !  YZCL H3O PES in Jacobi Coordinates for H + H2O
         else if (icoord .eq. 2) then
            hopilab=6  !  YZCL H3O PES in Jacobi Cartesian Coordinates with COM at O
         endif
      call h3oyzclinput
!
! -- Read data for the surfaces --

!-----------------------------------------------------------------------
! Gly2H^+-H2O Emperical Surface, Clary et al JCP 116(v), 9829(2002) 
!-----------------------------------------------------------------------
      else if (label .eq. 'gly2pot') then

!  -  Set default options
         icoord = 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)
           
   93       continue
           
            if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=1
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 93
         endif


         if (icoord .eq. 1) then
            hopilab=7  !  Gly2H^+-H2O in Jacobian Coordinates
         endif

!-----------------------------------------------------------------------
! i2ar ground state surface
!-----------------------------------------------------------------------
      else if (label .eq. 'i2ar_0') then
         hopilab=8


!-----------------------------------------------------------------------
! Ozone surfaces from Woywood, Domcke and Schinke
!-----------------------------------------------------------------------
      else if (label(1:5) .eq. 'o3woy') then

!  -  Set default options

         icoord = 2
         tfac  = 0.5_dop
         lfirst=.false.

!  -  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)

  95       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 95
         endif

!
! get which function
!
         call defo3woy(hopipar(1),label(1:ll))

         hopilab=9      ! Ozone
         hopipar(2)=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

!-----------------------------------------------------------------------
! NH3 Surface from POTLIB site by Truhlar Group 
! reduced to 2D
! Z H Li, R Valero and D G Truhlar, Theor CHem Acc 118, 9-24(2007)
!-----------------------------------------------------------------------
      else if (label(1:7) .eq. 'nh3_tru') then

!  -  Set default options

         icoord = 1
         lfirst=.false.

!  -  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)

  96       continue

            icoord=1

            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 96
         endif


!
! get which function
!
         call defnh3_truhlar(hopipar(1),label(1:ll))

         hopilab=10     ! Truhlar NH3 (2D)
         hopipar(2)=icoord
! set default PES energy bounds
         lcutpes=.true.
         if (pesmax .gt.  99999.0_dop) pesmax=12.0_dop/ev

!-----------------------------------------------------------------------
! IBr Dipole surfaces
!-----------------------------------------------------------------------
      else if (label(1:8) .eq. 'ibr_css_') then
!
! get which function
!
         call defibr_css_dz(hopipar(1),label(1:ll))
         hopilab=11     ! IBr Dipole surfaces

!-----------------------------------------------------------------------
! IBr ENBO surfaces
!-----------------------------------------------------------------------
      else if (label(1:7) .eq. 'ibr_css') then
!
! get which function
!
         call defibr_css(hopipar(1),label(1:ll))
         hopilab=12     ! IBr ENBO surfaces

!-----------------------------------------------------------------------
! Sorbie & Murrell ground state H2O, Mol. Phys. 29 (1975) 1387 & 31 (1976) 905
! Options: Jacobi [default, (OH)-H, first coord is O-H, third is O-X-H]
!          Cartesian, 9D
!-----------------------------------------------------------------------
      elseif (label(1:ll)=='SorbieH2O') then
!  -  Set default options
        hopilab=13     ! Sorbie & Murrell H2O
        tfac=1.0079_dop/(15.999_dop+1.0079_dop)
        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)
   11     continue
          if (pesopts(np)(ilbl:jlbl)=='Jacobi') then
            hopipar(1) = 1
          elseif ((pesopts(np)(ilbl:jlbl)=='Cartesian').or.  &
                  (pesopts(np)(ilbl:jlbl)=='9D')) then
            hopipar(1) = 2
          endif
          nopt=nopt+1
          ilbl=jlbl+1
          jlbl=zpesopts(nopt,np)
          if (jlbl .gt. 0) goto 11
        endif
        hoprpar(1)=tfac

!-----------------------------------------------------------------------
! H4 surface. Aguado et al. J.C.P. 101, 4004, (1994)
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h4srf') then
         hopilab=14

!-----------------------------------------------------------------------
! NOCl S1 surface from R. Schinke
! icoord : 1 curvilinear Jacobian coordinates (3D)
!        : 2 BF Cartesian coordinates of the Jacobi vectors (3D)
!        : 3 BF-constrained Cartesian coordinates of the nuclei (3D)
!        : 4 2D/1D-constrained Cartesian coordinates of the nuclei (5D)
!        : 5 SF free Cartesian coordinates of the nuclei (9D)
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'nocl1sch') then
!
! -- 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)
   48       continue
            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)
            if (jlbl .gt. 0) go to 48
         endif

         hopilab=15
         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(opdatadef)
!
! -- 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

!-----------------------------------------------------------------------
! i2ar excited state surfaces
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'i2ar_dim') then

!     extract surface parameters
!     determine an index attached to the label by the input routine.
!     via this index we access the proper surface parameters.

         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then

            read(label(ilbl:jlbl),*,iostat=ierr) np
            if (ierr .ne. 0) then
               message='Cannot read np in i2ar_dim{np}'
               go to 999
            endif
         else
            go to 990
         endif

!     here we expect the following 4 surface parameters:
!     (1) dimension of the dim matrix (integer)
!     (2) matrixelement to be picked up: row (string)
!     (3) matrixelement to be picked up: column (string)
!     (4) rare gas species (string)

!     in (2) and (3) you may append the string "_RKR" or "_rkr" to the state
!     label in order to select the experimental RKR curve instead of the
!     default calculated curve.

         if (zpesopts(1,np).gt.0) then
            ilbl=1
            jlbl=zpesopts(1,np)
            read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) matdim
            if (ierr .ne. 0) then
               message='Cannot read dimension parameter '// &
                  'for i2ar_dim surface'
               go to 999
            endif
            nofpars=nofpars+1
         else
            go to 990
         endif
         if (zpesopts(2,np).gt.0) then
            ilbl=zpesopts(1,np)+1
            jlbl=zpesopts(2,np)
            matrow=pesopts(np)(ilbl:jlbl)
            nofpars=nofpars+1
         else
            go to 990
         endif
         if (zpesopts(3,np).gt.0) then
            ilbl=zpesopts(2,np)+1
            jlbl=zpesopts(3,np)
            matcol=pesopts(np)(ilbl:jlbl)
            nofpars=nofpars+1
         else
            go to 990
         endif
         if (zpesopts(4,np).gt.0) then
            ilbl=zpesopts(3,np)+1
            jlbl=zpesopts(4,np)
            raregas=pesopts(np)(ilbl:jlbl)
            nofpars=nofpars+1
         else
            go to 990
         endif

!     initialize common blocks for the i2ar dim calculations
         call i2ar_dim_init(matdim,matrow,matcol,raregas,ierr,message)
         if (ierr.ne.0) then
            message = 'Error initialising I2Ar surface'
            goto 999
         endif

!     everything is correct so far: assign label
         hopilab=16

  990    if (nofpars.lt.4) then
            message='Less than 4 parameters for i2ar_dim surface'
            ierr=1
            goto 999
         endif

!-----------------------------------------------------------------------
! HOO surface. DMBE IV of Varandas (J.Phys.Chem. 94 8073 (1990))
! icoord : 1 Jacobian  coordinates
!        : 2 Cartesian coordinates
!        : 3 Binding   coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'hoosrf') then
!........set default options
         icoord=1
         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)
 12         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)
                  message='Cannot read argument to HOO tfac keyword'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'cartesian') then
               icoord=2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
               icoord=3
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 12
         endif
         hoprpar(1)=tfac
!
! now assign surface
!
         if (icoord .eq. 1) then
            hopilab=17    !  Jacobian Coordinates
         elseif (icoord .eq. 2) then
            hopilab=18    !  Cartesian coordinates
         elseif (icoord .eq. 3) then
            hopilab=19    !  Binding coordinates
         endif

!-----------------------------------------------------------------------
! BKMP2 surface
! imod: 1 collinear
!     : 2 3D
! icoord : 1 binding coordinates
!        : 2 Jacobian coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'bkmp2') 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)
   20       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)
                  message='Cannot read argument to BKMP2 tfac keyword'
                  go to 999
               endif
            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. 'colinear') then
               imod=1
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 20
         endif
!
! now assign surface
!
         if (imod .eq. 1 .and. icoord .eq. 2) then
            hopilab=20  !  collinear model in Jacobian Coordinates.
         elseif (imod .eq. 1 .and. icoord .eq. 1) then
            hopilab=21  !  collinear model in Binding Coordinates.
         elseif (imod .eq. 2 .and. icoord .eq. 2) then
            hopilab=22  !  3D model in Jacobian Coordinates.
         elseif (imod .eq. 2 .and. icoord .eq. 1) then
            hopilab=23  !  3D model in Binding Coordinates.
         endif

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

         hoprpar(1)=tfac

!-----------------------------------------------------------------------
! C2H 2A'' surface from Martial Boggio-Pasqua
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'c2hasec') then
         hopilab=24
         lcutpes=.true.
         if (pesmax .gt.  99999.0_dop) pesmax=5.0_dop/ev

!-----------------------------------------------------------------------
! C2H 1A' surface from Martial. Jacobi coordinates.
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'c2ha1') then
         lcutpes=.true.
         if (pesmax .gt.  99999.0_dop) pesmax=25.0_dop/ev

         if (ndof .eq. 2) then
            hopilab=25  ! Collinear model, 2D
         else
            hopilab=26  !  3D
         endif

!-----------------------------------------------------------------------
! WSLFH Surface for H3O from Wu,Schatz,Lendvay,Fang and Harding
!       (JCP, 113 3150 (2000)).
!
! icoord : 1 - Binding coordinates
!          2 - Jacobian coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'wslfh') then

! - 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)
   30       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
                  message='Cannot read argument to WSLFH 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
                  message='Cannot read argument to WSLFH 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
            else if (pesopts(np)(ilbl:jlbl) .eq. '5D') then
               icoord=3
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 30
         endif
! - Assign surface

         if (icoord.eq.1) then
            hopilab=28  !  WSLFH PES in Binding Coordinates
         else if (icoord.eq.2) then
            hopilab=29  !  WSLFH PES in Jacobian Coordinates
         else if (icoord.eq.3) then
            hopilab=27  !  WSLFH PES 5D in Jacobian Coordinates
         endif
         if (tfac1 .gt. 1.0_dop .or. tfac1 .le. 0.0_dop) then
            message='tfac1 in WSLFH 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
            message='tfac2 in WSLFH must be in interval (0,1) .'
            ierr=1
            go to 999
         endif

! -- Read data for the surfaces --
         call wslfhinput(opdatadef)

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

!-----------------------------------------------------------------------
! Tully surface.  
!-----------------------------------------------------------------------
        else if (label(1:ll) .eq. 'tully') then
           hopilab=30
           hopipar(1) = 1
           hopipar(2) = 1
           hopipar(3) = 4       ! default: all six DOF
           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)
 50           continue
              if (pesopts(np)(ilbl:jlbl) .eq. 'nxy') then
                 ilbl=zpesopts(nopt+1,np)+1
                 jlbl=zpesopts(nopt+2,np)
                 nopt=nopt+2
                 read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hopipar(1)
                 if (ierr .ne. 0) then
                    write(6,'(i3,a)') np, pesopts(np)(1:72)
                    write(6,'(3i5,a)') ilbl,jlbl,pesopts(np)(ilbl:jlbl)
                    message='Cannot read argument to tully nxy keyword'
                    go to 999
                 endif
              else if (pesopts(np)(ilbl:jlbl) .eq. 'nz') then
                 ilbl=zpesopts(nopt+1,np)+1
                 jlbl=zpesopts(nopt+2,np)
                 nopt=nopt+2
                 read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) hopipar(2)
                 if (ierr .ne. 0) then
                    write(6,'(i3,a)') np, pesopts(np)(1:72)
                    write(6,'(3i5,a)') ilbl,jlbl,pesopts(np)(ilbl:jlbl)
                    message='Cannot read argument to tully nz keyword'
                    go to 999
                 endif
              else if (pesopts(np)(ilbl:jlbl).eq.'2D') then
                 hopipar(3)=1   ! 2D system
              else if (pesopts(np)(ilbl:jlbl).eq.'4Da') then
                 hopipar(3)=2   ! 4D system with plane
              else if (pesopts(np)(ilbl:jlbl).eq.'4Db') then
                 hopipar(3)=3   ! 4D system with angles
              else if (pesopts(np)(ilbl:jlbl).eq.'6D') then
                 hopipar(3)=4   ! 6D system                 
              endif
              nopt=nopt+1
              ilbl=jlbl+1
              jlbl=zpesopts(nopt,np)
              if (jlbl.gt.0) go to 50
           endif

           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' Numerical details:                 '
           write(ilog,'(a)') ' -----------------------------------'
           write(ilog,'(a,i6)') ' Grid size xy-plane : ',hopipar(1)
           write(ilog,'(a,i6)') '            z-plane : ',hopipar(2)
           write(ilog,'(a)') '                                    '
           if (hopipar(3).eq.1) then
              write(ilog,'(a)') ' 2D system: using z,r '
           else if (hopipar(3).eq.2) then
              write(ilog,'(a)') ' 4D system: using x,y,z,r '
           else if (hopipar(3).eq.3) then
              write(ilog,'(a)') ' 4D system: using z,r,theta,phi '
           else if (hopipar(3).eq.4) then
              write(ilog,'(a)') ' 6D system: using x,y,z,r,theta,phi '
           endif

!-----------------------------------------------------------------------
! H2O (PJT2) Surface from Polyansky, Jensen and Tennyson
!  (JCP 105, 6490-6497 (1996); JCP 101, 7651 (1994))
! icoord : 1 binding coordinates
!        : 2 Jacobian coordinates - first variant:  (OH)H
!        : 3 Jacobian coordinates - second variant: (HH)O
!-----------------------------------------------------------------------
      else if (label(1:ll).eq.'pjt2') then

          icoord=0
          tfac=0.0_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)
   60       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)
                  message='Cannot read argument to PJT2 tfac keyword'
                  go to 999
               endif
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binding') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian1') then
               icoord=2
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian2') then
               icoord=3
            else if (pesopts(np)(ilbl:jlbl) .eq. 'binding2') then
               icoord=4
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if(jlbl .gt. 0) go to 60
         endif

! .........set default options
         if (icoord.eq.0) icoord=1
         if (tfac.eq.0.0_dop) then
            if (icoord.eq.2) tfac=1.0_dop/16.9994_dop
            if (icoord.eq.3) tfac=0.5_dop
         endif
!
! now assign surface
!
         if (icoord .eq. 1) then
            hopilab=31  !  PJT2 PES in Binding Coordinates
         elseif (icoord .eq. 2) then
            hopilab=32  !  PJT2 PES in Jacobian Coordinates: H-(OH)
         elseif (icoord .eq. 3) then
            hopilab=33  !  PJT2 PES in Jacobian Coordinates: O-(HH)
         elseif (icoord .eq. 4) then
            hopilab=34  !  PJT2 PES in Binding Coord. with u=cos(theta)
         endif

         lcutpes=.true.
!         if (pesmax .gt.  99999.0_dop) pesmax=0.0_dop
         hoprpar(1)=tfac

!-----------------------------------------------------------------------
! WDSE Surface for H3O from Walch, Dunning, Schatz and Elgersma.
!       (JCP 72, 1303 (1980), JCP 73, 21 (1980))
!
! icoord : 1 - Binding coordinates
!          2 - Jacobian coordinates
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'wdse') then

! - 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)
   70       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
                  message='Cannot read argument to WDSE 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
                  message='Cannot read argument to WDSE 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 70
         endif
! - Assign surface

         if (icoord.eq.1) then
            hopilab=35  !  WDSE PES in Binding Coordinates
         else if (icoord.eq.2) then
            hopilab=36  !  WDSE PES in Jacobian Coordinates
         endif
         if (tfac1 .gt. 1.0_dop .or. tfac1 .le. 0.0_dop) then
            message='tfac1 in WDSE 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
            message='tfac2 in WDSE must be in interval (0,1) .'
            ierr=1
            go to 999
         endif

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

!------------------------------------------------------------------------
!        Potential surface for H2+H2 BMKP
!        see Boothroyd et al. J. Chem. Phys. 116 (2002) 666
!------------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h4bmkp') then

!  -  Set default options

         icoord = 2
         tfac1  = 0.5_dop
         tfac2  = 0.5_dop
         interact = 0

!  -  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)

   80       continue

            if (pesopts(np)(ilbl:jlbl) .eq. 'tfac1') then
                ilbl=zpesopts(2,np)+1
                jlbl=zpesopts(3,np)
                read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) tfac1
                if (ierr .ne.  0) then
                    message='Cannot read argument to BMKB tfac1 keyword'
                    go to 999
                endif
            else if (pesopts(np)(ilbl:jlbl) .eq.  'tfac2') then
                 ilbl=zpesopts(4,np)+1
                 jlbl=zpesopts(5,np)
                 read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) tfac2
                 if (ierr .ne.  0) then
                    message='Cannot read argument to BMKB tfac2 keyword'
                    go to 999
                 endif
              jlbl=zpesopts(1,np)
            endif
            if (pesopts(np)(ilbl:jlbl) .eq. 'cartesian') then
               icoord=1
            else if (pesopts(np)(ilbl:jlbl) .eq. 'jacobian') then
               icoord=2
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'jacobian2') then
               icoord=3
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'interaction') then
               interact=1
            endif
            nopt=nopt+1
            ilbl=jlbl+1
            jlbl=zpesopts(nopt,np)
            if (jlbl .gt. 0) go to 80
         endif
!
! now assign surface
!
         if (icoord .eq. 1) then
            if (interact.eq.1) then
               message='option "interaction" not implemented for ' &
                  //'Cartesian coordinates'
               ierr=1
               goto 999
            endif
            hopilab=37 !  bmkp PES in Cartesian Coordinates
         else if (icoord .eq. 2) then
            hopilab=38 !  bmkp PES in Jacobian Coordinates
         else if (icoord .eq. 3) then
            hopilab=39 !  bmkp PES in reordered Jacobian Coordinates
         endif

         if (tfac1 .gt. 1.0_dop .or. tfac1 .le. 0.0_dop) then
            message='tfac1 in bmkp must be in interval (0,1) .'
            ierr=1
            goto 999
         endif
         if (tfac2 .gt.  1.0_dop .or.  tfac2 .le.  0.0_dop) then
            message='tfac2 in bmkp must be in interval (0,1) .'
            ierr=1
            goto 999
         endif

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

!-----------------------------------------------------------------------
!  HFCO surface from Kato, J.Chem.Phys. 107 (1997),6114-6122
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'hfco') then
         hopilab=40

!-----------------------------------------------------------------------
!  H2-H2 4D surface from Diep-Johnson, JCP 112 (2000) 4465.
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h4dj') then
         hopilab=41

!-----------------------------------------------------------------------
!  H2-H2 BMKP/rigid rotor (4D). Option r0 sets rotor length.
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'h4bmkprr') then
!        set defaults      
         bmkpr0=1.449
!        parse options
         ilbl=index(label,'{')+1
         jlbl=index(label,'}')-1
         if (ilbl .gt. 1) then
            read(label(ilbl:jlbl),*) np
            jlbl=zpesopts(1,np)
            if (pesopts(np)(ilbl:jlbl) .eq. 'r0') then
               ilbl=zpesopts(2,np)+1
               jlbl=zpesopts(3,np)
               read(pesopts(np)(ilbl:jlbl),*,iostat=ierr) bmkpr0
               if (ierr .ne.  0) then
                   message='Cannot read argument to h4bmkprr keyword r0'
                   go to 999
               endif
            endif
         endif
!        store data
         hoprpar(1)=bmkpr0
         hopipar(1)=0 ! not interaction-potential
         hopilab   =42

!---------------------------------------------------------------------------
! The gas-surface model of Frankcombe, Collins, Zhang, JCP 137 (2012) 144701
!---------------------------------------------------------------------------
      else if (label(1:ll)=='gsmodel') then
         hopilab = 43

!------------------------------------------------------------------------------
! The malonaldehyde model of Mizukami, Habershon, Tew, JCP, 141, (2014), 144310
!------------------------------------------------------------------------------
      else if (label(1:5)=='malon') then
         hopilab = 44

!----------------------------------------------------------------------------
! The 13-mode salicylaldimine model of Polyak, Allen, Worth, JCP, 143, 
! (2015), 084121
! Modes in order: v1 v5 v7 v9 v10 v11 v13 v16 v22 v23 v24 v32 v36
!----------------------------------------------------------------------------
      else if (label(1:7)=='salicyl') then
         hopilab = 45

!------------------------------------------------------------------------------
! 4th-order, 1D double well based on the malonaldehyde model of Mizukami, 
! Habershon, Tew, JCP, 141, (2014), 144310
!------------------------------------------------------------------------------
      else if (label(1:5)=='double') then
         hopilab = 46

!------------------------------------------------------------------------------
! 2nd-order, 1D well for testing
!------------------------------------------------------------------------------
      else if (label(1:5)=='single') then
         hopilab = 47

!------------------------------------------------------------------------------
! ArAr surface
!-----------------------------------------------------------------------
      else if (label(1:4) .eq. 'arar') then
         hopilab=48 

!-----------------------------------------------------------------------
! H5O2+ Surface by Bowman
!-----------------------------------------------------------------------
      else if (label(1:6)  .eq. 'zundel') then
         hopilab = 49
         if (label(1:10)  .eq. 'zundel_dip')  hopilab = 50
         ! defaults
         hoprpar(1) = 1.d0  ! H_mass
         hoprpar(2) = 1.5d0 ! d
         zlog = .false.     ! z -> z/(R-2d)
         rlog = .false.     ! true: valence, false: jacobi
         dipcomponent = 0
         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)
 87         continue
            if (pesopts(np)(ilbl:jlbl) .eq. 'm') 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
                  message='Cannot read argument "m" of zundel'
                  go to 999
               endif
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'd0') 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
                  message='Cannot read argument "d0" of zundel'
                  go to 999
               endif
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'x') then
               if ( hopilab .eq. 45) then
                  message="Diple component must only be used with " &
                       //"the dipole routine, not the potential."
                  go to 999
               end if
               dipcomponent = 1
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'y') then
               if ( hopilab .eq. 45) then
                  message="Diple component must only be used with " &
                       //"the dipole routine, not the potential."
                  go to 999
               end if
               dipcomponent = 2
            elseif (pesopts(np)(ilbl:jlbl) .eq. "z") then
               if ( hopilab .eq. 45) then
                  message="Diple component must only be used with " &
                       //"the dipole routine, not the potential."
                  go to 999
               end if
               dipcomponent = 3

            elseif (pesopts(np)(ilbl:jlbl) .eq. 'jacobi') then
               rlog=.false.
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'valence') then
               rlog=.true.
            elseif (pesopts(np)(ilbl:jlbl) .eq. 'ztrafo') then
               zlog = .true.
            else
               message='Unknown parameter for zundel (zundel_dip) PES: ' &
                    //pesopts(np)(ilbl:jlbl)
               go to 999
            end if

            ilbl=jlbl+1
            nopt = nopt + 1
            jlbl=zpesopts(nopt,np)

            if (jlbl.gt.0) go to 87
            if ( hopilab .eq. 46) then
               if (dipcomponent.eq.0) then
                  message = 'Dipole component for zundel not set.'
                  go to 999
               end if
               call  h5o2_set_dip_component(dipcomponent,ilog)
            end if
            call h5o2_init(zlog,rlog,hoprpar(1),hoprpar(2),ilog,mpirank)

         end if

!-----------------------------------------------------------------------
!     Potential surface of HCN
!     J.M. Bowman, B. Gazdy, J.A. Bentley, T.J. Lee, C.E. Dateo
!     J. Chem. Phys. 99 (1993) 308
!-----------------------------------------------------------------------
      else if (label(1:ll) .eq. 'hcn') then
         hopilab = 51
         call inhcn(24,opdatadef)

!------------------------------------------------------------------------------
!     Formaldehyde PES of X. Wang, P.L. Houston, J.M. Bowman,
!     Phil. Trans. R. Soc. A, 375, 20160194, 
!     http://dx.doi.org/10.1098/rsta.2016.0194
!------------------------------------------------------------------------------
      else if (label(1:4)=='form') then
         hopilab = 52

!-----------------------------------------------------------------------
! !!! next free hopilab is 53 !!!
!-----------------------------------------------------------------------

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

!-----------------------------------------------------------------------
! if label is assigned, assign ifile now.
!-----------------------------------------------------------------------
      if (hopilab .gt. 0) ifile=204

      return

  999 continue
      write(ilog,'(a)') 'Error in Defsrflib (funcsrflib.F90)'
      write(ilog,'(a)') message
      ierr = 1
      return

      end subroutine defsrflib

!#######################################################################
! 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 ipes 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)
! in usersrf.F  
!
!#######################################################################

      subroutine vlibpoint (v,gpoint,ndof1,ifunc,hopipar,hoprpar,label,&
                            ilog,opdatadef)

      implicit none

      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4) 
      

      integer(long), intent(in)              :: ifunc,ndof1,ilog
      integer(long), dimension(*),intent(in) :: hopipar
      integer(long)                          :: i,j,n,icoord,ideriv
      real(dop)                              :: tfac,rd,rv,mh,mo,&
                                                q1,q2,q3,theta,gamma,phi,&
                                                r1,r2,r3

      real(dop)                               :: r,th1,th2,tfac1,tfac2,rr,&
                                                 RAB,RAC,RBC,RAD,RBD,RCD,x,y,z, &
                                                 Xrv,Yrv,Xrd,Yrd,&
                                                 cosdi,AL,BE,GA,Derr, &
                                                 Vfit,rch,rcf,rco,ohco,ofco,&
                                                 xO,yO,zO,xN,yN,zN,xCl,yCl,&
                                                 zCl,rdz,rvx,rvz,v0
!
      real(dop), dimension(3)                :: h1,h2,h3,o
      real(dop), dimension(3)                :: rcoo,drv,a,b,c,d,dvtot
      real(dop), dimension(6)                :: r4,dVdr
      real(dop), intent(out)                 :: v
      real(dop), dimension(*), intent(in)    :: hoprpar
      real(dop), dimension(ndof1), intent(in):: gpoint
      real(dop), dimension(4,3)              :: cc,dVdcc
      character(len=*), intent(in)           :: label
      character(len=80)                      :: message
      real(dop), parameter :: pi = 4.0_dop*atan(1.0_dop)

      character(len=*), intent(in) :: opdatadef

      external :: vh2bmkp,h4dj
      real(dop) :: vh2bmkp,h4dj

     real(kind=dop)::SorbieMurrellH2O

     select case (ifunc)
     
     case default
      write(message,'(a,i5)') 'Function not found :',ifunc
      write(ilog,'(/,a)') 'Error in Vlibpoint (funcsrflib.F90)'
      write(ilog,'(a)') message
      stop

     case (1)
      q1=gpoint(1)
      q2=gpoint(2)
      q3=gpoint(3)
      theta=gpoint(4)
      gamma=gpoint(5)
      phi=gpoint(6)

      call arh2o_potfn(q1,q2,q3,theta,gamma,phi,v)
      return

     case (2)
      r1=gpoint(1)
      r2=gpoint(2)
      theta=gpoint(3)

      call H2O_POTS(v,r1,r2,theta)
      return

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

      tfac=hoprpar(1)
      r1=rd**2+(tfac*rv)**2-2_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)

      call SchwenkeH2O(r1,r3,r2,v)
      return

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

      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)

      call SchwenkeH2O(r1,r2,r3,v)
      return

     case (5)
      r1    = gpoint(1)
      r2    = gpoint(2)
      r3    = gpoint(3)
      theta = gpoint(4)
      gamma = gpoint(5)
      phi   = gpoint(6)

      call pot_h3o(r1,r2,r3,theta,gamma,phi,v)
      v = v/27.2113845_dop
      return

     case (6)
      h1(1) = gpoint(1)
      h2(1) = gpoint(2)
      h2(2) = gpoint(3)
      h3(1) = gpoint(4)
      h3(2) = gpoint(5)
      h3(3) = gpoint(6)

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

      h1(2) = 0.0_dop
      h1(3) = 0.0_dop
      h2(3) = 0.0_dop

      mh = 1.007825_dop
      mo = 15.9994_dop
      h3(1) = h3(1) + (mh /(mh + mo))*h1(1) + (mh/(2.0_dop*mh + mo))*h2(1)
      h3(2) = h3(2) + (mh /(2.0_dop*mh + mo))*h2(2)

      h2(1) = h2(1) + (mh /(mh + mo))*h1(1)

      call yzcl_h3o(o,h1,h2,h3,v)
      v = v + 76.835881382_dop + 0.13162e-04_dop + 0.215736e-03_dop
      v = v/27.2113845_dop
      return

     case (7)
      r1    = gpoint(1)
      r2    = gpoint(2)
      r3    = gpoint(3)
      theta = gpoint(4)
      phi   = gpoint(5)
      gamma = gpoint(6)

      call gly2_wat_empr(r1,r2,r3,theta,phi,gamma,v)
      return
 
     case (8)
      r1=gpoint(1)
      theta=gpoint(2)
      call i2ar_0_2d(r1,theta,v)
      return

     case (9)! O3
      call func_o3woy(gpoint,v,hopipar(1),hopipar(2))
      return
 
     case (10)! 2D NH3
      call func_nh3_truhlar(gpoint,v,hopipar(1),hopipar(2))
      return

     case (11)! IBr
      call func_ibr_css(v,gpoint,hopipar(1),hopipar(2))
      return

     case (12)! IBr dipoles
      call func_ibr_css_dz(v,gpoint,hopipar(1),hopipar(2))
      return
 
     case (13)   ! Sorbie & Murrell H2O
       if (hopipar(1)==1) then
         rv=gpoint(1)
         rd=gpoint(2)
         theta=gpoint(3)
         tfac=hoprpar(1)
         r1=rv
         r2=sqrt(rd**2+(tfac*rv)**2-2*tfac*rd*rv*cos(theta))
         r3=sqrt(rd**2+((1-tfac)*rv)**2+2*(1-tfac)*rd*rv*cos(theta))
       else
         r1=sqrt(sum((gpoint(4:6)-gpoint(1:3))**2))
         r2=sqrt(sum((gpoint(7:9)-gpoint(1:3))**2))
         r3=sqrt(sum((gpoint(7:9)-gpoint(4:6))**2))
       endif
        ! (Units below are identical to those in source/lib/utilities/keyunits.f)
        ! Convert to Angstrom
       r1=r1*0.52917720859_8
       r2=r2*0.52917720859_8
       r3=r3*0.52917720859_8
       v=SorbieMurrellH2O(r1,r2,r3)
        ! Convert to hartree
       v=v/27.21138386_8
       return

     case (14) ! H4 Aguado et al
      R=gpoint(1)
      RAB=gpoint(2)
      RCD=gpoint(3)
      AL=gpoint(4)
      BE=gpoint(5)
      GA=gpoint(6)
      call h4srf(R,RAB,RCD,AL,BE,GA,v)
      return

     case (15) ! NOCl Schinke
      icoord=hopipar(1)
      if (icoord .eq. 1) then
         rd=gpoint(1)
         rv=gpoint(2)
         theta=gpoint(3)

      else if (icoord .eq. 2) then
!.....NOCl S1 PES from R. Schinke 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 S1 PES from R. Schinke 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 S1 PES from R. Schinke 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 S1 PES from R. Schinke 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)

      endif

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

     case (16) ! I2Ar excited state
      r=gpoint(1)
      rr=gpoint(2)
      theta=gpoint(3)
      call i2ar_dim(r,rr,theta,v)
      return

     case (17) ! HOO Varandas
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)
      tfac=hoprpar(1)
      r1=rd**2+(tfac*rv)**2-2.0_dop*tfac*rd*rv*cos(theta)
      r3=rd**2+((1.0_dop-tfac)*rv)**2+2.0_dop*(1.0_dop-tfac)*rd*rv*cos(theta)
      rcoo(1)=rv
      rcoo(2)=sqrt(r1)
      rcoo(3)=sqrt(r3)
      call HO2SUR(rcoo,v)
      v = v + 0.19157005_dop    ! NORMIERUNG AUF H + O2 ALS ASYMPTOTIK
      return

     case (18) ! HOO Varandas
      Xrd=gpoint(1)
      Yrd=gpoint(2)
      Xrv=gpoint(3)
      Yrv=gpoint(4)
      tfac=hoprpar(1)
!.....Transformation cartesian -> Jacobian
      rd = SQRT(Xrd**2 + Yrd**2)
      rv = SQRT(Xrv**2 + Yrv**2)
      theta = ACOS((Xrd*Xrv+Yrd*Yrv)/(rd*rv))
!.....Transformation Jacobian -> binding
      r1=rd**2+(tfac*rv)**2-2.0_dop*tfac*rd*rv*cos(theta)
      r3=rd**2+((1.0_dop-tfac)*rv)**2+2.0_dop*(1.0_dop-tfac)*rd*rv*cos(theta)
      rcoo(1)=rv
      rcoo(2)=sqrt(r1)
      rcoo(3)=sqrt(r3)
      call HO2SUR(rcoo,v)
      v = v + 0.19157005_dop    ! NORMIERUNG AUF H + O2 ALS ASYMPTOTIK
      return

     case (19) ! HOO Varandas
      rcoo(1)=gpoint(1)       ! O-O distance
      rcoo(2)=gpoint(2)       ! O-H distance
      rcoo(3)=gpoint(3)       ! O-O-H angle
      call HO2SUR(rcoo,v)
      v = v + 0.19157005_dop    ! NORMIERUNG AUF H + O2 ALS ASYMPTOTIK
      return

     case (20) ! BKMP2
      ideriv=0
      rd=gpoint(1)
      rv=gpoint(2)

      tfac=hoprpar(1)
      r1 = max(0.2_dop,rd-tfac*rv)
      r2 = max(0.2_dop,rv)
      r3 = r1+r2
      call bkmp2(r1,r2,r3,v,dvtot,ideriv)
      return

     case (21) ! BKMP2
      ideriv=0
      r1=gpoint(1)
      r2=gpoint(2)

      r1 = max(0.2_dop,r1)
      r2 = max(0.2_dop,r2)
      r3 = r1+r2
      call bkmp2(r1,r2,r3,v,dvtot,ideriv)
      return

     case (22) ! BKMP2
      ideriv=0
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      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.2_dop,r1)
      r2 = max(0.2_dop,r2)
      r3 = max(0.2_dop,r3)

      call bkmp2(r1,r2,r3,v,dvtot,ideriv)
! add energy to make H + H2 limit (H2 at equilib r=1.405au) zero of energy:
! this is 4.748280454 eV or 0.174496 hartree
      v=v+0.174496_dop

      return

     case (23) ! BKMP2
      ideriv=0
      r1=gpoint(1)
      r2=gpoint(2)
      phi=gpoint(3)

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

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

      call bkmp2(r1,r2,r3,v,dvtot,ideriv)
      return

     case (24) ! C2H
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      tfac=1.0_dop/13.0_dop
      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)
      rcoo(1)=r1
      rcoo(2)=r2
      rcoo(3)=r3

      call c2hasec(0,rcoo,v,drv,opdatadef)
      return

     case (25) ! C2H
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      tfac=1.0_dop/13.0_dop
      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)
      rcoo(1)=r1
      rcoo(2)=r2
      rcoo(3)=r3

      call c2ha1(0,rcoo,v,drv,opdatadef)
      return

     case (26) ! C2H
      rd=gpoint(1)
      rv=gpoint(2)
      theta=(20.0_dop/180.0_dop)*pi

      tfac=1.0_dop/13.0_dop
      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)
      rcoo(1)=r1
      rcoo(2)=r2
      rcoo(3)=r3

      call c2ha1(0,rcoo,v,drv,opdatadef)
      return

     case (27) ! WSLFH
      r1 =gpoint(1)
      r2 = 1.81008_dop
      r3 =gpoint(2)
      th1=gpoint(3)
      th2=gpoint(4)
      phi=gpoint(5)

      call coortransb(r1,r2,r3,th1,th2,phi,c,b,a,d)
      call wslfh(a,b,c,d,v)
      return

     case (28) ! WSLFH
      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,c,b,a,d)
      call wslfh(a,b,c,d,v)
      return

     case (29) ! WSLFH
      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,c,b,a,d)
      call wslfh(a,b,c,d,v)
      return

     case (30) ! TULLY
      continue
      if (hopipar(3).eq.4) then
         x = gpoint(1)
         y = gpoint(2)
         z = gpoint(3) + 4.8053134_dop
         r = gpoint(4) + 2.1259463_dop
         theta = gpoint(5)
         phi   = gpoint(6)
      else if (hopipar(3).eq.1) then
         x = 0.0_dop
         y = 0.0_dop
         z = gpoint(1) + 4.8053134_dop
         r = gpoint(2) + 2.1259463_dop
         theta = 0.0_dop
         phi = 0.0_dop
      else if (hopipar(3).eq.2) then
         x = gpoint(1)
         y = gpoint(2)
         z = gpoint(3) + 4.8053134_dop
         r = gpoint(4) + 2.1259463_dop
         theta = 0.0_dop
         phi = 0.0_dop
      else if (hopipar(3).eq.3) then
         x = 0.0_dop
         y = 0.0_dop
         z = gpoint(1) + 4.8053134_dop
         r = gpoint(2) + 2.1259463_dop
         theta = gpoint(3)
         phi   = gpoint(4)
      endif

      call tully(x,y,z,r,theta,phi,v, &
              hopipar(1),hopipar(2))
      return

     case (31) ! PJT2 H2O
      r1=gpoint(1)
      r2=gpoint(2)
      theta=gpoint(3)

      call pjt2(r1,r2,theta,v)
      return

     case (32) ! PJT2 H2O
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      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)

      call surfh2o(r1,r3,r2,v)
      return

     case (33) ! PJT2 H2O
      rd=gpoint(1)
      rv=gpoint(2)
      theta=gpoint(3)

      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)

      call surfh2o(r1,r2,r3,v)
      return

    case (34) ! PJT2 H2O
      r1=gpoint(1)
      r2=gpoint(2)
      theta=acos(gpoint(3))

      call pjt2(r1,r2,theta,v)
      return

     case (35) ! WDSE H3O
      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,a,b,d,c)
      call coortrans(a,b,c,d,rab,rac,rad,rbc,rbd,rcd,cosdi)
      call wdse(rab,rbc,rbd,rac,rad,rcd,v)
      return

     case (36) ! WDSE H3O
      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,a,b,c,d)
      call coortrans(a,b,c,d,rab,rac,rad,rbc,rbd,rcd,cosdi)
      call wdse(rab,rbc,rbd,rac,rad,rcd,v)
      return

     case (37)!.....bmkp PES in Cartesian Coordinates
      n=1
      do j=1,3
         do i=1,3
         cc(j,i)=gpoint(n)
         n=n+1
         enddo
         cc(4,j)=0.0_dop
      enddo
      ideriv=0
      call h4bmkp(cc,Vfit,r4,dVdcc,dVdr,Derr,ideriv )
      v=Vfit
      return

     case (38)!.....bmkp PES in Jacobian Coordinates
      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,c,b,d,a)
      do i=1,3
         cc(1,i)=c(i)
         cc(2,i)=b(i)
         cc(3,i)=d(i)
         cc(4,i)=a(i)
      enddo
      ideriv=0
      call h4bmkp( cc, Vfit, r4, dVdcc, dVdr, Derr, ideriv )

!     if option 'interaction' is set, subtract the contribution
!     of the two isolated H2; notice that
!       vh2bmkp(r) = bmkp(R=inf,r,r0,pi/2,pi/2,0)
!     where r0 = 1.4483_dop
      if (hopipar(1).ne.0) then
         v = vh2bmkp(r1,tfac1,tfac2)
         Vfit = Vfit - v
         v = vh2bmkp(r2,tfac2,tfac1)
         Vfit = Vfit - v
!     and add v = vh2bmkp(r0,tfac1,tfac2)
         Vfit = Vfit - 0.3482224535461817_dop
      endif

      v = Vfit
      return

     case (39)!.....bmkp PES in reordered Jacobian Coordinates
      R  =gpoint(1)
      r1 =gpoint(2)
      th1=gpoint(3)
      r2 =gpoint(4)
      th2=gpoint(5)
      phi=gpoint(6)
!     This is the same as ifunc=39 but with reordered coordinates.
      tfac1=hoprpar(1)
      tfac2=hoprpar(2)
      call coortransj(R,r1,r2,th1,th2,phi,tfac1,tfac2,c,b,d,a)
      do i=1,3
         cc(1,i)=c(i)
         cc(2,i)=b(i)
         cc(3,i)=d(i)
         cc(4,i)=a(i)
      enddo
      ideriv=0
      call h4bmkp( cc, Vfit, r4, dVdcc, dVdr, Derr, ideriv )

!     if option 'interaction' is set, subtract the contribution
!     of the two isolated H2; notice that
!       vh2bmkp(r) = bmkp(R=inf,r,r0,pi/2,pi/2,0)
!     where r0 = 1.4483_dop
      if (hopipar(1).ne.0) then
         v = vh2bmkp(r1,tfac1,tfac2)
         Vfit = Vfit - v
         v = vh2bmkp(r2,tfac2,tfac1)
         Vfit = Vfit - v
!     and add v = vh2bmkp(r0,tfac1,tfac2)
         Vfit = Vfit - 0.3482224535461817_dop
      endif

      v = Vfit
      return

     case (40)!.....hfco
      rch=gpoint(1)
      rcf=gpoint(2)
      rco=gpoint(3)
      ohco=gpoint(4)
      ofco=gpoint(5)
      phi=gpoint(6)
      call hfco(rch,rcf,rco,ohco,ofco,phi,v)
      return

     case (41)!.....DJ PES (H2-H2 4D)
      R  =gpoint(1)
      th1=gpoint(2)
      th2=gpoint(3)
      phi=gpoint(4)
      v=h4dj(R,th1,th2,phi)
      return

     case (42)!.....bmkp PES in Jacobian Coordinates / rigid rotors
      R  =gpoint(1)
      r1 =hoprpar(1)
      r2 =hoprpar(1)
      th1=gpoint(2)
      th2=gpoint(3)
      phi=gpoint(4)
!     This is the same as ifunc=39 but with fixed r1/r2.
      tfac1=0.5
      tfac2=0.5
      call coortransj(R,r1,r2,th1,th2,phi,tfac1,tfac2,c,b,d,a)
      do i=1,3
         cc(1,i)=c(i)
         cc(2,i)=b(i)
         cc(3,i)=d(i)
         cc(4,i)=a(i)
      enddo
      ideriv=0
      call h4bmkp( cc, Vfit, r4, dVdcc, dVdr, Derr, ideriv )

!     if option 'interaction' is set, subtract the contribution
!     of the two isolated H2; notice that
!       vh2bmkp(r) = bmkp(R=inf,r,r0,pi/2,pi/2,0)
!     where r0 = 1.4483_dop
      if (hopipar(1).ne.0) then
         v = vh2bmkp(r1,tfac1,tfac2)
         Vfit = Vfit - v
         v = vh2bmkp(r2,tfac2,tfac1)
         Vfit = Vfit - v
!     and add v = vh2bmkp(r0,tfac1,tfac2)
         Vfit = Vfit - 0.3482224535461817_dop
      endif

      v = Vfit
      return

     case (43)!.....gas-surface free diatomic model
       call gsmodel(gpoint,v)
       return

     case (44) ! Malonaldehyde potential of Mizukami, Habershon, Tew.
        call malonmodel(gpoint,ndof1,v)
        return

     case (45) ! 13 mode salicylaldimine potential of Polyak, Allen, Worth.
        call salicyl(gpoint,v)
        return

     case (46) ! 1 mode approximate malonaldehyde potential of Mizukamie et al.
        call doublewell(gpoint,v)
        return

     case (47) ! 1 mode harmonic potential.
        call singlewell(gpoint,v)
        return

     case (48)! ArAr
      call arpotv(v,gpoint)
      return

     case (49) ! zundel (h5o2+)
      call zundel(gpoint,v)
      return

     case (50) ! zundel (h5o2+) dipole
      call dipzundel(gpoint,v)
      return

     case (51) ! HCN
      v0=-93.2862240766843627
      call hcnpot(gpoint(1),gpoint(2),gpoint(3),v)
      v=v-v0
      return

     case (52) ! Formaldehyde PES of Wang, Houston, Bowman, 2017.
        call formpot(gpoint,v)
        return

     end select

      end subroutine vlibpoint

!##########################################################################
!      SURFINF
! 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 srflibinfo (hopilab,hopipar,hoprpar,mlabel,jj,ilog)

      implicit none

      integer, parameter :: dop = selected_real_kind(15,307)
      integer, parameter  :: long = selected_int_kind(9)  ! >= 9-digit integers
                                                          !  (close to integer*4) 
      

      integer(long), intent(in)                   :: hopilab,ilog
      integer(long), dimension(*), intent(in)     :: hopipar
      real(dop), dimension(*), intent(in)         :: hoprpar
      integer(long), intent(out)                  :: jj
      character(len=*), dimension(*),  intent(in) :: mlabel
      character(len=80)                           :: message
      integer(long),external                      :: slen
      integer(long)                               :: icoord,i

     select case (hopilab)

     case default
      write(message,'(a,i5)') 'Surface not found :',hopilab
      write(ilog,'(/,a)') 'Error in Srflibinfo (funcsrflib.F90)'
      write(ilog,'(a)') message
      stop

     case (1)
      write(ilog,'(a)') 'ArH2O 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)
      jj=6

     case (2)
      write(ilog,'(a)') 'SCHWENKE PES : in Binding Coordinates'
      write(ilog,'(2a)') 'O-H1 Bond     : ',mlabel(1)
      write(ilog,'(2a)') 'O-H2 Bond     : ',mlabel(2)
      write(ilog,'(2a)') 'H1-O-H2 Angle : ',mlabel(3)
      jj=3

     case (3)
      write(ilog,'(2a,f9.4)')'SCHWENKE PES: in Jacobian Coordinates', &
                           ', tfac =',hoprpar(1)
      write(ilog,'(2a)')'Dissociative Coordinate H-(OH) : ',mlabel(1)
      write(ilog,'(2a)')'Vibrational Coordinate  (O-H)  : ',mlabel(2)
      write(ilog,'(2a)')'Angular Coordinate             : ',mlabel(3)
      jj=3

     case (4)
      write(ilog,'(2a,f9.4)')'SCHWENKE PES: in Jacobi Coordinates', &
                           ', tfac =',hoprpar(1)
      write(ilog,'(2a)')'Dissociative Coordinate O-(HH) : ',mlabel(1)
      write(ilog,'(2a)')'Vibrational Coordinate  (H-H)  : ',mlabel(2)
      write(ilog,'(2a)')'Angular Coordinate             : ',mlabel(3)
      jj=3

     case (5)
      write(ilog,'(2a,f9.4)')'H3O PES: in Jacobi Coordinates', &
                        ', tfac =',hoprpar(1)
      write(ilog,'(2a)')'Dissociative Coordinate H-(OH) :',mlabel(1)
      write(ilog,'(2a)')'Vibrational Coordinate  (OH-H) :',mlabel(2)
      write(ilog,'(2a)')'Vibrational Coordinate  (O-H)  :',mlabel(3)
      write(ilog,'(2a)') 'Angular Coordinate            :',mlabel(4)
      write(ilog,'(2a)') 'Angular Coordinate            :',mlabel(5)
      write(ilog,'(2a)') 'Angular Coordinate            :',mlabel(6)
      jj=6

     case (6)
         write(ilog,'(2a,f9.4)')'H3O PES: in Jacobi Cartesian Coordinates'
         write(ilog,'(2a)')'Cartesian Coordinate of H1(x) (O-H1):',mlabel(1)
         write(ilog,'(2a)')'Cartesian Coordinate of H2(x) (H2-OH1):',mlabel(2)
         write(ilog,'(2a)')'Cartesian Coordinate of H2(y) (H2-OH1):',mlabel(3)
         write(ilog,'(2a)')'Cartesian Coordinate of H3(x) (H3-OH1H2):',mlabel(4)
         write(ilog,'(2a)')'Cartesian Coordinate of H3(y) (H3-OH1H2):',mlabel(5)
         write(ilog,'(2a)')'Cartesian Coordinate of H3(z) (H3-OH1H2):',mlabel(6)
         jj=6

     case (7)
      write(ilog,'(a)') 'Gly2Hplus-H2O 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)
      jj=6

     case (8)
      write(ilog,'(2a)') ' I2Ar ground state surface, ', &
                         '2D, I2 distance fixed at equilibrium.'
      write(ilog,'(2a)')'Dissociative Coordinate R : ',mlabel(1)
      write(ilog,'(2a)')'Angular Coordinate  theta : ',mlabel(2)
      jj=3

     case (9)
      write(ilog,'(a)') &
             'O3 surfaces from Woywood et al. JCP (97) 107: 7282'
      write(ilog,'(a)') &
             '                            and JCP (97) 107: 7296'
      if (hopipar(2) .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)
      else if (hopipar(2) .eq. 1) then
         write(ilog,'(a)') 'Binding Coordinates '
         write(ilog,'(2a)') 'Bond 1     : ',mlabel(1)
         write(ilog,'(2a)') 'Bond 2     : ',mlabel(2)
         write(ilog,'(2a)') 'Bond Angle : ',mlabel(3)
      endif
      if (hopipar(1) .eq. 1) then
         write(ilog,'(a)') &
               'O3 S1 diabatic surface from Woywood et al. '
      else if (hopipar(1) .eq. 2) then
         write(ilog,'(a)') &
               'O3 S2 diabatic surface from Woywood et al. '
      else if (hopipar(1) .eq. 3) then
         write(ilog,'(a)') &
               'O3 S12 diabatic coupling from Woywood et al. '
      else if (hopipar(1) .eq. 4) then
         write(ilog,'(a)') &
               'O3 S1 dipole surface from Woywood et al. '
      else if (hopipar(1) .eq. 5) then
         write(ilog,'(a)') &
               'O3 S2 dipole surface from Woywood et al. '
      else if (hopipar(1) .eq. 6) then
         write(ilog,'(a)') &
               'O3 S0 diabatic surface after Woywood et al. '
      endif
      jj=4

     case (10)
       write(ilog,'(a)') &
             '2D NH3 surfaces from Truhlar et al.'
      if (hopipar(2) .eq. 1) then
         write(ilog,'(a)') '2D Jacobian Coordinates '
         write(ilog,'(2a)') 'Dissociative Coordinate : ',mlabel(1)
         write(ilog,'(2a)') 'Angular Coordinate      : ',mlabel(2)
      endif
      if (hopipar(1) .eq. 1) then
         write(ilog,'(a)') &
               'NH3 S0 diabatic surface. '
      else if (hopipar(1) .eq. 2) then
         write(ilog,'(a)') &
               'NH3 S1 diabatic surface. '
      else if (hopipar(1) .eq. 3) then
         write(ilog,'(a)') &
               'NH3 S01 diabatic coupling. '
      else if (hopipar(1) .eq. 4) then
         write(ilog,'(a)') &
               'NH3 S0 adiabatic surface. '
      else if (hopipar(1) .eq. 5) then
         write(ilog,'(a)') &
               'NH3 S1 adiabatic surface. '
      endif
      jj=4

     case (11)
      write(ilog,'(a)') 'IBr Surfaces'
      jj = 3

     case (12)
      write(ilog,'(a)') 'IBr Dipole Surfaces'
      jj = 3

     case (13)
      write(ilog,'(a)',advance='no')'Sorbie & Murrell H2O PES:'
      if (hopipar(1)==1) then
        write(ilog,'(2a,f9.4)')' in Jacobian Coordinates', &
                           ', tfac =',hoprpar(1)
        write(ilog,'(2a)')'Vibrational Coordinate  (O-H)  : ',mlabel(1)
        write(ilog,'(2a)')'Dissociative Coordinate H-(OH) : ',mlabel(2)
        write(ilog,'(2a)')'Angular Coordinate             : ',mlabel(3)
        jj=3
      else
        write(ilog,'(a)')' in Cartesian coordinates'
        write(ilog,'(2a,2('' '',a))')'O coordinates   : ',mlabel(1:3)
        write(ilog,'(2a,2('' '',a))')'H1 coordinates  : ',mlabel(4:6)
        write(ilog,'(2a,2('' '',a))')'H2 coordinates  : ',mlabel(7:9)
        jj=9
      endif

     case (14)
      write(ilog,'(a)') &
      'H4 surface. Aguado et al. J.C.P. 101, 4004, (1994)'
      write(ilog,'(a)') 'Six Jacobi coordinates'
      write(ilog,'(2a)') 'R   : ',mlabel(1)
      write(ilog,'(2a)') 'RAB : ',mlabel(2)
      write(ilog,'(2a)') 'RCD : ',mlabel(3)
      write(ilog,'(2a)') 'AL  : ',mlabel(4)
      write(ilog,'(2a)') 'BE  : ',mlabel(5)
      write(ilog,'(2a)') 'GA  : ',mlabel(6)
      jj=7

     case (15)
      icoord=hopipar(1)
      write(ilog,'(a)') 'NOCl S1 PES from R. Schinke. '
      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)
         jj=4
      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)
         jj=4
      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)
         jj=4
      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)
         jj=6
      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)
         jj=10
      endif

     case (16)
      write(ilog,'(a)') '  I2Ar excited state surface, 3D'
      write(ilog,'(2a)')'I2 Vibrational Coord.   r : ',mlabel(1)
      write(ilog,'(2a)')'Dissociative Coordinate R : ',mlabel(2)
      write(ilog,'(2a)')'Angular Coordinate  theta : ',mlabel(3)
      jj=4

     case (17)
      write(ilog,'(a)') &
     'HOO surface. DMBE IV of Varandas (J.Phys.Chem. 94 8073 (1990))'
      write(ilog,'(a,f9.4)') 'Jacobian Coordinates. tfac =', &
                                 hoprpar(1)
      write(ilog,'(2a)') 'Dissociative Coordinate : ',mlabel(1)
      write(ilog,'(2a)') 'Vibrational Coordinate  : ',mlabel(2)
      write(ilog,'(2a)') 'Angular Coordinate      : ',mlabel(3)
      jj=4

     case (18)
      write(ilog,'(a)')&
     'HOO surface. DMBE IV of Varandas (J.Phys.Chem. 94 8073 (1990))'
      write(ilog,'(a,f9.4)') 'Cartesian coordinates. tfac =',&
                                 hoprpar(1)
      write(ilog,'(2a)')'X-Dissociative Coordinate: ',mlabel(1)
      write(ilog,'(2a)')'Y-Dissociative Coordinate: ',mlabel(2)
      write(ilog,'(2a)')'X-Vibrational Coordinate : ',mlabel(3)
      write(ilog,'(2a)')'Y-Vibrational Coordinate : ',mlabel(4)
      jj=5

     case (19)
      write(ilog,'(a)')&
      'HOO surface. DMBE IV of Varandas (J.Phys.Chem. 94 8073 (1990))'
      write(ilog,'(a)') 'Binding coordinates'
      write(ilog,'(2a)')'O-O Coordinate: ',mlabel(1)
      write(ilog,'(2a)')'O-H Coordinate: ',mlabel(2)
      write(ilog,'(2a)')'O-O-H angle   : ',mlabel(3)
      jj=4

     case (20)
      write(ilog,'(2a,f9.4)') 'BKMP2 PES : Collinear model ',&
               'in Jacobian Coordinates. tfac =',hoprpar(1)
      write(ilog,'(2a)')'Dissociative Coordinate : ',mlabel(1)
      write(ilog,'(2a)')'Vibrational Coordinate  : ',mlabel(2)
      jj=3

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

     case (22)
      write(ilog,'(a)')&
           'BKMP2 PES : 3D model in Jacobian Coordinates'
      write(ilog,'(2a)')'Dissociative Coordinate : ',mlabel(1)
      write(ilog,'(2a)')'Vibrational Coordinate  : ',mlabel(2)
      write(ilog,'(2a)')'Angular Coordinate      : ',mlabel(3)
      jj=4

     case (23)
      write(ilog,'(a)')&
           'BKMP2 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)
      jj=4

     case (24)
      write(ilog,'(a)') 'C2H Asec from Martial Boggio-Pasqua'
      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)
      jj=4

     case (25)
      write(ilog,'(a)') 'C2H A1 from Martial Boggio-Pasqua'
      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)
      jj=4

     case (26)
      write(ilog,'(a)') 'C2H A1 from Martial Boggio-Pasqua'
      write(ilog,'(a)') 'Jacobian Coordinates, Collinear model'
      write(ilog,'(2a)') 'Dissociative Coordinate : ',mlabel(1)
      write(ilog,'(2a)') 'Vibrational Coordinate  : ',mlabel(2)
      jj=3

     case (27)
      write(ilog,'(2a)') &
           'WSLFH PESi 5D in Jacobian Coordinates ', &
           '(R,r1,r2,th1,th2,phi)'
      write(ilog,'(2a)')'Scattering coordinate  (R)         : ', &
                           mlabel(1)
      write(ilog,'(2a)')'First Vibrational Coordinate (r1) ', &
                          'set to fixed value: 1.81008 au'
      write(ilog,'(2a)')'Second Vibrational Coordinate (r2) : ', &
                           mlabel(2)
      write(ilog,'(2a)')'First Angular Coordinate  (th1)    : ', &
                           mlabel(3)
      write(ilog,'(2a)')'Second Angular Coordinate  (th1)   : ', &
                           mlabel(4)
      write(ilog,'(2a)')'Torsional Angular Coordinate (phi) : ', &
                           mlabel(5)
      write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1), &
                                  '  tfac2 = ',hoprpar(2)
      jj=6

     case (28)
      write(ilog,'(a)') &
        'WSLFH 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)
      jj=7

     case (29)
      write(ilog,'(2a)') &
        'WSLFH 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  (th1)   : ', &
                           mlabel(5)
      write(ilog,'(2a)')'Torsional Angular Coordinate (phi) : ', &
                           mlabel(6)
      write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1), &
                                  '  tfac2 = ',hoprpar(2)
      jj=7

     case (30)
      write(ilog,'(a)') ' Tully surface, taken from:'
      write(ilog,'(a)') &
              ' Tully et al., J. Vac. Sci. Technol. A 11(4), '
      write(ilog,'(a)') ' Jul/Aug 1993, 1914-1920 '
      write(ilog,'(a)') ' Numerical details:                 '
      write(ilog,'(a)') ' -----------------------------------'
      write(ilog,'(a,i4)') ' Grid size xy-plane : ',hopipar(1)
      write(ilog,'(a,i4)') '            z-plane : ',hopipar(2)
      if (hopipar(3).eq.1) then
         write(ilog,'(a)') ' 2D system: using z,r '
         write(ilog,'(2a)') 'Dissociative Coordinate : ',mlabel(1)
         write(ilog,'(2a)') 'Vibrational Coordinate  : ',mlabel(2)
         jj=3
      else if (hopipar(3).eq.2) then
         write(ilog,'(a)') ' 4D system: using x,y,z,r '
         write(ilog,'(2a)') 'Surface Coordinate      x: ',mlabel(1)
         write(ilog,'(2a)') 'Surface Coordinate      y: ',mlabel(2)
         write(ilog,'(2a)') 'Dissociative Coordinate z: ',mlabel(3)
         write(ilog,'(2a)') 'Vibrational Coordinate  r: ',mlabel(4)
         jj=5
      else if (hopipar(3).eq.3) then
         write(ilog,'(a)') ' 4D system: using z,r,theta,phi '
         write(ilog,'(2a)') 'Dissociative Coordinate z: ',mlabel(1)
         write(ilog,'(2a)') 'Vibrational Coordinate  r: ',mlabel(2)
         write(ilog,'(2a)') 'Angular Coordinate  theta: ',mlabel(3)
         write(ilog,'(2a)') 'Angulat Coordinate    phi: ',mlabel(4)
         jj=5
      else if (hopipar(3).eq.4) then
         write(ilog,'(a)') ' 6D system: using x,y,z,r,theta,phi '
         write(ilog,'(2a)') 'Surface Coordinate      x: ',mlabel(1)
         write(ilog,'(2a)') 'Surface Coordinate      y: ',mlabel(2)
         write(ilog,'(2a)') 'Dissociative Coordinate z: ',mlabel(3)
         write(ilog,'(2a)') 'Vibrational Coordinate  r: ',mlabel(4)
         write(ilog,'(2a)') 'Angular Coordinate  theta: ',mlabel(5)
         write(ilog,'(2a)') 'Angulat Coordinate    phi: ',mlabel(6)
         jj=7
      endif

      case (31)
       write(ilog,'(a)') 'PJT2 PES : in Binding Coordinates'
       write(ilog,'(2a)') 'O-H1 Bond     : ',mlabel(1)
       write(ilog,'(2a)') 'O-H2 Bond     : ',mlabel(2)
       write(ilog,'(2a)') 'H1-O-H2 Angle : ',mlabel(3)
       jj=4
          
      case (32)
       write(ilog,'(2a,f9.4)') 'PJT2 PES : in Jacobian Coordinates', &
                            ', tfac =',hoprpar(1)  
       write(ilog,'(2a)')'Dissociative Coordinate H-(OH) : ',mlabel(1)
       write(ilog,'(2a)')'Vibrational Coordinate  (O-H)  : ',mlabel(2)
       write(ilog,'(2a)')'Angular Coordinate             : ',mlabel(3)
       jj=4

      case (33)
       write(ilog,'(2a,f9.4)') 'PJT2 PES : in Jacobian Coordinates', &
                            ', tfac =',hoprpar(1)
       write(ilog,'(2a)')'Dissociative Coordinate O-(HH) : ',mlabel(1)
       write(ilog,'(2a)')'Vibrational Coordinate  (H-H)  : ',mlabel(2)
       write(ilog,'(2a)')'Angular Coordinate             : ',mlabel(3)
       jj=4

      case (34)
       write(ilog,'(a)') 'PJT2 PES in Binding Coordinates cos(theta)'
       write(ilog,'(2a)') 'O-H1 Bond     : ',mlabel(1)
       write(ilog,'(2a)') 'O-H2 Bond     : ',mlabel(2)
       write(ilog,'(2a)') 'cos(H1-O-H2)  : ',mlabel(3)
       jj=4
          
     case (35)
      write(ilog,'(a)') &
           'WDSE PES in Binding Coordinates (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)
      jj=7

     case (36)
      write(ilog,'(a)') &
           'WDSE 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  (th1)   : ', &
                           mlabel(5)
      write(ilog,'(2a)')'Torsional Angular Coordinate (phi) : ', &
                           mlabel(6)
      write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1), &
                                  '  tfac2 = ',hoprpar(2)
      jj=7

     case (37)
      write(ilog,'(2a)') &
             'bmkp PES in Cartesian Coordinates ', &
             '(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4)'
      write(ilog,'(2a)')'Cartesian coordinate X1 : ',mlabel(1)
      write(ilog,'(2a)')'Cartesian coordinate Y1 : ',mlabel(2)
      write(ilog,'(2a)')'Cartesian coordinate Z1 : ',mlabel(3)
      write(ilog,'(2a)')'Cartesian coordinate X2 : ',mlabel(4)
      write(ilog,'(2a)')'Cartesian coordinate Y2 : ',mlabel(5)
      write(ilog,'(2a)')'Cartesian coordinate Z2 : ',mlabel(6)
      write(ilog,'(2a)')'Cartesian coordinate X3 : ',mlabel(7)
      write(ilog,'(2a)')'Cartesian coordinate Y3 : ',mlabel(8)
      write(ilog,'(2a)')'Cartesian coordinate Z3 : ',mlabel(9)
      write(ilog,'(a)')'Coordinate X4,Y4,Z3 are zero by definition.'
      jj=10

     case (38)
      write(ilog,'(a)') &
        'bmkp PES in Jacobian Coordinates (R,r1,r2,th1,th2,phi)'
      if (hopipar(1).ne.0) then
         write(ilog,'(a)') 'interaction part only'
      endif
      write(ilog,'(2a)')'Scattering Coordinate R = G1-G2: ',mlabel(1)
      write(ilog,'(2a)')'1. Vibrational Coord. r1 = H-H : ',mlabel(2)
      write(ilog,'(2a)')'2. Vibrational Coord. r2 = H-H : ',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 Coord.   phi : ',mlabel(6)
      write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1), &
                                  '  tfac2 = ',hoprpar(2)
      jj=7

     case (39)
      write(ilog,'(a)') &
        'bmkp PES in Jacobian Coordinates (R,r1,th1,r2,th2,phi)'
      if (hopipar(1).ne.0) then
         write(ilog,'(a)') 'interaction part only'
      endif
      write(ilog,'(2a)')'Scattering Coordinate R = G1-G2: ',mlabel(1)
      write(ilog,'(2a)')'1. Vibrational Coord. r1 = H-H : ',mlabel(2)
      write(ilog,'(2a)')'First  Angular Coordinate  th1 : ',mlabel(3)
      write(ilog,'(2a)')'2. Vibrational Coord. r2 = H-H : ',mlabel(4)
      write(ilog,'(2a)')'Second Angular Coordinate  th2 : ',mlabel(5)
      write(ilog,'(2a)')'Torsional Angular Coord.   phi : ',mlabel(6)
      write(ilog,'(a,f9.4,a,f9.4)') 'tfac1 = ',hoprpar(1), &
                                  '  tfac2 = ',hoprpar(2)
      jj=7

     case (40)
      write(ilog,'(a)')'HFCO surface (Kato) in valence coordinates'
      write(ilog,'(2a)') 'Distance CH (rch) : ',mlabel(1)
      write(ilog,'(2a)') 'Distance CF (rcf) : ',mlabel(2)
      write(ilog,'(2a)') 'Distance CO (rco) : ',mlabel(3)
      write(ilog,'(2a)') 'Angle HCO (ohco)  : ',mlabel(4)
      write(ilog,'(2a)') 'Angle FCO (ofco)  : ',mlabel(5)
      write(ilog,'(2a)') 'Dihedral Angle phi: ',mlabel(6)
      jj=7

     case (41)
      write(ilog,'(a)') 'DJ PES for H2-H2'
      write(ilog,'(2a)') 'center-of-mass distance : ',mlabel(1)
      write(ilog,'(2a)') '1st in-plane angle      : ',mlabel(2)
      write(ilog,'(2a)') '2nd in-plane angle      : ',mlabel(3)
      write(ilog,'(2a)') 'torsional angle         : ',mlabel(4)
      jj=5

     case (42)
      write(ilog,'(a)') 'BMKP PES for H2-H2 / rigid rotors'
      write(ilog,'(2a)') 'center-of-mass distance : ',mlabel(1)
      write(ilog,'(a,f8.6)') &
         'fixed rotor length      : ',hoprpar(1)
      write(ilog,'(2a)') '1st in-plane angle      : ',mlabel(2)
      write(ilog,'(2a)') '2nd in-plane angle      : ',mlabel(3)
      write(ilog,'(2a)') 'torsional angle         : ',mlabel(4)
      jj=5

     case(43)
      write(ilog,'(a,/a,/a)') &
      'The dissociative diatomic gas-surface potential of Frankcombe, Collins,',&
      ' Zhang, JCP 137 (2012) 144701',&
      'Cartesians in Bohr'

     case(44)
      write(ilog,'(a)') 'Malonaldehyde PES of Mizukami, Habershon, Tew,', &
       ' JCP, 141, (2014), 144310'

     case(45)
        write(ilog,'(a)') 'Salicylaldimine potential of Polyak, Allen, Worth, JCP, 143, (2015), 084121'
           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' 13 normal mode coordinates'
           write(ilog,'(a)') ' v1 v5 v7 v9 v10 v11 v13 v16 v22 v23 v24 v32 v36 '
           write(ilog,'(a)') '                                    '
           write(ilog,'(a)') ' -----------------------------------'


     case(46)
      write(ilog,'(a)') '4th-order, 1D PES based on malonaldehyde PES of Mizukami, Habershon, Tew, JCP, 141, (2014), 144310'

     case(47)
      write(ilog,'(a)') 'Harmonic test potential'

     case(48)
      write(ilog,'(a)') 'Argon dimer including highly repulsive region.'
      write(ilog,'(a)') &
            'K. Patkowski, G. Murdachaew, C-M. Fou and K.  Szalewicz'
      write(ilog,'(a)') &
            'Mol. Phys. 103 15-16 2005 2031-2045'

     case(49)
      write(ilog,'(2a,/a,/a)') &
                 '15D Potential surface of H5O2+ in polyspherical ', &
                 'coordinates. ', &
                 'X. Huang, B. J. Braams, and J. M. Bowman, ', &
                 'J. Chem. Phys. 122, 044308 (2005)'

     case(50)
      write(ilog,'(2a,/a,/a)') &
                 '15D dipole surface of H5O2+ in polyspherical ', &
                 'coordinates. ', &
                 'X. Huang, B. J. Braams, and J. M. Bowman, ', &
                 'J. Chem. Phys. 122, 044308 (2005)'

     case(51)
      write(ilog,'(a,/a,/a)') &
        'Potential surface of HCN',&
        'J.M. Bowman, B. Gazdy, J.A. Bentley, T.J. Lee, C.E. Dateo',&
        'J. Chem. Phys. 99 (1993) 308' 
      write(ilog,'(a)') 'Jacobi Coordinates in au '
      write(ilog,'(2a)') 'R (H--CN)              : ',mlabel(1)
      write(ilog,'(2a)') 'Gamma (in-plane angle) : ',mlabel(2)
      write(ilog,'(2a)') 'r (C--N)               : ',mlabel(3)

     case(52)
        write(ilog,'(a)') 'Formaldehyde PES of Wang, Houston, Bowman,', &
             ' Phil. Trans. R. Soc. A, 375, 20160194, http://dx.doi.org/10.1098/rsta.2016.0194'

     end select

      close(50)
      end subroutine srflibinfo


