c
c  Given r1(OH), r2(BrO), and r3(HBr), the potential energy and gradient 
c  is returned in a.u. (E relative to the HOBr minimum). 
c
c  This PES corresponds to that of Peterson, J. Chem. Phys. 113, 4598 (2000)
c  for the X^1A' state of HOBr based on icMRCI+Q/CBS(DZ-QZ) calculations, 3/00
c
c...Notes:
c          1)  the flag "nadj" in subroutine vhobr can be set to 0 or 1 to choose 
c              the ab initio (0) or adjusted (1) potentials
c          2)  the flag "ng" also in subroutine vhobr can be set to 0 or 1 to
c              turn on (1) and off (0) the calculation of analytical gradients
c
c   Note: gradients are currently experimental
c
c************************************************************************************
c.. sample driver program
       program hobr
       implicit real*8 (a-h,o-z)
       dimension grad(3)
       common/stuff/nadj,ngrad
       external cosd
c
c.. nadj=1 (use the adjusted potential)
       nadj=1
c.. ng=1 (compute gradients with energy)
       ngrad=0
c
       r1=1.8
       r2=3.5
       theta=103
       r3=sqrt(r1*r1+r2*r2-2.d0*r1*r2*cosd(theta))
c
       call vhobr(r1,r2,r3,energy,grad)
c
       print *, r1,r2,r3,theta
       print *, 'energy= ',energy
       if(ngrad.ne.0) then
        print *, 'grad(1)=',grad(1)
        print *, 'grad(2)=',grad(2)
        print *, 'grad(3)=',grad(3)
       endif
c
c  for the ab initio PES (nadj=0)
c
c  energy=  3.230567717570892E-04
c  grad(1)= -1.187233102898588E-02
c  grad(2)= 7.840759556519744E-03
c  grad(3)= 1.904006844479320E-03
c
c  for the adjusted PES (nadj=1)
c
c  energy=  3.568167597522186E-04
c  grad(1)= -1.313307777347717E-02
c  grad(2)= 8.267287657428473E-03
c  grad(3)= 1.374323323963971E-03
c
       end
c************************************************************************************
c
       subroutine vhobr(r1,r2,r3,energy,g)
       implicit double precision  (a-h,o-z)
       parameter (hkcal=627.5096d0)
       common/gradient/ng,rg(5),grad(3)
       common/stuff/nadj,ngrad
       dimension g(3)
       data ncall/1/
       save ncall
c
c.. nadj=1 (use the adjusted potential)
c      nadj=1
c.. ng=1 (compute gradients with energy)
c      ng=ngrad
c
       if(ncall.eq.1) then
        write(6,*) 'loading potential function'
        call loadp
        ncall=2
       endif
       
       if((r1.lt.0.94).or.(r2.lt.1.94).or.(r3.lt.1.59)) then
         energy=1000.d0/hkcal
         g(1)=0.d0
         g(2)=0.d0
         g(3)=0.d0
         return
       endif
c
       if(nadj.lt.1) then
        call calcpes(r1,r2,r3,e,nadj)
        ehobr=-160.1913d0
c
       else
        call adjust(r1,r2,r3,r1o,r2o,r3o)
        call calcpes(r1o,r2o,r3o,e,nadj)
        ehobr=-160.1917d0
       endif
c
       energy=(e-ehobr)/hkcal
c
       if(ng.eq.1) then
        if(nadj.eq.1) then
         g(1)=(grad(1)*rg(1)+grad(2)*rg(4)+grad(3)*rg(5))/hkcal
         g(2)=grad(2)*rg(2)/hkcal
         g(3)=grad(3)*rg(3)/hkcal
       else
         g(1)=grad(1)/hkcal
         g(2)=grad(2)/hkcal
         g(3)=grad(3)/hkcal
        endif
       endif
c
       return 
       end
c
       subroutine calcpes(r1,r2,r3,energy,nadj)
       implicit real*8 (a-h,o-z)
       common/gradient/ng,rg(5),grad(3)
c
       y2a=func_mljc(1,r1)
       y2b=func_mljc(2,r2)
       y2c=func_mljc(3,r3)
       y3=func_3b(4,r1,r2,r3)
