program xqq c crystal spectra implicit none integer*4 i,nat,iccub,n,ice,ix,nat1,np,is,n3,ni,nj,nk, 1n27,ierr,qi,qj,qk,iq,it,IWR real*8,allocatable::r(:),P(:,:,:),A(:,:,:),f(:,:),Q1(:,:,:), 1P1(:,:),A1(:,:),al(:,:,:),gp(:,:,:),aa(:,:,:,:), 1al1(:,:,:),gp1(:,:,:),aa1(:,:,:,:),Dms(:,:),Dmc(:,:),Sqr(:,:), 1Sqi(:,:),Eq(:),fv1(:),fv2(:),fm1(:,:),m(:),sd(:,:),mu(:), 1qq(:,:) integer*4,allocatable::z(:) real*8 vi(3),vj(3),vk(3),sa,ecm,c1,wsmax,wsmin,DSX(3),QSX(3), 1MSX(3),ram(13),ENM,temp,ds,rs,pi,po,fwhm,ri,rj,rk,delq integer*4,allocatable::ip(:,:) logical ltab,lttt,lqpt,ldip,lg,lcub,lone character*8 QF(8),SF(16) data QF/'DOGX.TAB','DOGY.TAB','DOGZ.TAB',' DOG.TAB', 1 'ROAX.TAB','ROAY.TAB','ROAZ.TAB',' ROA.TAB'/ data SF/' d.prn',' dx.prn',' dy.prn',' dz.prn', 1 ' r.prn',' rx.prn',' ry.prn',' rz.prn', 1 ' ram.prn','ramx.prn','ramy.prn','ramz.prn', 1 ' roa.prn','roax.prn','roay.prn','roaz.prn'/ c1=1302.828d0 pi=4.0d0*atan(1.0d0) po=dsqrt(2.0d0*pi) write(6,601) 601 format(/, 1' Propagate periodic crystal-like 1D structures',/, 1' Input: FILE.X = CRYSTAL.X',/, 1 ' FILE.FC',/, 1 ' FILE.TEN',/, 1 ' FILE.QEN',/, 1 ' FILE.TTT',/, 1 ' XQQ.OPT',/, 1 ' CELL.TXT',/, 1 ' FTRY.INP',/,/, 1' Output: DOG.TAB DOGX.TAB DOGY.TAB DOGZ.TAB',/, 1' - isotropic and oriented IR/VCD intenisites',/, 1' ROA.TAB ROAX.TAB ROAY.TAB ROAZ.TAB',/, 1' - isotropic and oriented RAM/ROA intenisites',/, 1' d.prn dx.prn dy.prn dz,prn',/, 1' - isotropic and oriented IR spectra',/, 1' r.prn rx.prn ry.prn rz,prn',/, 1' - isotropic and oriented VCD spectra',/, 1' ram.prn ramx.prn ramy.prn ramz,prn',/, 1' - isotropic and oriented Raman spectra',/, 1' roa.prn roax.prn roay.prn roaz,prn',/, 1' - isotropic and oriented ROA spectra',/, 1' BIGC.X geometry',/, 1' DIP.TXT transition dipoles',/, 1' Q.PRNi mode contributions',/) call readopt(NI,NJ,NK,ltab,lttt,lqpt,ldip,np,n27,wsmin,wsmax, 1temp,fwhm,lg,lcub,lone,delq,ENM,IWR) allocate(z(1),r(3*1)) CALL READX(0,nat,z,r) deallocate(z,r) nat1=nat/n27 if(mod(nat,n27).ne.0)call report('mod nat n27 <> 0') write(6,6002)n27,nat1 6002 format(i6,' cells',i6,' atoms in a cell') n=3*nat n3=3*nat1 allocate(r(n),f(n,n),P(nat,3,3),A(nat,3,3),Eq(n3), 1al(n,3,3),gp(n,3,3),aa(n,3,3,3),Q1(n3,3,3), 1P1(n3,3),A1(n3,3),Dms(n3,n3),Dmc(n3,n3),fv1(n3),fv2(n3), 1al1(n3,3,3),gp1(n3,3,3),aa1(n3,3,3,3),ip(n27,3),fm1(2,n3), 1Sqr(n3,n3),Sqi(n3,n3),m(nat),mu(nat1),qq(n,6),z(nat)) ip(1,1)=0 CALL READX(1,nat,z,R) open(9,file='CELL.TXT',status='old') read(9,*)vi read(9,*)vj read(9,*)vk close(9) c determine central cube: ice=iccub(nat,r,n27,ip,vi,vj,vk) c read masses, force filed, dipole, quadrupole and pol. derivatives: CALL READM(nat,m) CALL READFF(n,f) CALL MASSFF(n,f,m) CALL READTEN(nat,P,A) call reqd(lqpt,qq,nat,R,P) if(lttt)CALL READTTT(al,aa,gp,NAT) c rewrite central cube to special variables: call massa(nat1,m,mu) call rwr(n,ice,nat1,nat,P1,A1,Q1,P,A,qq,al1,gp1,aa1, 1al,gp,aa) if(ldip)OPEN(29,FILE='DIP.TXT') c initialize tables: if(ltab)then OPEN(18,FILE='DOG.TAB') WRITE (18,501) 501 FORMAT (/,' FREQ. ', 1 ' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2 /,'-----------------------------------------------------------') do 59 i=1,3 open(18+i,file=QF(i)) 59 write(18+i,501) if(lttt)then do 62 i=5,8 call setis(QF(I),is) open(18+i,file=QF(i)(is:len(QF(i)))) 62 write(18+i,501) endif endif c initialize spectra: allocate(sd(16,np)) call vz2(sd,16,np) c spectra: 1 absorbtion c 2 abs x c 3 abs y c 4 abs z c 5 vcd c 6 vcd x c 7 vcd y c 8 vcd z c 9 ram c 10 ramx c 11 ramy c 12 ramz c 13 roa c 14 roax c 15 roay c 16 roaz c loop over wave vectors: it=0 c qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq ri=-delq do 1000 qi=0,NI-1 ri=ri+delq rj=-delq do 1000 qj=0,NJ-1 rj=rj+delq rk=-delq do 1000 qk=0,NK-1 rk=rk+delq c real write(6,6001)qi,qj,qk,ri,rj,rk 6001 format(3i5,3f8.3) c make dynamic matrix: call maked(ice,n,n3,f,Dmc,Dms,ri,rj,rk,NI,NJ,NK,n27,ip,lcub) c diagonalize: call ch(n3,Dmc,Dms,Eq,1,Sqr,Sqi,fv1,fv2,fm1,ierr) if(ierr.ne.0)call report('diagonalization not successfull') c now sum(Sqr(k,i)*Sqr(k,j)+Sqi(k,i)*Sqi(k,j)) = delta(i,j) c make mass-weighted S-matrix: call sm(Sqr,Sqi,mu,n3) do 11 iq=1,n3 sa=1.0d0 if(Eq(iq).lt.0.0d0)sa=-sa 11 Eq(iq)=sa*sqrt(dabs(Eq(iq)))*c1 if(IWR.eq.1)write(6,602)Eq 602 format(6f12.2) c loope over normal modes of this wave vector: do 1000 iq=1,n3 ecm=Eq(iq) if(ecm.gt.1.0d0)then it=it+1 call getdip(iq,ri,rj,rk,n3,P1,A1,Q1,Sqr,Sqi,NI,NJ,NK,ecm, 1 ds,rs,DSX,MSX,QSX,ldip,lttt,al1,gp1,aa1,ENM,ram,vi,vj,vk,lone) if(ltab)then WRITE(18,523)it,ecm,ds,rs,rs,po*Eq(iq)*ds do 60 ix=1,3 60 WRITE(18+ix,523)it,ecm,DSX(ix),QSX(ix)+MSX(ix) 523 FORMAT(I5,f16.4,4G12.4) if(lttt)then do 67 ix=5,7 67 WRITE(18+ix,3000)it,ecm,ram(1),ram(2),ram(3+ix),ram(4),ram(4), 1 ram(6+ix),ram(4),ram(6+ix),ram(6),ram(7) WRITE(18+8,3000)it,ecm,ram(1),ram(2),ram(3),ram(4),ram(4), 1 ram(5),ram(4),ram(5),ram(6),ram(7) 3000 FORMAT(I5,f9.2,6g12.3,f6.3,4g12.4) endif endif call sspread(ecm,ds,DSX(1),DSX(2),DSX(3),rs,QSX(1)+MSX(1), 1 QSX(2)+MSX(2),QSX(3)+MSX(3),sd,np,wsmin,wsmax,fwhm,lg) if(lttt)call spreadr(ecm,ram,sd,np,wsmin,wsmax,fwhm,lg,temp) endif 1000 continue c qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq if(ltab)then do 61 ix=0,3 WRITE(18+ix,1819) 1819 FORMAT('-------------------------------------------------') 61 close(18+ix) write(6,*)'DOG.TAB DOGX.TAB DOGY.TAB DOGZ.TAB written' if(lttt)then do 63 ix=5,8 WRITE(18+ix,1819) 63 close(18+ix) write(6,*)'ROA.TAB ROAX.TAB ROAY.TAB ROAZ.TAB written' endif endif do 64 ix=1,8 call setis(SF(ix),is) 64 call wspectra(SF(ix)(is:len(SF(ix))),wsmin,wsmax,sd,np,ix) if(lttt)then do 65 ix=9,16 call setis(SF(ix),is) 65 call wspectra(SF(ix)(is:len(SF(ix))),wsmin,wsmax,sd,np,ix) endif if(ldip)close(29) end subroutine trpa(iq,nu3,p,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3),ur(3), 1p(nu3,3) integer*4 nu3,iq,k,l do 1 k=1,3 ur(k)=0.0d0 ui(k)=0.0d0 do 1 l=1,nu3 ur(k)=ur(k)+p(l,k)*sr(l,iq) 1 ui(k)=ui(k)+p(l,k)*si(l,iq) return end SUBROUTINE READX(ic,nat,z,r) IMPLICIT none real*8 r(*) integer*4 ic,nat,i,ix,z(*) open(2,file='FILE.X',status='old') READ(2,*) read(2,*)nat if(ic.eq.1)then do 1 i=1,nat 1 READ(2,*)z(i),(r(ix+3*(i-1)),ix=1,3) endif close(2) RETURN END subroutine trpq(iq,nu3,q,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3,3),ur(3,3), 1q(nu3,3,3) integer*4 nu3,iq,k,l,m do 1 k=1,3 do 1 m=1,3 ur(k,m)=0.0d0 ui(k,m)=0.0d0 do 1 l=1,nu3 ur(k,m)=ur(k,m)+q(l,k,m)*sr(l,iq) 1 ui(k,m)=ui(k,m)+q(l,k,m)*si(l,iq) return end subroutine traa(iq,nu3,p,sr,si,ur,ui) implicit none real*8 sr(nu3,nu3),si(nu3,nu3),ui(3,3,3),ur(3,3,3), 1p(nu3,3,3,3) integer*4 nu3,iq,k,l,m,n do 1 k=1,3 do 1 m=1,3 do 1 n=1,3 ur(k,m,n)=0.0d0 ui(k,m,n)=0.0d0 do 1 l=1,nu3 ur(k,m,n)=ur(k,m,n)+p(l,k,m,n)*sr(l,iq) 1 ui(k,m,n)=ui(k,m,n)+p(l,k,m,n)*si(l,iq) return end subroutine getdip(iqt,qi,qj,qk,nu3,P1,A1,Q1,Sqr,Sqi,NI,NJ,NK,ecm, 1ds,rs,DSX,MSX,QSX,ldip,lttt,al1,gp1,aa1,ENM,ram,vi,vj,vk,lone) implicit none integer*4 nu3,iq,NI,NJ,NK,i,j,k,iqt,l,d real*8 P1(nu3,3),A1(nu3,3),sIr,sIi,sJr,sJi,sKr,sKi,Q1(nu3,3,3), 1Sqr(nu3,nu3),Sqi(nu3,nu3),ur(3),sr,si,ui(3),mr(3),mi(3),qr(3,3), 1qdi(3,3),pi,ANI,ANJ,ANK,urt(3),uit(3),mrt(3),mit(3),qrt(3,3), 1qit(3,3),ecm,debye,CM,AMU,clight,FD,FC,BOHR,WR,aulambda,auk,ds,rs, 1DSX(*),MSX(*),QSX(*),sp,fi,ENM,ram(*),vi(3),vj(3),vk(3),qi,qj,qk, 1al1(nu3,3,3),gp1(nu3,3,3),fj,fk,tIr,tIi,tJr,tJi,tKr,tKi, 1aa1(nu3,3,3,3),ali(3,3),alr(3,3),gi(3,3),gr(3,3),aai(3,3,3), 1aar(3,3,3),alti(3,3),altr(3,3),gti(3,3),gtr(3,3),aati(3,3,3), 1aatr(3,3,3),rr(3),ri(3),s2ii,s2ir,s2ji,s2jr,s2ki,s2kr logical ldip,lttt,lone debye=2.541732d0 CM=219474.63d0 AMU=1822.887D0 clight=137.03604d0 BOHR=0.52917705993d0 iq=iqt ANI=dble(NI) ANJ=dble(NJ) ANK=dble(NK) pi=4.0d0*datan(1.0d0) fi=2.0d0*pi*qi/ANI fj=2.0d0*pi*qj/ANJ fk=2.0d0*pi*qk/ANK if(lone)then c one cell only sIr=1.0d0 sIi=0.0d0 sJr=1.0d0 sJi=0.0d0 sKr=1.0d0 sKi=0.0d0 s2Ir=0.0d0 s2Ii=0.0d0 s2Jr=0.0d0 s2Ji=0.0d0 s2Kr=0.0d0 s2Ki=0.0d0 else c sum of phase factors over 0...NI-1 cells, etc.: call s1s(fi,sIr,sIi,ANI) call s1s(fj,sJr,sJi,ANJ) call s1s(fk,sKr,sKi,ANK) call s2s(fi,s2Ir,s2Ii,ANI) call s2s(fj,s2Jr,s2Ji,ANJ) call s2s(fk,s2Kr,s2Ki,ANK) endif call m3(sIr,sIi,sJr,sJi,sKr,sKi,sr,si) call m3(s2Ir,s2Ii,sJr,sJi,sKr,sKi,tIr,tIi) call m3(sIr,sIi,s2Jr,s2Ji,sKr,sKi,tJr,tJi) call m3(sIr,sIi,sJr,sJi,s2Kr,s2Ki,tKr,tKi) do 6 i=1,3 rr(i)=(tIr*vi(i)+tJr*vj(i)+tKr*vk(i))/BOHR 6 ri(i)=(tIi*vi(i)+tJi*vj(i)+tKi*vk(i))/BOHR FD =debye*dsqrt(CM/AMU/2.0d0) FC=debye*dsqrt(2.0d0/CM/AMU)/clight WR=dsqrt(ecm) c wavelength in atomic units aulambda=1.0d0/ecm*1.0d8/BOHR c wave vector in atomic units auk=2.0d0*pi/aulambda c transform complex dipolar derivatives to normal mode coordinates c du/dQ = dU/dx . S call trpa(iq,nu3,P1,Sqr,Sqi,ur,ui) call trpa(iq,nu3,A1,Sqr,Sqi,mr,mi) call trpq(iq,nu3,Q1,Sqr,Sqi,qr,qdi) if(lttt)then call trpq(iq,nu3,al1,Sqr,Sqi,alr,ali) call trpq(iq,nu3,gp1,Sqr,Sqi, gr,gi ) call traa(iq,nu3,aa1,Sqr,Sqi,aar,aai) endif c c intrinsic, dipoles: do 1 k=1,3 urt(k)=(ur(k)*sr-ui(k)*si)/WR*FD uit(k)=(ui(k)*sr+ur(k)*si)/WR*FD mrt(k)=(mr(k)*sr-mi(k)*si)*WR*FC mit(k)=(mi(k)*sr+mr(k)*si)*WR*FC do 1 j=1,3 qrt(k,j)=(qr(k,j)*sr-qdi(k,j)*si)/WR*FD*0.5d0 1 qit(k,j)=(qdi(k,j)*sr+qr(k,j)*si)/WR*FD*0.5d0 c intrinsic, polarizabilities: if(lttt)then do 11 k=1,3 do 11 j=1,3 altr(k,j)=(alr(k,j)*sr-ali(k,j)*si) alti(k,j)=(ali(k,j)*sr+alr(k,j)*si) gtr( k,j)=(gr( k,j)*sr-gi( k,j)*si) gti( k,j)=(gi( k,j)*sr+gr( k,j)*si) do 11 l=1,3 aatr(k,j,l)=(aar(k,j,l)*sr-aai(k,j,l)*si) 11 aati(k,j,l)=(aai(k,j,l)*sr+aar(k,j,l)*si) endif c non-local parts: do 4 i=1,3 call setjk(i,j,k) mrt(i)=mrt(i)+(rr(j)*ur(k)-ri(j)*ui(k) 1 -rr(k)*ur(j)+ri(k)*ui(j))*WR*FC/4.0d0 4 mit(i)=mit(i)+(ri(j)*ur(k)+rr(j)*ui(k) 1 -ri(k)*ur(j)-rr(k)*ui(j))*WR*FC/4.0d0 do 5 i=1,3 do 5 j=1,3 qrt(i,j)=qrt(i,j)+(rr(i)*ur(j)-ri(i)*ui(j) 1 + ur(i)*rr(j)-ui(i)*ri(j))/WR*FD/2.0d0 5 qit(i,j)=qit(i,j)+(rr(i)*ui(j)+ri(i)*ur(j) 1 + ui(i)*rr(j)+ur(i)*ri(j))/WR*FD/2.0d0 if(lttt)then do 2 d=1,3 do 2 i=1,3 call setjk(i,j,k) gtr(d,i)=gtr(d,i) -( rr(j)*alr(d,k)-ri(j)*ali(d,k) 1 - rr(k)*alr(d,j)+ri(k)*ali(d,j))/2.0d0 2 gti(d,i)=gti(d,i) -( ri(j)*alr(d,k)+rr(j)*ali(d,k) 1 - ri(k)*alr(d,j)-rr(k)*ali(d,j))/2.0d0 do 3 d=1,3 do 3 i=1,3 do 3 j=1,3 if(i.eq.j)then aatr(d,i,j)=aatr(d,i,j) 1 -alr(d,1)*rr(1)-alr(d,2)*rr(2)-alr(d,3)*rr(3) 1 +ali(d,1)*ri(1)+ali(d,2)*ri(2)+ali(d,3)*ri(3) aati(d,i,j)=aati(d,i,j) 1 -alr(d,1)*ri(1)-alr(d,2)*ri(2)-alr(d,3)*ri(3) 1 -alr(d,1)*ri(1)-alr(d,2)*ri(2)-alr(d,3)*ri(3) endif aatr(d,i,j)=aatr(d,i,j)+1.5d0*(alr(d,i)*rr(j)+rr(j)*alr(d,i) 1 -ali(d,i)*ri(j)-ri(j)*ali(d,i)) 3 aati(d,i,j)=aati(d,i,j)+1.5d0*(alr(d,i)*ri(j)+rr(j)*ali(d,i) 1 +ali(d,i)*rr(j)+ri(j)*alr(d,i)) endif if(ldip)write(29,619)iq,ecm, 1(urt(i),i=1,3),(uit(i),i=1,3), 1(mrt(i),i=1,3),(mit(i),i=1,3), 1((qrt(i,j)/10.0d0,i=1,j),j=1,3), 1((qit(i,j)/10.0d0,i=1,j),j=1,3) 619 format(i3,f8.1,' u:',6f10.4,/,11x,' m:',6f10.4,/, 111x,' q:',6f10.4,/,11x,' q:',6f10.4) c q: xx xy yy xz yz zz ds=sp(urt,urt)+sp(uit,uit) rs=sp(urt,mrt)+sp(uit,mit) c for i, j and k beam directions: do 58 i=1,3 j=i+1 if(j.eq.4)j=1 k=j+1 if(k.eq.4)k=1 DSX(i)= 1.5d0*(urt(j)*urt(j) + urt(k)*urt(k) 1 +uit(j)*uit(j) + uit(k)*uit(k) ) QSX(i)=-auk*1.5d0*(urt(j)*qrt(k,i)-urt(k)*qrt(i,j) 1 +uit(j)*qit(k,i)-uit(k)*qit(i,j)) 58 MSX(i)= 1.5d0*(urt(j)*mrt(j) + urt(k)*mrt(k) 1 +uit(j)*mit(j) + uit(k)*mit(k) ) if(lttt)call groa(ram,altr,alti,gtr,gti,aatr,aati,ENM,0) return end subroutine s1s(k,r,i,N) c 1 ikI c (r,i)= ------ sum(I=0..N-1) e c sqrt(N) implicit none real*8 k,r,i,N,down,l l=1.0d0 if(dabs(k).gt.1.0d-5)then down=2.0d0*dsqrt(N)*(l-cos(k)) r=(l+cos(k*(N-l))-cos(k)-cos(k*N))/down i=( sin(k*(N-l))+sin(k)-sin(k*N))/down else r=dsqrt(N) i=0.0d0 endif return end subroutine s2s(k,r,i,N) c 1 ikI c (r,i)= ------ sum(I=0..N-1) I e c sqrt(N) implicit none real*8 k,r,i,N,down,l l=1.0d0 if(dabs(k).gt.1.0d-5)then down=4.0d0*dsqrt(N)*(l-cos(k))**2 r=( (N-l)*cos(k*(N+l))-(3.0d0*N-2.0d0)*cos(k*N) 1 +(3.0d0*N-l)*cos(k*(N-l))-N*cos(k*(N-2.0d0)) 1 +2.0d0*cos(k)-2.0d0)/down i=( (N-l)*sin(k*(N+l))-(3.0d0*N-2.0d0)*sin(k*N) 1 +(3.0d0*N-l)*sin(k*(N-l))-N*sin(k*(N-2.0d0)))/down else r=dsqrt(N)*(N-l)/2.0d0 i=0.0d0 endif return end subroutine setjk(i,j,k) integer*4 i,j,k j=i+1 if(j.gt.3)j=1 k=j+1 if(k.gt.3)k=1 return end subroutine groa(ram,alr,ali,gr,gi,ar,ai,ENM,iwr) implicit none integer*4 I,J,K,L,iwr real*8 alr(3,3),ali(3,3),gr(3,3),gi(3,3),ar(3,3,3),ai(3,3,3), 1ram(*),SAL0,SAL1,SAG0,SAG1,SA1,SPALr,SPALi,EPS, 1S180,D180,D,YDX,YDY,P1,doc,AMU,BOHR,SpA3r,SpA3i,SpA32, 1gpisvejc,beta2,betaa,betag,roa3,roa1,ram1, 1EXCA,ENM EXCA=ENM*10.0d0 AMU=1822.0d0 BOHR=0.529177d0 c gpisvejc=980.0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 SPALr=0.0d0 SPALi=0.0d0 DO 3 I=1,3 SPALr=SPALr+alr(I,I) SPALi=SPALi+ali(I,I) DO 3 J=1,3 SAL0=SAL0+alr(I,I)*alr(J,J)+ali(I,I)*ali(J,J) SAL1=SAL1+alr(I,J)*alr(I,J)+ali(I,J)*ali(I,J) SAG0=SAG0+alr(I,I)*gr( J,J)+ali(I,I)*gi( J,J) SAG1=SAG1+alr(I,J)*gr( I,J)+ali(I,J)*gi( I,J) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*(alr(I,J)*ar(K,L,J)+ali(I,J)*ai(K,L,J))/3.0d0 if(iwr.ge.2)write(6,609)SAL0,SAL1,SAG0,SAG1,SA1 609 format(' SAL0 SAL1 SAG0 SAG1 SA1:',5g11.4) c IL+/-IR, backscattering S180=7.0d0*SAL1+SAL0 c write(6,*)ECM,SAL0,SAL1,SAG0,SAG1,SA1 D180=8.0d0*(3.0d0*SAG1-SAG0+SA1) D=D180 IF(S180.gt.1.0d-9)D=D/S180 c c polarized backscattering components (YDY+YDX=S180): YDY=6.0d0*(3.0d0*SAL1- SAL0)*(AMU*BOHR**4) YDX=6.0d0*(4.0d0*SAL1+2.0d0*SAL0)*(AMU*BOHR**4) if(YDX.gt.1.0d-9)then P1=YDY/YDX else P1=0.0d0 endif SpA3r=SPALr/3.0d0*dsqrt((AMU*BOHR**4)) SpA3i=SPALi/3.0d0*dsqrt((AMU*BOHR**4)) SpA32=SpA3r**2+SpA3i**2 beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) c c Barron's tensor invariants: c alpha G = [(1/3)Sp alpha] [(1/3)Sp G'] etc c beta(A)^2=Delta2 in Gaussian: betaa=SA1*3.0d0/2.0d0 *gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0 *gpisvejc c c DCP 180 c ram3=24.0d0*beta2 roa3=8.0d0*(12.0d0*betag+4.0d0*betaa) c c ICP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA32+7.0d0*beta2) doc=(3*P1-1)/(P1+1) c Raman, parallel polarization: ram(1)=YDX c Raman, perpendicular poalrization: ram(2)=YDY c Raman, total: ram(3)=ram1 c Depolarization ratio: ram(4)=P1 c ROA SCP 180: ram(5)=roa1 c DOC: ram(6)=doc c DCP ROA 180: ram(7)=roa3 do 6 i=1,3 call setjk(i,j,k) c S0(x)=(1/2)(azz azz + ayy ayy + azy azy + azy azy) c S0(y)=(1/2)(axx axx + azz azz + axz axz + axz axz) c S0(z)=(1/2)(axx axx + ayy ayy + axy axy + axy axy) c jj jj kk kk jk jk jk jk c 8-10 oriented ram1 ram(7+i)=180.0d0*0.5d0*AMU*BOHR**4*( 1alr(j,j)*alr(j,j)+alr(k,k)*alr(k,k)+ 1alr(j,k)*alr(j,k)+alr(j,k)*alr(j,k)+ 1ali(j,j)*ali(j,j)+ali(k,k)*ali(k,k)+ 1ali(j,k)*ali(j,k)+ali(j,k)*ali(j,k)) c -2ayx(Gxy+Gyx)+(ayy-axx)(Gxx-Gyy) c kj jk kj kk jj jj kk c +(2w/3)(ayyAx zy-axxAy zx+ayxAx zx-axyAy zy) c kk j ik jj k ij kj j ij jk k ik c 11-13 oriented roa1 6 ram(10+i)=-180.0d0*gpisvejc*( 1-2.0d0*alr(j,k)*(gr(j,k)+gr(k,j)) 1-2.0d0*ali(j,k)*(gi(j,k)+gi(k,j)) 1+(alr(k,k)-alr(j,j))*(gr(j,j)-gr(k,k)) 1+(ali(k,k)-ali(j,j))*(gi(j,j)-gi(k,k)) 1+2.0d0/3.0d0* 1(alr(k,k)*Ar(j,i,k)-alr(j,j)*Ar(k,i,j) 1+alr(k,j)*Ar(j,i,j)-alr(j,k)*Ar(k,i,k) 1+ali(k,k)*Ai(j,i,k)-ali(j,j)*Ai(k,i,j) 1+ali(k,j)*Ai(j,i,j)-ali(j,k)*Ai(k,i,k))) RETURN END subroutine sspread(ecm,d1,d2,d3,d4,d5,d6,d7,d8, 1s,np,wsmin,wsmax,fwhm,lg) implicit none real*8 s(16,np),wsmin,wsmax,w,dw,ecm,pd,fwhm,sf,d1, 1d2,d3,d4,d5,d6,d7,d8,dd(8),prod integer*4 i,np,j logical lg dd(1)=d1*108.7d0 dd(2)=d2*108.7d0 dd(3)=d3*108.7d0 dd(4)=d4*108.7d0 dd(5)=d5*435.0d0 dd(6)=d6*435.0d0 dd(7)=d7*435.0d0 dd(8)=d8*435.0d0 dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw if(lg)then pd=fwhm/2.0d0/dsqrt(log(2.0d0)) else pd=fwhm/2.0d0 endif do 1 i=1,np w=w+dw prod=w*sf(pd,lg,w,ecm,1) do 1 j=1,8 1 s(j,i)=s(j,i)+prod*dd(j) return end subroutine tq(qq,nat,x,PP,is) c is=1 : set quadrupole derivatives to local system c is=-1: put them to common origin implicit none integer*4 nat,is,L,I,J,LJ real*8 qq(3*nat,6),x(*),BOHR,PP(nat,3,3),R(3),s BOHR=0.52917705993d0 if(is.eq.1)then s=1.0d0 write(6,*)' QD set to common system' else if(is.eq.-1)then s=-1.0d0 write(6,*)' QD set to local systems' else write(6,*)is call report('tq unknown option') endif endif DO 2204 L=1,nat DO 3 I=1,3 3 R(I)=X(I+3*(L-1))/BOHR DO 2204 J=1,3 LJ=J+3*(L-1) s=1.0d0 qq(LJ,1)=qq(LJ,1)+(PP(L,J,1)*R(1)+PP(L,J,1)*R(1))*s qq(LJ,2)=qq(LJ,2)+(PP(L,J,2)*R(2)+PP(L,J,2)*R(2))*s qq(LJ,3)=qq(LJ,3)+(PP(L,J,3)*R(3)+PP(L,J,3)*R(3))*s qq(LJ,4)=qq(LJ,4)+(PP(L,J,1)*R(2)+PP(L,J,2)*R(1))*s qq(LJ,5)=qq(LJ,5)+(PP(L,J,1)*R(3)+PP(L,J,3)*R(1))*s 2204 qq(LJ,6)=qq(LJ,6)+(PP(L,J,2)*R(3)+PP(L,J,3)*R(2))*s return end subroutine reqd(lqpt,qq,nat,x,PP) implicit none integer*4 i,ia,ix,nat real*8 qq(3*nat,6),x(*),PP(nat,3,3) logical lqpt if(lqpt)then open(9,file='FILE.QEN',status='old') do 1 i=1,9 1 read(9,*) do 2 ia=1,nat do 2 ix=1,3 2 read(9,93)(qq(ix+3*(ia-1),i),i=1,6) 93 format(9x,3e13.4) close(9) write(6,*)'QD read' call tq(qq,nat,x,PP,-1) else call vz2(qq,3*nat,6) write(6,*)' local quadrupole derivatives initiated' endif return end subroutine wspectra(f,wsmin,wsmax,s,np,j) implicit none real*8 s(16,np),wsmin,wsmax,w,dw integer*4 i,np,j character*(*) f OPEN(30,FILE=f) dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw do 1 i=1,np w=w+dw 1 write(30,300)w,s(j,i) 300 format(f12.2,' ',g12.4) close(30) write(6,*)f//' written' return end subroutine sm(sr,si,m,n) implicit none integer*4 n,i,j,ia real*8 sr(n,n),si(n,n),m(*) do 1 i=1,n ia=(i+2)/3 do 1 j=1,n si(i,j)=si(i,j)/dsqrt(m(ia)) 1 sr(i,j)=sr(i,j)/dsqrt(m(ia)) return end SUBROUTINE MASSFF(N,f,m) IMPLICIT none integer*4 i,j,N REAL*8 m(*),F(N,N) do 1 i=1,N do 1 j=1,N 1 f(i,j)=f(i,j)/dsqrt(m((i+2)/3)*m((j+2)/3)) WRITE(*,*)' Mass-scaled FF done.' return end subroutine setis(s,is) character*(*) s integer*4 is do 1 is=1,len(s) 1 if(s(is:is).ne.' ')return end subroutine massa(natu,m,mu) implicit none integer*4 natu,ia real*8 m(*),mu(*) do 1 ia=1,natu 1 mu(ia)=m(ia) return end subroutine maked(ice,n,n3,f,Dc,Ds,ri,rj,rk,NI,NJ,NK,n27,ip,lcub) implicit none real*8 f(n,n),Dc(n3,n3),Ds(n3,n3),pi,fi,fj,fk,sJ,cJ,mJ,ri,rj,rk integer*4 ice,n,n3,NI,NJ,NK,n27,ip(n27,3),J,iof,jof,nat1,i1,j1,k1, 1ei,ej,ek logical lcub nat1=n3/3 pi=4.0d0*datan(1.0d0) fi=2.0d0*pi*ri/dble(NI) fj=2.0d0*pi*rj/dble(NJ) fk=2.0d0*pi*rk/dble(NK) ei=ip(ice,1) ej=ip(ice,2) ek=ip(ice,3) iof=nat1*(ice-1) call vz2(Dc,n3,n3) call vz2(Ds,n3,n3) c loop over all available elementary cells: do 1 J=1,n27 jof=nat1*(J -1) c relative position to ice: i1=ip(J,1)-ei j1=ip(J,2)-ej k1=ip(J,3)-ek mJ=fi*dble(i1)+fj*dble(j1)+fk*dble(k1) sJ=sin(mJ) cJ=cos(mJ) c exp( i 2 pi q I / N )/sqrt(N) on I = 0 ... N-1 1 if(.not.lcub.or.(i1.ge.0.and.j1.ge.0.and.k1.ge.0)) 1call addd(sJ,cJ,n,n3,nat1,iof,jof,Dc,Ds,f) return end subroutine addd(sJ,cJ,n,n3,nat1,iof,jof,Dc,Ds,f) implicit none integer*4 n,n3,nat1,iof,jof,ia,ix,iu,ii,ja,jx,ju,jj real*8 sJ,cJ,Dc(n3,n3),Ds(n3,n3),f(n,n) do 1 ia=1,nat1 do 1 ix=1,3 c in D: iu=ix+3*(ia -1) c in original ice: ii=ix+3*(ia+iof-1) do 1 ja=1,nat1 do 1 jx=1,3 ju=jx+3*(ja -1) jj=jx+3*(ja+jof-1) Dc(iu,ju)=Dc(iu,ju)+f(ii,jj)*cJ 1 Ds(iu,ju)=Ds(iu,ju)+f(ii,jj)*sJ return end SUBROUTINE READFF(N,FCAR) IMPLICIT none integer*4 N,N1,N3,LN,I,J real*8 FCAR(N,N),CONST OPEN(20,FILE='FILE.FC',STATUS='OLD') N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in.' CLOSE(20) RETURN END function iccub(nat,r,n27,ip,wi,wj,wk) c determine central cube implicit none integer*4 nat,nat1,ix,iu,ii,iccub,ice,i,n27,ip(n27,3) real*8 cet(3),r(*),xi,yi,zi,d,wi(3),wj(3),wk(3),down, 1vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,sp,vu(3),a,b,c real*8,allocatable::cei(:,:) allocate(cei(n27,3)) nat1=nat/n27 c center of all and of each one: do 40 ix=1,3 cet(ix)=0.0d0 do 42 iu=1,n27 cei(iu,ix)=0.0d0 ii=nat1*(iu-1) do 41 i=ii+1,ii+nat1 cet( ix)=cet( ix)+r(3*(i-1)+ix) 41 cei(iu,ix)=cei(iu,ix)+r(3*(i-1)+ix) 42 cei(iu,ix)=cei(iu,ix)/dble(nat1) 40 cet(ix)=cet(ix)/dble(nat) write(6,600)cet 600 format(/,' CRYSTAL.X center (A):',7x,3f7.2,/,' Cube centers (A):') write(6,601)(iu,(cei(iu,ix),ix=1,3),iu=1,n27) 601 format(3(i3,':',3f7.2)) d=1.0d10 ice=0 do 43 iu=1,n27 xi=cei(iu,1)-cet(1) yi=cei(iu,2)-cet(2) zi=cei(iu,3)-cet(3) if(dsqrt(xi**2+yi**2+zi**2).lt.d)then d=dsqrt(xi**2+yi**2+zi**2) ice=iu endif 43 continue if(ice.eq.0)then call report('Central cube not found') else write(6,602)ice,d 602 format(' Central cube:',i3,f10.2,' A from the center') iccub=ice endif c relative positions to the central cell do 1 iu=1,n27 do 2 ix=1,3 2 vu(ix)=cei(iu,ix)-cei(ice,ix) vi=sp(vu,wi) vj=sp(vu,wj) vk=sp(vu,wk) vii=sp(wi,wi) vij=sp(wi,wj) vik=sp(wi,wk) vjj=sp(wj,wj) vjk=sp(wj,wk) vkk=sp(wk,wk) down=(vkk*vii-vik**2)*(vjj*vii-vij**2) 1-(vjk*vii-vij*vik)*(vjk*vii-vik*vij) c write(6,*)vi,vj,vk,vii,vij,vik,vjj,vjk,vkk,down c=(vi*(vjk*vij-vjj*vik)+vj*(vij*vik-vjk*vii)+vk*(vjj*vii-vij**2))* 1vii/down b=(vj*vii-vi*vij-c*(vjk*vii-vik*vij))/(vjj*vii-vij**2) a=(vi-b*vij-c*vik)/vii write(6,603)iu,a,b,c 603 format(i3,3f5.1,$) if(mod(iu,4).eq.0)write(6,*) ip(iu,1)=nint(a) ip(iu,2)=nint(b) 1 ip(iu,3)=nint(c) if(mod(n27,4).ne.0)write(6,*) return end function sp(v1,v2) real*8 v1(*),v2(*),sp sp=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3) return end C SUBROUTINE READTEN(nat,P,A) IMPLICIT none integer*4 nat,L,I,J REAL*8 P(nat,3,3),A(nat,3,3) open(15,file='FILE.TEN',status='old') READ (15,*) C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) C DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) WRITE(6,*)' VCD and IR parameters read in ...' C close(15) RETURN END C SUBROUTINE report(s) character*(*) s write(6,*)s stop end FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END subroutine spreadr(ecm,ram,s,np,wsmin,wsmax,fwhm,lg,temp) implicit none real*8 s(16,np),wsmin,wsmax,w,dw,ecm,pd,fwhm,sf,ram(*), 1dd(8),prod,temp,c integer*4 i,np,j logical lg c=1.0d0/ecm/(1.0d0-exp(-(1.44d0*ecm/temp))) dd(1)=ram( 3)*c dd(2)=ram( 8)*c dd(3)=ram( 9)*c dd(4)=ram(10)*c dd(5)=ram( 5)*c dd(6)=ram(11)*c dd(7)=ram(12)*c dd(8)=ram(13)*c dw=(wsmax-wsmin)/dble(np-1) w=wsmin-dw if(lg)then pd=fwhm/2.0d0/dsqrt(log(2.0d0)) else pd=fwhm/2.0d0 endif do 1 i=1,np w=w+dw prod=sf(pd,lg,w,ecm,1) do 1 j=1,8 1 s(8+j,i)=s(8+j,i)+prod*dd(j) return end 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 SUBROUTINE READTTT(ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) open(15,file='FILE.TTT',status='old') READ(15,*) READ(15,*)NATC if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file') READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) RETURN END subroutine m3(ar,ai,br,bi,cr,ci,sr,si) c multiply 3 complex numbers a b c = s implicit none real*8 ar,ai,br,bi,cr,ci,sr,si,tr,ti c a b: tr=ar*br-ai*bi ti=ai*br+ar*bi c a b c: sr=tr*cr-ti*ci si=ti*cr+tr*ci return end subroutine readopt(NI,NJ,NK,ltab,lttt,lqpt,ldip,np,n27,wsmin, 1wsmax,temp,fwhm,lg,lcub,lone,delq,ENM,IWR) implicit none logical lex,ltab,lttt,lqpt,lg,ldip,lcub,lone integer*4 NI,NJ,NK,np,n27,IWR real*8 wsmin,wsmax,temp,fwhm,delq,ENM character*4 key NI=10 NJ=10 NK=10 delq=1.0d0 np=3901 n27=27 ENM=532.0d0 IWR=0 wsmin=100.0d0 fwhm=10.0d0 wsmax=4000.0d0 temp=293.0d0 c delete border interactions: lcub=.false. c signal from central cube only: lone=.false. ltab=.true. lttt=.true. lqpt=.false. ldip=.false. lg=.false. inquire(file='XQQ.OPT',exist=lex) if(lex)then open(9,file='XQQ.OPT') 1 read(9,90,end=99,err=99)key 90 format(a4) if(key.eq.'DELQ')read(9,*)delq if(key.eq.'EXCN')read(9,*)ENM if(key.eq.'FWHM')read(9,*)fwhm if(key.eq.'NINJ')read(9,*)NI,NJ,NK if(key.eq.'WRIT')read(9,*)IWR if(key.eq.'LCUB')read(9,*)lcub if(key.eq.'LGAU')read(9,*)lg if(key.eq.'LDIP')read(9,*)ldip if(key.eq.'LONE')read(9,*)lone if(key.eq.'LQPT')read(9,*)lqpt if(key.eq.'LTAB')read(9,*)ltab if(key.eq.'LTTT')read(9,*)lttt if(key.eq.'NPOI')read(9,*)np if(key.eq.'TEMP')read(9,*)temp if(key.eq.'WSMI')read(9,*)wsmin if(key.eq.'WSMA')read(9,*)wsmax if(key(1:3).eq.'N27')read(9,*)n27 goto 1 99 close(9) endif return end subroutine READM(nat,m) implicit none integer*4 i,nat real*8 m(*) open(9,file='FTRY.INP', status='old') read(9,*) read(9,*) read(9,90)(m(i),i=1,nat) 90 format(6F12.6) read(9,*) read(9,90)(m(i),i=1,nat) close(9) WRITE(*,*)' Masses read in.' return end subroutine vz2(v,n,m) implicit none real*8 v(n,m) integer*4 n,m,i,j do 1 i=1,n do 1 j=1,m 1 v(i,j)=0.0d0 return end subroutine rwr(n,ice,nat1,nat,P1,A1,Q1,P,A,qq,al1,gp1,aa1, 1al,gp,aa) implicit none integer*4 ice,nat1,nat,iof,i,ix,ii1,ii,jx,kx,lx,n real*8 P(nat,3,3),A(nat,3,3),al(n,3,3),gp(n,3,3), 1aa(n,3,3,3),P1(nat1*3,3),A1(nat1*3,3), 1al1(3*nat1,3,3),gp1(3*nat1,3,3),aa1(3*nat1,3,3,3), 1Q1(nat1*3,3,3),qq(3*nat,6) iof=nat1*(ice-1) do 1 i=1,nat1 do 1 ix=1,3 ii1=ix+3*(i -1) ii =ix+3*(i+iof-1) Q1(ii1,1,1)=qq(i+iof,1) Q1(ii1,2,2)=qq(i+iof,2) Q1(ii1,3,3)=qq(i+iof,3) Q1(ii1,1,2)=qq(i+iof,4) Q1(ii1,2,1)=qq(i+iof,4) Q1(ii1,1,3)=qq(i+iof,5) Q1(ii1,3,1)=qq(i+iof,5) Q1(ii1,2,3)=qq(i+iof,6) Q1(ii1,3,2)=qq(i+iof,6) do 1 jx=1,3 P1(ii1,jx)=P(i+iof,ix,jx) A1(ii1,jx)=A(i+iof,ix,jx) do 1 kx=1,3 al1(ii1,jx,kx)=al(ii,jx,kx) gp1( ii1,jx,kx)=gp( ii,jx,kx) do 1 lx=1,3 1 aa1( ii1,jx,kx,lx)=aa( ii,jx,kx,lx) return end c ============================================================ subroutine tql2(nm,n,d,e,z,ierr) c integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr double precision d(n),e(n),z(nm,n) double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure tql2, c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and c wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a symmetric tridiagonal matrix by the ql method. c the eigenvectors of a full symmetric matrix can also c be found if tred2 has been used to reduce this c full matrix to tridiagonal form. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c z contains the transformation matrix produced in the c reduction by tred2, if performed. if the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,ierr-1. c c e has been destroyed. c c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. if an error exit is made, c z contains the eigenvectors associated with the stored c eigenvalues. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c f = 0.0d0 tst1 = 0.0d0 e(n) = 0.0d0 c do 240 l = 1, n j = 0 h = dabs(d(l)) + dabs(e(l)) if (tst1 .lt. h) tst1 = h c .......... look for small sub-diagonal element .......... do 110 m = l, n tst2 = tst1 + dabs(e(m)) if (tst2 .eq. tst1) go to 120 c .......... e(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 220 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 l2 = l1 + 1 g = d(l) p = (d(l1) - g) / (2.0d0 * e(l)) r = pythag(p,1.0d0) d(l) = e(l) / (p + dsign(r,p)) d(l1) = e(l) * (p + dsign(r,p)) dl1 = d(l1) h = g - d(l) if (l2 .gt. n) go to 145 c do 140 i = l2, n 140 d(i) = d(i) - h c 145 f = f + h c .......... ql transformation .......... p = d(m) c = 1.0d0 c2 = c el1 = e(l1) s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml c3 = c2 c2 = c s2 = s i = m - ii g = c * e(i) h = c * p r = pythag(p,e(i)) e(i+1) = s * r s = e(i) / r c = p / r p = c * d(i) - s * g d(i+1) = h + s * (c * g + s * d(i)) c .......... form vector .......... do 180 k = 1, n h = z(k,i+1) z(k,i+1) = s * z(k,i) + c * h z(k,i) = c * z(k,i) - s * h 180 continue c 200 continue c p = -s * s2 * c3 * el1 * e(l) / dl1 e(l) = s * p d(l) = c * p tst2 = tst1 + dabs(e(l)) if (tst2 .gt. tst1) go to 130 220 d(l) = d(l) + f 240 continue c .......... order eigenvalues and eigenvectors .......... do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p c do 280 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 280 continue c 300 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end c ============================================================ subroutine htridi(nm,n,ar,ai,d,e,e2,tau) c integer i,j,k,l,n,ii,nm,jp1 double precision ar(nm,n),ai(nm,n),d(n),e(n),e2(n),tau(2,n) double precision f,g,h,fi,gi,hh,si,scale,pythag c c this subroutine is a translation of a complex analogue of c the algol procedure tred1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine reduces a complex hermitian matrix c to a real symmetric tridiagonal matrix using c unitary similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian input matrix. c only the lower triangle of the matrix need be supplied. c c on output c c ar and ai contain information about the unitary trans- c formations used in the reduction in their full lower c triangles. their strict upper triangles and the c diagonal of ar are unaltered. c c d contains the diagonal elements of the the tridiagonal matrix. c c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero. c c e2 contains the squares of the corresponding elements of e. c e2 may coincide with e if the squares are not needed. c c tau contains further information about the transformations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c tau(1,n) = 1.0d0 tau(2,n) = 0.0d0 c do 100 i = 1, n 100 d(i) = ar(i,i) c .......... for i=n step -1 until 1 do -- .......... do 300 ii = 1, n i = n + 1 - ii l = i - 1 h = 0.0d0 scale = 0.0d0 if (l .lt. 1) go to 130 c .......... scale row (algol tol then not needed) .......... do 120 k = 1, l 120 scale = scale + dabs(ar(i,k)) + dabs(ai(i,k)) c if (scale .ne. 0.0d0) go to 140 tau(1,l) = 1.0d0 tau(2,l) = 0.0d0 130 e(i) = 0.0d0 e2(i) = 0.0d0 go to 290 c 140 do 150 k = 1, l ar(i,k) = ar(i,k) / scale ai(i,k) = ai(i,k) / scale h = h + ar(i,k) * ar(i,k) + ai(i,k) * ai(i,k) 150 continue c e2(i) = scale * scale * h g = dsqrt(h) e(i) = scale * g f = pythag(ar(i,l),ai(i,l)) c .......... form next diagonal element of matrix t .......... if (f .eq. 0.0d0) go to 160 tau(1,l) = (ai(i,l) * tau(2,i) - ar(i,l) * tau(1,i)) / f si = (ar(i,l) * tau(2,i) + ai(i,l) * tau(1,i)) / f h = h + f * g g = 1.0d0 + g / f ar(i,l) = g * ar(i,l) ai(i,l) = g * ai(i,l) if (l .eq. 1) go to 270 go to 170 160 tau(1,l) = -tau(1,i) si = tau(2,i) ar(i,l) = g 170 f = 0.0d0 c do 240 j = 1, l g = 0.0d0 gi = 0.0d0 c .......... form element of a*u .......... do 180 k = 1, j g = g + ar(j,k) * ar(i,k) + ai(j,k) * ai(i,k) gi = gi - ar(j,k) * ai(i,k) + ai(j,k) * ar(i,k) 180 continue c jp1 = j + 1 if (l .lt. jp1) go to 220 c do 200 k = jp1, l g = g + ar(k,j) * ar(i,k) - ai(k,j) * ai(i,k) gi = gi - ar(k,j) * ai(i,k) - ai(k,j) * ar(i,k) 200 continue c .......... form element of p .......... 220 e(j) = g / h tau(2,j) = gi / h f = f + e(j) * ar(i,j) - tau(2,j) * ai(i,j) 240 continue c hh = f / (h + h) c .......... form reduced a .......... do 260 j = 1, l f = ar(i,j) g = e(j) - hh * f e(j) = g fi = -ai(i,j) gi = tau(2,j) - hh * fi tau(2,j) = -gi c do 260 k = 1, j ar(j,k) = ar(j,k) - f * e(k) - g * ar(i,k) x + fi * tau(2,k) + gi * ai(i,k) ai(j,k) = ai(j,k) - f * tau(2,k) - g * ai(i,k) x - fi * e(k) - gi * ar(i,k) 260 continue c 270 do 280 k = 1, l ar(i,k) = scale * ar(i,k) ai(i,k) = scale * ai(i,k) 280 continue c tau(2,l) = -si 290 hh = d(i) d(i) = ar(i,i) ar(i,i) = hh ai(i,i) = scale * dsqrt(h) 300 continue c return end c ============================================================ subroutine tqlrat(n,d,e2,ierr) c integer i,j,l,m,n,ii,l1,mml,ierr double precision d(n),e2(n) double precision b,c,f,g,h,p,r,s,t,epslon,pythag c c this subroutine is a translation of the algol procedure tqlrat, c algorithm 464, comm. acm 16, 689(1973) by reinsch. c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the rational ql method. c c on input c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e2 contains the squares of the subdiagonal elements of the c input matrix in its last n-1 positions. e2(1) is arbitrary. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c e2 has been destroyed. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e2(i-1) = e2(i) c f = 0.0d0 t = 0.0d0 e2(n) = 0.0d0 c do 290 l = 1, n j = 0 h = dabs(d(l)) + dsqrt(e2(l)) if (t .gt. h) go to 105 t = h b = epslon(t) c = b * b c .......... look for small squared sub-diagonal element .......... 105 do 110 m = l, n if (e2(m) .le. c) go to 120 c .......... e2(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 210 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 s = dsqrt(e2(l)) g = d(l) p = (d(l1) - g) / (2.0d0 * s) r = pythag(p,1.0d0) d(l) = s / (p + dsign(r,p)) h = g - d(l) c do 140 i = l1, n 140 d(i) = d(i) - h c f = f + h c .......... rational ql transformation .......... g = d(m) if (g .eq. 0.0d0) g = b h = g s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii p = g * h r = p + e2(i) e2(i+1) = s * r s = e2(i) / r d(i+1) = h + s * (h + d(i)) g = d(i) - e2(i) / g if (g .eq. 0.0d0) g = b h = g * p / r 200 continue c e2(l) = s * g d(l) = h c .......... guard against underflow in convergence test .......... if (h .eq. 0.0d0) go to 210 if (dabs(e2(l)) .le. dabs(c/h)) go to 210 e2(l) = h * e2(l) if (e2(l) .ne. 0.0d0) go to 130 210 p = d(l) + f c .......... order eigenvalues .......... if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. d(i-1)) go to 270 d(i) = d(i-1) 230 continue c 250 i = 1 270 d(i) = p 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end c ============================================================ subroutine htribk(nm,n,ar,ai,tau,m,zr,zi) c integer i,j,k,l,m,n,nm double precision ar(nm,n),ai(nm,n),tau(2,n),zr(nm,m),zi(nm,m) double precision h,s,si c c this subroutine is a translation of a complex analogue of c the algol procedure trbak1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine forms the eigenvectors of a complex hermitian c matrix by back transforming those of the corresponding c real symmetric tridiagonal matrix determined by htridi. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain information about the unitary trans- c formations used in the reduction by htridi in their c full lower triangles except for the diagonal of ar. c c tau contains further information about the transformations. c c m is the number of eigenvectors to be back transformed. c c zr contains the eigenvectors to be back transformed c in its first m columns. c c on output c c zr and zi contain the real and imaginary parts, c respectively, of the transformed eigenvectors c in their first m columns. c c note that the last component of each returned vector c is real and that vector euclidean norms are preserved. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (m .eq. 0) go to 200 c .......... transform the eigenvectors of the real symmetric c tridiagonal matrix to those of the hermitian c tridiagonal matrix. .......... do 50 k = 1, n c do 50 j = 1, m zi(k,j) = -zr(k,j) * tau(2,k) zr(k,j) = zr(k,j) * tau(1,k) 50 continue c if (n .eq. 1) go to 200 c .......... recover and apply the householder matrices .......... do 140 i = 2, n l = i - 1 h = ai(i,i) if (h .eq. 0.0d0) go to 140 c do 130 j = 1, m s = 0.0d0 si = 0.0d0 c do 110 k = 1, l s = s + ar(i,k) * zr(k,j) - ai(i,k) * zi(k,j) si = si + ar(i,k) * zi(k,j) + ai(i,k) * zr(k,j) 110 continue c .......... double divisions avoid possible underflow .......... s = (s / h) / h si = (si / h) / h c do 120 k = 1, l zr(k,j) = zr(k,j) - s * ar(i,k) - si * ai(i,k) zi(k,j) = zi(k,j) - si * ar(i,k) + s * ai(i,k) 120 continue c 130 continue c 140 continue c 200 return end c ============================================================ subroutine ch(n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) c integer i,j,n,nm,ierr,matz double precision ar(n,n),ai(n,n),w(n),zr(n,n),zi(n,n), x fv1(n),fv2(n),fm1(2,n) c c this subroutine calls the recommended sequence of c subroutines from the eigensystem subroutine package (eispack) c to find the eigenvalues and eigenvectors (if desired) c of a complex hermitian matrix. c c on input c c Petr Bour Change: nm=0 by default! c c nm must be set to the row dimension of the two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix a=(ar,ai). c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian matrix. c c matz is an integer variable set equal to zero if c only eigenvalues are desired. otherwise it is set to c any non-zero integer for both eigenvalues and eigenvectors. c c on output c c w contains the eigenvalues in ascending order. c c zr and zi contain the real and imaginary parts, c respectively, of the eigenvectors if matz is not zero. c c ierr is an integer output variable set equal to an error c completion code described in the documentation for tqlrat c and tql2. the normal completion code is zero. c c fv1, fv2, and fm1 are temporary storage arrays. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c nm=n if (n .le. nm) go to 10 ierr = 10 * n go to 50 c 10 call htridi(nm,n,ar,ai,w,fv1,fv2,fm1) if (matz .ne. 0) go to 20 c .......... find eigenvalues only .......... call tqlrat(n,w,fv2,ierr) go to 50 c .......... find both eigenvalues and eigenvectors .......... 20 do 40 i = 1, n c do 30 j = 1, n zr(j,i) = 0.0d0 30 continue c zr(i,i) = 1.0d0 40 continue c call tql2(nm,n,w,fv1,zr,ierr) if (ierr .ne. 0) go to 50 call htribk(nm,n,ar,ai,fm1,n,zr,zi) 50 return end double precision function epslon (x) double precision x c c estimate unit roundoff in quantities of size x. c double precision a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c a = 4.0d0/3.0d0 10 b = a - 1.0d0 c = b + b + b eps = dabs(c-1.0d0) if (eps .eq. 0.0d0) go to 10 epslon = eps*dabs(x) return end c ============================================================ double precision function pythag(a,b) double precision a,b c c finds dsqrt(a**2+b**2) without overflow or destructive underflow c double precision p,r,s,t,u p = dmax1(dabs(a),dabs(b)) if (p .eq. 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r if (t .eq. 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p r = (s/u)**2 * r go to 10 20 pythag = p return end