! ! Defining potential energy function ! subroutine poten(f,r12,r32,theta) ! implicit none ! ! r12,r32 are in Angstrom, theta in rad and f is the potential energy function value in cm-1 ! integer,parameter :: N = 99 ! double precision,intent(in) :: r12,r32,theta double precision,intent(out):: f double precision :: xcos,v,v1,v2,v3,v4,v5,v6,v7,v8,alpha double precision :: force(N) double precision :: aa1,re12,alphae,xst,y1,y2,y3,xs1,xs2,v0,vp1,vp2,vp3 double precision :: g1,g2,b1,b2,rhh,vhh,th1,th2 force( 1)= 8.000000000000E+05 ! b1 force( 2)= 8.000000000000E+04 ! b2 force( 3)= 1.300000000000E+01 ! g1 force( 4)= 5.500000000000E+00 ! g2 force( 5)= 0.000000000000E+00 ! f000 force( 6)= 2.529872472830E+00 ! f001 force( 7)= 7.600144603465E+00 ! f100 force( 8)= 1.911986956197E+04 ! f002 force( 9)= -2.633794236952E+03 ! f101 force( 10)= -3.464138372894E+02 ! f110 force( 11)= 3.714664033508E+04 ! f200 force( 12)= 1.037245694612E+03 ! f003 force( 13)= -4.836681196118E+03 ! f102 force( 14)= 3.117897942342E+03 ! f111 force( 15)= -1.271117618211E+03 ! f201 force( 16)= -1.872621560986E+02 ! f210 force( 17)= -1.179614898388E+03 ! f300 force( 18)= 4.747376342454E+03 ! f004 force( 19)= -1.361035046859E+03 ! f103 force( 20)= 2.914438002564E+03 ! f112 force( 21)= -5.065335793793E+03 ! f202 force( 22)= -1.109714671423E+03 ! f211 force( 23)= 1.037528298464E+02 ! f220 force( 24)= 4.804732872230E+03 ! f301 force( 25)= -1.406001811057E+02 ! f310 force( 26)= 2.242535326340E+03 ! f400 force( 27)= 1.791450002966E+03 ! f005 force( 28)= -9.818420229357E+02 ! f104 force( 29)= 6.129430506375E+01 ! f113 force( 30)= -2.715169671955E+03 ! f203 force( 31)= 2.611901366784E+03 ! f212 force( 32)= -6.170307353071E+03 ! f221 force( 33)= -4.883373618403E+03 ! f302 force( 34)= 1.004928914445E+03 ! f311 force( 35)= 2.870806133885E+02 ! f320 force( 36)= 6.785201750979E+03 ! f401 force( 37)= -1.048970863165E+03 ! f410 force( 38)= 7.472738249864E+02 ! f500 force( 39)= 2.049652341338E+03 ! f006 force( 40)= 2.480751796682E+02 ! f105 force( 41)= -1.957836278953E+03 ! f114 force( 42)= 7.016663326509E+01 ! f204 force( 43)= 5.182367675542E+03 ! f213 force( 44)= -2.973660382416E+03 ! f222 force( 45)= -3.906888774989E+03 ! f303 force( 46)= 2.270780941808E+03 ! f312 force( 47)= 2.379355196621E+04 ! f321 force( 48)= -3.493231488295E+03 ! f330 force( 49)= 8.592320739255E+02 ! f402 force( 50)= -1.132798726796E+04 ! f411 force( 51)= 3.766857130534E+02 ! f420 force( 52)= -1.595354970709E+04 ! f501 force( 53)= 1.489038647146E+03 ! f510 force( 54)= 5.148440268568E+02 ! f600 force( 55)= -5.461701474687E+03 ! f007 force( 56)= -6.994351599264E+01 ! f106 force( 57)= 1.742770552275E+03 ! f115 force( 58)= 2.914756012052E+02 ! f205 force( 59)= -3.364600432792E+03 ! f214 force( 60)= -1.120227311718E+03 ! f223 force( 61)= -3.832333563319E+02 ! f304 force( 62)= 2.644115181725E+03 ! f313 force( 63)= 1.936604957725E+03 ! f322 force( 64)= 3.817571881160E+03 ! f331 force( 65)= 3.077836781136E+03 ! f403 force( 66)= -2.635388174753E+03 ! f412 force( 67)= 8.399950057523E+03 ! f421 force( 68)= -5.245095152812E+03 ! f430 force( 69)= -3.472941663041E+02 ! f502 force( 70)= -5.005068446256E+03 ! f511 force( 71)= -2.821190432034E+03 ! f520 force( 72)= -1.058069608744E+04 ! f601 force( 73)= 5.378299620976E+03 ! f610 force( 74)= -9.973475090825E+02 ! f700 force( 75)= -7.404683471843E+02 ! f008 force( 76)= 0.000000000000E+00 ! f107 force( 77)= 0.000000000000E+00 ! f116 force( 78)= 0.000000000000E+00 ! f206 force( 79)= 0.000000000000E+00 ! f215 force( 80)= 0.000000000000E+00 ! f224 force( 81)= 0.000000000000E+00 ! f305 force( 82)= 0.000000000000E+00 ! f314 force( 83)= 0.000000000000E+00 ! f323 force( 84)= 0.000000000000E+00 ! f332 force( 85)= 0.000000000000E+00 ! f404 force( 86)= 0.000000000000E+00 ! f413 force( 87)= 0.000000000000E+00 ! f422 force( 88)= 0.000000000000E+00 ! f431 force( 89)= 0.000000000000E+00 ! f440 force( 90)= 0.000000000000E+00 ! f503 force( 91)= 0.000000000000E+00 ! f512 force( 92)= 0.000000000000E+00 ! f521 force( 93)= 0.000000000000E+00 ! f530 force( 94)= 0.000000000000E+00 ! f602 force( 95)= 0.000000000000E+00 ! f611 force( 96)= 0.000000000000E+00 ! f620 force( 97)= 0.000000000000E+00 ! f701 force( 98)= 0.000000000000E+00 ! f710 force( 99)= 0.000000000000E+00 ! f800 ! alpha = theta ! !write(6,"('r12,r32,alpha = ',3f20.8)") r12,r32,alpha ! re12 = 1.3359007d0 ! Ang alphae = 92.265883d0*3.1415926535897932385/180.d0 ! rad aa1 = 1.704d0 ! 1/Ang ! b1 = force(1) b2 = force(2) g1 = force(3) g2 = force(4) ! rhh=sqrt(r12**2+r32**2-2.d0*r12*r32*cos(alpha)) vhh=b1*exp(-g1*rhh)+b2*exp(-g2*rhh**2) ! ! calculate potential energy function values ! y1=1.0d+00-exp(-aa1*(r12-re12)) y2=1.0d+00-exp(-aa1*(r32-re12)) ! y3=(cos(alpha)-cos(alphae)) ! *exp(-1.0*(alpha-alphae)**2) ! v4 = 0 ; v5 = 0 ; v6 = 0 ; v7 = 0 ; v8 = 0 ! v0 = force(5)*y1**0*y2**0*y3**0 v1 = force(6)*y1**0*y2**0*y3**1& + force(7)*y1**1*y2**0*y3**0& + force(7)*y1**0*y2**1*y3**0 v2 = force(8)*y1**0*y2**0*y3**2& + force(9)*y1**1*y2**0*y3**1& + force(9)*y1**0*y2**1*y3**1& + force(10)*y1**1*y2**1*y3**0& + force(11)*y1**2*y2**0*y3**0& + force(11)*y1**0*y2**2*y3**0 v3 = force(12)*y1**0*y2**0*y3**3& + force(13)*y1**1*y2**0*y3**2& + force(13)*y1**0*y2**1*y3**2& + force(14)*y1**1*y2**1*y3**1& + force(15)*y1**2*y2**0*y3**1& + force(15)*y1**0*y2**2*y3**1& + force(16)*y1**2*y2**1*y3**0& + force(16)*y1**1*y2**2*y3**0& + force(17)*y1**3*y2**0*y3**0& + force(17)*y1**0*y2**3*y3**0 if (N>18) then v4 = force(18)*y1**0*y2**0*y3**4& + force(19)*y1**1*y2**0*y3**3& + force(19)*y1**0*y2**1*y3**3& + force(20)*y1**1*y2**1*y3**2& + force(21)*y1**2*y2**0*y3**2& + force(21)*y1**0*y2**2*y3**2& + force(22)*y1**2*y2**1*y3**1& + force(22)*y1**1*y2**2*y3**1& + force(23)*y1**2*y2**2*y3**0& + force(24)*y1**3*y2**0*y3**1& + force(24)*y1**0*y2**3*y3**1& + force(25)*y1**3*y2**1*y3**0& + force(25)*y1**1*y2**3*y3**0& + force(26)*y1**4*y2**0*y3**0& + force(26)*y1**0*y2**4*y3**0 endif if (N>26) then v5 = force(27)*y1**0*y2**0*y3**5& + force(28)*y1**1*y2**0*y3**4& + force(28)*y1**0*y2**1*y3**4& + force(29)*y1**1*y2**1*y3**3& + force(30)*y1**2*y2**0*y3**3& + force(30)*y1**0*y2**2*y3**3& + force(31)*y1**2*y2**1*y3**2& + force(31)*y1**1*y2**2*y3**2& + force(32)*y1**2*y2**2*y3**1& + force(33)*y1**3*y2**0*y3**2& + force(33)*y1**0*y2**3*y3**2& + force(34)*y1**3*y2**1*y3**1& + force(34)*y1**1*y2**3*y3**1& + force(35)*y1**3*y2**2*y3**0& + force(35)*y1**2*y2**3*y3**0& + force(36)*y1**4*y2**0*y3**1& + force(36)*y1**0*y2**4*y3**1& + force(37)*y1**4*y2**1*y3**0& + force(37)*y1**1*y2**4*y3**0& + force(38)*y1**5*y2**0*y3**0& + force(38)*y1**0*y2**5*y3**0 endif if (N>38) then v6 = force(39)*y1**0*y2**0*y3**6& + force(40)*y1**1*y2**0*y3**5& + force(40)*y1**0*y2**1*y3**5& + force(41)*y1**1*y2**1*y3**4& + force(42)*y1**2*y2**0*y3**4& + force(42)*y1**0*y2**2*y3**4& + force(43)*y1**2*y2**1*y3**3& + force(43)*y1**1*y2**2*y3**3& + force(44)*y1**2*y2**2*y3**2& + force(45)*y1**3*y2**0*y3**3& + force(45)*y1**0*y2**3*y3**3& + force(46)*y1**3*y2**1*y3**2& + force(46)*y1**1*y2**3*y3**2& + force(47)*y1**3*y2**2*y3**1& + force(47)*y1**2*y2**3*y3**1& + force(48)*y1**3*y2**3*y3**0& + force(49)*y1**4*y2**0*y3**2& + force(49)*y1**0*y2**4*y3**2& + force(50)*y1**4*y2**1*y3**1& + force(50)*y1**1*y2**4*y3**1& + force(51)*y1**4*y2**2*y3**0& + force(51)*y1**2*y2**4*y3**0& + force(52)*y1**5*y2**0*y3**1& + force(52)*y1**0*y2**5*y3**1& + force(53)*y1**5*y2**1*y3**0& + force(53)*y1**1*y2**5*y3**0& + force(54)*y1**6*y2**0*y3**0& + force(54)*y1**0*y2**6*y3**0 endif if (N>54) then v7 = force(55)*y1**0*y2**0*y3**7& + force(56)*y1**1*y2**0*y3**6& + force(56)*y1**0*y2**1*y3**6& + force(57)*y1**1*y2**1*y3**5& + force(58)*y1**2*y2**0*y3**5& + force(58)*y1**0*y2**2*y3**5& + force(59)*y1**2*y2**1*y3**4& + force(59)*y1**1*y2**2*y3**4& + force(60)*y1**2*y2**2*y3**3& + force(61)*y1**3*y2**0*y3**4& + force(61)*y1**0*y2**3*y3**4& + force(62)*y1**3*y2**1*y3**3& + force(62)*y1**1*y2**3*y3**3& + force(63)*y1**3*y2**2*y3**2& + force(63)*y1**2*y2**3*y3**2& + force(64)*y1**3*y2**3*y3**1& + force(65)*y1**4*y2**0*y3**3& + force(65)*y1**0*y2**4*y3**3& + force(66)*y1**4*y2**1*y3**2& + force(66)*y1**1*y2**4*y3**2& + force(67)*y1**4*y2**2*y3**1& + force(67)*y1**2*y2**4*y3**1& + force(68)*y1**4*y2**3*y3**0& + force(68)*y1**3*y2**4*y3**0& + force(69)*y1**5*y2**0*y3**2& + force(69)*y1**0*y2**5*y3**2& + force(70)*y1**5*y2**1*y3**1& + force(70)*y1**1*y2**5*y3**1& + force(71)*y1**5*y2**2*y3**0& + force(71)*y1**2*y2**5*y3**0& + force(72)*y1**6*y2**0*y3**1& + force(72)*y1**0*y2**6*y3**1& + force(73)*y1**6*y2**1*y3**0& + force(73)*y1**1*y2**6*y3**0& + force(74)*y1**7*y2**0*y3**0& + force(74)*y1**0*y2**7*y3**0 endif if (N>74) then v8 = force(75)*y1**0*y2**0*y3**8& + force(76)*y1**1*y2**0*y3**7& + force(76)*y1**0*y2**1*y3**7& + force(77)*y1**1*y2**1*y3**6& + force(78)*y1**2*y2**0*y3**6& + force(78)*y1**0*y2**2*y3**6& + force(79)*y1**2*y2**1*y3**5& + force(79)*y1**1*y2**2*y3**5& + force(80)*y1**2*y2**2*y3**4& + force(81)*y1**3*y2**0*y3**5& + force(81)*y1**0*y2**3*y3**5& + force(82)*y1**3*y2**1*y3**4& + force(82)*y1**1*y2**3*y3**4& + force(83)*y1**3*y2**2*y3**3& + force(83)*y1**2*y2**3*y3**3& + force(84)*y1**3*y2**3*y3**2& + force(85)*y1**4*y2**0*y3**4& + force(85)*y1**0*y2**4*y3**4& + force(86)*y1**4*y2**1*y3**3& + force(86)*y1**1*y2**4*y3**3& + force(87)*y1**4*y2**2*y3**2& + force(87)*y1**2*y2**4*y3**2& + force(88)*y1**4*y2**3*y3**1& + force(88)*y1**3*y2**4*y3**1& + force(89)*y1**4*y2**4*y3**0& + force(90)*y1**5*y2**0*y3**3& + force(90)*y1**0*y2**5*y3**3& + force(91)*y1**5*y2**1*y3**2& + force(91)*y1**1*y2**5*y3**2& + force(92)*y1**5*y2**2*y3**1& + force(92)*y1**2*y2**5*y3**1& + force(93)*y1**5*y2**3*y3**0& + force(93)*y1**3*y2**5*y3**0& + force(94)*y1**6*y2**0*y3**2& + force(94)*y1**0*y2**6*y3**2& + force(95)*y1**6*y2**1*y3**1& + force(95)*y1**1*y2**6*y3**1& + force(96)*y1**6*y2**2*y3**0& + force(96)*y1**2*y2**6*y3**0& + force(97)*y1**7*y2**0*y3**1& + force(97)*y1**0*y2**7*y3**1& + force(98)*y1**7*y2**1*y3**0& + force(98)*y1**1*y2**7*y3**0& + force(99)*y1**8*y2**0*y3**0& + force(99)*y1**0*y2**8*y3**0 endif f=v0+v1+v2+v3+v4+v5+v6+v7+v8+vhh ! if (v0+v1+v2>50000.0) f = v0+v1+v2+vhh end subroutine poten