c
       if(nadj.eq.1) then
C**Change OBr on asymptote
        de=2.04d0
        denom1=0.4780d0*exp(-4.370d0*(r2-5.5957d0))
        denom2=1.1522d0*exp(-0.8725d0*(r2-5.1402d0))
        denom=denom1+denom2+1.d0
        corr=-de/denom
        cut1=0.d0
        if(r3.gt.3.401d0) then
         cut1=exp(-0.1d0/(r3-3.4d0)**3)
        endif
        cut=1.d0-cut1
        y2b=y2b+corr-corr*cut
c
        if(ng.eq.1) then
         fac1=de/(denom**2)
         fac2=-4.37d0*denom1 - 0.8725d0*denom2
         grad(2)=grad(2)+(1.d0-cut)*fac1*fac2
         grad(3)=grad(3)+corr*cut1*0.3d0/((r3-3.4d0)**4)
        endif
       endif
c
       energy=y2a+y2b+y2c+y3
c
       return
       end
c
       subroutine adjust(r1in,r2in,r3in,r1,r2,r3)
       implicit real*8 (a-h,o-z)
       common/gradient/ng,rg(5),grad(3)
       data r1e,r2e,the,r1exp,r2exp,thexp,r3e,r3exp,sc1,sc2,sc3/
     > 1.819134d0,3.455719d0,102.7886d0  ,
     > 1.82226d0 , 3.45423d0 , 103.05d0  ,
     > 4.246681d0, 4.253816d0,
     > 0.99180d0,0.99235d0,0.99160d0 /
c
      cut=0.d0
      dcut=0.d0
      if(r1in.gt.3.001d0) then
       cut=exp(-0.1d0/(r1in-3.d0)**3)
      endif
C
C**ADJUSTMENTS*************************************
C**Shift ground state to experimental one**********
      dr1=r1exp-r1e
      dr2=r2exp-r2e
      dr3=r3exp-r3e
      ddr1=dr1-dr1*cut
      ddr2=dr2-dr2*cut
      ddr3=dr3-dr3*cut
c
      r1a=r1in-ddr1
      r2a=r2in-ddr2
      r3a=r3in-ddr3
c
      dsc1=1.d0-sc1
      dsc2=1.d0-sc2
      dsc3=1.d0-sc3
      ddsc1=dsc1-dsc1*cut
      ddsc2=dsc2-dsc2*cut
      ddsc3=dsc3-dsc3*cut
      sc1a=1.d0-ddsc1
      sc2a=1.d0-ddsc2
      sc3a=1.d0-ddsc3
c
C**Scaling of HOBr coordinates****************
      r1=sc1a*(r1a-r1exp)+r1exp
      r2=sc2a*(r2a-r2exp)+r2exp
      r3=sc3a*(r3a-r3exp)+r3exp
c
c.. contributions to the gradient computed here
      if(ng.eq.1) then
c
       rg(1)=sc1a
       rg(2)=sc2a
       rg(3)=sc3a
       rg(4)=0.d0
       rg(5)=0.d0
c
       if(r1in.gt.3.001d0) then
        dcut=cut*0.3d0/(r1in-3.d0)**4
        dsc1a=dsc1*dcut
        drp=1.d0+dr1*dcut
        dr2dr1=(sc2a*dr2+(r2a-r2exp)*dsc2)*dcut
        dr3dr1=(sc3a*dr3+(r3a-r3exp)*dsc3)*dcut
        rg(1)=sc1a*drp + (r1a-r1exp)*dsc1a
        rg(4)=dr2dr1
        rg(5)=dr3dr1
       endif
c
      endif
      return
      end
c       
      subroutine loadp
      implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(4,maxnp),indi(maxnp),indj(maxnp),
     >               indk(maxnp),indl(maxnp),e1,mci,nlp
c
        open(unit=11,file='hobr.pot',status='OLD',form='FORMATTED')
        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_mljc(nr,r)
       implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(4,maxnp),indi(maxnp),indj(maxnp),
     >              indk(maxnp),indl(maxnp),e1,mci,nlp
       common/gradient/ng,rg(5),grad(3)
       dimension coef(20)
