c
c   Skokov et al., Chem. Phys. Lett. 312, 494 (1999)
c
       real*8 function v(x2,x1,x3t)
       implicit double precision  (a-h,o-z)
       parameter (hkcal=627.5096d0,hcm=219474.63d0)
       dimension xs(50),ys(50),zs(50)
       data ncall/1/
       common/svdcoef/cf(1116) 
       common/svdbas/ns1,ns2,ns3,nsvd,ind(3,1116)
       common/svdpar/xs,ys,zs,cx,cy,cz
c
c  x1- OH, x2- Cl-OH, x3t - cos(theta): Jacob. coord.
c  IMPORTANT: zero Jacobi angle corresponds to linear H-O-Cl 
c  return energy in hartrees, zero is HOCl minimum
c
       pi = 4.0d0 * atan (1.0)
C**load svd.data and hocl.2.pot
         if(ncall.eq.1) then
       call loadp
       open(77, file='svd.data',status='old')
       read(77,*)ns1,ns2,ns3,nsvd
       read(77,*)s1min,s1max     
       read(77,*)s2min,s2max     
       read(77,*)s3min,s3max
       do i=1,nsvd
       read(77,*)(ind(m,i),m=1,3),cf(i)
       enddo
       dx=(s1max-s1min)/(ns1-1)
       dy=(s2max-s2min)/(ns2-1)
       dz=(s3max-s3min)/(ns3-1)
       cx=-1.0d0/dx**2
       cy=-1.0d0/dy**2
       cz=-1.0d0/dz**2     
       do i=1,ns1
       xs(i)=s1min+(i-1)*dx
       enddo
       do i=1,ns2
       ys(i)=s2min+(i-1)*dy
       enddo
       do i=1,ns3
       zs(i)=s3min+(i-1)*dz
       enddo
       ncall=2
         endif
C**for ncall >= 1
       pi = 4.0d0 * atan (1.0)
       x3=dacos(x3t)*180.d0/pi
       call jzmat(x1,x2,x3,r1in,r2in,r3in) !transform coord.
       call adjust(r1in,r2in,r3in,r1,r2,r3)  !coord. scaling
       call calcpes(r1,r2,r3,e)             
       ehocl=-57430.07  !global minimum in cm-1
c**compute svd potential**
       dv=0.d0
       do is=1,nsvd
       i=ind(1,is)
       j=ind(2,is)
       k=ind(3,is)
       dv=dv+cf(is)*dexp(cx*(xs(i)-x2)**2)*
     &              dexp(cy*(ys(j)-x3)**2)*
     &              dexp(cz*(zs(k)-x1)**2)
       enddo
C**add up all pieces 
       v=(e-ehocl+dv)/hcm 
       return 
       end

       subroutine adjust(r1in,r2in,r3in,r1,r2,r3)
C**r1in, r2in,r3in - input values
C**r1,r2,r3 - corrected values
       implicit double precision  (a-h,o-z)
       parameter (hkcal=627.5096d0)
       pi = 4.0 * atan (1.0)
C**ADJUSTMENTS*************************************
C**Shift ground state to experimental one**********
      r1exp=1.822395 !exper. ground state
      r2exp=3.191692
      r3exp=4.0147

      r1e=1.81898 !minimum of current PES
      r2e=3.19687
      r3e=4.01328

      dr1=r1exp-r1e
      dr2=r2exp-r2e
      dr3=r3exp-r3e

      r1=r1in-dr1
      r2=r2in-dr2
      r3=r3in-dr3
C**Scaling of HOCl coordinates****************
      sc1=1.d0-0.0065*exp(-0.9*abs(r1-r1exp))
      sc2=1.d0+0.00162
      sc3=1.d0-0.0069*exp(-0.6*abs(r3-r3exp))

      r1=sc1*(r1-r1exp)+r1exp
      r2=sc2*(r2-r2exp)+r2exp
      r3=sc3*(r3-r3exp)+r3exp
C**check triangle rule
       if(r1.ge.r2+r3) r1=r2+r3-1.d-6
       if(r2.ge.r1+r3) r2=r1+r3-1.d-6
       if(r3.ge.r1+r2) r3=r1+r2-1.d-6
       return 
       end

      subroutine jzmat(x1,x2,x3,r1,r2,r3)
C****Transformation from jacobian (x1,x2,x3) for HOCL**
C****to z-matrix (r1,r2,r3) for HOCL                 **
C*   x1 - R(OH), x2 - R(HO-CL), x3 - jacobian angle   *
C*   r1 - R(OH), r2 - R(OCl), r3 - r(HCl)             *
C******************************************************
       implicit double precision (a-h,o-z)
       amh=1.007825035
       amo=15.99491463
       pi = 4.0 * atan (1.0)
       x3=x3*pi/180.0d0
       dx1=x1*amh/(amh+amo)
       r1=x1
       r2=sqrt(dx1**2+x2**2-2.d0*dx1*x2*cos(x3))
       ar=(dx1**2+r2**2-x2**2)/(2*dx1*r2)
       if (ar.gt.1.0d0) ar=1.0d0
       if (ar.lt.-1.0d0) ar=-1.0d0
       th=180.0d0*acos(ar)/pi
       r3=dsqrt(r1**2+r2**2-2.d0*r1*r2*ar)
       x3=x3*180.0d0/pi
       return
       end


       subroutine calcpes(r1,r2,r3,e)
       implicit real*8 (a-h,o-z)
