      program h3tst
c--------------------------------------------------------------
c oct24/91  wkeogh
c  A simple test programme to verify that the H3 set of routines
c  produces the correct result on various computers.
c  This programme calculates the H3 energy for 6 geometries 
c  for which the known energy is supplied.
c  The derivatives are then check by comparing numerical derivatives
c  to the analytical derivatives for one geometry only.   
c--------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension x(18),Eknown(6),r(3),dVtot(3)
      dimension dVdr(3),dVdrn(3)
      parameter( npts=6,  dr=1.d-08 )
      data x/ 1.757d0, 1.757d0, 3.514d0,
     .        1.d0   , 2.d0   , 3.d0   ,
     .        3.d0   , 4.d0   , 5.d0   ,
     .        2.d0   , 2.d0   , 2.d0   ,
     .        1.d0   , 1.d0   , 1.d0   ,
     .        0.7d0  , 0.7d0  , 0.7d0  /
      data Eknown/ -0.1592929999, -0.0928127174,  -0.0561796335,
     .             -0.0739012770,  0.2122503352,   0.8559610809 /
      call idh3(h3id)
      write(6,*) ' H3 surface identification number: ',h3id
      id   = 1
      write(6,6000)
      do ipt=1,npts
	 ishift = (ipt-1)*3
	 r(1) = x(ishift+1)
	 r(2) = x(ishift+2)
	 r(3) = x(ishift+3)
         call h3tot(r,Vtot,dVtot,id)
	 diff = Vtot - Eknown(ipt)
	 write(6,6005) r(1),r(2),r(3),Vtot,Eknown(ipt),diff
      end do
c    check analytical versus numerical derivatives for one geometry:
c    save energy and analytical derivatives from last call
      r(1) = 0.95d0
      r(2) = 1.05d0
      r(3) = 1.85d0
      call h3tot(r,Vref,dVdr,id)
      do i=1,3
	 r(i) = r(i) + dr
	 call h3tot(r,Vtot,dVtot,id)
	 dVdrn(i) = (Vtot-Vref)/dr
	 r(i) = r(i) - dr
      end do
      write(6,6015) r
      write(6,6020) dVdr
      write(6,6025) dVdrn
 6000 format('    r1       r2      r3   ',
     .       '      Vtot           Eknown           diff  ')
 6005 format(3(f8.4,1x),2f16.10,2x,g12.6)
 6015 format('               r = ',3(1x,f12.6))
 6020 format('analytical dV/dr = ',3(1x,e12.6))
 6025 format('numerical  dV/dr = ',3(1x,e12.6))
      stop
      end
C
c   H3 surface evaluation routines
c   ==============================
c
      block data
c----------------------------------------------------------------------
c Note vector 'ipass' is just used to pass debugging print flags around:
c  if ipr.gt.0 print out some details of various contributions to Vtot
c  if ipr.gt.1  also print out some details about derivatives
c where ipr is set from various elements of ipass in various subroutines.
c----------------------------------------------------------------------
      common/pass/ipass(100)
      data ipass/100*0/
      end
c
      subroutine idh3(h3id)
c-----------------------------------------------------------------
c  h3 surface identification number is 706
c-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
	h3id = 706.0
      return
      end
c
      subroutine h3tot(r,Vtot,dVtot,id)
c-----------------------------------------------------------------
c Calculate total H3 potential from all of its parts
c if (id.gt.0) also calculate the derivative 
c all distances are in bohrs and all energies are in hartrees
c
c For a discussion of this surface, see:
c    A.I.Boothroyd, W.J.Keogh, P.G.Martin, M.R.Peterson 
c    Journal of Chemical Physics 95 pp4343-4359 (Sept15/91)
c Note: the surface parameter values as published lead to an anomolously
c       deep (several millihartree) van der Waals well for a very compact
c       H2 molecule (i.e., r < 0.8 bohr).  After that paper was submitted,
c       this problem was fixed and the corrected Cbend coefficients are
c       used in this version of the surface (version 706); this change
c       should have negligible effect, since such small H2 molecules are
c       unbound.  (The old coefficients are still in "subroutine vbcb"
c       but have been commented out.)
c
c any QUESTIONS/PROBLEMS/COMMENTS concerning this programme can be
c addressed to:      wkeogh@alchemy.chem.utoronto.ca
c (if necessary, via   boothroy@cita.utoronto.ca    boothroy@utordop.bitnet
c                 or   pgmartin@cita.utoronto.ca    pgmartin@utordop.bitnet )
c
c----------------------------------------------------------------------
c version of july27/91 ... surf706d.out Cbend values put in
c----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),dVtot(3)
      dimension dVlon(3),dVas(3),dVbnda(3),dVbndb(3),
     .          dCal(3), dCas(3),dCbnda(3),dCbndb(3)
      dimension dB1a(3),dB1b(3),dT(3),T(3)
      dimension dB2(3),dB3(6),expass(4)
c							{debugging:}
      dimension vb(2,25),cb(2,25),vbp(2,25),cbp(2,25)
      common/pass/ipass(100)
c
      ipr = ipass(1)
c
c FIRST, CHECK THAT IT IS A VALID H3 GEOMETRY, WITHIN TOLERANCE "DERR":
c
      derr = 1.d-5
      rhi = dmax1(r(1),r(2))
      if(r(3).ge.rhi)then
         rmid = rhi
         rhi = r(3)
      else
         rmid = dmax1( r(3) , dmin1(r(1),r(2)) )
      endif
      rlo = dmin1(r(1),r(2),r(3))
      if(rlo+rmid.gt.rhi+derr)
     .     stop ' STOP -- H3TOT Error: Not a valid H3 geometry. '
      if(rlo.lt.0.2d0)
     .     stop ' STOP -- H3TOT Error: shortest distance < 0.2 bohr. '
c
c  zero everything to avoid any 'funny' values:
      Vlon  = 0.d0
      Vas   = 0.d0
      Vbnda = 0.d0
      Vbndb = 0.d0
      Cal   = 0.d0
      Cas   = 0.d0
      cbnda = 0.d0
      cbndb = 0.d0
      do i=1,3
	 dVlon(i) = 0.d0
	 dVas(i)  = 0.d0
	 dVbnda(i)= 0.d0
	 dVbndb(i)= 0.d0
	 dCal(i)  = 0.d0
	 dCas(i)  = 0.d0
	 dCbnda(i)= 0.d0
	 dCbndb(i)= 0.d0
      end do  
c  zero the vb and cb arrays (used only for numerical derivatives)
      do i=1,25
	 vb(1,i) = 0.d0
	 vb(2,i) = 0.d0
	 cb(1,i) = 0.d0
	 cb(2,i) = 0.d0
      end do  
      call h3lond( r, Vlon, dVlon )
      call vascal( r, Vas, dVas )
c    now do any corrections required for compact geometries:
      call compac(r,icompc,T,dT)
c    oct.3/90 compact routines only called for compact geometries:
      if(icompc.ge.1)then
         call csym  ( r, cal, dCal )
         call casym ( r, Cas, dCas,  1 ,T,dT)
      end if
      call VbCb(r,icompc,T,dT,ipr,
     .          Vbnda,Vbndb,dVbnda,dVbndb,
     .          Cbnda,Cbndb,dCbnda,dCbndb)
c    add up the various parts of the potential:
      Vtot = Vlon + Vas + Vbnda + Vbndb + Cal + Cas + Cbnda + Cbndb