c
        np=p9(nr,1)-5
        de=p9(nr,2)
        rn=p9(nr,3)
        rs=p9(nr,4)
        as=p9(nr,5)
        re=p9(nr,6)
        ve=-de
c
        z=(r-re)/(r+re)
c
        ebsw=exp(as*(r-rs))
        bsw=1.d0 + ebsw
c        
        binf=0.d0
        do j=1,np
         binf=binf+p9(nr,j+6)
        enddo
c
        coef(1)=1.d0
        bsum=p9(nr,7)
        do j=2,np
         coef(j)=coef(j-1)*z
         bsum=bsum+p9(nr,j+6)*coef(j)
        enddo
        bsum=bsum-binf
c
        beta=binf+bsum/bsw
        
        ebz=exp(-beta*z)
        rmsix=(re/r)**rn
        s=1.0d0 - rmsix *ebz
        
        func_mljc= ve + de*s*s
c
c..compute gradient
        if(ng.eq.1) then
c
        coef(2)=1.d0
        dbsum=p9(nr,8)
        do j=3,np
         coef(j)=coef(j-1)*z * float(j-1)/float(j-2)
         dbsum=dbsum+p9(nr,j+6)*coef(j)
        enddo
        dzdr=2*re/(r+re)**2
        dbsum=dbsum*dzdr
        dbeta=dbsum/bsw - bsum*as*ebsw/(bsw*bsw)
        
        dbz=beta*dzdr + z*dbeta
        ds=ebz*(6.d0*rmsix/r) + rmsix*ebz*dbz
        grad(nr)=2.d0*de*s*ds
c
        endif
        return
        end
c
       function func_3b(ns,r1,r2,r3)
       implicit double precision (a-h,o-z)
       parameter (maxnp=700)
       common/genpef/p9(4,maxnp),indi(maxnp),indj(maxnp),
     >              indk(maxnp),indl(maxnp),e1,mci,nlp
       common/switches/tans1,tans2
       common/gradient/ng,rg(5),grad(3)
       common/legendre/pl1(20),pl2(20)
       dimension xexa(0:20),xeya(0:20),xeza(0:20)
       dimension xexb(0:20),xeyb(0:20),xezb(0:20),pl(maxnp)
       dimension dxexa(0:20),dxeya(0:20),dxeza(0:20)
       dimension dxexb(0:20),dxeyb(0:20),dxezb(0:20)
       external cosd,acosd
c
        nplin=p9(ns,1)-nlp
        r1e=p9(ns,2)
        r2e=p9(ns,3)
        r3e=p9(ns,4)
        b1a=p9(ns,5)
        b2a=p9(ns,6)
        b2b=p9(ns,7)
        b3b=p9(ns,8)
        b3a=p9(ns,9)
        b1b=p9(ns,10)
        ncoord=p9(ns,11)        
        mxpol=p9(ns,12)
        alph=p9(ns,13)
c
       dn1=p9(ns,14)
       dn2=p9(ns,15)
       dn3=p9(ns,16)
       dcut1=p9(ns,17)
       dcut2=p9(ns,18)
       dcut3=p9(ns,19)
c
       b4a=p9(ns,20)
       b4b=p9(ns,21)
       tans1=p9(ns,22)
       tans2=p9(ns,23)
c
       dr1=r1-r1e
       dr2=r2-r2e
       dr3=r3-r3e
c
       expon1a=b1a*dr1
       expon2a=b2a*dr2
       expon2b=b2b*dr2
       expon3b=b3b*dr3
       rho1a=r1*dexp(-expon1a)
       rho2a=r2*dexp(-expon2a)
       rho2b=r2*dexp(-expon2b)
       rho3b=r3*dexp(-expon3b)
c
       bcos1=(r3*r3-r1*r1-r2*r2)/(-2.d0*r1*r2)
       bcos3=(r1*r1-r2*r2-r3*r3)/(-2.d0*r2*r3)
       if(bcos1.gt.1.d0) bcos1=1.d0
       if(bcos3.gt.1.d0) bcos3=1.d0
       if(bcos1.lt.-1.d0) bcos1=-1.d0
       if(bcos3.lt.-1.d0) bcos3=-1.d0

       theta1=acosd(bcos1)
       theta2=acosd(bcos3)
