program dcroa c analytical CD for a row of twisted cigars implicit none real*8 V,w,pi,fi,u0,k,dist,dk,ac,as,magic,am,ap,wcm,distA,u0D, 1l_ss,l_sc,l_cs,l_cc, 2s_ss,s_sc,s_cs,s_cc,a,bohr,cm,debye,sm,sp,fiD, 3alt(9),ali(9),gi(9),gt(9),ai(27),at(27),al0(9),sf,cf, 4aj,sjk,cjk,sma,spa integer*4 N,np,i,j,iargc,ii real*8,allocatable::r(:),d(:),lambda(:),ks(:),ram(:),roa(:) character*8 s8 pi=4.0d0*atan(1.0d0) bohr=0.529177d0 cm=219474.0d0 debye=2.541746230211d0 magic=pi*bohr*cm*debye**2*1.0d-8 call iopt(wcm,distA,N,u0D,am,ap,fiD) if(iargc().ge.1)then call getarg(1,s8) read(s8,*)N endif if(iargc().ge.2)then call getarg(2,s8) read(s8,*)u0D endif if(iargc().ge.3)then call getarg(3,s8) read(s8,*)fiD endif if(iargc().ge.4)then call getarg(4,s8) read(s8,*)am endif if(iargc().ge.5)then call getarg(5,s8) read(s8,*)ap endif w=wcm/cm dist=distA/bohr u0=u0D/debye fi=fiD*pi/180.0d0 c potential: V=-u0**2*cos(fi)/dist**3 c number of energies: np=N write(6,6000)N,u0D,am,ap,fiD,distA,wcm,V*cm 6000 format(' Number of dipoles ',i12,/, 1 ' Dipole size ',f12.6,' debye',/, 1 ' polar || ',f12.6,' au',/, 1 ' polar _|_ ',f12.6,' au',/, 1 ' Mutual twist ',f12.6,' deg',/, 1 ' Distance ',f12.6,' A',/, 1 ' Frequency ',f12.6,' cm-1',/, 1 ' Potential element ',f12.6,' cm-1',/) c write(6,*) 0.111183288d0,s_ss(3.5d0,4.7d0,17) c write(6,*)-0.270666095d0,s_sc(3.5d0,4.7d0,17) c write(6,*) 0.945204956d0,s_cs(3.5d0,4.7d0,17) c write(6,*) 0.129688748d0,s_cc(3.5d0,4.7d0,17) c write(6,*) 1.619861744d0,l_ss(3.5d0,4.7d0,17) c write(6,*)-1.594762699d0,l_sc(3.5d0,4.7d0,17) c write(6,*) 7.434605983d0,l_cs(3.5d0,4.7d0,17) c write(6,*)10.206782d0 ,l_cc(3.5d0,4.7d0,17) call vz(al0,9) al0(1+3*(1-1))=am al0(2+3*(2-1))=ap al0(3+3*(3-1))=ap if(N.eq.2)then cf=cos(fi) sf=sin(fi) allocate(lambda(2),r(2),d(2),ks(2),ram(2),roa(2)) lambda(1)=w-V d(1)=u0**2*(1.0d0+cf) r(1)=-dist*w*u0**2/2.0d0*sin(fi) ks(1)=1 c <0|a|1>=sum(j) cj aj call vz(alt,9) call vz(gt,9) call vz(at,27) call mkga(0.0d0,0.0d0,0.0d0,al0,gi,ai) call suma(alt,gt,at,al0,gi,ai,1.0d0/dsqrt(2.0d0)) call mkali(ali,am,ap,sf,cf) call mkga(0.0d0,0.0d0,dist,ali,gi,ai) call suma(alt,gt,at,ali,gi,ai,1.0d0/dsqrt(2.0d0)) call doram(ram(1),roa(1),alt,gt,at) c call wr('al-1:',alt,9) c call wr('G -1:', gt,9) c call wr('A -1:', at,27) lambda(2)=w+V d(2)=u0**2*(1.0d0-cf) r(2)=dist*w*u0**2/2.0d0*sin(fi) ks(2)=-1 call vz(alt,9) call vz(gt,9) call vz(at,27) call mkga(0.0d0,0.0d0,0.0d0,al0,gi,ai) call suma(alt,gt,at,al0,gi,ai,1.0d0/dsqrt(2.0d0)) call mkga(0.0d0,0.0d0,dist,ali,gi,ai) call suma(alt,gt,at,ali,gi,ai,-1.0d0/dsqrt(2.0d0)) call doram(ram(2),roa(2),alt,gt,at) c call wr('al-2:',alt,9) c call wr('G -2:', gt,9) c call wr('A -2:', at,27) ii=2 else c transition index: ii=0 allocate(lambda(2*np+1),r(2*np+1),d(2*np+1),ks(2*np+1)) allocate(ram(2*np+1),roa(2*np+1)) c sums of positive and negative rotational strengths: sp=0.0d0 sm=0.0d0 spa=0.0d0 sma=0.0d0 c k-0 ... pi dk=pi/np k=-dk do 1 i=0,np k=k+dk as=0.0d0 ac=0.0d0 aj=0.0d0 do 2 j=1,N aj=aj+1.0d0 as=as+(sin(k*aj))**2 2 ac=ac+(cos(k*aj))**2 if(dabs(as).lt.1.0d-10)as=0.0d0 if(dabs(ac).lt.1.0d-10)ac=0.0d0 c sinus wave: if(as.ne.0.0d0)then ii=ii+1 a=1.0d0/dsqrt(as) lambda(ii)=w+2.0d0*V*cos(k) d(ii) = a**2*u0**2*s_ss(k,fi,N)**2 r(ii) =-a**2*u0**2*w*dist/2.0d0* 1 (s_ss(k,fi,N)*l_sc(k,fi,N)-s_sc(k,fi,N)*l_ss(k,fi,N)) if(r(ii).ge.0.0d0)then sp=sp+r(ii) else sm=sm+r(ii) endif ks(ii)=-k call vz(alt,9) call vz(gt,9) call vz(at,27) aj=0.0d0 do 301 j=1,N aj=aj+1.0d0 cf=cos(fi*aj) sf=sin(fi*aj) sjk=sin(k*aj) call mkali(ali,am,ap,sf,cf) call mkga(0.0d0,0.0d0,dist*aj,ali,gi,ai) 301 call suma(alt,gt,at,ali,gi,ai,sjk) call doram(ram(ii),roa(ii),alt,gt,at) ram(ii)=a**2*ram(ii) roa(ii)=a**2*roa(ii) if(roa(ii).ge.0.0d0)then spa=spa+roa(ii) else sma=sma+roa(ii) endif endif c cosinus wave: if(ac.ne.0.0d0)then ii=ii+1 a=1.0d0/dsqrt(ac) lambda(ii)=w+2.0d0*V*cos(k) d(ii)= a**2*u0**2*s_cc(k,fi,N)**2 r(ii)=-a**2*u0**2*w*dist/2.0d0* 1 (s_cs(k,fi,N)*l_cc(k,fi,N)-s_cc(k,fi,N)*l_cs(k,fi,N)) if(r(ii).ge.0.0d0)then sp=sp+r(ii) else sm=sm+r(ii) endif ks(ii)=k call vz(alt,9) call vz(gt,9) call vz(at,27) aj=0.0d0 do 201 j=1,N aj=aj+1.0d0 cf=cos(fi*aj) sf=sin(fi*aj) cjk=cos(k*aj) call mkali(ali,am,ap,sf,cf) call mkga(0.0d0,0.0d0,dist*aj,ali,gi,ai) 201 call suma(alt,gt,at,ali,gi,ai,cjk) call doram(ram(ii),roa(ii),alt,gt,at) ram(ii)=a**2*ram(ii) roa(ii)=a**2*roa(ii) if(roa(ii).ge.0.0d0)then spa=spa+roa(ii) else sma=sma+roa(ii) endif endif 1 continue write(6,605)' cd',sp,sm,sp+sm write(6,605)'roa',spa,sma,spa+sma 605 format(a3,/,' R+ = ',e20.6,/,' R- = ',e20.6,/,' R = ',e20.6) if(sp.ne.0.0d0)then do 3 i=1,ii 3 if(r(i).gt.0.0d0)r(i)=-r(i)*sm/sp endif if(spa.ne.0.0d0)then do 31 i=1,ii 31 if(roa(i).gt.0.0d0)roa(i)=-roa(i)*sma/spa endif endif do 32 i=1,ii d(i)=d(i)/dble(N) r(i)=r(i)/dble(N) roa(i)=roa(i)/dble(N) 32 ram(i)=ram(i)/dble(N) write(6,*)'Spectra normalized to one dipole!' open(90,file='s.tab') open(91,file='a.tab') write(90,900)N write(91,900)N 900 format(' n = ',i4,' normalized spectra'/,/,60(1h-)) do 4 i=1,ii write(90,600)i,lambda(i)*219474.0d0,d(i)*debye**2,r(i)*magic,ks(i) 4 write(91,601)i,lambda(i)*219474.0d0,ram(i),roa(i),roa(i),ks(i) 600 format(i4,f12.6,2g14.6,f8.3) 601 format(i4,f10.4,' 1.0 2.0 ',g14.6,' 4.0 5.0 ',g14.6,' 7.0 ', 1g14.6,f8.3) write(90,901) write(91,901) 901 format(60(1h-)) close(90) close(91) write(6,*)' s.tab a.tab written' end function s_ss(k,f,N) c sum(j=1...N,sin(kj)*sin(fj)) implicit none real*8 s_ss,k,f,m,p,a,m2,p2,ma,pa,s,aj integer*4 N,j s=0.0d0 aj=0.0d0 do 1 j=1,N aj=aj+1.0d0 1 s=s+sin(k*aj)*sin(f*aj) s_ss=s return end function s_sc(k,f,N) c sum(j=1...N,sin(kj)*cos(fj)) implicit none real*8 s_sc,k,f,m,p,a,m2,p2,ma,pa,s,aj integer*4 N,j s=0.0d0 aj=0.0d0 do 1 j=1,N aj=aj+1.0d0 1 s=s+sin(k*aj)*cos(f*aj) s_sc=s return end function s_cs(k,f,N) c sum(j=1...N,cos(kj)*sin(fj)) implicit none real*8 s_cs,k,f,s_sc,sj,s,aj integer*4 N,j s=0.0d0 aj=0.0d0 do 1 j=1,N aj=aj+1.0d0 1 s=s+cos(k*aj)*sin(f*aj) s_cs=s return end function s_cc(k,f,N) c sum(j=1...N,cos(kj)*cos(fj)) implicit none real*8 s_cc,k,f,m,p,a,m2,p2,ma,pa,s,aj integer*4 N,j s=0.0d0 aj=0.0d0 do 1 j=1,N aj=aj+1.0d0 1 s=s+cos(k*aj)*cos(f*aj) s_cc=s return end function l_sc(k,f,N) c sum(j=1...N,j*sin(kj)*cos(fj)) implicit none real*8 l_sc,k,f,m,p,a,m2,p2,ma,pa,h,s,aj integer*4 N,j s=0.0d0 aj=0 do 1 j=1,N aj=aj+1.0d0 1 s=s+aj*sin(k*aj)*cos(f*aj) l_sc=s return end function l_ss(k,f,N) c sum(j=1...N,j*sin(kj)*sin(fj)) implicit none real*8 l_ss,k,f,m,p,a,m2,p2,ma,pa,h,s,aj integer*4 N,j s=0.0d0 aj=0 do 1 j=1,N aj=aj+1.0d0 1 s=s+aj*sin(k*aj)*sin(f*aj) l_ss=s return end function l_cs(k,f,N) c sum(j=1...N,j*cos(kj)*sin(fj)) implicit none real*8 l_cs,k,f,m,p,a,m2,p2,ma,pa,h,s,aj integer*4 N,j s=0.0d0 aj=0 do 1 j=1,N aj=aj+1.0d0 1 s=s+aj*cos(k*aj)*sin(f*aj) l_cs=s return end function l_cc(k,f,N) c sum(j=1...N,j*cos(kj)*cos(fj)) implicit none real*8 l_cc,k,f,m,p,a,m2,p2,ma,pa,h,s,aj integer*4 N,j s=0.0d0 aj=0 do 1 j=1,N aj=aj+1.0d0 1 s=s+aj*cos(k*aj)*cos(f*aj) l_cc=s return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine vz(v,n) real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end subroutine mkga(x,y,z,al,g,a) c makes tensors G and A fromalpha implicit none real*8 x,y,z,al(*),g(*),a(*),s integer*4 i do 1 i=1,3 g(i+3*(1-1))=-0.5d0*(y*al(i+3*(3-1))-z*al(i+3*(2-1))) g(i+3*(2-1))=-0.5d0*(z*al(i+3*(1-1))-x*al(i+3*(3-1))) 1 g(i+3*(3-1))=-0.5d0*(x*al(i+3*(2-1))-y*al(i+3*(1-1))) do 2 i=1,3 s=-x*al(1+3*(i-1))-y*al(2+3*(i-1))-z*al(3+3*(i-1)) a(i+3*(1-1)+9*(1-1))=1.5d0*(x*al(i+3*(1-1))+x*al(i+3*(1-1))) a(i+3*(1-1)+9*(2-1))=1.5d0*(x*al(i+3*(2-1))+y*al(i+3*(1-1))) a(i+3*(1-1)+9*(3-1))=1.5d0*(x*al(i+3*(3-1))+z*al(i+3*(1-1))) a(i+3*(2-1)+9*(1-1))=1.5d0*(y*al(i+3*(1-1))+x*al(i+3*(2-1))) a(i+3*(2-1)+9*(2-1))=1.5d0*(y*al(i+3*(2-1))+y*al(i+3*(2-1))) a(i+3*(2-1)+9*(3-1))=1.5d0*(y*al(i+3*(3-1))+z*al(i+3*(2-1))) a(i+3*(3-1)+9*(1-1))=1.5d0*(z*al(i+3*(1-1))+x*al(i+3*(3-1))) a(i+3*(3-1)+9*(2-1))=1.5d0*(z*al(i+3*(2-1))+y*al(i+3*(3-1))) a(i+3*(3-1)+9*(3-1))=1.5d0*(z*al(i+3*(3-1))+z*al(i+3*(3-1))) a(i+3*(1-1)+9*(1-1))=a(i+3*(1-1)+9*(1-1))-s a(i+3*(2-1)+9*(2-1))=a(i+3*(2-1)+9*(2-1))-s 2 a(i+3*(3-1)+9*(3-1))=a(i+3*(3-1)+9*(3-1))-s return end subroutine suma(alt,gt,at,ali,gi,ai,c) c add i-polarizabilities to the totals implicit none real*8 alt(*),gt(*),at(*),ali(*),gi(*),ai(*),c integer*4 i do 101 i=1,9 alt(i)=alt(i)+ali(i)*c gt(i) =gt(i) +gi(i)*c at(i ) =at(i ) +ai(i )*c at(i+ 9) =at(i+ 9) +ai(i+ 9)*c 101 at(i+18) =at(i+18) +ai(i+18)*c return end FUNCTION EPS(I,J,K) IMPLICIT INTEGER*4 (I-N) REAL*8 EPS EPS=0.0D0 IF(I.EQ.1)then if(J.EQ.2.AND.K.EQ.3)EPS= 1.0D0 IF(J.EQ.3.AND.K.EQ.2)EPS=-1.0D0 endif IF(I.EQ.2)then IF(J.EQ.3.AND.K.EQ.1)EPS= 1.0D0 IF(J.EQ.1.AND.K.EQ.3)EPS=-1.0D0 endif IF(I.EQ.3)then IF(J.EQ.1.AND.K.EQ.2)EPS= 1.0D0 IF(J.EQ.2.AND.K.EQ.1)EPS=-1.0D0 endif RETURN END subroutine doram(ram,roa,al,g,a) implicit none real*8 ram,roa,al(*),g(*),a(*),SAL0,SAL1,SAG0,SAG1,SA1,SPAL, 1EPS,AMU,BOHR,SpA3,beta2,gpisvejc,EXCA,betaa,betag integer*4 i,j,k,l EXCA=5320.0d0 AMU=1822.0d0 BOHR=0.529177d0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPAL=0.0d0 DO 3 I=1,3 SPAL=SPAL+al(I+3*(I-1)) DO 3 J=1,3 SAL0=SAL0+al(I+3*(I-1))*al(J+3*(J-1)) SAL1=SAL1+al(I+3*(J-1))*al(I+3*(J-1)) SAG0=SAG0+al(I+3*(I-1))* g(J+3*(J-1)) SAG1=SAG1+al(I+3*(J-1))* g(I+3*(J-1)) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*al(I+3*(J-1))*a(K+3*(L-1)+9*(J-1))/3.0d0 SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) betaa=SA1*3.0d0/2.0d0*gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc roa=96.0d0*(betag+betaa/3.0d0) ram=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) return end subroutine mkali(ali,am,ap,sf,cf) implicit none real*8 ali(*),am,ap,sf,cf call vz(ali,9) ali(1+3*(1-1))= am*cf**2+ap*sf**2 ali(1+3*(2-1))= (am-ap)*sf*cf ali(2+3*(1-1))= (am-ap)*sf*cf ali(2+3*(2-1))= ap*cf**2+am*sf**2 ali(3+3*(3-1))= ap return end subroutine iopt(w,dist,N,u0,am,ap,fi) implicit none integer*4 N character*4 k real*8 w,dist,u0,am,ap,fi logical lex N=2 u0=0.1d0 dist=3.0d0 am=0.1d0 ap=0.2d0 fi=20.0d0 w=1650.0d0 inquire(file='DCROA.OPT',exist=lex) if(lex)then open(9,file='DCROA.OPT') 1 read(9,900,end=99,err=99)k 900 format(a4) if(k.eq.'WEXC')read(9,*)w if(k.eq.'DIST')read(9,*)dist if(k.eq.'NDIP')read(9,*)N if(k.eq.'PEXC')read(9,*)u0 if(k.eq.'AMEX')read(9,*)am if(k.eq.'APEX')read(9,*)ap if(k.eq.'ANGL')read(9,*)fi goto 1 99 close(9) endif return end subroutine wr(s,a,n) implicit none character*(*) s real*8 a(*) integer*4 n,i write(6,'(a)')s write(6,600)(a(i),i=1,n) 600 format(3F12.6) return end