program dcanal c analytical CD for a row of twisted dipoles implicit none real*8 V,w,pi,fi,u0,k,dist,dk,ac,as,magic, 1l_ss,l_sc,l_cs,l_cc, 2s_ss,s_sc,s_cs,s_cc,a,bohr,cm,debye,sm,sp integer*4 N,np,i,j,iargc,ii real*8,allocatable::r(:),d(:),lambda(:),ks(:) 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 c frequency: w=1650.0d0/cm c distance 3 A dist=3.0d0/bohr if(iargc().ne.3)then write(6,*)' usage dcanal N u0 fi' stop endif c read number of dipoles: call getarg(1,s8) read(s8,*)N c read transition dipole in debye call getarg(2,s8) read(s8,*)u0 c read mutual twist: call getarg(3,s8) read(s8,*)fi c potential: V=-(u0/debye)**2*cos(fi*pi/180.0d0)/dist**3 c number of energies: np=N write(6,6000)N,u0,fi,dist*bohr,w*cm,V*cm 6000 format(' Number of dipoles ',i12,/, 1 ' Dipole size ',f12.6,' debye',/, 1 ' Mutual twist ',f12.6,' deg',/, 1 ' Distance ',f12.6,' A',/, 1 ' Frequency ',f12.6,' cm-1',/, 1 ' Potential element ',f12.6,' cm-1',/) fi=fi*pi/180.0d0 u0=u0/debye 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) if(N.eq.2)then allocate(lambda(2),r(2),d(2),ks(2)) lambda(1)=w-V d(1)=u0**2*(1.0d0+cos(fi)) r(1)=-dist*w*u0**2/2.0d0*sin(fi) ks(1)=1 lambda(2)=w+V d(2)=u0**2*(1.0d0-cos(fi)) r(2)= dist*w*u0**2/2.0d0*sin(fi) ks(2)=-1 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)) c sums of positive and negative rotational strengths: sp=0.0d0 sm=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 do 2 j=1,N as=as+(sin(k*j))**2 2 ac=ac+(cos(k*j))**2 c sinus wave: if(i.ne.0)then a=1.0d0/dsqrt(as) ii=ii+1 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 endif c cosinus wave: a=1.0d0/dsqrt(ac) ii=ii+1 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 1 ks(ii)=k write(6,605)sp,sm,sp+sm 605 format(' 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 endif open(90,file='s.tab') write(90,900)N 900 format(' n = ',i4,/,/,60(1h-)) do 4 i=1,ii 4 write(90,600)i,lambda(i)*219474.0d0,d(i)*debye**2,r(i)*magic,ks(i) 600 format(i4,f12.6,2g14.6,f8.3) write(90,901) 901 format(60(1h-)) close(90) write(6,*)' s.tab written' end function s_ss(k,f,N) c sum(j=1...N,sin(kj)*sin(kj)) implicit none real*8 s_ss,k,f,m,p,a,m2,p2,ma,pa integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 a=dble(N)+0.5d0 ma=m*a pa=p*a s_ss=0.25d0*(sin(ma)/sin(m2)-sin(pa)/sin(p2)) return end function s_sc(k,f,N) c sum(j=1...N,sin(kj)*cos(kj)) implicit none real*8 s_sc,k,f,m,p,a,m2,p2,ma,pa integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 a=dble(N)+0.5d0 ma=m*a pa=p*a s_sc=0.25d0*((cos(ma)-cos(m2))/sin(m2)-(cos(pa)-cos(p2))/sin(p2)) return end function s_cs(k,f,N) c sum(j=1...N,sin(kj)*cos(kj)) implicit none real*8 s_cs,k,f,s_sc integer*4 N s_cs=s_sc(f,k,N) return end function s_cc(k,f,N) c sum(j=1...N,sin(kj)*cos(kj)) implicit none real*8 s_cc,k,f,m,p,a,m2,p2,ma,pa integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 a=dble(N)+0.5d0 ma=m*a pa=p*a s_cc=0.25d0*((sin(ma)-sin(m2))/sin(m2)+(sin(pa)-sin(p2))/sin(p2)) return end function l_sc(k,f,N) c sum(j=1...N,j*sin(kj)*cos(kj)) implicit none real*8 l_sc,k,f,m,p,a,m2,p2,ma,pa,h integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 h=0.5d0 a=dble(N)+h ma=m*a pa=p*a l_sc=-0.25d0*((-a*cos(ma)+h*cos(m2))/sin(m2) 1 +( a*cos(pa)-h*cos(p2))/sin(p2) 2+(sin(ma)-sin(m2))*cos(m2)*h/(sin(m2))**2 2-(sin(pa)-sin(p2))*cos(p2)*h/(sin(p2))**2) return end function l_ss(k,f,N) c sum(j=1...N,j*sin(kj)*sin(kj)) implicit none real*8 l_ss,k,f,m,p,a,m2,p2,ma,pa,h integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 h=0.5d0 a=dble(N)+h ma=m*a pa=p*a l_ss=-0.25d0*(-a*sin(ma)/sin(m2) 1 +a*sin(pa)/sin(p2) 2+(cos(m2)-cos(ma))*cos(m2)*h/(sin(m2))**2 2-(cos(p2)-cos(pa))*cos(p2)*h/(sin(p2))**2) return end function l_cs(k,f,N) c sum(j=1...N,j*sin(kj)*sin(kj)) implicit none real*8 l_cs,k,f,m,p,a,m2,p2,ma,pa,h integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 h=0.5d0 a=dble(N)+h ma=m*a pa=p*a l_cs= 0.25d0*(-a*cos(ma)/sin(m2) 1 -a*cos(pa)/sin(p2) 2+cos(m2)*sin(ma)*h/(sin(m2))**2 2+cos(p2)*sin(pa)*h/(sin(p2))**2) return end function l_cc(k,f,N) c sum(j=1...N,j*sin(kj)*cos(kj)) implicit none real*8 l_cc,k,f,m,p,a,m2,p2,ma,pa,h integer*4 N m=f-k p=f+k m2=m/2.0d0 p2=p/2.0d0 h=0.5d0 a=dble(N)+h ma=m*a pa=p*a l_cc=0.25d0*((a*sin(ma)-h*sin(m2))/sin(m2) 1 +(a*sin(pa)-h*sin(p2))/sin(p2) 2+(cos(ma)-cos(m2))*cos(m2)*h/(sin(m2))**2 2+(cos(pa)-cos(p2))*cos(p2)*h/(sin(p2))**2) return end