*dk surface
subroutine surface(v,rij)
use pes_module
c-----------------------------------------------------------------------
c This routine supplies the potential energy surface for the
c reactive scattering codes.
c
c On entering: The arrary rij contains the internuclear distances
c in units of bohr (NOT mass scaled).
c
c The common /pes/ must contain a character string
c of length 25 that properly describes the potential.
c
c On Exit: The variable v is the value of the potential energy
c surface in Hartree atomic units.
c
c
c-----------------------------------------------------------------------
!implicit none
implicit real*8 (a-h,o-z)
!real*8 v, rij
dimension rij(3)
character(LEN=25) potname
!common /pes/ pesname
c-----------------------------------------------------------------------
c The character string potname must contain a string that uniquely
c specifies the potential energy surface to be used.
c-----------------------------------------------------------------------
data potname/'DMBE'/
c-----------------------------------------------------------------------
c Check for correct potential energy surface specification.
c-----------------------------------------------------------------------
if(pesname.ne.potname)then
write(6,*)'The potential energy surface name does not match',
> ' the name in the potential energy surface routine.'
write(6,304)pesname,potname
304 format(1h ,2a25)
write(6,*)'***error***: Execution stopped in routine surface'
stop
end if
c-----------------------------------------------------------------------
c This call to surfac supplies the potential energy surface
c for the particular system of interest.
c-----------------------------------------------------------------------
call dmbe(v,rij)
return
end
*dk dmbe
subroutine dmbe(v,r1)
implicit real*8 (a-h,o-z)
dimension r1(3)
common/potcm/r(3),pote1,dedr(3)
common/pot2cm/pote2,dedr2(3),coup(3)
r(1)=r1(1)
r(2)=r1(2)
r(3)=r1(3)
call h3dmbe
v=pote1+.1744741119d0
return
end
subroutine h3dmbe
c***********************************************************************
c calculates double many-body expansion of the h3 potential
c pote1 is energy of surface 1.
c dedr contains derivatives for surface 1.
c pote2 is energy of surface 2.
c dedr2 contains derivatives for surface 2.
c coup contains nonadiabatic coupling.
c note: dedr,dedr2,coup not yet coded.
c***********************************************************************
c refs.: b. liu & p. siegbahn, j. chem. phys., 68, 2457, 1978;
c d.g. truhlar & c.j. horowitz, j. chem. phys., 68, 2466, 1978;
c a.j.c. varandas, mol. phys., 53, 1303, 1984;
c a.j.c. varandas, j. mol. struct. (theochem), 21, 401, 1985;
c a.j.c. varandas, d.g. truhlar, c.a. mead, & n.c. blais, to be
c published.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /acdcom/alph2,az2,ce0,cd1
common /crdcom/ r12,r22,r32,qcoord,rho,s,cphi3,s2,s3
common /h2com/rm,r0,rmt
common /fixcom/alpha5,beta1,beta2,beta3
common /picom/pi,sqrt3
common /potcm/r(3),pote1,dedr(3)
common /pot2cm/ pote2,dedr2(3),coup(3)
common /tercom/ vhf,hflinc,bendt,hfadd,vhfc,ce2,ce3,corr
common /lepcom/ vleps2,dlep(3),dlep2(3)
common /xcocom/x(15)
common /dampp/damps(10,3,2)
common /corr/corr2s(3,2)
external h3data
data rhol /2.4848d0/
r1=r(1)
r2=r(2)
r3=r(3)
call coord(r1,r2,r3)
per=r1+r2+r3
per2=per*per
per3=per2*per
exp1=exp(-beta1*per)
exp2=exp(-beta2*per2)
exp3=exp(-beta3*per)
hflinc=0.0d0
do 1 j=1,4
hflinc=hflinc+x(j)*alin(r1,r2,r3)**(2*j)
1 continue
hflinc=hflinc*exp(-alpha5*per3)
c construct damping function matrix damps(power,distance, state)
call damp(r)
c construct 2-body correlation fcn matrix corr2s(distance, state)
call corr2(r)
vhf=vleps(r1,r2,r3)
b1=b1fn(r1,r2,r3)
b12=b1**2
b13=b12*b1
b14=b13*b1
bend1=b1*(x(5)+x(6)*per)*exp1
bend2=(x(7)*b12+x(8)*b13+x(9)*b14)*exp2
bend3=(1.d0/r1+1.d0/r2+1.d0/r3)*(x(10)*b1*exp1+x(11)*
1 b12*exp2)
bend4=((r2-r1)**2+(r3-r1)**2+(r3-r2)**2)*
1 b1*(x(12)*exp1+x(13)*exp2)
bend5=b1*(x(14)+x(15)*per2)*exp3
bendt=bend1+bend2+bend3+bend4+bend5
11 hfadd=s2*(1.d0+s3*cphi3)*(ce0+cd1*rho)*exp(-alph2*
* (rho-rhol)**2)
45 vhfc=vhf+hflinc+bend1+bend2+bend3+bend4+bend5+hfadd
call corr3(r1,r2,r3,rm,r0,ce3)
ce2=corr2s(1,1)+corr2s(2,1)+corr2s(3,1)
corr=ce2+ce3
pote1=vhfc+corr
pote2=pote1-vhf+vleps2
return
end
subroutine prepot
implicit real*8 (a-h,o-z)
c dummy subprogram
return
end
function alin(r1,r2,r3)
c***********************************************************************
c computes base function for the second fit to the hartree-fock
c energy at the linear symmetrical geometries.
c***********************************************************************
implicit real*8 (a-h,o-z)
alin=(r1-r2)*(r2-r3)*(r3-r1)
return
end
function b1fn(r1,r2,r3)
c***********************************************************************
c computes b1 function.
c***********************************************************************
implicit real*8 (a-h,o-z)
t1=(r1**2-r2**2-r3**2)/(2.d0*r2*r3)
t2=(r2**2-r1**2-r3**2)/(2.d0*r1*r3)
t3=(r3**2-r1**2-r2**2)/(2.d0*r1*r2)
b1fn=1.d0+t1+t2+t3
return
end
subroutine corr2(r)
c***********************************************************************
c calculates 2-body correlation energy.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /noncom/chh(10),chhe(10),chh2(10)
common /dampp/damps(10,3,2)
common /corr/corr2s(3,2)
dimension ri(10),r(3)
do 1 id=1,3
ri(2)=1.0d0/(r(id)*r(id))
ri(4)=ri(2)*ri(2)
ri(6)=ri(4)*ri(2)
ri(8)=ri(6)*ri(2)
ri(10)=ri(8)*ri(2)
do 1 is=1,2
corr2s(id,is)=0.0d0
do 1 i=6,10,2
corr2s(id,is)=corr2s(id,is)-damps(i,id,is)*chh(i)*ri(i)
1 continue
return
end
function corr2d(r,rm,r0)
c***********************************************************************
c calculates 1st derivative of the 2-body correlation energy.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /noncom/chh(10),chhe(10),chh2(10)
corr2d=0.0d0
do 1 i=6,10,2
corr2d=corr2d-(chh(i)/r**i)*(ddamp(r,i,rm,r0)-float(i)*
* damp(r,i,rm,r0)/r)
1 continue
return
end
subroutine damp(r)
c***********************************************************************
c calculates matrix damps(power,dist,state) of dispersion damping
c functions.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /discom/alph0,alph1,bet0,bet1
common /dampp/damps(10,3,2)
common /h2com/rm,r0,rmt
dimension r(3), a(10),b(10)
do 1 i=6,10,2
a(i)=alph0/float(i)**alph1
b(i)=bet0*exp(-bet1*float(i))
1 continue
do 2 is=1,2
xm=rm
if(is.eq.2)xm=rmt
do 2 id=1,3
x=2.0d0*r(id)/(xm+2.5d0*r0)
xs=x*x
do 2 n=6,10,2
pol=a(n)*x+b(n)*xs
2 damps(n,id,is)=(1.0d0-exp(-pol))**n
return
end
function ddamp(r,n,rm,r0)
c***********************************************************************
c calculates 1st derivative of damping function.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /discom/alph0,alph1,bet0,bet1
a=alph0/float(n)**alph1
b=bet0*exp(-bet1*float(n))
denom=rm+2.5d0*r0
x=2.0d0*r/denom
pol=a*x+b*x**2
dpol=a+2.d0*b*x
ddamp=(float(2*n)/denom)*dpol*exp(-pol)*(1.d0-exp(-pol))**(n-1)
return
end
function g(r,ck0,ck1)
c***********************************************************************
c calculates g function.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /h2com/rm,r0,rmt
g=1.d0+ck0*exp(-ck1*(r-rm))
return
end
function gd(r,ck0,ck1)
c***********************************************************************
c calculates 1st derivative of g function.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /h2com/rm,r0,rmt
gd=-ck1*ck0*exp(-ck1*(r-rm))
return
end
function h(r,et)
c***********************************************************************
c calculates h function.
c***********************************************************************
implicit real*8 (a-h,o-z)
h=(tanh(et*r))**6
return
end
function hd(r,et)
c***********************************************************************
c calculates 1st derivative of h function.
c***********************************************************************
implicit real*8 (a-h,o-z)
hd=6.0d0*tanh(et*r)**5*et/cosh(et*r)**2
return
end
function h2hf(r)
c***********************************************************************
c calculates extended-hartree-fock curve for ground singlet state
c of h2 [a.j.c. varandas & j.d. silva, j. chem. soc. faraday ii
c (submitted)].
c***********************************************************************
implicit real*8 (a-h,o-z)
common /diacom/d,a1,a2,a3,gamma
common /h2com/rm,r0,rmt
dr=r-rm
h2hf=-d*(1.0d0+a1*dr+a2*dr**2+a3*dr**3)*exp(-gamma*dr)
return
end
function h2hfd(r)
c***********************************************************************
c calculates 1st derivative of 2-body extended-hartree-fock curve
c for the ground-singlet state of h2.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /diacom/d,a1,a2,a3,gamma
common /h2com/rm,r0,rmt
dr=r-rm
h2hfd=-gamma*h2hf(r)-d*(a1+2*a2*dr+3.0d0*a3*dr**2)*exp(-gamma*dr)
return
end
function switch(r1,r2,r3)
c***********************************************************************
c calculates switching term for the hartree-fock component of the diatomic
c triplet state function. the notation of thompson et al. (j. chem. phys.,
c 82,5597,1985; and references therein) is used throughout.
c***********************************************************************
implicit real*8 (a-h,o-z)
common/acdcom/alph2,az2,ce0,cd1
common /picom/ pi,sqrt3
common /crdcom/ r12,r22,r32,qcoord,rho,s,cphi3,s2,s3
common /swfcom/sf
if(s.ne.0.0d0)go to 1
switch=exp(-az2*qcoord**2)
go to 11
1 switch=exp(-az2*(1.d0+s3*cphi3)*qcoord**2)
11 sf=switch
return
end
function trip(r,rsquar,al0,al1,al2,al3)
c***********************************************************************
c calculates effective diatomic triplet state curve.
c***********************************************************************
implicit real*8 (a-h,o-z)
trip=(al0+al1*r+al2*rsquar)*exp(-al3*r)
return
end
function vleps(r1,r2,r3)
c***********************************************************************
c calculates 3-body extended-hartree-fock energy defined by a
c leps-type function.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /crdcom/ r12,r22,r32,qcoord,rho,s,cphi3,s2,s3
common /h2com/rm,r0,rmt
common /trncom/al0,al1,al2,al3
common /lepcom/ vleps2,dlep(3),dlep2(3)
common /corr/corr2s(3,2)
c vleps2 is leps upper surface.
c dlep(3) are leps lower surface derivatives (not yet coded)
c dlep2(3) " " upper " " " " "
sng1=h2hf(r1)
sng2=h2hf(r2)
sng3=h2hf(r3)
if(1)11,1,11
11 f=switch(r1,r2,r3)
call acct(r1,hft1,cet1,at1,1)
cesng1=corr2s(1,1)
t1=trip(r1,r12,al0,al1,al2,al3)*f+(1.d0-f)*(at1-cesng1)
call acct(r2,hft2,cet2,at2,2)
cesng2=corr2s(2,1)
t2=trip(r2,r22,al0,al1,al2,al3)*f+(1.d0-f)*(at2-cesng2)
call acct(r3,hft3,cet3,at3,3)
cesng3=corr2s(3,1)
t3=trip(r3,r32,al0,al1,al2,al3)*f+(1.d0-f)*(at3-cesng3)
go to 2
1 t1=trip(r1,r12,al0,al1,al2,al3)
t2=trip(r2,r22,al0,al1,al2,al3)
t3=trip(r3,r32,al0,al1,al2,al3)
2 q1=0.5d0*(sng1+t1)
q2=0.5d0*(sng2+t2)
q3=0.5d0*(sng3+t3)
ex1=0.5d0*(sng1-t1)
ex2=0.5d0*(sng2-t2)
ex3=0.5d0*(sng3-t3)
f1=(ex1-ex2)**2
f2=(ex2-ex3)**2
f3=(ex3-ex1)**2
qsum=q1+q2+q3
exch=sqrt(0.5d0*(f1+f2+f3))
vleps=qsum-exch
vleps2=qsum+exch
return
end
subroutine acct(r,hft,cet,t,i)
c***********************************************************************
c computes the hface (hartree-fock-approximate correlation energy)
c potential for the lowest triplet state of h2 [a.j.c. varandas &
c j. brandao mol. phys.,45,1982,857].
c***********************************************************************
implicit real*8 (a-h,o-z)
common /h2com/rm,r0,rmt
common /corr/corr2s(3,2)
data a,b1,b2,b3,b4/0.448467d0,-0.056687d0,0.246831d0,-0.018419d0,
* 0.000598d0/
r2=r*r
r3=r2*r
r4=r3*r
hft=a*exp(-(b1*r+b2*r2+b3*r3+b4*r4))/r
cet=corr2s(i,2)
t=hft+cet
return
end
subroutine coord(r1,r2,r3)
c**********************************************************************
c calculates d3h symmetry coordinates
c**********************************************************************
implicit real*8 (a-h,o-z)
common /crdcom/ r12,r22,r32,qcoord,rho,s,cphi3,s2,s3
common /picom/pi,sqrt3
r12=r1**2
r22=r2**2
r32=r3**2
qcoord=r12+r22+r32
rho=sqrt(qcoord/3.d0)
gamma=2.d0*r12-r22-r32
beta=sqrt3*(r22-r32)
s2=(beta**2+gamma**2)/qcoord**2
s=sqrt(s2)
s3=s2*s
if(s.ne.0.d0)go to 11
c if(s.eq.0.0d0) arg,theta,phi3,cphi3 are not computed.
go to 45
11 arg=gamma/qcoord/s
if(arg.le.-1.d0)arg=-1.d0
if(arg.ge.1.d0)arg=1.d0
theta=acos(arg)
phi3=3.d0*(pi-theta)
cphi3=cos(phi3)
45 return
end
subroutine corr3(r1,r2,r3,rm,r0,ce3)
c***********************************************************************
c calculates 3-body correlation energy.
c***********************************************************************
implicit real*8 (a-h,o-z)
dimension t1(10),t2(10),t3(10),co3(10)
common /coecom/ck0(10),ck1(10)
common /noncom/chh(10),chhe(10),chh2(10)
common /dampp/damps(10,3,2)
c define two statement functions
g(x,c0,c1)=1.d0+c0*exp(-c1*(x-rm))
h(x,c1)=(tanh(x*c1))**6
do 1 i=6,10,2
g1=g(r1,ck0(i),ck1(i))
g2=g(r2,ck0(i),ck1(i))
g3=g(r3,ck0(i),ck1(i))
h1=h(r1,ck1(i))
h2=h(r2,ck1(i))
h3=h(r3,ck1(i))
t1(i)=chh(i)*(1.0d0-0.5d0*(g2*h3+g3*h2))*damps(i,1,1)/r1**i
t2(i)=chh(i)*(1.0d0-0.5d0*(g3*h1+g1*h3))*damps(i,2,1)/r2**i
t3(i)=chh(i)*(1.0d0-0.5d0*(g1*h2+g2*h1))*damps(i,3,1)/r3**i
co3(i)=t1(i)+t2(i)+t3(i)
1 continue
ce3=co3(6)+co3(8)+co3(10)
return
end
subroutine corr3d(r1,r2,r3,rm,r0,dc3r1,dc3r2,dc3r3)
c***********************************************************************
c calculates first-order partial derivatives of 3-body correlation
c energy with respect to interatomic distances.
c***********************************************************************
implicit real*8 (a-h,o-z)
dimension d1(10),d2(10),d3(10)
common /coecom/ck0(10),ck1(10)
common /noncom/chh(10),chhe(10),chh2(10)
dc3r1=0.d0
dc3r2=0.d0
dc3r3=0.d0
do 1 i=6,10,2
g1=g(r1,ck0(i),ck1(i))
g2=g(r2,ck0(i),ck1(i))
g3=g(r3,ck0(i),ck1(i))
gd1=gd(r1,ck0(i),ck1(i))
gd2=gd(r2,ck0(i),ck1(i))
gd3=gd(r3,ck0(i),ck1(i))
h1=h(r1,ck1(i))
h2=h(r2,ck1(i))
h3=h(r3,ck1(i))
hd1=hd(r1,ck1(i))
hd2=hd(r2,ck1(i))
hd3=hd(r3,ck1(i))
dd1=ddamp(r1,i,rm,r0)
dd2=ddamp(r2,i,rm,r0)
dd3=ddamp(r3,i,rm,r0)
gh23=g2*h3+g3*h2
gh31=g3*h1+g1*h3
gh12=g1*h2+g2*h1
d1(i)=chh(i)*dd1*(1.d0-0.5d0*gh23)/r1**i-
1 float(i)*chh(i)*damp(r1,i,rm,r0)*(1.d0-0.5d0*gh23
2 )/r1**(i+1)+
3 chh(i)*damp(r2,i,rm,r0)*(-0.5d0*(g3*hd1+gd1*h3))/r2**i+
4 chh(i)*damp(r3,i,rm,r0)*(-0.5d0*(gd1*h2+g2*hd1))/r3**i
d2(i)=chh(i)*dd2*(1.d0-0.5d0*gh31)/r2**i-
1 float(i)*chh(i)*damp(r2,i,rm,r0)*(1.d0-0.5d0*gh31
2 )/r2**(i+1)+
3 chh(i)*damp(r1,i,rm,r0)*(-0.5d0*(gd2*h3+g3*hd2))/r1**i+
4 chh(i)*damp(r3,i,rm,r0)*(-0.5d0*(g1*hd2+h1*gd2))/r3**i
d3(i)=chh(i)*dd3*(1.d0-0.5d0*gh12)/r3**i-
1 float(i)*chh(i)*damp(r3,i,rm,r0)*(1.d0-0.5d0*gh12
2 )/r3**(i+1)+
3 chh(i)*damp(r1,i,rm,r0)*(-0.5d0*(g2*hd3+gd3*h2))/r1**i+
4 chh(i)*damp(r2,i,rm,r0)*(-0.5d0*(gd3*h1+g1*hd3))/r2**i
dc3r1=dc3r1+d1(i)
dc3r2=dc3r2+d2(i)
dc3r3=dc3r3+d3(i)
1 continue
return
end
block data h3data
c***********************************************************************
c data for the dmbe h3 surface.
c***********************************************************************
implicit real*8 (a-h,o-z)
common /acdcom/alph2,az2,ce0,cd1
common /barcom/rsp,vspkc
common /diacom/d,a1,a2,a3,gamma
common /discom/alph0,alph1,bet0,bet1
common /fixcom/alpha5,beta1,beta2,beta3
common /h2com/rm,r0,rmt
common /noncom/chh(10),chhe(10),chh2(10)
common /coecom/ck0(10),ck1(10)
common /picom/pi,sqrt3
common /trncom/al0,al1,al2,al3
common /xcocom/x(15)
data alph0,alph1,bet0,bet1/25.9528d0,1.1868d0,15.7381d0,0.09729d0/
data alpha5,beta1,beta2,beta3/ 8.2433033d-03,
2 0.53302897d+00,0.39156612d-01,0.69996945d+00/
data alph2/4.735364d-1/
data al0,al1,al2,al3/
2 0.45024472d+01,-0.62467617d+01,0.40966542d+01,0.21813012d+01/
data az2/4.453649d-4/
data ce0,cd1/6.333404d-03,-1.726839d-03/
data chh(6),chh(8),chh(10)/6.499027d0,124.3991d0,3285.8d0/
data chhe(6),chhe(8),chhe(10)/2.82d0,41.8d0,865.6d0/
data chh2(6),chh2(8),chh2(10)/8.813d0,163.87d0,4023.2d0/
data ck0,ck1/ 0.000000d+00, 0.000000d+00,
* 0.000000d+00, 0.000000d+00, 0.000000d+00,
* -1.372843d-01, 0.000000d+00, -1.638459d-01,
* 0.000000d+00, -1.973814d-01, 0.000000d+00,
* 0.000000d+00, 0.000000d+00, 0.000000d+00,
* 0.000000d+00, 1.011204d+00, 0.000000d+00,
* 9.988099d-01, 0.000000d+00, 9.399411d-01/
data d,a1,a2,a3,gamma/0.15796326d0,2.1977034d0,1.2932502d0,
* 0.64375666d0,2.1835071d0/
data pi/3.14159265358979d0/
data rm,r0,rmt/1.401d0,6.928203d0,7.82d0/
data rsp,vspkc/1.757d0,9.65d0/
data sqrt3/1.73205080756887d0/
data x/-0.9286579d-02, 0.2811592d-03, -0.4665659d-05,
2 0.2069800d-07, 0.2903613d+02, -0.2934824d+01,
3 0.7181886d+00, -0.3753218d+00, -0.1114538d+01,
4 0.2134221d+01, -0.4164343d+00, 0.2022584d+00,
5 -0.4662687d-01, -0.4818623d+02, 0.2988468d+00/
end