      module ccidata
      save
      double precision ::  eshift= 0.1744759989649372d0
      double precision ::  beta=0.9246250191375499d0    
      double precision ::  rnought=2.d0                         
      double precision, dimension(89) ::  xlin1=(/-0.1682626669809615d0,        
     &     1.697952225005269d0,      3.959367441888326d0,   
     &    -7.159854944676235d0,     -2.462488887191634d0,  
     &     7.998434264390882d0,     -3.665712084901418d0, 
     &    -0.2155463656483861d0,     0.5620418518154215d0, 
     &    -0.1369812453279634d0,     1.0063614579308723d-2,
     &    41.07358514471853d0,      67.02217912830267d0, 
     &  -100.0014758720550d0,        2.481875721304641d0, 
     &    32.46244740624601d0,     -15.32460432587821d0, 
     &     3.814631583515083d0,     -0.3155985331636701d0, 
     &    -4.2198516063816239d-3, -157.7203098415152d0, 
     &    17.11141015275416d0,     -61.48638359701442d0, 
     &     4.244841908917134d0,      5.300393723900312d0, 
     &    -0.6170558358973661d0,    -6.2143674437616274d-3,
     &     5.785766830955630d0,     -0.1133784523542150d0, 
     &    -0.2539145873400511d0,    -4.3406559735048131d-2,
     &     2.4548851428702549d-2,   -1.092841685240530d0, 
     &     0.8131178021738225d0,    -7.7933784638615652d-2,
     &    -0.1484455564689179d0,    -7.351078313940695d0, 
     &   128.7480567271992d0,      199.6058760392980d0, 
     &   -50.18444648711906d0,    -178.2184990028168d0, 
     &   129.6832433574917d0,      -37.57036725437851d0, 
     &     5.747838042162251d0,      0.1434298979682676d0, 
     &    -8.3110343245055174d-2,  428.7635440340811d0, 
     &   495.9564113906064d0,    -1742.123837168904d0, 
     &  -259.0564024322437d0,      250.2352551214438d0, 
     &   -44.15232947384926d0,       0.9963955068997190d0, 
     &     0.3136096547683901d0, -1032.921631265462d0, 
     &   199.9708129167335d0,     -135.4683186684037d0, 
     &   -14.50534672808297d0,       5.650612567917722d0, 
     &    -0.1473535040913008d0,     1.351710421967908d0, 
     &    29.60972395839821d0,      -7.224575522867537d0, 
     &     0.3582134126835411d0,     0.2589923362572861d0, 
     &     0.1819402887127863d0,  1312.633253770123d0, 
     &  5189.281073844522d0,     -2727.020597875299d0, 
     &  -122.2193440638810d0,      152.4859537357477d0, 
     &    -5.562127541501408d0,     -1.549418235204051d0, 
     &  2647.710764337011d0,       502.1991025406292d0, 
     &  -152.3069488960497d0,       -5.740041181528469d0, 
     &     1.362209926084833d0,     17.96504193139269d0, 
     &    10.13930464010770d0,      -2.287463377365149d0, 
     &     1.932492687974671d0,  -3337.181754475553d0, 
     &   182.9671187483376d0,       20.29606065582956d0, 
     &     1.192500349837268d0,    -22.69079850990068d0, 
     &    -0.5967747962499963d0,     1.091111370927579d0/)
      double precision, dimension(89) :: xlin2=(/3.2289532152874986d-2,   
     &    -8.6686506736620789d-2,   -0.3266024449171150d0,   
     &    -0.7972798579742405d0,     2.812854025206053d0,   
     &    -2.374871328104482d0,      0.6089164699935086d0,   
     &     0.1635658322754857d0,    -0.1263080841790984d0,   
     &     2.5221027611065153d-2,   -1.6677760633386995d-3,
     &    -0.3372069966080046d0,    -8.377537389356892d0, 
     &     3.694668699246826d0,     11.36868357221171d0, 
     &   -11.88889428297186d0,       4.270226804928180d0, 
     &    -0.4526863369836274d0,    -5.2518975083660345d-2,
     &     9.0223742463199624d-3,   38.03127084604298d0, 
     &   -12.82065963469141d0,      -5.663705994784566d0, 
     &     7.300080908068341d0,     -0.5157232521089076d0, 
     &    -0.2200954472805637d0,     2.0686337667851093d-2,
     &    14.19614073502447d0,       1.726255538125508d0, 
     &    -2.526391344985463d0,      6.0598667264186397d-2,
     &     3.8560006141770319d-2,    2.306725001955670d0, 
     &    -0.3197660708712612d0,    -1.9982819587147187d-2,
     &     0.1036787525350017d0,    -2.030231243513758d0, 
     &    25.16474735933574d0,     -14.13898387128746d0, 
     &   -62.01219930765056d0,      75.56471352787402d0, 
     &   -30.25154027581496d0,       4.257663341071418d0, 
     &     0.6142213624561507d0,    -0.2458345005137628d0, 
     &     1.9135101060142316d-2,  168.2482649287792d0, 
     &  -270.3587962365308d0,       72.46467150857366d0, 
     &    32.12278113625301d0,      11.89435258275814d0, 
     &    -5.521268211450904d0,      0.8934653854937415d0, 
     &    -7.6805320500971247d-2, -422.2672942154591d0, 
     &   -35.92489169679082d0,      14.74665733902546d0, 
     &     6.020706922914545d0,      2.4538599994799737d-3,
     &    -0.1370624616096560d0,    44.47654636843973d0,
     &   -11.18166310447657d0,       3.6110526441458540d-2,
     &     0.1539535582849851d0,     2.065981296572923d0, 
     &    -0.1624439939014451d0,  -658.3723053145102d0, 
     &  -883.4947620258348d0,     -263.8479489759550d0, 
     &    44.68359598895833d0,      29.88970883430300d0, 
     &    -5.663152586850906d0,      0.2413124485342995d0, 
     &   656.9813540494517d0,      -91.56710742448996d0, 
     &   -30.38175110815412d0,       2.066013906655882d0, 
     &     0.5314632792073760d0,    54.99422880821747d0, 
     &    -0.7887249891070589d0,    -0.9164068632403487d0, 
     &     0.6576346543435602d0,   292.2665706493461d0, 
     &   -32.04152511166628d0,      -2.104943871500433d0, 
     &   -0.7251129002870176d0,      4.709122797258667d0,   
     &    1.119971732384046d0,      -4.066519840824071d0/)
      double precision, dimension(71) ::  xlin3=(/6.3998636332917941d-2,    
     &    -0.5826083150212633d0,    -3.294968443823759d0, 
     &     5.751817019955082d0,     -1.486176533347650d0, 
     &    -0.6845411439568367d0,    -0.9083627566591157d0, 
     &     1.305678536333854d0,     -0.4821780067143284d0, 
     &     5.5182580699196751d-2,  -22.57213192201081d0, 
     &   -61.37436407473290d0,      52.64511100643983d0, 
     &    29.99521825023724d0,     -38.31341936937510d0, 
     &    14.86681566870566d0,      -3.460352146284503d0, 
     &     0.1343449644795501d0,    53.21459427503272d0, 
     &    64.60699757442232d0,      30.89466488573955d0, 
     &     8.412732049211240d0,     -4.274560634850636d0, 
     &    -0.1788031871862021d0,   -46.47154313325355d0, 
     &    24.88995914944731d0,      -2.757433396373881d0, 
     &    -0.3635530438220282d0,     1.049428472748707d0, 
     &    -0.9909460260393127d0,     7.623381850055868d0, 
     &   -75.76387696233419d0,    -162.8879929453359d0, 
     &   -29.98317728892977d0,     168.7972614497797d0, 
     &  -102.0864921135087d0,       31.38363768116054d0, 
     &    -6.084136544334489d0,     -2.9400115678577254d-2,
     &  -417.4808742201730d0,     -601.3246381855070d0,
     &  1111.017825017578d0,       749.0346863822200d0, 
     &  -254.9798923160399d0,       37.02896288679148d0, 
     &    -0.7010161055223685d0,   928.6100310886904d0, 
     &   513.3554862860957d0,      138.4248255502364d0, 
     &    -0.9347938923678116d0,     3.867048618649862d0, 
     &  -109.6794304656507d0,       10.30110710850976d0, 
     &     0.5118243051520034d0,     0.2069543381912151d0, 
     &  -796.4010440692600d0,    -4693.456970734772d0, 
     &  2597.246850944603d0,       630.4845628327919d0, 
     &  -221.0729177223342d0,        9.839396652075672d0, 
     & -5227.790162633692d0,      -529.4879950032479d0, 
     &    30.53701394166368d0,     -14.22864486129778d0, 
     &   -35.14731691232372d0,       2.140914121987492d0,  
     &  2250.367423678786d0,       138.5994085512685d0,  
     &    70.30416282155014d0,     -62.74844809925895d0/) 
      double precision, dimension(71) :: xlin4=(/-1.2612782515309912d-2,
     &     1.5724446750473936d-2,     0.2653887372495336d0, 
     &     0.1216019078797281d0,     -0.6102844930670878d0, 
     &     4.6420446436398755d-3,     0.5563651571664677d0, 
     &    -0.3890198684837317d0,      0.1052569138798340d0, 
     &    -1.0080347996402828d-2,     0.2409741384095758d0, 
     &     5.745106082351191d0,       1.047295525278513d0, 
     &   -12.03507114814484d0,        9.035186991905373d0, 
     &    -2.267198360038944d0,       2.0836099694222385d-2,
     &     4.7445938282087222d-2,   -18.26537244396756d0, 
     &    -4.809068965727743d0,       5.658647992475390d0, 
     &    -1.284093508478602d0,      -2.164733987255563d0, 
     &     0.3630073042947094d0,      3.760490137844927d0, 
     &    -5.927286374087972d0,       0.1278839922444694d0, 
     &     0.5381095462111573d0,     -2.179370013756247d0, 
     &     0.4127566904703107d0,      0.7168269469196163d0, 
     &    -9.502943939736065d0,      -9.597638311349975d0, 
     &    65.60860814240448d0,      -56.04003254284753d0, 
     &    15.38423202570607d0,       -0.3167914651468724d0, 
     &    -0.7643266917502949d0,      0.1110705439245336d0, 
     &  -146.5638710976931d0,       151.8079548350792d0, 
     &    15.86486188082157d0,      -65.56160562330304d0, 
     &    -1.126379669058077d0,      -1.759058937005427d0, 
     &    -3.1582037336166588d-4,   410.9371301349821d0, 
     &   142.2842073695382d0,       -15.14809406881419d0, 
     &    -9.858662824266736d0,       0.2622049408092991d0, 
     &   -72.67984057864001d0,       11.28350883179339d0, 
     &    -0.5323962789000201d0,     -2.447540583751350d0, 
     &   148.9575982135905d0,      1249.939675860206d0, 
     &   268.6478693121586d0,        33.67689052007242d0, 
     &   -35.65384689269180d0,        3.464021415311834d0, 
     &  -103.4190791219546d0,        76.76961997409407d0, 
     &    11.59182442022548d0,       -0.5254501536114633d0, 
     &   -32.07157614897808d0,        1.102243506822152d0, 
     &  -952.9757997893194d0,        52.11461967121677d0, 
     &     6.929202177070909d0,      -6.322557284152793d0/)
      end module
      subroutine prepot
      implicit none
      integer          :: ivp,nt
      double precision :: rvp(nt,3),evp(nt),r1x,r2x,r3x
      call prepotcbs
      return

      entry pot(rvp, evp, nt)