c
c r1-OH, r2-OCl,r3-HCl, e in a.u. 
c
       y2a=func_ap(1,r1)        !OH 2-body term
       y2b=func_ap(2,r2)        !ClO 2-body term
       y2c=func_mlj(3,r3)       !HCl 2-body term
       y3=func_3b(4,r1,r2,r3)   !HOCl 3-body term
C**Change OCl on asymptote
       de=0.23931d0    !kcal/mol
       if(r2.le.5.001) y2b=y2b-de
       if(r2.gt.5.001) y2b=y2b-de*(1.d0-exp(-1.d0/(r2-5.0)**2)) 
       e=(y2a+y2b+y2c+y3)*349.75504d0
       return
       end


c*****************************
      subroutine loadp
      implicit double precision (a-h,o-z)
      parameter (maxnp=700)
      common/genpef/p9(5,maxnp),indi(maxnp),indj(maxnp),
     >              indk(maxnp),indl(maxnp),e1,mci,nlp
c
        open(unit=11,file='hocl.2.pot',status='OLD',form='FORMATTED')
        read(11,*)
        read(11,*) e1
        do i=1,3
         read(11,*)
         read(11,*) p9(i,1)
         read(11,*) (p9(i,j),j=2,p9(i,1)+1)
        enddo
        read(11,*)
        read(11,*) p9(4,1), nlp, mci
        read(11,*) (p9(4,i),i=2,nlp+1)
        do j=1,p9(4,1)-nlp
          read(11,*) indi(j),indj(j),indk(j),indl(j),
     >              p9(4,j+nlp+1)
        enddo
        return
        end
c*****************************
       function func_mlj(nr,r)
       implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(5,maxnp),indi(maxnp),indj(maxnp),
     >               indk(maxnp),indl(maxnp),e1,mci,nlp
c
        np=p9(nr,1)
        de=p9(nr,2)
        rnexp=p9(nr,3)
        ve=-de
        re=p9(nr,4)
        z=(r-re)/(r+re)
        beta=p9(nr,5)
        do j=1,np-4
         beta=beta+p9(nr,j+5)*(z**j)
        enddo
        s=1.0d0 - (re/r)**rnexp *exp(-beta*z)
        func_mlj=ve + de*s*s
        return
        end
c*****************************
       function func_ap(nr,r)
       implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(5,maxnp),indi(maxnp),indj(maxnp),
     >               indk(maxnp),indl(maxnp),e1,mci,nlp
c
        np=p9(nr,1)
        re=p9(nr,2)
        alpha=p9(nr,3)
        beta=p9(nr,4)
        dr=(r-re)
        rho=r*exp(-beta*dr)
        sum=0.d0
        do j=1,np-4
         sum=sum+p9(nr,j+4)*rho**(j)
        enddo
        func_ap=sum + abs(p9(nr,np+1))*exp(-alpha*dr)/r
        return
        end
c*****************************
       function func_3b(ns,r1,r2,r3)
       implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(5,maxnp),indi(maxnp),indj(maxnp),
     >               indk(maxnp),indl(maxnp),e1,mci,nlp
       dimension r(3)
       external cosd,acosd
c
       b4=0.60d0
       c1=0.15d0
       c2=0.21d0
       c3=0.11d0
c
        nplin=p9(ns,1)-nlp
        r(1)=r1
        r(2)=r2
        r(3)=r3
        r1e=p9(ns,2)
        r2e=p9(ns,3)
        r3e=p9(ns,4)
        b1=p9(ns,5)
        b2=p9(ns,6)
        b3=p9(ns,7)
c
       dn1=p9(ns,8)
       dn2=p9(ns,9)
       dn3=p9(ns,10)
       dcut1=p9(ns,11)
       dcut2=p9(ns,12)
       dcut3=p9(ns,13)
       g1=p9(ns,14)
       g2=p9(ns,15)
       g3=p9(ns,16)
       gg1=p9(ns,17)
       gg2=p9(ns,18)
       gg3=p9(ns,19)
c
        dr1=r1-r1e
        dr2=r2-r2e
        dr3=r3-r3e
c
        rho1=exp(-b1*dr1)
        rho2=exp(-b2*dr2)
        rho3=exp(-b3*dr3)