c
       if((r1.gt.dcut1).and.(r2.gt.dcut2).and.(r3.gt.dcut3)) then
         damp=exp(dn1/(dcut1-r1))*exp(dn2/(dcut2-r2))*
     >        exp(dn3/(dcut3-r3))
       else
         damp=0.d0
       endif
c
       cuta=switch(r1-r3,alph)
       cutb=switchp(r1-r3,alph)
       cut3a=switch((r3-2.0d0*r3e),b3a)
       cut1b=switch((r1-2.0d0*r1e),b1b)
c
       tha=cosd(180.d0-theta1)
       thb=cosd(180.d0-theta2)
       xexa(0)=1.d0
       xeya(0)=1.d0
       xeza(0)=1.d0
       xexa(1)=rho1a
       xeya(1)=rho2a
       xeza(1)=fleg1(tha,1)
       xexb(0)=1.d0
       xeyb(0)=1.d0
       xezb(0)=1.d0
       xexb(1)=rho3b
       xeyb(1)=rho2b
       xezb(1)=fleg2(thb,1)
       do m=2,mxpol
        xexa(m)=xexa(m-1)*rho1a
        xeya(m)=xeya(m-1)*rho2a
        xeza(m)=fleg1(tha,m)
        xexb(m)=xexb(m-1)*rho3b
        xeyb(m)=xeyb(m-1)*rho2b
        xezb(m)=fleg2(thb,m)
       enddo
c
      suma=0.d0
      sumb=0.d0
      do m=1,mci
       pl(m)=p9(ns,m+nlp+1)
       i=indi(m)
       j=indj(m)
       k=indk(m)
       if(indl(m).eq.1) then
         eta=xeza(k)*cut3a+1.d0
         suma=suma+xexa(i)*xeya(j)*eta*pl(m)
       else
         eta=xezb(k)*cut1b+1.d0
         sumb=sumb+xexb(i)*xeyb(j)*eta*pl(m)
       endif
      enddo

c conical intersection stuff
       call citerm(r1,r2,r3,b4a,bcos1,ci,dci1,dci2,dci3)
       call citerm0(r1,r2,r3,b4b,bcos3,ci0,dci01,dci02,dci03)
c      call citerm0(r3,r2,b4b,bcos3,ci0)
c
      sum1a=0.d0
      sum1b=0.d0
      do m=mci+1,nplin
       pl(m)=p9(ns,m+nlp+1)
       i=indi(m)
       j=indj(m)
       k=indk(m)
       if(indl(m).eq.1) then
        eta=xeza(k)*cut3a+1.d0
        sum1a=sum1a+xexa(i)*xeya(j)*eta*pl(m)
       else
        eta=xezb(k)*cut1b+1.d0
        sum1b=sum1b+xexb(i)*xeyb(j)*eta*pl(m)
       endif
      enddo
c
        sumv=(suma+sum1a*ci)*cuta + (sumb+sum1b*ci0)*cutb
        func_3b=sumv*damp
c
c..calculate gradients if requested
        if(ng.eq.1) then
c
         drho1a=rho1a*(1/r1 - b1a)
         drho2a=rho2a*(1/r2 - b2a)
         drho2b=rho2b*(1/r2 - b2b)
         drho3b=rho3b*(1/r3 - b3b)
c
         dcosa1=r2/(2.d0*r1*r1) - 1.d0/(2.d0*r2) - r3*r3/(2.d0*r1*r1*r2)
         dcosa2=r1/(2.d0*r2*r2) - 1.d0/(2.d0*r1) - r3*r3/(2.d0*r1*r2*r2)
         dcosa3=r3/(r1*r2)
         dcosb1=r1/(r2*r3)
         dcosb2=r3/(2.d0*r2*r2) - 1.d0/(2.d0*r3) - r1*r1/(2.d0*r2*r2*r3)
         dcosb3=r2/(2.d0*r3*r3) - 1.d0/(2.d0*r2) - r1*r1/(2.d0*r2*r3*r3)
     
         dgam1=dswitch((r3-2.d0*r3e),b3a)
         dgam2=dswitch((r1-2.d0*r1e),b1b)

         domeg1r1=dswitch((r1-r3),alph)
         domeg1r3=-domeg1r1
         domeg2r1=-domeg1r1
         domeg2r3=domeg1r1
