!------------------------------------------------------------------- ! function : syspot2 ! ! package : li3 Potentials ! ! Language : Fortran 77 ! ! author : V. Venturi (vanessa.venturi@nist.gov) ! ! revision : F. Colavecchia (flavioc@lanl.gov) ! ! purpose : Get the potential of the singlet and triplet Li2. ! ! ! input : r -> interatomic distance in bohr. ! nv -> singlet (nv=1) or triplet (nv=2) potential ! lsetup -> sets up the reading of data or calculation ! of potential ! ! output : v -> li2 potential in a.u. ! ! ! modules : ! ! common : ! ! notes : The original code included RKR data from table 3 ! (Linton et al J. Mol. Spect. 196, 20 (1999)). ! Now we make use of new RKR data provided by A. Ross ! (ross@in2p3.fr), read from rkr.dat. See also ! readme_li2.txt for the values of De, Do, a, rmin. ! Comments noted by CFC were included by F. Colavecchia ! !------------------------------------------------------------------- DOUBLE PRECISION FUNCTION SYS_SS_pot_3(R,V,NV,LSETUP) IMPLICIT NONE DOUBLE PRECISION R,V INTEGER NV LOGICAL LSETUP C C ... Local Variables C double precision c6,c8,c10,CA,beta,gamma,as,bs,at,bt, > re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat save c6,c8,c10,CA,beta,gamma,as,bs,at,bt, > re1,re3,shift1,shift3,rdatas,vdatas,rdatat,vdatat double precision hartree,bohra save hartree,bohra double precision correction_ss double precision rs1,vs1,rs2,vs2,rs3,vs3,rt0,vt0,rt1,vt1,rt2,vt2, > rt3,vt3,res,des,ret,det,rdat,vdats,vdatt double precision vs,vt,vav,vdiff,vex,vdisp,amp,vshift, > rfit,rjoin,rfinal,r0,rl,ri,r1,r2,v1,v2,delta integer np,ns,nt,ndatas,ndatat save np,ns,nt,ndatas,ndatat integer i,j,k,ni,nsl,ntl,ndat,nt2off C integer ns1,ns2,ns3,nt1,nt2,nt3,nt0 integer nmax parameter (nmax=500) parameter (ns1=3,ns2=42,ns3=41,nt1=11,nt0=11,nt2=27,nt3=26) C real*8 vdatatp(nmax),vdatasp(nmax) dimension rs1(ns1),vs1(ns1),rs2(ns2),vs2(ns2),rs3(ns3),vs3(ns3), > rt1(nt1),vt1(nt1),rt2(nt2),vt2(nt2),rt3(nt3),vt3(nt3) dimension rdatas(nmax),vdatas(nmax),rdatat(nmax),vdatat(nmax), > rdat(nmax),vdats(nmax),vdatt(nmax) dimension vt0(nt0),rt0(nt0) logical do_setup data do_setup /.true./ save do_setup logical order,norepeat C C output flags for potential routine: C output_data - output raw potential data C output_fitex - exchange fitting constant A as a function of R C output_v - output final potentials to be used in calculations C output_vmesh - output interpolated mesh of final potentials + C correction included C logical output_data, output_fitex, output_v, output_vmesh data output_data /.false./ data output_fitex /.false./ data output_v /.false./ data output_vmesh /.false./ C C=====Data for ^7Li2==================================================== C C Physical constants used in potential routines C DATA HARTREE /219474.630680D0/, BOHRA /0.5291772490D0/ C C np is the degree of polynomial used in the akima spline interpolation C np=3 C C rfit= point at which the exchange form is fit to the data and at C which is calculated the shift in the data required to match the C short and long range parts C rjoin= point at which the short range data is smoothly connected C to the long range form C rfinal= last point for short range data set used in the calculations C DATA rfit /15.5D0/, rjoin /17.D0/,rfinal /20.0d0/ C C --- Long range: C Yan, Babb, Dalgarno and Drake, Phys. Rev. A 54, 2824 (1996) C (in atomic units) C DATA C6 /1393.39D0/, C8 /83425.8D0/, C10 /7372100.0D0/ C C --- Long range exchange: C Cote, Dalgarno, and Jamieson, Phys. Rev. A 50, 399 (1994) C (in atomic units) C DATA beta /1.25902D0/, gamma /4.55988D0/ C C --- Short-range singlet potential data C ---------------------------------- C De determined by Zemke and Stwalley, J. Phys. Chem. 97, 2053 (1993) C (in units of cm-1) C DATA Des /8516.61D0/ C C Section with E>0 : small R C 1st 2 pts from Konowalow and Olson J. Chem. Phys. 71, 450 (1979) C 3rd pt from Schmidt-Mink et al Chem. Phys. 92, 263 (1985) C (in atomic units) C DATA rs1 / 2.75, 3.00, 3.25 / C DATA vs1 / 0.061866D0, 0.036411D0, 0.013546D0 / C C Section with E<0 : RKR potential C from Barakat et al Chem. Phys. 102, 215 (1986) C note: misprints in data corrected, see Cote et al C Phys. Rev A 50 399 (1994) C (in units of angstroms and cm-1) C DATA res / 2.673247D0 / C DATA rs2 / 1.823201D0, 1.823462D0, 1.824057D0, & 1.825118D0, 1.826727D0, 1.828918D0, & 1.831708D0, 1.834867D0, 1.838851D0, & 1.843462D0, 1.848639D0, 1.854376D0, & 1.860716D0, 1.867431D0, 1.874904D0, & 1.882864D0, 1.891368D0, 1.900497D0, & 1.910173D0, 1.920282D0, 1.931173D0, & 1.942717D0, 1.954984D0, 1.967947D0, & 1.981710D0, 1.996354D0, 2.011906D0, & 2.028431D0, 2.046073D0, 2.064963D0, & 2.085263D0, 2.107177D0, 2.130967D0, & 2.156977D0, 2.185675D0, 2.217724D0, & 2.254114D0, 2.296439D0, 2.347580D0, & 2.413867D0, 2.517207D0, 2.673247D0 / C DATA vs2 / 8514.7777D0, 8507.8421D0, 8492.0440D0, & 8463.9645D0, 8421.6125D0, 8364.3068D0, & 8292.0294D0, 8205.2323D0, 8104.4730D0, & 7990.4162D0, 7863.7083D0, 7724.9765D0, & 7574.8736D0, 7413.8431D0, 7242.5556D0, & 7061.4199D0, 6870.8931D0, 6671.3979D0, & 6463.3140D0, 6246.9482D0, 6022.6578D0, & 5790.7056D0, 5551.3992D0, 5304.9322D0, & 5051.5343D0, 4791.4274D0, 4524.7756D0, & 4251.7309D0, 3972.4624D0, 3687.1094D0, & 3395.7978D0, 3090.6412D0, 2795.7419D0, & 2487.1914D0, 2173.0721D0, 1853.4573D0, & 1528.4128D0, 1197.9974D0, 862.2642D0, & 521.2611D0, 175.0320D0, 0.0000D0 / C DATA rs3 / 12.639382D0, 9.863785D0, 8.495553D0, & 7.657218D0, 7.090991D0, 6.675227D0, & 6.352341D0, 6.089904D0, 5.870023D0, & 5.680896D0, 5.514781D0, 5.366746D0, CX & 5.580896D0, 5.514781D0, 5.366746D0, & 5.232506D0, 5.109851D0, 4.996566D0, & 4.890860D0, 4.791595D0, 4.697788D0, CX & 4.490860D0, 4.791595D0, 4.697788D0, & 4.608426D0, 4.522846D0, 4.440679D0, & 4.361370D0, 4.284400D0, 4.209387D0, & 4.135091D0, 4.064160D0, 3.993263D0, & 3.923163D0, 3.853594D0, 3.784295D0, & 3.714965D0, 3.645290D0, 3.574908D0, & 3.503385D0, 3.430172D0, 3.354534D0, & 3.275416D0, 3.191166D0, 3.098852D0, & 2.992096D0, 2.848951D0 / C DATA vs3 / 8514.7777D0, 8507.8421D0, 8492.0440D0, & 8463.9645D0, 8421.6125D0, 8364.3068D0, & 8292.0294D0, 8205.2323D0, 8104.4730D0, & 7990.4162D0, 7863.7083D0, 7724.9765D0, & 7574.8736D0, 7413.8431D0, 7242.5556D0, & 7061.4199D0, 6870.8931D0, 6671.3979D0, & 6463.3140D0, 6246.9482D0, 6022.6578D0, & 5790.7056D0, 5551.3992D0, 5304.9322D0, & 5051.5343D0, 4791.4274D0, 4524.7756D0, & 4251.7309D0, 3972.4624D0, 3687.1094D0, & 3395.7978D0, 3090.6412D0, 2795.7419D0, & 2487.1914D0, 2173.0721D0, 1853.4573D0, & 1528.4128D0, 1197.9974D0, 862.2642D0, & 521.2611D0, 175.0320D0 / C C --- Short-range triplet potential data C ---------------------------------- C De determined by Linton et al, J. Mol. Spect. 196, 20 (1999) C (in units of cm-1) C DATA Det /333.70D0/ C C Section with E>0 : small R C 1st pt from Konowalow and Olson J. Chem. Phys. 81, 4534 (1984) C 7 pts from Schmidt-Mink et al Chem. Phys. 92, 263 (1985) C (in atomic units) C c DATA rt1 / 3.00D0 , 3.25D0, 3.50D0, c & 4.00D0 , 4.50D0, 5.00D0, c & 5.50D0 , 6.00D0/ C c DATA vt1 / 0.104471D0, 0.082074D0, 0.063725D0, c & 0.036778D0, 0.020006D0, 0.010053D0, c & 0.004370D0, 0.001258D0 / C C 11 pts from Halls et al, Chem. Phys. Lett. 339, 427 (2001) C distances in angstroms c DATA rt1 / 2.00d0, 2.125d0, 2.25D0, 2.375D0, & 2.50d0, 2.625d0, 2.75d0, & 2.875d0, 3.00d0, 3.125d0, 3.25d0/ DATA vt1 / 0.0472751972d0, & 0.0361788060d0, 0.0273239836d0, 0.0203412490d0, & 0.0148869076d0, 0.0106605276d0, 0.0074111640d0, & 0.0049337380d0, 0.0030628654d0, 0.0016656470d0, & 0.0006356372d0 / C C Section with E<0 : RKR potential C from Linton et al J. Mol. Spect. 196, 20 (1999) C (in units of angstroms and cm-1) C DATA ret / 4.173D0 / C CFC DATA rt2 / 3.356D0, 3.358D0, 3.365D0, CFC & 3.377D0, 3.395D0, 3.422D0, CFC & 3.458D0, 3.505D0, 3.571D0, CFC & 3.668D0, 3.846D0, 4.159D0,4.173D0 / C CFC DATA vt2 / 333.269D0, 330.170D0, 322.155D0, CFC & 308.098D0, 287.665D0, 260.837D0, CFC & 227.679D0, 188.240D0, 142.523D0, CFC & 90.453D0, 31.857D0, 0.024D0, 0.000D0 / C CFC DATA rt3 / 16.052D0, 11.392D0, 9.441D0, CFC & 8.297D0, 7.501D0, 6.885D0, CFC & 6.373D0, 5.922D0, 5.503D0, CFC & 5.092D0, 4.630D0 , 4.185D0/ C CFC DATA vt3 / 333.269D0, 330.170D0, 322.155D0, CFC & 308.098D0, 287.665D0, 260.837D0, CFC & 227.679D0, 188.240D0, 142.523D0, CFC & 90.453D0, 31.857D0, 0.024D0 / C C=============Initial setup============================================= C IF(LSETUP) THEN if (do_setup) then shift1=0.0d0 shift3=0.13d-5 c6=1393.39d0 CFC ! ! Read new RKR points ! open(unit=1200,file='./Potentials/li3web/rkr.dat', +status='old') do i=1,nt2 read(1200,*) rt2(i),vt2(i) !write(*,*) i,rt2(i),vt2(i) end do do i=1,nt3 read(1200,*) rt3(i),vt3(i) !write(*,*) i,rt3(i),vt3(i) end do close(1200) nt2off = 4 ! skip the first nt2off values of rkr data C C --- Calculate total number of potential data points C ns=ns1+ns2+ns3 nt=nt1+nt2+nt3-nt2off if (ns.gt.nmax.or.nt.gt.nmax) then write(*,1000) ns,nt,nmax endif C C --- Order each set of data, in ascending order, into an array C and smooth data by multiplying by R^6 for interpolation. C Convert data to atomic units. C re1=res/bohra re3=ret/bohra C do i=1,ns1 rdatas(i)=rs1(i) vdatas(i)=vs1(i) vdatas(i)=vdatas(i)*rdatas(i)**6 enddo do i=1,ns2 j=ns1+i rdatas(j)=rs2(i)/bohra vdatas(j)=(vs2(i)-des)/hartree vdatas(j)=vdatas(j)*rdatas(j)**6 enddo do i=1,ns3 j=ns1+ns2+i rdatas(j)=rs3(i)/bohra vdatas(j)=(vs3(i)-des)/hartree vdatas(j)=vdatas(j)*rdatas(j)**6 enddo C order=.true. do i=2,ns if (rdatas(i).eq.rdatas(i-1)) then write(*,1010) i-1,i stop elseif (rdatas(i).lt.rdatas(i-1)) then order=.false. endif enddo if (.not. order) then call sort (rdatas,vdatas,ns) endif C c do i=1,nt0 c rdatat(i)=rt0(i)/bohra c vdatat(i)=vt0(i) c vdatat(i)=vdatat(i)*rdatat(i)**6 c enddo do i=1,nt1 j=i rdatat(j)=rt1(i)/bohra vdatat(j)=vt1(i) vdatat(j)=vdatat(j)*rdatat(j)**6 enddo C C Note: 1st and 2nd points of the RKR data had to be omitted C to prevent oscillations in the spline fit between C the two sets of data CFC CFC Now we omit the first nt2off points of the left side CFC of RKR data. CFC do i=nt2off+1,nt2 j=nt1-nt2off+i rdatat(j)=rt2(i)/bohra vdatat(j)=(vt2(i)-det)/hartree vdatat(j)=vdatat(j)*rdatat(j)**6 enddo do i=1,nt3 j=nt1+nt2-nt2off+i rdatat(j)=rt3(i)/bohra vdatat(j)=(vt3(i)-det)/hartree vdatat(j)=vdatat(j)*rdatat(j)**6 enddo C order=.false. do i=2,nt if (rdatat(i).eq.rdatat(i-1)) then write(*,*) 'Error in RKR data' write(*,1010) i-1,i write(*,'(2f20.8)') rdatat(i-1),vdatat(i-1) write(*,'(2f20.8)') rdatat(i),vdatat(i) stop elseif (rdatat(i).lt.rdatat(i-1)) then order=.false. endif enddo if (.not. order) then call sort (rdatat,vdatat,nt) endif C if (output_data) then do i=1,ns write(50,1100) rdatas(i),vdatas(i)/rdatas(i)**6 enddo do i=1,nt write(51,1100) rdatat(i),vdatat(i)/rdatat(i)**6 enddo C do i=1,ns rdat(i)=rdatas(i) vdats(i)=vdatas(i) do j=1,nt if (rdatas(i).eq.rdatat(j)) then vdatt(i)=vdatat(j) else call akima(np,nt,rdatat,vdatat,1,rdat(i),vt) vdatt(i)=vt endif enddo enddo C k=ns do i=1,nt norepeat=.true. do j=1,ns if (rdatat(i).eq.rdatas(j)) then norepeat=.false. endif enddo if (norepeat) then k=k+1 rdat(k)=rdatat(i) vdatt(k)=vdatat(i) call akima(np,ns,rdatas,vdatas,1,rdat(k),vs) vdats(k)=vs endif enddo ndat=k call sort2 (rdat,vdats,vdatt,ndat) C do i=1,ndat write(52,1110) rdat(i),vdats(i),vdatt(i) enddo C endif C C --- Fit Smirnov and Chibisov exchange form: A exp(-beta R) R^gamma C to the exchange energy determined from the potential data to C obtain the amplitude A C C Calculate A at the chosen fitting point C C Interpolate potential data using akima spline interpolation C call akima(np,ns,rdatas,vdatas,1,rfit,vs) call akima(np,nt,rdatat,vdatat,1,rfit,vt) C C Determine fitting constant A C vt=vt/rfit**6 vs=vs/rfit**6 vex=log((vt-vs)/2.0d0) CA=vex+beta*rfit-gamma*log(rfit) CA=exp(CA) write(*,*) CA C C --- Calculate the average of vs and vt, which should correspond C with vdisp C vav=(vt+vs)/2.0d0 vdisp=-c6/rfit**6-c8/rfit**8-c10/rfit**10 vshift=vdisp-vav write(*,*) vshift write(*,*) vdisp write(*,*) exp(vex) C C --- Shift short range data by calculated difference to match the C long range form C ! vshift = 0d0 nsl=0 do i=1,ns ! vdatas(i)=vdatas(i)+vshift*rdatas(i)**6 if(nsl.eq.0.and.rdatas(i).ge.rfit) nsl=i enddo ntl=0 do i=1,nt ! vdatat(i)=vdatat(i)+vshift*rdatat(i)**6 if(ntl.eq.0.and.rdatat(i).ge.rfit) ntl=i enddo if(nsl.eq.0.or.ntl.eq.0) then write(*,1020) rfit stop endif C C --- Connect short range data onto long range form using akima spline C to smooth across the connection C ! call akima(np,ns,rdatas,vdatas,1,rfit,vs) ! call akima(np,nt,rdatat,vdatat,1,rfit,vt) rdatas(nsl)=rfit rdatat(ntl)=rfit vdatas(nsl)=(-(vt-vs)/2.0d0+vdisp)*rfit**6 vdatat(ntl)=( (vt-vs)/2.0d0+vdisp)*rfit**6 write(*,*) 'ntl =',ntl write(*,*) 'nt =',nt ! vdatas(nsl)=vs ! vdatat(ntl)=vt C delta = 0.5d0 ni=int((rfinal-rjoin)/delta)+1 ! ni = 2 ndatas=nsl+ni ndatat=ntl+ni ri=rjoin do i=1,ni vex=CA*exp(-beta*ri)*ri**gamma vdisp=-c6/ri**6-c8/ri**8-c10/ri**10 C rdatas(nsl+i)=ri rdatat(ntl+i)=ri vdatas(nsl+i)=(vdisp-vex)*ri**6 vdatat(ntl+i)=(vdisp+vex)*ri**6 C ri=ri+delta enddo rfinal=rdatas(ndatas) C C Extrapolate short range data using the form: aexp(-br) C r1=rdatas(1) r2=rdatas(1)+1d-6 call akima(np,ndatas,rdatas,vdatas,1,r1,v1) call akima(np,ndatas,rdatas,vdatas,1,r2,v2) if (v1*v2.lt.0.0d0) then write(*,1030) stop endif v1=v1/r1**6 v2=v2/r2**6 bs=log(v1/v2)/(r2-r1) as=v2*exp(bs*r2) C r1=rdatat(1) r2=rdatat(1)+1d-6 !call akima(np,ndatat,rdatat,vdatat,1,r1,v1) call akima(np,ndatat,rdatat,vdatat,1,r2,v2) if (v1*v2.lt.0.0d0) then write(*,1030) stop endif v1=vdatat(1)/r1**6 v2=v2/r2**6 bt=log(v1/v2)/(r2-r1) at=v2*exp(bt*r2) C C Output potential details/parameters C write(*,1040) c6,c8,c10 write(*,1050) rfit,CA,beta,gamma write(*,1060) vshift,vshift*hartree write(*,1070) rjoin,rfinal write(*,1080) as,bs,at,bt endif C do_setup=.not.do_setup SYS_SS_pot_3=0.0d0 C RETURN ENDIF C C======================================================================= C IF((NV.LE.0).OR.(NV.GT.2)) THEN WRITE(*,900) STOP END IF C GO TO (10,20),NV C C======================================================================= C C --- SECTION FOR NV=1 -- Calculate singlet potential C 10 call mpot(np,ndatas,rdatas,vdatas,c6,c8,c10, > CA,beta,gamma,as,bs,nv,r,v) C v=v+correction_ss(r,shift1,re1) SYS_SS_pot_3=V C RETURN C C======================================================================= C C --- SECTION FOR NV=2 -- Calculate triplet potential C 20 call mpot(np,ndatat,rdatat,vdatat,c6,c8,c10, > CA,beta,gamma,at,bt,nv,r,v) C v=v+correction_ss(r,shift3,re3) SYS_SS_pot_3=V C RETURN C C======================================================================= C C --- FORMAT STATEMENTS ------------------------------------------------ C 900 FORMAT(//,' Program HALTED in Potential Subroutine -- NV OUT ', X 'OF RANGE') C 990 FORMAT(//,'MOLECULAR POTENTIAL DATA') 1000 format(//,'ERROR in Potential Subroutine -- no. of data points', X 'exceeds array size:',3i6) 1010 format(//,'ERROR in Potential Subroutine -- identical data for', X 'pair:',2i6) 1020 format(//,'ERROR in Potential Subroutine -- rfit (au)=',d14.6, X 'exceeds short range data') 1030 format(//,'ERROR in Potential Subroutine -- data points used in', X 'short range extrapolation of data cross a zero') 1040 format(//,'Long range dispersion coefficients (au):',/, X ' C6=',d14.6,' C8=',d14.6,' C10=',d14.6) 1050 format(//,'Long range exchange, Aexp(-beta R)R^gamma (au):',/, X ' fit to short range data at R=',d14.6,/, X ' A=',d14.6,' beta=',d14.6,' gamma=',d14.6) 1060 format(/,'Shift applied to short range data to match to', X ' long range form:',/, X d14.6,' au',' or',d14.6,' cm-1') 1070 format(//,'Long range form assumed at R (au)=',d14.6,/, X 'Short range data extended to R (au)=',d14.6) 1080 format(//,'Short range extrapolation, aexp(-bR) (au)',/, X 'Singlet: a=',d14.6,' b=',d14.6,/, X 'Triplet: a=',d14.6,' b=',d14.6) 1100 format(2d16.8) 1110 format(3d16.8) 1120 format(4d16.8) C END