c
       bcos1=(r(3)*r(3)-r(1)*r(1)-r(2)*r(2))/(-2.d0*r(1)*r(2))
       bcos2=(r(2)*r(2)-r(1)*r(1)-r(3)*r(3))/(-2.d0*r(1)*r(3))
       bcos3=(r(1)*r(1)-r(2)*r(2)-r(3)*r(3))/(-2.d0*r(2)*r(3))
       sumcos=1.d0+bcos1+bcos2+bcos3
       theta1=acosd(bcos1)
       theta2=acosd(bcos3)
       w1=sumcos
       if((r(1).ne.0.d0).and.(r(2).ne.0.d0).and.(r(3).ne.0.d0)) then
         w2=1.d0/(1.d0/r(1) + 1.d0/r(2) + 1.d0/r(3))
       else
         w2=0.d0
       endif
c
       if((r(1).gt.dcut1).and.(r(2).gt.dcut2).and.(r(3).gt.dcut3)) then
         damp=exp(dn1/(dcut1-r(1)))*exp(dn2/(dcut2-r(2)))*
     >        exp(dn3/(dcut3-r(3)))
     >   *switch(((r(1)-gg1*r1e)/2.d0),g1)
     >   *switch(((r(2)-gg2*r2e)/2.d0),g2)
     >   *switch(((r(3)-gg3*r3e)/2.d0),g3)
       else
         damp=0.d0
       endif
c..
        sum=0.d0
        do m=1,mci
          i=indi(m)
          j=indj(m)
          k=indk(m)
          l=indl(m)
          if(l.eq.0) then
           sum=sum+p9(ns,m+nlp+1)*rho1**i *rho2**j *rho3**k 
          elseif(l.eq.1) then
           sum=sum+p9(ns,m+nlp+1)*rho1**i *rho2**j *rho3**k *w1
          else
           sum=sum+p9(ns,m+nlp+1)*rho1**i *rho2**j *rho3**k *w2
          endif
        enddo
c
       call citerm(r(1),r(2),b4,theta1,ci)
       call citerm0(r(3),r(2),r(1),b4,theta2,ci0)
c
       rho1=exp(-c1*dr1)
       rho2=exp(-c2*dr2)
       rho3=exp(-c3*dr3)
c
        do m=mci+1,nplin
          i=indi(m)
          j=indj(m)
          k=indk(m)
          l=indl(m)
          if(l.eq.1) then
           sum=sum+p9(ns,m+nlp+1)*rho1**i *rho2**j *rho3**k * ci
          else
           sum=sum+p9(ns,m+nlp+1)*rho1**i *rho2**j *rho3**k * ci0
          endif
        enddo
        func_3b=damp*sum
        return
        end
c****************************************
      function switch(x,a)
      implicit double precision (a-h,o-z)
      switch=0.5d0*(1.d0-tanh(a*x))
      return
      end
c****************************************
      subroutine citerm(r1,r2,b4,theta,ci)
      implicit real*8 (a-h,o-z)
      external cosd,sind
c
      bcos=cosd(theta)
c
       seam=3.247d0+0.9162d0*r1-0.37491d0*r1**2+0.02955d0*r1**3
c
       ridge= r2*bcos+seam
c
       tans=5.d0
c
       dampci=switch(3.0d0-r1,tans)
c
       ci= exp(-sqrt(ridge**2 + dampci/(r2+0.01d0) +
     >        (b4*r2*sind(theta))**2 + 0.0001d0))
     >      * switch((bcos-cosd(120.0d0)),tans)
     >      * switch((r1-3.0d0),tans/2.d0)
     >      * switch((1.2d0-r1),tans)
c
       return
       end
c************************
      subroutine citerm0(r3,r2,r1,b4,theta,ci)
      implicit real*8 (a-h,o-z)
      external cosd,sind
c
      bcos=cosd(theta)
      tans=5.d0
c
       dampci=switch(4.0d0-r3,tans)+switch(5.0d0-r2,tans)
c
        seam=-6.12d0 + 3.50d0*r3
c
       ridge= r2*bcos+seam
c
       ci= exp(-sqrt(ridge**2 + dampci +
     >        (b4*r2*sind(theta))**2 + 0.0001d0))
     >      * switch((bcos-cosd(170.0d0)),tans)
c
       return
       end
c-----------------------------------------------------------------------------
        double precision function cosd(x)

        implicit real*8 (a-h,o-z)
        data pi/3.141592653589d0/

        cosd = dcos(x*pi/180.0d0)

        return
        end
c-----------------------------------------------------------------------------
        double precision function sind(x)

        implicit real*8 (a-h,o-z)
        data pi/3.141592653589d0/

        sind = sin(x*pi/180.0d0)

        return
        end
c-----------------------------------------------------------------------------
        double precision function acosd(x)

        implicit real*8 (a-h,o-z)
        data pi/3.141592653589d0/

        acosd = acos(x)*180.d0/pi

        return
        end