c
         if(damp.ne.0.d0) then
          dchi1=dn1/((dcut1-r1)**2) * damp
          dchi2=dn2/((dcut2-r2)**2) * damp
          dchi3=dn3/((dcut3-r3)**2) * damp
         else 
          dchi1=0.d0
          dchi2=0.d0
          dchi3=0.d0
         endif

       dxexa(0)=0.d0
       dxeya(0)=0.d0
       dxeza(0)=0.d0
       dxexa(1)=1.d0
       dxeya(1)=1.d0
       dxeza(1)=1.d0
       dxexb(0)=0.d0
       dxeyb(0)=0.d0
       dxezb(0)=0.d0
       dxexb(1)=1.d0
       dxeyb(1)=1.d0
       dxezb(1)=1.d0
       do m=2,mxpol
        dxexa(m)=float(m)*xexa(m-1)
        dxeya(m)=float(m)*xeya(m-1)
        dxeza(m)=dfleg1(tha,m)
        dxexb(m)=float(m)*xexb(m-1)
        dxeyb(m)=float(m)*xeyb(m-1)
        dxezb(m)=dfleg2(thb,m)
       enddo
c
      dsuma1=0.d0
      dsuma2=0.d0
      dsuma3=0.d0
      dsumb1=0.d0
      dsumb2=0.d0
      dsumb3=0.d0
      do m=1,mci
       pl(m)=p9(ns,m+nlp+1)
       i=indi(m)
       j=indj(m)
       k=indk(m)
       if(indl(m).eq.1) then
         eta=xeza(k)*cut3a+1.d0
         deta=cut3a*dxeza(k)
         dsuma1=dsuma1+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa1) +
     >          eta*(xeya(j)*dxexa(i)*drho1a) )
         dsuma2=dsuma2+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa2) +
     >          eta*(xexa(i)*dxeya(j)*drho2a) )
         dsuma3=dsuma3+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa3 +
     >          xeza(k)*dgam1) )
       else
         eta=xezb(k)*cut1b+1.d0
         deta=cut1b*dxezb(k)
         dsumb1=dsumb1+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb1 +
     >          xezb(k)*dgam2) )
         dsumb2=dsumb2+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb2) +
     >          eta*(xexb(i)*dxeyb(j)*drho2b) )
         dsumb3=dsumb3+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb3) +
     >          eta*(xeyb(j)*dxexb(i)*drho3b) )
       endif
      enddo
c.. conical intersection stuff
      dsum1a1=0.d0
      dsum1a2=0.d0
      dsum1a3=0.d0
      dsum1b1=0.d0
      dsum1b2=0.d0
      dsum1b3=0.d0
      do m=mci+1,nplin
       pl(m)=p9(ns,m+nlp+1)
       i=indi(m)
       j=indj(m)
       k=indk(m)
       if(indl(m).eq.1) then
         eta=xeza(k)*cut3a+1.d0
         deta=cut3a*dxeza(k)
         dsum1a1=dsum1a1+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa1) +
     >          eta*(xeya(j)*dxexa(i)*drho1a) )
         dsum1a2=dsum1a2+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa2) +
     >          eta*(xexa(i)*dxeya(j)*drho2a) )
         dsum1a3=dsum1a3+pl(m)*( xexa(i)*xeya(j)*(deta*dcosa3 +
     >          xeza(k)*dgam1) )
       else
         eta=xezb(k)*cut1b+1.d0
         deta=cut1b*dxezb(k)
         dsum1b1=dsum1b1+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb1 +
     >          xezb(k)*dgam2) )
         dsum1b2=dsum1b2+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb2) +
     >          eta*(xexb(i)*dxeyb(j)*drho2b) )
         dsum1b3=dsum1b3+pl(m)*( xexb(i)*xeyb(j)*(deta*dcosb3) +
     >          eta*(xeyb(j)*dxexb(i)*drho3b) )
      endif
      enddo
