PROGRAM DC_PSI c c transition dipole coupling//short wavelengths extension c IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 nsmax,gj REAL*8 alt(9),gt(9),at(27),ali(9),gi(9),ai(27),gc(9),ac(27),XT(3), 1U(3,3),pi,tau,bohr,debye,cm,Q(3,3),ram(4),roa(4),gau,Qim(3,3),t, 1cid,wau,roa1,CLIGHT,raman logical*4 lpsi,lnm,lor,lnorm,lnum,lcor,lttt,lamid,lw0,la3,l12, 1usebuf,lrt,llook,lcid,lah,ldd,lquad,lpol,ppol integer*4,allocatable::NS(:),nsit(:),which(:) real*8,allocatable::al(:,:,:),G(:,:,:),A(:,:,:,:),E(:),cr3(:), 1sr3(:),EN0(:),EVAL(:),UM(:,:),UA(:,:),R(:,:),RT(:,:),COEFI(:,:), 1DM(:,:),MM(:,:),MI(:,:),COEF(:,:),V(:,:),buf(:),bufr(:), 1rb(:,:),ub(:,:),umb(:,:),eb(:),VB(:,:),rs(:,:,:),we(:,:),VI(:,:), 1DMim(:,:),MIim(:,:),MMim(:,:),V0(:,:) logical,allocatable::ldistr(:) character*8 QF(3),RF(4) data QF/'DCQX.TAB','DCQY.TAB','DCQZ.TAB'/ data RF/'ROAX.TAB','ROAY.TAB','ROAZ.TAB','ROAI.TAB'/ pi=4.0d0*atan(1.0d0) bohr=0.529177d0 debye=2.541765d0 cm=219470.0d0 clight=137.03604d0 c c R in debye^2=1e-36 cgs - this is magic c con=2.0d0*pi*bohr*cm*debye**2*1.0d-8 c 2.0d0*pi*bohr*cm* 1.0d-8 = 1/clight con= debye**2 / clight WRITE(6,66666)2018 66666 FORMAT(72(1H-),/, 115X,'TRANSITION DIPOLE COUPLING - FOR ANY LIGHT WAVELENGTH',/, 272(1H-),/,I15, ' version, UOCHB, Prague',/,/, 35X,'INPUT : DC.SC ... dipole definition',/, 35X,' DC.PAR ... options (Optional for uvcd)',/, 35X,'OUTPUT: DC.TAB ... spectrum',/, 55X,' DIP.MOM ... dipole moments (optional)',/, 65X,' EIGEN.VEC ... eigenvectors',/) c c Read DC.SC: NG=1 nsmax=1 nsdim=1 allocate(NS(1),e(1),UA(1,3),UM(1,3),R(1,3),ldistr(1), 1rs(1,1,3),we(1,1),nsit(1),which(1)) N=1 call rdsc(NG,N,XT,NS,E,UA,UM,R,0,ldistr,nsmax,nsdim,rs,we,nsit, 1which) WRITE(6,*)' ',N,' OSCILLATORS ',NG,' GROUPS' WRITE(6,*) c NS number of oscillators on each group c E oscillator energies c cr3,cr3 - something for psi c EN0 buffer for diagonalization c EVAL eigenvalues c UA,UM magnetic/electric dipoles c R positions c RT magnetic dipole corrected positions c DM,MM,MI resultant transition el and internal and dog mag dipoles c C,V coefficients, hamiltonian deallocate(NS,e,UA,UM,R,ldistr,rs,we,nsit,which) nsdim=nsmax allocate(E(N),UM(N,3),UA(N,3),R(N,3),RT(N,3),V(N,N),NS(NG), 1ldistr(NG),rs(NG,nsdim,3),we(NG,nsdim),nsit(NG),which(N),VI(N,N)) call rdsc(NG,N,XT,NS,E,UA,UM,R,1,ldistr,nsmax,nsdim,rs,we,nsit, 1which) c read some options: call rdpar(EPS0,iwre,NI,lpsi,lnm,lor,lnorm,lnum,lcor,lttt,lamid, 1lw0,la3,usebuf,lrt,NT,NPER,ir,tau,tx,ty,tz,llook,lcid,lah, 1dsf,l12,lquad,lpol,gau,ppol) C SC INPUT FILE USES NANOMETERS, E in eV,debyes c C WORKING UNITS : ATOMIC UNITS !!!!! c c OUTPUT UNITS : cm-1 or nm, debyes c c Rewrite input and put everything into aus: c ============================================================ DO 4 I=1,N c eV -> hartrees: E(I)=E(I)*8065.54093D0/cm DO 4 K=1,3 c nm -> bohrs: R(I,K)=R(I,K)*10.0D0/bohr c debyes -> aus, dsf:scale factor: UA(I,K)=UA(I,K)/debye*dsf 4 UM(I,K)=UM(I,K)/debye tx=tx/bohr ty=ty/bohr tz=tz/bohr ldd=.false. do 6 I=1,NG if(ldistr(I))ldd=.true. if(ldistr(I))then do 7 K=1,nsit(I) do 7 ix=1,3 7 rs(I,K,ix)=10.0d0*rs(I,K,ix)/bohr endif 6 continue c ============================================================ DO 214 I = 1,N xi=R(I,1) yi=R(I,2) zi=R(I,3) if(lcor)then c correct position for magnetic dipole call correct(xi,yi,zi,UA(I,1),UA(I,2),UA(I,3), 1 UM(I,1),UM(I,2),UM(I,3),E(I)) RT(I,1)=xi RT(I,2)=yi RT(I,3)=zi else RT(I,1)=R(I,1) RT(I,2)=R(I,2) RT(I,3)=R(I,3) endif 214 continue c c setup potential matrix call setpg(NG,N,V,E,UA,RT,EPS0,which,ldistr,nsdim,nsit,we,rs, 1lpol,gau,VI) c if ppol, safe dipole-dipole interactions for later if(ppol)then allocate(V0(N,N)) V0=V endif R=RT if(iwre.ne.0)then open(56,file='V.TXT') do 101 I=1,N write(56,609)I 609 format(' Vector ',i5,':') 101 write(56,619)(V(I,J)*cm,J=1,I) 619 format(10f8.0) close(56) endif c correction to amide-amide covalent bonding, AI and AII together if(l12)call vcorrect12(V,N,N) c correction to amide-amide covalent bonding: if(lamid)call vcorrect(V,N,N) c correction to amide-amide interraction through H-bond: if(lah)call hcorrect(V,N,N) c correcting amide to fi1 psi1 fi2 psi2 if(lw0)call wcorrect(V,N,N) c correcting 1-3 amide interaction to fi1 psi1 fi2 psi2 if(la3)call tcorrect(V,N,N) c extending the system by rotation/translation: if(lrt)then M=NT*N write(6,6005)NT,NPER,M,tx*bohr,ty*bohr,tz*bohr,ir,tau 6005 format(' Extending ',i3,' times, periodic unit of ',i5,/, 1 ' New dimension: ',i5,/, 1 ' Translational vector (A): ',3f10.2,/, 1 ' Rotational axis: ',i1,/, 1 ' Angle (deg): ',f10.2) if((N/NPER)*NPER.ne.N)call report('mod N Nper <> 0') c rotational martix: allocate(rb(M,3),umb(M,3),ub(M,3),eb(M)) do 1 ii=1,NT call srm(U,tau*dble(ii-1),ir) call rotv(R,rb,N,M,ii,U,tx,ty,tz) c dipoles: rotate only: call rotv(UA,ub,N,M,ii,U,0.0d0,0.0d0,0.0d0) call rotv(UM,umb,N,M,ii,U,0.0d0,0.0d0,0.0d0) do 1 i=1,N 1 eb(i+N*(ii-1))=E(i) open(8,file='ext.x') write(8,*)'extended system' write(8,*)2*M do 2 i=1,M ubn=dsqrt(ub(i,1)**2+ub(i,2)**2+ub(i,3)**2) x=rb(i,1)*bohr+ub(i,1)/ubn y=rb(i,2)*bohr+ub(i,2)/ubn z=rb(i,3)*bohr+ub(i,3)/ubn write(8,80)6,(rb(i,ix)*bohr,ix=1,3) 2 write(8,80)8,x,y,z 80 format(i5,3f12.6,' 0 0 0 0 0 0 0 0.0') close(8) write(6,*)'final geometry written into ext.x' c setup potential matrix, dipole-dipole interaction: deallocate(E) allocate(VB(M,M),E(M)) call setp(M,VB,eb,ub,rb,rb,EPS0) write(6,*)'new potential set, dipole-dipole' c propagating small V by NPER increments: do 3 ishift=0,M/NPER-1,NPER do 3 i=1,N do 3 j=1,N 3 VB(i+ishift,j+ishift)=V(i,j) write(6,*)'old potential implanted ',M/NPER,' times' deallocate(V,R,UA,UM) allocate(V(M,M),R(M,3),UA(M,3),UM(M,3)) V=VB R=rb E=eb UA=ub UM=umb N=M deallocate(VB,eb,rb,ub,umb) endif c c manual inspection of Hamiltonian: if(llook)call lookatme(V,N) allocate(MI(N,3),DM(N,3),cr3(N),sr3(N),EN0(N),EVAL(N), 1MM(N,3),COEF(N,N),COEFI(N,N),MMim(N,3),DMim(N,3),MIim(N,3)) COEFI=0.0d0 IERR=0 ICODE=2 WRITE(6,*)' DIAGONALIZING THE HAMILTONIAN' if(lpol)then call ch(N,V,VI,EVAL,1,COEF,COEFI) else CALL TRED12(N,V,COEF,EVAL,ICODE,IERR,EN0) IF(IERR.NE.0)call report('DIAGONALIZATION ERROR') endif write(6,6767)EVAL(1)*cm,EVAL(N)*cm 6767 format(f12.2,' .. ',F12.2) if(ppol)call redial(N,EVAL,V0,COEF,COEFI,E,gau) c NN=1 if(lrt)then IF(lnorm)NN=NG*NT if(lpsi)then OPEN(23,FILE='DCE.PSI.TAB') else OPEN(23,FILE='DCE.TAB') endif else IF(lnorm)NN=NG if(lpsi)then OPEN(23,FILE='DC.PSI.TAB') else OPEN(23,FILE='DC.TAB') endif endif Z=1.0d0/dble(NN) if(lquad)then do 10 i=1,3 OPEN( 23+i,FILE=QF(i)) 10 WRITE(23+i,2323) endif OPEN(43,FILE='COS.LST') WRITE(23,2323) 2323 FORMAT( 1 ' MODE D RTOT RINT', 2/,' E DEBYE^2 debye^2=1e-36 cgs (esu**2 cm**2)', 3/,80(1H-)) DTOT=0.0d0 RTOT=0.0d0 dth=pi/dble(NI) dpsi=2.0d0*pi/dble(NI) c psipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsi if(usebuf.and.lpsi)then ibuff=N*(N-1)/2 allocate(buf(ibuff),bufr(ibuff)) if(.not.lor)then write(6,*)'rotation averaging: ' do 801 ib=1,ibuff bufr(ib)=0.0d0 801 buf(ib)=0.0d0 cmax=0.0d0 if(lnum)then write(6,*)' numerical' c preintegrate average: th=-dth/2.0d0 do 285 ith=1,NI th=th+dth st=sin(th) ct=cos(th) psi=-dpsi/2.0d0 do 285 ipsi=1,NI psi=psi+dpsi sp=sin(psi) cp=cos(psi) XxXx = pi*(cp**2+ct**2*sp**2) XyXx = pi*(cp*sp-ct**2*sp*cp) XzXx =-pi*st*ct*sp XyXy = pi*(sp**2+ct**2*cp**2) XzXy = pi*st*cp*ct XzXz = pi*st**2 YxXx=0.0d0 YxXy=-cp**2*ct*pi-ct*sp**2*pi YxXz=-st*cp*pi YyXx= cp**2*ct*pi+ct*sp**2*pi YyXy=0.0d0 YyXz=-st*sp*pi YzXx= st*cp*pi YzXy= st*sp*pi YzXz=0.0d0 Zx=sp*st Zy=-cp*st Zz=ct do 86 k=1,N c 2pi/lambda in bohr^-1: EIb=2.0d0*pi*E(K)*cm*1.0d-8*bohr r3k=(Zx*R(k,1)+Zy*R(k,2)+Zz*R(k,3))*EIb if(abs(r3k).gt.cmax)cmax=abs(r3k) cr3(k)=cos(r3k) 86 sr3(k)=sin(r3k) ib=0 do 285 k=1,N ukx=UA(k,1)*st uky=UA(k,2)*st ukz=UA(k,3)*st ck=cr3(k) sk=sr3(k) do 285 kp=k+1,N ib=ib+1 upx=UA(kp,1) upy=UA(kp,2) upz=UA(kp,3) buf(ib) =buf(ib) +6.0d0*(ck*cr3(kp)+sk*sr3(kp))*( 1 ukx*(upx*XxXx+upy*XyXx+upz*XzXx)+ 1 uky*(upx*XyXx+upy*XyXy+upz*XzXy)+ 1 ukz*(upx*XzXx+upy*XzXy+upz*XzXz)) 285 bufr(ib)=bufr(ib)+6.0d0*(sk*cr3(kp)-ck*sr3(kp))*( 1 ukx*(upx*YxXx+upy*YxXy+upz*YxXz)+ 1 uky*(upx*YyXx+upy*YyXy+upz*YyXz)+ 1 ukz*(upx*YzXx+upy*YzXy+upz*YzXz)) else write(6,*)' analytical' ib=0 do 1862 k=1,N ux=UA(k,1) uy=UA(k,2) uz=UA(k,3) do 1862 kp=k+1,N ib=ib+1 EIb=2.0d0*pi*(E(K)+E(kp))/2.0d0*cm*1.0d-8*bohr wx=UA(kp,1) wy=UA(kp,2) wz=UA(kp,3) xij=R(k,1)-R(kp,1) yij=R(k,2)-R(kp,2) zij=R(k,3)-R(kp,3) r2=xij**2+yij**2+zij**2 if(r2.gt.1.0d-3)then rij=dsqrt(r2) rkkp=rij*EIb if(abs(rkkp).gt.cmax)cmax=abs(rkkp) uu=ux*wx +uy*wy +uz*wz ur=ux*xij+uy*yij+uz*zij wr=wx*xij+wy*yij+wz*zij buf(ib)=3.0d0*((uu-ur*wr/r2)*bessel(0,rkkp) 1 -(uu-3.0d0*ur*wr/r2)*bessel(1,rkkp)/rkkp) bufr(ib)=1.5d0*(wx*(uy*zij-uz*yij)+ 1 wy*(uz*xij-ux*zij)+wz*(ux*yij-uy*xij))*bessel(1,rkkp)/rij else buf(ib)=2.0d0*(ux*wx+uy*wy+uz*wz) endif 1862 continue endif write(6,*)' maximal gomiometric argument: ',cmax endif endif c psipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsi DO 8 I=1,N cosf=0.0d0 c c psipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsi if(lpsi)THEN pf=1.0d0/(2.0d0*pi*cm*1.0d-8*bohr) D=0.0d0 S=0.0d0 if(lnum)then c psi, numerical formulae if(lor)then c fixed orientation do 862 k=1,N c 2pi/lambda in bohr^-1: EIb=2.0d0*pi*E(K)*cm*1.0d-8*bohr r3k=R(k,3)*EIb cr3(k)=cos(r3k) 862 sr3(k)=sin(r3k) do 87 k=1,N cik=COEF(k,I) ukx=UA(k,1) uky=UA(k,2) d =d +1.5d0*cik**2*(ukx*ukx+uky*uky) do 87 kp=k+1,N cp=COEF(kp,I)*cik upx=UA(kp,1) upy=UA(kp,2) d=d+3.0d0*cp*(ukx*upx+uky*upy)*(cr3(k)*cr3(kp)+sr3(k)*sr3(kp)) 87 S=S+1.5d0*cp*(ukx*upy-uky*upx)*(sr3(k)*cr3(kp)-cr3(k)*sr3(kp)) S=S*pf RI=0.0d0 else c rotation average if(.not.usebuf)then c integration not be stored in buffer th=-dth/2.0d0 do 85 ith=1,NI th=th+dth psi=-dpsi/2.0d0 st=sin(th) ct=cos(th) do 85 ipsi=1,NI psi=psi+dpsi sp=sin(psi) cp=cos(psi) XxXx = pi*(cp**2+ct**2*sp**2) XyXx = pi*(cp*sp-ct**2*sp*cp) XzXx =-pi*st*ct*sp XyXy = pi*(sp**2+ct**2*cp**2) XzXy = pi*st*cp*ct XzXz = pi*st**2 YxXx=0.0d0 YxXy=-cp**2*ct*pi-ct*sp**2*pi YxXz=-st*cp*pi YyXx= cp**2*ct*pi+ct*sp**2*pi YyXy=0.0d0 YyXz=-st*sp*pi YzXx= st*cp*pi YzXy= st*sp*pi YzXz=0.0d0 Zx=sp*st Zy=-cp*st Zz=ct do 861 k=1,N c 2pi/lambda in bohr^-1: EIb=2.0d0*pi*E(K)*cm*1.0d-8*bohr r3k=(Zx*R(k,1)+Zy*R(k,2)+Zz*R(k,3))*EIb cr3(k)=cos(r3k) 861 sr3(k)=sin(r3k) do 85 k=1,N cik=COEF(k,I)*st ukx=UA(k,1) uky=UA(k,2) ukz=UA(k,3) do 85 kp=k+1,N cp=COEF(kp,I)*cik upx=UA(kp,1) upy=UA(kp,2) upz=UA(kp,3) d =d +6.0d0*cp*( 1 ukx*(upx*XxXx+upy*XyXx+upz*XzXx)+ 1 uky*(upx*XyXx+upy*XyXy+upz*XzXy)+ 1 ukz*(upx*XzXx+upy*XzXy+upz*XzXz))* 2 (cr3(k)*cr3(kp)+sr3(k)*sr3(kp)) 85 S=S+6.0d0*cp*( 1 ukx*(upx*YxXx+upy*YxXy+upz*YxXz)+ 1 uky*(upx*YyXx+upy*YyXy+upz*YyXz)+ 1 ukz*(upx*YzXx+upy*YzXy+upz*YzXz))* 2 (sr3(k)*cr3(kp)-cr3(k)*sr3(kp)) else c integration on the fly ib=0 do 185 k=1,N cik=COEF(k,I) do 185 kp=k+1,N ib=ib+1 cp=COEF(kp,I)*cik d =d +cp*buf(ib) 185 S=S+cp*bufr(ib) endif D=d*dth*dpsi/8.0d0/pi/pi S=S*dth*dpsi/8.0d0/pi/pi*pf/(-2.0d0) RI=0.0d0 endif else c psi, analytical formulae if(.not.usebuf)then c psi, analytical formulae - on the fly do 186 k=1,N cik=COEF(k,I) ux=UA(k,1) uy=UA(k,2) uz=UA(k,3) do 186 kp=k+1,N wx=UA(kp,1) wy=UA(kp,2) wz=UA(kp,3) cp=COEF(kp,I)*cik xij=R(k,1)-R(kp,1) yij=R(k,2)-R(kp,2) zij=R(k,3)-R(kp,3) r2=xij**2+yij**2+zij**2 if(r2.gt.1.0d-3)then EIb=2.0d0*pi*(E(K)+E(kp))/2.0d0*cm*1.0d-8*bohr rij=dsqrt(r2) rkkp=rij*EIb uu=ux*wx+uy*wy+uz*wz ur=ux*xij+uy*yij+uz*zij wr=wx*xij+wy*yij+wz*zij D=D+3.0d0*cp*((uu-ur*wr/r2)*bessel(0,rkkp) 1 -(uu-3.0d0*ur*wr/r2)*bessel(1,rkkp)/rkkp) S=S+1.5d0*cp*(wx*(uy*zij-uz*yij)+ 1 wy*(uz*xij-ux*zij)+wz*(ux*yij-uy*xij))*bessel(1,rkkp)/rij else D=D+2.0d0*(ux*wx+uy*wy+uz*wz) endif 186 continue else c psi, analytical formulae - in memory ib=0 do 1861 k=1,N cik=COEF(k,I) do 1861 kp=k+1,N ib=ib+1 D=D+COEF(kp,I)*cik*buf(ib) 1861 S=S+COEF(kp,I)*cik*bufr(ib) endif S=S*pf/(-1.0d0) RI=0.0d0 endif if(.not.lor)then do 853 k=1,N 853 D=D+COEF(k,I)**2*(UA(k,1)**2+UA(k,2)**2+UA(k,3)**2) endif c psipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsipsi ELSE c dipolar approximation do 81 ix=1,3 MI( I,ix)=0.0d0 MIim(I,ix)=0.0d0 DM( I,ix)=0.0d0 81 DMim(I,ix)=0.0d0 ax=0.0d0 ay=0.0d0 az=0.0d0 aix=0.0d0 aiy=0.0d0 aiz=0.0d0 if(lquad)then Q=0.0d0 Qim=0.0d0 endif DO 370 j = 1,N cij=COEF(j,I) cijim=COEFI(j,I) if(ldd)then c distributed dipoles gj=which(j) do 9 jj=1,nsit(gj) rx=rs(gj,jj,1) ry=rs(gj,jj,2) rz=rs(gj,jj,3) t=0.5d0*E(j)*(ry*ua(j,3)-rz*ua(j,2))*we(gj,jj) ax=ax+cij*t aix=aix+cijim*t t=0.5d0*E(j)*(rz*ua(j,1)-rx*ua(j,3))*we(gj,jj) ay=ay+cij*y aiy=aiy+cijim*t t=0.5d0*E(j)*(rx*ua(j,2)-ry*ua(j,1))*we(gj,jj) az=az+cij*t 9 aiz=aiz+cijim*t else c Magnetic moment = (1/5) Ej R x u t=0.5d0*E(j)*(R(j,2)*ua(j,3)-R(j,3)*ua(j,2)) ax=ax+cij*t aix=aix+cijim*t t=0.5d0*E(j)*(R(j,3)*ua(j,1)-R(j,1)*ua(j,3)) ay=ay+cij*t aiy=aiy+cijim*t t=0.5d0*E(j)*(R(j,1)*ua(j,2)-R(j,2)*ua(j,1)) az=az+cij*t aiz=aiz+cijim*t if(lquad)then c Quadrupole moment = Ru + uR do 12 ix=1,3 do 12 jx=1,3 t=(R(j,ix)*ua(j,jx)+R(j,jx)*ua(j,ix))/2.0d0 Q( ix,jx)=Q( ix,jx)+cij *t 12 Qim(ix,jx)=Qim(ix,jx)+cijim*t endif endif DO 370 ix=1,3 MI( I,ix)=MI( I,ix)+cij *UM(j,ix) MIim(I,ix)=MIim(I,ix)+cijim*UM(j,ix) DM( I,ix)=DM( I,ix)+cij *UA(j,ix) 370 DMim(I,ix)=DMim(I,ix)+cijim*UA(j,ix) MM(I,1)=ax MM(I,2)=ay MM(I,3)=az MMim(I,1)=aix MMim(I,2)=aiy MMim(I,3)=aiz c real parts: c dipole strength <0|u|1><1|u|0>: D =DM( I,1)* DM(I,1)+DM( I,2)* DM(I,2)+DM( I,3)* DM(I,3) 1 +DMim(I,1)*DMim(I,1)+DMim(I,2)*DMim(I,2)+DMim(I,3)*DMim(I,3) c internal rotational strength <0|u|1><1|m|0>: RI=DM( I,1)* MI(I,1)+DM( I,2)* MI(I,2)+DM( I,3)* MI(I,3) 1 +DMim(I,1)*MIim(I,1)+DMim(I,2)*MIim(I,2)+DMim(I,3)*MIim(I,3) c <0|m|1><1|m|0>: AMR=MM( I,1)* MM(I,1)+MM( I,2)* MM(I,2)+MM( I,3)* MM(I,3) 1 +MMim(I,1)*MMim(I,1)+MMim(I,2)*MMim(I,2)+MMim(I,3)*MMim(I,3) c rotational strength from dipoles <0|u|1><1|m|0>: S =DM( I,1)* MM(I,1)+DM( I,2)* MM(I,2)+DM( I,3)*MM( I,3) 1 +DMim(I,1)*MMim(I,1)+DMim(I,2)*MMim(I,2)+DMim(I,3)*MMim(I,3) if(dabs(D).gt.0.0d0.and.dabs(AMR).gt.0.0d0)then cosf=S/dsqrt(D*AMR) else cosf=0.0d0 endif ENDIF c D in debyes^2 D=D*debye**2 S=S*con RI=RI*con DTOT=DTOT+D*Z RTOT=RTOT+(S+RI)*Z c EI=EVAL(I)*cm if(lnm)EI=1.0d7/EI if(I.eq.1.or.I.EQ.N)WRITE(6,2324)I,EI,D*Z,Z*(S+RI),Z*RI if(I.eq.2)WRITE(6,*)' ... ' c c cosinus between the electric and magnetic moment WRITE(43,2324)I,EI,cosf c if(lquad)then c wavelength in atomic units aulambda=1.0d0/(EVAL(I)*cm)*1.0d8/bohr c wave vector in atomic units auk=2.0d0*pi/aulambda c x-polarization: c ux^2 + (2) (k ux Qyz -(1/c) ux mx); the (2) is somehow hidden in c along x,y,z: do 431 ix=1,3 iy=ix+1 if(iy.gt.3)iy=1 iz=iy+1 if(iz.gt.3)iz=1 dip= (DM(I,iy)**2 +DM(I,iz)**2) *debye**2 quad=auk*(DM(I,iy)*Q(iz,ix)-DM(I,iz)*Q(iy,ix)) *debye**2 amag=- (DM(I,iy)*MM(I,iy)+DM(I,iz)*MM(I,iz))/clight*debye**2 431 WRITE(23+ix,2324)I,EI,1.5d0*dip*Z,-1.5d0*Z*(quad+amag),Z*amag endif 8 WRITE(23,2324)I,EI,D*Z,Z*(S+RI),Z*RI 2324 FORMAT(I4,F20.6,4G16.6) WRITE(23,2325) 2325 FORMAT(80(1H-)) WRITE (23,450) DTOT,RTOT,NN 450 FORMAT ('DSUM=',G16.6,' ','RSUM=',G16.6,/, 1'SPECTRUM OF ALL SYSTEM DEVIDED BY ',I4) CLOSE(23) if(lquad)then do 13 i=24,26 WRITE(i,2325) CLOSE(i) 13 if(lrt)write(6,*)QF(i-23)//' written' endif CLOSE(43) if(lrt)then if(lpsi)then WRITE(6,*)' DCE.PSI.TAB written' else WRITE(6,*)' DCE.TAB written' endif else if(lpsi)then WRITE(6,*)' DC.PSI.TAB written' else WRITE(6,*)' DC.TAB written' endif endif if(lcid)then write(6,*)'Approximate RROA CID' raman=0.0d0 roa1=0.0d0 do 5 i=1,N uxi=UA(I,1) uyi=UA(I,2) uzi=UA(I,3) xi=R(i,1) yi=R(i,2) zi=R(i,3) do 5 j=1,N uxj=UA(j,1) uyj=UA(j,2) uzj=UA(j,3) xj=R(j,1) yj=R(j,2) zj=R(j,3) x=yj*uzj-zj*uyj y=zj*uxj-xj*uzj z=xj*uyj-yj*uxj raman=raman+(uxi*uxi+uyi*uyi+uzi*uzi)*(uxj*uxj+uyj*uyj+uzj*uzj) 1 +7.0d0*(uxi*uxj+uyi*uyj+uzi*uzj)*(uxi*uxj+uyi*uyj+uzi*uzj) 5 roa1=roa1 +(uxi*uxj+uyi*uyj+uzi*uzj)*(uxi* x +uyi* y +uzi* z ) c 532nm in au: wau=0.085645646d0 CLIGHT=137.053d0 cid=(48.0d0*wau*roa1/CLIGHT)/(6.0d0*raman) write(6,504)cid 504 format(' CID = ',e13.4) endif if(lttt)then write(6,*)'Raman and ROA spectra requested' inquire(file='DC.TT',exist=lttt) if(lttt)then open(33,file='DC.TT') read(33,*)NT if(NG.eq.NT)then allocate(al(N,3,3),G(N,3,3),A(N,3,3,3)) I=0 do 82 IG=1,NG read(33,*) do 82 IS=2,NS(IG) I=I+1 read(33,*) read(33,*) do 821 jx=1,3 821 read(33,*)(al(I,ix,jx ),ix=1,3) read(33,*) do 822 jx=1,3 822 read(33,*)(G( I,ix,jx ),ix=1,3) read(33,*) do 823 kx=1,3 do 823 jx=1,3 823 read(33,*)(A( I,ix,jx,kx),ix=1,3) read(33,*) 82 read(33,*) close(33) write(6,*)'DC.TT read' do 90 i=1,4 OPEN(23+i,FILE=RF(i)) 90 write(23+i,900)N 900 format(' n = ',i4,' Raman ROA spectra'/,/,60(1h-)) DO 83 I=1,N EI=EVAL(I)*cm if(lnm)EI=1.0d7/EI call vz(alt,9) call vz(gt,9) call vz(at,27) c rewrite local polarizabilities: DO 84 j = 1,N do 88 ix=1,3 do 88 jx=1,3 ali(ix+3*(jx-1))=al(j,ix,jx) Gi( ix+3*(jx-1))=G( j,ix,jx) do 88 kx=1,3 88 Ai( ix+3*(jx-1)+9*(kx-1))=A( j,ix,jx,kx) c make non-local part of G and A: call mkga(R(j,1),R(j,2),R(j,3),ali,gc,ac) c add to the local part: do 89 ix=1,9 Gi(ix)=Gi(ix)+gc(ix) Ai(ix )=Ai(ix )+ac(ix ) Ai(ix+ 9)=Ai(ix+ 9)+ac(ix+ 9) 89 Ai(ix+18)=Ai(ix+18)+ac(ix+18) c add all to total polarizability: 84 call suma(alt,gt,at,ali,gi,ai,COEF(j,I)) c get the intensity and write it: call doram(ram,roa,alt,gt,at) do 83 io=1,4 83 write(23+io,601)I,EI,ram(io),roa(io),roa(io) 601 format(i4,f10.4,' 1.0 2.0 ',g14.6,' 4.0 5.0 ',g14.6,' 7.0 ', 1 g14.6) do 91 i=1,4 write(23+i,901) 901 format(60(1h-)) write(6,*)RF(i)//' written' 91 close(23+i) else write(6,*)'DC.TT found' call report('NT <> NG') endif else call report('DC.TT not found') endif endif if(iwre.eq.1)CALL WREI(N,COEF,DM,MM,MI,EVAL,R,UA) WRITE(6,*)' >>> PROGRAM TERMINATED <<<' STOP END subroutine redial(N,EVAL,V0,COEF,COEFI,E,gau) implicit none integer*4 N,i,ii,k,J real*8 EVAL(N),V0(N,N),E(N),gau,w,wk,dw,down,cm, 1COEF(N,N),COEFI(N,N) real*8,allocatable::VR(:,:),VI(:,:),CIk(:,:),cnr(:,:),cni(:,:), 1ctr(:,:),cti(:,:) cm=219470.0d0 write(6,*)'Post-diagonalization polarization correction' allocate(VR(N,N),VI(N,N),CIk(N,N),cnr(N,N),cni(N,N), 1ctr(N,N),cti(N,N)) VR=0.0d0 VI=0.0d0 CIk=0.0d0 do 1 I=1,N VR(I,I)=EVAL(I) do 1 k=1,N do 11 ii=1,N 11 CIk(I,k)=CIk(I,k)+COEF(ii,I)*V0(ii,k) 1 CIk(I,k)=CIk(I,k)-COEF( k,I)*V0( k,k) do 2 I=1,N do 2 J=I,N w=dsqrt(EVAL(I)*EVAL(J)) do 22 k=1,N wk=E(k) dw=w**2-wk**2 down=dw**2+(gau*w)**2 VI(I,J)=VI(I,J)+2.0d0*w*w*gau/down*CIk(I,k)*CIK(J,k) 22 VR(I,J)=VR(I,J)+2.0d0*w*dw /down*CIk(I,k)*CIK(J,k) VR(J,I)= VR(I,J) 2 VI(J,I)=-VI(I,J) c diagonalize, replace old energies and coefficients: write(6,*)'Re-diagonalizing' call ch(N,VR,VI,EVAL,1,cnr,cni) do 3 I=1,N do 3 ii=1,N ctr(ii,I)=0.0d0 cti(ii,I)=0.0d0 do 3 J=1,N ctr(ii,I)=ctr(ii,I)+cnr(J,I)*COEF(ii,J) 3 cti(ii,I)=cti(ii,I)+cni(J,I)*COEF(ii,J) COEF=ctr COEFI=cti write(6,6767)EVAL(1)*cm,EVAL(N)*cm 6767 format(f12.2,' .. ',F12.2) return end SUBROUTINE WREI(N,COEF,DM,AMM,AMMI,EVAL,X,UA) IMPLICIT REAL*8(A-H,O-Z) implicit integer*4(i-n) DIMENSION COEF(N,N),DM(N,3),EVAL(N),AMM(N,3),AMMI(N,3), 1X(N,3),UA(N,3) real*8,allocatable::R(:,:) REAL*8 ME(3),MM(3),AI(3,3),TI(3) OPEN(26,FILE='EIGEN.VEC') WRITE(26,*)' THE NORMALIZED EIGENVECTORS (BASE MODE):' M=8 MI=-7 207 MI=MI+8 MJ=M-7 IF(M.GT.N)M=N WRITE(26,28)(I,I=MJ,M) 28 FORMAT(8X,8(I7,2X)) WRITE(26,29)(EVAL(I)*219474.0d0,I=MJ,M) 29 FORMAT(8X,8(f9.2)) DO 206 J=1,N 206 WRITE(26,27)J,(COEF(J,I),I=MJ,M) 27 FORMAT(I5,3X,8F9.5) IF (M.LT.N)THEN M=M+8 GOTO 207 ENDIF CLOSE(26) write(6,*)'EIGEN.VEC written' open(26,file='DIPOLES.TXT') write(26,*)N,' dipoles and positions' do 1 I=1,N 1 write(26,2626)I,(UA(I,ix),ix=1,3),(X(I,ix),ix=1,3) 2626 format(i6,3f10.4,3f12.3) close(26) write(6,*)'DIPOLES.TXT written' DO 281 I=1,N DO 281 J=1,3 281 AMM(I,J)=AMM(I,J)+AMMI(I,J) allocate(R(3,N)) do 2 IA=1,N do 2 ix=1,3 2 R(ix,IA)=X(IA,ix) CALL INERTIA(R,AI,TI,N) OPEN(26,FILE='DIP.MOM') WRITE(26,*)' PROJECTIONS OF TRANSITION DIPOLES:' WRITE(26,2627)(TI(I),I=1,3) 2627 FORMAT(3F15.5,/, 1'MODE DX DY DZ [%] / MX MY MZ [%]') DO 2280 I = 1,N EI=EVAL(I)*219470.0d0 DO 39 II=1,3 ME(II)=0.0d0 MM(II)=0.0d0 DO 39 J=1,3 MM(II)=MM(II)+AI(J,II)*AMM(I,J) 39 ME(II)=ME(II)+AI(J,II)*DM (I,J) ANORM=ME(1)**2+ME(2)**2+ME(3)**2 BNORM=MM(1)**2+MM(2)**2+MM(3)**2 DO 404 IX=1,3 IF(BNORM.GT.1.0D-38)MM(IX)=MM(IX)**2/BNORM*100.0d0 404 IF(ANORM.GT.1.0D-38)ME(IX)=ME(IX)**2/ANORM*100.0d0 2280 WRITE (26,1626)I,EI,(ME(IX),IX=1,3),(MM(IX),IX=1,3) 1626 FORMAT(I4,F8.1,3F6.1,' /',3F6.1) CLOSE(26) write(6,*)'DIP.MOM written' c electric dipoles RETURN END SUBROUTINE INERTIA(R,A,TI,N) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(i-n) DIMENSION R(3,N),A(3,3),T(3,3),TI(3),CM(3),EE(3) DO 2 I=1,3 CM(I)=0.0d0 DO 3 J=1,N 3 CM(I)=CM(I)+R(I,J) 2 CM(I)=CM(I)/N DO 1 I=1,3 DO 1 J=1,3 T(I,J)=0.0d0 DO 1 IS=1,N 1 T(I,J)=T(I,J)+(R(I,IS)-CM(I))*(R(J,IS)-CM(J)) CALL TRED(3,3,T,A,TI,2,IERR,EE) RETURN C DEL(I,J)=SUM(OVER JS) OF A(I,JS)*A(J,JS) C DIAGONAL TD(I,J) = A(IS,I)*T(IS,JS)*A(JS,j) END SUBROUTINE TRED(N0,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(i-n) DIMENSION A(N0,N0),Z(N0,N0),D(N0),E(N0) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 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 DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) JP1 = J + 1 IF (L .LT. JP1) GO TO 220 DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE HH = F / (H + H) DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) 290 D(I) = H 300 CONTINUE 320 D(1) = 0.0D0 E(1) = 0.0D0 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 DO 360 J = 1, L G = 0.0D0 DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE 500 CONTINUE 620 CALL TQL (N0,N,Z,D,IERR,ICODE,E) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END SUBROUTINE TQL(N0,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT integer*4(i-n) DIMENSION Z(N0,N0),D(N0),E(N0) 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 IERR = 0 IF (N .EQ. 1) GO TO 1001 DO 100 I = 2, N 100 E(I-1) = E(I) F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 110 CONTINUE 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 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) DO 140 I = L1, N 140 D(I) = D(I) - H F = F + H P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L 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)) 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 200 CONTINUE 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 GO TO 1001 1000 IERR = L 1001 RETURN END subroutine report(s) character*(*) s write(6,60) write(6,*)s write(6,60) 60 format(60(1H*)) stop end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(N),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(N),E(N) 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 ===================================================================== function bessel(n,x) implicit none integer*4 n real*8 x,bessel if(n.eq.0)then bessel=sin(x)/x else if(n.eq.1)then bessel=(sin(x)/x-cos(x))/x else call report('unknown Bessel function') endif endif return end c ===================================================================== c correct position for magnetic dipole subroutine correct(x,y,z,ux,uy,uz,mx,my,mz,e) implicit none real*8 x,y,z,ux,uy,uz,mx,my,mz,e,bx,by,bz,mb,dr,u,ax,ay,az ax=uy*mz-uz*my ay=uz*mx-ux*mz az=ux*my-uy*mx call norm(ax,ay,az) bx=ay*uz-az*uy by=az*ux-ax*uz bz=ax*uy-ay*ux call norm(bx,by,bz) mb=bx*mx+by*my+bz*mz u=dsqrt(ux**2+uy**2+uz**2) if(u.lt.1.0d-5)return if(e.lt.1.0d-15)return dr=-2.0d0*mb/(e*219470.0d0*3.3248d-8)/u x=x+dr*ax y=y+dr*ay z=z+dr*az return end c ===================================================================== subroutine norm(ax,ay,az) implicit none real*8 ax,ay,az,an an=dsqrt(ax**2+ay**2+az**2) if(dabs(an).lt.1.0d-10)return ax=ax/an ay=ay/an az=az/an return end c ===================================================================== subroutine vz(v,n) real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end c ===================================================================== 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 c ===================================================================== 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 c ===================================================================== subroutine doram(ram,roa,al,g,a) implicit none real*8 ram(4),roa(4),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) c xyz directional scattering: do 1 i=1,3 call setjk(i,j,k) ram(i)=180.0d0*0.5d0*AMU*BOHR**4*( 1al(j+3*(j-1))**2+al(k+3*(k-1))**2+2.0d0*al(j+3*(k-1))**2) 1 roa(i)=-180.0d0*gpisvejc*( 1-2.0d0*al(j+3*(k-1))*(g(j+3*(k-1))+g(k+3*(j-1))) 1+(al(k+3*(k-1))-al(j+3*(j-1)))*(g(j+3*(j-1))-g(k+3*(k-1))) 1+2.0d0/3.0d0* 1(al(k+3*(k-1))*a(j+3*(i-1)+9*(k-1)) 1-al(j+3*(j-1))*a(k+3*(i-1)+9*(j-1)) 1+al(k+3*(j-1))*a(j+3*(i-1)+9*(j-1)) 1-al(j+3*(k-1))*a(k+3*(i-1)+9*(k-1)))) c isotropic average: ram(4)=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) return end c ===================================================================== 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 c ===================================================================== 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 c ===================================================================== subroutine rdpar(EPS0,iwre,NI,lpsi,lnm,lor,lnorm,lnum,lcor,lttt, 1lamid,lw0,la3,usebuf,lrt,NT,NPER,ir,tau,tx,ty,tz,llook,lcid,lah, 1dsf,l12,lquad,lpol,gau,ppol) implicit none reaL*8 EPS0,tau,tx,ty,tz,dsf,gau,gcm integer*4 iwre,NI,NT,ir,NPER logical lpsi,lnm,lor,lnorm,lnum,lcor,lttt,lamid,lex,lw0,la3, 1usebuf,lrt,llook,lcid,lah,l12,lquad,lpol,ppol CHARACTER*80 s80,s81 EPS0=1.0d0 gcm=20.0d0 iwre=0 lpsi=.false. lquad=.false. lnm=.false. lor=.false. lnorm=.false. lnum=.false. lcor=.false. lttt=.false. lamid=.false. l12=.false. lw0=.false. la3=.false. lah=.false. usebuf=.true. lrt=.false. llook=.false. lpol=.false. NI=20 NPER=1 tx=0.0d0 ty=0.0d0 tz=0.0d0 tau=0.0d0 lcid=.true. ppol=.false. NT=1 dsf=1.0d0 c NI - number of integration point inquire(file='DC.PAR',exist=lex) if(lex)then open(45,file='DC.PAR') 1011 read(45,4545)s80 4545 format(a80) c realtive permittivity: if(s80(1:3).eq.'EPS')read(45,*)EPS0 c psi-integrate through buffer: if(s80(1:3).eq.'BUF')read(45,*)usebuf if(s80(1:4).eq.'IWRE')read(45,*)iwre c use the general formula: if(s80(1:4).eq.'LPSI')read(45,*)lpsi c extend by rotation and translation: if(s80(1:2).eq.'RT')read(45,*)lrt,NT,NPER,ir,tau,tx,ty,tz if(s80(1:2).eq.'NI')read(45,*)NI if(s80(1:3).eq.'LNM')read(45,*)lnm c use fixed orientation only: if(s80(1:3).eq.'LOR')read(45,*)lor if(s80(1:5).eq.'LNORM')read(45,*)lnorm if(s80(1:4).eq.'LNUM')read(45,*)lnum if(s80(1:4).eq.'LCOR')read(45,*)lcor if(s80(1:4).eq.'LTTT')read(45,*)lttt c manually inspect H: if(s80(1:4).eq.'LOOK')read(45,*)llook c replace amide I - amide II coupling by pre-fitted dependence if(s80(1:3).eq.'L12')read(45,*)l12 c replace amide I - amide I coupling by pre-fitted dependence if(s80(1:4).eq.'LAMI')read(45,*)lamid c replace diagonal amide freq by pre-fitted coupling: if(s80(1:3).eq.'LW0')read(45,*)lw0 c replace 1-3 amide interaction: if(s80(1:3).eq.'LA3')read(45,*)la3 c replace amide interactions across hydrogen bond: if(s80(1:3).eq.'LAH')read(45,*)lah c the approximate model for CID: if(s80(1:4).eq.'LCID')read(45,*)lcid c dipole scale factor" if(s80(1:3).eq.'DSF')read(45,*)dsf c Quadrupolar interaction for oriented samples: if(s80(1:4).eq.'LQUA')read(45,*)lquad c Add the polarizability interaction if(s80(1:4).eq.'LPOL')read(45,*)lpol c post-diagonalization polarizability interaction if(s80(1:4).eq.'PPOL')read(45,*)ppol c bandwidth for polarizability in cm-1: if(s80(1:4).eq.'GAMM')read(45,*)gcm if(s80(1:3).ne.'END')then backspace 45 read(45,4545)s81 write(6,4546)s80(1:40),s81(1:10) 4546 format(a40,' ',a10) goto 1011 endif close(45) WRITE(6,*) if(.not.lnum.and.lor)then write(6,*)'Incompatible options lor without lnum' write(6,*)'Numerical integration switched on' lnum=.true. endif else WRITE(6,*) WRITE(6,*)' Relative permitivity:' READ(5,*)EPS0 endif gau=gcm/219474.0d0 return end c ===================================================================== subroutine rdsc(NG,N,XT,NS,E,UA,UM,R,ic,ldistr,nsmax,nsdim, 1rs,we,nsit,which) c ic = 0 ... read dimensions c 1 ... read variables implicit none real*8 XT(*),OMEGA,UA(N,3),UM(N,3),R(N,3),E(*),rs(NG,nsdim,3), 1we(NG,nsdim) integer*4 nsmax,NG,N,I,NS(*),ii,J,K,ic,nt,nsit(*),nsi, 1nsdim,which(*),IG,ix logical ldistr(*),ltem OPEN(22,FILE='DC.SC',STATUS='OLD') READ(22,*)NG N=0 DO 1 I=1,NG READ(22,*)IG ltem=IG.lt.0 if(ic.ne.0)ldistr(I)=ltem if(ic.ne.0)then READ(22,*)(XT(j),j=1,3) if(ldistr(I))then read(22,*)nsit(I) do 3 j=1,nsit(I) 3 read(22,*)(rs(I,j,ix),ix=1,3),we(I,j) endif READ(22,*)NS(I) nt=NS(I) else READ(22,*) if(ltem)then read(22,*)nsi if(nsi.gt.nsmax)nsmax=nsi do 4 j=1,nsi 4 read(22,*) endif READ(22,*)nt endif READ(22,*) READ(22,*)OMEGA,OMEGA ii=N DO 2 J=2,nt N=N+1 if(ic.ne.0)then which(N)=I READ(22,*)E(N),E(N) E(N)=E(N)-OMEGA else READ(22,*) endif 2 continue N=ii READ(22,*) DO 1 J=2,nt N=N+1 if(ic.ne.0)then READ(22,*)UA(N,1),(UA(N,K),K=1,3),(UM(N,K),K=1,3) DO 10 K=1,3 10 R(N,K)=XT(K) else READ(22,*) endif 1 continue CLOSE(22) return end c ===================================================================== subroutine hcorrect(V,N,N0) implicit none integer*4 N,N0,na,nat,ia,ib,ix,ip,nb,i real*8 V(N0,N0),da,vold,cm,dab real*8,allocatable::r(:) integer*4,allocatable::m(:) real*8,allocatable::c(:) character*80 s80 cm=219470.0d0 WRITE(6,605) 605 format(/,'Correcting amide I - amide I H-bond neighbors') open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat allocate(r(3*nat)) do 1 ia=1,nat 1 read(9,*)r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' open(9,file='FFITH.PAR',status='old') read(9,*) read(9,*)nb if(nb.ne.1)call report('hcorrect for nb<>1 not implemented') allocate(c(nb)) do 4 i=1,nb read(9,909)s80 909 format(a80) 4 read(s80(38:50),*)c(i) close(9) write(6,*)nb,' basis coefficients in FFITH.PAR' open(9,file='AMIDE.LST',status='old') read(9,*)na allocate(m(6*na)) do 2 ia=1,na 2 read(9,*)m(1+6*(ia-1)),(m(ix+6*(ia-1)),ix=1,6) close(9) write(6,*)na,' amides in AMIDE.LST' if(N.ne.na)call report('N <> na') c order C - CO - NH-C amide groups: c 6 21 34 5 c find neighboring amides H-bond connected ip=0 do 3 ia=1,na do 3 ib=1,na dab=da(m(1+6*(ia-1)),m(4+6*(ib-1)),r) if(ia.ne.ib.and.dab.lt.2.20d0)then c 6 21 34 5 c C - C - NH-C ia c \\ c O c . c H c / c C - CO - N -C ib c 6 21 34 5 vold=V(ia,ib) V(ia,ib)=C(1)/cm V(ib,ia)=C(1)/cm ip=ip+1 write(6,609)ip,ia,ib,dab,vold*cm,v(ia,ib)*cm 609 format(3i5,4f12.3) endif 3 continue write(6,*)ip,' elements replaced' return end c ===================================================================== subroutine vcorrect(V,N,N0) implicit none integer*4 N,N0,na,nat,ia,ib,ix,ip,nb,i real*8 V(N0,N0),ta,va,phi,psi,vold,cm real*8,allocatable::r(:) integer*4,allocatable::m(:) logical replace real*8,allocatable::c(:) character*80 s80 cm=219470.0d0 WRITE(6,605) 605 format(/,'Correcting amide I - amide I neighbors') open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat allocate(r(3*nat)) do 1 ia=1,nat 1 read(9,*)r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' open(9,file='FFIT.PAR',status='old') read(9,*) read(9,*)nb allocate(c(nb)) do 4 i=1,nb read(9,909)s80 909 format(a80) 4 read(s80(38:50),*)c(i) close(9) write(6,*)nb,' basis coefficients in FFIT.PAR' open(9,file='AMIDE.LST',status='old') read(9,*)na allocate(m(6*na)) do 2 ia=1,na 2 read(9,*)m(1+6*(ia-1)),(m(ix+6*(ia-1)),ix=1,6) close(9) write(6,*)na,' amides in AMIDE.LST' if(N.ne.na)call report('N <> na') c order C - CO - NH-C amide groups: c 6 21 34 5 c find neighboring amides ip=0 do 3 ia=1,na do 3 ib=ia+1,na replace=.false. if(m(6+6*(ia-1)).eq.m(5+6*(ib-1)))then c phi:2 3 5 2 c psi: 3 5 2 3 c C - CO - NH-(C = C) - CO - NH-C c 6 21 34 (5 = 6) 21 34 5 c ib | ia phi=ta(m(2+6*(ib-1)),m(3+6*(ib-1)),m(5+6*(ib-1)),m(2+6*(ia-1)),r) psi=ta(m(3+6*(ib-1)),m(5+6*(ib-1)),m(2+6*(ia-1)),m(3+6*(ia-1)),r) replace=.true. endif if(m(6+6*(ib-1)).eq.m(5+6*(ia-1)))then c C - CO - NH-(C = C) - CO - NH-C c 6 21 34 (5 = 6) 21 34 5 c ia | ib phi=ta(m(2+6*(ia-1)),m(3+6*(ia-1)),m(5+6*(ia-1)),m(2+6*(ib-1)),r) psi=ta(m(3+6*(ia-1)),m(5+6*(ia-1)),m(2+6*(ib-1)),m(3+6*(ib-1)),r) replace=.true. endif if(replace)then vold=V(ia,ib) V(ia,ib)=va(phi,psi,nb,c)/cm V(ib,ia)=V(ia,ib) ip=ip+1 write(6,609)ip,ia,ib,phi,psi,vold*cm,v(ia,ib)*cm 609 format(3i5,4f12.3) endif 3 continue write(6,*)ip,' elements replaced' return end c ===================================================================== subroutine wcorrect(V,N,N0) implicit none integer*4 N,N0,na,nat,ia,ib,ix,ip,nb,nd,cac,can,ccc,cbn,ic real*8 V(N0,N0),ta,v0,phi1,psi1,phi2,psi2,vold,cm real*8,allocatable::r(:) integer*4,allocatable::m(:) integer*4,allocatable::state(:,:) real*8,allocatable::c(:) cm=219470.0d0 WRITE(6,605) 605 format(/,'Correcting amide I diagonal frequency') open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat allocate(r(3*nat)) do 1 ia=1,nat 1 read(9,*)r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' open(9,file='FFIT4.PAR',status='old') read(9,*)nd if(nd.ne.4)call report('nd <> 4') read(9,*)nb allocate(c(nb),state(nb,nd)) call rdp(nb,nd,state,c) write(6,*)nb,' basis coefficients in FFIT.PAR' open(9,file='AMIDE.LST',status='old') read(9,*)na allocate(m(6*na)) do 2 ia=1,na 2 read(9,*)m(1+6*(ia-1)),(m(ix+6*(ia-1)),ix=1,6) close(9) write(6,*)na,' amides in AMIDE.LST' if(N.ne.na)call report('N <> na') c order C - CO - NH-C amide groups: c 6 21 34 5 c find three neighboring amides ip=0 do 3 ia=1,na cac=m(6+6*(ia-1)) can=m(5+6*(ia-1)) do 3 ib=1,na cbn=m(5+6*(ib-1)) do 3 ic=1,na ccc=m(6+6*(ic-1)) if(cac.eq.cbn.and.can.eq.ccc)then c phi1 :2 3 5 2 c psi1 : 3 5 2 3 c phi2 : 2 3 5 2 c psi2 : 3 5 2 3 c C - CO - NH-(C = C ) - CO - NH-(C = C ) - CO - NH-C c 6 21 34 (5 = 6 ) 21 34 (5 = 6 21 34 5 c ib | ia | ic c cbc cbn cac can ccc cnc phi1=ta(m(2+6*(ib-1)),m(3+6*(ib-1)),cbn,m(2+6*(ia-1)),r) psi1=ta(m(3+6*(ib-1)),cbn,m(2+6*(ia-1)),m(3+6*(ia-1)),r) phi2=ta(m(2+6*(ia-1)),m(3+6*(ia-1)),can,m(2+6*(ic-1)),r) psi2=ta(m(3+6*(ia-1)),can,m(2+6*(ic-1)),m(3+6*(ic-1)),r) vold=V(ia,ia) c deviation from the average from the fit V(ia,ia)=vold+(v0(phi1,psi1,phi2,psi2,nb,c,state,nd)-c(1))/cm ip=ip+1 write(6,609)ip,ia,ib,ic,phi1,psi1,phi2,psi2,vold*cm,v(ia,ia)*cm 609 format(4i5,6f12.3) endif 3 continue write(6,*)ip,' elements replaced' return end c ===================================================================== subroutine rdp(nb,nd,state,c) implicit none integer*4 nb,nd,state(nb,nd),i,j,k,kk real*8 c(*) character*16 s16 do 4 i=1,nb read(9,9091) 9091 format(5x,$) do 41 j=1,nd read(9,9092)s16 9092 format(a16,$) if(s16(1:3).eq.'one')then k=1 else read(s16(5:7),*)kk if(s16(1:3).eq.'sin')then k=kk else k=kk+1 endif endif 41 state(i,j)=k 4 read(9,9096)c(i) 9096 format(f12.4) close(9) return end c ===================================================================== subroutine tcorrect(V,N,N0) implicit none integer*4 N,N0,na,nat,ia,ib,ix,ip,nb,nd,cac,can,ccc,cbn,ic real*8 V(N0,N0),ta,v0,phi1,psi1,phi2,psi2,vold,cm real*8,allocatable::r(:) integer*4,allocatable::m(:) integer*4,allocatable::state(:,:) real*8,allocatable::c(:) cm=219470.0d0 WRITE(6,605) 605 format(/,'Correcting 1-3 amide interactions') open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat allocate(r(3*nat)) do 1 ia=1,nat 1 read(9,*)r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' open(9,file='FFIT4.13.PAR',status='old') read(9,*)nd if(nd.ne.4)call report('nd <> 4') read(9,*)nb allocate(c(nb),state(nb,nd)) call rdp(nb,nd,state,c) write(6,*)nb,' basis coefficients in FFIT4.13.PAR' open(9,file='AMIDE.LST',status='old') read(9,*)na allocate(m(6*na)) do 2 ia=1,na 2 read(9,*)m(1+6*(ia-1)),(m(ix+6*(ia-1)),ix=1,6) close(9) write(6,*)na,' amides in AMIDE.LST' if(N.ne.na)call report('N <> na') c order C - CO - NH-C amide groups: c 6 21 34 5 c find three neighboring amides ip=0 do 3 ia=1,na cac=m(6+6*(ia-1)) can=m(5+6*(ia-1)) do 3 ib=1,na cbn=m(5+6*(ib-1)) do 3 ic=1,na ccc=m(6+6*(ic-1)) if(cac.eq.cbn.and.can.eq.ccc)then c phi1 :2 3 5 2 c psi1 : 3 5 2 3 c phi2 : 2 3 5 2 c psi2 : 3 5 2 3 c C - CO - NH-(C = C ) - CO - NH-(C = C ) - CO - NH-C c 6 21 34 (5 = 6 ) 21 34 (5 = 6 21 34 5 c ib | ia | ic c cbc cbn cac can ccc cnc phi1=ta(m(2+6*(ib-1)),m(3+6*(ib-1)),cbn,m(2+6*(ia-1)),r) psi1=ta(m(3+6*(ib-1)),cbn,m(2+6*(ia-1)),m(3+6*(ia-1)),r) phi2=ta(m(2+6*(ia-1)),m(3+6*(ia-1)),can,m(2+6*(ic-1)),r) psi2=ta(m(3+6*(ia-1)),can,m(2+6*(ic-1)),m(3+6*(ic-1)),r) vold=V(ib,ic) V(ib,ic)=v0(phi1,psi1,phi2,psi2,nb,c,state,nd)/cm V(ic,ib)=V(ib,ic) ip=ip+1 write(6,609)ip,ia,ib,ic,phi1,psi1,phi2,psi2,vold*cm,v(ib,ic)*cm 609 format(4i5,6f12.3) endif 3 continue write(6,*)ip,' elements replaced' return end function ta(l1,l2,l3,l4,r) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) DIMENSION XX(3),YY(3),A(3),P(3),E0(3),B0(3), 1AP(3),r(*) PI=3.141592653589D0 do 23 ix=1,3 A(ix) =r(ix+3*(l3-1))-r(ix+3*(l4-1)) P(ix) =r(ix+3*(l1-1))-r(ix+3*(l2-1)) 23 E0(ix)=r(ix+3*(l2-1))-r(ix+3*(l3-1)) call vp(B0,E0,A) SSB0B0=sp(B0,B0) SSE0E0=sp(E0,E0) DO 8 KK=1,3 B0(KK)=B0(KK)/DSQRT(SSB0B0) 8 E0(KK)=E0(KK)/DSQRT(SSE0E0) SSAE0=sp(A,E0) SSPE0=sp(P,E0) DO 9 KK=1,3 XX(KK)=A(KK)-E0(KK)*SSAE0 9 YY(KK)=P(KK)-E0(KK)*SSPE0 SSXXXX=sp(XX,XX) SSYYYY=sp(YY,YY) SSXXYY=sp(XX,YY) call vp(AP,XX,YY) SSE0AP=sp(AP,E0) SIGN=1.0D0 IF (dabs(SSE0AP).gt.1.0D-9)SIGN=-SSE0AP/ABS(SSE0AP) CA=0.0D0 IF (dabs(SSXXXX).gt.1.0D-9.AND.dabs(SSYYYY).gt.1.0D-9) 1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY) ANGLE=0.0D0 IF (CA.LE.-1.0D0) ANGLE=PI IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA) ANGLE=180.0D0* ANGLE/PI ANGLE=SIGN*ANGLE if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE ta=ANGLE return end function da(l1,l2,r) c distance of two atoms IMPLICIT none INTEGER*4 l1,l2 real*8 r(*),da da=dsqrt( (r(1+3*(l1-1))-r(1+3*(l2-1)))**2+ 1 (r(2+3*(l1-1))-r(2+3*(l2-1)))**2+ 1 (r(3+3*(l1-1))-r(3+3*(l2-1)))**2) return end function sp(a,b) IMPLICIT none real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(A,X,Y) implicit none real*8 A(*),X(*),Y(*) A(1)=X(2)*Y(3)-X(3)*Y(2) A(2)=X(3)*Y(1)-X(1)*Y(3) A(3)=X(1)*Y(2)-X(2)*Y(1) return end function va(x,y,nb,c) implicit none integer*4 nb,nb1,k1,k2,ib real*8 va,x,y,c(*),s,a1,fk a1=sqrt(dble(nb)) nb1=nint(a1) ib=0 s=0.0d0 do 1 k1=1,nb1 do 1 k2=1,nb1 ib=ib+1 1 s=s+c(ib)*fk(k1,x,0.0d0,360.0d0)*fk(k2,y,0.0d0,360.0d0) va=s return end function v0(x1,x2,x3,x4,nb,c,state,nd) implicit none integer*4 nb,nd,state(nb,nd),i real*8 v0,x1,x2,x3,x4,c(*),s,fk s=0.0d0 do 1 i=1,nb 1 s=s+c(i)*fk(state(i,1),x1,0.0d0,360.0d0) 1 *fk(state(i,2),x2,0.0d0,360.0d0) 1 *fk(state(i,3),x3,0.0d0,360.0d0) 1 *fk(state(i,4),x4,0.0d0,360.0d0) v0=s return end function fk(k,x,x0,L) c K: 1 2 3 4 5 6 c fik: 1 sin(2pix/L) cos(2pix/L) sin(4pix/L) cos(4pix/L) sin(6pix/L) implicit none integer*4 k real*8 x,x0,L,y,fk,pi pi=4.0d0*datan(1.0d0) y=x-x0 if(k.eq.1)then fk=1.0d0 else if(mod(k,2).ne.0)then fk=cos(dble(2*((k-1)/2))*pi*y/L) else fk=sin(dble(2*(k/2))*pi*y/L) endif endif return end subroutine setpg(NG,N,V,E,UA,RT,EPS0,which,ldistr,nsdim,nsit, 1we,rs,lpol,gau,VI) implicit none integer*4 N,I,J,which(*),nsit(*),nsdim,NG,gi,gj,ii,jj real*8 V(N,N),E(*),UA(N,3),xi,yi,zi,gau,RT(N,3),uix,uiy,uiz, 1xj,yj,zj,VIJ,EPS0,ujx,ujy,ujz,we(NG,nsdim), 1rs(NG,nsdim,3),vijf,VI(N,N) logical ldistr(*),lpol DO 2101 I = 1,N 2101 V(I,I)=E(I) c DO 210 I = 1,N gi=which(I) xi=RT(I,1) yi=RT(I,2) zi=RT(I,3) uix=UA(I,1) uiy=UA(I,2) uiz=UA(I,3) DO 210 J = I+1,N gj=which(J) xj=RT(J,1) yj=RT(J,2) zj=RT(J,3) ujx=UA(J,1) ujy=UA(J,2) ujz=UA(J,3) VIJ=0.0d0 IF(gi.ne.gj)then if(ldistr(gi).and.ldistr(gj))then do 1 ii=1,nsit(gi) uix=UA(I,1)*we(gi,ii) uiy=UA(I,2)*we(gi,ii) uiz=UA(I,3)*we(gi,ii) xi=rs(gi,ii,1) yi=rs(gi,ii,2) zi=rs(gi,ii,3) do 1 jj=1,nsit(gj) ujx=UA(J,1)*we(gj,jj) ujy=UA(J,2)*we(gj,jj) ujz=UA(J,3)*we(gj,jj) xj=rs(gj,jj,1) yj=rs(gj,jj,2) zj=rs(gj,jj,3) 1 VIJ=VIJ+vijf(xi,yi,zi,xj,yj,zj,uix,uiy,uiz,ujx,ujy,ujz,EPS0) else VIJ=vijf(xi,yi,zi,xj,yj,zj,uix,uiy,uiz,ujx,ujy,ujz,EPS0) endif ENDIF V(I,J)=VIJ 210 V(J,I)=VIJ call minimax(V,N,'Vij:') if(lpol)call addpol(N,V,VI,gau) return end subroutine addpol(N,V,VI,gau) implicit none integer*4 N,I,J,K real*8 V(N,N),VI(N,N),gau do 1 I=1,N VI(I,I)=0.0d0 do 1 J=I+1,N VI(I,J)=0.0d0 do 2 K=1,N 2 if(K.ne.I.and.K.ne.J)VI(I,J)=VI(I,J)-4.0d0/gau*V(I,K)*V(J,K) 1 VI(J,I)=VI(I,J) call minimax(VI,N,'Imaginary Vij:') return 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 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,vimax,ijimax,ijjmax,vijmax 659 format( 'Maximal diagonal element ',i5,5x, 1f12.6,/,'Maximal off-diagonal element ',i5,i5,f12.6) return end subroutine setp(N,V,E,UA,R,RT,EPS0) implicit none integer*4 N,I,J,imax,ijimax,ijjmax real*8 V(N,N),E(*),UA(N,3),R(N,3),xio,yio,zio,xi,yi,zi, 1RT(N,3),uix,uiy,uiz,xjo,yjo,zjo,xj,yj,zj,xji,yji,zji,D2,D2o, 1rij,rijo,VIJ,uij,uir,ujr,EPSR,EPS0,vimax,vijmax,ujx,ujy,ujz,D5 DO 2101 I = 1,N 2101 V(I,I)=E(I) c DO 210 I = 1,N xio=R(I,1) yio=R(I,2) zio=R(I,3) xi=RT(I,1) yi=RT(I,2) zi=RT(I,3) uix=UA(I,1) uiy=UA(I,2) uiz=UA(I,3) DO 210 J = I+1,N xjo=R(J,1) yjo=R(J,2) zjo=R(J,3) xj=RT(J,1) yj=RT(J,2) zj=RT(J,3) ujx=UA(J,1) ujy=UA(J,2) ujz=UA(J,3) xji=xj-xi yji=yj-yi zji=zj-zi D2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 D2o=(xio-xjo)**2+(yio-yjo)**2+(zio-zjo)**2 rij=dsqrt(D2) rijo=dsqrt(D2o) D5=D2*D2*rij c c allow to interact only oscillators in different group (rij>0): VIJ=0.0d0 IF(rijo.GT.1.0d-4)THEN uij=uix*ujx+uiy*ujy+uiz*ujz uir=uix*xji+uiy*yji+uiz*zji ujr=ujx*xji+ujy*yji+ujz*zji EPSR=EPS0 IF (EPS0.LT.0.0D0) EPSR=-EPS0*rij VIJ=(uij*D2-3.0d0*uir*ujr)/D5/EPSR ENDIF V(I,J)=VIJ 210 V(J,I)=VIJ vimax=abs(v(1,1)) imax=1 vijmax=abs(v(1,2)) ijimax=1 ijjmax=2 do i=1,N if(vimax.lt.dabs(v(i,i)))then vimax=abs(v(i,i)) imax=i endif do j=i+1,N if(vijmax.lt.dabs(v(i,j)))then vijmax=dabs(v(i,j)) ijimax=i ijjmax=j endif enddo enddo write(6,659)imax,vimax,ijimax,ijjmax,vijmax 659 format( 'Maximal diagonal element ',i5,5x, 1f12.6,/,'Maximal off-diagonal element ',i5,i5,f12.6) return end subroutine lookatme(V,N) implicit none integer*4 N,I,J real*8 V(N,N),cm cm=219470.0d0 212 WRITE(6,*)' I,J FOR THE H(I,J) {[0,0] TO CONTINUE, [-100,0]', 1',[-200,0] OPTIONS}:' READ(5,*)I,J IF(I.EQ.-100)THEN WRITE(6,*)' 1-2,3-4,5-6 ETC. INTERACTIONS CANCELED' DO 67 I=1,N-1,2 V(I,I+1)=0.0d0 67 V(I+1,I)=0.0d0 GOTO 212 ENDIF IF(I.EQ.-200)THEN WRITE(6,*)' ODD-EVEN AND EVEN-ODD INTERACTIONS CANCELED' DO 68 I=1,N DO 68 J=1,I IF(MOD(I-J,2).NE.0)THEN V(I,J)=0.0d0 V(J,I)=0.0d0 ENDIF 68 CONTINUE GOTO 212 ENDIF IF(I*J.EQ.0)GOTO 211 WRITE(6,666)V(IABS(I),IABS(J)),' au' WRITE(6,666)V(IABS(I),IABS(J))*cm,' cm-1' if(dabs(V(IABS(I),IABS(J))).gt.1.0d-8) 1WRITE(6,666)1.0d7/(V(IABS(I),IABS(J))*cm),' nm' 666 format(f15.6,a5) IF (I.LT.0) THEN I=IABS(I) J=IABS(J) WRITE(6,*)' NEW V(I,J):' READ(5,*)V(I,J) V(J,I)=V(I,J) ENDIF GOTO 212 211 return end subroutine srm(A,tau,ix) c set rotational matrix implicit none integer*4 ix,i,j,iy,iz real*8 A(3,3),tau,s,c,pi pi=4.0d0*datan(1.0d0) s=sin(tau*pi/180.0d0) c=cos(tau*pi/180.0d0) do 1 i=1,3 do 2 j=1,3 2 A(i,j)=0.0d0 1 A(i,i)=1.0d0 iy=ix+1 if(iy.gt.3)iy=1 iz=iy+1 if(iz.gt.3)iz=1 A(iy,iy)=c A(iz,iz)=c A(iy,iz)=s A(iz,iy)=-s return end subroutine rotv(R,rb,N,M,ii,U,tx,ty,tz) c rotate and translate ensemble of coordinates implicit none integer*4 N,M,ii,i real*8 R(N,3),rb(M,3),ai,x,y,z,U(3,3),tx,ty,tz ai=dble(ii-1) do 1 i=1,N x=R(i,1) y=R(i,2) z=R(i,3) rb(i+N*(ii-1),1)=U(1,1)*x+U(1,2)*y+U(1,3)*z+tx*ai rb(i+N*(ii-1),2)=U(2,1)*x+U(2,2)*y+U(2,3)*z+ty*ai 1 rb(i+N*(ii-1),3)=U(3,1)*x+U(3,2)*y+U(3,3)*z+tz*ai return end function vijf(xi,yi,zi,xj,yj,zj,uix,uiy,uiz,ujx,ujy,ujz,EPS0) implicit none real*8 vijf,xi,yi,zi,xj,yj,zj,uix,uiy,uiz,ujx,ujy,ujz, 1uij,uir,ujr,EPS0,D2,rij,D5,v,EPSR,xji,yji,zji xji=xj-xi yji=yj-yi zji=zj-zi uij=uix*ujx+uiy*ujy+uiz*ujz uir=uix*xji+uiy*yji+uiz*zji ujr=ujx*xji+ujy*yji+ujz*zji D2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 rij=dsqrt(D2) D5=D2*D2*rij v=(uij*D2-3.0d0*uir*ujr)/D5 EPSR=EPS0 IF(EPS0.LT.0.0D0) EPSR=-EPS0*rij vijf=v/EPSR return end c ===================================================================== subroutine vcorrect12(V,N,N0) implicit none integer*4 N,N0,na,nat,ia,ib,ix,ip,nb,i,ic,j real*8 V(N0,N0),ta,va,phi,psi,vold1,cm,vnew1, 1vold2,vnew2,vold3,vnew3,vold4,vnew4 real*8,allocatable::r(:) integer*4,allocatable::m(:) real*8,allocatable::c(:,:),d(:) character*80 s80 cm=219470.0d0 WRITE(6,605) 605 format(/,'Correcting amide I - amide I neighbors') open(9,file='FILE.X',status='old') read(9,*) read(9,*)nat allocate(r(3*nat)) do 1 ia=1,nat 1 read(9,*)r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) close(9) write(6,*)nat,' atoms in FILE.X' open(9,file='FFIT.PAR',status='old') do 4 ic=1,4 read(9,*)nb if(ic.eq.1)allocate(c(4,nb),d(nb)) do 4 i=1,nb read(9,909)s80 909 format(a80) 4 read(s80(38:50),*)c(ic,i) close(9) write(6,*)nb,'AI AII basis coefficients in FFIT.PAR' open(9,file='AMIDE.LST',status='old') read(9,*)na allocate(m(6*na)) do 2 ia=1,na 2 read(9,*)m(1+6*(ia-1)),(m(ix+6*(ia-1)),ix=1,6) close(9) write(6,*)na,' amides in AMIDE.LST' if(N.ne.2*na)call report('N <> 2na') c order C - CO - NH-C amide groups: c 6 21 34 5 c find neighboring amides ip=0 do 3 ia=1,na do 3 ib=1,na if(m(6+6*(ib-1)).eq.m(5+6*(ia-1)))then c C - CO - NH-(C = C) - CO - NH-C c 6 21 34 (5 = 6) 21 34 5 c ia | ib c N-end(1) C-end(2) phi=ta(m(2+6*(ia-1)),m(3+6*(ia-1)),m(5+6*(ia-1)),m(2+6*(ib-1)),r) psi=ta(m(3+6*(ia-1)),m(5+6*(ia-1)),m(2+6*(ib-1)),m(3+6*(ib-1)),r) c c AI AI do 101 i=1,nb 101 d(i)=c(1,i) i=2*ia-1 j=2*ib-1 vold1=V(i,j) vnew1=va(phi,psi,nb,d)/cm V(i,j)=vnew1 V(j,i)=V(i,j) c c AII AII do 102 i=1,nb 102 d(i)=c(2,i) i=2*ia j=2*ib vold2=V(i,j) vnew2=va(phi,psi,nb,d)/cm V(i,j)=vnew2 V(j,i)=V(i,j) c c AI AII 1 2 do 103 i=1,nb 103 d(i)=c(3,i) i=2*ia-1 j=2*ib vold3=V(i,j) vnew3=va(phi,psi,nb,d)/cm V(i,j)=vnew3 V(j,i)=V(i,j) c c AI AII 2 1 do 104 i=1,nb 104 d(i)=c(4,i) i=2*ia j=2*ib-1 vold4=V(i,j) vnew4=va(phi,psi,nb,d)/cm V(i,j)=vnew4 V(j,i)=V(i,j) ip=ip+1 write(6,609)ip,ia,ib,phi,psi,vold1*cm,vnew1*cm, 1 vold2*cm,vnew2*cm,vold3*cm,vnew3*cm,vold4*cm,vnew4*cm 609 format(3i5,10f12.3) endif 3 continue write(6,*)ip,' x 4 elements replaced' 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) 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