c    add up the various parts of the derivative:
      dVtot(1) =   dVlon(1) + dVas(1) + dVbnda(1) + dVbndb(1) 
     .           +  dCal(1) + dCas(1) + dCbnda(1) + dCbndb(1)
      dVtot(2) =   dVlon(2) + dVas(2) + dVbnda(2) + dVbndb(2) 
     .           +  dCal(2) + dCas(2) + dCbnda(2) + dCbndb(2)
      dVtot(3) =   dVlon(3) + dVas(3) + dVbnda(3) + dVbndb(3) 
     .           +  dCal(3) + dCas(3) + dCbnda(3) + dCbndb(3)

      if(ipr.gt.0) then
	  write(7,*)
	  write(7,*) 'h3tot enter --------------------------------'
          write(7,7000) r,icompc,vtot,
     .             vlon,vas,cal,cas,vbnda,vbndb,cbnda,cbndb
	  if(ipr.gt.1)then
		   write(7,7100) dVtot,dVlon,dVas,dVbnda,dVbndb,
     .                           dCal,dCas,dCbnda,dCbndb
	  end if
	  write(7,7400) (vb(1,i),i=1,5),(vb(2,i),i=1,5)
	  write(7,7500) (cb(1,i),i=1,9),(cb(2,i),i=1,9)
          write(7,7999) Vbnda,Vbndb,Cbnda,Cbndb,  
     .                  dVbnda,dvbndb,dCbnda,dCbndb
	  write(7,*) 'exiting subr.h3tot'
	  write(7,*) 'h3tot exit ---------------------------------'
      end if
      return
 7400 format('Vba values: ',5(1x,g12.6),/,'Vbb values: ',5(1x,g12.6))
 7500 format('Cba values: ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       'Cbb values: ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8),/,
     .       '            ',3(1x,f16.8))  
 7000 format(5x,'    r =',3(1x,f16.10),' icompac = ',i1,/,
     .       5x,'Vtot  = ',f18.12,/,
     .       5x,'Vlon  = ',f18.12,'       Vas   = ',g18.12,/,
     .       5x,'Cal   = ',g18.12,'       Cas   = ',g18.12,/,
     .       5x,'Vbnda = ',g18.12,'       Vbndb = ',g18.12,/,
     .       5x,'Cbnda = ',g18.12,'       Cbndb = ',g18.12)
 7100 format(' dVtot   = ',3(1x,g18.12),/,
     .       ' dVlon   = ',3(1x,g18.12),/,
     .       ' dVas    = ',3(1x,g18.12),/,
     .       ' dVbnda  = ',3(1x,g18.12),/,
     .       ' dVbndb  = ',3(1x,g18.12),/,
     .       ' dCal    = ',3(1x,g18.12),/,
     .       ' dCas    = ',3(1x,g18.12),/,
     .       ' dCbnda  = ',3(1x,g18.12),/,
     .       ' dCbndb  = ',3(1x,g18.12))
 7750 format(3(1x,g16.10))
 7751 format(/,3(1x,g16.10))
 7999 format('   Vbnda        Vbndb        Cbnda        Cbndb ',/,
     .        4(e12.6,2x),/,
     .       'dVbnda: ',3(e16.10,1x),/,
     .       'dVbndb: ',3(e16.10,1x),/,
     .       'dCbnda: ',3(e16.10,1x),/,
     .       'dCbndb: ',3(e16.10,1x))
      end
c
      subroutine triplet(r,E3,id)
C--------------------------------------------------------------------
c  version of oct4/90  surface626 values
C  H2 triplet curve and derivatives:
c  calculates triplet potential and first derivative
c  uses truhlar horowitz equation with our extension
c  uses the Johnson correction at short distances (r < rr)
c     if r .ge. rr         use modified t/h triplet equation
c     if r .le. rl         use the Johnson correction
c     in between           use the transition equation
C--------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension E3(3),E1(3)
      common/pass/ipass(100)
      parameter( rl=0.95d0, rr=1.15d0 )
      parameter( z1=1.d0, z2=2.d0, z4=4.d0)
c  triplet and johnson values from fit621:
      parameter(a1=-0.0253496194,a2=-29.2302126444,
     .          a3=-50.7225015503,a4= 2.0452676876 ,
     .          a5=-12.2408908509,a6= 1.6733157383 )
      parameter( C1=-0.4170298146519658,
     .           C2=-0.0746027774843370,C3= 0.4297899952237434 )
c surface parameters from fit601.out
c     parameter(a1=-0.66129429,a2=-1.99434198,a3=-2.37604328,
c    .          a4= 2.08107802,a5=-0.0313032510,a6=3.76546699,
c    .          C1=-0.4222590135447196,C2=-0.0731117796738824,
c    .          C3= 0.4295918082189010 )
      ipr = ipass(9)
      E3(3) = 0.d0
      if(r.ge.rr )then
c       modified truhlar/horowitz triplet equation:
         exdr = dexp( -a4*r )
         rsq  = r*r
         ra6  = r**(-a6)
         E3(1) = a1* ( a2 + r + a3*rsq + a5*ra6 )*exdr 
c       first derivative of triplet curve:
         ra61 = r**(-a6-z1)
         E3(2) = a1*exdr*
     .   ( z1 -a2*a4 +(z2*a3-a4)*r -a3*a4*rsq -a5*a6*ra61 -a4*a5*ra6 )
      end if
      if(r.lt.rr )then
          dr = r- rl  
          call vh2opt(r,e1,2)
	  if( r.le.rl ) then
c          Johnson triplet equation:
            E3(1) = E1(1) + c2*dr + c3
	    E3(2) = E1(2) + c2
          else                 
c          Transition equation:    
            E3(1) = E1(1) + c1*dr*dr*dr + c2*dr + c3
	    E3(2) = E1(2) + 3.d0*c1*dr*dr + c2
	  end if
      end if
      if(ipr.gt.0) write(7,7100) e3
 7100 format('  triplet:  E3 = ',3(1x,f12.8))
      return
      end
C
c
      subroutine vH2opt(r,E,ideriv)  
C-----------------------------------------------------------------
c  version of jul 09/90 ... super duper speedy version
C self-contained version of schwenke's H2 potential
C all distances in bohrs and all energies in Hartrees
C (1st deriv added on May 2 1989; 2nd deriv on May 28 1989)
C-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension E(3)
      common/pass/ipass(100)
      parameter(a0=0.03537359271649620,    a1= 2.013977588700072    ,
     .  a2= -2.827452449964767        ,    a3= 2.713257715593500    ,
     .  a4= -2.792039234205731        ,    a5= 2.166542078766724    ,
     .  a6= -1.272679684173909        ,    a7= 0.5630423099212294   ,
     .  a8= -0.1879397372273814       ,    a9= 0.04719891893374140  ,
     . a10= -0.008851622656489644     ,   a11= 0.001224998776243630 ,
     . a12= -1.227820520228028d-04    ,   a13= 8.638783190083473d-06,
     . a14= -4.036967926499151d-07    ,   a15= 1.123286608335365d-08,
     . a16= -1.406619156782167d-10 )
      parameter( r0=3.5284882d0,  DD=0.160979391d0,
     .           c6=6.499027d0,   c8=124.3991d0,    c10=3285.828d0)
c     Ediss = 0.174445d0
      E(1)  = -999.d0
      E(2)  = -999.d0
      E(3)  = -999.d0
      r2  = r  *r
      r3  = r2 *r
      r4  = r3 *r
      r5  = r4 *r
      r6  = r5 *r
      r7  = r6 *r
      r8  = r7 *r
      r9  = r8 *r
      r10 = r9 *r
      r11 = r10*r
      r12 = r11*r
      r13 = r12*r
      r14 = r13*r
      r15 = r14*r
      r02 = r0*r0
      r04 = r02*r02
      r06 = r04*r02
      rr2  = r2 + r02
      rr4  = r4 + r04
      rr6  = r6 + r06
      rr25 = rr2*rr2*rr2*rr2*rr2
c     general term:  a(i)*r(i-1),  i=0,16
      alphaR =  a0/r + a1
     .            + a2 *r   + a3 *r2  + a4 *r3  + a5 *r4  + a6 *r5
     .            + a7 *r6  + a8 *r7  + a9 *r8  + a10*r9  + a11*r10
     .            + a12*r11 + a13*r12 + a14*r13 + a15*r14 + a16*r15
      exalph = dexp(alphaR)
      vsr = DD*(exalph-1.d0)*(exalph-1.d0) - DD
      vlr =  -C6/rr6   -C8/(rr4*rr4)   -C10/rr25
      E(1)=  vsr +  vlr 
c
C  calculate first derivative if required:
      if(ideriv.ge.1)then        
         r3 = r2*r
         r5 = r4*r
         rr26 = rr25*rr2
         rr43 = rr4*rr4*rr4
c        general term:  (i-1)*a(i)*r**(i-2)  ,  i=0,16
	 dalphR = -a0/r2 + a2 + 2.d0*a3*r
     .          +  3.d0 *a4 *r2  +  4.d0 *a5 *r3  +  5.d0 *a6 *r4 
     .          +  6.d0 *a7 *r5  +  7.d0 *a8 *r6  +  8.d0 *a9 *r7 
     .          +  9.d0 *a10*r8  + 10.d0 *a11*r9  + 11.d0 *a12*r10
     .          + 12.d0 *a13*r11 + 13.d0 *a14*r12 + 14.d0 *a15*r13
     .          + 15.d0 *a16*r14
	 dvsr = 2.d0 *DD *(exalph-1.d0) *exalph *dalphR
	 dvlr =    6.d0*C6*r5 / (rr6*rr6)
     .          +  8.d0*C8*r3 / rr43  + 10.d0*C10*r / rr26
	 E(2) = dvsr + dvlr
      end if
c
C  calculate second derivative if required:
      if(ideriv.ge.2)then
         r10  = r6*r4
         rr27 = rr26*rr2
         rr44 = rr43*rr4
         rr62 = rr6*rr6
         rr63 = rr62*rr6
