SUBROUTINE potv(V,R1,R2,xcos) IMPLICIT DOUBLE PRECISION (A-H,O-Y),logical(z) COMMON /MASS/ XMASS(3),G1,G2,xmassr(3) dimension pv(32),u(32) dimension ht(50) DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/, CMTOAU/219474.624/ DATA A/2.22600/ DATA X1/1.0D0/,X0/0.0D0/,TINY/9.0D-15/,X2/2.0D0/ DATA pv/0.999963011541667,-0.000831585900216,0.000440024019720, *-0.000790390661341,0.000682751133154,-0.000787433811192, *0.000137727815420,-0.000051565403109,-0.000122703443403, *0.003387614188353,-0.000271677251731,-0.004052838735161, *0.001707147193552,-0.001805888500169,-0.006788069703931, *0.001208782633980,0.001108778944947,0.009632702226484, *-0.009093431064531,-0.000836179586411,0.001775570609778, *0.000051292593874,0.006623947296453,-0.001827850224422, *-0.012626322010510,-0.003304654368357,0.010437389017074, *0.000227788649284,0.008030675091008, *-0.000723724819511,0.002063490878186,0.001356187797029/ pi=dacos(-1.d0) ! Equilibrium IF (G1 .EQ. X0) THEN ! BONDLENGTH BONDANGLE COORDINATES: ATOM 1 = ATOM 2 Q1 = R1 Q2 = R2 THETA = ACOS(XCOS) ELSE IF (G2 .EQ. X0) THEN ! SCATTERING COORDINATES: ATOM 2 = ATOM 3 XX = R1 * G1 YY = R1 * (X1 - G1) IF (R2 .EQ. X0 .OR. XCOS .GE. (X1 - TINY)) THEN Q1 = ABS(XX - R2) Q2 = (YY + R2) COST = -X1 ELSE IF (XCOS .LE. (TINY - X1)) THEN Q1 = (XX + R2) Q2 = ABS(YY + R2) COST = X1 ELSE Q1 = SQRT(XX*XX + R2*R2 - X2*XX*R2*XCOS) Q2 = SQRT(YY*YY + R2*R2 + X2*YY*R2*XCOS) COST = (Q1**2 + Q2**2 - R1**2) / (X2 * Q1 * Q2) ENDIF THETA = ACOS(COST) ELSE ! GENERAL COORDINATES (INCLUDING RADAU): ATOM 1 = ATOM 2 F1= X1/G1 F2= X1/G2 F12= X1 - F1*F2 P1= R1*(X1-F1)/(G2*F12) P2= R2*(X1-F2)/(G1*F12) S1= R1-P1 S2= R2-P2 Q1= SQRT(P1*P1 + S2*S2 + X2*P1*S2*XCOS)/(X1-G1) Q2= SQRT(P2*P2 + S1*S1 + X2*P2*S1*XCOS)/(X1-G2) Q3= SQRT(P1*P1 + P2*P2 - X2*P1*P2*XCOS) COST = (Q1*Q1 + Q2*Q2 - Q3*Q3)/(X2*Q1*Q2) THETA = ACOS(COST) ENDIF THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/TOANG Y1 = Q1 - Req Y2 = Cos(THETA) - Cos(THETAeq) Y3 = Q2 - Req S1 = (Y1 + Y3)/2.0d0 S2 = Y2 S3 = (Y1 - Y3)/2.0d0 ! Potential symmetrical in S1,S2,S3. Full set of ! derivatives for each order (1-4) and upper terms for others. HT(1)=1.0d0 ht(2)=s1*s2 ht(3)=s1**2 ht(4)=s1**3 ht(5)=s3**2 ht(6)=s1**4 ht(7)=s2**2 ht(8)=s2**3 ht(9)=s2**4 ht(10)=s1**2*s2 ht(11)=s2**2*s1 ht(12)=s3**2*s1 ht(13)=s3**2*s2 ht(14)=s3**4 ht(15)=s1**3*S2 ht(16)=S2**3*s1 ht(17)=s1**2*s2**2 ht(18)=s1**2*s3**2 ht(19)=s1*s2*s3**2 ht(20)=s2**2*s3**2 ht(21)=s1**5 ht(22)=s2**5 ht(23)=s1**4*s2 ht(24)=s1**3*s2**2 ht(25)=s1**3*s3**2 ht(26)=s1**2*s2**3 ht(27)=s1**2*s2*s3**2 ht(28)=s1*s2**4 ht(29)=s1*s3**4 ht(30)=s2*s3**4 ht(31)=s1*s2**2*s3**2 ht(32)=s2**3*s3**2 vt = 0.0d0 do 89 i = 1,32 u(i) = (ht(i) * pv(i)) vt = vt + u(i) 89 continue call pants(vn,r1,r2,xcos,vbn) ! This is where the ab initio and empirical potentials are defined ! vn: ab initio ! vbn: empirically determined barrier to linearity ! vt: empirical morphing of the ab initio PES v = vn*vt ! Adiabatic correction for H2-18O xmaso16=15.990526 xmaso17=16.9947425 xmaso18=17.9947714 xmasd=2.013553214 xmash=1.00727647 CALL bodc16(V16,Q1,Q2,THETA) CALL bodc18(V18,Q1,Q2,THETA) Vad16=(xmaso18*xmaso16*(v16-v18)/(xmaso18-xmaso16))/xmaso16+ *(2*xmash*(xmaso18*v18-xmaso16*v16)/(2*(xmaso18-xmaso16)))/xmash Vad18=(xmaso18*xmaso16*(v16-v18)/(xmaso18-xmaso16))/xmaso18+ *(2*xmash*(xmaso18*v18-xmaso16*v16)/(2*(xmaso18-xmaso16)))/xmash v=v-Vad16+Vad18 end SUBROUTINE pants(V,R1,R2,xcos,vbn) implicit double precision(a-h,o-y),logical(z) common /mass/ xmass(3),g1,g2,xmassr(3) dimension hp(84),hc(84),jp(84),pp(84),uz(84) dimension rij(3) DATA X1/1.0D0/,X0/0.0D0/,TINY/9.0D-15/,X2/2.0D0/ IF (G1 .EQ. X0) THEN C BONDLENGTH BONDANGLE COORDINATES: ATOM 1 = ATOM 2 Q1 = R1 Q2 = R2 THETA = ACOS(XCOS) ELSE IF (G2 .EQ. X0) THEN C SCATTERING COORDINATES: ATOM 2 = ATOM 3 XX = R1 * G1 YY = R1 * (X1 - G1) IF (R2 .EQ. X0 .OR. XCOS .GE. (X1 - TINY)) THEN Q1 = ABS(XX - R2) Q2 = (YY + R2) COST = -X1 ELSE IF (XCOS .LE. (TINY - X1)) THEN Q1 = (XX + R2) Q2 = ABS(YY + R2) COST = X1 ELSE Q1 = SQRT(XX*XX + R2*R2 - X2*XX*R2*XCOS) Q2 = SQRT(YY*YY + R2*R2 + X2*YY*R2*XCOS) COST = (Q1**2 + Q2**2 - R1**2) / (X2 * Q1 * Q2) ENDIF THETA = ACOS(COST) ELSE C GENERAL COORDINATES (INCLUDING RADAU): ATOM 1 = ATOM 2 F1= X1/G1 F2= X1/G2 F12= X1 - F1*F2 P1= R1*(X1-F1)/(G2*F12) P2= R2*(X1-F2)/(G1*F12) S1= R1-P1 S2= R2-P2 Q1= SQRT(P1*P1 + S2*S2 + X2*P1*S2*XCOS)/(X1-G1) Q2= SQRT(P2*P2 + S1*S1 + X2*P2*S1*XCOS)/(X1-G2) Q3= SQRT(P1*P1 + P2*P2 - X2*P1*P2*XCOS) COST = (Q1*Q1 + Q2*Q2 - Q3*Q3)/(X2*Q1*Q2) THETA = ACOS(COST) ENDIF vbn=1.0 if(q1.gt.4.0) q1=4.0 if(q2.gt.4.0) q2=4.0 CALL BREITB3lin(Vbr,Q1,Q2,THETA) CALL PESleq6(VQ,Q1,Q2,THETA) CALL pesd2x(Vd,Q1,Q2,THETA) CALL bodc16(V16,Q1,Q2,THETA) c rij -- bohrs, THETA -- rad rij(1) = Q1 rij(2) = Q2 rij(3) = THETA Call rel(Vr,rij(1),rij(2),rij(3)) call pots(vp,rij(1),rij(2),rij(3)) c rij -- angstrems rij(1) = Q1*0.5291772 rij(2) = Q2*0.5291772 rij(3) = THETA Call wifin(Vai2,rij(1),rij(2),rij(3)) c v -- hartries, vai2 -- 1/cm v=vai2/219474.624+vp+vr+vbr+vd+v16+vq RETURN END c Subroutine for calculation of BODC for H2-16O with Schwenke c and SCF data. SUBROUTINE BODC16(V,R1,R2,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FT(119), cv(119) DATA C0/0.0D0/,SCALE/1.0D-6/ data nv/119/ DATA ZERO/0.0D0/ DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/ DATA cv/2785.4109260806,-90.0463360612,48.5226385100, *201.6994168988,-19.4625405248,217.8639531832, *-48.9272565544,-76.0803927294,25.9675092231,-65.2206704082, *60.7284555026,-386.5173198637,-2.0352745279,3.5405842322, *-32.1830475213,126.8885218953,60.2606780157,134.6300248947, *-14.0991224631,560.8792812623,27.8464660146,-67.1464124186, *-181.5695164733,43.7287769150,-150.7468776877,-103.2652779984, *-159.6432245839,-24.8956181482,185.4825335893,-231.4775497546, *-51.9208759775,3.7066140140,-212.4422129941,207.1884490703, *383.1659239100,50.8641728660,35.1754939408,127.2280510484, *-154.4544312699,55.5787967758,282.6216945851,116.5606405651, *5.4433254144,-107.1461094167,-173.8895717556,-26.5859674990, *-560.8756697840,-237.7109212157,143.9462552048,-592.3478209334, *0.0,-198.6468835005,-19.9674473372,-14.1731270087,193.4510720304, *4.6347021028,32.9502486772,-221.1685318724,26.4090449111, *-268.432837385,-147.1422366151,133.5465868568,363.9554096142, *673.9006484856,214.9454642643,40.7735822438,65.2246188257, *173.0708970426,1.9795259929,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/ c Equilibrium THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/toang S1 = R1 - Req S2 = DCOS(THETA) - DCOS(THETAeq) S3 = R2 - Req Y1 = (S1 + S3)/2.0d0 Y2 = S2 Y3 = (S1 - S3)/2.0d0 c Potential in Y1,Y2,Y3. Full set of c derivatives for each order (1-4) and upper terms for others. DO 88 I5=1,NV 88 FT(I5)=0.0D0 FT(1)=1.0d0 FT(2)=Y1 FT(3)=Y2 FT(4)=Y1**2 FT(5)=Y2**2 FT(6)=Y3**2 FT(7)=Y1*Y2 FT(8)=Y1**3 FT(9)=Y2**3 FT(10)=Y1**2*Y2 FT(11)=Y2**2*Y1 FT(12)=Y3**2*Y1 FT(13)=Y3**2*Y2 FT(14)=Y1**4 FT(15)=Y2**4 FT(16)=Y3**4 FT(17)=Y1**3*Y2 FT(18)=Y2**3*Y1 FT(19)=Y1**2*Y2**2 FT(20)=Y1**2*Y3**2 FT(21)=Y2**2*Y3**2 FT(22)=Y3**2*Y1*Y2 FT(23)=Y1**5 FT(24)=Y2**5 FT(25)=Y1**4*Y2 FT(26)=Y2**4*Y1 FT(27)=Y3**4*Y1 FT(28)=Y3**4*Y2 FT(29)=Y1**3*Y2**2 FT(30)=Y1**3*Y3**2 FT(31)=Y2**3*Y1**2 FT(32)=Y2**3*Y3**2 FT(33)=Y1**2*Y2*Y3**2 FT(34)=Y1*Y2**2*Y3**2 FT(35)=Y1**6 FT(36)=Y2**6 FT(37)=Y3**6 FT(38)=Y1**5*Y2 FT(39)=Y2**5*Y1 FT(40)=Y1**4*Y2**2 FT(41)=Y2**4*Y1**2 FT(42)=Y2**4*Y3**2 FT(43)=Y3**4*Y2**2 FT(44)=Y1**4*Y3**2 FT(45)=Y3**4*Y1**2 FT(46)=Y3**4*Y1*Y2 FT(47)=Y1**3*Y2**3 FT(48)=Y2**3*Y1**2*Y3**2 FT(49)=Y1**3*Y3**2*Y2 FT(50)=Y2**3*Y3**2*Y1 FT(51)=Y1**2*Y2**2*Y1**2 FT(52)=Y1**7 FT(53)=Y2**7 FT(54)=Y1**6*Y2 FT(55)=Y2**6*Y1 FT(56)=Y3**6*Y1 FT(57)=Y3**6*Y2 FT(58)=Y1**5*Y2**2 FT(59)=Y1**5*Y3**2 FT(60)=Y2**5*Y1**2 FT(61)=Y2**5*Y3**2 FT(62)=Y1**4*Y2*Y3**2 FT(63)=Y1**4*Y2**3 FT(64)=Y2**4*Y1*Y3**2 FT(65)=Y2**4*Y1**3 FT(66)=Y3**4*Y1*Y2**2 FT(67)=Y3**4*Y2*Y1**2 FT(68)=Y3**4*Y1**3 FT(69)=Y3**4*Y2**3 FT(70)=Y1**3*Y3**2*Y2**2 FT(71)=Y1**8 FT(72)=Y2**8 FT(73)=Y3**8 FT(74)=Y1**7*Y2 FT(75)=Y2**7*Y1 FT(76)=Y1**6*Y2**2 FT(77)=Y1**6*Y3**2 FT(78)=Y2**6*Y1**2 FT(79)=Y2**6*Y3**2 FT(80)=Y3**6*Y1**2 FT(81)=Y3**6*Y2**2 FT(82)=Y3**6*Y2*Y1 FT(83)=Y1**5*Y3**2*Y2 FT(84)=Y1**5*Y2**3 FT(85)=Y2**5*Y1**3 FT(86)=Y2**5*Y3**2*Y1 FT(87)=Y1**4*Y2**4 FT(88)=Y1**4*Y3**4 FT(89)=Y1**4*Y2**2*Y3**2 FT(90)=Y2**4*Y3**4 FT(91)=Y2**4*Y3**2*Y1**2 FT(92)=Y3**4*Y1**2*Y2**2 FT(93)=Y3**4*Y1*Y2**3 FT(94)=Y3**4*Y1**3*Y2 FT(95)=Y1**3*Y2**3*Y3**2 FT(96)=Y2**9 FT(97)=Y2**9*Y1*Y3**4 FT(98)=Y2**9*Y1**5 FT(99)=Y2**10 FT(100)=Y2**10*Y1**2*Y3**2 FT(101)=Y2**10*Y3**4 FT(102)=Y2**10*Y1**4 FT(103)=Y2**11 FT(104)=Y2**11*Y1*Y3**2 FT(105)=Y2**11*Y1**3 FT(106)=Y2**12 FT(107)=Y2**12*Y1**2 FT(108)=Y2**12*Y3**2 FT(109)=Y2**13 FT(110)=Y2**13*Y1 FT(111)=Y2**14 FT(112)=Y2**8*Y1**6 FT(113)=Y2**8*Y3**6 FT(114)=Y2**8*Y1**4*Y3**2 FT(115)=Y2**8*Y3**4*Y1**2 FT(116)=Y2**7*Y1**7 FT(117)=Y2**7*Y1**5*Y3**2 FT(118)=Y2**7*Y1**3*Y3**4 FT(119)=Y2**7*Y1*Y3**6 V=ZERO DO 40 I=1,NV 40 V=V+CV(I)*FT(I) C SCALE AND SHIFT THE ZERO V=C0+SCALE*V RETURN END SUBROUTINE BODC18(V,R1,R2,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FT(117), cv(117) DATA C0/0.0D0/,SCALE/1.0D-6/ data nv/117/ DATA ZERO/0.0D0/ DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/ DATA cv/2520.8349350538,-91.1209324309,49.9198880161, *201.5986877584,-27.0761946057,213.0632778654,-42.3461775982, *-89.5081758955,3.1590606555,-82.9610483797,91.2462747271, *-345.6845385214,2.9125633339,-9.2743637228,19.1249589270, *126.8439173732,102.4748093797,29.1206937761,-0.4903235370, *548.5647508851,57.4425043137,-65.3115364766,-32.7831814238, *111.9709105628,-278.4048997815,-231.4930911427,-323.4608284981, *-98.0602324965,181.1676886771,-282.2900662875,167.5274873935, *47.8266275957,-213.7402930366,20.8555467237,169.6474418937, *-62.1263131255,71.9040085206,72.1511749317,123.9432784777, *173.6777915640,58.5107154936,-27.8257872910,99.6002434065, *-135.2045098416,145.2475743775,-23.3633794048,-599.5375919902, *-540.4989859645,275.5466454040,-502.8009780288,167.2307803194, *-109.3097362364,-77.1423552840,492.0803282409,313.6574722992, *-22.4885664677,141.5175472800,-1134.9499411146,88.6368503668, *-587.9451428738,-247.5329792460,-4.9333140627,322.2004639448, *916.1768044740,479.0204491692,93.9437933859,107.5055433254, *0.0,63.1640462956,191.6712842623,1.6032078252,78.8696064123, *1.0174252672,-377.6585243760,-192.1212462644,653.6712620983, *0.0,261.5131272695,146.4973633837,0.0,-137.1867814770, *-98.7015480745,0.0,0.0,0.0,0.0,-97.1438825386,0.0,0.0,0.0,0.0, *0.0,22.7027546667,0.0,0.0,0.0,-341.9933632553,0.0,0.0, *-111.8520186501,0.0,0.0,0.0,0.0,-48.2716070600,0.0,0.0,0.0,0.0, *-5.0447226741,0.0,0.0,0.0,0.0,0.0,0.0,144.1528855579/ THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/toang S1 = R1 - Req S2 = DCOS(THETA) - DCOS(THETAeq) S3 = R2 - Req Y1 = (S1 + S3)/2.0d0 Y2 = S2 Y3 = (S1 - S3)/2.0d0 c Potential in Y1,Y2,Y3. Full set of c derivatives for each order (1-4) and upper terms for others. DO 88 I5=1,NV 88 FT(I5)=0.0D0 FT(1)=1.0d0 FT(2)=Y1 FT(3)=Y2 FT(4)=Y1**2 FT(5)=Y2**2 FT(6)=Y3**2 FT(7)=Y1*Y2 FT(8)=Y1**3 FT(9)=Y2**3 FT(10)=Y1**2*Y2 FT(11)=Y2**2*Y1 FT(12)=Y3**2*Y1 FT(13)=Y3**2*Y2 FT(14)=Y1**4 FT(15)=Y2**4 FT(16)=Y3**4 FT(17)=Y1**3*Y2 FT(18)=Y2**3*Y1 FT(19)=Y1**2*Y2**2 FT(20)=Y1**2*Y3**2 FT(21)=Y2**2*Y3**2 FT(22)=Y3**2*Y1*Y2 FT(23)=Y1**5 FT(24)=Y2**5 FT(25)=Y1**4*Y2 FT(26)=Y2**4*Y1 FT(27)=Y3**4*Y1 FT(28)=Y3**4*Y2 FT(29)=Y1**3*Y2**2 FT(30)=Y1**3*Y3**2 FT(31)=Y2**3*Y1**2 FT(32)=Y2**3*Y3**2 FT(33)=Y1**2*Y2*Y3**2 FT(34)=Y1*Y2**2*Y3**2 FT(35)=Y1**6 FT(36)=Y2**6 FT(37)=Y3**6 FT(38)=Y1**5*Y2 FT(39)=Y2**5*Y1 FT(40)=Y1**4*Y2**2 FT(41)=Y2**4*Y1**2 FT(42)=Y2**4*Y3**2 FT(43)=Y3**4*Y2**2 FT(44)=Y1**4*Y3**2 FT(45)=Y3**4*Y1**2 FT(46)=Y3**4*Y1*Y2 FT(47)=Y1**3*Y2**3 FT(48)=Y2**3*Y1**2*Y3**2 FT(49)=Y1**3*Y3**2*Y2 FT(50)=Y2**3*Y3**2*Y1 FT(51)=Y1**2*Y2**2*Y1**2 FT(52)=Y1**7 FT(53)=Y2**7 FT(54)=Y1**6*Y2 FT(55)=Y2**6*Y1 FT(56)=Y3**6*Y1 FT(57)=Y3**6*Y2 FT(58)=Y1**5*Y2**2 FT(59)=Y1**5*Y3**2 FT(60)=Y2**5*Y1**2 FT(61)=Y2**5*Y3**2 FT(62)=Y1**4*Y2*Y3**2 FT(63)=Y1**4*Y2**3 FT(64)=Y2**4*Y1*Y3**2 FT(65)=Y2**4*Y1**3 FT(66)=Y3**4*Y1*Y2**2 FT(67)=Y3**4*Y2*Y1**2 FT(68)=Y3**4*Y1**3 FT(69)=Y3**4*Y2**3 FT(70)=Y1**3*Y3**2*Y2**2 FT(71)=Y1**8 FT(72)=Y2**8 FT(73)=Y3**8 FT(74)=Y1**7*Y2 FT(75)=Y2**7*Y1 FT(76)=Y1**6*Y2**2 FT(77)=Y1**6*Y3**2 FT(78)=Y2**6*Y1**2 FT(79)=Y2**6*Y3**2 FT(80)=Y3**6*Y1**2 FT(81)=Y3**6*Y2**2 FT(82)=Y3**6*Y2*Y1 FT(83)=Y1**5*Y3**2*Y2 FT(84)=Y1**5*Y2**3 FT(85)=Y2**5*Y1**3 FT(86)=Y2**5*Y3**2*Y1 FT(87)=Y1**4*Y2**4 FT(88)=Y1**4*Y3**4 FT(89)=Y1**4*Y2**2*Y3**2 FT(90)=Y2**4*Y3**4 FT(91)=Y2**4*Y3**2*Y1**2 FT(92)=Y3**4*Y1**2*Y2**2 FT(93)=Y3**4*Y1*Y2**3 FT(94)=Y3**4*Y1**3*Y2 FT(95)=Y1**3*Y2**3*Y3**2 FT(96)=Y2**9 FT(97)=Y2**9*Y1*Y3**4 FT(98)=Y2**9*Y1**5 FT(99)=Y2**10 FT(100)=Y2**10*Y1**2*Y3**2 FT(101)=Y2**10*Y3**4 FT(102)=Y2**10*Y1**4 FT(103)=Y2**11 FT(104)=Y2**11*Y1*Y3**2 FT(105)=Y2**11*Y1**3 FT(106)=Y2**12 FT(107)=Y2**12*Y1**2 FT(108)=Y2**12*Y3**2 FT(109)=Y2**13 FT(110)=Y2**13*Y1 FT(111)=Y2**14 FT(112)=Y2**8*Y1**6 FT(113)=Y2**8*Y3**6 FT(114)=Y2**8*Y1**4*Y3**2 FT(115)=Y2**8*Y3**4*Y1**2 FT(116)=Y2**7*Y1**7 FT(117)=Y2**7*Y1**5*Y3**2 V=ZERO DO 40 I=1,NV 40 V=V+CV(I)*FT(I) C SCALE AND SHIFT THE ZERO V=C0+SCALE*V RETURN END c Subroutine for calculation of MVD1 for H2-16O SUBROUTINE rel(V,R1,R2,THETA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FT(200), cv(200) DATA C0/0.0D0/,SCALE/1.0D-6/ data nv/200/ DATA ZERO/0.0D0/ DATA RZ/0.9576257 /,RHO/75.48992362 / DATA TOANG/0.5291772/ DATA cv/6.7010435494,24.4217638356,0.1926915998, *1.3694841332,1.7989449050,0.0,0.0,0.0,0.0,0.0,-0.0519929240, *-0.0000574005,-0.0001860844,-0.0006152872,0.0000879624, *-0.0004932809,-0.0000565213,0.0006886853,-0.0001027470, *0.0003081795,0.0001078370,0.0019925277,0.0001286772, *-0.0004511478,0.0000072157,-0.0004439080,-0.0001416903, *0.0000984644,-0.0003092456,-0.0029118771,-0.0000562847, *0.0000377062,0.0002610717,0.0000367337,0.0,-0.0002565517, *0.0008262448,0.0,-0.0001036317,0.0025090209,0.0002105614, *-0.0000204982,0.0,-0.0000414402,-0.0000676532,0.0000918240, *-0.0000400594,0.0,-0.0006468368,0.0001397619,0.0005356017, *0.0001585601,-0.0001899817,-0.0015914516,-0.0002822918, *0.0,0.0004567782,0.0,0.0001064879,0.0,0.0,0.0,-0.0001250236, *0.0000077559,0.0007064063,0.0,0.0,0.0,0.0006102666,-0.0004987995, *0.0,-0.0001943536,-0.0002855510,0.0,-0.0002176976,0.0002410759, *-0.0001644075,-0.0001182853,0.0,0.0,-0.0000272857,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, *0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/ c Equilibrium THETAeq = (180.d0 - RHO)*3.141592654/180.d0 Req= RZ/toang S1 = R1 - Req S2 = DCOS(THETA) - DCOS(THETAeq) S3 = R2 - Req Y1 = (S1 + S3)/2.0d0 Y2 = S2 Y3 = (S1 - S3)/2.0d0 c Potential in Y1,Y2,Y3. Full set of c derivatives for each order (1-4) and upper terms for others. DO 88 I5=1,NV 88 FT(I5)=0.0D0 FT(11)=1.0d0 FT(12)=Y1 FT(13)=Y2 FT(14)=Y1**2 FT(15)=Y2**2 FT(16)=Y3**2 FT(17)=Y1*Y2 FT(18)=Y1**3 FT(19)=Y2**3 FT(20)=Y1**2*Y2 FT(21)=Y2**2*Y1 FT(22)=Y3**2*Y1 FT(23)=Y3**2*Y2 FT(24)=Y1**4 FT(25)=Y2**4 FT(26)=Y3**4 FT(27)=Y1**3*Y2 FT(28)=Y2**3*Y1 FT(29)=Y1**2*Y2**2 FT(30)=Y1**2*Y3**2 FT(31)=Y2**2*Y3**2 FT(32)=Y3**2*Y1*Y2 FT(33)=Y1**5 FT(34)=Y2**5 FT(35)=Y1**4*Y2 FT(36)=Y2**4*Y1 FT(37)=Y3**4*Y1 FT(38)=Y3**4*Y2 FT(39)=Y1**3*Y2**2 FT(40)=Y1**3*Y3**2 FT(41)=Y2**3*Y1**2 FT(42)=Y2**3*Y3**2 FT(43)=Y1**2*Y2*Y3**2 FT(44)=Y1*Y2**2*Y3**2 FT(45)=Y1**6 FT(46)=Y2**6 FT(47)=Y3**6 FT(48)=Y1**5*Y2 FT(49)=Y2**5*Y1 FT(50)=Y1**4*Y2**2 FT(51)=Y2**4*Y1**2 FT(52)=Y2**4*Y3**2 FT(53)=Y3**4*Y2**2 FT(54)=Y1**4*Y3**2 FT(55)=Y3**4*Y1**2 FT(56)=Y3**4*Y1*Y2 FT(57)=Y1**3*Y2**3 FT(58)=Y2**3*Y1**3 FT(59)=Y1**3*Y3**2*Y2 FT(60)=Y2**3*Y3**2*Y1 FT(61)=Y1**2*Y2**2*Y1**2 FT(62)=Y1**7 FT(63)=Y2**7 FT(64)=Y1**6*Y2 FT(65)=Y2**6*Y1 FT(66)=Y3**6*Y1 FT(67)=Y3**6*Y2 FT(68)=Y1**5*Y2**2 FT(69)=Y1**5*Y3**2 FT(70)=Y2**5*Y1**2 FT(71)=Y2**5*Y3**2 FT(72)=Y1**4*Y2*Y3**2 FT(73)=Y1**4*Y2**3 FT(74)=Y2**4*Y1*Y3**2 FT(75)=Y2**4*Y1**3 FT(76)=Y3**4*Y1*Y2**2 FT(77)=Y3**4*Y2*Y1**2 FT(78)=Y3**4*Y1**3 FT(79)=Y3**4*Y2**3 FT(80)=Y1**3*Y3**2*Y2**2 FT(81)=Y1**8 FT(82)=Y2**8 FT(83)=Y3**8 FT(84)=Y1**7*Y2 FT(85)=Y2**7*Y1 FT(86)=Y1**6*Y2**2 FT(87)=Y1**6*Y3**2 FT(88)=Y2**6*Y1**2 FT(89)=Y2**6*Y3**2 FT(90)=Y3**6*Y1**2 FT(91)=Y3**6*Y2**2 FT(92)=Y3**6*Y2*Y1 FT(93)=Y1**5*Y3**2*Y2 FT(94)=Y1**5*Y2**3 FT(95)=Y2**5*Y1**3 FT(96)=Y2**5*Y3**2*Y1 FT(97)=Y1**4*Y2**4 FT(98)=Y1**4*Y3**4 FT(99)=Y1**4*Y2**2*Y3**2 FT(100)=Y2**4*Y3**4 FT(101)=Y2**4*Y3**2*Y1**2 FT(102)=Y3**4*Y1**2*Y2**2 FT(103)=Y3**4*Y1*Y2**3 FT(104)=Y3**4*Y1**3*Y2 FT(105)=Y1**3*Y2**3*Y3**2 FT(106)=Y2**9 FT(107)=Y2**9*Y1*Y3**4 FT(108)=Y2**9*Y1**5 FT(109)=Y2**10 FT(110)=Y2**10*Y1**2*Y3**2 FT(111)=Y2**10*Y3**4 FT(112)=Y2**10*Y1**4 FT(113)=Y2**11 FT(114)=Y2**11*Y1*Y3**2 FT(115)=Y2**11*Y1**3 FT(116)=Y2**12 FT(117)=Y2**12*Y1**2 FT(118)=Y2**12*Y3**2 FT(119)=Y2**13 FT(120)=Y2**13*Y1 FT(121)=Y2**14 FT(122)=Y2**8*Y1**6 FT(123)=Y2**8*Y3**6 FT(124)=Y2**8*Y1**4*Y3**2 FT(125)=Y2**8*Y3**4*Y1**2 FT(126)=Y2**7*Y1**7 FT(127)=Y2**7*Y1**5*Y3**2 FT(128)=Y2**7*Y1**3*Y3**4 FT(129)=Y2**7*Y1*Y3**6 FT(130)=Y2**6*Y1**8 FT(131)=Y2**6*Y3**8 FT(132)=Y2**6*Y1**6*Y3**2 FT(133)=Y2**6*Y1**4*Y3**4 FT(134)=Y2**6*Y1**2*Y3**6 FT(135)=Y2**5*Y1**9 FT(136)=Y2**5*Y1**7*Y3**2 FT(137)=Y2**5*Y1**5*Y3**4 FT(138)=Y2**5*Y1**3*Y3**6 FT(139)=Y2**5*Y1*Y3**8 FT(140)=Y2**4*Y1**10 FT(141)=Y2**4*Y3**10 FT(142)=Y2**4*Y1**8*Y3**2 FT(143)=Y2**4*Y1**6*Y3**4 FT(144)=Y2**4*Y1**4*Y3**6 FT(145)=Y2**4*Y1**2*Y3**8 FT(146)=Y2**3*Y1**11 FT(147)=Y2**3*Y1**9*Y3**2 FT(148)=Y2**3*Y1**7*Y3**4 FT(149)=Y2**3*Y1**5*Y3**6 FT(150)=Y2**3*Y1**3*Y3**8 FT(151)=Y2**3*Y1*Y3**10 FT(152)=Y2**2*Y1**12 FT(153)=Y2**2*Y3**12 FT(154)=Y2**2*Y1**10*Y3**2 FT(155)=Y2**2*Y1**8*Y3**4 FT(156)=Y2**2*Y1**6*Y3**6 FT(157)=Y2**2*Y1**4*Y3**8 FT(158)=Y2**2*Y1**2*Y3**10 FT(159)=Y2*Y1**12 FT(160)=Y2*Y1**11*Y3**2 FT(161)=Y2*Y1**9*Y3**4 FT(162)=Y2*Y1**7*Y3**6 FT(163)=Y2*Y1**5*Y3**8 FT(164)=Y2*Y1**3*Y3**10 FT(164)=Y2*Y1*Y3**12 FT(165)=Y2**8*Y1 FT(166)=Y2**7*Y1**2 FT(167)=Y2**7*Y3**2 FT(168)=Y2**6*Y1**3 FT(169)=Y2**6*Y3**2*Y1 FT(170)=Y2**5*Y1**4 FT(171)=Y2**5*Y3**4 FT(172)=Y2**5*Y1**2*Y3**2 FT(173)=Y2**4*Y1**5 FT(174)=Y2**4*Y1**3*Y3**2 FT(175)=Y2**4*Y1*Y3**4 FT(176)=Y2**3*Y1**6 FT(177)=Y2**3*Y3**6 FT(178)=Y2**3*Y1**4*Y3**2 FT(179)=Y2**3*Y1**2*Y3**4 FT(180)=Y2**2*Y1**7 FT(181)=Y2**2*Y1**5*Y3**2 FT(182)=Y2**2*Y1**3*Y3**4 FT(183)=Y2**2*Y1*Y3**6 FT(184)=Y2*Y1**8 FT(185)=Y2*Y3**8 FT(186)=Y2*Y1**6*Y3**2 FT(187)=Y2*Y1**4*Y3**4 FT(188)=Y2*Y1**2*Y3**6 FT(189)=Y2**9*Y1 FT(190)=Y2**8*Y1**2 FT(191)=Y2**8*Y3**2 FT(192)=Y2**7*Y1**3 FT(193)=Y2**7*Y3**2*Y1 FT(194)=Y2**6*Y1**4 FT(195)=Y2**6*Y3**4 FT(196)=Y2**6*Y1**2*Y3**2 FT(197)=Y2**5*Y1**5 FT(198)=Y2**5*Y1*Y3**4 FT(199)=Y2**5*Y1**3*Y3**2 FT(200)=Y2**4*Y1**6 V=ZERO DO 40 I=12,NV 40 V=V+CV(I)*FT(I) RETURN END SUBROUTINE BREITB3lin(VR,R1,R2,Z) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA SCALE/1.0D-6/ c This subroutine contains corrections to the water NBO PES due to the BREIT c term. c [see for istance HM Quiney et al Chem. Phys. Lett. 290 (1998) 473 , c Bethe and Salpheter, "Quantum mechanics of one and two-electron atoms"] c The corrections have been computed on a grid based on the 325 point grid c from P&S (see J. chem. phys., 106 (1997) 4618). c Moreover, a few extra points have been added , as well as a cut in the radial c coordinates (lines 1-2), in order to account for high bending modes. c The final grid used contains 293 points. c Then the points have been fitted with a polynomial in X,Y and Z, using a c Mathematica script. c The basis set used for the electronic calculations is the set called 'B' c provided by H.M. Quiney (see Chem. Phys. Lett. ). c c Those corrections have been computed by Paolo c email: paolo@theory.phys.ucl.ac.uk . c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c The inputs are in u.a., and X,Y are the distances of the H atoms from c the oxygen, and Z is the angle HOH in radiants. The final result is in c Hartree. c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c Potential symmetrical in X,Y,Z. Full set of c derivatives for each order (1-4) and upper terms for others Seq=1.8240445d0 e=1.809645d0 h=0.5d0 x=r1 y=r2 1 if (x.gt.3.5d0) x=3.5d0 2 if (y.gt.3.5d0) y=3.5d0 xm=x-y xp=x+y xp2=-e+h*xp zm=-Seq+z v1=7577.371974812879+54.16425693884002*xm**2+ y 6.774810680339735*xm**4-0.7125792752605449*xm**6- y 212.3502737780579*xp2- y 119.8400228527018*xm**2*xp2-29.18501086262345*xm**4*xp2+ y 221.1271759652365*xp2**2+72.97844564936518*xm**2*xp2**2+ y 89.0829000401668*xm**4*xp2**2-180.6736318109525*xp2**3+ y 321.1679099895901*xm**2*xp2**3-86.9772033725647*xm**4*xp2**3+ y 81.467133012329*xp2**4-890.357807854604*xm**2*xp2**4+ y 57.52112099193129*xp2**5+676.0925261740228*xm**2*xp2**5- y 58.228375665496*xp2**6+1.323061103513093*zm- y 0.5270121336782832*xm**2*zm+0.01985434592556156*xp2* y zm+7.715533877300036*xm**2*xp2*zm-6.229144779356836*xm**4*xp2* y zm-0.1710351457311522*xp2**2*zm+9.92240371274571*xm**4*xp2**2* y zm+10.47917821735698*xp2**3*zm-27.67291046310705*xm**2*xp2**4* y zm+15.42094531691969*xp2**6*zm+7.749223096520598*zm**2- y 0.2147443261198608*xm**2*zm**2-1.164069372530965*xm**4*zm**2+ y 10.24148015693046*xp2*zm**2-2.364801830233726*xm**2*xp2*zm**2 v2=3.405435470980123*xm**4*xp2*zm**2+11.53659470986242*xp2**2* y zm**2+40.71096970108562*xp2**3*zm**2-65.49114275587444*xp2**4* y zm**2+0.5246601257334035*zm**3+1.025008298074623*xm**2*zm**3+ y 13.57824254391274*xp2*zm**3-7.469419914001589*xp2**2* y zm**3-33.70757112970705*xp2**3*zm**3+30.20514216972833*xp2**4* y zm**3-10.53913543447923*zm**4-0.6159136295163627*xm**2*zm**4- y 19.56431274355461*xp2*zm**4-20.81965238209867*xp2**2*zm**4+ y 5.998958874987758*xp2**3*zm**4-9.44711265431818*zm**5- y 22.55622148750276*xp2*zm**5+16.30440168684215*xp2**2* y zm**5+19.20675957512514*zm**6+19.78080962673524*xp2* y zm**6+8.08923849773384*zm**7-10.68490632273025*zm**8 VR=V1+V2 C SCALE AND SHIFT THE ZERO VR=VR*1.0d-6 RETURN END SUBROUTINE PESd2x(VR,x,y,z) IMPLICIT REAL*8 (A-H,O-Z) c This subroutine contains corrections to the water NBO PES due to the Darwin c 2 electrons term. Those corrections have been computed by Gyorgy c tarczay@para.chem.elte.hu c Those corrections are computed on the P&S grid of 325 points. c (see J. chem. phys., 106 (1997) 4618), to which a few x points have been c added in irder to account to high bending modes. The final grid containd 341 c points. c The input are in u.a., and X,Y are the distances of the H atoms from c the oxygen, and Z is the angle HOH in radiants. The final result is in c Hartree. Seq=1.8240445d0 pi = dacos(-1.0d0) e=1.809645d0 h=0.5d0 x1=x y1=y 1 if (x.gt.3.5d0) x1=3.5d0 2 if (y.gt.3.5d0) y1=3.5d0 xm=x1-y1 xp=x1+y1 xp2=-e+h*xp zm=-Seq+z c Potential symmetrical in S1,S2,S3. Full set of c derivatives for each order (1-4) and upper terms for others C v1=-3263.067522028298-9.17958228916955*xm**2- y 0.5165611032611932*xm**4-0.5157212949525876*xm**6+ y 18.31161578215203*xp2+ y 30.14193751791963*xm**2*xp2-13.62543868575853*xm**4*xp2- y 43.37232019119388*xp2**2-50.50353364079294*xm**2*xp2**2+ y 70.36441193443143*xm**4*xp2**2+27.43935454999898*xp2**3+ y 123.751990625258*xm**2*xp2**3-76.80240321256033*xm**4*xp2**3- y 9.50017804016001*xp2**4-363.4487347625543*xm**2*xp2**4+ y 113.1940248029751*xp2**5+376.6560011408163*xm**2*xp2**5- y 164.6523756673548*xp2**6+9.16256842998227*zm y -1.22230095639504*xm**2*zm+1.33032571356463*xp2* y zm-0.94822119654751*xm**2*xp2*zm+0.7645470802285307*xm**4*xp2* y zm-11.77270680473595*xp2**2*zm-0.4065994514809928*xm**4*xp2**2* y zm-2.113651214829342*xp2**3*zm-3.653921741665064*xm**2*xp2**4* y zm+26.53983199106825*xp2**6*zm+3.099164302936567*zm**2 v2=-0.4668245990549825*xm**2*zm**2+0.05845413180128318*xm**4*zm**2 y +2.708722250876111*xp2*zm**2-2.578482367020144*xm**2*xp2*zm**2- y 0.1605392233404811*xm**4*xp2*zm**2-10.57780429022803*xp2**2* y zm**2-3.496293826717189*xp2**3*zm**2+23.46280699747645*xp2**4* y zm**2+1.8547816858377*zm**3-0.4003662844685243*xm**2*zm**3+ y 3.040229985315839*xp2*zm**3-4.955739113923876*xp2**2* y zm**3+14.05364889791468*xp2**3*zm**3-21.6926320924828*xp2**4* y zm**3-1.321464834042384*zm**4+2.298844571392118*xm**2*zm**4- y 2.633405421645483*xp2*zm**4+20.97178840867901*xp2**2* y zm**4-32.18658937476802*xp2**3*zm**4-0.5992225949734171*zm**5+ y 2.059827452250273*xp2*zm**5+0.6453850286056735*xp2**2* y zm**5-0.4620689505336259*zm**6+0.7465042626807512*xp2* y zm**6-0.1254018119377959*zm**7+0.01947721364782498*zm**8 VR=V1+V2 C SCALE AND SHIFT THE ZERO VR=VR*1.0d-6 RETURN END SUBROUTINE PESleq6(V1,x,y,z) IMPLICIT DOUBLE PRECISION (A-H,O-Z) V1=-3514.850376005703+1.189315215135138*(x-y)**2+ - 0.824157459531989*(x-y)**4+0.03853108456851828*(x-y)**6+ - 12.83265590340491*(-1.809659+0.5*(x+y))- - 9.51736455454466*(x-y)**2*(-1.809659+0.5*(x+y))- - 3.027576695974858*(x-y)**4*(-1.809659+0.5*(x+y))+ - 10.94033338777717*(-1.809659+0.5*(x+y))**2+ - 15.53332877554612*(x-y)**2*(-1.809659+0.5*(x+y))**2+ - 6.063907309056958*(x-y)**4*(-1.809659+0.5*(x+y))**2- - 13.79644533708051*(-1.809659+0.5*(x+y))**3- - 26.67549601926293*(x-y)**2*(-1.809659+0.5*(x+y))**3- - 6.200836894255189*(x-y)**4*(-1.809659+0.5*(x+y))**3+ - 5.688103460541242*(-1.809659+0.5*(x+y))**4+ - 52.9835500898771*(x-y)**2*(-1.809659+0.5*(x+y))**4- - 11.88910471926647*(-1.809659+0.5*(x+y))**5- - 43.99657824332825*(x-y)**2*(-1.809659+0.5*(x+y))**5+ - 15.874512160015*(-1.809659+0.5*(x+y))**6- - 8.60706112101134*(-1.824045+z)+ - 1.264485336667462*(x-y)**2*(-1.824045+z)- - 0.915127202947929*(-1.809659+0.5*(x+y))* - (-1.824045+z)+ - 0.6556566908758441*(x-y)**2*(-1.809659+0.5*(x+y))* - (-1.824045+z)- - 0.813078328219753*(x-y)**4*(-1.809659+0.5*(x+y))* - (-1.824045+z)+ - 12.42234678727481*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)+ - 2.805488560712774*(x-y)**4*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)- - 4.937250627623143*(-1.809659+0.5*(x+y))**3* - (-1.824045+z)- - 3.095201035295474*(x-y)**2*(-1.809659+0.5*(x+y))**4* - (-1.824045+z)+ - 18.85309150691318*(-1.809659+0.5*(x+y))**6* - (-1.824045+z)- - 3.209703208057476*(-1.824045+z)**2+ - 0.5360421552708203*(x-y)**2*(-1.824045+z)**2- - 0.263467844585989*(x-y)**4*(-1.824045+z)**2- - 1.13298516075929*(-1.809659+0.5*(x+y))* - (-1.824045+z)**2- - 0.06909229322445753*(x-y)**2*(-1.809659+0.5*(x+y))* - (-1.824045+z)**2+ - 1.649063526503709*(x-y)**4*(-1.809659+0.5*(x+y))* - (-1.824045+z)**2+ - 3.603611347474725*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)**2+ - 3.757418764813337*(-1.809659+0.5*(x+y))**3* - (-1.824045+z)**2+ - 4.607672502246032*(-1.809659+0.5*(x+y))**4* - (-1.824045+z)**2- - 0.7490414640610651*(-1.824045+z)**3- - 0.0888181500794012*(x-y)**2*(-1.824045+z)**3- - 5.334303151299991*(-1.809659+0.5*(x+y))* - (-1.824045+z)**3+ - 1.37948603262339*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)**3+ - 11.24395154910416*(-1.809659+0.5*(x+y))**3* - (-1.824045+z)**3 - - 17.85690001161674*(-1.809659+0.5*(x+y))**4* - (-1.824045+z)**3+ - 0.7694433624551493*(-1.824045+z)**4- - 0.939662303404418*(x-y)**2*(-1.824045+z)**4- - 2.296000209594694*(-1.809659+0.5*(x+y))* - (-1.824045+z)**4- - 4.514249057965571*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)**4- - 2.324765391545952*(-1.809659+0.5*(x+y))**3* - (-1.824045+z)**4+ - 0.223711667169141*(-1.824045+z)**5+ - 1.164515013150094*(-1.809659+0.5*(x+y))* - (-1.824045+z)**5- - 2.825913168656484*(-1.809659+0.5*(x+y))**2* - (-1.824045+z)**5+ - 0.4811142779617512*(-1.824045+z)**6+ - 1.292817090808966*(-1.809659+0.5*(x+y))* - (-1.824045+z)**6+ - 0.1657130839026308*(-1.824045+z)**7- - 0.02192338698614548*(-1.824045+z)**8 v1=v1/1000000.0d0 RETURN END SUBROUTINE POTS(V,rij1,rij2,rij3) implicit double precision (a-h,o-y),logical(z) c c pes for h2o, c Harry Partridge and David W. Schwenke, J. Chem. Phys., c submitted Nov. 8, 1996. c rij(i,1)& rij(i,2) are oh distances in au c rij(i,3) is hoh angle in rad c v(i) is pes in au c n is number of geometries c mass dependent factors are included. the nuclear masses c should be passed to this program using the array xm in c common potmcm. xm(1) is the c mass of the hydrogen associated with rij(i,1), and xm(2) c is the mass of the hydrogen associated with rij(i,2). c all masses are in au. c dimension c5z(245),cbasis(245),ccore(245), $ crest(245),idx(245,3),fmat(15,3),cmass(9),idxm(9,3) c c expansion indicies c data (idx(i,1),i=1,245)/ $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, $ 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, $ 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, $ 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4, $ 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, $ 6, 6, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, $ 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 5, 5, $ 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, $ 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, $ 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 9, 9, $ 9, 9, 9, 9, 9/ data (idx(i,2),i=1,245)/ $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, $ 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, $ 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, $ 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, $ 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, $ 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, $ 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, $ 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, $ 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, $ 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, $ 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 1, 1, $ 1, 1, 1, 1, 1/ data (idx(i,3),i=1,245)/ $ 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15, 1, 2, 3, 4, 5, $ 6, 7, 8, 9,10,11,12,13,14, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, $12,13, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13, 1, 2, 3, 4, 5, $ 6, 7, 8, 9,10,11,12, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12, 1, $ 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, $11, 1, 2, 3, 4, 5, 6, 7, 8, 9,10,11, 1, 2, 3, 4, 5, 6, 7, 8, $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, 1, 2, 3, 4, 5, 6, 7, 8, $ 9,10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, $ 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, $ 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, $ 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, $ 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7, 1, 2, $ 3, 4, 5, 6, 7/ c c expansion coefficients for 5z ab initio data c data (c5z(i),i=1,45)/ $ 4.2278462684916D+04, 4.5859382909906D-02, 9.4804986183058D+03, $ 7.5485566680955D+02, 1.9865052511496D+03, 4.3768071560862D+02, $ 1.4466054104131D+03, 1.3591924557890D+02,-1.4299027252645D+03, $ 6.6966329416373D+02, 3.8065088734195D+03,-5.0582552618154D+02, $-3.2067534385604D+03, 6.9673382568135D+02, 1.6789085874578D+03, $-3.5387509130093D+03,-1.2902326455736D+04,-6.4271125232353D+03, $-6.9346876863641D+03,-4.9765266152649D+02,-3.4380943579627D+03, $ 3.9925274973255D+03,-1.2703668547457D+04,-1.5831591056092D+04, $ 2.9431777405339D+04, 2.5071411925779D+04,-4.8518811956397D+04, $-1.4430705306580D+04, 2.5844109323395D+04,-2.3371683301770D+03, $ 1.2333872678202D+04, 6.6525207018832D+03,-2.0884209672231D+03, $-6.3008463062877D+03, 4.2548148298119D+04, 2.1561445953347D+04, $-1.5517277060400D+05, 2.9277086555691D+04, 2.6154026873478D+05, $-1.3093666159230D+05,-1.6260425387088D+05, 1.2311652217133D+05, $-5.1764697159603D+04, 2.5287599662992D+03, 3.0114701659513D+04/ data (c5z(i),i=46,90)/ $-2.0580084492150D+03, 3.3617940269402D+04, 1.3503379582016D+04, $-1.0401149481887D+05,-6.3248258344140D+04, 2.4576697811922D+05, $ 8.9685253338525D+04,-2.3910076031416D+05,-6.5265145723160D+04, $ 8.9184290973880D+04,-8.0850272976101D+03,-3.1054961140464D+04, $-1.3684354599285D+04, 9.3754012976495D+03,-7.4676475789329D+04, $-1.8122270942076D+05, 2.6987309391410D+05, 4.0582251904706D+05, $-4.7103517814752D+05,-3.6115503974010D+05, 3.2284775325099D+05, $ 1.3264691929787D+04, 1.8025253924335D+05,-1.2235925565102D+04, $-9.1363898120735D+03,-4.1294242946858D+04,-3.4995730900098D+04, $ 3.1769893347165D+05, 2.8395605362570D+05,-1.0784536354219D+06, $-5.9451106980882D+05, 1.5215430060937D+06, 4.5943167339298D+05, $-7.9957883936866D+05,-9.2432840622294D+04, 5.5825423140341D+03, $ 3.0673594098716D+03, 8.7439532014842D+04, 1.9113438435651D+05, $-3.4306742659939D+05,-3.0711488132651D+05, 6.2118702580693D+05, $-1.5805976377422D+04,-4.2038045404190D+05, 3.4847108834282D+05/ data (c5z(i),i=91,135)/ $-1.3486811106770D+04, 3.1256632170871D+04, 5.3344700235019D+03, $ 2.6384242145376D+04, 1.2917121516510D+05,-1.3160848301195D+05, $-4.5853998051192D+05, 3.5760105069089D+05, 6.4570143281747D+05, $-3.6980075904167D+05,-3.2941029518332D+05,-3.5042507366553D+05, $ 2.1513919629391D+03, 6.3403845616538D+04, 6.2152822008047D+04, $-4.8805335375295D+05,-6.3261951398766D+05, 1.8433340786742D+06, $ 1.4650263449690D+06,-2.9204939728308D+06,-1.1011338105757D+06, $ 1.7270664922758D+06, 3.4925947462024D+05,-1.9526251371308D+04, $-3.2271030511683D+04,-3.7601575719875D+05, 1.8295007005531D+05, $ 1.5005699079799D+06,-1.2350076538617D+06,-1.8221938812193D+06, $ 1.5438780841786D+06,-3.2729150692367D+03, 1.0546285883943D+04, $-4.7118461673723D+04,-1.1458551385925D+05, 2.7704588008958D+05, $ 7.4145816862032D+05,-6.6864945408289D+05,-1.6992324545166D+06, $ 6.7487333473248D+05, 1.4361670430046D+06,-2.0837555267331D+05, $ 4.7678355561019D+05,-1.5194821786066D+04,-1.1987249931134D+05/ data (c5z(i),i=136,180)/ $ 1.3007675671713D+05, 9.6641544907323D+05,-5.3379849922258D+05, $-2.4303858824867D+06, 1.5261649025605D+06, 2.0186755858342D+06, $-1.6429544469130D+06,-1.7921520714752D+04, 1.4125624734639D+04, $-2.5345006031695D+04, 1.7853375909076D+05,-5.4318156343922D+04, $-3.6889685715963D+05, 4.2449670705837D+05, 3.5020329799394D+05, $ 9.3825886484788D+03,-8.0012127425648D+05, 9.8554789856472D+04, $ 4.9210554266522D+05,-6.4038493953446D+05,-2.8398085766046D+06, $ 2.1390360019254D+06, 6.3452935017176D+06,-2.3677386290925D+06, $-3.9697874352050D+06,-1.9490691547041D+04, 4.4213579019433D+04, $ 1.6113884156437D+05,-7.1247665213713D+05,-1.1808376404616D+06, $ 3.0815171952564D+06, 1.3519809705593D+06,-3.4457898745450D+06, $ 2.0705775494050D+05,-4.3778169926622D+05, 8.7041260169714D+03, $ 1.8982512628535D+05,-2.9708215504578D+05,-8.8213012222074D+05, $ 8.6031109049755D+05, 1.0968800857081D+06,-1.0114716732602D+06, $ 1.9367263614108D+05, 2.8678295007137D+05,-9.4347729862989D+04/ data (c5z(i),i=181,225)/ $ 4.4154039394108D+04, 5.3686756196439D+05, 1.7254041770855D+05, $-2.5310674462399D+06,-2.0381171865455D+06, 3.3780796258176D+06, $ 7.8836220768478D+05,-1.5307728782887D+05,-3.7573362053757D+05, $ 1.0124501604626D+06, 2.0929686545723D+06,-5.7305706586465D+06, $-2.6200352535413D+06, 7.1543745536691D+06,-1.9733601879064D+04, $ 8.5273008477607D+04, 6.1062454495045D+04,-2.2642508675984D+05, $ 2.4581653864150D+05,-9.0376851105383D+05,-4.4367930945690D+05, $ 1.5740351463593D+06, 2.4563041445249D+05,-3.4697646046367D+03, $-2.1391370322552D+05, 4.2358948404842D+05, 5.6270081955003D+05, $-8.5007851251980D+05,-6.1182429537130D+05, 5.6690751824341D+05, $-3.5617502919487D+05,-8.1875263381402D+02,-2.4506258140060D+05, $ 2.5830513731509D+05, 6.0646114465433D+05,-6.9676584616955D+05, $ 5.1937406389690D+05, 1.7261913546007D+05,-1.7405787307472D+04, $-3.8301842660567D+05, 5.4227693205154D+05, 2.5442083515211D+06, $-1.1837755702370D+06,-1.9381959088092D+06,-4.0642141553575D+05/ data (c5z(i),i=226,245)/ $ 1.1840693827934D+04,-1.5334500255967D+05, 4.9098619510989D+05, $ 6.1688992640977D+05, 2.2351144690009D+05,-1.8550462739570D+06, $ 9.6815110649918D+03,-8.1526584681055D+04,-8.0810433155289D+04, $ 3.4520506615177D+05, 2.5509863381419D+05,-1.3331224992157D+05, $-4.3119301071653D+05,-5.9818343115856D+04, 1.7863692414573D+03, $ 8.9440694919836D+04,-2.5558967650731D+05,-2.2130423988459D+04, $ 4.4973674518316D+05,-2.2094939343618D+05/ c c expansion coefficients for basis correction c data (cbasis(i),i=1,45)/ $ 6.9770019624764D-04,-2.4209870001642D+01, 1.8113927151562D+01, $ 3.5107416275981D+01,-5.4600021126735D+00,-4.8731149608386D+01, $ 3.6007189184766D+01, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $-7.7178474355102D+01,-3.8460795013977D+01,-4.6622480912340D+01, $ 5.5684951167513D+01, 1.2274939911242D+02,-1.4325154752086D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00,-6.0800589055949D+00, $ 8.6171499453475D+01,-8.4066835441327D+01,-5.8228085624620D+01, $ 2.0237393793875D+02, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 3.3525582670313D+02, 7.0056962392208D+01,-4.5312502936708D+01/ data (cbasis(i),i=46,90)/ $-3.0441141194247D+02, 2.8111438108965D+02, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-1.2983583774779D+02, 3.9781671212935D+01, $-6.6793945229609D+01,-1.9259805675433D+02, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-8.2855757669957D+02,-5.7003072730941D+01, $-3.5604806670066D+01, 9.6277766002709D+01, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 8.8645622149112D+02,-7.6908409772041D+01, $ 6.8111763314154D+01, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (cbasis(i),i=91,135)/ $ 2.5090493428062D+02,-2.3622141780572D+02, 5.8155647658455D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 2.8919570295095D+03, $-1.7871014635921D+02,-1.3515667622500D+02, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-3.6965613754734D+03, 2.1148158286617D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00,-1.4795670139431D+03, $ 3.6210798138768D+02, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $-5.3552886800881D+03, 3.1006384016202D+02, 0.0000000000000D+00/ data (cbasis(i),i=136,180)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 1.6241824368764D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 4.3764909606382D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 1.0940849243716D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 3.0743267832931D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (cbasis(i),i=181,225)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (cbasis(i),i=226,245)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00/ c c expansion coefficients for core correction c data (ccore(i),i=1,45)/ $ 2.4332191647159D-02,-2.9749090113656D+01, 1.8638980892831D+01, $-6.1272361746520D+00, 2.1567487597605D+00,-1.5552044084945D+01, $ 8.9752150543954D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $-3.5693557878741D+02,-3.0398393196894D+00,-6.5936553294576D+00, $ 1.6056619388911D+01, 7.8061422868204D+01,-8.6270891686359D+01, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00,-3.1688002530217D+01, $ 3.7586725583944D+01,-3.2725765966657D+01,-5.6458213299259D+00, $ 2.1502613314595D+01, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 5.2789943583277D+02,-4.2461079404962D+00,-2.4937638543122D+01/ data (ccore(i),i=46,90)/ $-1.1963809321312D+02, 2.0240663228078D+02, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-6.2574211352272D+02,-6.9617539465382D+00, $-5.9440243471241D+01, 1.4944220180218D+01, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-1.2851139918332D+03,-6.5043516710835D+00, $ 4.0410829440249D+01,-6.7162452402027D+01, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 1.0031942127832D+03, 7.6137226541944D+01, $-2.7279242226902D+01, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (ccore(i),i=91,135)/ $-3.3059000871075D+01, 2.4384498749480D+01,-1.4597931874215D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 1.6559579606045D+03, $ 1.5038996611400D+02,-7.3865347730818D+01, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-1.9738401290808D+03,-1.4149993809415D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00,-1.2756627454888D+02, $ 4.1487702227579D+01, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $-1.7406770966429D+03,-9.3812204399266D+01, 0.0000000000000D+00/ data (ccore(i),i=136,180)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-1.1890301282216D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 2.3723447727360D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-1.0279968223292D+03, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 5.7153838472603D+02, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (ccore(i),i=181,225)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (ccore(i),i=226,245)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00/ c c expansion coefficients for v rest c data (crest(i),i=1,45)/ $ 0.0000000000000D+00,-4.7430930170000D+00,-1.4422132560000D+01, $-1.8061146510000D+01, 7.5186735000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $-2.7962099800000D+02, 1.7616414260000D+01,-9.9741392630000D+01, $ 7.1402447000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00,-7.8571336480000D+01, $ 5.2434353250000D+01, 7.7696745000000D+01, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 1.7799123760000D+02, 1.4564532380000D+02, 2.2347226000000D+02/ data (crest(i),i=46,90)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-4.3823284100000D+02,-7.2846553000000D+02, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00,-2.6752313750000D+02, 3.6170310000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (crest(i),i=91,135)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (crest(i),i=136,180)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (crest(i),i=181,225)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00/ data (crest(i),i=226,245)/ $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00, 0.0000000000000D+00, $ 0.0000000000000D+00, 0.0000000000000D+00/ c c expansion indicies for mass correction c data idxm/1,2,1,1,3,2,1,2,1, $ 2,1,1,3,1,2,2,1,1, $ 1,1,2,1,1,1,2,2,3/ c c expansion coefficients for mass correction c data cmass/ -8.3554183D+00,3.7036552D+01,-5.2722136D+00, $ 1.6843857D+01,-7.0929741D+01,5.5380337D+00,-2.9962997D+01, $ 1.3637682D+02,-3.0530195d+00/ c c two body parameters c data reoh,thetae,b1,roh,alphaoh,deoh,phh1,phh2/0.958649d0, $ 104.3475d0,2.0d0,0.9519607159623009d0,2.587949757553683d0, $ 42290.92019288289d0,16.94879431193463d0,12.66426998162947d0/ c c scaling factors for contributions to emperical potential c data f5z,fbasis,fcore,frest/0.0d0,0.0d0,-1.0d0,0.0d0/ save data ifirst/0/ if(ifirst.eq.0)then ifirst=1 1 format(/1x,'pes for h2o', $ /1x,'by Harry Partridge and David W. Schwenke', $ /1x,'submitted to J. Chem. Phys. Nov. 8, 1996') 56 format(/1x,'parameters before adjustment') 55 format(/1x,'two body potential parameters:', $ /1x,'hh: phh1 = ',f10.1,' phh2 = ',f5.2, $ /1x,'oh: deoh = ',f10.1,' alpha = ',f7.4, $ ' re = ',f7.4) 4 format(/1x,'three body parameters:', $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, $ /1x,'betaoh = ',f10.4, $ /1x,' i j k',7x,'c5z',9x,'cbasis',10x,'ccore', $ 10x,'crest') do 2 i=1,245 2 continue c c remove mass correction from vrest c xmh=1836.152697d0 xmhi=1d0/xmh xmd=3670.483031d0 fact=1.d0/((1.d0/xmd)-(1.d0/xmh)) 65 format(/1x,'parameters for delta v hdo ', $ /1x,' i j k') do 60 i=1,9 cmass(i)=cmass(i)*fact corr=cmass(i)*xmhi if(idxm(i,1).eq.idxm(i,2))corr=corr*0.5d0 do 61 j=1,245 if(idx(j,1).eq.idxm(i,1).and.idx(j,2).eq.idxm(i,2).and. $ idx(j,3).eq.idxm(i,3))then crest(j)=crest(j)-corr go to 62 end if 61 continue 62 continue do 63 j=1,245 if(idx(j,2).eq.idxm(i,1).and.idx(j,1).eq.idxm(i,2).and. $ idx(j,3).eq.idxm(i,3))then crest(j)=crest(j)-corr go to 64 end if 63 continue 64 continue 60 continue xm1=1.d0/1836.152697d0 xm2=1.d0/1836.152697d0 c c adjust parameters using scale factors c 57 format(/1x,'adjusting parameters using scale factors ', $ /1x,'f5z = ',f11.8, $ /1x,'fbasis = ',f11.8, $ /1x,'fcore = ',f11.8, $ /1x,'frest = ',f11.8) phh1=phh1*f5z deoh=deoh*f5z do 59 i=1,245 c5z(i)=f5z*c5z(i)+fbasis*cbasis(i)+fcore*ccore(i) $ +frest*crest(i) 59 continue 58 format(/1x,'three body parameters:', $ /1x,'reoh = ',f10.4,' thetae = ',f10.4, $ /1x,'betaoh = ',f10.4, $ /1x,' i j k cijk', $ /(1x,3i5,1pe15.7)) do 66 i=1,9 cmass(i)=cmass(i)*frest 66 continue 76 format(/1x,'mass correction factors ', $ /1x,' i j k cijk', $ /(1x,3i5,1pe15.7)) c c convert parameters from 1/cm, angstrom to a.u. c reoh=reoh/0.529177249d0 b1=b1*0.529177249d0*0.529177249d0 do 3 i=1,245 c5z(i)=c5z(i)*4.556335d-6 3 continue do 67 i=1,9 cmass(i)=cmass(i)*4.556335d-6 67 continue rad=acos(-1d0)/1.8d2 ce=cos(thetae*rad) phh1=phh1*exp(phh2) phh1=phh1*4.556335d-6 phh2=phh2*0.529177249d0 deoh=deoh*4.556335d-6 roh=roh/0.529177249d0 alphaoh=alphaoh*0.529177249d0 c5z(1)=c5z(1)*2d0 end if x1=(rij1-reoh)/reoh x2=(rij2-reoh)/reoh x3=cos(rij3)-ce rhh=sqrt(rij1**2+rij2**2 $ -2d0*rij1*rij2*cos(rij3)) vhh=phh1*exp(-phh2*rhh) ex=exp(-alphaoh*(rij1-roh)) voh1=deoh*ex*(ex-2d0) ex=exp(-alphaoh*(rij2-roh)) voh2=deoh*ex*(ex-2d0) fmat(1,1)=1d0 fmat(1,2)=1d0 fmat(1,3)=1d0 do 10 j=2,15 fmat(j,1)=fmat(j-1,1)*x1 fmat(j,2)=fmat(j-1,2)*x2 fmat(j,3)=fmat(j-1,3)*x3 10 continue v=0d0 do 12 j=2,245 term=c5z(j)*(fmat(idx(j,1),1)*fmat(idx(j,2),2) $ +fmat(idx(j,2),1)*fmat(idx(j,1),2)) $ *fmat(idx(j,3),3) v=v+term 12 continue v1=0d0 v2=0d0 do 13 j=1,9 v1=v1+cmass(j)*fmat(idxm(j,1),1)*fmat(idxm(j,2),2) $ *fmat(idxm(j,3),3) v2=v2+cmass(j)*fmat(idxm(j,2),1)*fmat(idxm(j,1),2) $ *fmat(idxm(j,3),3) 13 continue v=v+xm1*v1+xm2*v2 v=v*exp(-b1*((rij1-reoh)**2+(rij2-reoh)**2)) $ +c5z(1) $ +voh1+voh2+vhh return end SUBROUTINE wifin(v,r1,r2,th) c c In here almost all the points have been included, 342 in total c c implicit real*8 (a-h,o-z) reoh=0.958649d0 thetae=104.3475d0 b1=2.0d0 roh=0.951961d0 alphaoh=2.587949757553683d0 phh2=6.70164303995d0 t0=0.01d0 ut=20.d0 x0=2.5d0 ut2=20.d0 x02=2.5d0 thetae=thetae*.314159265358979312d01*.00555555555555555555d0 xs1=(r1+r2)*0.5d0-reoh xs2=(r1-r2)*0.5d0 xst=dcos(th)-dcos(thetae) rs=dsqrt(0.5d0)*(r1+r2) rm=dsqrt(0.5d0)*(r1-r2) rr1=r1-roh rr2=r2-roh xep1=dexp(-2.d0*alphaoh*rr1)-2.d0*dexp(-alphaoh*rr1)+1.d0 xep2=dexp(-2.d0*alphaoh*rr2)-2.d0*dexp(-alphaoh*rr2)+1.d0 xep3=dexp(-b1*(rr1**2+rr2**2)) rhh=dsqrt(r1**2+r2**2-2.d0*r1*r2*dcos(th)) vhh=0.900642240911285975d6*dexp(-phh2*rhh) vpb1=0.518622556959170834d5*xep1 vpb2=0.518622556959170834d5*xep2 v0=-0.69833729793860E+02*xs1**0*xs2**0*xst**0 vp1=-0.88931405097186E+02*xs1**0*xs2**0*xst**1 *-0.87065167087729E+04*xs1**1*xs2**0*xst**0 *+0.18561833211448E+05*xs1**0*xs2**0*xst**2 *-0.22925714045631E+06*xs1**0*xs2**2*xst**0 *-0.25473078391522E+05*xs1**1*xs2**0*xst**1 *-0.24221902469321E+06*xs1**2*xs2**0*xst**0 *+0.91104690191314E+03*xs1**0*xs2**0*xst**3 *-0.20846258713362E+05*xs1**0*xs2**2*xst**1 *-0.97259666765996E+04*xs1**1*xs2**0*xst**2 *+0.21895229646374E+07*xs1**1*xs2**2*xst**0 *+0.23321398548491E+05*xs1**2*xs2**0*xst**1 *+0.70135119326976E+06*xs1**3*xs2**0*xst**0 *+0.35115546540196E+04*xs1**0*xs2**0*xst**4 *+0.56216397977173E+05*xs1**0*xs2**2*xst**2 *-0.20038134529604E+07*xs1**0*xs2**4*xst**0 *-0.10382409164897E+05*xs1**1*xs2**0*xst**3 *-0.15791163967577E+05*xs1**1*xs2**2*xst**1 *+0.69131950382604E+05*xs1**2*xs2**0*xst**2 *-0.83355495109342E+07*xs1**2*xs2**2*xst**0 *-0.72835013274781E+05*xs1**3*xs2**0*xst**1 *-0.20429129369404E+07*xs1**4*xs2**0*xst**0 *+0.55792697169003E+03*xs1**0*xs2**0*xst**5 *+0.24482126719682E+04*xs1**0*xs2**2*xst**3 *-0.27399255063928E+05*xs1**0*xs2**4*xst**1 *+0.30011447505506E+04*xs1**1*xs2**0*xst**4 *-0.40810705479272E+05*xs1**1*xs2**2*xst**2 *+0.13148302968514E+08*xs1**1*xs2**4*xst**0 *-0.59362784536737E+04*xs1**2*xs2**0*xst**3 *-0.86061422326869E+05*xs1**2*xs2**2*xst**1 *-0.32156522858267E+05*xs1**3*xs2**0*xst**2 *+0.20293315578616E+08*xs1**3*xs2**2*xst**0 *+0.19710424146155E+05*xs1**4*xs2**0*xst**1 *+0.37828666711928E+07*xs1**5*xs2**0*xst**0 *+0.53407786571380E+03*xs1**0*xs2**0*xst**6 *+0.24159882636341E+05*xs1**0*xs2**2*xst**4 *+0.12125087320692E+06*xs1**0*xs2**4*xst**2 *-0.63007628874696E+07*xs1**0*xs2**6*xst**0 *-0.30411762363758E+03*xs1**1*xs2**0*xst**5 *-0.78383097659001E+05*xs1**1*xs2**2*xst**3 *+0.58781557929045E+05*xs1**1*xs2**4*xst**1 *+0.13375434109158E+03*xs1**2*xs2**0*xst**4 *+0.33893910748383E+06*xs1**2*xs2**2*xst**2 *-0.37826889997246E+08*xs1**2*xs2**4*xst**0 *+0.62361486440787E+04*xs1**3*xs2**0*xst**3 *-0.10596388642131E+06*xs1**3*xs2**2*xst**1 *+0.28702790231058E+05*xs1**4*xs2**0*xst**2 *-0.37680124674448E+08*xs1**4*xs2**2*xst**0 *+0.76453502516977E+05*xs1**5*xs2**0*xst**1 *-0.65507833178799E+07*xs1**6*xs2**0*xst**0 *-0.19940650191011E+04*xs1**0*xs2**0*xst**7 *-0.53454319818075E+05*xs1**0*xs2**2*xst**5 *-0.18482082317071E+06*xs1**0*xs2**4*xst**3 *-0.32482357937925E+06*xs1**0*xs2**6*xst**1 *+0.14750414212517E+05*xs1**1*xs2**0*xst**6 *+0.10992166354556E+06*xs1**1*xs2**2*xst**4 *-0.42127450223572E+06*xs1**1*xs2**4*xst**2 *+0.31697813513995E+08*xs1**1*xs2**6*xst**0 *-0.76319535999048E+05*xs1**2*xs2**0*xst**5 *-0.39996355835188E+04*xs1**2*xs2**2*xst**3 *-0.10848426889348E+06*xs1**2*xs2**4*xst**1 *+0.16031984625558E+06*xs1**3*xs2**0*xst**4 *-0.99440969626463E+06*xs1**3*xs2**2*xst**2 *+0.68374298978030E+08*xs1**3*xs2**4*xst**0 *-0.71454371570541E+05*xs1**4*xs2**0*xst**3 *+0.18244007909156E+07*xs1**4*xs2**2*xst**1 *+0.12228399205468E+06*xs1**5*xs2**0*xst**2 *+0.50353404665938E+08*xs1**5*xs2**2*xst**0 *-0.20016358839030E+06*xs1**6*xs2**0*xst**1 *+0.90962316465138E+07*xs1**7*xs2**0*xst**0 *+0.10149096018617E+04*xs1**0*xs2**0*xst**8 *+0.32566540304056E+05*xs1**0*xs2**2*xst**6 *+0.17026834226150E+05*xs1**0*xs2**4*xst**4 *-0.20338737696627E+06*xs1**0*xs2**6*xst**2 *-0.11419910468431E+08*xs1**0*xs2**8*xst**0 *-0.20146404554912E+05*xs1**1*xs2**0*xst**7 *+0.20490886901277E+05*xs1**1*xs2**2*xst**5 *+0.11627387113119E+07*xs1**1*xs2**4*xst**3 *-0.26762513146301E+07*xs1**1*xs2**6*xst**1 *+0.63330899794849E+05*xs1**2*xs2**0*xst**6 *-0.29823655377502E+06*xs1**2*xs2**2*xst**4 *-0.11572269239294E+04*xs1**2*xs2**4*xst**2 *-0.53010926363911E+08*xs1**2*xs2**6*xst**0 *-0.58084107449611E+05*xs1**3*xs2**0*xst**5 *+0.10772705558304E+07*xs1**3*xs2**2*xst**3 *-0.47494756934226E+05*xs1**3*xs2**4*xst**1 *+0.21811000394523E+05*xs1**4*xs2**0*xst**4 *+0.14998583534579E+07*xs1**4*xs2**2*xst**2 *-0.64501462536687E+08*xs1**4*xs2**4*xst**0 *-0.14570114276761E+06*xs1**5*xs2**0*xst**3 *-0.66198424020631E+05*xs1**5*xs2**2*xst**1 *+0.32017806514222E+06*xs1**6*xs2**0*xst**2 vp2=-0.46302044457094E+08*xs1**6*xs2**2*xst**0 *+0.13193687870007E+06*xs1**7*xs2**0*xst**1 *-0.11157788556702E+08*xs1**8*xs2**0*xst**0 *+0.20805390750196E+04*xs1**0*xs2**0*xst**9 *+0.10977717065665E+06*xs1**0*xs2**2*xst**7 *+0.40544491766050E+06*xs1**0*xs2**4*xst**5 *-0.24296852330663E+05*xs1**0*xs2**6*xst**3 *+0.38177105134606E+07*xs1**0*xs2**8*xst**1 *-0.16390212024836E+05*xs1**1*xs2**0*xst**8 *-0.35305243614455E+06*xs1**1*xs2**2*xst**6 *+0.70273286813606E+06*xs1**1*xs2**4*xst**4 *+0.16135959231803E+05*xs1**1*xs2**6*xst**2 *+0.76303520865603E+04*xs1**1*xs2**8*xst**0 *+0.14575708768624E+06*xs1**2*xs2**0*xst**7 *+0.85864519904820E+06*xs1**2*xs2**2*xst**5 *-0.53581284988223E+07*xs1**2*xs2**4*xst**3 *+0.12582250565590E+08*xs1**2*xs2**6*xst**1 *-0.37156076835264E+06*xs1**3*xs2**0*xst**6 *-0.47107091715944E+06*xs1**3*xs2**2*xst**4 *+0.10422174650334E+08*xs1**3*xs2**4*xst**2 *-0.73996618932697E+08*xs1**3*xs2**6*xst**0 *+0.25959907471094E+06*xs1**4*xs2**0*xst**5 *-0.34057704907398E+07*xs1**4*xs2**2*xst**3 *-0.17848038023478E+05*xs1**4*xs2**4*xst**1 *-0.29021508562549E+06*xs1**5*xs2**0*xst**4 *+0.10456945081877E+08*xs1**5*xs2**2*xst**2 *-0.71012605744803E+08*xs1**5*xs2**4*xst**0 *-0.64745699046240E+06*xs1**6*xs2**0*xst**3 *-0.34346142250957E+08*xs1**6*xs2**2*xst**1 *+0.12578165297978E+06*xs1**7*xs2**0*xst**2 *+0.18473904269184E+08*xs1**7*xs2**2*xst**0 *-0.35458625565554E+06*xs1**8*xs2**0*xst**1 *+0.97231092585162E+07*xs1**9*xs2**0*xst**0 *-0.97725063032968E+03*xs1**0*xs2**0*xst**10 *-0.58810552328534E+05*xs1**0*xs2**2*xst**8 *-0.15097015389582E+06*xs1**0*xs2**4*xst**6 *+0.28619491235330E+05*xs1**0*xs2**6*xst**4 *+0.64423406036478E+07*xs1**0*xs2**8*xst**2 *+0.13304147310031E+05*xs1**0*xs2**10*xst**0 *+0.24818228591588E+05*xs1**1*xs2**0*xst**9 *+0.82145427857098E+05*xs1**1*xs2**2*xst**7 *-0.28560408228532E+07*xs1**1*xs2**4*xst**5 *+0.17758056545809E+07*xs1**1*xs2**6*xst**3 *+0.47618501940831E+08*xs1**1*xs2**8*xst**1 *-0.14362715517470E+06*xs1**2*xs2**0*xst**8 *+0.87939534136632E+06*xs1**2*xs2**2*xst**6 *+0.43542816466615E+05*xs1**2*xs2**4*xst**4 *+0.40855419184907E+08*xs1**2*xs2**6*xst**2 *+0.15590054335119E+09*xs1**2*xs2**8*xst**0 *+0.23847828237922E+06*xs1**3*xs2**0*xst**7 *-0.30988931738056E+07*xs1**3*xs2**2*xst**5 *+0.81912153955218E+07*xs1**3*xs2**4*xst**3 *-0.14524552518400E+09*xs1**3*xs2**6*xst**1 *-0.27542471524966E+06*xs1**4*xs2**0*xst**6 *+0.27085838951528E+07*xs1**4*xs2**2*xst**4 *-0.48766162840619E+08*xs1**4*xs2**4*xst**2 *+0.64814050166329E+09*xs1**4*xs2**6*xst**0 *+0.10272921778939E+07*xs1**5*xs2**0*xst**5 *-0.61301714047235E+07*xs1**5*xs2**2*xst**3 *+0.48355468765097E+07*xs1**5*xs2**4*xst**1 *-0.16950012608773E+07*xs1**6*xs2**0*xst**4 *-0.36875319697041E+08*xs1**6*xs2**2*xst**2 *+0.32375440102365E+09*xs1**6*xs2**4*xst**0 *+0.20890611536445E+07*xs1**7*xs2**0*xst**3 *+0.85634774262523E+08*xs1**7*xs2**2*xst**1 *-0.28453560365440E+07*xs1**8*xs2**0*xst**2 *-0.70110715851666E+05*xs1**8*xs2**2*xst**0 *+0.16593108812323E+06*xs1**9*xs2**0*xst**1 *-0.46099243127663E+07*xs1**10*xs2**0*xst**0 *-0.12157030586644E+04*xs1**0*xs2**0*xst**11 *-0.82857116167377E+05*xs1**0*xs2**2*xst**9 *-0.58451292296010E+06*xs1**0*xs2**4*xst**7 *+0.11599886074607E+07*xs1**0*xs2**6*xst**5 *-0.81818449411873E+07*xs1**0*xs2**8*xst**3 *-0.38544362571641E+08*xs1**0*xs2**10*xst**1 *+0.13772894688675E+05*xs1**1*xs2**0*xst**10 *+0.45365215594235E+06*xs1**1*xs2**2*xst**8 *+0.79988903784077E+06*xs1**1*xs2**4*xst**6 *+0.35926515007257E+05*xs1**1*xs2**6*xst**4 *-0.92442595978415E+08*xs1**1*xs2**8*xst**2 *+0.43271322021484E+08*xs1**1*xs2**10*xst**0 *-0.17400692445498E+06*xs1**2*xs2**0*xst**9 vp3=-0.14268059220845E+07*xs1**2*xs2**2*xst**7 *+0.40545378149964E+07*xs1**2*xs2**4*xst**5 *+0.17585280543815E+05*xs1**2*xs2**6*xst**3 *-0.10276063972665E+09*xs1**2*xs2**8*xst**1 *+0.56807964581015E+06*xs1**3*xs2**0*xst**8 *+0.18437826095465E+07*xs1**3*xs2**2*xst**6 *-0.12864944253983E+08*xs1**3*xs2**4*xst**4 *-0.10451574046722E+09*xs1**3*xs2**6*xst**2 *-0.43651678669326E+09*xs1**3*xs2**8*xst**0 *-0.73083658333429E+06*xs1**4*xs2**0*xst**7 *-0.17056814638312E+05*xs1**4*xs2**2*xst**5 *+0.59610184323397E+08*xs1**4*xs2**4*xst**3 *+0.49938980132233E+09*xs1**4*xs2**6*xst**1 *+0.68167414002953E+06*xs1**5*xs2**0*xst**6 *-0.59994275353066E+07*xs1**5*xs2**2*xst**4 *+0.58536371084931E+04*xs1**5*xs2**4*xst**2 *-0.14136428052005E+10*xs1**5*xs2**6*xst**0 *+0.95044891506254E+05*xs1**6*xs2**0*xst**5 *+0.50803779370631E+08*xs1**6*xs2**2*xst**3 *-0.17596700209434E+05*xs1**6*xs2**4*xst**1 *+0.24656515735182E+07*xs1**7*xs2**0*xst**4 *+0.27067874384387E+08*xs1**7*xs2**2*xst**2 *-0.30519783492194E+09*xs1**7*xs2**4*xst**0 *+0.99859651431732E+05*xs1**8*xs2**0*xst**3 *-0.58311937002508E+08*xs1**8*xs2**2*xst**1 *+0.31058183048947E+07*xs1**9*xs2**0*xst**2 *+0.33566859753505E+05*xs1**9*xs2**2*xst**0 *+0.28534012133708E+03*xs1**0*xs2**0*xst**12 *+0.73181259963006E+05*xs1**0*xs2**2*xst**10 *+0.30712308016695E+06*xs1**0*xs2**4*xst**8 *-0.38988311746447E+06*xs1**0*xs2**6*xst**6 *+0.65071010834709E+04*xs1**0*xs2**8*xst**4 *+0.33025369570808E+08*xs1**0*xs2**10*xst**2 *-0.17164906136195E+08*xs1**0*xs2**12*xst**0 *-0.13913272952356E+05*xs1**1*xs2**0*xst**11 *-0.33399040011605E+06*xs1**1*xs2**2*xst**9 *+0.10266852525562E+07*xs1**1*xs2**4*xst**7 *-0.39836520885698E+07*xs1**1*xs2**6*xst**5 *+0.34049616900886E+08*xs1**1*xs2**8*xst**3 *-0.21561729380190E+04*xs1**1*xs2**10*xst**1 *+0.14552323333248E+06*xs1**2*xs2**0*xst**10 *+0.59410384696083E+06*xs1**2*xs2**2*xst**8 *-0.23886244165805E+07*xs1**2*xs2**4*xst**6 *+0.50371572986824E+07*xs1**2*xs2**6*xst**4 *+0.96234635335528E+08*xs1**2*xs2**8*xst**2 *-0.88806690050628E+08*xs1**2*xs2**10*xst**0 *-0.43369177090809E+06*xs1**3*xs2**0*xst**9 *+0.33593315021506E+05*xs1**3*xs2**2*xst**7 *+0.88172686594762E+05*xs1**3*xs2**4*xst**5 *-0.44813807185195E+08*xs1**3*xs2**6*xst**3 *+0.11963770299610E+09*xs1**3*xs2**8*xst**1 *+0.59120204513323E+06*xs1**4*xs2**0*xst**8 *+0.12298199958256E+06*xs1**4*xs2**2*xst**6 *+0.18693214990721E+08*xs1**4*xs2**4*xst**4 *+0.13670141790264E+09*xs1**4*xs2**6*xst**2 *+0.41735524210803E+09*xs1**4*xs2**8*xst**0 *-0.11245306922000E+07*xs1**5*xs2**0*xst**7 *+0.17460051518372E+07*xs1**5*xs2**2*xst**5 *-0.11089446429012E+09*xs1**5*xs2**4*xst**3 *-0.58155496949524E+09*xs1**5*xs2**6*xst**1 *+0.14312315939245E+07*xs1**6*xs2**0*xst**6 *+0.26916166982085E+06*xs1**6*xs2**2*xst**4 *+0.10382617150293E+09*xs1**6*xs2**4*xst**2 *+0.10618975583625E+10*xs1**6*xs2**6*xst**0 *-0.34896597462504E+07*xs1**7*xs2**0*xst**5 *-0.50144885137780E+08*xs1**7*xs2**2*xst**3 vp=vp1+vp2+vp3 vps1=42395.535333d0*xep1 vps2=42395.535333d0*xep2 y1=1.d0/(1.d0+dexp(ut*(x0-r1))) y2=1.d0/(1.d0+dexp(ut*(x0-r2))) y12=1.d0/(1.d0+dexp(ut2*(x02-r1))) y22=1.d0/(1.d0+dexp(ut2*(x02-r2))) vp=vp*xep3*(1-y12)*(1-y22) voh1=vpb1*(1-y1)+y1*vps1 voh2=vpb2*(1-y2)+y2*vps2 v=v0+vp+voh1+voh2+vhh return end ! wifin from ai 10.2Dr subroutine poten(xp,v,r1,r2,th) integer,parameter :: npropin = 246 double precision,intent(in) :: r1,r2,th double precision,intent(out):: v double precision,intent(in) :: xp(1:npropin) double precision :: reoh,thetae,b1,roh,alphaoh double precision :: phh2,t0,ut,x0,ut2,x02,xs1,xs2,xst,rs,rm,rr1,rr2,xep1,xep2,xep3 double precision :: rhh,vhh,vpb1,vpb2,v0,vp1,vp2,vp3,vp,vps1,vps2,y1,y2,y12,y22,voh1,voh2 integer(4) :: i ! reoh=0.95860d0 ! thetae=104.3475d0 thetae=104.4800d0 b1=2.15d0 ! roh=0.9578d0 alphaoh=0.0d0 phh2=6.20164303995d0 thetae=thetae*.314159265358979312d01*.00555555555555555555d0 xs1=(r1+r2)*0.5d0-reoh xs2=(r1-r2)*0.5d0 xst=dcos(th)-dcos(thetae) rr1=r1-reoh rr2=r2-reoh alpha=2.2668d0 ! alpha=1.5d0 xep1=(dexp(-2.d0*alpha*rr1)-2.d0*dexp(-alpha*rr1))* & (0.4389830771267D+05 & -0.8606164763375D+03*rr1 & +0.1746549380317D+05*rr1**2 & -0.3326053761996D+05*rr1**3 & +0.4988743744461D+04*rr1**4 & +0.8751979626215D+04*rr1**5 & -0.3770596827415D+04*rr1**6 & +0.4270854122634D+03*rr1**7) + 0.4389830771267D+05 xep2=(dexp(-2.d0*alpha*rr2)-2.d0*dexp(-alpha*rr2))* & (0.4389830771267D+05 & -0.8606164763375D+03*rr2 & +0.1746549380317D+05*rr2**2 & -0.3326053761996D+05*rr2**3 & +0.4988743744461D+04*rr2**4 & +0.8751979626215D+04*rr2**5 & -0.3770596827415D+04*rr2**6 & +0.4270854122634D+03*rr2**7) + 0.4389830771267D+05 xep3=dexp(-b1*((r1-reoh)**2+(r2-reoh)**2)) rhh=dsqrt(r1**2+r2**2-2.d0*r1*r2*dcos(th)) vhh=0.820894739261131734D+06*dexp(-phh2*rhh) v0 = xp( 1) *xs1**0 *xs2**0 *xst**0 vp1=+xp( 2) *xs1**0 *xs2**0 *xst**1 & +xp( 3) *xs1**1 *xs2**0 *xst**0 & +xp( 4) *xs1**0 *xs2**0 *xst**2 & +xp( 5) *xs1**0 *xs2**2 *xst**0 & +xp( 6) *xs1**1 *xs2**0 *xst**1 & +xp( 7) *xs1**2 *xs2**0 *xst**0 & +xp( 8) *xs1**0 *xs2**0 *xst**3 & +xp( 9) *xs1**0 *xs2**2 *xst**1 & +xp( 10) *xs1**1 *xs2**0 *xst**2 & +xp( 11) *xs1**1 *xs2**2 *xst**0 & +xp( 12) *xs1**2 *xs2**0 *xst**1 & +xp( 13) *xs1**3 *xs2**0 *xst**0 & +xp( 14) *xs1**0 *xs2**0 *xst**4 & +xp( 15) *xs1**0 *xs2**2 *xst**2 & +xp( 16) *xs1**0 *xs2**4 *xst**0 & +xp( 17) *xs1**1 *xs2**0 *xst**3 & +xp( 18) *xs1**1 *xs2**2 *xst**1 & +xp( 19) *xs1**2 *xs2**0 *xst**2 & +xp( 20) *xs1**2 *xs2**2 *xst**0 & +xp( 21) *xs1**3 *xs2**0 *xst**1 & +xp( 22) *xs1**4 *xs2**0 *xst**0 & +xp( 23) *xs1**0 *xs2**0 *xst**5 & +xp( 24) *xs1**0 *xs2**2 *xst**3 & +xp( 25) *xs1**0 *xs2**4 *xst**1 & +xp( 26) *xs1**1 *xs2**0 *xst**4 & +xp( 27) *xs1**1 *xs2**2 *xst**2 & +xp( 28) *xs1**1 *xs2**4 *xst**0 & +xp( 29) *xs1**2 *xs2**0 *xst**3 & +xp( 30) *xs1**2 *xs2**2 *xst**1 & +xp( 31) *xs1**3 *xs2**0 *xst**2 & +xp( 32) *xs1**3 *xs2**2 *xst**0 & +xp( 33) *xs1**4 *xs2**0 *xst**1 & +xp( 34) *xs1**5 *xs2**0 *xst**0 & +xp( 35) *xs1**0 *xs2**0 *xst**6 & +xp( 36) *xs1**0 *xs2**2 *xst**4 & +xp( 37) *xs1**0 *xs2**4 *xst**2 & +xp( 38) *xs1**0 *xs2**6 *xst**0 & +xp( 39) *xs1**1 *xs2**0 *xst**5 & +xp( 40) *xs1**1 *xs2**2 *xst**3 & +xp( 41) *xs1**1 *xs2**4 *xst**1 & +xp( 42) *xs1**2 *xs2**0 *xst**4 & +xp( 43) *xs1**2 *xs2**2 *xst**2 & +xp( 44) *xs1**2 *xs2**4 *xst**0 & +xp( 45) *xs1**3 *xs2**0 *xst**3 & +xp( 46) *xs1**3 *xs2**2 *xst**1 & +xp( 47) *xs1**4 *xs2**0 *xst**2 & +xp( 48) *xs1**4 *xs2**2 *xst**0 & +xp( 49) *xs1**5 *xs2**0 *xst**1 & +xp( 50) *xs1**6 *xs2**0 *xst**0 & +xp( 51) *xs1**0 *xs2**0 *xst**7 & +xp( 52) *xs1**0 *xs2**2 *xst**5 & +xp( 53) *xs1**0 *xs2**4 *xst**3 & +xp( 54) *xs1**0 *xs2**6 *xst**1 & +xp( 55) *xs1**1 *xs2**0 *xst**6 & +xp( 56) *xs1**1 *xs2**2 *xst**4 & +xp( 57) *xs1**1 *xs2**4 *xst**2 & +xp( 58) *xs1**1 *xs2**6 *xst**0 & +xp( 59) *xs1**2 *xs2**0 *xst**5 & +xp( 60) *xs1**2 *xs2**2 *xst**3 & +xp( 61) *xs1**2 *xs2**4 *xst**1 & +xp( 62) *xs1**3 *xs2**0 *xst**4 & +xp( 63) *xs1**3 *xs2**2 *xst**2 & +xp( 64) *xs1**3 *xs2**4 *xst**0 & +xp( 65) *xs1**4 *xs2**0 *xst**3 & +xp( 66) *xs1**4 *xs2**2 *xst**1 & +xp( 67) *xs1**5 *xs2**0 *xst**2 & +xp( 68) *xs1**5 *xs2**2 *xst**0 & +xp( 69) *xs1**6 *xs2**0 *xst**1 & +xp( 70) *xs1**7 *xs2**0 *xst**0 & +xp( 71) *xs1**0 *xs2**0 *xst**8 & +xp( 72) *xs1**0 *xs2**2 *xst**6 & +xp( 73) *xs1**0 *xs2**4 *xst**4 & +xp( 74) *xs1**0 *xs2**6 *xst**2 & +xp( 75) *xs1**0 *xs2**8 *xst**0 & +xp( 76) *xs1**1 *xs2**0 *xst**7 & +xp( 77) *xs1**1 *xs2**2 *xst**5 & +xp( 78) *xs1**1 *xs2**4 *xst**3 & +xp( 79) *xs1**1 *xs2**6 *xst**1 & +xp( 80) *xs1**2 *xs2**0 *xst**6 & +xp( 81) *xs1**2 *xs2**2 *xst**4 & +xp( 82) *xs1**2 *xs2**4 *xst**2 & +xp( 83) *xs1**2 *xs2**6 *xst**0 & +xp( 84) *xs1**3 *xs2**0 *xst**5 & +xp( 85) *xs1**3 *xs2**2 *xst**3 & +xp( 86) *xs1**3 *xs2**4 *xst**1 & +xp( 87) *xs1**4 *xs2**0 *xst**4 & +xp( 88) *xs1**4 *xs2**2 *xst**2 & +xp( 89) *xs1**4 *xs2**4 *xst**0 & +xp( 90) *xs1**5 *xs2**0 *xst**3 & +xp( 91) *xs1**5 *xs2**2 *xst**1 & +xp( 92) *xs1**6 *xs2**0 *xst**2 vp2=+xp( 93) *xs1**6 *xs2**2 *xst**0 & +xp( 94) *xs1**7 *xs2**0 *xst**1 & +xp( 95) *xs1**8 *xs2**0 *xst**0 & +xp( 96) *xs1**0 *xs2**0 *xst**9 & +xp( 97) *xs1**0 *xs2**2 *xst**7 & +xp( 98) *xs1**0 *xs2**4 *xst**5 & +xp( 99) *xs1**0 *xs2**6 *xst**3 & +xp(100) *xs1**0 *xs2**8 *xst**1 & +xp(101) *xs1**1 *xs2**0 *xst**8 & +xp(102) *xs1**1 *xs2**2 *xst**6 & +xp(103) *xs1**1 *xs2**4 *xst**4 & +xp(104) *xs1**1 *xs2**6 *xst**2 & +xp(105) *xs1**1 *xs2**8 *xst**0 & +xp(106) *xs1**2 *xs2**0 *xst**7 & +xp(107) *xs1**2 *xs2**2 *xst**5 & +xp(108) *xs1**2 *xs2**4 *xst**3 & +xp(109) *xs1**2 *xs2**6 *xst**1 & +xp(110) *xs1**3 *xs2**0 *xst**6 & +xp(111) *xs1**3 *xs2**2 *xst**4 & +xp(112) *xs1**3 *xs2**4 *xst**2 & +xp(113) *xs1**3 *xs2**6 *xst**0 & +xp(114) *xs1**4 *xs2**0 *xst**5 & +xp(115) *xs1**4 *xs2**2 *xst**3 & +xp(116) *xs1**4 *xs2**4 *xst**1 & +xp(117) *xs1**5 *xs2**0 *xst**4 & +xp(118) *xs1**5 *xs2**2 *xst**2 & +xp(119) *xs1**5 *xs2**4 *xst**0 & +xp(120) *xs1**6 *xs2**0 *xst**3 & +xp(121) *xs1**6 *xs2**2 *xst**1 & +xp(122) *xs1**7 *xs2**0 *xst**2 & +xp(123) *xs1**7 *xs2**2 *xst**0 & +xp(124) *xs1**8 *xs2**0 *xst**1 & +xp(125) *xs1**9 *xs2**0 *xst**0 & +xp(126) *xs1**0 *xs2**2 *xst**8 & +xp(127) *xs1**0 *xs2**4 *xst**6 & +xp(128) *xs1**0 *xs2**6 *xst**4 & +xp(129) *xs1**0 *xs2**8 *xst**2 & +xp(130) *xs1**1 *xs2**0 *xst**9 & +xp(131) *xs1**1 *xs2**2 *xst**7 & +xp(132) *xs1**1 *xs2**4 *xst**5 & +xp(133) *xs1**1 *xs2**6 *xst**3 & +xp(134) *xs1**1 *xs2**8 *xst**1 & +xp(135) *xs1**2 *xs2**0 *xst**8 & +xp(136) *xs1**2 *xs2**2 *xst**6 & +xp(137) *xs1**2 *xs2**4 *xst**4 & +xp(138) *xs1**2 *xs2**6 *xst**2 & +xp(139) *xs1**3 *xs2**0 *xst**7 & +xp(140) *xs1**3 *xs2**2 *xst**5 & +xp(141) *xs1**3 *xs2**4 *xst**3 & +xp(142) *xs1**3 *xs2**6 *xst**1 & +xp(143) *xs1**4 *xs2**0 *xst**6 & +xp(144) *xs1**4 *xs2**2 *xst**4 & +xp(145) *xs1**4 *xs2**4 *xst**2 & +xp(146) *xs1**5 *xs2**0 *xst**5 & +xp(147) *xs1**5 *xs2**2 *xst**3 & +xp(148) *xs1**5 *xs2**4 *xst**1 & +xp(149) *xs1**6 *xs2**0 *xst**4 & +xp(150) *xs1**6 *xs2**2 *xst**2 & +xp(151) *xs1**6 *xs2**4 *xst**0 & +xp(152) *xs1**7 *xs2**0 *xst**3 & +xp(153) *xs1**7 *xs2**2 *xst**1 & +xp(154) *xs1**8 *xs2**0 *xst**2 & +xp(155) *xs1**8 *xs2**2 *xst**0 & +xp(156) *xs1**9 *xs2**0 *xst**1 & +xp(157) *xs1**10*xs2**0 *xst**0 & +xp(158) *xs1**0 *xs2**2 *xst**9 & +xp(159) *xs1**0 *xs2**4 *xst**7 & +xp(160) *xs1**0 *xs2**6 *xst**5 & +xp(161) *xs1**0 *xs2**8 *xst**3 & +xp(162) *xs1**0 *xs2**10*xst**1 & +xp(163) *xs1**1 *xs2**0 *xst**10 & +xp(164) *xs1**1 *xs2**2 *xst**8 & +xp(165) *xs1**1 *xs2**4 *xst**6 & +xp(166) *xs1**1 *xs2**6 *xst**4 & +xp(167) *xs1**1 *xs2**8 *xst**2 & +xp(168) *xs1**1 *xs2**10*xst**0 & +xp(169) *xs1**2 *xs2**0 *xst**9 & +xp(170) *xs1**2 *xs2**2 *xst**7 & +xp(171) *xs1**2 *xs2**4 *xst**5 & +xp(172) *xs1**2 *xs2**6 *xst**3 & +xp(173) *xs1**2 *xs2**8 *xst**1 & +xp(174) *xs1**3 *xs2**0 *xst**8 vp3=+xp(175) *xs1**3 *xs2**2 *xst**6 & +xp(176) *xs1**3 *xs2**4 *xst**4 & +xp(177) *xs1**3 *xs2**6 *xst**2 & +xp(178) *xs1**3 *xs2**8 *xst**0 & +xp(179) *xs1**4 *xs2**0 *xst**7 & +xp(180) *xs1**4 *xs2**2 *xst**5 & +xp(181) *xs1**4 *xs2**4 *xst**3 & +xp(182) *xs1**4 *xs2**6 *xst**1 & +xp(183) *xs1**5 *xs2**0 *xst**6 & +xp(184) *xs1**5 *xs2**2 *xst**4 & +xp(185) *xs1**5 *xs2**4 *xst**2 & +xp(186) *xs1**5 *xs2**6 *xst**0 & +xp(187) *xs1**6 *xs2**0 *xst**5 & +xp(188) *xs1**6 *xs2**2 *xst**3 & +xp(189) *xs1**6 *xs2**4 *xst**1 & +xp(190) *xs1**7 *xs2**0 *xst**4 & +xp(191) *xs1**7 *xs2**2 *xst**2 & +xp(192) *xs1**8 *xs2**0 *xst**3 & +xp(193) *xs1**8 *xs2**2 *xst**1 & +xp(194) *xs1**9 *xs2**0 *xst**2 & +xp(195) *xs1**9 *xs2**2 *xst**0 & +xp(196) *xs1**10*xs2**0 *xst**1 & +xp(197) *xs1**11*xs2**0 *xst**0 & +xp(198) *xs1**0 *xs2**0 *xst**2 & +xp(199) *xs1**0 *xs2**2 *xst**0 & +xp(200) *xs1**0 *xs2**4 *xst**8 & +xp(201) *xs1**0 *xs2**6 *xst**6 & +xp(202) *xs1**0 *xs2**8 *xst**4 & +xp(203) *xs1**0 *xs2**10*xst**2 & +xp(204) *xs1**0 *xs2**12*xst**0 & +xp(205) *xs1**1 *xs2**0 *xst**11 & +xp(206) *xs1**1 *xs2**2 *xst**9 & +xp(207) *xs1**1 *xs2**4 *xst**7 & +xp(208) *xs1**1 *xs2**6 *xst**5 & +xp(209) *xs1**1 *xs2**8 *xst**3 & +xp(210) *xs1**1 *xs2**10*xst**1 & +xp(211) *xs1**2 *xs2**0 *xst**10 & +xp(212) *xs1**2 *xs2**2 *xst**8 & +xp(213) *xs1**2 *xs2**4 *xst**6 & +xp(214) *xs1**2 *xs2**6 *xst**4 & +xp(215) *xs1**2 *xs2**8 *xst**2 & +xp(216) *xs1**2 *xs2**10*xst**0 & +xp(217) *xs1**3 *xs2**0 *xst**9 & +xp(218) *xs1**3 *xs2**2 *xst**7 & +xp(219) *xs1**3 *xs2**4 *xst**5 & +xp(220) *xs1**3 *xs2**6 *xst**3 & +xp(221) *xs1**3 *xs2**8 *xst**1 & +xp(222) *xs1**4 *xs2**0 *xst**8 & +xp(223) *xs1**4 *xs2**2 *xst**6 & +xp(224) *xs1**4 *xs2**4 *xst**4 & +xp(225) *xs1**4 *xs2**6 *xst**2 & +xp(226) *xs1**4 *xs2**8 *xst**0 & +xp(227) *xs1**5 *xs2**0 *xst**7 & +xp(228) *xs1**5 *xs2**2 *xst**5 & +xp(229) *xs1**5 *xs2**4 *xst**3 & +xp(230) *xs1**5 *xs2**6 *xst**1 & +xp(231) *xs1**6 *xs2**0 *xst**6 & +xp(232) *xs1**6 *xs2**2 *xst**4 & +xp(233) *xs1**6 *xs2**4 *xst**2 & +xp(234) *xs1**6 *xs2**6 *xst**0 & +xp(235) *xs1**7 *xs2**0 *xst**5 & +xp(236) *xs1**7 *xs2**2 *xst**3 & +xp(237) *xs1**7 *xs2**4 *xst**1 & +xp(238) *xs1**8 *xs2**0 *xst**4 & +xp(239) *xs1**8 *xs2**2 *xst**2 & +xp(240) *xs1**8 *xs2**4 *xst**0 & +xp(241) *xs1**9 *xs2**0 *xst**3 & +xp(242) *xs1**9 *xs2**2 *xst**1 & +xp(243) *xs1**10*xs2**0 *xst**2 & +xp(244) *xs1**10*xs2**2 *xst**0 & +xp(245) *xs1**11*xs2**0 *xst**1 & +xp(246) *xs1**12*xs2**0 *xst**0 vp=vp1+vp2+vp3 v=v0+vp*xep3+vhh+xep1+xep2 ! v=v0+vp*xep3+xep1+xep2 end subroutine poten