c
      do ivp = 1, nt
         r1x=rvp(ivp,1)  
         r2x=rvp(ivp,2)  
         r3x=rvp(ivp,3)  
         call potcbs(r1x,r2x,r3x,evp(ivp))
      enddo    
c
      return
      end

      subroutine prepotcbs
      use ccidata
      implicit double precision (a-h,o-z)
      save
      parameter (n2=5000)
      common /indy/iind(n2),jind(n2),kind(n2),ipar(n2),maxi,maxi2
      double precision :: rho1t(0:12),rho2t(0:12),rho3t(0:12)

      write(6,*)' Using CCI H3 PES of 8/15/01'
      write(6,*)' S. L. Mielke, B. C. Garrett, and K. A. Peterson, J. Ch
     &em. Phys. 116, 4142 (2002).'

      epslon=1.d-12
      epscem=1.d-12
      ald=0.02d0
      betad=0.72d0

      maxi2=0
      maxi=0
      call indexa3(12,12)     
      maxi2a=maxi2

      maxi=0
      call indexa3(11,11)        
      maxi2b=maxi2-maxi2a
      return

      entry potcbs(r1x,r2x,r3x,eee)
      call trip(r1x,etrip1)
      call trip(r2x,etrip2)
      call trip(r3x,etrip3)
      call singlet(r1x,es1)
      call singlet(r2x,es2)
      call singlet(r3x,es3)
      q1=0.5d0*(es1+etrip1)
      q2=0.5d0*(es2+etrip2)
      q3=0.5d0*(es3+etrip3)
      xj1=0.5d0*(es1-etrip1)
      xj2=0.5d0*(es2-etrip2)
      xj3=0.5d0*(es3-etrip3)
      xj=sqrt(epslon+0.5d0*((xj3-xj1)**2+(xj2-xj1)**2+(xj3-xj2)**2))    
      xjs=0.25d0*(xj1+xj2+xj3)
      vlondon=q1+q2+q3-xj

      rho1=exp(-beta*(r1x-rnought))
      rho2=exp(-beta*(r2x-rnought))
      rho3=exp(-beta*(r3x-rnought))

      if((r1x.le.betad).or.(r2x.le.betad).or.(r3x.le.betad))then
        damper=0.d0
      else
        rdamp1=ald/(betad-r1x)
        rdamp2=ald/(betad-r2x)
        rdamp3=ald/(betad-r3x)
        damper=exp(rdamp1+rdamp2+rdamp3)
      endif

      bb=sqrt(epscem+(r1x-r3x)**2+(r1x-r2x)**2+(r2x-r3x)**2)
      warp=1.d0/(1.d0/r1x+1.d0/r2x+1.d0/r3x)

      v=eshift+vlondon

      sum1=0.d0
      sum2=0.d0
      sum3=0.d0
      sum4=0.d0
      
      rho1t(0)=1.d0
      rho2t(0)=1.d0
      rho3t(0)=1.d0

      do it=1,12
        rho1t(it)=rho1t(it-1)*rho1
        rho2t(it)=rho2t(it-1)*rho2
        rho3t(it)=rho3t(it-1)*rho3
      enddo    

      do ii=1,maxi2a
        rprod=rho1t(kind(ii))*rho2t(jind(ii))*rho3t(iind(ii))
        sum1=sum1+xlin1(ipar(ii))*rprod
        sum2=sum2+xlin2(ipar(ii))*rprod
      enddo     

      do ii=maxi2a+1,maxi2a+maxi2b
        rprod=rho1t(kind(ii))*rho2t(jind(ii))*rho3t(iind(ii))
        sum3=sum3+xlin3(ipar(ii))*rprod
        sum4=sum4+xlin4(ipar(ii))*rprod
      enddo     

      eee=v+damper*(sum1+sum2*bb+sum3*warp+sum4*bb*warp)

      return
      end
      
      subroutine indexa3(mbig,mtop)            
      parameter (n2=5000)