c        general term: (i-1)*(i-2)*a(i)*r**(i-3),    i=0,16
	 ddalph = 2.d0*a0/r3 +  2.d0*a3      +  6.d0*a4 *r
     .     + 12.d0*a5 *r2    + 20.d0*a6 *r3  + 30.d0*a7 *r4
     .     + 42.d0*a8 *r5    + 56.d0*a9 *r6  + 72.d0*a10*r7
     .     + 90.d0*a11*r8    +110.d0*a12*r9  +132.d0*a13*r10
     .     +156.d0*a14*r11   +182.d0*a15*r12 +210.d0*a16*r13
	 ddvsr = 2.d0 *DD *exalph 
     .   *( (2.d0*exalph-1.d0)*dalphR*dalphR + (exalph-1.d0)*ddalph )
	 ddvlr =- 72.d0* C6 *r10 / rr63 -96.d0 *C8 *r6 / rr44  
     .		-120.d0*C10 *r2  / rr27 +30.d0 *C6 *r4 / rr62
     .		+ 24.d0* C8 *r2  / rr43 +10.d0 *C10    / rr26
	 E(3) = ddvsr + ddvlr
      end if
      return
      end
c
      subroutine h3lond(r,Vlon,dVlon)
c----------------------------------------------------------------------
c version of may 12/90  ... derivatives corrected (0.5 changed to 0.25)
c calculates the h3 london terms and derivatives 
c modified oct 7/89 to include eps**2 term which rounds off the
c cusp in the h3 potential which occurs at equilateral triangle 
c configurations
c----------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      real*8 Q(3),J(3),Jt
      dimension r(3),e1(3),e3(3),esing(3),etrip(3),dVlon(3)
      dimension de1(3),de3(3),dJ1(3),dJ2(3),dJ3(3)
      common/pass/ipass(100)
      parameter( half=0.5d0, two=2.d0 , eps2=1.d-12 )
      ipr=ipass(2)
      do i=1,3
         call vh2opt(r(i),esing,2)
	  e1(i) = esing(1)
	 de1(i) = esing(2)
	 call triplet(r(i),etrip,2)
	 if(ipr.gt.0) write(7,7400) i,esing,etrip
	  e3(i) = etrip(1)
	 de3(i) = etrip(2)
         q(i)  = half*(E1(i) + E3(i))  
         j(i)  = half*(E1(i) - E3(i)) 
      end do   
      sumQ  =   q(1) + q(2) + q(3)
      sumJ  =   dabs( j(2)-j(1) )**2 
     .        + dabs( j(3)-j(2) )**2 
     .        + dabs( j(3)-j(1) )**2 
      Jt     = half*sumJ + eps2
      rootJt = dsqrt(Jt)
      Vlon   = sumQ - rootJt   
      if(ipr.gt.0) then
	 write(7,7410) sumq,sumj
         write(7,7420) vlon,rootJt
      end if
c  calculate the derivatives with respect to r(i):
      dVlon(1) = half*(dE1(1)+de3(1))
     .         - 0.25d0*(two*j(1)-j(2)-j(3))*(dE1(1)-de3(1))/rootJt
      dVlon(2) = half*(dE1(2)+de3(2))
     .         - 0.25d0*(two*j(2)-j(3)-j(1))*(dE1(2)-de3(2))/rootJt
      dVlon(3) = half*(dE1(3)+de3(3))
     .         - 0.25d0*(two*j(3)-j(1)-j(2))*(dE1(3)-de3(3))/rootJt
      if(ipr.gt.0) then
         write(7,7000) r,e1,e3,vlon
	 write(7,7100) q,j
	 write(7,7200) dVlon
      end if
 7000 format('             r = ',3(1x,f12.6),/,
     .       '            E1 = ',3(1x,f12.8),/,
     .       '            E3 = ',3(1x,f12.8),/,
     .       '          vlon = ',1x,f12.8)
 7100 format(13x,'q = ',3(1x,e12.6),/,13x,'j = ',3(1x,e12.6))
 7200 format('         dVlon = ',3(1x,g12.6))
 7400 format('from subr.london: ',/,
     .       '  using r',i1,':   Esinglet=',3(1x,f12.8),/,
     .       '         ',1x,'    Etriplet=',3(1x,f12.8))
 7410 format('         sumQ = ',g12.6,'        sumJ = ',g12.6)
 7420 format('         Vlon = ',f12.8,'      rootJt = ',g12.6)
      return
      end
c
      subroutine Vascal(rpass,Vas,dVas)
c------------------------------------------------------------------
c  version of oct11/90  fit632c.out values
c  version of oct5/90   surf626 values
c  calculate the asymmetric correction term and its derivatives
c  see equations [14] to [16] of truhlar/horowitz 1978 paper
c------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      common/pass/ipass(100)
      dimension dVas(3),dA(3),dS(3),rpass(3)
c--- Vasym values from fit632c ----------------------
      parameter(
     . aa1=0.3438222224E-02,aa2=0.1398145763E-02,
     . aa3=-.1923999449E-03,aa4=0.9712737075E-05,
     . aa5=-.1263794562E-06,aa6=0.5181432712E+00,
     . aa7=-.9487002995E-03 )
      ipr = ipass(3)
      r1 = rpass(1)
      r2 = rpass(2)
      r3 = rpass(3)
      R  = r1 + r2 + r3
      Rsq = R*R
      Rcu = Rsq*R
C   calculate the Vas term first (eq.14 of truhlar/horowitz)
      call acalc(r1,r2,r3,A,dA)
        A2 = A *A
        A3 = A2*A
        A4 = A3*A
        A5 = A4*A
        exp1 = dexp(-aa1*Rcu)
	exp6 = dexp(-aa6*R)
        S = aa2*A2 + aa3*A3 + aa4*A4 + aa5*A5
      Vas = S*exp1 +  aa7*A2 *exp6 / R
         dS(1)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(1)
         dVas(1) =  -3.d0*aa1*Rsq*S*exp1 + dS(1)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(1)*exp6/R
     .             -aa6*aa7*A2*exp6/R
         dS(2)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(2)
         dVas(2) =  -3.d0*aa1*Rsq*S*exp1 + dS(2)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(2)*exp6/R
     .             -aa6*aa7*A2*exp6/R
         dS(3)  = ( 2.d0*aa2*A  +3.d0*aa3*A2 
     .             +4.d0*aa4*A3 +5.d0*aa5*A4) * dA(3)
         dVas(3) =  -3.d0*aa1*Rsq*S*exp1 + dS(3)*exp1
     .             -aa7*A2*exp6/Rsq  + 2.d0*aa7*A*dA(3)*exp6/R
     .             -aa6*aa7*A2*exp6/R
      if(ipr.gt.0) write(7,7100) Vas,dS,dVas
 7100 format('    from subr.vascal:   Vas  = ',  1x,g12.6, /,
     .       '                        dS   = ',3(1x,g12.6),/,
     .       '                        dVas = ',3(1x,g12.6))
      return
      end
c
      subroutine acalc(r1,r2,r3,A,dA)
c---------------------------------------------------------------------
c  version of may 1, 1990
c  based on equations from notes date april 4, 1990
c  calculate the A value and its derivatives
c           A  = dabs[ (r1-r2)*(r2-r3)*(r3-r1) ]
c---------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension dA(3)
      common/pass/ipass(100)
      ipr = ipass(11)
      A = (r1-r2)*(r2-r3)*(r3-r1)
      dA(1) = ( -2.d0*r1 + r2 + r3 )*(r2-r3)
      dA(2) = ( -2.d0*r2 + r3 + r1 )*(r3-r1)
      dA(3) = ( -2.d0*r3 + r1 + r2 )*(r1-r2)
      if(A.lt.0.d0)then
         A = -A
         dA(1) = -dA(1)
         dA(2) = -dA(2)
         dA(3) = -dA(3)
      end if
      if(ipr.gt.0) write(7,7100) A,dA
 7100 format('Acalc enter ---------------------------------------',/,
     .       '     A = ',f15.10,'   dA = ',3(1x,g12.6),/,
     .       'Acalc exit ----------------------------------------')
      return
      end 
c
      subroutine compac(r,icompc,T,dT)
C---------------------------------------------------------------
c  version of may 14/90 ... calculates T and dT values also
c  decide whether or not this particular geometry is compact,
c  that is, are any of the three distances smaller than the
c  rr value from the johnson correction.
c-----------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),T(3),dT(3)
      common/pass/ipass(100)
      parameter( rr = 1.15d0, rp = 1.25d0 )
c  first see if this is a compact geometry:
      do i=1,3
	 T(i) = 0.d0
	 dT(i)= 0.d0
      end do   
      icompc = 0
      if(r(1).lt.rr) icompc = icompc + 1
      if(r(2).lt.rr) icompc = icompc + 1
      if(r(3).lt.rr) icompc = icompc + 1
      if(icompc.eq.0) return
c calculate the T(i) values:
      do i=1,3
      if(r(i).lt.rr) then
	 top = rr-r(i)
	 bot = rp-r(i)
	 top2 = top  * top
	 top3 = top2 * top
	 bot2 = bot  * bot
         T(i)  = top3/bot
         dT(i) = -3.d0*top2/bot + top3/bot2
      end if
      end do   
      return
      end
C
      subroutine casym ( r,Cas,dCas,id,T,dT )
