program ssdn c supporting program for molecules with (two) ch3 rotations c normal mode phase correction during the scan IMPLICIT none integer*4 N0,nf0,nw0,nd0 parameter (N0=100,nf0=1500,nw0=10,nd0=2) real*8 w(nf0,N0*3),APT(nf0,N0*3,3),AAT(nf0,N0*3,3), 1AL(nf0,N0*3,3,3),G(nf0,N0*3,3,3),A(nf0,N0*3,3,3,3), 1wc(nd0,nw0),ws(nd0,nw0), 2APTu(nw0**nd0,N0*3,3),AATu(nw0**nd0,N0*3,3), 3wu(nw0**nd0,3*N0),Au(nw0**nd0,N0*3,3,3,3), 3ALu(nw0**nd0,N0*3,3,3),Gu(nw0**nd0,N0*3,3,3), 1fi(nd0),fis(nd0),fie(nd0),dfi(nd0) integer*4 ncos(nd0),nsin(nd0),ir,i,nf,NAT,L,II,IND,J,K, 1idum,ic,l0,l1,ix,in,jn,injn,ndf,ie character*80 fn1(nf0),fn2(nf0) character*9 s9,num real*8 Z0,pi,s1,s2,s3,s4,dd,a1,cs Z0=0.0d0 pi=4.0d0*atan(1.0d0) open(9,file='SSDN.OPT') read(9,*)ndf if(ndf.gt.2)call report('too many rotations') do 100 ir=1,ndf read(9,*)ncos(ir) if(ncos(ir).gt.nw0)call report('too many cosinus freqs') read(9,*)(wc(ir,i),i=1,ncos(ir)) read(9,*)nsin(ir) if(nsin(ir)+ncos(ir).gt.nw0)call report('too many c and s freqs') read(9,*)(ws(ir,i),i=1,nsin(ir)) 100 read(9,*)fis(ir),fie(ir),dfi(ir) read(9,*)nf if(nf.gt.nf0)call report('too many points') c c correct the phase factors: 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,*)IDUM,IDUM,(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,*)IDUM,IDUM,(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,*)IDUM,IDUM,(A(i,IND,II,J,K),K=1,3) 1 CLOSE(2) open(88,file='SW.TXT.SCR') 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) close(88) c c Initialize Fourier expansions do 10 i=1,(ncos(1)+nsin(1))*(ncos(2)+nsin(2)) 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 c and collect Fourier expansions open(9,file='LIST.TXT') open(61,file='AL11.TXT') do 9 i=1,nf read(9,*)(fi(ix),ix=1,ndf) c c angle ranges must agree to ensure orthonormality, otherwise c all is screwed up!: do 91 ix=1,ndf if(fi(ix).lt.fis(ix))call report('error in angle ranges') 91 if(fi(ix).gt.fie(ix))call report('error in angle ranges') c c correct periodicity, set integration intervals c start and end points extra: dd=1.0d0 do 93 ix=1,ndf dd=dd*dabs(dfi(ix)/(fie(ix)-fis(ix))) if(dabs(fi(ix)-fis(ix)).lt.dabs(0.2d0*dfi(ix)))dd=dd/2.0d0 93 if(dabs(fi(ix)-fie(ix)).lt.dabs(0.2d0*dfi(ix)))dd=dd/2.0d0 do 92 ix=1,ndf 92 fi(ix)=fi(ix)*pi/180.0d0 do 102 in=1,ncos(1)+nsin(1) if(in.le.ncos(1))then a1=cos(wc(1,in )*fi(1))*2.0d0*dd if(in.eq.1.and.abs(wc(1,in)).lt.1.0d-4)a1=a1/2.0d0 else a1=sin(ws(1,in-ncos(1))*fi(1))*2.0d0*dd endif do 102 jn=1,ncos(2)+nsin(2) injn=in+(jn-1)*(ncos(1)+nsin(1)) if(jn.le.ncos(2))then cs=a1*cos(wc(2,jn )*fi(2))*2.0d0 if(jn.eq.1.and.abs(wc(2,jn)).lt.1.0d-4)cs=cs/2.0d0 else cs=a1*sin(ws(2,jn-ncos(2))*fi(2))*2.0d0 endif DO 102 L=1,NAT*3 wu(injn,L)=wu(injn,L)+w(i,L)*cs c if(injn.eq.1.and.L.eq.8)write(6,*)i,fi(1),fi(2),w(i,L),cs do 102 II=1,3 APTu(injn,L,II)=APTu(injn,L,II)+APT(i,L,II)*cs AATu(injn,L,II)=AATu(injn,L,II)+AAT(i,L,II)*cs do 102 J=1,3 c if(injn.eq.1.and.L.eq.1.and.II.eq.1.and.J.eq.1)write(6,*) c 1AL(i,L,II,J),cs c stop ALu(injn,L,II,J)=ALu(injn,L,II,J)+AL(i,L,II,J)*cs Gu(injn,L,II,J)=Gu(injn,L,II,J)+G(i,L,II,J)*cs do 102 K=1,3 102 Au(injn,L,II,J,K)=Au(injn,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) write(61,*)fi(1)*180.0d0/pi,fi(2)*180.0d0/pi,AL(i,1,1,1) 9 close(22) write(6,*)nf,' points' close(9) close(91) c save the transforms to disk do 103 in=1,ncos(1)+nsin(1) do 103 jn=1,ncos(2)+nsin(2) injn=in+(jn-1)*(ncos(1)+nsin(1)) write(s9,999)injn 999 format(i9) open(15,file=s9//'.u.ten') WRITE(15,1504) NAT,NAT-6,0,in,jn 1504 FORMAT(3I5,' normal derivatives, fourier',2i3) DO 104 L=1,NAT*3 104 WRITE(15,1501) (APTu(injn,L,II),II=1,3),L DO 2244 L=1,NAT*3 2244 WRITE(15,1501) (AATu(injn,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(injn,L,II),II=1,3),L CLOSE(15) open(22,file=s9//'.u.ttt') WRITE(22,2002)in,jn,NAT 2002 FORMAT(' ROA tensors, fourier',2i3/, 1I4/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') CALL WRITETTT(N0*3,NAT,ALu,Au,Gu,injn,nw0**nd0) 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 2255 format(i4) do 106 jn=1,ncos(2)+nsin(2) 106 write(22,2256)(wu(in+(jn-1)*(ncos(1)+nsin(1)),L), 1in=1,ncos(1)+nsin(1)) 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(88,90)ia,IND 90 format('Angle ',i6,': 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