c
        dsumv1=(dsuma1+dsum1a1*ci+sum1a*dci1)*cuta + 
     >         (dsumb1+dsum1b1*ci0+sum1b*dci01)*cutb +
     >         (suma+sum1a*ci)*domeg1r1 + (sumb+sum1b*ci0)*domeg2r1
c
        dsumv2=(dsuma2+dsum1a2*ci+sum1a*dci2)*cuta + 
     >         (dsumb2+dsum1b2*ci0+sum1b*dci02)*cutb
c
        dsumv3=(dsuma3+dsum1a3*ci+sum1a*dci3)*cuta + 
     >         (dsumb3+dsum1b3*ci0+sum1b*dci03)*cutb +
     >         (suma+sum1a*ci)*domeg1r3 + (sumb+sum1b*ci0)*domeg2r3

         grad(1)=grad(1) + damp*dsumv1 + sumv*dchi1
         grad(2)=grad(2) + damp*dsumv2 + sumv*dchi2
         grad(3)=grad(3) + damp*dsumv3 + sumv*dchi3
c
        endif
        return
        end
c
      function fleg1(x,n)
      implicit real*8 (a-h,o-z)
      common/legendre/pl1(20),pl2(20)
      nl=n+1 
      fleg1=0.d0
      if(n.eq.0) then
       fleg1=1.d0
       return
      elseif(n.eq.1) then
       fleg1=1.d0+x
       return
      else
       pl1(1)=1.d0
       pl1(2)=x
       twox=2.d0*x
       f2=x
       d=1.d0
       do j=3,nl
        f1=d
        f2=f2+twox
        d=d+1.d0
        pl1(j)=(f2*pl1(j-1)-f1*pl1(j-2))/d
       enddo
      endif
      do j=1,nl
       fleg1=fleg1+pl1(j)
      enddo
      return
      end
c
      function fleg2(x,n)
      implicit real*8 (a-h,o-z)
      common/legendre/pl1(20),pl2(20)
      nl=n+1 
      fleg2=0.d0
      if(n.eq.0) then
       fleg2=1.d0
       return
      elseif(n.eq.1) then
       fleg2=1.d0+x
       return
      else
       pl2(1)=1.d0
       pl2(2)=x
       twox=2.d0*x
       f2=x
       d=1.d0
       do j=3,nl
        f1=d
        f2=f2+twox
        d=d+1.d0
        pl2(j)=(f2*pl2(j-1)-f1*pl2(j-2))/d
       enddo
      endif
      do j=1,nl
       fleg2=fleg2+pl2(j)
      enddo
      return
      end
c
      function dfleg1(x,n)
      implicit real*8 (a-h,o-z)
      common/legendre/pl1(20),pl2(20)
      nl=n+1 
      dfleg1=1.d0
      if(n.eq.0) then
       dfleg1=0.d0
       return
      elseif(n.eq.1) then
       return
      else
      do j=3,nl
       nn=j-1
       dfleg1=dfleg1+float(nn)*(x*pl1(j)-pl1(j-1))/(x**2 -1)
      enddo
      endif
      return
      end
c
      function dfleg2(x,n)
      implicit real*8 (a-h,o-z)
      common/legendre/pl1(20),pl2(20)
      nl=n+1 
      dfleg2=1.d0
      if(n.eq.0) then
       dfleg2=0.d0
       return
      elseif(n.eq.1) then
       return
      else
      do j=3,nl
       nn=j-1
       dfleg2=dfleg2+float(nn)*(x*pl2(j)-pl2(j-1))/(x**2 -1)
      enddo
      endif
      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
      function switchp(x,a)
      implicit double precision (a-h,o-z)
      switchp=0.5d0*(1.d0+tanh(a*x))
      return
      end
c
      subroutine citerm(r1,r2,r3,b4,bcos,ci,dci1,dci2,dci3)
      implicit real*8 (a-h,o-z)
      common/switches/tans,dummy
      common/gradient/ng,rg(5),grad(3)
      data a1,b1,c1,a2,b2,c2/3.9519d0,0.075273d0,0.0384d0,
     >                       3.1367d0,0.660d0,0.1667d0/
      data tiny/1.d-14/
      external cosd
      two=2.d0