c     calculate the fitting indices for A3 symmetry
      common /indy/iind(n2),jind(n2),kind(n2),ipar(n2),maxi,maxi2

      l=maxi2
      l2=maxi

      do  i=0,min(mtop,mbig)
       do  j=i,min(mtop,mbig)
        do  k=j,min(mtop,mbig)
         isum=i+j+k
          if(isum.le.mbig)then
           if(((isum-i).gt.0).and.((isum-j).gt.0).and.
     &        ((isum-k).gt.0))then     
            l2=l2+1
            l=l+1
            ipar(l)=l2
            iind(l)=i
            jind(l)=j
            kind(l)=k
            if(i.ne.j)then
             if(j.ne.k)then
               l=l+1
               ipar(l)=l2
               iind(l)=i
               jind(l)=k
               kind(l)=j
               l=l+1
               ipar(l)=l2
               iind(l)=j
               jind(l)=i
               kind(l)=k
               l=l+1
               ipar(l)=l2
               iind(l)=j
               jind(l)=k
               kind(l)=i
               l=l+1
               ipar(l)=l2
               iind(l)=k
               jind(l)=i
               kind(l)=j
               l=l+1
               ipar(l)=l2
               iind(l)=k
               jind(l)=j
               kind(l)=i
              else
               l=l+1
               ipar(l)=l2
               iind(l)=j
               jind(l)=i
               kind(l)=k
               l=l+1
               ipar(l)=l2
               iind(l)=j
               jind(l)=k
               kind(l)=i
             endif
            elseif (j.ne.k)then
               l=l+1
               ipar(l)=l2
               iind(l)=j
               jind(l)=k
               kind(l)=i
               l=l+1
               ipar(l)=l2
               iind(l)=k
               jind(l)=j
               kind(l)=i
            endif
           endif
          endif
        enddo     
       enddo    
      enddo    
      maxi=l2
      maxi2=l
      return
      end

      subroutine trip(r,v)
