PROGRAM new5 implicit none integer*4 A,B,AP,NOAT,NQ,I,J,L,III,IR,LP,IG,ip,ia,ID,ix,iy,iz REAL*8 j1,ELD(3),MGD(3),MDN(3),MDD(3),MDT(3),debye, 3 ME(3),MM(3),AI(3,3),TI(3),aulambda,auk,clight,CM,AMU, 4BOHR,pi,FD,FC,WLIM,W,WR,cmax,ABSW,EPS, 5x,y,z,r2,rk,sd,rpx,rpy,rpz,Q(3,3), 6ssp,spm,sdm,QLD(6),x0,y0,z0,pr, 7prp,pa,pap,ot,rkc,xl,yl,zl,po, 8qaxx,qaxy,qaxz,qayx,qayy,qayz,qazx,qazy,qazz, 8qpxx,qpxy,qpxz,qpyx,qpyy,qpyz,qpzx,qpzy,qpzz, 8olxx,olxy,olxz,olyx,olyy,olyz,olzx,olzy,olzz, 8opxx,opxy,opxz,opyx,opyy,opyz,opzx,opzy,opzz logical lex,LPSI,lexq,LAPT,lexe,LQPT character*20 s20 REAL*8,allocatable::SV(:,:,:),R(:,:),FRQ(:),qq(:,:), 1PP(:,:,:),II(:,:,:),JJ(:,:,:),PV(:,:,:),PQ(:,:),QSX(:,:), 1AQ(:,:),j0(:),j1rk(:),RDO(:),RGO(:),DS(:),DSX(:,:),MSX(:,:) character*8 QF(3) data QF/'DOGX.TAB','DOGY.TAB','DOGZ.TAB'/ c pi=4.0d0*atan(1.0d0) debye=2.541732d0 po=dsqrt(2.0d0*pi) clight=137.03604d0 CM=219474.63d0 AMU=1822.887D0 c FD =SQRT(388.9219d0) FD =debye*dsqrt(CM/AMU/2.0d0) c FC=131.14D-8 FC=debye*dsqrt(2.0d0/CM/AMU)/clight BOHR=0.52917715d0 rkc=2.0d0*pi*1.0d-8*BOHR WLIM=0.000001d0 LPSI=.false. LAPT=.false. LQPT=.false. WRITE(*,*)' COMBINATION OF POLARIZABILITY TENSORS' WRITE(*,*)' AND THE S-VECTORS' WRITE(*,*) WRITE(*,6666) 6666 FORMAT(/, 1' INPUT FILES: F .INP ... VIBRATIONAL PARAMETERS',/, 1' FILE.TEN ... TENSORS',/, 1' NEW5.OPT ... options',/, 3' OUTPUT FILE: DOG.TAB',/, 4' DIP.MOM',/) c Options: oooooooooooooooooooooooooooooooooooooooooooooooooo OPEN(19,FILE='DIP.MOM') inquire(file='NEW5.OPT',exist=lex) if(lex)then open(17,file='NEW5.OPT') 1 read(17,170,end=179,err=179)s20 170 format(a20) write(6,*)s20 if(s20(1:4).eq.'LPSI')read(17,*)LPSI if(s20(1:4).eq.'LAPT')read(17,*)LAPT if(s20(1:4).eq.'LQPT')read(17,*)LQPT goto 1 179 close(17) endif c c deimensions from S-matrix: call readsd(NQ,NOAT) c S-matrix: SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS allocate(SV(NOAT,NQ,3),R(3,NOAT),FRQ(NQ)) call reads(NQ,NOAT,SV,FRQ,R) C READ DIPOLE DERIVATIVES = TENSOR P(A,B): DDDDDDDDDDDDDDDDDD allocate(PP(NOAT,3,3),JJ(NOAT,3,3),PV(NOAT,3,3), 1II(NOAT,3,3)) OPEN(15,FILE='FILE.TEN',STATUS='OLD') READ (15,*)NOAT DO 10 L=1,NOAT DO 10 J=1,3 10 READ (15,*) (PP(L,J,I),I=1,3) if(LAPT)then write(6,*)' AAT from APT ("APT model") requested' DO 2203 L=1,NOAT DO 2203 J=1,3 DO 2203 I=1,3 II(L,J,I)=0.0d0 JJ(L,J,I)=0.0d0 DO 2203 ID=1,3 DO 2203 IG=1,3 2203 II(L,J,I)=II(L,J,I)+0.25d0*EPS(I,IG,ID)*R(IG,L)*PP(L,J,ID) else C C READ TENSORS I(A,B) AND J(A,B) FOR THE 0 0 0 ORIGIN: C (=ELECTRONIC & NUCLEAR TERMS) C DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (II(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (JJ(L,J,I),I=1,3) C C READ VELOCITY REPRESENTATION OF DIPOLE DERIVATIVES C DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) c Distributed origin gauge (kind of obsolete): CALL VCDD0(II,PV,PP,R,NOAT) WRITE(*,*)' Tensors transformed' endif CLOSE(15) WRITE(*,*)' Dipole derivatives read in' c quadrupole derivatives: QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ allocate(qq(3*NOAT,6)) c some older format of quadrupoles, probably never used: inquire(file='FILEQQ.TEN',exist=lexq) if(lexq)then write(6,*)'quadrupole detivatives found' open(15,file='FILEQQ.TEN') read(15,*) read(15,*) write(6,*)'atoms' do 130 i=1,3*NOAT 130 read(15,*)(qq(i,j),j=1,2),(qq(i,j),j=1,6) c xx yy zz xy xz yz, in atomic units close(15) endif c quadrupoles for anisotropic spectra: inquire(file='FILE.QEN',exist=lexe) if(lexe)call reqd('FILE.QEN',qq,NOAT) if(LQPT)then write(6,*)' Quadrupole derivs from APT requested' DO 2204 L=1,NOAT DO 2204 J=1,3 qq(J+3*(L-1),1)=PP(L,J,1)*R(1,L)+PP(L,J,1)*R(1,L) qq(J+3*(L-1),2)=PP(L,J,2)*R(2,L)+PP(L,J,2)*R(2,L) qq(J+3*(L-1),3)=PP(L,J,3)*R(3,L)+PP(L,J,3)*R(3,L) qq(J+3*(L-1),4)=PP(L,J,1)*R(2,L)+PP(L,J,2)*R(1,L) qq(J+3*(L-1),5)=PP(L,J,1)*R(3,L)+PP(L,J,3)*R(1,L) 2204 qq(J+3*(L-1),6)=PP(L,J,2)*R(3,L)+PP(L,J,3)*R(2,L) call wrqd0('LQPT.QEN',qq,NOAT) endif CALL INERTIA(R,AI,TI,NOAT) WRITE (19,5019)(TI(I),I=1,3) 5019 FORMAT( 1' POLARIZATIONS ACCORDING TO THE INERTIA TENSOR AXES',/, 2' ( ',3F15.4, ')',/, 3' MODE ELECTRIC ', 4' MAGNETIC [%] ',/,80(1H-),/) C allocate(AQ(3*NOAT,3),PQ(3*NOAT,3)) DO 411 I=1,3*NOAT DO 411 A=1,3 AQ(I,A)=0.0d0 411 PQ(I,A)=0.0d0 cmax=0.0d0 ip=-1 if(LPSI)then C Put axial tensor origins to atoms: do 2201 L=1,NOAT DO 2201 J=1,3 DO 2201 A=1,3 DO 2201 B=1,3 DO 2201 IG=1,3 2201 II(L,J,A)=II(L,J,A)-0.25d0*EPS(A,IG,B)*R(IG,L)*PP(L,J,B) c c Put quadrupole derivatives to atoms: if(lexq)call TRQDB(qq,R,NOAT,1,PP) endif allocate(DS(NQ),RDO(NQ),RGO(NQ)) if(lexe.or.LQPT)allocate(DSX(NQ,3),QSX(NQ,3),MSX(NQ,3)) DO 413 I=1,NQ DS(I)=0.0d0 RDO(I)=0.0d0 413 RGO(I)=0.0d0 if(LPSI)then write(*,*)'Arbitrary wavelength calculation' c c one-atomic contributions: do 21 L=1,NOAT do 21 B=1,3 DO 21 I=1,NQ DS(I)=DS(I)+ 1 (PP(L,1,B)*SV(L,I,1)+PP(L,2,B)*SV(L,I,2)+PP(L,3,B)*SV(L,I,3))**2 21 RDO(I)=RDO(I)+ 1 (II(L,1,B)*SV(L,I,1)+II(L,2,B)*SV(L,I,2)+II(L,3,B)*SV(L,I,3))* 1 (PP(L,1,B)*SV(L,I,1)+PP(L,2,B)*SV(L,I,2)+PP(L,3,B)*SV(L,I,3)) c c two-atomic contributions: do 2 L=1,NOAT open(80,file='progress') write(80,*)L,' of ',NOAT,' atoms' close(80) xl=R(1,L) yl=R(2,L) zl=R(3,L) do 2 LP=L+1,NOAT x=xl-R(1,LP) y=yl-R(2,LP) z=zl-R(3,LP) r2=dsqrt(x**2+y**2+z**2) x0=x/r2 y0=y/r2 z0=z/r2 c c get Bessel functions here, to save time allocate(j0(NQ),j1rk(NQ)) do 23 I=1,NQ rk=r2*rkc*FRQ(I) if(rk.gt.cmax)cmax=rk j0(I)=sin(rk)/rk j1=(j0(I)-cos(rk))/rk 23 j1rk(I)=j1/rk c DO 2 AP=1,3 IP=3*(LP-1)+AP rpx=y*PP(LP,AP,3)-z*PP(LP,AP,2) rpy=z*PP(LP,AP,1)-x*PP(LP,AP,3) rpz=x*PP(LP,AP,2)-y*PP(LP,AP,1) if(lexq)then qpxx=qq(IP,1) qpyy=qq(IP,2) qpzz=qq(IP,3) qpxy=qq(IP,4) qpyx=qq(IP,4) qpxz=qq(IP,5) qpzx=qq(IP,5) qpyz=qq(IP,6) qpzy=qq(IP,6) c opxx=qpyx*z0-qpzx*y0 opxy=qpyy*z0-qpzy*y0 opxz=qpyz*z0-qpzz*y0 opyx=qpzx*x0-qpxx*z0 opyy=qpzy*x0-qpxy*z0 opyz=qpzz*x0-qpxz*z0 opzx=qpxx*y0-qpyx*x0 opzy=qpxy*y0-qpyy*x0 opzz=qpxz*y0-qpyz*x0 endif c DO 2 A=1,3 IA=3*(L -1)+A spm=rpx*PP(L,A,1)+rpy*PP(L,A,2)+rpz*PP(L,A,3) pr =PP(L ,A ,1)*x0+PP(L ,A ,2)*y0+PP(L , A,3)*z0 prp=PP(LP,AP,1)*x0+PP(LP,AP,2)*y0+PP(LP,AP,3)*z0 pa =II(L ,A ,1)*x0+II(L ,A ,2)*y0+II(L , A,3)*z0 pap=II(LP,AP,1)*x0+II(LP,AP,2)*y0+II(LP,AP,3)*z0 sd= 1 PP(L,A,1)*PP(LP,AP,1)+PP(L,A,2)*PP(LP,AP,2)+PP(L,A,3)*PP(LP,AP,3) sdm= 1 PP(L,A,1)*II(LP,AP,1)+PP(L,A,2)*II(LP,AP,2)+PP(L,A,3)*II(LP,AP,3) 1+II(L,A,1)*PP(LP,AP,1)+II(L,A,2)*PP(LP,AP,2)+II(L,A,3)*PP(LP,AP,3) c c quadrupolar contributions ot=0.0d0 c xx yy zz xy xz yz, in atomic units if(lexq)then qaxx=qq(IA,1) qayy=qq(IA,2) qazz=qq(IA,3) qaxy=qq(IA,4) qayx=qq(IA,4) qaxz=qq(IA,5) qazx=qq(IA,5) qayz=qq(IA,6) qazy=qq(IA,6) c olxx=qayx*z0-qazx*y0 olxy=qayy*z0-qazy*y0 olxz=qayz*z0-qazz*y0 olyx=qazx*x0-qaxx*z0 olyy=qazy*x0-qaxy*z0 olyz=qazz*x0-qaxz*z0 olzx=qaxx*y0-qayx*x0 olzy=qaxy*y0-qayy*x0 olzz=qaxz*y0-qayz*x0 c ot=((olxx*PP(LP,AP,1)+opxx*PP(L,A,1))*x0+ 1 (olxy*PP(LP,AP,1)+opxy*PP(L,A,1))*y0+ 1 (olxz*PP(LP,AP,1)+opxz*PP(L,A,1))*z0+ 1 (olyx*PP(LP,AP,2)+opyx*PP(L,A,2))*x0+ 1 (olyy*PP(LP,AP,2)+opyy*PP(L,A,2))*y0+ 1 (olyz*PP(LP,AP,2)+opyz*PP(L,A,2))*z0+ 1 (olzx*PP(LP,AP,3)+opzx*PP(L,A,3))*x0+ 1 (olzy*PP(LP,AP,3)+opzy*PP(L,A,3))*y0+ 1 (olzz*PP(LP,AP,3)+opzz*PP(L,A,3))*z0)/r2 endif do 2 I=1,NQ ssp=SV(L,I,A)*SV(LP,I,AP) j1=j1rk(I) RDO(I)=RDO(I)+0.75d0*ssp*(-spm*j1+2.0d0*sdm*(1.0d0-j1) 1+(ot+2.0d0*(pap*pr+prp*pa))*(3.0d0*j1-j0(I))) 2 DS(I)=DS(I)+3.0d0*ssp*(sd*(j0(I)-j1)-pr*prp*(j0(I)-3.0d0*j1)) c DO 25 I=1,NQ W=FRQ(I) ABSW=DABS(W) WR=DSQRT(ABSW) IF(ABSW.GT.WLIM)DS(I)=DS(I)/WR**2*FD**2 RDO(I)=RDO(I)*FD*FC 25 RGO(I)=RDO(I) else write(*,*)'Dipolar approximation route' DO 41 I=1,NQ W=FRQ(I) ABSW=DABS(W) IF(ABSW.GT.WLIM)then DO 391 III=1,3 QLD(III )=0.0d0 QLD(III+3)=0.0d0 ME(III)=0.0d0 391 MM(III)=0.0d0 WR=DSQRT(ABSW) if(lexq.or.lexe.or.LQPT)then C QUADRUPOLE MOMENT TOTAL: DO 57 B=1,6 DO 57 L=1,NOAT DO 57 A=1,3 57 QLD(B) = QLD(B) + qq(3*(L-1)+A,B)*SV(L,I,A) Q(1,1)=QLD(1) Q(2,2)=QLD(2) Q(3,3)=QLD(3) Q(1,2)=QLD(4) Q(2,1)=QLD(4) Q(1,3)=QLD(5) Q(3,1)=QLD(5) Q(2,3)=QLD(6) Q(3,2)=QLD(6) endif DO 56 B=1,3 MGD(B)=0.0d0 ELD(B)=0.0d0 MDN(B)=0.0d0 MDD(B)=0.0d0 MDT(B)=0.0d0 DO 56 L=1,NOAT DO 56 A=1,3 C ELECTRONIC MOMENT TOTAL: ELD(B) = ELD(B) + PP(L,A,B)*SV(L,I,A) C MAGNETIC DIPOLE NUCLEAR SAME FOR DO/GO MDN(B) = MDN(B) + JJ(L,A,B) *SV(L,I,A) C MAGNETIC DIPOLE ELECTRONIC DO MDD(B) = MDD(B) + II(L,A,B)*SV(L,I,A) C MAGNETIC DIPOLE TOTAL DO MDT(B) = MDT(B) + (JJ(L,A,B)+ II(L,A,B))*SV(L,I,A) C MAGNETIC DIPOLE TOTAL GO 56 MGD(B) = MGD(B) + (II(L,A,B)+JJ(L,A,B)) *SV(L,I,A) c c APT,AAT in normal modes: do 561 B=1,3 AQ(I,B)=MDT(B) 561 PQ(I,B)=ELD(B) C DO 51 B=1,3 IF(ABSW.GT.WLIM)ELD(B)=ELD(B)/WR*FD MDT(B)=MDT(B)*FC*WR MDN(B)=MDN(B)*FC*WR MDD(B)=MDD(B)*FC*WR 51 MGD(B)=MGD(B)*FC*WR if((lexe.or.LQPT).and.ABSW.GT.WLIM)then DO 59 A=1,3 DO 59 B=1,3 59 Q(A,B)=Q(A,B)*FD/WR*0.5d0 c the half is because Q is sum(i) qi ri ri c but we need quadrupole as (1/2) of it endif C DS(I) = ELD(1)*ELD(1)+ELD(2)*ELD(2)+ELD(3)*ELD(3) RDO(I)= MDT(1)*ELD(1)+MDT(2)*ELD(2)+MDT(3)*ELD(3) RGO(I)= MGD(1)*ELD(1)+MGD(2)*ELD(2)+MGD(3)*ELD(3) c c X, Y and Z beam directions if quadrupole present: if(lexe.or.LQPT)then c wavelength in atomic units aulambda=1.0d0/FRQ(I)*1.0d8/BOHR c wave vector in atomic units auk=2.0d0*pi/aulambda do 58 ix=1,3 iy=ix+1 if(iy.eq.4)iy=1 iz=iy+1 if(iz.eq.4)iz=1 DSX(I,ix)= 1.5d0*(ELD(iy)*ELD(iy) +ELD(iz)*ELD(iz) ) QSX(I,ix)=-auk*1.5d0*(ELD(iy)*Q(iz,ix)-ELD(iz)*Q(ix,iy)) 58 MSX(I,ix)= 1.5d0*(ELD(iy)*MDT(iy) +ELD(iz)*MDT(iz) ) endif C DO 39 III=1,3 DO 39 J=1,3 MM(III)=MM(III)+AI(J,III)*MDT(J) 39 ME(III)=ME(III)+AI(J,III)*ELD(J) C CALL NORM(ME) CALL NORM(MM) C DO 4 III=1,3 MM(III)=MM(III)*MM(III)*100.0d0 4 ME(III)=ME(III)*ME(III)*100.0d0 C if(lexq)then WRITE(19,1920)I,W,(ME(IR),IR=1,3),(MM(IR),IR=1,3), 1 (QLD(IR),IR=1,6) WRITE(17,522)(ELD(IR),IR=1,3),(MDT(IR),IR=1,3), 1 (QLD(IR),IR=1,6) else WRITE(19,1920)I,W,(ME(IR),IR=1,3),(MM(IR),IR=1,3) 1920 FORMAT(I4,F7.1,4(3F10.0,' ')) WRITE(17,522)(ELD(IR),IR=1,3),(MDT(IR),IR=1,3) 522 FORMAT(12G13.5) endif endif 41 continue endif CLOSE(17) c OPEN(18,FILE='DOG.TAB') WRITE (18,501) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') DO 412 I=1,NQ 412 WRITE (18,523) I,FRQ(I),DS(I),RDO(I),RGO(I),po*FRQ(I)*DS(I) 523 FORMAT(I5,f16.7,4G15.6) WRITE(18,1819) 1819 FORMAT('-------------------------------------------------') CLOSE(18) if(lexe.or.LQPT)then do 60 ix=1,3 OPEN(18,FILE=QF(ix)) WRITE (18,501) DO 61 I=1,NQ 61 WRITE (18,523) I,FRQ(I),DSX(I,ix),QSX(I,ix)+MSX(I,ix) WRITE(18,1819) write(6,*)QF(ix)//' written' 60 CLOSE(18) endif C WRITE(19,1819) CLOSE(19) if(LPSI)write(6,*)'cmax: ',cmax call WRITEQTEN(PQ,AQ,NOAT) WRITE(*,*)' >> PROGRAM TERMINATED <<' END SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT none INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA,NAT real*8 VCD(NAT,3,3),VEL(NAT,3,3),DIP(NAT,3,3), 1C(3,NAT),SKEW,EPS DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END FUNCTION EPS(I,J,K) integer*4 I,J,K REAL*8 EPS EPS=0.0D0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0D0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0D0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0D0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0D0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0D0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0D0 RETURN END SUBROUTINE INERTIA(R,A,TI,N) IMPLICIT integer*4(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION R(3,N),A(3,3),T(3,3),TI(3),CM(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 TRED12(3,T,A,TI,2,IERR) 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 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(N) real*8,allocatable::E(:) allocate(E(N)) 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 TQL12 (N,Z,D,IERR,ICODE,E) RETURN 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END 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 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 NORM(A) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3) AN=DSQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3)) IF(ABS(AN).LT.1.0d-30)RETURN A(1)=A(1)/AN A(2)=A(2)/AN A(3)=A(3)/AN RETURN END SUBROUTINE WRITEQTEN(P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(NAT*3,3),A(NAT*3,3) OPEN(15,FILE='FILE.Q.TEN') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5,' normal mode AAT AAP derivatives') DO 10 L=1,NAT*3 10 WRITE(15,1501) (P(L,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 221 L=1,NAT*3 221 WRITE(15,1501) (A(L,I),I=1,3),L DO 230 L=1,NAT*3 230 WRITE(15,1501) (0.0d0,I=1,3),L DO 100 L=1,NAT*3 100 WRITE(15,1501) (P(L,I),I=1,3),L CLOSE(15) RETURN END SUBROUTINE TRQDB(Q,X,NAT,ic,P) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(3*NAT,6),X(3,*),P(NAT,3,3) if(ic.eq.0)then write(6,*)'no change in quadrupole derivatives' return else if(ic.eq.1)then sign= 0.5d0 write(6,*)'quadrupole derivatives put on atoms' else sign=-0.5d0 write(6,*)'quadrupole derivatives shifted back' endif endif DO 10 L=1,NAT DO 10 J=1,3 LL=3*(L-1)+J Q(LL,1)=Q(LL,1)+sign*(X(1,L)*P(L,J,1)+P(L,J,1)*X(1,L)) Q(LL,2)=Q(LL,2)+sign*(X(2,L)*P(L,J,2)+P(L,J,2)*X(2,L)) Q(LL,3)=Q(LL,3)+sign*(X(3,L)*P(L,J,3)+P(L,J,3)*X(3,L)) Q(LL,4)=Q(LL,4)+sign*(X(1,L)*P(L,J,2)+P(L,J,1)*X(2,L)) Q(LL,5)=Q(LL,5)+sign*(X(1,L)*P(L,J,3)+P(L,J,1)*X(3,L)) 10 Q(LL,6)=Q(LL,6)+sign*(X(2,L)*P(L,J,3)+P(L,J,2)*X(3,L)) c xx yy zz xy xz yz in atomic units write(6,*)'ok' return end subroutine reads(NQ,NOAT,SV,FRQ,R) implicit none integer*4 NQ,NOAT,K,I,J real*8 SV(NOAT,NQ,3),R(3,NOAT),FRQ(*),BOHR BOHR=0.52917715d0 OPEN(17,FILE='F.INP',STATUS='OLD') READ(17,*)NQ,NOAT,NOAT DO 3 K=1,NOAT READ(17,*)R(1,K),(R(I,K),I=1,3) DO 3 I=1,3 3 R(I,K)=R(I,K)/BOHR WRITE(*,*)' Geometry read in' READ(17,*) DO 4 K=1,NOAT DO 4 I=1,NQ 4 READ (17,*)(SV(K,I,J),J=1,2),(SV(K,I,J),J=1,3) READ(17,*) READ(17,*)(FRQ(I),I=1,NQ) c1717 FORMAT(3F15.7) WRITE(*,*)' S-matrix read in' return end subroutine readsd(NQ,NOAT) implicit none integer*4 NQ,NOAT OPEN(17,FILE='F.INP',STATUS='OLD') READ(17,*)NQ,NOAT,NOAT close(17) return end subroutine reqd(fn,q,nat) implicit none character*(*) fn integer*4 i,ia,ix,nat real*8 q(3*nat,6) open(9,file=fn) do 1 i=1,9 1 read(9,*) do 2 ia=1,nat do 2 ix=1,3 2 read(9,93)(q(ix+3*(ia-1),i),i=1,6) 93 format(9x,3e13.4) close(9) write(6,*)fn//' read' return end subroutine wrqd0(fn,d,nat) implicit none character*(*) fn integer*4 i,ia,ix,nat real*8 d(3*nat,6),q(6) character*1 xyz(3) data xyz/'x','y','z'/ do 1 i=1,6 1 q(i)=0.0d0 open(9,file=fn) write(9,90)nat,'XX',q(1),'YY',q(2),'ZZ',q(3), 1 'XY',q(4),'XZ',q(5),'YZ',q(6) 90 format(i6,' atoms',/, 1 ' Equilibrium quadrupole, sum(i)qi ri ri, atomic units:',/, 12(3(3x,a2,2x,f19.4),/)) write(9,91) 91 format(' Quadrupole derivatives, atomic units: ',/, 1 ' Atom Coord XX YY ZZ ',/, 1 ' XY XZ YZ ',/, 1 ' ------------------------------------------------') do 2 ia=1,nat do 2 ix=1,3 if(ix.eq.1)then write(9,92)ia,xyz(ix),(d(ix+3*(ia-1),i),i=1,6) 92 format(i5,3x,a1,3e13.4,/,9x,3e13.4) else write(9,93)xyz(ix),(d(ix+3*(ia-1),i),i=1,6) 93 format(8x,a1,3e13.4,/,9x,3e13.4) endif 2 continue close(9) write(6,*)fn//' written' return end