program correctpol implicit none integer*4 nat,N,i,j,IERR,ia,is real*8 epsr,CM,gau,wmin,wmax,v(4),ux(4),uy(4),rmax real*8,allocatable::r(:),p(:,:),f(:,:),ZM(:),S(:,:),W(:),W2(:), 1Hr(:,:),Hi(:,:),Cr(:,:),Ci(:,:),Eq(:),u(:,:) integer*4,allocatable::q(:),coa(:,:) parameter (CM=219474.0d0) write(6,800) 800 format(/,' Corrects FILE.FC by vibrational polarizabilities',/,/, 1' Input: FILE.FC',/, 1' FILE.TEN',/, 1' FILE.X ',/, 1' Output: NEW.FC',/) call readopt(gau,epsr,wmin,wmax,v,ux,uy,is,rmax) allocate(r(1),q(1)) call rx(0,'FILE.X',nat,r,q) deallocate(r,q) N=3*nat allocate(q(nat),r(N),p(2*N,3),F(N,N),ZM(N),S(N,N),W(N),W2(N), 1Hr(N,N),Hi(N,N),Cr(N,N),Ci(N,N),Eq(N),u(2*N,3)) call rx(1,'FILE.X',nat,r,q) write(6,*)'Geometry read' CALL READFF(N,F) call readM(nat,ZM) CALL FAST(N,F,ZM) CALL TRED12(N,F,S,W2,2,IERR) IF(IERR.NE.0)call report('diagonalization not ok') do 1 i=1,N if(W2(i).gt.0)then W(i)= sqrt(W2(i)) else W(i)=0.0d0 endif 1 continue write(6,601) 601 format(' Vibrational frequencies:') write(6,602)(W(I)*CM,I=1,N) 602 format(8f10.2) C Calculate the true mass-weighted S-matrix: DO 210 I=1,N DO 210 J=1,N 210 S(I,J)=S(I,J)*ZM(I) c read APT and AAT to one variable p: call rdp(N,p) c make transition dipole moments: call new66(N,p,S,u,W) c recognize amides call ra(ia,r,q,nat) allocate(coa(ia,6)) call rda(ia,coa) c set potential call ddm(N,ia,coa,S,Hr,Hi,W,r,p,gau,wmin,wmax, 1epsr,v,ux,uy,is,rmax) call minimax(Hr,N,'Dr:') call minimax(Hi,N,'Di:') c diagonalizc diagonalize: write(6,605)N,N 605 format(' Diagonalizing',i6,' x',i6) call ch(N,Hr,Hi,Eq,1,Cr,Ci) call new6c(N,Cr,Ci,Eq,u) end subroutine minimax(V,N,s) implicit none character*(*) s integer*4 N,I,J,imax,ijimax,ijjmax real*8 V(N,N),vimax,vijmax,CM CM=219474.0d0 vimax=abs(v(1,1)) imax=1 vijmax=abs(v(1,2)) ijimax=1 ijjmax=2 do 1 i=1,N if(vimax.lt.dabs(v(i,i)))then vimax=abs(v(i,i)) imax=i endif do 1 j=i+1,N if(vijmax.lt.dabs(v(i,j)))then vijmax=dabs(v(i,j)) ijimax=i ijjmax=j endif 1 continue write(6,*)s write(6,659)imax,imax,vimax*CM,ijimax,ijjmax,vijmax*CM 659 format(' Maximal diagonal element ',2i6,f12.2,' cm-1',/, 1 ' Maximal off-diagonal element ',2i6,f12.2,' cm-1') return end subroutine readopt(gau,epsr,wmin,wmax,v,ux,uy,is,rmax) implicit none integer*4 is,i real*8 gau,epsr,CM,gcm,wmin,wmax,v(4),ux(4),uy(4), 1vcm(4),rmax logical ldeut,lex character*4 key CM=219474.0d0 epsr=3.0d0 ldeut=.false. gcm=20.0d0 wmin=1400.0d0 wmax=1800.0d0 c Amide I H: c vibrational frequency in cm-1: vcm(1)=1675.0d0 c projection of dipole to x ~ C->N bond direction c z ~ CN x CO, y ~ z x x ux(1)=-0.101d0 c projection of dipole to y-direction: uy(1)=+0.103d0 c Amide II H: vcm(2)=1552.0d0 ux(2)=-0.117d0 uy(2)=-0.024d0 c Amide I D: vcm(3)=1667.0d0 ux(3)=-0.116d0 uy(3)=+0.100d0 c Amide II D: vcm(4)=1482.0d0 ux(4)= 0.058d0 uy(4)=+0.054d0 c up to where not to add it: rmax=4.0d0 inquire(file='CORRECTPOL.OPT',exist=lex) if(lex)then open(9,file='CORRECTPOL.OPT') 1 read(9,90,end=99,err=99)key 90 format(a4) if(key.eq.'EPSR')read(9,*)epsr if(key.eq.'GAMM')read(9,*)gcm if(key.eq.'DEUT')read(9,*)ldeut if(key.eq.'WMIN')read(9,*)wmin if(key.eq.'WMAX')read(9,*)wmax if(key.eq.'RMAX')read(9,*)rmax goto 1 99 close(9) endif do 2 i=1,4 2 v(i)=vcm(i)/CM gau=gcm/CM write(6,609)epsr 609 format(/,' Epsr = ',f6.2) if(ldeut)then is=2 else is=0 endif rmax=rmax/0.529177d0 return end subroutine new6c(N,Cr,Ci,E,u) implicit none integer*4 N,i,ix,ii real*8 Cr(N,N),Ci(N,N),u(2*N,3),CM,debye,E(*),mr,mi,dr,di,R,D parameter (CM=219474.0d0,debye=2.541580253d0) open(9,file='DOGC.TAB') write(9,900) 900 format(' Dipole and rotational strengths from correctpol',/,/,/) do 1 i=1,N if(E(i).gt.0.0d0)then D=0.0d0 R=0.0d0 do 2 ix=1,3 dr=0.0d0 mr=0.0d0 di=0.0d0 mi=0.0d0 do 3 ii=1,N mr=mr+Cr(ii,i)*u(N+ii,ix) mi=mi+Ci(ii,i)*u(N+ii,ix) dr=dr+Cr(ii,i)*u( ii,ix) 3 di=di+Ci(ii,i)*u( ii,ix) D=D+dr*dr+di*di 2 R=R+dr*mr+di*mi write(9,902)i,E(i)*CM,D*debye**2,R*debye**2 endif 1 continue 902 format(i6,f10.2,2e15.4) write(9,903) 903 format(60(1h-)) close(9) return end subroutine new66(N,p,S,u,W) implicit none integer*4 N,i,ix,ii real*8 S(N,N),p(2*N,3),W(*),CM,debye,u(2*N,3),D,R,clight parameter (CM=219474.0d0,debye=2.541580253d0,clight=137.03604d0) u=0.0d0 open(9,file='DOG0.TAB') write(9,900) 900 format(' Dipole and rotational strengths from correctpol',/,/,/) do 1 i=1,N if(W(i).gt.0.0d0)then D=0.0d0 R=0.0d0 do 2 ix=1,3 u( i,ix)=0.0d0 u(N+i,ix)=0.0d0 do 3 ii=1,N u(N+i,ix)=u(N+i,ix)+S(ii,i)*p(N+ii,ix) 3 u( i,ix)=u( i,ix)+S(ii,i)*p( ii,ix) u( i,ix)=u( i,ix)/dsqrt(2.0d0*W(i)) u(N+i,ix)=u(N+i,ix)*dsqrt(2.0d0*W(i))/clight D=D+u(i,ix)*u( i,ix) 2 R=R+u(i,ix)*u(N+i,ix) write(9,902)i,W(i)*CM,D*debye**2,R*debye**2 902 format(i6,f10.2,2e15.4) endif 1 continue write(9,903) 903 format(60(1h-)) close(9) return end subroutine ddm(N,ia,coa,S,Dr,Di,W,r,p,gau,wmin,wmax, 1epsr,v,ux,uy,is,rmax) implicit none integer*4 N,ia,coa(ia,6),i,j,ii,i6,ai,ix,jx,is,kk, 1kmax,imax,iumax,kumax,ipmax,jpmax,kpmax real*8 W(*),r(*),CM,x(3),y(3),z(3),gau,wmin,wmax,tt,wij, 1co(3),v(4),alr(3,3),ali(3,3),rk(3),tik(3,3), 1p(2*N,3),S(N,N),epsr,Dr(N,N),Di(N,N),ri(3),ux(4),uy(4), 1tmax,umax,dip,pmax,dik,rmax,alimax(3,3),alrmax(3,3) real*8,allocatable::tau(:,:),ui(:,:),u(:,:,:) Dr=0.0d0 Di=0.0d0 CM=219474.0d0 do 1 i=1,N 1 Dr(i,i)=W(i) write(6,*)'seting Hamiltonian matrix' c pre-calculate tau: allocate(tau(N,ia*3),ui(ia,3),u(ia,4,3)) tau=0.0d0 tmax=0.0d0 umax=0.0d0 kmax=0 imax=0 iumax=0 kumax=0 alimax=0.0d0 alrmax=0.0d0 c loop over normal modes: do 8 I=1,N if(W(I)*CM.le.wmax.and.W(I)*CM.ge.wmin)then c loop over amides ii do 111 ii=1,ia do 11 ix=1,3 ui(ii,ix)=0.0d0 c loop over amide atoms do 112 i6=1,6 ai=coa(ii,i6) 112 if(ai.ne.0)ui(ii,ix)=ui(ii,ix)+S(3*ai-2,I)*p(3*ai-2,ix) 1 +S(3*ai-1,I)*p(3*ai-1,ix) 1 +S(3*ai ,I)*p(3*ai ,ix) 11 ui(ii,ix)=ui(ii,ix)/dsqrt(2.0d0*W(I)) dip=dsqrt(ui(ii,1)**2+ui(ii,2)**2+ui(ii,3)**2) if(dip.gt.umax)then umax=dip iumax=I kumax=ii endif 111 continue c loop over amides kk: do 9 kk=1,ia c position of amide kk: call asx(rk,r,coa(kk,1)) c loop over amides ii: do 9 ii=1,ia if(ii.ne.kk)then call asx(ri,r,coa(ii,1)) c position tensor: call dtik(tik,ri,rk,epsr,dik) if(dik.gt.rmax)then do 901 ix=1,3 do 901 jx=1,3 if(abs(tik(ix,jx)).gt.tmax)then tmax=abs(tik(ix,jx)) kmax=kk imax=ii endif 901 continue do 903 ix=1,3 903 tau(I,3*(kk-1)+ix)=tau(I,3*(kk-1)+ix) 1 +ui(ii,1)*tik(1,ix)+ui(ii,2)*tik(2,ix)+ui(ii,3)*tik(3,ix) endif endif 9 continue endif 8 continue c set vibrational transition dipole moments for each amide do 4 kk=1,ia c setup local axes x ~ CN, z ~ CN x CO, y ~ z x x call asx(rk,r,coa(kk,1)) do 3 ix=1,3 x( ix)=r(ix+3*(coa(kk,4)-1))-rk(ix) 3 co(ix)=r(ix+3*(coa(kk,3)-1))-rk(ix) call nv(x) call vp(x ,co,z) call nv(z) call vp(z , x,y) c transitional vibrational dipoles: do 4 I=1,4 do 4 ix=1,3 4 u(kk,I,ix)=ux(I)*x(ix)+uy(I)*y(ix) pmax=0.0d0 ipmax=0 jpmax=0 c loop over single excited states I do 2 I=1,N if(W(I)*CM.le.wmax.and.W(I)*CM.ge.wmin)then c loop over single excited states J do 201 J=I+1,N if(W(J)*CM.le.wmax.and.W(J)*CM.ge.wmin)then c effective frequency wij=dsqrt(W(I)*W(J)) c loop over "polarizabilities"=amide groups do 7 kk=1,ia c set effective polarizability k for this frequency: call seta(wij,ia,is,kk,u,v,alr,ali,gau) do 7 ix=1,3 do 7 jx=1,3 tt=2.0d0*tau(I,ix+3*(kk-1))*tau(J,jx+3*(kk-1)) if(ali(ix,jx)**2+alr(ix,jx)**2.gt.pmax)then pmax=ali(ix,jx)**2+alr(ix,jx)**2 alimax=ali alrmax=alr kpmax=kk ipmax=I jpmax=J endif Di(I,J)=Di(I,J)-tt*ali(ix,jx) 7 Dr(I,J)=Dr(I,J)-tt*alr(ix,jx) Di(J,I)=-Di(I,J) Dr(J,I)= Dr(I,J) endif 201 continue endif 2 continue write(6,*) call asx(ri,r,coa(imax,1)) call asx(rk,r,coa(kmax,1)) write(6,607)tmax,imax,kmax, 1dsqrt((ri(1)-rk(1))**2+(ri(2)-rk(2))**2+(ri(3)-rk(3))**2) 1*0.529177d0,umax,iumax,kumax,pmax,Ipmax,Jpmax,kpmax, 1((alrmax(ix,jx),ix=1,jx),jx=1,3),((alimax(ix,jx),ix=1,jx),jx=1,3) 607 format(' tmax kmax imax r :',e12.4,2i5,f10.2,'A',/, 1 ' umax Imax Kmax :',e12.4,2i5,/, 1 ' pmax Imax Jmax kmax:',e12.4,3i5,/, 16f10.2,/,6f10.2) return end subroutine dtik(tk,ri,rk,epsr,dik) implicit none integer*4 ix,jx real*8 tk(3,3),ri(3),rk(3),dik,rik(3),epsr do 1 ix=1,3 1 rik(ix)=ri(ix)-rk(ix) dik=dsqrt(rik(1)**2+rik(2)**2+rik(3)**2) do 2 ix=1,3 do 2 jx=1,3 2 tk(ix,jx)=3.0d0*rik(ix)*rik(jx) do 3 ix=1,3 3 tk(ix,ix)=tk(ix,ix)-dik**2 do 4 ix=1,3 do 4 jx=1,3 4 tk(ix,jx)=tk(ix,jx)/dik**5/epsr return end subroutine asx(ri,r,ai) integer*4 ai real*8 ri(3),r(*) ri(1)=r(1+3*(ai-1)) ri(2)=r(2+3*(ai-1)) ri(3)=r(3+3*(ai-1)) return end subroutine seta(w,ia,is,k,u,v,ar,ai,gau) implicit none integer*4 is,ix,jx,iv,k,ia real*8 w,u(ia,4,3),v(4),ar(3,3),ai(3,3),ww,t,gau do 5 ix=1,3 do 5 jx=1,3 ar(ix,jx)=0.0d0 ai(ix,jx)=0.0d0 do 5 iv=is+1,is+2 ww=v(iv)**2-w**2 t=2.0d0*v(iv)*u(k,iv,ix)*u(k,iv,jx)/(ww**2+(v(iv)*gau)**2) ar(ix,jx)=ar(ix,jx) + ww * t 5 ai(ix,jx)=ai(ix,jx) + w*gau * t c wv^2-w^2 + iwG c a_ij = 2 wy ----------------------- ui uj c (wv^2-w^2)^2 + wv^2G^2) return end subroutine vp(a,b,c) real*8 a(*),b(*),c(*) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end subroutine nv(v) real*8 v(*),a,sp a=dsqrt(sp(v,v)) v(1)=v(1)/a v(2)=v(2)/a v(3)=v(3)/a return end function sp(a,b) real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine rdp(N,p) implicit none integer*4 N,i,j real*8 p(2*N,3) open(9,file='FILE.TEN',status='old') read(9,*) do 2 i=1,2*N 2 read(9,*)(p(i,j),j=1,3) close(9) write(6,*)'APT read' return end subroutine rda(ia,coa) implicit none integer*4 ia,coa(ia,6),i,j open(9,file='AM.SCR') do 1 i=1,ia 1 read(9,*)(coa(i,j),j=1,6) close(9) write(6,608)ia 608 format(i6,' amides') return end subroutine ra(ia,r,q,nat) implicit none integer*4 ia,q(*),nat,i,j,io,ic,nn,k,jc(6),icn,ihn real*8 r(*),xi,yi,zi,xj,yj,zj,dij,bohr,xk,yk,zk,dkj bohr=0.529177d0 ia=0 open(9,file='AM.SCR') do 1 i=1,nat xi=r(1+3*(i-1))*bohr yi=r(2+3*(i-1))*bohr zi=r(3+3*(i-1))*bohr do 5 j=1,6 5 jc(j)=0 if(q(i).eq.6)then jc(1)=i io=0 ic=0 nn=0 icn=0 ihn=0 do 2 j=1,nat if(j.ne.i)then xj=r(1+3*(j-1))*bohr yj=r(2+3*(j-1))*bohr zj=r(3+3*(j-1))*bohr dij=dsqrt((xi-xj)**2+(yi-yj)**2+(zi-zj)**2) if(dij.lt.1.65d0.and.q(j).eq.6)then ic=ic+1 jc(2)=j endif if(dij.lt.1.35d0.and.q(j).eq.8)then io=io+1 jc(3)=j endif if(dij.lt.1.55d0.and.q(j).eq.7)then nn=nn+1 jc(4)=j do 3 k=1,nat if(k.ne.i)then xk=r(1+3*(k-1))*bohr yk=r(2+3*(k-1))*bohr zk=r(3+3*(k-1))*bohr dkj=dsqrt((xk-xj)**2+(yk-yj)**2+(zk-zj)**2) if(dkj.lt.1.11d0.and.q(k).eq.1)then ihn=ihn+1 jc(5)=k endif if(dkj.lt.1.55d0.and.q(k).eq.6)then icn=icn+1 jc(6)=k endif endif 3 continue endif endif 2 continue if(io.eq.1.and.ic.eq.1.and.nn.eq.1)then ia=ia+1 write(6,607)ia,(jc(j),j=1,6) 607 format(i6,':',i6,' CO',i6,' C',i6,' O',i6,' N',i6,' H',i6,' CN') c 3O c || c 1C 6C c / \ / c 2C 4N c | c 5H write(9,608)(jc(j),j=1,6) 608 format(6i6) endif endif 1 continue close(9) return end SUBROUTINE FAST(N,F,ZM) IMPLICIT none INTEGER*4 N,I,J real*8 F(N,N),ZM(N),t DO 214 I=1,N t=ZM(I) DO 214 J=I,N F(I,J)=F(I,J)*ZM(J)*t 214 F(J,I)=F(I,J) RETURN END C SUBROUTINE READFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) 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) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) CLOSE(20) WRITE(*,*)' Cartesian FF read in ... ' RETURN END C ============================================================ subroutine readM(nat,ZM) implicit none integer*4 nat,I,nosc real*8 ZM(*),BM,AMU real*8,allocatable::m(:) AMU=1822.88d0 allocate(m(nat)) OPEN(15,FILE='FTRY.INP') READ(15,*) READ(15,*)nosc if(nosc.gt.0)READ(15,*) C Read-in atomic charges and forget: READ(15,*) (m(I),I=1,nat) C Read-in atomic masses: READ(15,*) READ(15,1040)(m(I),I=1,nat) 1040 FORMAT(6G12.6) BM=0.0d0 DO 102 I=1,nat BM=BM+m(I) m(I)=m(I)*AMU ZM(3*I-2)=1.0d0/SQRT(m(I)) ZM(3*I-1)=1.0d0/SQRT(m(I)) 102 ZM(3*I )=1.0d0/SQRT(m(I)) CLOSE(15) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS',F12.2) return end C ============================================================ CD subroutine report(s) character*(*) s write(6,*)s stop end subroutine rx(ic,s,nat,r,q) implicit none integer*4 ic,nat,i,ix,q(*) real*8 r(*),bohr character*(*) s bohr=0.529177d0 open(9,file=s,status='old') read(9,*) read(9,*)nat if(ic.eq.1)then do 1 i=1,nat 1 read(9,*)q(i),(r(ix+3*(i-1)),ix=1,3) do 3 i=1,3*nat 3 r(i)=r(i)/bohr endif close(9) return end c SUBROUTINE SMATOUT(S,e,N,r,q) Cc dX = S dQ S(nat3xnint) c IMPLICIT none c INTEGER*4 N,I,J,q(*) c real*8 S(N,N),e(N),X,Y,Z,r(*),bohr c bohr=0.529177d0 c OPEN(34,FILE='F.INP') c WRITE(34,10)N,N,N/3 c10 FORMAT(3I7) c DO 3 I=1,N/3 c x=r(1+3*(I-1))*bohr c y=r(2+3*(I-1))*bohr c z=r(3+3*(I-1))*bohr c3 WRITE(34,11)q(I),X,Y,Z c11 FORMAT(I7,3F12.6) c WRITE(34,14) c14 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') c DO 1 I=1,N/3 c DO 1 J=1,N c1 WRITE(34,15)I,J,S(3*(I-1)+1,J),S(3*(I-1)+2,J),S(3*(I-1)+3,J) c15 FORMAT(2I7,3F11.6) c WRITE(34,11)N c do 4 I=1,N c WRITE(34,13)e(I) c13 format(F11.3,$) c4 if(mod(I,6).eq.0)write(34,*) c write(34,*) c CLOSE(34) c RETURN c END SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) 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 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) 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 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (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 :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 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 E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) 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) .LE. 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 IF (ICODE .NE. 2) GO TO 300 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 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) c integer i,j,n,nm,ierr,matz double precision ar(n,n),ai(n,n),w(n),zr(n,n),zi(n,n) real*8,allocatable:: fv1(:),fv2(:),fm1(:,:) 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 allocate(fv1(n),fv2(n),fm1(2,n)) 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 if(ierr.ne.0)call report('Diagonalization error') 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