program testvheo2 implicit none c real*8 R, z, theta, r1, v, vrccsdt, pi, abserr, + relerr1, relerr2, tmp, vOO integer n, i c pi=dacos(-1d0) abserr=0d0 relerr1=0d0 relerr2=0d0 n=0 write(*,*)'Using atomic units, theta in degrees' write(*,*)'Input format: nr R theta rOO vHeO2' write(*,*)'The potential is printed when "nr" equals zero or' write(*,*)'a multiple of 150. "vHeO2" may be set to zero.' 10 continue read(*,*,end=100)i,R,theta,r1,vrccsdt n=n+1 z=cos(theta*pi/180d0) call vheo2(R,z,r1,v) if (vrccsdt.lt.0d0) then tmp=dabs(v-vrccsdt) if (tmp.gt.abserr) abserr=tmp endif if (vrccsdt.gt.7.65d-4) then tmp=dabs(v/vrccsdt-1d0) if (tmp.gt.relerr1) relerr1=tmp endif if (R.gt.7d0) then tmp=dabs(v/vrccsdt-1d0) if (tmp.gt.relerr2) relerr2=tmp endif if (i.eq.150*(i/150)) then call xsigma(r1,vOO) write(*,'(1a,i3,3(1a,f10.6))')'nr=',i,' R=',R,' theta=', + theta,' rOO=',r1 write(*,200)'vHeO2 = ',vrccsdt,' (ab initio)' write(*,200)'vHeO2 = ',v,' (fit)' write(*,200)'vO2 = ',voo write(*,*) 200 format(1a,f18.12,1a) endif go to 10 100 continue write(*,*)'number of points ',n write(*,*)'Largest relative error for V>7.65e-4: ',relerr1*100d0, +' %' write(*,*)'Largest relative error for R>7: ',relerr2*100d0, +' %' write(*,*)'Largest absolute error for V<0: ',abserr end c ************************************************************************ c subroutine vheo2(R,z,r1,v) implicit none real*8 R, z, r1, v c c input: c R: distance from He to center of mass of O2 (bohr). c z=cos(theta), theta=0 corresponds to the linear geometry. c r1: O-O bond distance (bohr). c c output: c v: He + O2 interaction potential (hartree). c c This P.E.S. is based on RCCSD(T) calculations with the MOLPRO package, c P. J. Knowles, C. Hampel and H.-J. Werner, J. Chem. Phys. 99 (1993) 5219. c c The grid used includes 750 points for which v<0.1 hartree in the region c 3.05 <= R <= 15, 1.9 < r1 < 3.1 c For r1>4 and R<2.5 the potential may become unphysical c The c_n,l with n=6, 8 are fixed from a fit to points with R>10. c c This subroutine will automatically perform a self-test on the first call, c by checking the potential for one specific geometry. c c Gerrit C. Groenenboom, 21 August 1999. c logical ifirst/.true./ real*8 beta, betar, c(118) data beta/2.15d0/, betar/0.055d0/ data c/ + 4.09722719689086d+00, -5.87850086171064d+00, + 2.62860977876725d+01, 5.33977970603514d+03, + -7.03836915442398d+03, 3.08701418724401d+03, + -4.49162257820320d+02, -1.74812523863630d+04, + 2.29039582871963d+04, -9.93850262047042d+03, + 1.42335345746197d+03, 1.77736241396470d+04, + -2.33288030910082d+04, 1.01352371811363d+04, + -1.45652788815395d+03, -6.23418699018776d+03, + 8.01027388055036d+03, -3.40685290682507d+03, + 4.78878899906915d+02, -3.09385962612368d+03, + 4.11820103176701d+03, -1.82194380085292d+03, + 2.69674770177646d+02, 1.12016768336068d+04, + -1.46813538957560d+04, 6.38150949160146d+03, + -9.16271316952770d+02, -1.13367913344927d+04, + 1.48805598612998d+04, -6.47556198430401d+03, + 9.31287308474901d+02, 3.96349644539987d+03, + -5.09619063647082d+03, 2.17171045556329d+03, + -3.06482889284378d+02, 5.36526094399032d+02, + -7.10309442465364d+02, 3.14112121326586d+02, + -4.66097444499953d+01, -2.33873817730045d+03, + 3.05248585569763d+03, -1.32364829295424d+03, + 1.89906454609946d+02, 2.37889179406219d+03, + -3.11439979610650d+03, 1.35423641791217d+03, + -1.95078301958074d+02, -7.86519888471896d+02, + 1.00565170713353d+03, -4.27199326425664d+02, + 6.02789202579294d+01, -1.84135627390178d+01, + 2.42246635614066d+01, -1.11610909492054d+01, + 1.77480604349836d+00, 1.45465852414389d+02, + -1.87253678047965d+02, 8.10891496006745d+01, + -1.17030768297363d+01, -1.45733059813131d+02, + 1.89161213069541d+02, -8.24141135939573d+01, + 1.19684422084845d+01, 4.59828585364982d+01, + -5.79512931436981d+01, 2.44608084322934d+01, + -3.45233993959894d+00, -1.22239254764827d+01, + 6.99567706600854d+00, -4.94111744139816d-01, + -4.49529381223018d+00, 4.39827732944237d+00, + 6.69408361918515d-01, -3.90839419951740d+02, + -3.57488087850102d+02, 3.22042997461517d+02, + -1.39598770751104d+02, -4.18882366674752d+02, + 1.65799083233748d+02, 1.23172469401706d+01, + 1.46573401817661d+02, -6.64328930023274d+01, + -2.84334402957489d+04, 2.25544569447039d+05, + -1.93364408972475d+05, 5.12862150223910d+04, + -7.51711206190477d+04, 3.64392258273544d+05, + -3.03121711252478d+05, 8.73346907324946d+04, + -3.45519750882223d+04, 4.25719100787138d+05, + -3.75943308872289d+05, 9.45902804164792d+04, + -4.60944149675849d+03, 1.14916075693361d+05, + -1.04417118162543d+05, 2.54022051744389d+04, + 2.88402408686071d+06, -2.71785831822367d+07, + 2.52655446346514d+07, -7.03236130965375d+06, + 3.11298177116258d+06, -3.63253073945702d+07, + 3.45445576187005d+07, -9.27396845434792d+06, + 1.74021246580691d+06, -4.13399996654562d+07, + 3.66427868291820d+07, -8.60763269788007d+06, + -8.22733314352684d+04, -6.46296378998269d+06, + 5.79777234997689d+06, -1.26247393590976d+06, + -6.29973251125557d+04, 2.17132361189737d+05, + -1.87680229897387d+05, 6.28352247097358d+04/ c real*8 Ra, za, Rb, zb, pla(4), plb(4), pl(9), + Ran(4), Rbn(4), r1n(4), sum, Rn(4), f_long, + Rsave, zsave, r1save, v0 integer i, j, k, ind, kmax c if (ifirst) then write(*,'(1a)') +'=================================================', +'RCCSD(T) He+O2 interaction potential,', +'Gerrit C. Groenenboom and Izabela M. Struniewicz,', +'J. Chem. Phys. vol 113, December (2000)', +'E-mail: gerritg@theochem.kun.nl', +'=================================================' c c Save the point c Rsave=R zsave=z r1save=r1 c c ... and first check whether this point gives the right value: c R=7d0 z=0.5d0 r1=2.5d0 v0=-0.866042743602981517D-04 go to 20 endif 10 continue if (ifirst) then if (dabs((v/v0)-1d0).gt.1d-12) then write(*,*)'** ERROR in vHeO2: integrity check not passed' write(*,*)'** v(',R,',',z,',',r1,')=',v write(*,*)'** this should be ',v0 stop 230 endif c OK, check passed, do the required point R=Rsave z=zsave r1=r1save ifirst=.false. endif c 20 continue v=0d0 c if (dabs(z).gt.1d0) then write(*,*)'** ERROR in VHEO2: illegal cos(theta) = z = ',z stop 231 endif call jac2yu(R,z,r1,Ra,za,Rb,zb) call p_leg(pla,3,za) call p_leg(plb,3,zb) Ran(1)=dexp(-beta*Ra) Rbn(1)=dexp(-beta*Rb) r1n(1)=dexp(-betar*r1**3) do i=2,4 Ran(i)=Ra*Ran(i-1) Rbn(i)=Rb*Rbn(i-1) r1n(i)=r1*r1n(i-1) end do c c the r1 independent part of the short range c ind=0 do i=1,3 ind=ind+1 v=v+(Ran(1)*pla(i)+Rbn(1)*plb(i))*c(ind) end do c c the r1 dependent part of the short range c do i=1,4 do j=1,4 sum=0d0 do k=1,4 ind=ind+1 sum=sum+c(ind)*r1n(k) end do v=v+(Ran(i)*pla(j)+Rbn(i)*plb(j))*sum end do end do c c the long range c do i=1,4 Rn(i)=f_long(R,beta,4+2*i) end do call p_leg(pl,8,z) c kmax=2 do i=1,4 if (i.eq.3) kmax=kmax+1 do j=1,i+1 ind=ind+1 sum=c(ind) do k=1,kmax ind=ind+1 sum=sum+r1n(k)*c(ind) end do v=v+Rn(i)*pl(2*j-1)*sum end do end do if (ifirst) go to 10 end c ************************************************************************ c subroutine jac2yu(R,z,r1,Ra,za,Rb,zb) implicit none real*8 R, r1, z, Ra, za, Rb, zb real*8 alpha, tmpa, tmpb, x, y2 if (R.eq.0d0) then Ra=0.5d0*r1 za=-1d0 Rb=Ra zb=za else alpha=r1/(2d0*R) tmpa=dsqrt(1d0-2d0*alpha*z+alpha*alpha) tmpb=dsqrt(1d0+2d0*alpha*z+alpha*alpha) if (tmpa.eq.0d0) then za=1d0 else za=(z-alpha)/tmpa end if if (tmpb.eq.0d0) then zb=1d0 else zb=-(z+alpha)/tmpb end if x=R*z y2=R*R*(1d0-z*z) Ra=dsqrt( ( x-0.5d0*r1)**2 + y2 ) Rb=dsqrt( ( x+0.5d0*r1)**2 + y2 ) endif end c ************************************************************************ c real*8 function f_long(r,beta,n) implicit none real*8 r, beta integer n c c f_long(r,beta,n)=d_n(beta*r,n)/r^n c c Tang-Toennies damping function: c d_n(x,n) = 1 - exp(-x) * sum_k=0:n x^k / k! c c For small x (when d_n<1e-8) use: c d_n(x,n) = exp(-x)* sum_k=n+1:infty x^k/k! c real*8 sum, term, dex, rn, dn, x integer i c if (r.eq.0d0) then f_long=0d0 return endif c x=beta*r sum=1d0 term=1d0 rn=1d0 do i=1,n rn=rn*r term=term*x/dfloat(i) sum=sum+term end do dex=dexp(-x) dn=1d0-dex*sum if (dabs(dn).lt.1d-8) then c c alternative formula for small x c sum=0d0 do i=n+1,1000 term=term*x/dfloat(i) sum=sum+term if (term/sum .lt. 1d-8) go to 10 enddo write(*,*)'** ERROR in F_LONG: no convergence, x=',x,' n=',n stop 232 10 continue dn=sum*dex endif f_long=dn/rn end c ************************************************************************ c subroutine p_leg(p,lmax,z) implicit none real*8 p(0:lmax), z, di integer lmax, i c p(0) = 1.d0 if (lmax.eq.0) return p(1) = z di = 0d0 do i=1,lmax-1 di = di + 1d0 p(i+1) = ( (2.d0*di + 1.d0) * z * p(i) + - di * p(i-1) ) / (di + 1d0) end do end c ************************************************************************ c subroutine xsigma(r,v) parameter ( npts = 24000 ) implicit double precision(a-h,o-z) c r & v in atomic units. c c ground state x3sgm potential energy curve of 02 c rkr values from cosby and helms for 1.79480 < r < 4.0840 c dissociation energy 41260.2 cm-1 from brix and herzberg c zero point energy 788.1 cm-1 from cosby and helms c short range from ron friedman c long range -17.57/r**6 assumed valid for r > 5.0 c nine points 4.2< r < 5.0 added to join smoothly to long range c generated by arthur allison july 1991 using difference between c 4.5 and 5.0 from steve guberman c dimension ra(200),pr(200),c(4,200) data nc,yp0,ypm/148,-1.018436615d0,1.3494d-03/ data (ra(i),i=1,76)/ * 0.17948246d+01, 0.17963175d+01, 0.17978482d+01, 0.17993977d+01, * 0.18010229d+01, 0.18027425d+01, 0.18045189d+01, 0.18064086d+01, * 0.18083550d+01, 0.18103959d+01, 0.18125313d+01, 0.18147234d+01, * 0.18169911d+01, 0.18193154d+01, 0.18217154d+01, 0.18241720d+01, * 0.18267043d+01, 0.18292743d+01, 0.18319199d+01, 0.18346411d+01, * 0.18374001d+01, 0.18402347d+01, 0.18431449d+01, 0.18461118d+01, * 0.18491731d+01, 0.18522912d+01, 0.18554848d+01, 0.18587729d+01, * 0.18621556d+01, 0.18656138d+01, 0.18691664d+01, 0.18728136d+01, * 0.18765742d+01, 0.18804292d+01, 0.18843976d+01, 0.18884794d+01, * 0.18926746d+01, 0.18969832d+01, 0.19014241d+01, 0.19060161d+01, * 0.19107404d+01, 0.19156159d+01, 0.19206426d+01, 0.19258393d+01, * 0.19312251d+01, 0.19367998d+01, 0.19425823d+01, 0.19485916d+01, * 0.19548655d+01, 0.19613851d+01, 0.19682070d+01, 0.19753502d+01, * 0.19828524d+01, 0.19907515d+01, 0.19990851d+01, 0.20079102d+01, * 0.20173021d+01, 0.20273177d+01, 0.20380891d+01, 0.20497109d+01, * 0.20623721d+01, 0.20763183d+01, 0.20918518d+01, 0.21095208d+01, * 0.21301944d+01, 0.21555923d+01, 0.21902310d+01, 0.22157990d+01, * 0.22816991d+01, 0.23539947d+01, 0.23859878d+01, 0.24333632d+01, * 0.24715168d+01, 0.25049461d+01, 0.25353896d+01, 0.25637355d+01/ data (ra(i),i=77,148)/ * 0.25905507d+01, 0.26161565d+01, 0.26408174d+01, 0.26647036d+01, * 0.26879472d+01, 0.27106617d+01, 0.27329227d+01, 0.27547868d+01, * 0.27763297d+01, 0.27975891d+01, 0.28186029d+01, 0.28393899d+01, * 0.28599879d+01, 0.28804348d+01, 0.29007493d+01, 0.29209316d+01, * 0.29410005d+01, 0.29609938d+01, 0.29809115d+01, 0.30007726d+01, * 0.30205769d+01, 0.30403623d+01, 0.30601100d+01, 0.30798765d+01, * 0.30996431d+01, 0.31194285d+01, 0.31392328d+01, 0.31590939d+01, * 0.31790116d+01, 0.31990049d+01, 0.32190927d+01, 0.32392750d+01, * 0.32595517d+01, 0.32799797d+01, 0.33005399d+01, 0.33212513d+01, * 0.33421328d+01, 0.33632032d+01, 0.33845005d+01, 0.34060055d+01, * 0.34277752d+01, 0.34497905d+01, 0.34721271d+01, 0.34947660d+01, * 0.35177640d+01, 0.35411588d+01, 0.35649882d+01, 0.35892901d+01, * 0.36141211d+01, 0.36395569d+01, 0.36656540d+01, 0.36924881d+01, * 0.37201915d+01, 0.37488587d+01, 0.37786218d+01, 0.38096323d+01, * 0.38420978d+01, 0.38762451d+01, 0.39123200d+01, 0.39506815d+01, * 0.39917074d+01, 0.40359270d+01, 0.40839828d+01, 0.42d+01, * 0.43d+01, 0.44d+01, 0.45d+01, 0.46d+01, * 0.47d+01, 0.48d+01, 0.49d+01, 0.50d+01/ data (pr(i),i=1,76)/ * -0.12174598d-01,-0.13688451d-01,-0.15266193d-01,-0.16903654d-01, * -0.18597132d-01,-0.20343372d-01,-0.22139485d-01,-0.23982944d-01, * -0.25871496d-01,-0.27803170d-01,-0.29776206d-01,-0.31789060d-01, * -0.33840347d-01,-0.35928850d-01,-0.38053471d-01,-0.40213235d-01, * -0.42407253d-01,-0.44634729d-01,-0.46894947d-01,-0.49187246d-01, * -0.51511020d-01,-0.53865722d-01,-0.56250838d-01,-0.58665903d-01, * -0.61110484d-01,-0.63584175d-01,-0.66086604d-01,-0.68617436d-01, * -0.71176354d-01,-0.73763069d-01,-0.76377315d-01,-0.79018852d-01, * -0.81687464d-01,-0.84382962d-01,-0.87105179d-01,-0.89853955d-01, * -0.92629171d-01,-0.95430714d-01,-0.98258492d-01,-0.10111243d+00, * -0.10399248d+00,-0.10689857d+00,-0.10983069d+00,-0.11278879d+00, * -0.11577288d+00,-0.11878293d+00,-0.12181893d+00,-0.12488087d+00, * -0.12796875d+00,-0.13108257d+00,-0.13422229d+00,-0.13738792d+00, * -0.14057943d+00,-0.14379683d+00,-0.14704006d+00,-0.15030915d+00, * -0.15360407d+00,-0.15692484d+00,-0.16027149d+00,-0.16364407d+00, * -0.16704269d+00,-0.17046749d+00,-0.17391871d+00,-0.17739663d+00, * -0.18090168d+00,-0.18443440d+00,-0.18799547d+00,-0.18978691d+00, * -0.19158624d+00,-0.18978691d+00,-0.18799547d+00,-0.18443440d+00, * -0.18090168d+00,-0.17739663d+00,-0.17391871d+00,-0.17046749d+00/ data (pr(i),i=77,148)/ * -0.16704269d+00,-0.16364407d+00,-0.16027149d+00,-0.15692484d+00, * -0.15360407d+00,-0.15030915d+00,-0.14704006d+00,-0.14379683d+00, * -0.14057943d+00,-0.13738792d+00,-0.13422229d+00,-0.13108257d+00, * -0.12796875d+00,-0.12488087d+00,-0.12181893d+00,-0.11878293d+00, * -0.11577288d+00,-0.11278879d+00,-0.10983069d+00,-0.10689857d+00, * -0.10399248d+00,-0.10111243d+00,-0.98258492d-01,-0.95430714d-01, * -0.92629171d-01,-0.89853955d-01,-0.87105179d-01,-0.84382962d-01, * -0.81687464d-01,-0.79018852d-01,-0.76377315d-01,-0.73763069d-01, * -0.71176354d-01,-0.68617436d-01,-0.66086604d-01,-0.63584175d-01, * -0.61110484d-01,-0.58665903d-01,-0.56250838d-01,-0.53865722d-01, * -0.51511020d-01,-0.49187246d-01,-0.46894947d-01,-0.44634729d-01, * -0.42407253d-01,-0.40213235d-01,-0.38053471d-01,-0.35928850d-01, * -0.33840347d-01,-0.31789060d-01,-0.29776206d-01,-0.27803170d-01, * -0.25871496d-01,-0.23982944d-01,-0.22139485d-01,-0.20343372d-01, * -0.18597132d-01,-0.16903654d-01,-0.15266193d-01,-0.13688451d-01, * -0.12174598d-01,-0.10729373d-01,-0.93581104d-02,-0.690d-02, * -0.500d-02, -0.357d-02, -0.260d-02, -0.207d-02, * -0.170d-02, -0.145d-02, -0.127d-02, -0.11245d-02/ logical first/.true./ if (first) then write(*,'(1a)') +'=========================================', +'O2 potential, J. F. Babb and A. Dalgarno,', +'Phys. Rev. A, 51, 3021 (1995)', +'=========================================' first=.false. endif c c asym used only for short range fit a*exp(-alpha*r) rel v11 = 0 asym = 41260.2d0/219474.4354d0 nint = 1 call splicn(ra,pr,nc,c,yp0,ypm) c x state of o2 if ( r .gt. ra(1) ) then if ( r .lt. ra(nc) ) then c spline fit 70 if ( (r-ra(nint+1)) .lt. 0 ) goto 71 nint = nint + 1 goto 70 71 continue rl = r - ra(nint) ru = ra(nint+1) - r v = c(4,nint)*rl + c(3,nint)*ru * + c(2,nint)*rl**3 + c(1,nint)*ru**3 else c long range v = -17.57d0/r**6 endif else c short range v = 5757.005316d0*dexp(-5.792467035d0*r) - asym endif c write(21,*)r,v c ucm = 42048.3 + v*219474.4354 return end c subroutine splicn(x,y,m,c,ypo,ypm) implicit double precision(a-h,o-z) c c spline function subroutine c set of m pts,abscissae in x,ordinates in y,initial deriv yp0, c final deriv ypm, coeffs left in c max of 200 points c dimension x(m),y(m),d(200),p(200),e(200),c(4,m),a(200,3), * b(200),z(200) c mm = m-1 do 2 k = 1,mm d(k) = x(k+1) - x(k) p(k) = d(k)/6.0d0 e(k) = (y(k+1) - y(k))/d(k) 2 continue p(m) = 0.0d0 b(1) = (e(1) - ypo)/(2.0d0*p(1)) do 3 k = 2,mm b(k) = e(k) - e(k-1) 3 continue c b(m) = ypm - e(m-1) a(1,3) = 0.5d0 do 4 k = 2,m a(k,2) = 2.0d0*(p(k-1) + p(k)) - p(k-1)*a(k-1,3) b(k) = b(k) - p(k-1)*b(k-1) a(k,3) = p(k)/a(k,2) b(k) = b(k)/a(k,2) 4 continue c z(m) = b(m) do 6 i = 1,mm k = m-i z(k) = b(k) - a(k,3)*z(k+1) 6 continue c do 7 k = 1,mm q = 1.0d0/(6.0d0*d(k)) c(1,k) = z(k)*q c(2,k) = z(k+1)*q c(3,k) = y(k)/d(k) - z(k)*p(k) c(4,k) = y(k+1)/d(k) - z(k+1)*p(k) 7 continue return end