program tabprnf implicit none integer*4 np0,ic,npoint,i,i1,iargc,nm,column,iu PARAMETER (np0=2000000) real*8 wsmax,wsmin,s(np0),t(np0),x,dw,dd,hw,xx,const,wo,w(np0), 1sf,temp,cm,sfc,wm,dd0 character*80 fn,s80 character*1 pt logical LG,lnm,lwe,la c lwe - weather to use constant bandwidth in cm-1 scale write(6,6009) 6009 format(' Usage: tabprnf name type npoints min max G T width',/,/, 1' Arguments: 1 name - name of the .tab file ',/, 1' 2 type - type of spectrum (column # in name)',/, 1' 3 npoints - number of spectral points',/, 1' 4 min - wmin',/, 1' 5 max -wmax',/, 1' 6 G - Gaussian peaks instead of Lorentzians',/, 1' 7 T - temperature',/, 1' 8 width - full width at half height') lnm=.false. if(iargc().lt.1)then write(6,*)'.TAB filename:' read(5,500)fn 500 format(a80) else call getarg(1,fn) write(6,*)'generating spectra from' write(6,*)fn endif if(iargc().gt.1)then call getarg(2,s80) read(s80,*)ic if(ic.lt.0)then ic=-ic lnm=.true. write(6,*)'cm-1 converted to nm !!!' endif write(6,*)ic,' column after frequency' if(ic.eq.1 )write(6,*)'absorption' if(ic.eq.2 )write(6,*)'circular dichroism' if(ic.eq.3 )write(6,*)'Raman scattering' if(ic.eq.8 )write(6,*)'Raman optical activity' if(ic.eq.20)write(6,*)'magnetic circular dichroism' if(ic.eq.30)write(6,*)'ORD from Rs in the second column' if(ic.eq.40)write(6,*)'NMR from shifts in the second column' else write(6,*) 1 'column to read (1 abs,2 vcd, 3 ram, 8 roa, 20 mcd, 30 ord,' 2//' 40 nmr):' read(5,*)ic endif if(iargc().gt.2)then call getarg(3,s80) read(s80,*)npoint write(6,*)npoint,' points' else write(6,*)'number of points:' read(5,*)npoint endif if(npoint.gt.np0)call report('too many spectral points') if(iargc().gt.3)then call getarg(4,s80) read(s80,*)wsmin call getarg(5,s80) read(s80,*)wsmax write(6,6000)wsmin,wsmax 6000 format(' frequency interval: ',f12.2,' - ',f12.2) else write(6,*)'Wmin:' read(5,*)wsmin write(6,*)'Wmax:' read(5,*)wsmax endif if(iargc().gt.5)then call getarg(6,pt) if(pt.eq.'g'.or.pt.eq.'G')then write(6,*)'Gaussian peaks' LG=.true. else write(6,*)'Lorentzian peaks' LG=.false. endif else write(6,*)'peak type (G/L):' read(5,501)pt 501 format(a1) if(pt.eq.'G'.or.pt.eq.'g')then LG=.true. else LG=.false. endif endif if(iargc().gt.6)then call getarg(7,s80) read(s80,*)temp write(6,6001)temp 6001 format(' Temperature: ',f5.1,' K') else write(6,*)'temperature (K):' read(5,*)temp endif if(iargc().gt.7)then call getarg(8,s80) read(s80,*)hw write(6,6101)hw 6101 format(' FWHH: ',f12.2) else write(6,*)'peak width (FWHH):' read(5,*)hw endif c band width in cm-1: la=.true. c spectrum in cm-1: lwe=.true. if(hw.lt.0)then hw=-hw if(iargc().gt.8)then call getarg(8,s80) read(s80,*)iu else 988 write(6,607) 607 format(' Units: 1: peak width in nm , spectrum in nm ',/, 1 ' 2: peak width in cm-1, spectrum in cm-1',/, 1 ' 3: peak width in cm-1, spectrum in nn ',/, 1 ' ---------------------------------------',/, 1 ' Input choice: ',$) read(5,*)iu endif if(iu.eq.1)then write(6,503)hw 503 format(' FWHH ',f10.3,' nm') else if(iu.eq.2)then lwe=.true. write(6,502)hw 502 format(' FWHH ',f10.3,' cm-1') else if(iu.eq.3)then write(6,502)hw la=.true. else write(6,*)' Unknown option' goto 988 endif endif endif else iu=0 endif column=ic if(ic.ge.10)column=ic/10 if(ic.eq.30)column=2 if(ic.eq.40)column=1 if(ic.eq.13)column=3 if(ic.eq.18)column=8 open(9,file=fn) do 1 i=1,3 1 read(9,*) nm=0 2 read(9,*,end=99,err=99)wo,wo,(x,i=1,column) nm=nm+1 if(nm.gt.np0)call report('Too many transitions') if(lnm)wo=1.0d7/wo w(nm)=wo t(nm)=x goto 2 99 close(9) write(6,*)nm, ' modes' do 3 i=1,npoint 3 s(i)=0.0d0 dw=(wsmax-wsmin)/dble(npoint-1) if(LG)then dd0=hw/2.0d0/dsqrt(log(2.0d0)) else dd0=hw/2.0d0 endif dd=dd0 do 4 i=1,npoint xx=wsmin+dw*dble(i-1) do 4 i1=1,nm wo=w(i1) const=t(i1) if(iu.eq.0)then sfc=sf(dd,LG,xx,wo,ic) wm=wo else if(lwe)then c lambda in cm: cm=1.0d7/xx wm=1.0d7/wo sfc=sf(dd,LG,cm,wm,ic) else c lambda in nm: c if la then dd0 in cm: if(la)dd=dd0*xx**2/1.0d7 sfc=sf(dd,LG,xx,wo,ic) wm=wo endif endif c abs: if(ic.eq.1)const=108.7d0*t(i1)*wm c vcd: if(ic.eq.2)const=435.0d0*t(i1)*wm c mcd - B-terms: if(ic.eq.20)const=-5.98442d-3*t(i1)*wm c Raman: if(ic.eq.3)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp))) c roa: if(ic.eq.8)const=t(i1)/wo/(1.0d0-exp(-(1.44d0*wo/temp))) c Raman - absolute: if(ic.eq.13)const=t(i1) c roa - absolute: if(ic.eq.18)const=t(i1) c ord: if(ic.eq.30)const=t(i1)*8.35d18 c nmr: if(ic.eq.40)const=t(i1) 4 s(i)=s(i)+sfc*const c call wrspec2('S.PRN',npoint,wsmin,dw,s) end c ================================================================== subroutine wrspec2(sn,nw,wmin,dw,spec) implicit none character*(*) sn integer*4 nw,iw real*8 wmin,dw,spec(*),w open(8,file=sn) c format according to the interval between the points: if(dabs(dw).gt.1.0d-1)then do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w,spec(iw) 805 format(f12.2,' ',g25.12) else if(dabs(dw).gt.1.0d-3)then do 2 iw=1,nw w=wmin+dble(iw-1)*dw 2 write(8,806)w,spec(iw) 806 format(f14.4,' ',g25.12) else do 3 iw=1,nw w=wmin+dble(iw-1)*dw 3 write(8,807)w,spec(iw) 807 format(g25.12,' ',g25.12) endif endif close(8) write(6,*)sn,' written' return end c ================================================================== function sf(d,g,x,x0,ic) implicit none logical g integer*4 ic real*8 pi,spi,sf,d,x,x0,dd,ab pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(ic.eq.30)then c ORD, real parts, only Lorentz possible: ab=(x-x0)*(x+x0)*x**2*x0**2 sf=ab/(ab**2+(d**2)**2) else c dispersion spectra, imaginary part if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d10)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif endif return end c ================================================================== subroutine report(s) character*(*)s write(6,*)s stop end