program calchocl c************************************** c c Calculates HOCl X^1A' PES as described in: c c K.A. Peterson, S. Skokov, and J.M. Bowman, J. Chem. Phys. 111, 7446 (1999) c c this version corresponds to the "adjusted" PES c c************************************** implicit real*8 (a-h,o-z) c r1=1.8 r2=3.21 theta=102.8 e=V(r1,r2,theta) print *, r1,r2,theta,e c 1.8 3.21 102.8 -163.886520396843 - 0.23380 end c real*8 function V(r1in,r2in,thin) implicit double precision (a-h,o-z) parameter (hkcal=627.5096d0) data ncall/1/ c if(ncall.eq.1) then call loadp ncall=2 endif pi = 4.0 * atan (1.0) c x3=dacos(x3t)*180/pi c call jzmat(x1,x2,x3,r1in,r2in,thin) c call adjust(r1in,r2in,thin,r1,r2,theta) call calcpes(r1,r2,theta,e) c c call calcpes(r1in,r2in,thin,e) c ehocl=-57428.07/349.75504d0 ehclo=-38646.2/349.75504d0 c v=(e-ehocl)/hkcal v=e return end subroutine jzmat(x1,x2,x3,r1,r2,theta) C****Transformation from jacobian (x1,x2,x3) for HOCL** C****to z-matrix (r1,r2,theta) for HOCL ** C* x1 - R(H-OCl), x2 - R(OCL), x3 - jacobian angle * C* r1 - R(OH), r2 - R(OCL), theta - HOCL angle * c* distance in bohr, angle in degree * C****************************************************** implicit double precision (a-h,o-z) amcl=34.9688527 amo =15.99491463 pi = 4.0 * atan (1.0) x3=x3*pi/180.0d0 dx1=x2*amcl/(amcl+amo) dx2=x2*amo/(amcl+amo) r2=x2 r1=sqrt(dx1**2+x1**2-2.0*dx1*x1*cos(x3)) ar=(dx1**2+r1**2-x1**2)/(2*dx1*r1) if (ar.gt.1.0) ar=1.0 if (ar.lt.-1.0) ar=-1.0 theta=180.0d0*acos(ar)/pi x3=x3*180.0d0/pi return end subroutine adjust(r1in,r2in,thin,r1,r2,theta) C**r1in, r2in,thin - input values C**r1,r2,theta - corrected values implicit double precision (a-h,o-z) parameter (hkcal=627.5096d0) pi = 4.0 * atan (1.0) C**ADJUSTMENTS************************************* C**Shift ground state to experimental one********** r1exp=1.822395 !exper. ground state r2exp=3.191692 thexp=102.965 r1e=1.81830 !minimum of current PES r2e=3.1977 r3e=2.49 !HCl in HClO roh=4.36 !HO in HClO the=102.821 th1e=108.0 dr1=(r1exp-r1e) dr2=r2exp-r2e dth=thexp-the r1=r1in-dr1 r2=r2in-dr2 theta=thin-dth C**Scaling of HOCl coordinates**************** sc=425-25*exp(-4*(abs(r1-r1exp)**3.5)) if (r1.gt.3.0001d0) sc=sc-900*exp(-0.1/(r1-3.d0)**3) sc1=1.00000d0-1.d-5*sc sc2=1.00d0+0.0018*exp(-2.d0*abs(r1-r1exp)) sc3=0.9914d0+0.004*exp(-3.d0*abs(r1-r1exp)) r1=sc1*(r1-r1exp)+r1exp r2=sc2*(r2-r2exp)+r2exp theta=sc3*(theta-thexp)+thexp return end subroutine calcpes(r1,r2,theta,e) implicit real*8 (a-h,o-z) r3=sqrt(r1*r1+r2*r2-2.d0*r1*r2*cosd(theta)) c if(r1.lt.0.8.or.r2.lt.0.8.or.r3.lt.0.8) then e=500.0d0 return endif c y2a=func_ap(1,r1) !OH 2-body term y2b=func_ap(2,r2) !ClO 2-body term y2c=func_mlj(3,r3) !HCl 2-body term y3=func_3b(4,r1,r2,r3) !HOCl 3-body term C**Change OCl on asymptote de=0.23380d0 if(r2.le.5.001) y2b=y2b-de if(r2.gt.5.001) y2b=y2b-de*(1.d0-exp(-1.d0/(r2-5.0)**2)) e=y2a+y2b+y2c+y3 return end c function func_mlj(nr,r) implicit double precision (a-h,o-z) parameter (maxnp=700) common/genpef/nset,p9(5,maxnp),indi(maxnp),indj(maxnp), > indk(maxnp),indl(maxnp),e1,mci,nlp c np=p9(nr,1) de=p9(nr,2) rnexp=p9(nr,3) ve=-de re=p9(nr,4) z=(r-re)/(r+re) beta=p9(nr,5) do j=1,np-4 beta=beta+p9(nr,j+5)*(z**j) enddo s=1.0d0 - (re/r)**rnexp *exp(-beta*z) func_mlj=ve + de*s*s return end function func_ap(nr,r) implicit double precision (a-h,o-z) parameter (maxnp=700) common/genpef/nset,p9(5,maxnp),indi(maxnp),indj(maxnp), > indk(maxnp),indl(maxnp),e1,mci,nlp c np=p9(nr,1) re=p9(nr,2) alpha=p9(nr,3) beta=p9(nr,4) dr=(r-re) rho=r*exp(-beta*dr) sum=0.d0 do j=1,np-4 sum=sum+p9(nr,j+4)*rho**(j) enddo func_ap=sum + abs(p9(nr,np+1))*exp(-alpha*dr)/r return end c function func_3b(ns,r1,r2,r3) implicit double precision (a-h,o-z) parameter (maxnp=700) common/genpef/nset,p9(5,maxnp),indi(maxnp),indj(maxnp), > indk(maxnp),indl(maxnp),e1,mci,nlp dimension rho(3),r(3) data c1,c2,c3/.38156d0,.46172d0,.33730d0/ data b5,b6,b7,b8,b9/-8.64605d0,18.01025d0,-10.44204d0,2.85753d0, > -0.338997d0/ nplin=p9(ns,1)-nlp r(1)=r1 r(2)=r2 r(3)=r3 r1e=p9(ns,2) r2e=p9(ns,3) r3e=p9(ns,4) b1=p9(ns,5) b2=p9(ns,6) b3=p9(ns,7) b4=p9(ns,8) dr1=r1-r1e dr2=r2-r2e dr3=r3-r3e c rho1=exp(-r1e/r1)*exp(-b1*dr1) rho2=exp(-r2e/r2)*exp(-b2*dr2) rho3=exp(-r3e/r3)*exp(-b3*dr3) bcos1=(r(3)*r(3)-r(1)*r(1)-r(2)*r(2))/(-2.d0*r(1)*r(2)) bcos2=(r(2)*r(2)-r(1)*r(1)-r(3)*r(3))/(-2.d0*r(1)*r(3)) bcos3=(r(1)*r(1)-r(2)*r(2)-r(3)*r(3))/(-2.d0*r(2)*r(3)) sumcos=1.d0+bcos1+bcos2+bcos3 theta1=acosd(bcos1) c sum=0.d0 do m=1,mci i=indi(m) j=indj(m) k=indk(m) l=indl(m) sum=sum+p9(ns,m+nlp+1) > *rho2**j *rho1**i *rho3**k *sumcos**l enddo c c conical intersection stuff (theta=180) ridge=b5 + r(2)*bcos1 + b6*r(1) + b7*r(1)**2 + b8*r(1)**3 + > b9*r(1)**4 c tans=21.d0 ci=sqrt(ridge**2 + (b4*b4*r(2)*sind(theta1))**2 + > 0.0001d0) > * switch((bcos1-cosd(160.0d0)),tans/3.d0) > * switch((r(2)-4.0d0),tans/3.d0) > * switch((3.0d0-r(2)),tans/3.d0) > * switch((r(1)-2.8d0),tans/3.d0) > * switch((1.4d0-r(1)),tans/3.d0) c c.. definition change rho1=exp(-c1*c1*dr1) rho2=exp(-c2*c2*dr2) rho3=exp(-c3*c3*dr3) c do m=mci+1,nplin i=indi(m) j=indj(m) k=indk(m) l=indl(m) sum=sum+p9(ns,m+nlp+1)*rho1**i * rho2**j * rho3**k * ci enddo func_3b=sum return end c function switch(x,a) implicit double precision (a-h,o-z) switch=0.5d0*(1.d0-tanh(a*x)) return end c subroutine loadp implicit double precision (a-h,o-z) parameter (maxnp=700) common/genpef/nset,p9(5,maxnp),indi(maxnp),indj(maxnp), > indk(maxnp),indl(maxnp),e1,mci,nlp c open(unit=11,file='hocl1.pot',status='OLD',form='FORMATTED') nset=4 do i=1,3 read(11,*) read(11,*) p9(i,1) read(11,*) (p9(i,j),j=2,p9(i,1)+1) enddo read(11,*) read(11,*) p9(4,1), nlp, mci read(11,*) (p9(4,i),i=2,nlp+1) do j=1,p9(4,1)-nlp read(11,*) indi(j),indj(j),indk(j),indl(j), > p9(4,j+nlp+1) enddo return end