PROGRAM newa IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0) REAL*8 II DIMENSION kt(NA0),r(3,NA0),XM(NA0),iso(NA0),riso(NA0), 1wexp(NA0),FCAR(MX3,MX3),CHAR(NA0),SMAT(MX3,MX3),FRQ(MX3), 1PP(NA0,3,3),II(NA0,3,3), 2ALPHA(MX3,3,3),A(MX3,3,3,3),G(MX3,3,3) logical lex,lopt CHARACTER*5 s5 common/opts/lopt(300) BOHR=0.52917715d0 C do 3 i=1,300 3 lopt(300)=.false. c 1 . calculate frequencies c 2 . generate FTRY.INP c 3 . substitute acidic hydrogens by deuteria c 4 . read in FTRY.INP c 5 . project zero vibrations c 6 . write F.INP c 7 . do FPC c 8 . call new5/abs/vcd c 9 . write DOG.TAB c 10 . write DIP.MOM c 11 . read F.INP c 12 . write dipoles in F.INP c 13 . call new6/Raman/ROA c 14 . polar ROA c 15 . empty c 16 . write ROA.TAB c 17 . large Raman print c 18 . generate abs spectrum c 19 . write abs spectrum do d.prn c 20 . generate vcd spectrum c 21 . write vcd spectrum do d.prn c 22 . generate Raman spectrum c 23 . write Raman spectrum do d.prn c 24 . generate ROA spectrum c 25 . write ROA spectrum do d.prn c c WMIN, WMAX min, max frequency c EXCA excitation wavelength for Raman in A lopt(1)=.true. lopt(5)=.true. lopt(6)=.true. lopt(9)=.true. EXCA=4880.0d0 WRITE(6,6000) 6000 format(' Spectra generation',/,/, 1 ' Input: NEWA.OPT options',/, 2 ' FILE.X geometry',/, 2 ' FILE.FC force constants',/, 2 ' FILE.TEN tensors',/) inquire(file='NEWA.OPT',exist=lex) if(lex)then open(20,file='NEWA.OPT') 1 read(20,200,end=99,err=99)s5 if(s5.eq.'LFREQ')read(20,*)lopt(1) if(s5(1:4).eq.'FTRY')read(20,*)lopt(2) if(s5(1:4).eq.'DEUT')read(20,*)lopt(3) if(s5(1:4).eq.'RFTR')read(20,*)lopt(4) if(s5(1:4).eq.'PROJ')read(20,*)lopt(5) if(s5(1:5).eq.'F.INP')read(20,*)lopt(6) if(s5(1:3).eq.'FPC')read(20,*)lopt(7) if(s5(1:3).eq.'ABS')read(20,*)lopt(8) if(s5(1:3).eq.'DOG')read(20,*)lopt(9) if(s5(1:3).eq.'MOM')read(20,*)lopt(10) if(s5(1:3).eq.'R.F')read(20,*)lopt(11) if(s5(1:3).eq.'W.F')read(20,*)lopt(12) if(s5(1:3).eq.'RAM')read(20,*)lopt(13) if(s5(1:3).eq.'POL')read(20,*)lopt(14) if(s5(1:3).eq.'ROA')read(20,*)lopt(16) if(s5(1:2).eq.'##')then read(20,*)io read(20,*)lopt(io) endif if(s5(1:4).eq.'WMIN')read(20,*)wmin if(s5(1:4).eq.'WMAX')read(20,*)wmax if(s5(1:4).eq.'EXCA')read(20,*)EXCA if(s5(1:4).eq.'ISOT')then read(20,*)ISOT read(20,*)(iso(i),i=1,ISOT) read(20,*)(riso(i),i=1,ISOT) endif goto 1 200 format(a5) 99 close(20) endif OPEN(55,FILE='FILE.X',STATUS='OLD') READ(55,*) read(55,*)N if(N.gt.NA0)call report('too many atoms') do 2 I=1,N 2 READ(55,*)kt(I),(r(ix,I),ix=1,3) CLOSE(55) call INPTRY(XM,NA0,N,kt,r,iso,riso,ISOT,wexp,nm,CHAR) DO 35 K=1,N DO 35 I=1,3 35 R(I,K)=R(I,K)/BOHR C Read-in cartesian FF if(lopt(1))then OPEN(20,FILE='FILE.FC',STATUS='OLD') CALL READFF(MX3,N*3,FCAR) CLOSE(20) call new4(N,FCAR,CHAR,XM,wexp,nm,SMAT,FRQ) endif c if(lopt(11))then OPEN(17,FILE='F.INP') READ(17,*)nm DO 36 K=1,N 36 READ(17,*) READ(17,*) DO 45 LL=1,N DO 45 I=1,nm L3=3*(LL-1) 45 READ(17,*)idumm,idumm,(SMAT(L3+ix,I),ix=1,3) READ(17,*) READ(17,1717)(FRQ(I),I=1,nm) 1717 FORMAT(6F11.3) CLOSE(17) WRITE(*,*)' S-matrix read in' endif c c absorption/vcd spectra if(lopt(8))then 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) DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (II(L,J,I),I=1,3) CLOSE(15) write(6,*)'FILE.TEN read' call new5(PP,II,NOAT,r,SMAT,FRQ,nm) endif if(lopt(13))then CALL READTEN(MX3,NAT,ALPHA,G,A) call new6(SMAT,wmin,wmax,ALPHA,G,A,N,EXCA,nm,FRQ,r) endif stop end subroutine INPTRY(XM,NA0,N,kt,r,iso,riso,ISOT,wexp,nm,CHAR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MENDELEV=89) DIMENSION XM(NA0),CHAR(*),kt(*),r(3,NA0),iso(*),riso(*),wexp(*) CHARACTER*4 CODE dimension amas(MENDELEV) data amas/1.007825,4.003, 2 6.941, 9.012, 10.810,12.000,14.003074,15.994915,18.9984,20.179, 3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ logical lopt common/opts/lopt(300) C do 1 i=1,N CHAR(i)=0.0d0 1 XM(i)=amas(kt(i)) nm=3*N do 3 i=1,nm 3 wexp(i)=0.0d0 do 2 i=1,ISOT 2 XM(iso(i))=riso(i) c c Substitute all acidic Hs by Ds if(lopt(3))then ND=0 do 8 i=1,N a=r(1,i) y=r(2,i) z=r(3,i) if(kt(i).eq.1)then do 31 j=1,N it=kt(j) dd=sqrt((a-r(1,j))**2+(y-r(2,j))**2+(z-r(3,j))**2) if(it.eq.7.or.it.eq.8.or.it.eq.9.or.it.eq.16.or.it.eq.17.or. 1 it.eq.35.or.it.eq.53)then if(dd.lt.1.23d0)then ND=ND+1 XM(i)=2.014000d0 endif endif 31 continue endif 8 continue write(6,7009)ND 7009 format(i5,' acidic hydrogens have been substituted by deuteria') endif c c generate FTRY.INP if(lopt(2))then CODE='INTY' OPEN(7,FILE='FTRY.INP') WRITE(7,1111)CODE,0,1,0,1,0,1,0 1111 FORMAT(A4,7I4) WRITE(7,333)0 333 FORMAT(20I4) WRITE(7,444)(CHAR(i),I=1,N) 444 format(6f12.6) WRITE(7,333)1 WRITE(7,444)(XM(I),I=1,N) WRITE(7,*)3*N WRITE(7,444)(wexp(i),I=1,3*N) CLOSE(7) write(6,*)'FTRY.INP created' ENDIF c read FTRY.INP if(lopt(4))then OPEN(7,FILE='FTRY.INP') read(7,*) read(7,*) read(7,444)(CHAR(i),I=1,N) read(7,*) read(7,444)(XM(I),I=1,N) read(7,*)nm read(7,444)(wexp(I),I=1,nm) CLOSE(7) write(6,*)'FTRY.INP read' ENDIF return end subroutine new4(NOAT,FCAR,CHAR,ZIM,wexp,nm,SMAT,FRQ) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION ZIM(*),CHAR(*),FCAR(MX3,MX3),wexp(*),FRQ(*), 1SMAT(MX3,MX3) COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA logical lopt common/opts/lopt(300) C NA=3*NOAT do 130 i=1,3 do 130 j=1,NA 130 AS(i,j)=0.0d0 C CHATOT=0.0d0 DO 77 I=1,NOAT 77 CHATOT=CHATOT+CHAR(I) WRITE(*,777)CHATOT,NOAT 777 FORMAT(' TOTAL CHARGE',F20.6,/I6,' atoms') C BM=0.0d0 DO 102 IA=1,NOAT 102 BM=BM+ZIM(IA) WRITE(*,778)BM 778 FORMAT(' TOTAL MASS ',F20.6) DO 220 I=1,NOAT DO 220 J=1,3 220 ZM(3*I-3+J)=1/SQRT(ZIM(I)) C DO 211 I=1,NA 211 ZNU(I)=0.0d0 DO 212 I=1,nm 212 ZNU(I)=wexp(I) C C IF(lopt(5))CALL PROJECT(FCAR,MX3) C C ie=0 ADUM=SUMk(1,ie,FCAR,SMAT) C Calculate the true mass-weighted S-matrix: DO 210 IA=1,NA T=ZM(IA) DO 210 IM=1,NA 210 SMAT(IA,IM)=SMAT(IA,IM)*T DO 213 IM=1,NA 213 FRQ(IM)=AS(1,IM) if(lopt(6))CALL SMATOUT(SMAT,AS,NA) if(lopt(7))CALL FPC(Q,CHAR,SMAT,NOAT,3*NOAT) return END C SUBROUTINE SMATOUT(S,AS,NAT3) C dX = S dQ S(nat3xnint) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION S(MX3,NAT3),AS(3,NAT3) WRITE(*,*)' RANGE OF MODES FOR THE INP FILE (Use the bracket#): ' IRANGE=1 JRANGE=NAT3 NI=JRANGE-IRANGE+1 c c look at the phase and make the biggest element always positive: DO 101 J=IRANGE,JRANGE amx=S(1,J) DO 102 I=2,NAT3 102 if(abs(S(I,J)).gt.abs(amx))amx=S(I,J) if(amx.lt.0)then do 103 I=1,NAT3 103 S(I,J)=-S(I,J) endif 101 continue c OPEN(34,FILE='F.INP') WRITE(34,10)NI,NAT3,NAT3/3 10 FORMAT(3I5) OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*)NAT DO 3 I=1,NAT READ(35,*)IAT,X,Y,Z 3 WRITE(34,11)IAT,X,Y,Z 11 FORMAT(I4,3F12.6) CLOSE(35) WRITE(34,14) 14 FORMAT(' Atom Mode X-disp. Y-disp. Z-disp.') DO 1 I=1,NAT3/3 DO 1 J=IRANGE,JRANGE IU=3*(I-1) K=NAT3-J+1 1 WRITE(34,15)I,K,S(IU+1,J),S(IU+2,J),S(IU+3,J) 15 FORMAT(2I5,3F11.6) WRITE(34,11)NI WRITE(34,13)(AS(1,I),I=IRANGE,JRANGE) 13 format(6F11.3) CLOSE(34) WRITE(*,*)' S-matrix written out ...' RETURN END FUNCTION SUMk(iwr,iee,FCAR,SMAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NA0=1000,MX3=3*NA0) dimension FCAR(MX3,MX3),EWORK(MX3),SMAT(MX3,MX3),FNEW(MX3,MX3) COMMON/ITERA/ 1ZM(MX3), 2Q(MX3),ZNU(MX3), 2AS(3,MX3),NA c wtol=0.01d0 c1=1302.828d0 CALL FAST(NA,FCAR,FNEW,ZM) if(IWR.gt.0)WRITE(6,*)' New cartesian force field done ' IMOD=2 CALL TRED12(MX3,NA,FNEW,SMAT,Q,IMOD,IERR,EWORK) IF(IERR.NE.0)STOP if(iwr.gt.0)WRITE(6,*)' Diagonalization OK' c RMSW=0.0D0 ie=0 DO 310 I=1,NA T=Q(I) SIGN=1.0d0 C Make imaginary frequencies negative: IF(T.LT.0.0d0)SIGN=-1.0d0 T=SIGN*SQRT(ABS(T))*c1 Q(I)=T TP=ZNU(I) DEL=T-TP if(ABS(TP).le.wtol)then DEL=0.0d0 else ie=ie+1 endif RMSW=RMSW+DEL**2 AS(1,I)=T AS(2,I)=TP 310 AS(3,I)=T-TP iee=ie if(ie.gt.0)then RMSW=DSQRT(RMSW/DFLOAT(ie)) else if(NA.gt.0)RMSW=DSQRT(RMSW/DFLOAT(NA)) endif SO1=0.0D0 SO2=0.0D0 DO 5151 J=1,NA SO1=SO1+AS(1,J)*AS(2,J) 5151 SO2=SO2+AS(1,J)*AS(1,J) AK=SO1/SO2 C if(iwr.gt.0)WRITE(*,5555) 5555 FORMAT(1X,'************ NEW FREQUENCIES ******************') DO 311 J=NA,1,-1 a3=0.0d0 if(AS(2,NA-J+1).gt.wtol)a3=AS(3,NA-J+1) 311 if(iwr.gt.0)WRITE(*,6666)J,NA-J+1,(AS(I,NA-J+1),I=1,2),a3 6666 FORMAT(1X, 1 I5,' (',I4,') ...... ',F8.2,' -',F8.2,'/',F8.2) if(iwr.gt.0)WRITE(*,5551) 5551 FORMAT(1X,'************************************************') if(iwr.gt.0)WRITE(6,6667)ie,RMSW,AK 6667 FORMAT(i6,'experimental frequencies found',/, 1 ' >>>>><< >> ',F8.2,' <<<<<< ',F8.6) c SUMk=RMSW RETURN END SUBROUTINE FPC(V,CHAR,S,NOAT,NA3) C PETR BOUR: LAST UPDATE 4-28-92 IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0,MX=MX3/3) DIMENSION V(MX3), S(MX3,NA3), X(3,MX), 1 CHAR(MX), DS(3), SCRT(3) DCOR = 388891.4D0 RCOR = 12221.7D0 OPEN(35,FILE='FILE.X',STATUS='OLD') READ(35,*) READ(35,*)NOAT DO 51 J=1,NOAT 51 READ(35,*)IATDUM,(X(I,J),I=1,3) CLOSE(35) NA=3*NOAT NVIBS=NA C INCLUDE TRANSLATIONS AND ROTATIONS, TOO OPEN(13,FILE='FPC.TAB') WRITE(13,50) 50 FORMAT(3X,'FREQ',8X,'DIP. STRN.*E39',6X,'ROT. STRN.*E45'/ 1 3X,'CM-1',2X,2(6X,' (FR CM)**2 '),/,'---------------') DO 42 I=1,NVIBS DO 784 IKY=1,3 DS(IKY)=0.0D0 SCRT(IKY)=0.0D0 DO 784 JR=1,NOAT DS(IKY)=DS(IKY)+CHAR(JR)*S(3*(JR-1)+IKY,I) DO 784 IB=1,3 DO 784 IG=1,3 784 SCRT(IKY)=SCRT(IKY)+EPS(IKY,IB,IG)*X(IB,JR)*S(3*(JR-1)+IG,I) 1 *CHAR(JR) DT =DS(1)*DS (1)+DS(2)*DS (2)+DS(3)*DS (3) ROT=DS(1)*SCRT(1)+DS(2)*SCRT(2)+DS(3)*SCRT(3) ROTC =+ROT*RCOR DIPOL=0.0D0 IF (ABS(V(I)).GT.0) DIPOL = DT*DCOR/ABS(V(I)) 42 WRITE(13,85) I,V(I),DIPOL,ROTC,DT,ROT 85 FORMAT(I4,F8.2,' ',F13.6,' ',f13.6,' ',2G14.5) WRITE(13,86) 86 FORMAT('---------------------------------------------') CLOSE(13) RETURN END C SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) 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) C CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' RETURN END SUBROUTINE FAST(NA,FCAR,FNEW,ZM) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION FCAR(MX3,NA),FNEW(MX3,NA),ZM(NA) DO 214 I=1,NA UI=ZM(I) DO 214 J=I,NA FNEW(I,J)=FCAR(I,J)*ZM(J)*UI 214 FNEW(J,I)=FNEW(I,J) RETURN END C C SUBROUTINE PROJECT(F,N0) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION F(N0,N0),TEM(MX3,MX3),A(MX3,MX3) C C Using the common here to save memory for A and TEM COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA N=NA WRITE(*,*)'Projecting transl/rotations from the force field...' CALL DOMA(A) DO 1 I=1,N DO 1 J=1,N TEM(I,J)=0.0d0 DO 1 II=1,N 1 TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J) DO 2 I=1,N DO 2 J=1,N F(I,J)=0.0d0 DO 2 II=1,N 2 F(I,J)=F(I,J)+TEM(I,II)*A(J,II) RETURN END SUBROUTINE DOMA(A) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (NA0=1000,MX3=3*NA0,MX=MX3/3) DIMENSION C(3,MX),TEM(MX3,MX3),A(MX3,MX3) COMMON/ITERA/ZM(MX3),Q(MX3),ZNU(MX3),AS(3,MX3),NA N=NA AMACH=0.00000000000001d0 NAT=N/3 OPEN(4,FILE='FILE.X',STATUS='OLD') READ(4,*) READ(4,*)NAT DO 11 I=1,NAT 11 READ(4,*)ID,(C(II,I),II=1,3) CLOSE(4) DO 12 I=1,N DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) N6=6 CALL ORT(A,MX3,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(*,*)N6,' coordinates projected' DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT(A,NDIM,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(NDIM,N6) AMACH=0.00000000000001d0 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/SQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END SUBROUTINE TRED12(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) 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 (MX3,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(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) 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 subroutine new5(PP,II,NOAT,R,SV,FRQ,nm) implicit none integer*4 NA0,A,B,NOAT,I,J,L,III,IR, 1nm PARAMETER (NA0=1000) REAL*8 PP(NA0,3,3),II(NA0,3,3),SV(3*NA0,3*NA0), 1 R(3,NA0),FRQ(*), 2 ELD(3),MDT(3), 3 ME(3),MM(3),AI(3,3),TI(3), 4pi,FD,FC,Z0,WLIM,W,WR,ABSW, 5RDO, 6SVL,sumd logical lopt common/opts/lopt(300) c pi=4.0d0*atan(1.0d0) FD =SQRT(388.9219d0) FC=131.14D-8 Z0=0.0d0 WLIM=0.000001d0 if(lopt(12))then open(17,file='F.INP') 101 read(17,*,end=102,err=102) goto 101 102 continue endif if(lopt(9))OPEN(18,FILE='DOG.TAB') if(lopt(10))OPEN(19,FILE='DIP.MOM') if(lopt(9))WRITE (18,501) if(lopt(10))CALL INERTIA(R,AI,TI,NOAT) if(lopt(10))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-),/) 501 FORMAT (/,' FREQ. ', 1' D [D^2] ',' R [D^-2]DO R - CO I[km/mol]', 2/,'-----------------------------------------------------------') C DO 41 I=1,nm W=FRQ(I) ABSW=DABS(W) WR=DSQRT(ABSW) DO 56 B=1,3 ELD(B)=Z0 MDT(B)=Z0 DO 56 L=1,NOAT DO 56 A=1,3 SVL=SV(3*(L-1)+A,I) ELD(B) = ELD(B) + PP(L,A,B)*SVL 56 MDT(B) = MDT(B) + II(L,A,B)*SVL c DO 51 B=1,3 IF(ABSW.GT.WLIM)ELD(B)=ELD(B)/WR*FD MDT(B)=MDT(B)*FC*WR 51 continue C SUMD = ELD(1)*ELD(1)+ELD(2)*ELD(2)+ELD(3)*ELD(3) RDO = MDT(1)*ELD(1)+MDT(2)*ELD(2)+MDT(3)*ELD(3) C DO 39 III=1,3 ME(III)=Z0 MM(III)=Z0 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(lopt(10))WRITE(19,1920)I,W,(ME(IR),IR=1,3),(MM(IR),IR=1,3) 1920 FORMAT(I4,F7.1,4(3F10.0,' ')) if(lopt(12))WRITE(17,522)(ELD(IR),IR=1,3),(MDT(IR),IR=1,3) 522 FORMAT(12G13.5) IF(ABSW.LT.WLIM)SUMD=Z0 41 if(lopt(9))WRITE (18,523) I,W,SUMD,RDO,RDO,SQRT(PI+PI)*W*SUMD 523 FORMAT(I4,5G15.6) C if(lopt(9))WRITE(18,1819) if(lopt(10))WRITE(19,1819) 1819 FORMAT('-------------------------------------------------') if(lopt(12))CLOSE(17) if(lopt(9))CLOSE(18) if(lopt(10))CLOSE(19) 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 REAL*8(A-H,O-Z) PARAMETER (NA0=1000) DIMENSION R(3,NA0),A(3,3),T(3,3),TI(3),CM(3),E(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,3,T,A,TI,2,IERR,E) 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 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 report(s) character*(*) s write(6,*)s stop end subroutine new6(S,wmin,wmax,ALPHA,G,A,NAT,EXCA,nm,FRQ,r) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION S(MX3,MX3),ALPHA(MX3,3,3),A(MX3,3,3,3),G(MX3,3,3), 3R(3,MX3/3),FRQ(*) LOGICAL LPOLAR logical lopt common/opts/lopt(300) C LPOLAR=lopt(14) IF(LPOLAR)THEN c set A, G to zero: CALL TTTT(MX3,ALPHA,A,G,NAT,3,R) c make common A, G from alpha: CALL TTTT(MX3,ALPHA,A,G,NAT,2,R) write(6,*)' A, G made from alpha !!!' ENDIF CALL TRTEN(MX3,S,3*NAT,ALPHA,A,G) CALL DORAMAN(MX3,nm,ALPHA,G,A,EXCA,wmin,wmax,FRQ) return END C ========================================== SUBROUTINE READTEN(MX3,NAT,ALPHA,G,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(MX3,3,3),G(MX3,3,3),A(MX3,3,3,3) OPEN(2,FILE='FILE.TTT',STATUS='OLD') READ(2,*) READ(2,*)NAT READ(2,*) READ(2,*) DO 1 I=1,3 READ(2,*) DO 1 L=1,NAT DO 1 IX=1,3 IND=3*(L-1)+IX 1 READ(2,*)LDUM,IXDUM,(ALPHA(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 2 I=1,3 READ(2,*) DO 2 L=1,NAT DO 2 IX=1,3 IND=3*(L-1)+IX c last index magnetic, inner 2 READ(2,*)LDUM,IXDUM,(G(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 3 I=1,3 DO 3 J=1,3 READ(2,*) DO 3 L=1,NAT DO 3 IX=1,3 IND=3*(L-1)+IX 3 READ(2,*)LDUM,IXDUM,(A(IND,I,J,K),K=1,3) CLOSE(2) WRITE(6,*)NAT,' atoms, tensors read in from FILE.TTT' RETURN END C ============================================================== SUBROUTINE DORAMAN(NS,N,ALPHA,G,A,EXCA,wmin,wmax,E) c c supposedly G/W is supplied c IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) DIMENSION ALPHA(NS,3,3),G(NS,3,3),A(NS,3,3,3),E(*) logical lopt common/opts/lopt(300) c c if(lopt(16))OPEN(4,FILE='ROA.TAB') if(lopt(16))WRITE(4,3304) 3304 FORMAT( 1'MODE FREQ Ramx (180) Ramy Ramtot CID90', 2' CID-X CID180 DEP ROA180 RamanDCPI',/, 3' cm-1 A^4/AMU x10^4', 4' x10^4 x10^4 A^5/AMUx1e4 au.10^4') if(lopt(16))write(4,3002) 3002 FORMAT(80(1H-)) c DO 1 IQ=1,N ECM=E(IQ) IF(ABS(ECM).GT.0.01d0)THEN if(wmin.eq.0.0d0.or.ECM.gt.wmin)then if(wmax.eq.0.0d0.or.ECM.lt.wmax)then 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+ALPHA(IQ,I,I) DO 3 J=1,3 SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(IQ,J,J) SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(IQ,I,J) SAG0=SAG0+ALPHA(IQ,I,I)* G(IQ,J,J) SAG1=SAG1+ALPHA(IQ,I,J)* G(IQ,I,J) DO 3 K=1,3 DO 3 L=1,3 3 SA1=SA1+EPS(I,K,L)*ALPHA(IQ,I,J)*A(IQ,K,L,J)/3.0d0 c c IL+IR, backscattering, DCP_I c SDCP180=3.0d0*SAL1-SAL0 c c IL+/-IR, backscattering S180=7.0d0*SAL1+SAL0 D180=8.0d0*(3.0d0*SAG1-SAG0+SA1) D=D180 IF(S180.gt.1.0d-9)D=D/S180 c c IL+/-IR, 90o, x S90X=7.0d0*SAL1+SAL0 D90X=2.0d0*(7.0d0*SAG1+SAG0+SA1) DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c IL+/-IR, 90o, z S90Z=2.0d0*(3.0d0*SAL1-SAL0) D90Z=4.0d0*(3.0d0*SAG1-SAG0-SA1) DZ=D90Z IF(S90Z.gt.1.0d-9)DZ=DZ/S90Z c c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 P1=YDY/YDX c c Try to get Gaussian units: AMU=1822.0d0 BOHR=0.529177d0 SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) c gpisvejc=980.0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c c c beta(A)^2=Delta2 in Gaussian: betaa=SA1*3.0d0/2.0d0*gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0*gpisvejc c c DCP 180 c ram3=24.0d0*beta2 roa3=8.0d0*(12.0d0*betag+4.0d0*betaa) c c unknown: S90X=S90X*(AMU*BOHR**4) D90X=D90x*gpisvejc DX=D90X IF(S90X.gt.1.0d-9)DX=DX/S90X c c ICP(90)=90Z roa2=8.0d0*(3.0d0*betag-betaa) ram2=12.0d0*beta2 CID2=0.0d0 if(ram2.gt.1.0d-9)CID2=roa2/ram2 c c ICP(180): c 8/c*(12*beta(G')^2 + 4*beta(A)^2): roa1=96.0d0*(betag+betaa/3.0d0) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) CID1=0.0d0 if(ram1.gt.1.0d-9)CID1=roa1/ram1 if(lopt(16))WRITE(4,3300)IQ,ECM,YDX,YDY, 1 ram1,CID2,DX,CID1,P1,roa1,roa3 3300 FORMAT(I4,f9.2,6f9.2,f6.3,2f11.2) c endif endif ENDIF 1 CONTINUE if(lopt(16))WRITE(4,3002) if(lopt(16))CLOSE(4) RETURN END SUBROUTINE TRTEN(NS,S,N,ALPHA,A,G) IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) PARAMETER (NA0=1000,MX3=3*NA0) DIMENSION S(NS,N),A(NS,3,3,3), 1ALPHA(NS,3,3),G(NS,3,3),V(MX3) logical lopt common/opts/lopt(300) SFAC=0.0234280d0 DO 1 I=1,3 DO 1 J=1,3 DO 3 JQ=1,N 3 V(JQ)=ALPHA(JQ,I,J) DO 1 IQ=1,N ALPHAS=0.0d0 DO 2 IX=1,N 2 ALPHAS=ALPHAS+S(IX,IQ)*V(IX)*SFAC 1 ALPHA(IQ,I,J)=ALPHAS DO 11 I=1,3 DO 11 J=1,3 DO 31 JQ=1,N 31 V(JQ)=G(JQ,I,J) DO 11 IQ=1,N GS=0.0d0 DO 21 IX=1,N 21 GS=GS+S(IX,IQ)*V(IX)*SFAC 11 G(IQ,I,J)=GS DO 12 I=1,3 DO 12 J=1,3 DO 12 K=1,3 DO 32 JQ=1,N 32 V(JQ)=A(JQ,I,J,K) DO 12 IQ=1,N AS=0.0D0 DO 22 IX=1,N 22 AS=AS+S(IX,IQ)*V(IX)*SFAC 12 A(IQ,I,J,K)=AS if(lopt(17))WRITE(6,3001) 3001 FORMAT(1X,'TENSORS HAVE BEEN TRANSFORMED TO Q COORDINATES',/, 11X,'ALPHA DERIVATIVES LISTED:') DO 4 IQ=1,N SP=ALPHA(IQ,1,1)+ALPHA(IQ,2,2)+ALPHA(IQ,3,3) 4 if(lopt(17))WRITE(6,3000)IQ,SP,((ALPHA(IQ,I,J),I=1,3),J=1,3) 3000 FORMAT(1X,I4,' MODE; Trace: ',f15.8,/,3(3F15.8,/),/) if(lopt(17))WRITE(6,3002) 3002 FORMAT(1X, 11X,'G (magn index -->)') DO 5 IQ=1,N SP=G(IQ,1,1)+G(IQ,2,2)+G(IQ,3,3) 5 if(lopt(17))WRITE(6,3000)IQ,SP,((G(IQ,J,I),I=1,3),J=1,3) if(lopt(17))WRITE(6,3003) 3003 FORMAT(1X, 11X,'A DERIVATIVES LISTED:') DO 6 IQ=1,N 6 if(lopt(17))WRITE(6,3004)IQ,(((A(IQ,I,J,K),I=1,3),J=1,3),K=1,3) 3004 FORMAT(1X,I4,' MODE',/,3(9F8.5,/),/) RETURN END c ==================== SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R) c C Transforms A and G to local origins (ICO=1) or common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in AU C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=600) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J) C IF(ICO.EQ.1)S=1.0d0 IF(ICO.EQ.2)S=-1.0d0 IF(ICO.EQ.3)S=0.0d0 C DO 1 I=1,3 DO 1 J=1,3 DO 1 L=1,NAT DO 1 IE=1,3 II=3*(L-1)+IE DO 1 K=1,3 DO 1 M=1,3 c G-second index-magnetic 1 G(II,I,J)=G(II,I,J)*S*S+S*0.5d0*EPS(J,K,M)*X(K,L)*ALPHA(II,I,M) IF(ICO.EQ.1)WRITE(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' G set to zero' C DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 II=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(II,IA,ID) ENDIF A(II,IA,IB,IC)=A(II,IA,IB,IC)+SUM*S 1-1.5d0*S*(X(IB,L)*ALPHA(II,IA,IC)+X(IC,L)*ALPHA(II,IA,IB)) IF(ICO.EQ.3)A(II,IA,IB,IC)=0.0d0 3 CONTINUE c IF(ICO.EQ.1)WRITE(6,*)' A,G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' A,G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A,G set to zero' RETURN END