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