c---------------------------------------------------------------
c  version of sept14/90 
c  the compact asymmetric correction term and derivatives
c---------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),T(3)
      dimension dPr(3),dT(3),dsumT(3),dterm1(3),determ(3),dseries(3),
     .          dCas(3),dA(3)
      common/pass/ipass(100)
c---- Casym values from fit633a   
      parameter( 
     . u1=0.1537481166E+00, u2=0.2745950036E+00, u3=0.7501206780E-02,
     . u4=0.3119136023E+01, u5=0.9969170798E+00, u6=-.3373682823E+01,
     . u7=0.6807215913E+00, u8=-.4920325491E-01, u9=-.3919467989E+01,
     .u10=0.5085532326E+01,u11=-.1415264778E+01,u12=0.1138681785E+00,
     .u13=-.1525367566E-02)
      ipr = ipass(6)
      cas = 0.d0
      dCas(1) = 0.d0
      dCas(2) = 0.d0
      dCas(3) = 0.d0
      call acalc(r(1),r(2),r(3),A,dA)
      A2 = A*A
c     if(A.eq.0.d0) return
      sumT = T(1) + T(2) + T(3)
      Sr   = r(1) + r(2) + r(3)
      Pr   = r(1) * r(2) * r(3)
      Sr2  = Sr*Sr
      Sr3  = Sr2*Sr
      Pr2  = Pr*Pr
      Pr3  = Pr2*Pr
c    write out the series explicitly:
      series = 1.d0 + u4/Pr2 +  u5/Pr + u6 +   u7*Pr +  u8*Pr2
     .           + A*(u9/Pr2 + u10/Pr + u11 + u12*Pr + u13*Pr2)
      term1  = u1/Pr**u2
      eterm  = dexp(-u3*Sr3)
      Cas    = sumt * a2 * term1 * series * eterm
      if(ipr.gt.0) then
	    write(7,*)
            write(7,*) 'Casym: ------------------'
	    write(7,*) ' T(i) = ',T
	    write(7,*) 'dT(i) = ',dT
	    write(7,*) ' id = ',id
	    write(7,7400) series,term1,eterm,Cas
	    write(7,*) ' Cas = ',Cas
      end if
c			{calculate the derivative as well:}
      if(id.gt.0)then
	    dPr(1) = r(2)*r(3)
	    dPr(2) = r(3)*r(1)
	    dPr(3) = r(1)*r(2)
	 do i=1,3
	    dsumt(i) = dt(i)
	    dterm1(i)=  -1.d0*u1*u2*Pr**(-u2-1.d0)*dPr(i)
	    determ(i)= eterm *(-3.d0*u3*Sr2)
	    dseries(i)=
     .        dPr(i)*( -2.d0*u4/Pr3  -u5/Pr2 +u7  +2.d0*u8*Pr    )
     .        +dA(i)*( u9/Pr2 +u10/Pr +u11 +u12*Pr +u13*Pr2    )
     .        +A*dPr(i)*( -2.d0*u9/Pr3 -u10/Pr2 +u12 + 2.d0*u13*Pr )
	    dCas(i) = 
     .        dsumt(i)      * a2    * term1 * series * eterm
     .      + 2.d0*a*dA(i)  * sumt  * term1 * series * eterm
     .      + dterm1(i)     * sumt  * a2    * series * eterm
     .      + dseries(i)    * sumt  * a2    * term1  * eterm
     .      + determ(i)     * sumt  * a2    * term1  * series 
         end do   
	 if(ipr.gt.0) then
	    write(7,7500) dPr,dt,dterm1,determ,dseries,dCas
	    write(7,*) ' dCas = ',dCas
	 end if
      end if
      return
 7400 format('subr.Cas:   series = ',g12.6,'     term1 = ',g12.6,/,
     .       '             eterm = ',g12.6,'      Cas  = ',g12.6)
 7500 format('    derivatives: dPr = ',3(1x,g12.6),/,
     .       '                 dT  = ',3(1x,g12.6),/,
     .       '              dterm1 = ',3(1x,g12.6),/,
     .       '              determ = ',3(1x,g12.6),/,
     .       '             dseries = ',3(1x,g12.6),/,
     .       '              dCas   = ',3(1x,g12.6))  
      end
c
      subroutine csym(r,Cal,dCal)
c---------------------------------------------------------------
c  version of oct12/90  surf636 values
c  calculate the 'compact all' correction term and derivatives
c  a correction term (added sept 11/89), which adds a small
c  correction to all compact geometries
c---------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension r(3),t(3),dCal(3),G(3)
      dimension dG(3),sumv(3)
      common/pass/ipass(100)
      parameter( rr=1.15d0, rp=1.25d0 )
c---- Csym values from fit633a                                                  
      parameter( 
     . v1=-.2070776049E+00, v2=-.5350898737E+00, v3=0.1011861942E-01)
      cal   = 0.d0
      ipr   = ipass(5)
      Sr    = r(1)+r(2)+r(3)
      Sr2   = Sr*Sr
      Sr3   = Sr*Sr2  
      exp3  = dexp( -v3*Sr3 )
      dexp3 = -3.d0*v3*Sr2*exp3
      do i=1,3
         ri      = r(i)
	 rrri    = rr-ri
	 rrri2   = rrri*rrri
	 rrri3   = rrri*rrri2
	 rpri    = rp-ri
	 rpri2   = rpri*rpri
         G(i)    = 0.d0
	 dG(i)   = 0.d0
	 sumv(i) = v1+v1*v2*ri
         if(ri.lt.rr) then
	    G(i)  = (rrri3/rpri) *sumv(i)
  	    dG(i) = (rrri3/rpri2)*sumv(i) 
     .                -3.d0*(rrri2/rpri)*sumv(i)
     .                + (rrri3/rpri)*v1*v2
	 end if
      end do   
      sumG = G(1) + G(2) + G(3)
      cal  = sumG*exp3
      dCal(1) = dG(1)*exp3 + sumG*dexp3
      dCal(2) = dG(2)*exp3 + sumG*dexp3
      dCal(3) = dG(3)*exp3 + sumG*dexp3
      if(ipr.gt.0)then
	 write(7,*)
         write(7,7100) Cal
	 write(7,7000) G,dG,sumv,dCal
	 write(7,*) 'exiting subr.csym'
      end if
 7100 format('Csym enter ----------------------------------------',/,
     .       '    Cal = ',1x,g20.14)
 7000 format('      G = ',3(1x,g12.6),/,
     .       '     dG = ',3(1x,g12.6),/,
     .       '   sumv = ',3(1x,g12.6),/,
     .       '   dCal = ',3(1x,g16.10),/,
     .       'Csym exit -----------------------------------------')
      return
      end
c
      subroutine VbCb(rpass,icompc,T,dT,ipr,
     .                Vbnda,Vbndb,dVbnda,dVbndb,
     .                Cbnda,Cbndb,dCbnda,dCbndb)
C---------------------------------------------------------------
c  jul24/91 fit705 cbend values added
c  feb27/91 modified to match equation in h3 paper more closely
c  nov 4/90 Vbend coefficients now a,g  Cbend coeff's still c,d
c  in this version, the derivatives are always calculated
c  subroutines vbend, cbend and B1ab all combined into this one
c  module in order to improve efficiency by not having to pass
c  around B1a,B1b,B2,B3a,B3b functions and derivatives
c   B1a = 1 - sum of [ P1(cos(theta(i))) ]
c   B1b = 1 - sum of [ P3(cos(theta(i))) ]
c--------------------------------------------------------------
      implicit double precision (a-h,o-z)
c    arrays required for B1,B2,B3 calculations:
      dimension rpass(3),dB1a(3),dB1b(3),th(3),
     .          dB2(3),dB3(6)
c    arrays required for Vbend calculations:
      dimension dVbnda(3),dVbndb(3)
      dimension dVb1(3),dvb2(3),dvb3(3),dvb4(3),dvb5(3),dVbnd(3)
      dimension Vbnd(2)
c			{passed out for debugging only:}
      dimension vb(2,25)
c    arrays required for Cbend calculations:
      dimension Cbnds(2),dCbnds(2,3)
      dimension dCb1(3),dCb2(3),dCb3(3),dCb4(3),dCb5(3),dCb6(3),
     .          dCb7(3),dCb8(3)
      dimension T(3),dP(3)
      dimension dT(3),dCbnda(3),dCbndb(3)
c					{passed out for debugging only:}
      dimension cb(2,25)
c						{3/8 and 5/8:}
      parameter( z58=0.625d0, z38=0.375d0 )
c   parameters required for Vbend calculations:
      parameter(beta1=0.52d0, beta2=0.052d0, beta3=0.79d0 ) 
c-- Vbenda values from fit634b                                               
      parameter( 
     .a11=-.2918252280E+03,a12=0.2164569141E+02,a13=-.4005497543E+00,
     .a21=-.2890774947E+01,a22=0.1032032542E+02,a23=0.2681748056E+02,
     .a24=0.2633751055E+02,a31=0.6180641351E+01,a32=0.5037667041E+01,
     .a41=-.1125570079E+00,a42=-.3176529304E-02,a43=0.9068915355E+00,
     .a44=-.7228418516E+00,a51=0.2785898232E+03,a52=0.4764442446E+02,
     .a53=0.8621545662E+01)