c
      seam1=a1+b1*r1-c1*r1*r1
      seam2=a2+b2*r1-c2*r1*r1
      sw1=switch(r1-1.9d0,25.d0)
      sw2=switch(1.9d0-r1,25.d0)
c
      seam=(seam1*sw1 + seam2*sw2 )
c
       ridge= r2*bcos+seam
c
       if(r1.gt.3.0d0) then
         dampci=0.05d0*(1.d0-exp(-(r1-3.0d0)))
       else
         dampci=0.d0
       endif
c
       sw3=switch((bcos-cosd(160.0d0)),tans)
       sw4=switch((r1-3.0d0),tans)
       sw5=switch((1.2d0-r1),tans)
       sw6=switch((r2-4.2d0),tans)
       sw7=switch((3.6d0-r2),tans)
       fac=b4*b4*r2*r2
       sin2=1.d0-bcos*bcos
       part1=sqrt(ridge**2 + dampci + fac*sin2 + 0.0001d0)
       part2=exp(-part1)
       fullsw=sw3 * sw4 * sw5 * sw6 * sw7
       ci= part2 * fullsw
c
       if(ng.eq.1) then
c
        dseam=(b2-two*c2*r1)*sw2 + (b1-two*c1*r1)*sw1 - 
     >    seam2*dswitch(1.9-r1,25.d0) + seam1*dswitch(r1-1.9d0,25.d0)
c
        dcos1=-(r2/(two*r1*r1)-1.d0/(two*r2)-r3*r3/(two*r1*r1*r2))
        dcos2=-(r1/(two*r2*r2)-1.d0/(two*r1)-r3*r3/(two*r1*r2*r2))
        dcos3=-r3/(r1*r2)
c
        dridge1=dseam+r2*dcos1
        dridge2=bcos+r2*dcos2
        dridge3=r2*dcos3
c
        if(dampci.ne.0.d0) then
         ddampci=0.05d0-dampci 
        else
         ddampci=0.d0
        endif
        tp1=two*part1
c
        if(abs(sw3).lt.tiny) sw3=sign(tiny,sw3)
        if(abs(sw4).lt.tiny) sw4=sign(tiny,sw4)
        if(abs(sw5).lt.tiny) sw5=sign(tiny,sw5)
        if(abs(sw6).lt.tiny) sw6=sign(tiny,sw6)
        if(abs(sw7).lt.tiny) sw7=sign(tiny,sw7)
c
        dpart2r1=-(part2*(ddampci-two*fac*bcos*dcos1 + 
     >             two*ridge*dridge1))/tp1
        dpart2r2=-(part2*(two*b4*b4*sin2*r2 -
     >            two*fac*bcos*dcos2 +
     >            two*ridge*dridge2))/tp1
        dpart2r3=-(part2*(-two*fac*bcos*dcos3 +
     >            two*ridge*dridge3))/tp1
c
       dswth=dswitch((bcos-cosd(160.0d0)),tans)
c
       dci1=part2*(fullsw/sw4*dswitch((r1-3.0d0),tans) -
     >      fullsw/sw5*dswitch((1.2d0-r1),tans) +
     >      fullsw/sw3*dswth*dcos1) +
     >      fullsw*dpart2r1
       dci2=part2*(fullsw/sw6*dswitch((r2-4.2d0),tans) -
     >      fullsw/sw7*dswitch((3.6d0-r2),tans) +
     >      fullsw/sw3*dswth*dcos2) +
     >      fullsw*dpart2r2
       dci3=part2*fullsw/sw3*dswth*dcos3 +
     >      fullsw*dpart2r3
c
       endif
       return
       end
c
      subroutine citerm0(r1,r2,r3,b4,bcos,ci,dci01,dci02,dci03)
      implicit real*8 (a-h,o-z)
      common/switches/tans1,tans
      common/gradient/ng,rg(5),grad(3)
      data a1,b1/-12.358d0,5.785d0/
      data tiny/1.d-14/
      external cosd
      two=2.d0
