*deck potbo c c this deck potbo contains: c the interface routine pot c vbase adjusts the energy zero c convr converts internuclear distances to a c convv converts potential values to kcal.mol c r are internucl. distances in the order rac,rab,rbc c vv is the output potential value in desired units c the routine prepot for reading potential parameters c the routine potbo for calculating bo potentials c as a function of rab,rbc,rac (in a). it returns as v the c value of the potential in kcalmol. the energy zero is set c at total dissociation. c the soubroutine pol3 for calculating the bo polynomial c c subroutine prepot c implicit real*8(a-h,o-z) common/b21/ de1,de2,de3,re1,re2,re3 > ,b1,b2,b3 common/b23/cof(65),nt,ie(65),je(65),ke(65) dimension coef(65),iexp(65),jexp(65),kexp(65) common/b22/a11,a12,a13,a21,a22,a23,a31,a32,a33 > ,a14,a24,a34 common/b25/ pg common/btanh/ref1,ref2,ref3,sl1,sl2,sl3,r01,r02,r03 common/repuls/a,n1,n2,n3 common/refine/agb,cgbx,cgby,cgbz,agw,cgwx,cgwy,cgwz, > xn0b,yn0b,zn0b,xn0w,yn0w,zn0w c data coef / * 0.9224052d+03, 0.3272860d+03, 0.3381706d+03, -0.1048643d+04, * -0.2962569d+03, -0.3839140d+03, -0.7478068d+03, -0.7533124d+03, * -0.1715470d+03, -0.2134407d+01, 0.5508481d+03, 0.1134382d+04, * -0.1160042d+04, 0.3557390d+03, 0.1017086d+02, 0.6969160d+03, * 0.1258877d+04, 0.5197582d+02, 0.4328018d+03, -0.1062101d+03, * 0.5471942d+02, -0.3211568d+02, -0.1333515d+03, -0.4993582d+03, * -0.3928734d+03, 0.1143333d+04, -0.2252141d+03, 0.2641766d+03, * 0.9964729d+02, -0.3031654d+03, -0.9919408d+03, 0.7123589d+02, * -0.5026490d+03, -0.2019268d+03, 0.1085769d+03, -0.3689024d+02, * 0.1259192d+03, 0.1398510d+02, -0.5489551d+01, -0.2721741d+02, * -0.8672702d+01, 0.9243725d+02, 0.9455495d+02, 0.4378558d+02, * -0.3012895d+03, 0.5875433d+02, -0.7644942d+01, -0.1544242d+03, * 0.1174672d+01, 0.1881099d+02, 0.2699340d+03, -0.3803897d+02, * 0.2211474d+02, 0.1713199d+03, 0.6709845d+02, -0.2999395d+02, * 0.1285414d+02, 0.5457174d+02, -0.9745477d+02, -0.7197123d+01, * -0.1377454d+01, -0.2120028d+02, 0.2071087d+02, 0.7003032d+00, * 0.4229195d+01/ data iexp/ * 1, 1, 0, 2, 1, 2, 1, 0, 1, 0, 3, 2, 1, 3, 2, 1, 0, 2, 1, 0, * 1, 0, 4, 3, 2, 1, 4, 3, 2, 1, 0, 3, 2, 1, 0, 2, 1, 0, 1, 0, * 5, 4, 3, 2, 1, 5, 4, 3, 2, 1, 0, 4, 3, 2, 1, 0, 3, 2, 1, 0, * 2, 1, 0, 1, 0/ data jexp/ * 1, 0, 1, 1, 2, 0, 1, 2, 0, 1, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, * 0, 1, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 0, 1, 2, 0, 1, * 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 0, 1, 2, 3, * 0, 1, 2, 0, 1/ data kexp/ * 0, 1, 1, 0, 0, 1, 1, 1, 2, 2, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, * 3, 3, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, * 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, * 4, 4, 4, 5, 5/ c data vvre,vvr0,eps/1d-4,.01d0,.001d0/ data tailbx,tailby,tailbz/.111d0,.037d0,.078d0/ data tailwx,tailwy,tailwz/.088d0,.096d0,.042d0/ data x0b,y0b,angb/1.6523d0,1.3715d0,72.5235d0/ data x0w,y0w,angw/1.8848d0,0.9363d0,106.929d0/ agw = -1.8d0 agb = -0.1d0 c 2 format(i2) 3 format(e15.7,3i3) c pg=3.1415926536d0 c c dissociation energies de1=137.59d0 de2=141.20d0 de3=58.d0 c c b parameters c b1=.970755d0 b2=2.194185d0 b3=1.17086d0 c c equilibrium distances re1=1.563864d0 re2=0.916808d0 re3=1.5957d0 c c diatomic parameters a11=2.30439d0 a12=-2.30172d0 a13=1.69028d0 a14=-.69295d0 a21=2.07810d0 a22=-1.25667d0 a23=0.27905d0 a24=-0.10047d0 a31=2.15123d0 a32=-1.38219d0 a33=0.31067d0 a34=-.07972d0 c c finds the largest short range distance at which diatomic c potentials are zero (r01,r02,r03) do 1 i=1,3 if(i.ne.1)goto 63 b=b1 re=re1 de=de1 a1=a11 a2=a12 a3=a13 a4=a14 goto67 63 continue if(i.ne.2)goto65 b=b2 re=re2 de=de2 a1=a21 a2=a22 c-------------------------------------------------------------- c Make sure a3=0.0 and not a3=a23 c-------------------------------------------------------------- ccccc a3=a23 a3=0.d0 a4=a24 goto 67 65 continue b=b3 re=re3 de=de3 a1=a31 a2=a32 a3=a33 a4=a34 67 continue dx=-.1d0 x=re-dx 70 continue x=x+dx ex=exp(-b*(x-re)) v=-de*ex*(a1+ex*(a2+ex*(a3+ex*a4))) if(v.lt.0d0)goto70 if(abs(dx).lt.0.0001d0)goto110 x=x-dx dx=.1d0*dx goto70 110 continue if(i.eq.1)r01=x if(i.eq.2)r02=x if(i.eq.3)r03=x 1 continue c c determines the value at which tanh has a flex ref1=.49d0*r01+.51d0*re1 ref2=.49d0*r02+.51d0*re2 ref3=.49d0*r03+.51d0*re3 c c determines the slope of tanh at the flex alneps=log(1d0/eps-1d0) sl1=.5d0*alneps/(re1-ref1) sl2=.5d0*alneps/(re2-ref2) sl3=.5d0*alneps/(re3-ref3) c c determines the parameters for the additional repulsive inner c diatomic term the coefficient a and the powers n1,n2,n3 c a=vvre n1=log(vvr0/vvre)/(-b1*(r01-re1)) n2=log(vvr0/vvre)/(-b2*(r02-re2)) n3=log(vvr0/vvre)/(-b3*(r03-re3)) c c determines the parameters for the additional repulsive short range c terms of diatoms (given as D*a*BO**n) c a=vvre n1=log(vvr0/vvre)/(-b1*(r01-re1)) n2=log(vvr0/vvre)/(-b2*(r02-re2)) n3=log(vvr0/vvre)/(-b3*(r03-re3)) c c calculates the parameters for gaussian bond order corrections to c the barrier and the well c gamb=angb*pg/180.d0 gamw=angw*pg/180.d0 z0b=sqrt(x0b*x0b+y0b*y0b-2.d0*x0b*y0b*cos(gamb)) z0w=sqrt(x0w*x0w+y0w*y0w-2.d0*x0w*y0w*cos(gamw)) xn0b=exp(-b1*(x0b-re1)) yn0b=exp(-b2*(y0b-re2)) zn0b=exp(-b3*(z0b-re3)) xn0w=exp(-b1*(x0w-re1)) yn0w=exp(-b2*(y0w-re2)) zn0w=exp(-b3*(z0w-re3)) cl1=log(.5d0) cgbx=-cl1/(tailbx*tailbx) cl1=log(.5d0) cgby=-cl1/(tailby*tailby) cl1=log(.5d0) cgbz=-cl1/(tailbz*tailbz) cl1=log(.5d0) cgwx=-cl1/(tailwx*tailwx) cl1=log(.5d0) cgwy=-cl1/(tailwy*tailwy) cl1=log(.5d0) cgwz=-cl1/(tailwz*tailwz) c c coefficients for the three-body term c nt=65 do 10 i=1,nt cof(i)=coef(i) ie(i)=1+iexp(i) je(i)=1+jexp(i) ke(i)=1+kexp(i) c write(6,*)'prepot', cof(i),ie(i),je(i),ke(i) 10 continue c return end