c---------------------------------------------------------------------------- c c DJ - Pseudo-diatomic model for a linear dimer A...B of two c linear molecules A and B C c---------------------------------------------------------------------------- C C Input is from file DJ.INP C Output is to file DJ.OUT C C C C Version 21.X.2016 ----- Zbigniew KISIEL ----- C C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C c Sources of formulae: c C [BCKF] T.J.Balle, E.J.Campbell, M.R.Keenan, W.H.Flygare, C "A new method for observing the rotational spectra of weak C molecular complexes", C J.Chem.Phys. 72 (1980) 922-932. c C [RCH] W.G.Read, E.J.Campbell, G.Henderson, C "The rotational spectrum and molecular structure of the C benzene-hydrogen chloride complex", C J.Chem.Phys. 78 (1983) 3501-3508 c c [Millen] D.J.Millen, c "Determination of stretching force constants of weakly bound c dimers from centrifugal distortion constants", c Canad.J.Chem. 7 (1985) 1477-1479 c c [KC] S.G.Kukolich, E.J.Campbell, c "Microwave measurements of bromine quadrupole coupling constants c and the molecular structure of XeHBr", C Chem.Phys.Lett. 94 (1983) 73-76 c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c C Modification history: C C 4.11.01: stable version C 4.04.02: cyclic operation C 19.07.07: commenting in preparationj for web upload C 2.08.16: major facelift C 21.10.16: more facelifting/equations/testing C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c IMPLICIT REAL*8(A-H,O-Z) PARAMETER (CONV=505379.01D0, * convcm=29979.2458d0, * convd=57.2957795d0) character line*150,name*30,errmes*100 real*8 ma,mb,mu,ia,ib,iab,ie,ibab,icab,iea C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | DJ - Pseudo-diatomic model calculations for', * T79,'|'/ * ' | weakly bound dimers', * T79,'|'/ * ' |',76(1H_),'|'/' version 21.X.2016',T64,'Zbigniew KISIEL'/) c c open(1,file='dj.inp',status='old',err=1000) open(2,file='dj.out',status='unknown') c c c...Input c C !----------------------------------------------------------------------------------------------------------------------------------------- C ! Name of complex a..b B/MHz C/MHz DJ/kHz Ma/u Mb/u Ba/MHz Bb/MHz Theta.a Theta.b C !_____________________________---------------_______________----------____________------------____________------------__________---------- C ! C ! C benzene...H35Cl 1237.68362 1237.68362 1.223 78.046950 35.9766762 5687.7 312989.297 0.0 22.96 C 1 read(1,'(a)',err=8,end=8)line if(line(1:1).eq.'!')goto 1 c do 20 n=31,45 if(line(n:n).eq.'-')goto 8 20 continue c if(len_trim(line).lt.109)goto 1 c read(line,2,err=8)name,bab,cab,dj,ma,mb,ba,bb,thetaa,thetab 2 format(a30,2f15.8,f10.5,4f12.6,2f10.5) if(b.lt.0.0d0)goto 8 c if(bab.le.0.d0)then errmes='Bad value of B.ab' goto 1010 endif if(cab.le.0.d0)then errmes='Bad value of C.ab' goto 1010 endif if(ba.le.0.d0)then errmes='Bad value of B.a' goto 1010 endif if(bb.le.0.d0)then errmes='Bad value of B.b' goto 1010 endif if(bb.le.0.d0)then errmes='Bad value of B.ab' goto 1010 endif if(ma.le.0.d0)then errmes='Bad value of M.a' goto 1010 endif if(mb.le.0.d0)then errmes='Bad value of M.b' goto 1010 endif c C C..Reduced mass (u) and moments of inertia (u Angstr^2) C if(bab.eq.cab)then b=bab else b=0.5d0*(bab+cab) ibab=conv/bab icab=conv/cab endif ia=conv/ba ib=conv/bb iab=conv/b c mu=mb*ma/(mb+ma) C C...Echo of input C write(2,3)name 3 format(79(1h-)/1x,a/79(1h-)/) write(2,'(1x)') write(2,15)'INPUT DATA:' write(2,'(1x)') write(2,15)' Monomer A '// * ' Monomer B' write(2,'(1x)') write(2,15)' B.a /MHz ',ba, * ' B.b /MHz ',bb write(2,15)' Ib.a /uA^2 ',ia, * ' Ib.b /uA^2 ',ib write(2,15)' M.a /u ',ma, * ' M.b /u ',mb write(2,15)' Theta.a /deg ',thetaa, * ' Theta.b /deg ',thetab write(2,'(1x)') write(2,15)' Dimer AB' write(2,'(1x)') c if(bab.eq.cab)then write(2,15)' B.ab /MHz ',b write(2,15)' Ib.ab /uA^2 ',iab else write(2,15)' B.ab /MHz ',bab, * ' Ib.ab /uA^2 ',ibab write(2,15)' C.ab /MHz ',cab, * ' Ic.ab /uA^2 ',icab write(2,15)' (B+C)/2 ',b, * ' I_(B+C)/2 ',iab endif c write(2,15)' mu /u ',mu write(2,15)' DJ.ab /kHz ',dj c 15 format(a,f16.7,7x,a,f16.7,10x,a,f16.7) C C C...Straight diatomic intermolecular (stretching) force constant C fs_diat = 2.6221D-10 * B**3 * mu / DJ C C C..Stretching force constant (Millen) for C 1/ linear complex of two linear molecules eq.12 C 2/ symmetric top complex of symmetric+linear molecules eq.16 C (eq.12 and eq.16 are identical) C C fs = f.sigma /Nm-1 from eq.(12)/(16) of Millen where C for B in MHz, DJ in kHz, mu in u: C C 2.6222181D-10 = 16*Pi^2*10^3*10^12*1.660539040*10^(-27) C fs_millen=2.6222181D-10 * B**3 * mu/dj * (1.d0-b/ba-b/bb) C C C..Stretching force constant (Millen) for planar asymmetric top dimer C with C2v symmetry Eq.(21) C C bsmall=bab/ba+bab/bb csmall=cab/ba+cab/bb fs_millen21=1.3111091D-10 * mu/dj * * (bab**3 * (1-bsmall) + cab**3 * (1-csmall) - * 0.25d0*(Bab-Cab)**2*(Bab+Cab)*(2.d0-bsmall-csmall) ) C C C..Stretching force constant (Millen) for asymmetric top dimer with C hydrogen bond perpendicular to planar monomer Eq.(25) C C bsmall=bab/ba+bab/bb csmall=cab/ba+cab/bb fs_millen25=1.3111091D-10 * mu/dj * * (bab**3 * (1-bsmall) + cab**3 * (1-csmall)) C C c Diatomic stretching frequency: c c oms = omega.sigma /cm-1 from eq.(14) of [BCKF] c where for B in MHz, DJ in kHz: C C 2.109645*D-3 = sqrt(4*10^3) 10^6 / 2.99792458*10^10 C c oms= 2.109645D-3*dsqrt(B**3/dj) omsa=1.30277D+2*dsqrt(fs_millen/mu) C c c Equilibrium rotational constant Be /MHz is obtained by c solving eq.(21) of [BCKF]: c c B0 = Be - 18 ( Be^2/om), c c for which the quadratic formula solution (Be, om, Be all in MHz) is c c 2Be = om/18 - sqrt( (om/18)^2 - 4(om/18) B0 ) c c This solution is actually given explicitly in the last equation of [KC] c om= convcm*(oms) be=0.5d0* ( (om/18.d0) - dsqrt( (om/18.d0)**2 * - 4.d0 * (om/18.d0) * b) ) ie=conv/be c oma=convcm*(omsa) bea=0.5d0*( (oma/18.d0) - dsqrt( (oma/18.d0)**2 * - 4.d0 * (oma/18.d0) * b) ) iea=conv/bea c c...diatomic distances c r0=dsqrt( conv/(mu*b) ) re=dsqrt( conv/(mu*be) ) rea=dsqrt( conv/(mu*bea) ) c costa=dcos(thetaa/convd) costb=dcos(thetab/convd) c c...Centre of mass separation: a general expression is somewhat elusive to locate C but it seems to begin with eq.(2) of [RCH] where one of the clustered C molecules has a correction for large amplitude motion. c C The formula for two rods, both with angular large amplitude C motion is given, for example, for OC...HCN in eq.(6) of C E.J.Goodwin,A.C.Legon, Chem.Phys. 87 (1984)81-92 C C The formula for a symmetric top and a rod, both with angular large amplitude C motion is given, for example, for H3P...HCN, in eq.(4) of C A.C.Legon, L.C.Willoughby, Chem.Phys. 85(1984) 443-450. c c For no oscillations the basic formula reduces to: C Ib_dimer = Ib_A + Ib_B + mu rcm^2 C which is also eq.(17) of [Millen] c rcm=dsqrt( * (conv/b-0.5d0*ia*(1.d0+costa**2)-0.5d0*ib*(1.d0+costb**2)) * /mu ) rcmno=dsqrt( * (conv/b-0.5d0*ia*(1.d0+1.d0 )-0.5d0*ib*(1.d0+1.d0 )) * /mu ) C C C Stretching force constant for asymmetric complex consisting of C rigid rod and rigid body where rod is colinear to one of c inertial axes of the body, eq.(6) of [RCH] C C fs = 16 pi^4 (mu rcm)^2 [ (B^2+C^2)^2 + 2 (B^4+C^4) / (h DJ) C C For B,C in MHz, DJ in kHz, rcm in Angstr, mu in uAnstr^2: C C fs = 6.485771 * 10^(-17) * C (mu rcm)^2 [ (B^2+C^2)^2 + 2 (B^4+C^4) / DJ) c C where 6.485777*10^(-17) = 16 pi^4 u^2 10^(-20) 10^24/(h 10^3) C = 16 pi^4 1.66053905*10^(-19)/6.62607004 C fs_RCH = 6.485771D-17 * * (mu*rcm)**2 * ((BAB**2+CAB**2)**2 + * 2*(BAB**4+CAB**4)) / DJ c c c eps = Lennard-Jones well depth (epsilon) /cm-1 from eq.(19) [BCKF] where c for fs in N/m, re in Angstr: c c 503.4117 = 10^(-20) / ( 6.62607004*10^(-34) * 2.99792458*10^10 ) c eps_diat= 503.4117d0*fs_diat *re**2/72.d0 eps_mill= 503.4117d0*fs_millen *re**2/72.d0 eps_asym= 503.4117d0*fs_rch *re**2/72.d0 eps_mill25= 503.4117d0*fs_millen25 *re**2/72.d0 c c Output c write(*,'('' Processing: '',a)')name 7 format(1x,a,f15.6,3x,a,f15.6,3x,a,f15.6) c write(2,'(1x)') write(2,'(1x)') write(2,15)'CALCULATED:' write(2,'(1x)') write(2,115)' R.cm_0/Angstr',r0, ' diatomic, from '// * 'Ib_ab = mu R.cm_0^2' write(2,115)' R.cm/Angstr ',rcmno,' from Ib_ab '// * '= Ib_a + Ib_b + mu R_cm^2, e.g. eq.17 [Millen]' if(thetaa.ne.0.d0.or.thetab.ne.0.d0) * write(2,115)' R.cm/Angstr ',rcm,' as above but with '// * '(1+cos^2 theta) corrections' c write(2,'(1x)') write(2,115)' f.sigma/Nm-1 ',fs_diat, * ' diatomic from B^3 mu / DJ' write(2,115)' f.sigma/Nm-1 ',fs_millen, * ' linear/symmetric complex, eq.12=16 [Millen]' c if(bab.ne.cab)then write(2,115)' f.sigma/Nm-1 ',fs_RCH, * ' using B.ab and C.ab, eq.6 [RCH]' write(2,115)' f.sigma/Nm-1 ',fs_millen21, * ' using B.ab and C.ab, eq.21 [Millen], but monomer B=C' write(2,115)' f.sigma/Nm-1 ',fs_millen25, * ' using B.ab and C.ab, eq.25 [Millen], but monomer B=C' endif c 115 format(a,f14.5,a) write(2,'(1x)') write(2,115)' omega_s/cm-1 ',oms, * ' diatomic from DJ=4B^3/omega^2, eq.14 [BCKF]' write(2,115)' omega_s/cm-1 ',omsa, * ' from 1/(2pi)sqrt(f.sigma/mu), eq.15 [BCKF]' write(2,'(32x,a)')'and using f.sigma from eq.12=16 [Millen]' c write(2,'(1x)') write(2,115)' B.ab_e/MHz ',be, ' diatomic, last eq. [KC]' write(2,115)' I.ab_e/uA^2 ',ie, ' from B.ab_e' write(2,115)' R.cm_e/Angstr',re, ' diatomic, eq.22a [BCKF]' c write(2,'(1x)') write(2,115)' B.ab_e/MHz ',bea, * ' last eq. [KC] with omega_s '// * 'eq.15 [BCKF} + eq.12=16 [Millen]' write(2,115)' I.ab_e/uA^2 ',iea,' from B.ab_e above' write(2,115)' R.cm_e/Angstr',rea, * ' diatomic, eq.22a [BCKF] but from I.ab_e above' c write(2,'(1x)') write(2,115)' eps/cm-1 ',eps_diat, * ' eq.19 [BCKF] using diatomic f.sigma' write(2,115)' eps/cm-1 ',eps_mill, * ' eq.19 [BCKF] using f.sigma from eq.12=16 [Millen]' if(bab.ne.cab)then write(2,115)' eps/cm-1 ',eps_asym, * ' eq.19 [BCKF] using asymmetric f.sigma from eq.6 [RCH]' write(2,115)' eps/cm-1 ',eps_mill25, * ' eq.19 [BCKF] using asymmetric f.sigma from eq.25 [Millen]' endif write(2,'(1x)') c write(2,'(1x//)') goto 1 C 8 write(2,'(79(1H-)/)') c write(2,1001) * '[BCKF] T.J.Balle, E.J.Campbell, M.R.Keenan, W.H.Flygare,', * ' A new method for observing the rotational spectra'// * ' of weak', * ' molecular complexes,', * ' J.Chem.Phys. 72 (1980) 922-932.', * ' ', * '[RCH] W.G.Read, E.J.Campbell, G.Henderson,', * ' The rotational spectrum and molecular structure'// * ' of the', * ' benzene-hydrogen chloride complex,', * ' J.Chem.Phys. 78 (1983) 3501-3508.', * ' ', * '[Millen] D.J.Millen,', * ' Determination of stretching force constants of'// * ' weakly bound', * ' dimers from centrifugal distortion constants",', * ' Canad.J.Chem. 7 (1985) 1477-1479.', * ' ', * '[KC] S.G.Kukolich, E.J.Campbell,', * ' Microwave measurements of bromine quadrupole '// * 'coupling constants', * ' and the molecular structure of XeHBr",', * ' Chem.Phys.Lett. 94 (1983) 73-76.' 1001 format(50(1x,a/)) write(2,'(79(1H-))') c close(1) close(2) write(*,'(1x//)') C STOP C 1000 write(*,'(1x/5x,a//)') * '**** ERROR: Cannot open the input file DJ.INP' stop c 1010 write(*,'(1x/5x,a,'' when reading from''/a//)') * errmes(1:len_trim(errmes)), * line(1:len_trim(line)) c C STOP END c------------------------------------------------------------------------ c------------------------------------------------------------------------