c-- Vbendb values from fit634b                                               
      parameter( 
     .g11=-.4241912309E+02,g12=0.2951365281E+01,g13=-.4840201514E-01,
     .g21=-.1159168549E-03,g22=0.8688003567E-02,g23=-.3486923900E-02,
     .g24=0.8312212839E-03,g31=0.5621810473E-01,g32=-.9776930747E-02,
     .g41=-.1178456251E-01,g42=0.3491086729E-02,g43=0.7430516993E-01,
     .g44=-.9643636957E-01,g51=0.4735533782E+02,g52=0.3001038808E+01,
     .g53=0.1896630453E+01)
c parameters required for Cbend calculations:
c-- Cbenda values from fit702c -------------------                          
cx    parameter( 
cx   .c11=0.3408013485E+04,c12=-.2569675505E+04,c13=-.4171156109E+03,
cx   .c14=-.6179297444E+04,c15=0.3235223485E+04,c21=-.3409049347E+03,
cx   .c22=0.9612365597E+02,c23=0.4104361425E+03,c24=0.6238298637E+03,
cx   .c31=0.2655770314E+04,c32=0.5314085916E+03,c41=0.2222002477E+02,
cx   .c42=-.6611541362E+01,c43=-.2541255632E+03,c44=0.7752854954E+02,
cx   .c51=0.9613217362E+04,c52=-.7010043704E+03,c53=0.1406657509E+03,
cx   .c61=-.1677846467E+04,c62=0.7411196008E+04,c63=-.1104074171E+05,
cx   .c71=-.1077399844E+03,c72=-.6598646429E+01,c73=0.2035775472E+00,
cx   .c74=0.1097670442E+04,c75=0.4873456879E+02,c81=0.1005135589E+03,
cx   .c82=-.1972894789E+02,c83=-.8495015577E+04,c84=-.6106570657E+02)
c-- Cbendb values from fit702c -------------------                          
cx    parameter( 
cx   .d11=0.8450990097E+03,d12=0.1466587231E+03,d13=0.1239413655E+02,
cx   .d14=-.5649327722E+03,d15=0.3112693920E+03,d21=0.8029625257E-01,
cx   .d22=0.2379962852E-01,d23=0.1233566879E+00,d24=-.3582414079E-01,
cx   .d31=-.1012436386E+03,d32=0.8904152887E-02,d41=0.8740633605E+01,
cx   .d42=-.2937137124E+01,d43=-.5497174488E+02,d44=0.1656850593E+02,
cx   .d51=0.3859496208E+05,d52=0.8458336162E+02,d53=-.1431690434E+02,
cx   .d61=0.5487645288E+02,d62=-.1954058540E+03,d63=-.7226235345E+02,
cx   .d71=-.1895543609E+02,d72=0.1281231795E+00,d73=0.1291641223E+00,
cx   .d74=-.3052282568E+03,d75=-.6194562183E+01,d81=-.1746535079E+02,
cx   .d82=-.8402855300E+00,d83=-.3878264141E+05,d84=0.5655392575E+01)
c-- Cbenda values from fit705a.out -------------                             
cxx   parameter( 
cxx  .c11=0.7138332997E+04,c12=-.3941467854E+04,c13=0.1214839191E+03,
cxx  .c14=-.9715171923E+04,c15=0.4549112075E+04,c21=-.7153573233E+02,
cxx  .c22=0.1410801562E+03,c23=0.5310286989E+03,c24=0.6632669305E+03,
cxx  .c31=0.9648491179E+03,c32=0.1336994724E+03,c41=0.3986106227E+01,
cxx  .c42=-.2071394753E+01,c43=-.1304296435E+03,c44=0.4433517239E+02,
cxx  .c51=-.2160974730E+05,c52=-.5040404494E+03,c53=0.1150251071E+03,
cxx  .c61=-.1900650362E+04,c62=0.8065979973E+04,c63=-.1022333250E+05,
cxx  .c71=-.4023633922E+02,c72=-.5593869000E+01,c73=0.2901552041E+00,
cxx  .c74=0.1332299065E+04,c75=0.4830933142E+02,c81=0.1527546398E+03,
cxx  .c82=-.2277696550E+02,c83=0.2214880233E+05,c84=-.5951714583E+02)
c-- Cbendb values from fit705a.out -------------                             
cxx   parameter( 
cxx  .d11=0.1733880674E+04,d12=-.5105474641E+03,d13=0.5739180116E+02,
cxx  .d14=-.1549259263E+04,d15=0.6494865366E+03,d21=-.1180847496E+00,
cxx  .d22=0.3847529877E-01,d23=0.1035799123E+00,d24=-.2783059123E-01,
cxx  .d31=-.4049790116E+02,d32=0.3623102614E-01,d41=0.2132254354E+01,
cxx  .d42=-.1337841873E+01,d43=-.2911759497E+02,d44=0.9698950726E+01,
cxx  .d51=0.7851192036E+05,d52=0.1244750132E+03,d53=-.2157749426E+02,
cxx  .d61=0.1945719436E+02,d62=-.4336295493E+02,d63=0.8025442515E+02,
cxx  .d71=0.5232787178E+01,d72=-.2082854432E+01,d73=0.1740967673E+00,
cxx  .d74=-.3069559809E+02,d75=0.9099007208E+01,d81=-.6088433744E+01,
cxx  .d82=-.9840250626E+00,d83=-.7880550115E+05,d84=-.9084945347E+01)
c-- Cbenda values from fit706d.out 
      parameter( 
     .c11=0.7107064647E+04,c12=-.3728421267E+04,c13=0.1757654624E+03,
     .c14=-.9725998132E+04,c15=0.4665086074E+04,c21=-.4435165986E+02,
     .c22=0.1604477309E+03,c23=0.5805142925E+03,c24=0.6892349445E+03,
     .c31=0.6581730442E+03,c32=0.6078389739E+02,c41=0.2885182566E+01,
     .c42=-.1728169916E+01,c43=-.1119535503E+03,c44=0.4052536250E+02,
     .c51=0.2540505673E+03,c52=-.5762083627E+03,c53=0.1295901320E+03,
     .c61=-.2131706075E+04,c62=0.9084452020E+04,c63=-.1138253963E+05,
     .c71=-.3964298833E+02,c72=-.5019979693E+01,c73=0.2906541488E+00,
     .c74=0.1212943686E+04,c75=0.4140463415E+02,c81=0.1752855549E+03,
     .c82=-.2496320107E+02,c83=0.3765413052E+03,c84=-.5480488130E+02)
c-- Cbendb values from fit706d.out 
      parameter( 
     .d11=0.1917166552E+04,d12=-.6542563392E+03,d13=0.6793758367E+02,
     .d14=-.1694968847E+04,d15=0.6866649703E+03,d21=-.2137567948E+00,
     .d22=0.4975938228E-01,d23=0.9364998295E-01,d24=-.2444320779E-01,
     .d31=-.2863126914E+02,d32=0.5443219625E-01,d41=0.9673956120E+00,
     .d42=-.1160159706E+01,d43=-.2424199759E+02,d44=0.8569424490E+01,
     .d51=-.6517635862E+04,d52=0.1518098147E+03,d53=-.2706514366E+02,
     .d61=0.4308392956E+02,d62=-.1234851732E+03,d63=0.2320626055E+03,
     .d71=0.1049541418E+02,d72=-.2424169341E+01,d73=0.1745646946E+00,
     .d74=0.2603615561E+02,d75=0.1345799970E+02,d81=-.6653710513E+01,
     .d82=-.2576854447E+00,d83=0.6172608425E+04,d84=-.1328142473E+02)
c   feb27/91 new c51 = c51+c83;   new d51 = d51+d83
      cx1 = c51 + c83
      dx1 = d51 + d83
      r1 = rpass(1)
      r2 = rpass(2)
      r3 = rpass(3)
      t1 = r1*r1 -r2*r2 -r3*r3 
      t2 = r2*r2 -r3*r3 -r1*r1 
      t3 = r3*r3 -r1*r1 -r2*r2 
c  calculate the cosines of the three internal angles:
      c1 = t1 / (-2.d0*r2*r3)
      c2 = t2 / (-2.d0*r3*r1)
      c3 = t3 / (-2.d0*r1*r2)
      sum = c1 + c2 + c3
      B1a = 1.d0 - sum         
      c1cube = c1*c1*c1
      c2cube = c2*c2*c2
      c3cube = c3*c3*c3
      cos3t1 = 4.d0*c1cube - 3.d0*c1
      cos3t2 = 4.d0*c2cube - 3.d0*c2
      cos3t3 = 4.d0*c3cube - 3.d0*c3
      sumb   = cos3t1 + cos3t2 + cos3t3 
      B1b    = 1.d0 - ( z58*sumb + z38*sum )