c
       if((r2.gt.5.0d0).and.(r3.gt.3.0d0)) then
         dampci=0.05d0*(1.d0-exp(-(r2-5.0d0)))*(1.d0-exp(-(r3-3.0d0)))
       else
         dampci=0.d0
       endif
c
        seam=a1 + b1*r3
c
       ridge= r2*bcos+seam
c
       sw1=switch((bcos-cosd(175.0d0)),tans)
       sw2=switch((2.4d0-r2),tans)
       sw3=switch((r2-5.0d0),tans)
       sw4=switch((2.4d0-r3),tans)
       sw5=switch((r3-3.0d0),tans)
       fac=b4*b4*r2*r2
       sin2=1.d0-bcos*bcos
       part1=sqrt(ridge**2 + dampci + fac*sin2 + 0.0001d0)
       part2=exp(-part1)
       fullsw=sw1 * sw2 * sw3 * sw4 * sw5
       ci= part2 * fullsw
c
       if(ng.eq.1) then
c
        dseam=b1
        dcos1=-r1/(r2*r3)
        dcos2=-(r3/(2.d0*r2*r2)-1.d0/(2.d0*r3)-r1*r1/(2.d0*r2*r2*r3))
        dcos3=-(r2/(2.d0*r3*r3)-1.d0/(2.d0*r2)-r1*r1/(2.d0*r2*r3*r3))
c
        dridge1=r2*dcos1
        dridge2=bcos+r2*dcos2
        dridge3=dseam+r2*dcos3
c
        if(dampci.ne.0.d0) then
         ddampci2=0.05d0*exp(-(r2-5.0d0))*(1.d0-exp(-(r3-3.0d0)))
         ddampci3=0.05d0*exp(-(r3-3.0d0))*(1.d0-exp(-(r2-5.0d0)))
        else
         ddampci2=0.d0
         ddampci3=0.d0
        endif
c
        tp1=two*part1
c
        dpart2r1=-(part2*(-two*fac*bcos*dcos1 +
     >             two*ridge*dridge1))/tp1
        dpart2r2=-(part2*(two*b4*b4*sin2*r2 +
     >            ddampci2-two*fac*bcos*dcos2 +
     >            two*ridge*dridge2))/tp1
        dpart2r3=-(part2*(ddampci3-two*fac*bcos*dcos3 +
     >            two*ridge*dridge3))/tp1
c
        dswth=dswitch((bcos-cosd(175.0d0)),tans)
c
        if(abs(sw1).lt.tiny) sw1=sign(tiny,sw1)
        if(abs(sw2).lt.tiny) sw2=sign(tiny,sw2)
        if(abs(sw3).lt.tiny) sw3=sign(tiny,sw3)
        if(abs(sw4).lt.tiny) sw4=sign(tiny,sw4)
        if(abs(sw5).lt.tiny) sw5=sign(tiny,sw5)
c
        dci01=part2*fullsw/sw1*dswth*dcos1 +
     >        fullsw*dpart2r1
        dci02=part2*(-fullsw/sw2*dswitch((2.4d0-r2),tans) +
     >        fullsw/sw3*dswitch((r2-5.0d0),tans) +
     >        fullsw/sw1*dswth*dcos2) +
     >        fullsw*dpart2r2
        dci03=part2*(-fullsw/sw4*dswitch((2.4d0-r3),tans) +
     >        fullsw/sw5*dswitch((r3-3.0d0),tans) +
     >        fullsw/sw1*dswth*dcos3) +
     >        fullsw*dpart2r3
c
       endif
       return
       end
c
      function dswitch(x,a)
      implicit double precision (a-h,o-z)
      dswitch=-(a/cosh(a*x)**2)/2.d0
      return
      end

       double precision function cosd(x)
       implicit double precision (a-h,o-z)
       pi=acos(-1.d0)
       theta=x/180.d0 * pi
       cosd=cos(theta) 
       return
       end

       double precision function acosd(x)
       implicit double precision (a-h,o-z)
       y=acos(x)
       pi=acos(-1.d0)
       acosd=180.d0*y / pi
       return
       end