!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!System Information!!!!!!!!!! !!!!PKH3 Reactive Scattering!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! $Title_Labels Title1='Positron-Hydrogen reactive scattering Exact coupled-channel result' Title2='channel(1)=Positron+H, channel(2)=Positronium+Protron, channel(3)=closed' Title3='H3 Porter-Karplus PES' $END Title_Labels h-pos -------------------------------------------------------------------------- Main option parameters: LMESHER logical variable set equal to .true. if mesher is to be executed. LSFUNC logical variable set equal to .true. if sfunc is to be executed. LMATELM logical variable set equal to .true. if matelm is to be executed. LOVERLAP logical variable set equal to .true. if overlap is to be executed. LPLOT logical variable set equal to .true. for plotting wavefunctions and potential. LAPH3D logical variable set equal to .true. if aph3d is to be executed. LXSEC logical variable set equal to .true. if xsec is to be executed. CRAY logical set equal to .true. if running on a Cray -------------------------------------------------------------------------- $options lmesher=.true.,lsfunc=.true.,lmatelm=.true.,loverlap=.true., cray=.false.,lenlvls=.true. ,scheme=2 $end ----------------------------------------------------------------------- debug printing options nsubs ...... number of subroutines which will have their debug writes turned on The following variables must be defined for each subroutine for which you wish to turn on debug writes: subs(i) ...... name of the subroutine to turn on debug writes; must be a Hollerith character string of no more than eight characters iprt(i) ...... what level of debug writes: 0=none 1=minimum 2=moderate 3=maximum 4=combination (must read in iter(j,i), j=1 to 3) iter(j,i) ...... how many times to loop thru subroutine i with the jth debug level turned on (j=1 to 3) ----------------------------------------------------------------------- $debug nsubs=2 subs(1)=' debug',iprt(1)=3, subs(2)=' adini',iprt(2)=3, subs(3)=' adress',iprt(3)=3 $end -------------------------------------------------------------------------- Read in the masses (carbon 12 mass units) for each atom. MASSES: massa=mass(Lithium), massb=mass(Flourine), massc=mass(Hydrogen) -------------------------------------------------------------------------- $mass amass=5.4858026e-4,bmass=1.007276470,cmass=5.4858026e-4 $end ----------------------------------------------------------------------- Specify the following: JTOT total angular momentum IPAR parity (IPAR=0=even-parity, IPAR=1=odd-parity) NSYMC is true only if the diatom in that channel is homonuclear. NSYMC can be true for only one arrangement channel. ----------------------------------------------------------------------- $momentum jtot=0,parity=0,nsymc(1)=.false.,.false.,.false. $end ----------------------------------------------------------------------- c Set the zero point of the potential (in Hartree), if needed. c Set vzero=15.2eV for H+(e+) ----------------------------------------------------------------------- $setpot vzero=0.558585145 $end $savnam ithsav=500 $end ----------------------------------------------------------------------- Read in the mesher data: rhofirst Value of the first rho for calculating surface functions. (units of Bohr) rhomax Value of the last rho for calculating surface functions. (units of Bohr) deltarho Length of the first sector. (unit of Bohr) ithrho If ithrho is less than or equal to 1 a new calculation will be started. If ithrho is greater than 1 it is assumed that ithrho steps have already been calculated and the calculations will start at the ithrho+1 step. ----------------------------------------------------------------------- $sfundist rhomin=9.999,rhomax=10.0,deltarho1=0.001 iaccn=1,iset=.false.,nconsidr=16,ovrtol=0.00001,rhowrit=50.000, rhoend=37.2 $end ----------------------------------------------------------------------- NLT initial number of elements in the theta direction. NLC initial number of elements in the chi direction for each arrangement channel. ENGEV maximum total energy (units of eV) for scattering. NDIV number of times each element can be subdivided. NFREQ number of eigenenergies desired. NQ number of basis vectors used (must be larger than NFREQ). NODMAX maximum number of nodal points desired. ----------------------------------------------------------------------- $mshdat nlt=8,nlc(1)=8,8,16,ndiv=0,nfreq=16, nodmax=4450,tolfn1=0.000025,tolfn2=0.000025,xkmesh=1.012035,nlancz=2, nblock=4,maxj=200,tolexp=0.865,nxtra=27 $end $stst betaw(1)=0.85,betaw(2)=0.55,betaw(3)=30.0, nquad=13,seq(1)=0.0,seq(2)=0.0,seq(3)=0.0, chncof(1)=1.0,chncof(2)=1.0,chncof(3)=0.0, betafrst(1)=0.425,betafrst(2)=0.275,betafrst(3)=15.0 $end ----------------------------------------------------------------------- Read in the propagation data: EFIRST total energy for the first propagation DELTAENG energy increment NENERGY number of total energies --------------------------------------------------------------------------- $energy emin=0.3734800000d0,ejump=1.837451135d-3,numengs=1 $end c-----------------------------------------------------------------------readi197 c If INTEGRAT(i) = true, then channel i is propagated. readi198 c-----------------------------------------------------------------------readi199 $intchanl integrat(1)=.true.,.true.,.false. $end c-----------------------------------------------------------------------readi201 c If IRDINDEP = true, read the energy independent propagation info from readi202 c disk on the first energy. If false, do not read from readi203 c disk on the first energy. readi204 c If IWRINDEP = true, write the energy independent propagation info ontoreadi205 c disk on the first energy. If false, do not write readi206 c onto disk for the first energy. readi207 c-----------------------------------------------------------------------readi208 $store irdindep=.false., iwrindep=.true. $end --------------------------------------------------------------------------- Read in data to determine the asymptotic vibrational-rotational states of the isolated diatoms. MINVIB the minimum vibrational state for each arrangement channel. MAXVIB the maximum vibrational state for each arrangement channel. JMIN the minimum rotational state for each vibrational state and each arrangement channel. JMAX the maximum rotational state for each vibrational state and each arrangement channel. --------------------------------------------------------------------------- $quantum minvib(1)=1,1,0,maxvib(1)=2,1,0, jmin(1,1)=4*0,jmax(1,1)=1,0,0,0, jmin(1,2)=4*0,jmax(1,2)=0,0,0, jmin(1,3)=4*0,jmax(1,3)=4*0, $end --------------------------------------------------------------------------- Read in data necessary to generate the asymptotic basis and for expanding the interaction potential so that interaction matrix elements can be calculated analytically. NOSCIL number of harmonic oscillators basis functions. NGLEGN number of Gauss-Legendre points used to expand the interaction potential in a Legendre polynomial. NHERMT numer of Gauss-Hermite points used in calculating the Hamiltonian matrix elements. RE equilibrium position of the diatom for each arrangement channel. RX RX times RE is the position of the minimum of the parabola which defines the harmonic oscillator basis. This parameter can usually be set equal to 1.1 for each arrangement channel. This parameter should be greater than one because the long range part of the potential is softer than the short range part of the potential. RALFA RALFA times ALFA is curvature of the parabola which defines the harmonic oscillator basis. This number should be slightly less than one so that the harmonic oscillator basis will have amplitude over a large enough range to adequadely expand the exact wavefunction. This parameter can usually be set equal to 0.95 for each arrangement channel. NPOW the order of the polynomial used to expand the vibrational dependence of the interaction potential. The variable x=(r-re)/re is used as the independent variable. SCHEME if 1, use harmonic oscillators for asymptotic analysis if 2, use Hylleras functions for asymptotic analysis --------------------------------------------------------------------------- $gauss noscil(1)=10,10,10,nglegn(1)=6,6,6,nhermt(1)=6,6,6, ralpha(1)=1.0002722,1.0002722,0.0, rx(1)=0.00,0.00,0.00,npow(1)=10,10,10, re(1)=1.00000,re(2)=1.413829,re(3)=0.00000 $end --------------------------------------------------------------------------- Read in the spectroscopic constants to define the experimental vibrational rotational states of the isolated diatoms. These constants are obtained from Herzberg's book Spectra of Diatomic Molecules and have units of inverse centimeters. These are used to insure that one has calculated accurate asymptotic wavefunctions. The calculated energies are used and not the 'experimental' energies. --------------------------------------------------------------------------- $spectro wecm(1)=4464.788,wexecm(1)=130.17,becm(1)=60.82,alfecm(1)=2.301, weyecm(1)=0.0,decm(1)=0.0d0, wecm(2)=4464.788,wexecm(2)=130.17,becm(2)=60.82,alfecm(2)=2.301, weyecm(2)=0.0,decm(2)=0.0d0, wecm(3)=4464.788,wexecm(3)=130.17,becm(3)=60.82,alfecm(3)=2.301, weyecm(3)=0.0000,decm(3)=0.0d0 $end $etachanl eta(1)=0.48,0.53,0.0, scalagr=0.87 $end --------------------------------------------------------------------------- Read in the number of APH channels used in the propagation and the value of lambda (projection of the total angular momentum on the body-fixed z-axis) for each channel. --------------------------------------------------------------------------- $aphinp naph=15,lambda(1)=4*0 $end --------------------------------------------------------------------------- Read in the propagation distances. START starting distance used in the propagation. ENDAPH distance which signals the end of the APH region. ENDDELVE distance which signals the end of the Delves region. SWITCH distance to switch from the Log-derivative propagator to the VIVS propagator. FINISH maximum propagation distance. Distance at which asymptotic functions are used to obtain the S-matrix (scattering matrix). NOTE the following inequalities must be satisfied: START .le. ENDAPH .le. ENDDELVE .le. FINISH and START .le. SWITCH .le. FINISH The following schematic may be useful. Log-derivative VIVS --------------------- ------------------------- | | | | | start endaph switch enddelve finish --------- ------------------------ ------------ APH Delves Jacobi distance -----> --------------------------------------------------------------------------- $distance start=4.637245201715,endaph=7.104530507108,enddelve=0.5, switch=190.0,finish=0.30490 $end --------------------------------------------------------------------------- Read in nstep used to control the accuracy in the log-derivative method. NSTEPS number of steps per asymptotic wavelength at the highest total energy. NSTEP=20 is usually adequate and NSTEP=40 usually gives very accurate S-matrix elements. For energies less than the higest total energy the scattering wavefunction will be calculated more accurately. --------------------------------------------------------------------------- $logder nsteps=40 $end --------------------------------------------------------------------------- Read in parameters used in VIVS for control of the accuracy in the propagation. DRNOW initial interval size. TOFF tolerance parameter which requires that the perturbation corrections are less than this number. This parameter is roughly the accuracy of the resulting S-matrix elements at the lowest energy. S-matrix elements at higher energies are probably calculated more accurately because the perturbation corrections are smaller. IAPLHA number of step per interval. Usually set equal to 4. --------------------------------------------------------------------------- $vivs drnow=0.05,toff=.032,ialpha=4 $end --------------------------------------------------------------------------- Read in contour plotting parameters. For plotting contours only (not needed otherwise) --------------------------------------------------------------------------- $contour nx=100,xi=-1.0,dx=0.020202,ny=100,yi=-1.0, dy=0.020202,iaph=9 $end --------------------------------------------------------------------------- Read in 3-D plotting parameters. For plotting 3-dimensional graphs using grafic (not used otherwise) --------------------------------------------------------------------------- $plot irhostrt=1,irhoend=1,ipotntl=.false.,nsurf=2, nthmode(1)=1,2,3,4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,20, 21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40, 41,42,43,44,45,46,47,48,49,50,logpot=.true.,denplot=.false., iplt=2,nl=750,nil=30,kolb(1)=7,7,view(1)=0.,-50.,50.,cutoff=2100, win(1)=0.0,0.0,0.0,4.0,iproj=1,irdsurf=.false.,npltskp=0 $end