module a2data save double precision :: eshift=0.1662307680057728d0 double precision :: beta=0.900d0 double precision :: rnought=2.d0 double precision, dimension(89) :: xlin1=(/0.3122096268380428d0, & -0.5091675305404870d0, 17.09937094766591d0, & -47.55205835431389d0, 58.80549483490758d0, & -38.83145155729753d0, 11.88747689358562d0, & 0.2617252059369060d0, -1.142054614272080d0, & 0.2728967965513714d0, -1.9323083032991791d-2, & 57.04567478144305d0, 148.5636591372696d0, & -318.6179220492799d0, 182.3742028409573d0, & -144.0665991614994d0, 62.04257754042266d0, & -5.877331776092552d0, 0.4288520217861720d0, & -0.1036073697125561d0, 120.0475914921467d0, & 263.6164106497357d0, -462.6236113095622d0, & 116.0608660637904d0, -19.91101916036278d0, & 1.193693981250914d0, 0.1422177054289994d0, & -160.3114246747076d0, 89.72253917752164d0, & -9.071109638571208d0, -2.0135523158362437d-2, & 0.1243557501187861d0, -12.96200461060273d0, & -2.347887744857734d0, 0.4369816290614699d0, & 0.1861905441025299d0, -8.284586975531951d0, & -20.42241348042063d0, 958.1242541736171d0, &-1928.583970392355d0, 1119.594819408116d0, & -696.0896180464565d0, 272.2524874095656d0, & -17.44951726252673d0, -1.588843767916973d0, & -8.3560790211881619d-2, 2386.104049771720d0, & 3941.711173912721d0, -5088.409610323098d0, & 2147.437212950090d0, -1209.070332780972d0, & 228.5650959646700d0, -18.75623216539055d0, & 1.331518943744410d0, 1323.582901230102d0, & 1366.507728689508d0, -470.6376874755080d0, & -4.840542534010448d0, 4.418210908054986d0, & 0.2615031877550023d0, 365.2205662123347d0, & 56.60696871480308d0, 20.20287252717788d0, & -8.939307433105093d0, -7.579362430843242d0, & 1.945102121992439d0, 23893.31732811452d0, &14785.75023538304d0, -5590.061708848788d0, & 2659.366722765638d0, -1103.540790332176d0, & 133.2401706836243d0, -2.779352658653005d0, &-2487.131159960726d0, 223.7986868185370d0, & 940.7101507094950d0, -119.8581541479964d0, & -5.864853168420393d0, -1471.392587738388d0, & 20.46733470843232d0, 44.81137223542368d0, & -37.13875660767655d0, -4114.467297104680d0, & 435.2902050180833d0, 110.1713440583027d0, & -6.086172632709467d0, -58.60516318227123d0, & -30.37217546242515d0, 114.1717244633075d0/) double precision, dimension(89) :: xlin2=(/-2.3691917237142486d-2, & 0.6433593796852592d0, -3.191080817325861d0, & 5.082622923008778d0, -4.303785253384310d0, & 2.677917664616706d0, -1.482618068308200d0, & 0.6784912558836848d0, -0.2049921467505755d0, & 3.1560425722410100d-2, -1.9786106444049946d-3, & -1.510131880107122d0, -28.60321443064080d0, & 53.22691305234915d0, -40.11411939395515d0, & 20.81861842742526d0, -4.447977855272196d0, & 0.2301695113250538d0, -0.1295406658545045d0, & 2.4496342125756746d-2, 23.93815251101609d0, & -15.97857074145742d0, -19.75274167695325d0, & 15.86955351115030d0, 0.4841116429447808d0, & -1.3958200902611101d-2, -6.8392975615942639d-2, & -6.396616748692148d0, -0.5798087307661497d0, & 1.564469079111058d0, -9.8289900209482819d-2, & -4.9186020637759537d-2, 16.99519683721401d0, & -2.466567635268690d0, 1.0845662786683728d-2, & 0.2941408959952037d0, 1.777059184168285d0, & 46.60356704672536d0, -149.4923339330999d0, & 153.0203582840253d0, -129.5745301688355d0, & 83.01390198308063d0, -20.50411666044072d0, & 1.329739978962026d0, 8.1037167237220931d-2, & -5.6360423765908704d-3, -224.7409562531643d0, & -445.7702964503590d0, -182.9799944501726d0, & -133.0432327364785d0, 227.3077678979116d0, & -34.61637061865006d0, -8.9508185399590101d-2, & 7.1955361459616729d-2, 442.7525395172289d0, & 2.458625503871524d0, -202.8428507849074d0, & 32.31023899069420d0, -2.282723397214718d0, & 0.3055028229001248d0, 181.8487130166304d0, & -69.54292969634861d0, -6.889401006201262d0, & 2.218301674443548d0, 32.61450081881122d0, & -1.829579173006778d0, -3053.146122439747d0, & 2039.754739773428d0, -1925.218456128402d0, & 149.3339608613740d0, 28.43501684561451d0, & 0.3254831648262675d0, -0.3970879353555506d0, & 3299.500750582070d0, -594.7502330927563d0, & -65.96404318767495d0, 28.40550550560758d0, & -1.865692389204451d0, 448.8499645057911d0, & -30.33721822994903d0, -3.221780195307217d0, & 3.085642390961864d0, 498.9201490733377d0, & -151.2816443510426d0, -76.13951656899020d0, & 3.665378992703659d0, 72.17259482865131d0, & 6.152241900437081d0, -31.61340871373989d0/) double precision, dimension(71) :: xlin3=(/-0.1344994646896329d0, & -3.2141410415166161d-2, -7.517359288762402d0, & 18.44051637133763d0, -14.31587410374696d0, & -2.834545980827915d0, 12.12504845332679d0, & -8.087779341389851d0, 2.266439914858519d0, & -0.2458975078101396d0, -24.78954125497926d0, & -130.9696482818210d0, 208.5056810553900d0, & -102.7798259674824d0, 94.50458094288865d0, & -26.66242091873472d0, -3.796095988884443d0, & -0.2135440885386056d0, -101.9329951911512d0, & -237.0984558426526d0, 340.3675666528260d0, & -5.600055272704649d0, -6.626620457234992d0, & 2.250749424450074d0, 72.61834696241013d0, & 28.58895062725391d0, -25.96494974227168d0, & 1.775974157654598d0, -10.06825158236352d0, & 7.123101687824712d0, 2.773473562095701d0, & 34.25830296691608d0, -613.3162919232790d0, & 1003.318114503314d0, -293.8447014270225d0, & 238.0037203719179d0, -45.54821241600535d0, & -38.15460316807260d0, 2.152193113454373d0, & -1265.105998633106d0, -4182.774652447684d0, & 4323.977859927758d0, -1511.567227580535d0, & 1013.128818248927d0, -145.0171911214698d0, & 17.32391332654659d0, -2672.569843994831d0, & -494.6368772637625d0, -176.6910389060321d0, & 333.1709756767322d0, -44.43847568245701d0, & -715.3808858788194d0, -86.17184561227721d0, & -19.32116072356918d0, -16.23431382053576d0, & -17732.22229184206d0, -21836.37425726624d0, & 6027.463948377972d0, -3413.560409537637d0, & 1462.120516808686d0, -132.4258810442149d0, & -5461.158978972337d0, 1036.230524299047d0, & -1687.326937421790d0, 136.8271126588730d0, & 1563.706329850235d0, 55.53727444282163d0, & 6546.413397438990d0, 273.3238154461579d0, & -123.8048659677774d0, -69.99455204122073d0/) double precision, dimension(71) :: xlin4=(/1.2472601237000719d-2, & -0.2563237394846267d0, 1.141626638949766d0, & -1.268789070877211d0, 0.2669120728039442d0, & 0.4795337070923109d0, -0.4870693693992560d0, & 0.2180413266051965d0, -5.4940116209431157d-2, & 8.1898224212963791d-3, -0.9232881113457609d0, & 22.37373927583910d0, -31.43191815577488d0, & 15.83511958871624d0, -4.774566352751313d0, & -1.809376059889746d0, 0.8755187623690739d0, & 3.5850954944800173d-2, -7.088532812068095d0, & 5.757713051107022d0, 10.86426353414864d0, & -5.823082328751028d0, -3.528581474490192d0, & -0.1229368602139166d0, -10.79753576660658d0, & 24.25244965677145d0, -8.953867875122164d0, & 0.4615248048240077d0, -23.16730300552956d0, & 0.3228037267248082d0, -0.4960917683113312d0, & -29.27293769041133d0, 77.13556531104766d0, & -65.28389682361755d0, 48.11222198276875d0, & -29.05873944480265d0, -0.1063025548510765d0, & 2.448839826203129d0, -0.2802338471724751d0, & 51.92481774609280d0, 530.2135642610900d0, & -9.028255840313982d0, 346.5780338983027d0, & -240.1913305251466d0, 8.295317216645913d0, & -0.2022074738812011d0, -181.7378257010903d0, & 56.03276885734908d0, 155.3619996347792d0, & -45.12505868315263d0, 10.96008703215353d0, & -242.7878951933611d0, 101.8270659259155d0, & -0.4864797895058910d0, -7.874126046335711d0, & 2460.968493617616d0, -65.55982339327284d0, & 1112.649677115970d0, 300.2383741022979d0, & -150.4025249750507d0, 23.70714634112154d0, & -4536.070861321965d0, 451.5619825013044d0, & 167.3937904161444d0, -51.79228062736297d0, & -288.7604569805079d0, 10.58119538655392d0, & -2266.622734963641d0, 164.3048434953309d0, & 102.2707553778119d0, -71.61259882125256d0/) end module subroutine prepot implicit none integer :: ivp,nt double precision :: rvp(nt,3),evp(nt),r1x,r2x,r3x call prepota2 return entry pot(rvp, evp, nt) c do ivp = 1, nt r1x=rvp(ivp,1) r2x=rvp(ivp,2) r3x=rvp(ivp,3) call pota2(r1x,r2x,r3x,evp(ivp)) enddo c return end subroutine prepota2 use a2data 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 A2 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 pota2(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/aug-cc-pVDZ 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.9139661662872520d0 double precision,dimension(17) :: xlp=(/2.2349220271817220d-01, & 6.0356997147268380d-01, 8.9383829921944970d-01, & 9.0317491929985030d-01, 7.2775729411232600d-01, & 4.9697572693224850d-01, 2.5866961299917790d-01, & 6.8824590522579050d-02, 4.8273710815396330d-03, & 2.8416581979833760d-03, -5.7032138025732140d-03, & 2.3004916075318040d-03, 3.1740661346757730d-03, & -2.3028140070272670d-03, 6.4778830165694920d-04, & -8.5334197479958010d-05, 4.5529457450819180d-06/) 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/aug-cc-pVDZ 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.0578706387229880d0 double precision, dimension(17) :: xlp= (/-1.6244445106378250d-1, & -6.6081156795821980d-01, -1.1236729331029230d+00, & -1.1564282634392240d+00, -7.9003598933823580d-01, & -3.7515947232144770d-01, -1.3997385158833170d-01, & -6.5955219641305790d-02, -4.0788777170244650d-02, & -1.8988748368877420d-02, -2.2257283819080430d-03, & 2.5357629098565200d-03, 1.5612692551680850d-05, & -6.9340042488041930d-04, 2.3175123448115890d-04, & -3.3226980442299310d-05, 1.6041565480651940d-06/) 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