c  calculate derivatives if desired
         dc1dr1 = -r1/(r2*r3)
         dc2dr2 = -r2/(r1*r3)
         dc3dr3 = -r3/(r1*r2)
         dc1dr2 = ( t1/(r2*r2) + 2.d0 )/(2.d0*r3)
         dc1dr3 = ( t1/(r3*r3) + 2.d0 )/(2.d0*r2)
         dc2dr1 = ( t2/(r1*r1) + 2.d0 )/(2.d0*r3)
         dc2dr3 = ( t2/(r3*r3) + 2.d0 )/(2.d0*r1)
         dc3dr1 = ( t3/(r1*r1) + 2.d0 )/(2.d0*r2)
         dc3dr2 = ( t3/(r2*r2) + 2.d0 )/(2.d0*r1)
         db1a(1) = -1.d0*( dc1dr1 + dc2dr1 + dc3dr1 )
         db1a(2) = -1.d0*( dc1dr2 + dc2dr2 + dc3dr2 )
         db1a(3) = -1.d0*( dc1dr3 + dc2dr3 + dc3dr3 )
         d1 = 12.d0*c1*c1 - 3.d0
         d2 = 12.d0*c2*c2 - 3.d0
         d3 = 12.d0*c3*c3 - 3.d0
         db1b(1) = -z58*( d1*dc1dr1 + d2*dc2dr1 + d3*dc3dr1 )
     .             -z38*(    dc1dr1 +    dc2dr1 +    dc3dr1 )
         db1b(2) = -z58*( d1*dc1dr2 + d2*dc2dr2 + d3*dc3dr2 )
     .             -z38*(    dc1dr2 +    dc2dr2 +    dc3dr2 )
         db1b(3) = -z58*( d1*dc1dr3 + d2*dc2dr3 + d3*dc3dr3 )
     .             -z38*(    dc1dr3 +    dc2dr3 +    dc3dr3 )
d     if(ipr.gt.1)then
d         write(7,*)
d         write(7,*) ' from subr.vbcb -----------------------------'
d         write(7,6000) rpass
d         write(7,6010) b1a,b1b
d         write(7,6020) db1a,db1b
d         if(ipr.gt.2)then
d            write(7,7000) dc1dr1,dc1dr2,dc1dr3,
d    .                     dc2dr1,dc2dr2,dc2dr3,
d    .                     dc3dr1,dc3dr2,dc3dr3
d            write(7,7010) t1,t2,t3,c1,c2,c3,d1,d2,d3
d         end if
d     end if
c
c  calculate the quantities used by both Vbend and Cbend:
c
      R    = r1 + r2 + r3
      Rsq  = R*R
      B2   = 1.d0/r1 + 1.d0/r2 + 1.d0/r3
      B3   = (r2-r1)*(r2-r1) + (r3-r2)*(r3-r2) + (r1-r3)*(r1-r3)
      eps2 = 1.d-12
      B3b  = dsqrt( B3 + eps2 )
	 dB2(1) = -1.d0/(r1*r1)
	 dB2(2) = -1.d0/(r2*r2)
	 dB2(3) = -1.d0/(r3*r3)
	 dB3(1) = 4.d0*r1 -2.d0*r2 -2.d0*r3
	 dB3(2) = 4.d0*r2 -2.d0*r3 -2.d0*r1
	 dB3(3) = 4.d0*r3 -2.d0*r1 -2.d0*r2
	 dB3(4) = 0.5d0*dB3(1)/B3b
	 dB3(5) = 0.5d0*dB3(2)/B3b
	 dB3(6) = 0.5d0*dB3(3)/B3b
         exp1   = dexp(-beta1*R)
         exp2   = dexp(-beta2*Rsq)
         exp7   = dexp(-beta3*R)
	 dexp1  = -beta1*exp1
	 dexp2  = -2.d0*beta2*R*exp2
	 dexp7  = -beta3*exp7
c   do the vbnda calculations:
         B1    = B1a
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         asum  = a11 +  a12*R +  a13*Rsq 
         bsum  = a21*B12 +  a22*B13 +  a23*B14 +  a24*B15
	 csum  = a31*B1*exp1 +  a32*B12*exp2
         dsum1 = a41*exp1 +  a42*exp2
         dsum2 = a43*exp1 +  a44*exp2
         fsum  = a51 +  a52*R +  a53*Rsq 
         vb(1,1) = B1*asum*exp1 
         vb(1,2) = bsum * exp2
         vb(1,3) = B2*csum
         vb(1,4) = B1*B3*dsum1 + B1*B3b*dsum2
         vb(1,5) = B1* fsum  *exp7 
         Vbnd(1) = vb(1,1)+vb(1,2)+vb(1,3)+vb(1,4)+vb(1,5)