c     H2 triplet curve for FCI/CBS       
      implicit none
      integer          :: j
      double precision :: r,v,damp,xd,prefac
      double precision :: c6=-6.499027d0
      double precision :: c8=-124.3991d0
      double precision :: c10=-3285.828d0
      double precision :: beta=0.2d0
      double precision :: re=1.401d0

      double precision :: alpha=4.342902497495857d0    
      double precision,dimension(17) :: xlp=(/2.190254292077946d-1,
     &           6.903970886899540d-1, 1.163805164813647d+0,
     &           1.337274224278977d+0, 1.176131501286896d+0,
     &           9.093025316067311d-1, 5.908490740317524d-1,
     &           1.467271248985162d-1, 6.218667627370255d-2,
     &           1.105747186790804d-1, -7.112594221275334d-2,
     &          -8.866929977503379d-3, 3.823523348718375d-2,
     &          -2.055915301253834d-2, 5.419433399532862d-3,
     &          -7.263483303294653d-4, 4.154351972583585d-5/)
      damp=1.d0-exp(-beta*r**2)
      xd=damp/r  
      v=c6*xd**6+c8*xd**8+c10*xd**10

      prefac=exp(- alpha*(r-re))
      if((r-re).ne.0.d0)then
        do j=1,17
          v=v+xlp(j)*(r-re)**(j-1)*prefac
        enddo    
      else
        v=v+xlp(1)*prefac
      endif

      return
      end

      subroutine singlet(r,v)
