c Ground state potential energy surface He + H2+ ----> HeH+ + H c T. Joseph and Sathyamurthy, JCP 86,704(1987). subroutine heh2pes(xd,vxd) implicit real*8 (a-h,o-z) c calc. of potential and derivatives using carter-murrel fun dimension rx(3),parm(26),vdtm(3),reql(3),d(3),a(3,3) dimension rho(3),power(2),anr(3),dr(3),deriv(3) dimension title(18),line(13),secan(2),xd(3),r(3) c CMS-NON LS. FIT TO HeHH+ PES TRITM FIT (AB+3.3992) **BEST FIT** *** CI *** data m,n/520,26/ data parm/6.34650,9.34530,-3.63380,5.17180,-2.96930, & -0.08740,6.93820,0.61410,4.78700,5.88260,-9.74530, & -3.57910,4.15630,-1.95470,-0.45990,6.43060,-4.00110, & 4.28260,-3.80790,7.96950,-0.12490,0.20040,2.22250, & 1.37490,0.95250,0.95250/ data d/2.78650,2.02790,2.02790/ data reql/1.05360,0.77600,0.77600/ data ((a(i,j),j=1,3),i=1,3)/2.61800,1.57460,1.24510,4.60930, & 2.94020,2.61850,4.60930,2.94020,2.61850/ r(1)=sngl(xd(1)) r(2)=sngl(xd(2)) r(3)=sngl(xd(3)) c 1e1,dum1,q,p,t,w35,w37,w39,w11 data line(1)/'-'/ data ij,ifirst,iderv/1,1,0/ if(ifirst.ne.1)go to 99 6 format(18a4) c read(5,6)title write(9,8)(title(i),i=1,18) 8 format(35a4) c read(5,10)m,n 10 format(2x,2i3) write(9,11)m,n 11 format(2x,'no.of points = ',i3/,2x,'no. of parms = ',i3) c read(5,100)(parm(i),i=1,n) 100 format(2x,9f8.5) write(9,21)(parm(i),i=1,n) 21 format(2x,'potential parms'/,6(4x,6f12.5/)) c read(5,100)d(1),reql(1),(a(1,i),i=1,3) c read(5,100)d(2),reql(2),(a(2,i),i=1,3) c read(5,100)d(3),reql(3),(a(3,i),i=1,3) do 20 i=1,3 20 write(9,13)d(i),reql(i),(a(i,j),j=1,3) 13 format(5x,2f13.6,3f13.6) ifirst = ifirst+1 99 continue c x1 = exp((re(1)-r(1))*alph(1)) c c z1 = exp((re(2)-r(2))*alph(2)) c x2 = x1*x1 c y2 = y1*y1 c z2 = z1*z1 c e1ab = d1(1)*(x2 - 2.0*x1) c e1bc = d1(2)*(z2 - 2.0*z1) c e1ac = d1(3)*(y2 - 2.0*y1) c e1(1) = e1ab c e1(2) = e1bc c e1(3) = e1ac ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c in the main programme r(1)=r(heh),r(2)=r(h2+),r(3)=r(heh) c in this subroutine rx(1)=r(hh),rx(2)=r(heh),r(3)=r(heh) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc angs=0.529177 rx(1)=r(2)*angs rx(2)=r(1)*angs rx(3)=r(3)*angs c rx(3) = sqrt(rx(1)*rx(1)+rx(2)*rx(2)-2.*rx(1)*rx(2)*cosd(theta)) do 50 i=1,3 rho(i) =rx(i)-reql(i) vdtm(i) =-(d(i)*(1.+a(i,1)*rho(i)+a(i,2)*rho(i)*rho(i) 1 +a(i,3)*rho(i)*rho(i)*rho(i)))*exp(-(a(i,1)*rho(i))) 50 continue vdiat = vdtm(1)+vdtm(2)+vdtm(3) c vdiat = vdiat*27.2116 rho1 =rx(1)-parm(25) rho2 =rx(2)-parm(26) rho3 =rx(3)-parm(26) rho1s= rho1*rho1 rho2s= rho2*rho2 rho3s =rho3*rho3 rho12 =rho1*rho2 rho13 =rho1*rho3 rho23 =rho2*rho3 rhosum = rho2+rho3 term1 =parm(1)+parm(2)*rho1+parm(3)*(rho2+rho3)+parm(4)*rho1s term2 =parm(5)*(rho2s+rho3s)+parm(6)*(rho12+rho13)+parm(7)*rho23 term3 =parm(8)*rho1*rho1s+parm(9)*(rho2*rho2s+rho3*rho3s) term4 =parm(10)*rho1s*(rho2+rho3)+parm(11)*rho1*(rho2s+rho3s) term5 =parm(12)*rho2*(rho23+rho3s)+parm(13)*rho12*rho3 term6 =parm(14)*rho1s*rho1s+parm(15)*(rho2s*rho2s+rho3s*rho3s) term7 =parm(16)*rho1s*(rho12+rho13)+parm(17)*rho1s*(rho2s+rho3s) term8 =parm(18)*(rho12*rho2s+rho13*rho3s) term9 =parm(19)*rho2*(rho2s*rho3+rho3s*rho3) term10 =parm(20)*rho2s*rho3s+parm(21)*rho1s*rho23 term11 =parm(22)*(rho13*rho2s+rho12*rho3s) poly1 =term1+term2+term3+term4+term5+term6+term7 poly = poly1+term8+term9+term10+term11 power(1) =parm(23)*rho1/2. power(2) =parm(24)*(rho2+rho3)/2. do 70 i=1,2 anr(i) =exp(power(i))-exp(-power(i)) dr(i) =exp(power(i))+exp(-power(i)) 70 secan(i) = 1./dr(i)/dr(i) ta1 =1.-(anr(1)/dr(1)) ta2 =1.-(anr(2)/dr(2)) ta =ta1*ta2 ycalc =poly*ta e = ycalc+vdiat+d(1) vxd=e if(iderv.eq.0)return c******* calculation of derivatives ******** do 110 i=1,3 dterm = -vdtm(i)*a(i,1) part2 = a(i,1)+2.*a(i,2)*rho(i)+3.*a(i,3)*rho(i)**2 deriv(i)=dterm-(d(i)*exp(-a(i,1)*rho(i))*part2) 110 continue d1tm1 = parm(2)+2.*parm(4)*rho1+parm(6)*rhosum d1tm2 = 3.*parm(8)*rho1s+2.*parm(10)*rho1*rhosum d1tm3 = parm(11)*(rho2s+rho3s)+parm(13)*rho23 d1tm4 = 4.*parm(14)*rho1s*rho1+3.*parm(16)*rho1s*rhosum d1tm5 = 2.*parm(17)*rho1*(rho2s+rho3s) d1tm6 = parm(18)*(rho2s*rho2+rho3s*rho3) d1tm7 = 2.*parm(21)*rho1*rho23+parm(22)*rho2*(rho23+rho3s) d1tm = d1tm1+d1tm2+d1tm3+d1tm4+d1tm5+d1tm6+d1tm7 dwrt1 = deriv(1)+ta*d1tm-poly*ta2*2.*parm(23)*secan(1) d2tm1 = parm(3)+2.*parm(5)*rho2+parm(6)*rho1+parm(7)*rho3 d2tm2 = 3.*parm(9)*rho2s+parm(10)*rho1s+2.*parm(11)*rho12 d2tm3 = 2.*parm(12)*rho23+parm(12)*rho3s+parm(13)*rho13 d2tm4 = 4.*parm(15)*rho2s*rho2+parm(16)*rho1s*rho1 d2tm5 = 2.*parm(17)*rho1s*rho2+3.*parm(18)*rho1*rho2s d2tm6 = parm(19)*rho3*(3.*rho2s+rho3s)+2.*parm(20)*rho2*rho3s d2tm7 = parm(21)*rho1s*rho3+parm(22)*rho1*(2.*rho23+rho3s) d2tm = d2tm1+d2tm2+d2tm3+d2tm4+d2tm5+d2tm6+d2tm7 dwrt2 = deriv(2)+ta*d2tm-poly*ta1*2.*parm(24)*secan(2) d3tm1 = parm(3)+2.*parm(5)*rho3+parm(6)*rho1+parm(7)*rho2 d3tm2 = 3.*parm(9)*rho3s+parm(10)*rho1s+2.*parm(11)*rho13 d3tm3 = parm(12)*(rho2s+2.*rho23)+parm(13)*rho12 d3tm4 = 4.*parm(15)*rho3s*rho3+parm(16)*rho1s*rho1 d3tm5 = 2.*parm(17)*rho1s*rho3+3.*parm(18)*rho1*rho3s d3tm6 = parm(19)*rho2*(rho2s+3.*rho3s)+2.*parm(20)*rho2s*rho3 d3tm7 = parm(21)*rho1s*rho2+parm(22)*rho1*(rho2s+2.*rho23) d3tm = d3tm1+d3tm2+d3tm3+d3tm4+d3tm5+d3tm6+d3tm7 dwrt3 = deriv(3)+ta*d3tm-poly*ta1*2.*parm(24)*secan(2) c dedr2 =(dwrt1+dwrt3)*angs c dedr1 =(dwrt2+dwrt3)*angs c dedr3 =0.0 dedr1 = dwrt2*angs dedr2 = dwrt1*angs dedr3 = dwrt3*angs C write(9,210)r(1),r(2),r(3),e,dedr1,dedr2,dedr3 c210 format(2x,8e16.4) return end