cccc------------------------------------------------ SUBROUTINE h3ps(x,vcor) IMPLICIT REAL * 8 (A-H,O-Z) DIMENSION x(3) r1=x(1) r2=x(2) r3=x(3) CALL h3pot(r1,r2,r3,vcor) RETURN END ccc------------------------------------------------------------------- subroutine h3pot(r1,r2,r3,vcor) c J. Chem. Phys. 121, 9343(2004). C POTENTIAL ENERGY SURFACE FOR H_3^-. r1,r2,r3 represent C the bond distances in au for Ha-Hb, Hb-Hc and Hc-Ha respectively. C vcor is the energy including the long-range correction in au. implicit real * 8 (a-h,o-z) dimension der(3) smlr=r2 call h3pot_ANP(r1,r3,r2,e123,e12,e13,e23,der) c LONG RANGE POTENTIAL CALCULATION c quadrupole moment data q1,q2,q3,q4,q5,q6,q7/-52.0714,21.6823,-3.35413,2.09243, 1 60.5523,-32.8347,8.11269/ quad=(q7*smlr**5+q6*smlr**4+q5*smlr**3+q1*smlr**2 1 +q2*smlr+q3)*exp(-q4*smlr) c alpha0 data p1,p2,p3,p4,p5,p6/37.9853,-15.3491,4.52484, 1 1.25065,-25.5122,11.0013/ alpha01=(p6*smlr**4+p5*smlr**3+p1*smlr**2 1 +p2*smlr+p3)*exp(-p4*smlr) alpha02=9.14+2.53*exp(-2.6*(smlr-4.0)) if(smlr.le.4.0)then alpha0=alpha2 else alpha0=alpha1 endif c alpha2 data s1,s2/9.22521,0.5022/ alpha21=s1*exp(-s2*(smlr-3.2d0)*(smlr-3.2d0)) alpha2=alpha21*0.67 c conversion from bond-angle to jacobi coordinates capr=DSQRT(0.50D0*R1**2+0.50D0*R3**2+ & -0.25D0*R2**2) cgamma=(R3**2-R1**2)/R2/capr/2.0D0 plg=((3*cgamma*cgamma)-1)/2 vlr1=-0.5*(alpha0)/capr**4-0.5*alpha2*plg/capr**4 vlr2=-quad*plg/capr**3 vlr=vlr1+vlr2 xx=6.6d0/capr if(xx.lt.1)then swi=exp(1.0d0)*exp(-1.0d0/(1.0d0-xx)) etot=e123*(1.0d0-swi)+swi*vlr elseif(xx.gt.1.0d0)then swi=0.0d0 etot=e123 else endif vcor=etot+e12+e13+e23+0.17300816d0 return end subroutine h3pot_ANP(r12,r13,r23,e123,e12,e13,e23,der) implicit real * 8 (a-h,o-z) dimension der(3) call diat12(r12,e12,d12) call diat12(r13,e13,d13) call diat12(r23,e23,d23) call triaaa(r12,r13,r23,e123,der) der(1)=d12+der(1) der(2)=d13+der(2) der(3)=d23+der(3) return end ************************************************************************ subroutine diat12(r,ener,der) ************************************************************************ * This subroutine computes the energies of a diatomic potential * fitted to 17 points * rms = 0.00063876 kcal/mol * emax = 0.00153520 kcal/mol ************************************************************************ implicit real*8 (a-h,o-z) dimension cf( 10) data cf( 1)/0.613866431409E+00/ data cf( 2)/-.190101798306E+01/ data cf( 3)/-.125228951225E+02/ data cf( 4)/0.462035129873E+03/ data cf( 5)/-.719421334160E+04/ data cf( 6)/0.675855488723E+05/ data cf( 7)/-.402534224315E+06/ data cf( 8)/0.148601194147E+07/ data cf( 9)/-.310617990248E+07/ data cf( 10)/0.281672666657E+07/ e0= 0.0000000E+00 der=0.d0 vex1= 0.1499133E+01 vex2= 0.1000000E+01 aux = 1.d0/r bux = dexp(-vex2*r)*aux cux = dexp(-vex1*r) ener=e0+cf(1)*bux dux=1.d0 eux=r*cux do 1 i=2, 10 der=der+(i-1)*cf(i)*dux dux=dux*eux ener=ener+cf(i)*dux 1 continue der=der*(1.d0-vex1*r)*cux der=der-cf(1)*(vex2+aux)*bux return end ************************************************************* subroutine triaaa(r12,r13,r23,ener,der) ************************************************************* * This subroutine computes the energies of a 3D PES * for the AAA system class * rms = 0.10481571 kcal/mol * emax = 0.42928555 kcal/mol ************************************************************* implicit real*8(a-h,o-z) dimension i1( 43),i2( 43),i3( 43),i4( 43),cf( 43) dimension f12(0: 8),f13(0: 8),f23(0: 8) dimension der(3) data cf( 1)/-.1112208493660213E+02/ data i1( 1)/ 0/,i2( 1)/ 1/,i3( 1)/ 1/,i4( 1)/ 3/ data cf( 2)/0.8397428790889362E+03/ data i1( 2)/ 1/,i2( 2)/ 1/,i3( 2)/ 1/,i4( 2)/ 1/ data cf( 3)/0.1094798235180020E+03/ data i1( 3)/ 2/,i2( 3)/ 1/,i3( 3)/ 0/,i4( 3)/ 6/ data cf( 4)/-.5709190211594033E+04/ data i1( 4)/ 2/,i2( 4)/ 1/,i3( 4)/ 1/,i4( 4)/ 3/ data cf( 5)/-.1298244769404030E+04/ data i1( 5)/ 0/,i2( 5)/ 2/,i3( 5)/ 2/,i4( 5)/ 3/ data cf( 6)/-.5972616424494424E+03/ data i1( 6)/ 3/,i2( 6)/ 1/,i3( 6)/ 0/,i4( 6)/ 6/ data cf( 7)/0.3340252211548539E+05/ data i1( 7)/ 1/,i2( 7)/ 2/,i3( 7)/ 2/,i4( 7)/ 3/ data cf( 8)/0.1928328064734260E+05/ data i1( 8)/ 3/,i2( 8)/ 1/,i3( 8)/ 1/,i4( 8)/ 3/ data cf( 9)/0.6458200480851449E+04/ data i1( 9)/ 3/,i2( 9)/ 2/,i3( 9)/ 0/,i4( 9)/ 6/ data cf( 10)/0.2435060288684594E+04/ data i1( 10)/ 4/,i2( 10)/ 1/,i3( 10)/ 0/,i4( 10)/ 6/ data cf( 11)/-.1434637341827002E+06/ data i1( 11)/ 2/,i2( 11)/ 2/,i3( 11)/ 2/,i4( 11)/ 1/ data cf( 12)/-.8873294052910528E+05/ data i1( 12)/ 3/,i2( 12)/ 2/,i3( 12)/ 1/,i4( 12)/ 6/ data cf( 13)/-.2560409973248071E+05/ data i1( 13)/ 0/,i2( 13)/ 3/,i3( 13)/ 3/,i4( 13)/ 3/ data cf( 14)/-.4155091370361293E+05/ data i1( 14)/ 4/,i2( 14)/ 1/,i3( 14)/ 1/,i4( 14)/ 3/ data cf( 15)/-.1908724741559524E+05/ data i1( 15)/ 4/,i2( 15)/ 2/,i3( 15)/ 0/,i4( 15)/ 6/ data cf( 16)/-.6916964417547845E+04/ data i1( 16)/ 5/,i2( 16)/ 1/,i3( 16)/ 0/,i4( 16)/ 6/ data cf( 17)/0.2700864073926629E+06/ data i1( 17)/ 3/,i2( 17)/ 2/,i3( 17)/ 2/,i4( 17)/ 3/ data cf( 18)/0.1761420315708862E+06/ data i1( 18)/ 1/,i2( 18)/ 3/,i3( 18)/ 3/,i4( 18)/ 3/ data cf( 19)/0.1285760798314841E+06/ data i1( 19)/ 4/,i2( 19)/ 2/,i3( 19)/ 1/,i4( 19)/ 6/ data cf( 20)/0.5241468442401338E+05/ data i1( 20)/ 4/,i2( 20)/ 3/,i3( 20)/ 0/,i4( 20)/ 6/ data cf( 21)/0.6179297668523817E+05/ data i1( 21)/ 5/,i2( 21)/ 1/,i3( 21)/ 1/,i4( 21)/ 3/ data cf( 22)/0.3613316158567466E+05/ data i1( 22)/ 5/,i2( 22)/ 2/,i3( 22)/ 0/,i4( 22)/ 6/ data cf( 23)/0.1167502561479660E+05/ data i1( 23)/ 6/,i2( 23)/ 1/,i3( 23)/ 0/,i4( 23)/ 6/ data cf( 24)/-.3483578214904613E+06/ data i1( 24)/ 2/,i2( 24)/ 3/,i3( 24)/ 3/,i4( 24)/ 3/ data cf( 25)/-.1987372527942004E+06/ data i1( 25)/ 4/,i2( 25)/ 2/,i3( 25)/ 2/,i4( 25)/ 3/ data cf( 26)/-.1655884747032285E+06/ data i1( 26)/ 4/,i2( 26)/ 3/,i3( 26)/ 1/,i4( 26)/ 6/ data cf( 27)/-.7617059629974826E+05/ data i1( 27)/ 0/,i2( 27)/ 4/,i3( 27)/ 4/,i4( 27)/ 3/ data cf( 28)/-.1202657194761244E+06/ data i1( 28)/ 5/,i2( 28)/ 2/,i3( 28)/ 1/,i4( 28)/ 6/ data cf( 29)/-.4853231398644997E+05/ data i1( 29)/ 5/,i2( 29)/ 3/,i3( 29)/ 0/,i4( 29)/ 6/ data cf( 30)/-.5201212017261994E+05/ data i1( 30)/ 6/,i2( 30)/ 1/,i3( 30)/ 1/,i4( 30)/ 3/ data cf( 31)/-.4169368997779521E+05/ data i1( 31)/ 6/,i2( 31)/ 2/,i3( 31)/ 0/,i4( 31)/ 6/ data cf( 32)/-.9856384447393651E+04/ data i1( 32)/ 7/,i2( 32)/ 1/,i3( 32)/ 0/,i4( 32)/ 6/ data cf( 33)/0.4031304431979369E+06/ data i1( 33)/ 3/,i2( 33)/ 3/,i3( 33)/ 3/,i4( 33)/ 1/ data cf( 34)/0.9745826736095028E+05/ data i1( 34)/ 4/,i2( 34)/ 3/,i3( 34)/ 2/,i4( 34)/ 6/ data cf( 35)/0.9724244007877808E+05/ data i1( 35)/ 1/,i2( 35)/ 4/,i3( 35)/ 4/,i4( 35)/ 3/ data cf( 36)/0.2355094235748628E+05/ data i1( 36)/ 5/,i2( 36)/ 2/,i3( 36)/ 2/,i4( 36)/ 3/ data cf( 37)/0.7589289221958217E+05/ data i1( 37)/ 5/,i2( 37)/ 3/,i3( 37)/ 1/,i4( 37)/ 6/ data cf( 38)/0.3505234464517287E+05/ data i1( 38)/ 5/,i2( 38)/ 4/,i3( 38)/ 0/,i4( 38)/ 6/ data cf( 39)/0.6811389394620169E+05/ data i1( 39)/ 6/,i2( 39)/ 2/,i3( 39)/ 1/,i4( 39)/ 6/ data cf( 40)/0.1432035524851620E+05/ data i1( 40)/ 6/,i2( 40)/ 3/,i3( 40)/ 0/,i4( 40)/ 6/ data cf( 41)/0.1201363089340074E+05/ data i1( 41)/ 7/,i2( 41)/ 1/,i3( 41)/ 1/,i4( 41)/ 3/ data cf( 42)/0.2214599648186910E+05/ data i1( 42)/ 7/,i2( 42)/ 2/,i3( 42)/ 0/,i4( 42)/ 6/ data cf( 43)/0.2937776762252248E+04/ data i1( 43)/ 8/,i2( 43)/ 1/,i3( 43)/ 0/,i4( 43)/ 6/ vex1=0.9100010000000002E+00 f12(0)=1.d0 f13(0)=1.d0 f23(0)=1.d0 bux12=r12*dexp(-vex1*r12) bux13=r13*dexp(-vex1*r13) bux23=r23*dexp(-vex1*r23) do 1 i=1, 8 f12(i)=f12(i-1)*bux12 f13(i)=f13(i-1)*bux13 f23(i)=f23(i-1)*bux23 1 continue ener = 0.d0 der12 = 0.d0 der13 = 0.d0 der23 = 0.d0 do 2 l=1, 43 if (i4(l).eq.1) then aux=f12(i1(l))*f13(i2(l))*f23(i3(l)) dux12=i1(l)*f12(i1(l)-1)*f13(i2(l))*f23(i3(l)) dux13=i2(l)*f12(i1(l))*f13(i2(l)-1)*f23(i3(l)) dux23=i3(l)*f12(i1(l))*f13(i2(l))*f23(i3(l)-1) elseif (i4(l).eq.3) then aux1=f12(i1(l))*f13(i2(l))*f23(i3(l)) aux2=f12(i3(l))*f13(i1(l))*f23(i2(l)) aux3=f12(i2(l))*f13(i3(l))*f23(i1(l)) aux=aux1+aux2+aux3 dux1=i1(l)*f12(i1(l)-1)*f13(i2(l))*f23(i3(l)) dux2=i3(l)*f12(i3(l)-1)*f13(i1(l))*f23(i2(l)) dux3=i2(l)*f12(i2(l)-1)*f13(i3(l))*f23(i1(l)) dux12=dux1+dux2+dux3 dux1=i2(l)*f12(i1(l))*f13(i2(l)-1)*f23(i3(l)) dux2=i1(l)*f12(i3(l))*f13(i1(l)-1)*f23(i2(l)) dux3=i3(l)*f12(i2(l))*f13(i3(l)-1)*f23(i1(l)) dux13=dux1+dux2+dux3 dux1=i3(l)*f12(i1(l))*f13(i2(l))*f23(i3(l)-1) dux2=i2(l)*f12(i3(l))*f13(i1(l))*f23(i2(l)-1) dux3=i1(l)*f12(i2(l))*f13(i3(l))*f23(i1(l)-1) dux23=dux1+dux2+dux3 elseif (i4(l).eq.6) then aux1=f12(i1(l))*f13(i2(l))*f23(i3(l)) aux2=f12(i1(l))*f13(i3(l))*f23(i2(l)) aux3=f12(i2(l))*f13(i1(l))*f23(i3(l)) aux4=f12(i2(l))*f13(i3(l))*f23(i1(l)) aux5=f12(i3(l))*f13(i1(l))*f23(i2(l)) aux6=f12(i3(l))*f13(i2(l))*f23(i1(l)) aux=aux1+aux2+aux3+aux4+aux5+aux6 dux1=i1(l)*f12(i1(l)-1)*f13(i2(l))*f23(i3(l)) dux2=i1(l)*f12(i1(l)-1)*f13(i3(l))*f23(i2(l)) dux3=i2(l)*f12(i2(l)-1)*f13(i1(l))*f23(i3(l)) dux4=i2(l)*f12(i2(l)-1)*f13(i3(l))*f23(i1(l)) dux5=i3(l)*f12(i3(l)-1)*f13(i1(l))*f23(i2(l)) dux6=i3(l)*f12(i3(l)-1)*f13(i2(l))*f23(i1(l)) dux12=dux1+dux2+dux3+dux4+dux5+dux6 dux1=i2(l)*f12(i1(l))*f13(i2(l)-1)*f23(i3(l)) dux2=i3(l)*f12(i1(l))*f13(i3(l)-1)*f23(i2(l)) dux3=i1(l)*f12(i2(l))*f13(i1(l)-1)*f23(i3(l)) dux4=i3(l)*f12(i2(l))*f13(i3(l)-1)*f23(i1(l)) dux5=i1(l)*f12(i3(l))*f13(i1(l)-1)*f23(i2(l)) dux6=i2(l)*f12(i3(l))*f13(i2(l)-1)*f23(i1(l)) dux13=dux1+dux2+dux3+dux4+dux5+dux6 dux1=i3(l)*f12(i1(l))*f13(i2(l))*f23(i3(l)-1) dux2=i2(l)*f12(i1(l))*f13(i3(l))*f23(i2(l)-1) dux3=i3(l)*f12(i2(l))*f13(i1(l))*f23(i3(l)-1) dux4=i1(l)*f12(i2(l))*f13(i3(l))*f23(i1(l)-1) dux5=i2(l)*f12(i3(l))*f13(i1(l))*f23(i2(l)-1) dux6=i1(l)*f12(i3(l))*f13(i2(l))*f23(i1(l)-1) dux23=dux1+dux2+dux3+dux4+dux5+dux6 endif ener=ener+cf(l)*aux der12=der12+cf(l)*dux12 der13=der13+cf(l)*dux13 der23=der23+cf(l)*dux23 2 continue der(1)=der12*(1.d0-vex1*r12)*dexp(-vex1*r12) der(2)=der13*(1.d0-vex1*r13)*dexp(-vex1*r13) der(3)=der23*(1.d0-vex1*r23)*dexp(-vex1*r23) return end