program ssd c supporting program for toluen ch3 rotation c normal mode phase correction during the scan IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) parameter (N0=100,nf0=150,nw0=10) DIMENSION APT(nf0,N0*3,3),AAT(nf0,N0*3,3),wc(nw0),ws(nw0), 1AL(nf0,N0*3,3,3),G(nf0,N0*3,3,3),A(nf0,N0*3,3,3,3), 2APTu(nw0,N0*3,3),AATu(nw0,N0*3,3),w(nf0,N0*3),wu(nw0,3*N0), 3ALu(nw0,N0*3,3,3),Gu(nw0,N0*3,3,3),Au(nw0,N0*3,3,3,3) character*80 fn1(nf0),fn2(nf0) character*9 s9,num Z0=0.0d0 pi=4.0d0*atan(1.0d0) open(9,file='SSD.OPT') read(9,*)ncos if(ncos.gt.nw0)call report('too many cosinus freqs') read(9,*)(wc(i),i=1,ncos) read(9,*)nsin if(nsin+ncos.gt.nw0)call report('too many cos and sinus freqs') read(9,*)(ws(i),i=1,nsin) read(9,*)fis,dfi fis=fis*pi/180.0d0 dfi=dfi*pi/180.0d0 read(9,*)nf if(nf.gt.nf0)call report('too many angles') fie=fis+dble(nf-1)*dfi do 1 i=1,nf write(num,9001)i 9001 format(i9) c c read frequencies open(4,file=num//'/F.INP') call READS(nf0,3*N0,i,w) close(4) c c read normal mode dipole derivatives fn1(i)=num//'/FILE.Q.A.TEN' open(15,file=fn1(i)) read(15,*)NAT if(NAT.gt.N0)call report('too many atoms') DO 12 L=1,NAT*3 12 read(15,*) (APT(i,L,II),II=1,3) DO 221 L=1,NAT*3 221 read(15,*) (AAT(i,L,II),II=1,3) CLOSE(15) c read polarizability derivatives fn2(i)=num//'/FILE.Q.A.TTT' open(2,file=fn2(i)) READ(2,*) READ(2,*)NAT if(NAT.gt.N0)call report('too many atoms') READ(2,*) READ(2,*) DO 15 II=1,3 READ(2,*) DO 15 IND=1,NAT*3 15 READ(2,*)LDUM,IXDUM,(AL(i,IND,II,J),J=1,3) READ(2,*) READ(2,*) DO 13 II=1,3 READ(2,*) DO 13 IND=1,NAT*3 13 READ(2,*)LDUM,IXDUM,(G(i,IND,II,J),J=1,3) READ(2,*) READ(2,*) DO 14 II=1,3 DO 14 J=1,3 READ(2,*) DO 14 IND=1,NAT*3 14 READ(2,*)LDUM,IXDUM,(A(i,IND,II,J,K),K=1,3) 1 CLOSE(2) do 4 IND=1,3*NAT c c first three angles, investigate +++,++-,+-+,+-- for mode IND: s1=0.0d0 s2=0.0d0 s3=0.0d0 s4=0.0d0 do 3 K=1,3 do 3 J=1,3 s1=s1+(AL(1,IND,K,J)+AL(3,IND,K,J)-2.0d0*AL(2,IND,K,J))**2 s2=s2+(AL(1,IND,K,J)-AL(3,IND,K,J)-2.0d0*AL(2,IND,K,J))**2 s3=s3+(AL(1,IND,K,J)+AL(3,IND,K,J)+2.0d0*AL(2,IND,K,J))**2 3 s4=s4+(AL(1,IND,K,J)-AL(3,IND,K,J)+2.0d0*AL(2,IND,K,J))**2 if( s1.le.s2.and.s1.le.s3.and.s1.le.s4)ic=1 if(s2.le.s1.and. s2.le.s3.and.s2.le.s4)ic=2 if(s3.le.s1.and.s3.le.s2.and. s3.le.s4)ic=3 if(s4.le.s1.and.s4.le.s2.and.s4.le.s3 )ic=4 c switch the phase of mode IND in the third angle if(ic.eq.2.or.ic.eq.4)call ss(3,IND,G,AL,A,AAT,APT,N0,nf0) c switch the phase of mode IND in the second angle if(ic.eq.3.or.ic.eq.4)call ss(2,IND,G,AL,A,AAT,APT,N0,nf0) c c rest of angles, set phase according to the previous two: do 4 l=4,nf l0=l-1 l1=l-2 s1=0.0d0 s2=0.0d0 do 8 II=1,3 do 8 J=1,3 s1=s1+(AL(l1,IND,II,J)+AL(l,IND,II,J)-2.0d0*AL(l0,IND,II,J))**2 8 s2=s2+(AL(l1,IND,II,J)-AL(l,IND,II,J)-2.0d0*AL(l0,IND,II,J))**2 c switch the phase of mode IND in the l^th angle 4 if(s2.lt.s1)call ss(l,IND,G,AL,A,AAT,APT,N0,nf0) c c Initialize Fourier expansions do 10 i=1,ncos+nsin DO 10 L=1,NAT*3 wu(i,L)=0.0d0 do 10 II=1,3 APTu(i,L,II)=0.0d0 AATu(i,L,II)=0.0d0 do 10 J=1,3 ALu(i,L,II,J)=0.0d0 Gu(i,L,II,J)=0.0d0 do 10 K=1,3 10 Au(i,L,II,J,K)=0.0d0 c c save phase-corrected tensors fi=fis-dfi do 9 i=1,nf fi=fi+dfi c c collect Fourier expansions c start and end points extra: if(i.eq.1.or.i.eq.nf)then dd=dfi/2.0d0 else dd=dfi endif do 102 in=1,ncos+nsin if(in.eq.1.and.abs(wc(in)).lt.1.0d-4)then ff=dd/(fie-fis) else ff=2.0d0*dd/(fie-fis) endif if(in.le.ncos)then cs=ff*cos(wc(in)*fi) else cs=ff*sin(ws(in-ncos)*fi) endif DO 102 L=1,NAT*3 wu(in,L)=wu(in,L)+w(i,L)*cs do 102 II=1,3 APTu(in,L,II)=APTu(in,L,II)+APT(i,L,II)*cs AATu(in,L,II)=AATu(in,L,II)+AAT(i,L,II)*cs do 102 J=1,3 ALu(in,L,II,J)=ALu(in,L,II,J)+AL(i,L,II,J)*cs Gu(in,L,II,J)=Gu(in,L,II,J)+G(i,L,II,J)*cs do 102 K=1,3 102 Au(in,L,II,J,K)=Au(in,L,II,J,K)+A(i,L,II,J,K)*cs do 6 ie=80,1,-1 6 if(fn1(i)(ie:ie).ne.' ')goto 66 66 open(15,file=fn1(i)(1:ie)//'_') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5,' normal derivatives, phase corrected') DO 20 L=1,NAT*3 20 WRITE(15,1501) (APT(i,L,II),II=1,3),L 1501 FORMAT(3F14.8,I5) DO 224 L=1,NAT*3 224 WRITE(15,1501) (AAT(i,L,II),II=1,3),L DO 2310 L=1,NAT*3 2310 WRITE(15,1501) (Z0,II=1,3),L DO 101 L=1,NAT*3 101 WRITE(15,1501) (APT(i,L,II),II=1,3),L CLOSE(15) do 5 ie=80,1,-1 5 if(fn2(i)(ie:ie).ne.' ')goto 55 55 open(22,file=fn2(i)(1:ie)//'_') WRITE(22,2001)NAT 2001 FORMAT(' ROA tensors, normal modes derivatives phase corrected',/, 1I4/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') CALL WRITETTT(N0*3,NAT,AL,A,G,i,nf0) 9 close(22) c save the transforms to disk do 103 in=1,ncos+nsin write(s9,999)in 999 format(i9) open(15,file=s9//'.u.ten') WRITE(15,1504) NAT,NAT-6,0,in 1504 FORMAT(3I5,' normal derivatives, fourier',i3) DO 104 L=1,NAT*3 104 WRITE(15,1501) (APTu(in,L,II),II=1,3),L DO 2244 L=1,NAT*3 2244 WRITE(15,1501) (AATu(in,L,II),II=1,3),L DO 2311 L=1,NAT*3 2311 WRITE(15,1501) (Z0,II=1,3),L DO 105 L=1,NAT*3 105 WRITE(15,1501) (APTu(in,L,II),II=1,3),L CLOSE(15) open(22,file=s9//'.u.ttt') WRITE(22,2002)in,NAT 2002 FORMAT(' ROA tensors, fourier',i3/, 1I4/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') CALL WRITETTT(N0*3,NAT,ALu,Au,Gu,in,nw0) 103 close(22) open(22,file='u.finp') WRITE(22,2003) 2003 FORMAT(' Frequency Fourier expansions:') do 106 L=1,3*NAT write(22,2255)L 106 write(22,2256)(wu(in,L),in=1,ncos+nsin) 2255 format(i4,$) 2256 format(f8.2,20f8.3) close(22) stop end SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END SUBROUTINE WRITETTT(N30,NAT,AL,A,G,io,nf0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION AL(nf0,N30,3,3),G(nf0,N30,3,3),A(nf0,N30,3,3,3) c DO 1 I=1,3 WRITE(22,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 II=3*(L-1)+IX 1 WRITE(22,2001)L,IX,(AL(io,II,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F13.7) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(22,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 II=3*(L-1)+IX c last inner index magnetic 2 WRITE(22,2001)L,IX,(G(io,II,I,J),J=1,3) WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(22,2006)I,J 2006 FORMAT(' A(',I1,'(u),',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 II=3*(L-1)+IX 3 WRITE(22,2001)L,IX,(A(io,II,I,J,K),K=1,3) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' Atom/x vx vy vz') DO 11 I=1,3 WRITE(22,2012)I 2012 FORMAT(' AlphaV(',I1,',J):') DO 11 L=1,NAT DO 11 IX=1,3 II=3*(L-1)+IX 11 WRITE(22,2001)L,IX,(AL(io,II,I,J),J=1,3) RETURN END subroutine ss(ia,IND,G,AL,A,AAT,APT,N0,nf0) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION APT(nf0,N0*3,3),AAT(nf0,N0*3,3) DIMENSION AL(nf0,N0*3,3,3),G(nf0,N0*3,3,3),A(nf0,N0*3,3,3,3) do 5 II=1,3 APT(ia,IND,II)=-APT(ia,IND,II) AAT(ia,IND,II)=-AAT(ia,IND,II) do 5 J=1,3 G(ia,IND,II,J)=-G(ia,IND,II,J) AL(ia,IND,II,J)=-AL(ia,IND,II,J) do 5 K=1,3 5 A(ia,IND,II,J,K)=-A(ia,IND,II,J,K) write(6,90)ia,IND 90 format('Angle ',i3,': phase of mode ',i4,' switched') return end subroutine READS(nf0,N3,iu,w) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION W(nf0,N3) read(4,*)nint,N if(nint.gt.N3)call report('nint>N3') NAT=N/3 N7=N-NINT+1 do 1 i=1,NAT 1 read(4,*) read(4,*) DO 2 I=1,NAT DO 2 J=1,nint 2 read(4,*) read(4,*) READ(4,4000)(w(iu,I),I=N7,N) 4000 FORMAT(6F11.6) DO 7 I=1,N7-1 7 w(iu,I)=0.0d0 RETURN end