c
c  Given r1(OH), r2(ClO), and r3(HCl), the potential energy and gradient 
c  is returned in a.u. (E relative to the HOCl minimum). This PES
c  corresponds to that of Bittererova, Bowman, Peterson, JCP 113, 6186 (2000)
c  for the X^1A' state of HOCl based on icMRCI+Q/CBS(DZ-QZ) calculations, 3/00
c
c...Notes:
c            the flag "ng" also in subroutine vhocl can be set to 0 or 1 to
c            turn on (1) and off (0) the calculation of analytical gradients
c
c            The analytical gradients are experimental!
c
c            there are no empirical corrections to this PES
c
c************************************************************************************
c.. sample driver program
       program hocl
       implicit real*8 (a-h,o-z)
       dimension grad(3)
       external cosd
c
       r1=1.8
       r2=3.2
       theta=103
       r3=sqrt(r1*r1+r2*r2-2.d0*r1*r2*cosd(theta))
c
       call vhocl(r1,r2,r3,energy,grad)
c
       print *, r1,r2,r3
       print *, 'energy= ',energy
       print *, 'grad(1)=',grad(1)
       print *, 'grad(2)=',grad(2)
       print *, 'grad(3)=',grad(3)
c
c  energy=  9.712413646872547E-05
c  grad(1)= -1.046776564682647E-02
c  grad(2)= 6.998191228601866E-04
c  grad(3)= 5.710774841197419E-04
c
       end
c************************************************************************************
c
       subroutine vhocl(r1,r2,r3,energy,g)
       implicit double precision  (a-h,o-z)
       parameter (hkcal=627.5096d0)
       common/stuff/ncall
       common/gradient/ng,grad(3)
       dimension g(3)
       data ncall/1/
c
c.. ng=1 (compute gradients with energy)
       ng=0
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)=1.d0
         g(2)=1.d0
         g(3)=1.d0
         return
       endif
c
        call calcpes(r1,r2,r3,e)
        ehocl=-163.96197168d0
c
       energy=(e-ehocl)/hkcal
c
       g(1)=grad(1)/hkcal
       g(2)=grad(2)/hkcal
       g(3)=grad(3)/hkcal
c
       return 
       end
c
       subroutine calcpes(r1,r2,r3,energy)
       implicit real*8 (a-h,o-z)
       common/gradient/ng,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
       energy=y2a+y2b+y2c+y3
c
       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='hocl.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,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,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
      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
      common/gradient/ng,grad(3)
      data a1,b1,c1,d1/3.247d0,0.9162d0,0.37491d0,0.02955d0/
      data tiny/1.d-14/
      external cosd
      two=2.d0
c
      seam=a1+b1*r1-c1*r1*r1+d1*r1**3
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=b1-two*c1*r1+3.d0*d1*r1*r1
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

        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
        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
       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,grad(3)
      data a1,b1/-6.12d0,3.50d0/
      data tiny/1.d-14/
      external cosd
      two=2.d0
c
       if((r2.gt.5.0d0).and.(r3.gt.3.4d0)) then
         dampci=0.05d0*(1.d0-exp(-(r2-5.0d0)))*(1.d0-exp(-(r3-3.4d0)))
       else
         dampci=0.d0
       endif
c
        seam=a1 + b1*r3
c
       ridge= r2*bcos+seam
c
       sw1=switch((bcos-cosd(170.0d0)),tans)
       sw2=switch((2.4d0-r2),tans)
       sw3=switch((r2-5.0d0),tans)
       sw4=switch((2.2d0-r3),tans)
       sw5=switch((r3-3.4d0),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.4d0)))
         ddampci3=0.05d0*exp(-(r3-3.4d0))*(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(170.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.2d0-r3),tans) +
     >        fullsw/sw5*dswitch((r3-3.4d0),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 real*8 (a-h,o-z)
        data pi/3.141592653589d0/
       
        cosd = dcos(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