c
	 dasum   = a12+2.d0*a13*R
	 dbsum   =  2.d0* a21*B1  + 3.d0* a22*B12
     .            + 4.d0* a23*B13 + 5.d0* a24*B14
	 ddsum1  =  a41*dexp1 +  a42*dexp2
	 ddsum2  =  a43*dexp1 +  a44*dexp2
	 dfsum   =  a52 +2.d0*a53*R
         dVb1(1) = dB1a(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(2) = dB1a(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(3) = dB1a(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
c 
         dVb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2
         dVb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2
         dVb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2
c
c calculate the Vb3 derivatives:
         dVb3(1)=dB2(1)*csum + B2*(  a31*dB1a(1)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(1)*exp2 +  a32*B12*dexp2)
         dVb3(2)=dB2(2)*csum + B2*(  a31*dB1a(2)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(2)*exp2 +  a32*B12*dexp2)
         dVb3(3)=dB2(3)*csum + B2*(  a31*dB1a(3)*exp1 +  a31*B1*dexp1
     .                  +2.d0* a32*B1*dB1a(3)*exp2 +  a32*B12*dexp2)
c
c calculate the Vb4 derivatives (may 27/90):
         dVb4(1)= dB1a(1)*B3 *dsum1+  B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
         dVb4(2)= dB1a(2)*B3 *dsum1+  B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
         dVb4(3)= dB1a(3)*B3 *dsum1+  B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .          + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c
c calculate the Vb5 derivatives:
	 dVb5(1) = dB1a(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(2) = dB1a(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(3) = dB1a(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
c
c calculate the overall derivatives:
         dVbnda(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1)
         dVbnda(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2)
         dVbnda(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3)
	 Vbnda = Vbnd(1)
c
c---- now do the vbendb calculations --------
         B1    = B1b
         B12   = B1*B1
         B13   = B12*B1
         B14   = B13*B1
         B15   = B14*B1
         asum  = g11 +  g12*R +  g13*Rsq 
         bsum  = g21*B12 +  g22*B13 +  g23*B14 +  g24*B15
	 csum  = g31*B1*exp1 +  g32*B12*exp2
         dsum1 = g41*exp1 +  g42*exp2
         dsum2 = g43*exp1 +  g44*exp2
         fsum  = g51 +  g52*R +  g53*Rsq 
         vb(2,1) = B1*asum*exp1 
         vb(2,2) = bsum * exp2
         vb(2,3) = B2*csum
         vb(2,4) = B1*B3*dsum1 + B1*B3b*dsum2
         vb(2,5) = B1* fsum  *exp7 
         Vbnd(2) = vb(2,1)+vb(2,2)+vb(2,3)+vb(2,4)+vb(2,5)
c
	 dasum   = g12+2.d0*g13*R
	 dbsum   =  2.d0* g21*B1  + 3.d0* g22*B12
     .             +4.d0* g23*B13 + 5.d0* g24*B14
	 ddsum1  =  g41*dexp1 +  g42*dexp2
	 ddsum2  =  g43*dexp1 +  g44*dexp2
	 dfsum   =  g52 +2.d0*g53*R
         dVb1(1) = dB1b(1)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(2) = dB1b(2)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
         dVb1(3) = dB1b(3)*asum*exp1 + B1*dasum*exp1 +B1*asum*dexp1
c 
         dVb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2
         dVb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2
         dVb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2
c
c calculate the Vb3 derivatives:
         dVb3(1)=dB2(1)*csum + B2*(  g31*dB1b(1)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(1)*exp2 +  g32*B12*dexp2)
         dVb3(2)=dB2(2)*csum + B2*(  g31*dB1b(2)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(2)*exp2 +  g32*B12*dexp2)
         dVb3(3)=dB2(3)*csum + B2*(  g31*dB1b(3)*exp1 +  g31*B1*dexp1
     .                  +2.d0* g32*B1*dB1b(3)*exp2 +  g32*B12*dexp2)
c
c calculate the Vb4 derivatives (may 27/90):
         dVb4(1)= dB1b(1)*B3 *dsum1+  B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
         dVb4(2)= dB1b(2)*B3 *dsum1+  B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
         dVb4(3)= dB1b(3)*B3 *dsum1+  B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .          + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c
c calculate the Vb5 derivatives:
	 dVb5(1) = dB1b(1)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(2) = dB1b(2)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
	 dVb5(3) = dB1b(3)*fsum*exp7 + B1*dfsum*exp7 + B1*fsum*dexp7
c
c calculate the overall derivatives:
         dVbndb(1) = dvb1(1)+dvb2(1)+dvb3(1)+dvb4(1)+dvb5(1)
         dVbndb(2) = dvb1(2)+dvb2(2)+dvb3(2)+dvb4(2)+dvb5(2)
         dVbndb(3) = dvb1(3)+dvb2(3)+dvb3(3)+dvb4(3)+dvb5(3)
	 Vbndb = Vbnd(2)
c
c  now calculate the Cbend correction terms:
c  Cbnda uses the B1a formula and Cbndb uses the B1b formula
      if(icompc.eq.0) return
      sumT  = T(1) + T(2) + T(3)
      Rcu   = Rsq*R
      P     = r1*r2*r3
      Psq   = P*P
      Pcu   = Psq*P
      dP(1) = r2*r3
      dP(2) = r3*r1
      dP(3) = r1*r2
      Cbnds(1)  = 0.d0
      Cbnds(2)  = 0.d0
      dCbnda(1) = 0.d0
      dCbnda(2) = 0.d0
      dCbnda(3) = 0.d0
      dCbndb(1) = 0.d0
      dCbndb(2) = 0.d0
      dCbndb(3) = 0.d0
c  calculate exponentials and derivatives (common to cbnda and cbndb):
            exp7  = dexp( -beta2*Rcu ) 
	    dexp7 = -3.d0*beta2*Rsq*exp7
c ----- calculate the Cbnda correction term ----------------
       B1    = B1a
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       asum  = c11 + c12*R + c13*Rsq + c14/R + c15/Rsq 
       bsum  = c21*B12 + c22*B13 + c23*B14 + c24*B15
       csum  = c31*B1 *exp1  + c32*B12*exp2  
       dsum1 = c41*exp1 + c42*exp2
       dsum2 = c43*exp1 + c44*exp2
       fsum  = cx1 + c52*R + c53*Rsq 
       gsum  = c61 + c62/R + c63/Rsq 
       aasum = c71 + c72*P + c73*Psq + c74/P + c75/Psq 
       ffsum =       c81*P + c82*Psq + c84/Psq 
       dasum = c12 + 2.d0*c13*R -c14/Rsq -2.d0*c15/Rcu
       dbsum = 2.d0*c21*B1 +3.d0*c22*B12 +4.d0*c23*B13 +5.d0*c24*B14
       ddsum1= c41*dexp1 + c42*dexp2
       ddsum2= c43*dexp1 + c44*dexp2
       dfsum = c52 + 2.d0*c53*R
       dgsum = -c62/Rsq - 2.d0*c63/Rcu
       daasum = c72 +2.d0*c73*P -c74/Psq -2.d0*c75/Pcu
       dffsum = c81 + 2.d0*c82*P -2.d0*c84/Pcu
         cb(1,1)  = B1*asum * exp1  /P
         cb(1,2)  = bsum * exp2  
         cb(1,3)  = B2*csum 
         cb(1,4)  = B1*B3*dsum1 + B1*B3b*dsum2
         cb(1,5)  = B1 * fsum * exp7 /P
         cb(1,6)  = B1 * gsum * exp7  
         cb(1,7)  = B1 * aasum * exp2 
         cb(1,8)  = B1 * ffsum * exp7 
         Cbnds(1) = cb(1,1)+cb(1,2)+cb(1,3)+cb(1,4)+cb(1,5)
     .                     +cb(1,6)+cb(1,7)+cb(1,8)
c       calculate the derivatives:
	 dCb1(1) = dB1a(1)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq
	 dCb1(2) = dB1a(2)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq
	 dCb1(3) = dB1a(3)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq
c
	 dCb2(1) = dbsum*dB1a(1)*exp2 + bsum*dexp2
	 dCb2(2) = dbsum*dB1a(2)*exp2 + bsum*dexp2
	 dCb2(3) = dbsum*dB1a(3)*exp2 + bsum*dexp2
c
	 dCb3(1) = dB2(1)*csum + B2*( c31*dB1a(1)*exp1 +c31*B1*dexp1
     .            +c32 *2.d0*B1*dB1a(1)*exp2 + c32 *B12*dexp2 )
	 dCb3(2) = dB2(2)*csum + B2*( c31 *dB1a(2)*exp1 +c31 *B1*dexp1
     .            +c32 *2.d0*B1*dB1a(2)*exp2 + c32 *B12*dexp2 )
	 dCb3(3) = dB2(3)*csum + B2*( c31 *dB1a(3)*exp1 +c31 *B1*dexp1
     .            +c32 *2.d0*B1*dB1a(3)*exp2 + c32 *B12*dexp2 )
c          dCb4 equations may 27/90:
	 dCb4(1) = dB1a(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
	 dCb4(2) = dB1a(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
	 dCb4(3) = dB1a(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .           + dB1a(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c          dCb5 equations may 27/90: (corrected jun01)
	 dCb5(1) = dB1a(1)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq
	 dCb5(2) = dB1a(2)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(2)/Psq
	 dCb5(3) = dB1a(3)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*Fsum*exp7*dP(3)/Psq
c
	 dCb6(1) = dB1a(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(2) = dB1a(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(3) = dB1a(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
c
	 dCb7(1) = dB1a(1)*aasum*exp2 + B1*daasum*dP(1)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(2) = dB1a(2)*aasum*exp2 + B1*daasum*dP(2)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(3) = dB1a(3)*aasum*exp2 + B1*daasum*dP(3)*exp2
     .                               + B1* aasum*dexp2
	 dCb8(1) = dB1a(1)*ffsum*exp7  +B1*dffsum*dP(1)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(2) = dB1a(2)*ffsum*exp7  +B1*dffsum*dP(2)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(3) = dB1a(3)*ffsum*exp7  +B1*dffsum*dP(3)*exp7 
     .                                +B1* ffsum*dexp7 
           dCbnds(1,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1)
     .                  +dCb6(1) +dCb7(1) +dCb8(1) 
           dCbnds(1,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2)
     .                  +dCb6(2) +dCb7(2) +dCb8(2) 
           dCbnds(1,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3)
     .                  +dCb6(3) +dCb7(3) +dCb8(3) 
c
c ----- calculate the Cbndb correction term ----------------
       B1    = B1b
       B12   = B1*B1
       B13   = B12*B1
       B14   = B13*B1
       B15   = B14*B1
       asum  = d11 + d12*R + d13*Rsq + d14/R + d15/Rsq 
       bsum  = d21*B12 + d22*B13 + d23*B14 + d24*B15
       csum  = d31*B1 *exp1  + d32*B12*exp2  
       dsum1 = d41*exp1 + d42*exp2
       dsum2 = d43*exp1 + d44*exp2
       fsum  = dx1 + d52*R + d53*Rsq 
       gsum  = d61 + d62/R + d63/Rsq 
       aasum = d71 + d72*P + d73*Psq + d74/P + d75/Psq 
       ffsum =       d81*P + d82*Psq + d84/Psq 
       dasum = d12 + 2.d0*d13*R -d14/Rsq -2.d0*d15/Rcu
       dbsum = 2.d0*d21*B1 +3.d0*d22*B12 +4.d0*d23*B13 +5.d0*d24*B14
       ddsum1= d41*dexp1 + d42*dexp2
       ddsum2= d43*dexp1 + d44*dexp2
       dfsum = d52 + 2.d0*d53*R
       dgsum = -d62/Rsq - 2.d0*d63/Rcu
       daasum = d72 +2.d0*d73*P -d74/Psq -2.d0*d75/Pcu
       dffsum = d81 + 2.d0*d82*P -2.d0*d84/Pcu
         cb(2,1)  = B1*asum * exp1  /P
         cb(2,2)  = bsum * exp2  
         cb(2,3)  = B2*csum 
         cb(2,4)  = B1*B3*dsum1 + B1*B3b*dsum2
         cb(2,5)  = B1 * fsum * exp7 /P
         cb(2,6)  = B1 * gsum * exp7  
         cb(2,7)  = B1 * aasum * exp2 
         cb(2,8)  = B1 * ffsum * exp7  
         Cbnds(2) = cb(2,1)+cb(2,2)+cb(2,3)+cb(2,4)+cb(2,5)
     .                     +cb(2,6)+cb(2,7)+cb(2,8)
c       calculate the derivatives:
	 dCb1(1) = dB1b(1)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(1)/Psq
	 dCb1(2) = dB1b(2)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(2)/Psq
	 dCb1(3) = dB1b(3)*asum*exp1/P + B1*dasum*exp1/P
     .                +B1*asum*dexp1/P - B1*asum*exp1*dP(3)/Psq
c
	 dCb2(1) = dbsum*dB1b(1)*exp2 + bsum*dexp2
	 dCb2(2) = dbsum*dB1b(2)*exp2 + bsum*dexp2
	 dCb2(3) = dbsum*dB1b(3)*exp2 + bsum*dexp2
c
	 dCb3(1)= dB2(1)*csum + B2*( d31 *dB1b(1)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(1)*exp2 + d32 *B12*dexp2 )
	 dCb3(2)= dB2(2)*csum + B2*( d31 *dB1b(2)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(2)*exp2 + d32 *B12*dexp2 )
	 dCb3(3)= dB2(3)*csum + B2*( d31 *dB1b(3)*exp1 +d31 *B1*dexp1
     .           +d32 *2.d0*B1*dB1b(3)*exp2 + d32 *B12*dexp2 )
c          dCb4 equations may 27/90:
	 dCb4(1) = dB1b(1)*B3 *dsum1 + B1*dB3(1)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(1)*B3b*dsum2 + B1*dB3(4)*dsum2 + B1*B3b*ddsum2
	 dCb4(2) = dB1b(2)*B3 *dsum1 + B1*dB3(2)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(2)*B3b*dsum2 + B1*dB3(5)*dsum2 + B1*B3b*ddsum2
	 dCb4(3) = dB1b(3)*B3 *dsum1 + B1*dB3(3)*dsum1 + B1*B3 *ddsum1
     .           + dB1b(3)*B3b*dsum2 + B1*dB3(6)*dsum2 + B1*B3b*ddsum2
c          dCb5 equations may 27/90:
	 dCb5(1) = dB1b(1)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(1)/Psq
	 dCb5(2) = dB1b(2)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(2)/Psq
	 dCb5(3) = dB1b(3)*fsum*exp7/P + B1*dfsum*exp7/P
     .               + B1*fsum*dexp7/P - B1*fsum*exp7*dP(3)/Psq
c
	 dCb6(1) = dB1b(1)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(2) = dB1b(2)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
	 dCb6(3) = dB1b(3)*gsum*exp7 + B1*dgsum*exp7 + B1*gsum*dexp7
c
	 dCb7(1) = dB1b(1)*aasum*exp2 + B1*daasum*dP(1)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(2) = dB1b(2)*aasum*exp2 + B1*daasum*dP(2)*exp2
     .                               + B1* aasum*dexp2
	 dCb7(3) = dB1b(3)*aasum*exp2 + B1*daasum*dP(3)*exp2
     .                               + B1* aasum*dexp2
	 dCb8(1) = dB1b(1)*ffsum*exp7  +B1*dffsum*dP(1)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(2) = dB1b(2)*ffsum*exp7  +B1*dffsum*dP(2)*exp7 
     .                                +B1* ffsum*dexp7 
	 dCb8(3) = dB1b(3)*ffsum*exp7  +B1*dffsum*dP(3)*exp7 
     .                                +B1* ffsum*dexp7 
           dCbnds(2,1) = dCb1(1) +dCb2(1) +dCb3(1) +dCb4(1) +dCb5(1)
     .                  +dCb6(1) +dCb7(1) +dCb8(1) 
           dCbnds(2,2) = dCb1(2) +dCb2(2) +dCb3(2) +dCb4(2) +dCb5(2)
     .                  +dCb6(2) +dCb7(2) +dCb8(2) 
           dCbnds(2,3) = dCb1(3) +dCb2(3) +dCb3(3) +dCb4(3) +dCb5(3)
     .                  +dCb6(3) +dCb7(3) +dCb8(3) 
c
c  calculate the total derivative from the pieces:
         dCbnda(1) = dT(1) * Cbnds(1)    + sumT * dCbnds(1,1)
         dCbndb(1) = dT(1) * Cbnds(2)    + sumT * dCbnds(2,1) 
         dCbnda(2) = dT(2) * Cbnds(1)    + sumT * dCbnds(1,2)
         dCbndb(2) = dT(2) * Cbnds(2)    + sumT * dCbnds(2,2) 
         dCbnda(3) = dT(3) * Cbnds(1)    + sumT * dCbnds(1,3)
         dCbndb(3) = dT(3) * Cbnds(2)    + sumT * dCbnds(2,3) 
      Cbnda = sumT*Cbnds(1)
      Cbndb = sumT*Cbnds(2)
c  some debugging print statements:
d     if(ipr.gt.1)then
d        write(7,7650) B2,B3
d        write(7,7700) db2,db3
d        write(7,*)
d        write(7,*) 'Vbend values :'
d        k=1
d        write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbnda
d        write(7,7155) (vb(1,i),i=1,5)
d        k=2
d        write(7,7150) k,dvb1,dvb2,dvb3,dvb4,dvb5,dVbndb
d        write(7,7155) (vb(2,i),i=1,5)
d        write(7,7100) Vbnda,Vbndb
d        write(7,*)
d        write(7,*) 'Cbend values:'
d        write(7,*) 'c81,c81,c82=',c81,c81,c82
d        write(7,*) 'c91,c83,c84=',c91,c83,c84
d        write(7,*) 'd81,d81,d82=',d81,d81,d82
d        write(7,*) 'd91,d83,d84=',d91,d83,d84
d        k=1              
d           write(7,2020) k,(cb(1,j),j=1,9)
d           write(7,2033) sumT,Cbnds(1)
d        k=2              
d           write(7,2020) k,(cb(k,j),j=1,9)
d           write(7,2033) sumT,Cbnds(2)
d        k=1
d        write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6,
d    .               dcb7,dcb8,dcb9,(dCbnds(1,i),i=1,3)
d        k=2
d        write(7,7152) k,dcb1,dcb2,dcb3,dcb4,dcb5,dcb6,
d    .                 dcb7,dcb8,dcb9,(dCbnds(k,i),i=1,3)
d        write(7,7200) dT,dCbnda,dCbndb
d        write(7,*) 'VbCb  exit -------------------------'
d        write(7,*)
d     end if
      return
 2000 format(5x,'r1, r2, r3 = ',3(1x,f12.8))
 2010 format(5x,'B1a, B1b, B2, B3 = ',4(1x,g12.6))
 7100 format('    Vbnda = ',g12.6,'     Vbndb = ',g12.6)
 7650 format('    B2 =    ',g12.6,'       B3 =  ',g12.6)
 7700 format('    dB2 =   ',3(1x,g12.6),/,
     .       '    dB3 =   ',6(1x,g12.6))
 7150 format(' k=',i1,'     derivatives: ',/,
     .       ' dVb1  = ',3(1x,g16.10),/,
     .       ' dVb2  = ',3(1x,g16.10),/,
     .       ' dVb3  = ',3(1x,g16.10),/,
     .       ' dVb4  = ',3(1x,g16.10),/,
     .       ' dVb5  = ',3(1x,g16.10),/,
     .       ' dVbnd = ',3(1x,g16.10))
 7155 format('  Vb1 = ',g16.10,'  Vb2 = ',g16.10,/,
     .       '  Vb3 = ',g16.10,'  Vb4 = ',g16.10,/,
     .       '  Vb5 = ',g16.10)
 2020 format(5x,'k=',i1,' cb1,cb2,cb3= ',3(1x,g16.10),
     .     /,8x,' cb4,cb5,cb6= ',3(1x,g16.10),
     .     /,8x,' cb7,cb8,cb9= ',3(1x,g16.10))
 2033 format(5x,'sumT = ',g12.6 ,'  cbnds(k) = ',g12.6)
 7152 format('Subr.Cbend:   k=',i1,'     derivatives: ',/,
     .       '  dCb1  = ',3(1x,g18.12),/,
     .       '  dCb2  = ',3(1x,g18.12),/,
     .       '  dCb3  = ',3(1x,g18.12),/,
     .       '  dCb4  = ',3(1x,g18.12),/,
     .       '  dCb5  = ',3(1x,g18.12),/,
     .       '  dCb6  = ',3(1x,g18.12),/,
     .       '  dCb7  = ',3(1x,g18.12),/,
     .       '  dCb8  = ',3(1x,g18.12),/,
     .       '  dCb9  = ',3(1x,g18.12),/,
     .       '  dCbnd = ',3(1x,g18.12))
 7200 format('  dT    = ',3(1x,g18.12),/,
     .       '  dCbnda= ',3(1x,g18.12),/,
     .       '  dCbndb= ',3(1x,g18.12))
 6000 format(' from subr.b1ab:  r=',3(1x,f9.6))
 6010 format('    b1a = ',f16.10,'     b1b = ',f16.10)
 6020 format('    db1a: ',3(1x,e12.6),/,'    db1b: ',3(1x,e12.6))
 7000 format('    dc(j)/dr(i) derivatives: ',/,
     .    8x,3(1x,e12.6),/,8x,3(1x,e12.6),/,8x,3(1x,e12.6))
 7010 format('    ri2-rj2-rk2: ',3(1x,e12.6),/,
     .       '       cosines : ',3(1x,e12.6),/,
     .       '    12c**2 - 3 : ',3(1x,e12.6))
      end
c


