subroutine fit3d(r12,r13,r23,e,der) implicit real * 8 (a-h,o-z) dimension der(3) call diat12(r12,e12,d12) call diat12(r13,e13,d13) call diat23(r23,e23,d23) call triabb(r12,r13,r23,e123,der) e=e12+e13+e23+e123 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 38 points * rms = 0.01306248 kcal/mol * emax = 0.04102054 kcal/mol ************************************************************************ implicit real*8 (a-h,o-z) dimension cf( 6) data cf( 1)/0.227918348133E+01/ data cf( 2)/-.426241519342E-01/ data cf( 3)/-.277033983580E+00/ data cf( 4)/-.484220650845E+01/ data cf( 5)/0.137726092882E+02/ data cf( 6)/-.151873900316E+02/ e0= 0.0000000E+00 der=0.d0 vex1= 0.1073343E+01 vex2= 0.2926225E+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, 6 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 diat23(r,ener,der) ************************************************************************ * This subroutine computes the energies of a diatomic potential * fitted to 36 points * rms = 0.01785098 kcal/mol * emax = 0.04658124 kcal/mol ************************************************************************ implicit real*8 (a-h,o-z) dimension cf( 6) data cf( 1)/0.100944747431E+01/ data cf( 2)/-.704880722032E+00/ data cf( 3)/0.205248231684E+01/ data cf( 4)/-.656461128235E+01/ data cf( 5)/0.125212270626E+02/ data cf( 6)/-.998059288920E+01/ e0= 0.0000000E+00 der=0.d0 vex1= 0.9768039E+00 vex2= 0.1634212E+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, 6 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 triabb(r12,r13,r23,ener,der) ************************************************************* * This subroutine computes the energies of a 3D PES * for the ABB system class fitted to 1511 points * rms = 0.32654048 kcal/mol * emax = 2.90375166 kcal/mol ************************************************************* implicit real*8(a-h,o-z) dimension i1( 37),i2( 37),i3( 37),i4( 37),cf( 37) dimension f12(0: 5),f13(0: 5),f23(0: 5) dimension der(3) data cf( 1)/0.1147027902313981E+01/ data i1( 1)/ 1/,i2( 1)/ 1/,i3( 1)/ 0/,i4( 1)/ 2/ data cf( 2)/0.1372628580809719E+01/ data i1( 2)/ 1/,i2( 2)/ 0/,i3( 2)/ 1/,i4( 2)/ 1/ data cf( 3)/0.1980014275534632E+02/ data i1( 3)/ 1/,i2( 3)/ 1/,i3( 3)/ 1/,i4( 3)/ 1/ data cf( 4)/0.5379201080452515E+01/ data i1( 4)/ 2/,i2( 4)/ 1/,i3( 4)/ 0/,i4( 4)/ 2/ data cf( 5)/-.3832448627603355E+01/ data i1( 5)/ 2/,i2( 5)/ 0/,i3( 5)/ 1/,i4( 5)/ 2/ data cf( 6)/-.1545142894835985E+02/ data i1( 6)/ 0/,i2( 6)/ 2/,i3( 6)/ 1/,i4( 6)/ 2/ data cf( 7)/0.1637473417369147E+03/ data i1( 7)/ 2/,i2( 7)/ 1/,i3( 7)/ 1/,i4( 7)/ 2/ data cf( 8)/-.2645243308320998E+03/ data i1( 8)/ 1/,i2( 8)/ 2/,i3( 8)/ 1/,i4( 8)/ 1/ data cf( 9)/-.3079377355246953E+01/ data i1( 9)/ 2/,i2( 9)/ 2/,i3( 9)/ 0/,i4( 9)/ 2/ data cf( 10)/-.2996884803782916E+02/ data i1( 10)/ 2/,i2( 10)/ 0/,i3( 10)/ 2/,i4( 10)/ 1/ data cf( 11)/-.4997515477576375E+01/ data i1( 11)/ 3/,i2( 11)/ 1/,i3( 11)/ 0/,i4( 11)/ 2/ data cf( 12)/-.5428562192603896E+02/ data i1( 12)/ 3/,i2( 12)/ 0/,i3( 12)/ 1/,i4( 12)/ 2/ data cf( 13)/0.7030596418655885E+02/ data i1( 13)/ 0/,i2( 13)/ 3/,i3( 13)/ 1/,i4( 13)/ 2/ data cf( 14)/-.1829653141104358E+03/ data i1( 14)/ 2/,i2( 14)/ 2/,i3( 14)/ 1/,i4( 14)/ 2/ data cf( 15)/-.1815811115480261E+03/ data i1( 15)/ 2/,i2( 15)/ 1/,i3( 15)/ 2/,i4( 15)/ 1/ data cf( 16)/-.2688659034849295E+03/ data i1( 16)/ 3/,i2( 16)/ 1/,i3( 16)/ 1/,i4( 16)/ 2/ data cf( 17)/0.6996258916985428E+03/ data i1( 17)/ 1/,i2( 17)/ 3/,i3( 17)/ 1/,i4( 17)/ 1/ data cf( 18)/0.1529276303827616E+03/ data i1( 18)/ 3/,i2( 18)/ 2/,i3( 18)/ 0/,i4( 18)/ 2/ data cf( 19)/0.1128756109189617E+03/ data i1( 19)/ 3/,i2( 19)/ 0/,i3( 19)/ 2/,i4( 19)/ 2/ data cf( 20)/-.1232147843694094E+03/ data i1( 20)/ 0/,i2( 20)/ 3/,i3( 20)/ 2/,i4( 20)/ 2/ data cf( 21)/-.9399744188456941E+02/ data i1( 21)/ 4/,i2( 21)/ 1/,i3( 21)/ 0/,i4( 21)/ 2/ data cf( 22)/0.1503166731921269E+03/ data i1( 22)/ 4/,i2( 22)/ 0/,i3( 22)/ 1/,i4( 22)/ 2/ data cf( 23)/-.1168862115918631E+03/ data i1( 23)/ 0/,i2( 23)/ 4/,i3( 23)/ 1/,i4( 23)/ 2/ data cf( 24)/-.3129363925426480E+03/ data i1( 24)/ 2/,i2( 24)/ 2/,i3( 24)/ 2/,i4( 24)/ 1/ data cf( 25)/0.1466861135826718E+04/ data i1( 25)/ 3/,i2( 25)/ 2/,i3( 25)/ 1/,i4( 25)/ 2/ data cf( 26)/0.4552071379940093E+03/ data i1( 26)/ 3/,i2( 26)/ 1/,i3( 26)/ 2/,i4( 26)/ 2/ data cf( 27)/-.5473616801633061E+03/ data i1( 27)/ 1/,i2( 27)/ 3/,i3( 27)/ 2/,i4( 27)/ 2/ data cf( 28)/-.2245601795436739E+03/ data i1( 28)/ 3/,i2( 28)/ 3/,i3( 28)/ 0/,i4( 28)/ 2/ data cf( 29)/0.1616003964945219E+02/ data i1( 29)/ 3/,i2( 29)/ 0/,i3( 29)/ 3/,i4( 29)/ 1/ data cf( 30)/-.7661081485095968E+03/ data i1( 30)/ 4/,i2( 30)/ 1/,i3( 30)/ 1/,i4( 30)/ 2/ data cf( 31)/-.3896061458270143E+03/ data i1( 31)/ 1/,i2( 31)/ 4/,i3( 31)/ 1/,i4( 31)/ 1/ data cf( 32)/0.9145578840021393E+00/ data i1( 32)/ 4/,i2( 32)/ 2/,i3( 32)/ 0/,i4( 32)/ 2/ data cf( 33)/-.2509670273471583E+03/ data i1( 33)/ 4/,i2( 33)/ 0/,i3( 33)/ 2/,i4( 33)/ 2/ data cf( 34)/0.2220285111685380E+03/ data i1( 34)/ 0/,i2( 34)/ 4/,i3( 34)/ 2/,i4( 34)/ 2/ data cf( 35)/0.1500218991718269E+03/ data i1( 35)/ 5/,i2( 35)/ 1/,i3( 35)/ 0/,i4( 35)/ 2/ data cf( 36)/0.2021220548231449E+02/ data i1( 36)/ 5/,i2( 36)/ 0/,i3( 36)/ 1/,i4( 36)/ 2/ data cf( 37)/0.5143092110746637E+02/ data i1( 37)/ 0/,i2( 37)/ 5/,i3( 37)/ 1/,i4( 37)/ 2/ vex1=0.1053863292070094E+01 vex2=0.8545830504451079E+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(-vex2*r23) do 1 i=1, 5 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, 37 if (i4(l).eq.1) then aux=f12(i1(l))*f13(i3(l))*f23(i2(l)) dux12=i1(l)*f12(i1(l)-1)*f13(i3(l))*f23(i2(l)) dux13=i3(l)*f12(i1(l))*f13(i3(l)-1)*f23(i2(l)) dux23=i2(l)*f12(i1(l))*f13(i3(l))*f23(i2(l)-1) else aux1=f12(i1(l))*f13(i3(l)) aux2=f12(i3(l))*f13(i1(l)) aux=(aux1+aux2)*f23(i2(l)) dux23=(aux1+aux2)*i2(l)*f23(i2(l)-1) dux1=i1(l)*f12(i1(l)-1)*f13(i3(l)) dux2=i3(l)*f12(i3(l)-1)*f13(i1(l)) dux12=(dux1+dux2)*f23(i2(l)) dux1=i3(l)*f12(i1(l))*f13(i3(l)-1) dux2=i1(l)*f12(i3(l))*f13(i1(l)-1) dux13=(dux1+dux2)*f23(i2(l)) 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-vex2*r23)*dexp(-vex2*r23) return end