c     H2 singlet potential for FCI/CBS
      implicit none
      integer          :: j
      double precision :: r,v,damp,xd,prefac
      double precision :: c6=-6.499027d0
      double precision :: c8=-124.3991d0
      double precision :: c10=-3285.828d0
      double precision :: beta=0.2d0
      double precision :: re=1.401d0
      double precision :: alpha=3.980917850296971d0    
      double precision, dimension(17) :: xlp= (/ -1.709663318869429d-1,
     &          -6.675286769482015d-1, -1.101876055072129d+0, 
     &          -1.106460095658282d+0, -7.414724284525231d-1,
     &          -3.487882731923523d-1, -1.276736255756598d-1,
     &          -5.875965709151867d-2, -4.030128017933840d-2,
     &          -2.038653237221007d-2, -7.639198558462706d-4,
     &           2.912954920483885d-3, -2.628273116815280d-4,
     &          -4.622088855684211d-4, 1.278468948126147d-4,
     &          -1.157434070240206d-5, -2.609691840882097d-12/)

      damp=1.d0-exp(-beta*r**2)
      xd=damp/r 
      prefac=exp(- alpha*(r-re))
      v=c6*xd**6+c8*xd**8+c10*xd**10

      if((r-re).ne.0.d0)then
        do j=1,17
          v=v+xlp(j)*(r-re)**(j-1)*prefac
        enddo    
      else
        v=v+xlp(1)*prefac
      endif

      return
      end
