PROGRAM NEW6 IMPLICIT none INTEGER*4 NAT,N,ia,ja,ii,ix,M,ma,k,im,jj,jx,ij,ji real*8 rai,wmax,wmin,neg,dij real*8,allocatable::SS(:,:),S(:,:),ALPHA(:,:,:),A(:,:,:,:), 1G(:,:,:),R(:),E(:),iram(:),RAM(:,:,:),PM(:,:),AM(:,:) integer*4,allocatable::Q(:),mram(:),nram(:) logical lpos character*5 sk(4) data sk/'IR ','VCD ','Raman','ROA '/ C WRITE(6,6666) 6666 FORMAT(/,/, A' * * * * * * * * * * * * * * * * * * * * * *',/, 1' * * ',/, 2' * VIBRATIONAL OPTICAL ACTIVITY * ',/, 3' * * ',/, 4' * * * * * * * * * * * * * * * * * * * * * *',/, 8' ',/, 8' ATOMIC CONTRIBUTION GRAPHICS ',/, 8' Petr Bour 2021 ',/,/,/, A'INPUT FILES: FILE.TTT,FILE.TEN,F.INP',/, B'OUTPUT : GRX.TXT',/) OPEN(2,FILE='FILE.TTT',STATUS='OLD') READ(2,*) READ(2,*)NAT close(2) N=3*NAT allocate(S(N,N),ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),R(N),E(N),Q(NAT), 1SS(N,N),RAM(4,N,N),PM(N,3),AM(N,3)) c RAM(1 ... IR c 2 ... VCD c 3 ... Raman c 4 ... ROA c CALL READTEN(N,NAT,ALPHA,G,A) call READDD(N,PM,AM) CALL READS(N,S,E,M,R,Q) write(6,609)NAT,M 609 format(i6,' atoms,',i6,' modes') CALL DORAMAN(N,PM,AM,ALPHA,G,A,5320.0d0,RAM) C 11 write(6,799) 799 format(' MA WMIN WMAX (0 0 to quit): ',$) read(5,*)ma,wmin,wmax lpos=ma.lt.0 ma=abs(ma) ma=min(ma,NAT*(NAT+1)/2) if(ma.eq.0)stop c SS=S.S call dss(N,M,SS,S,E,WMIN,WMAX) c for each atom pair, get intenzity, record maximum c number of maxima: open(9,file='GRX.TXT') write(9,904)wmin,wmax 904 format(2f12.2) allocate(mram(ma),nram(ma),iram(ma)) do 4 k=1,4 iram=0.0d0 nram=0 mram=0 do 1 ia=1,NAT do 1 ja=ia,NAT if(lpos)then c separate positive and negative contributions rai=0.0d0 neg=0.0d0 do 6 im=1,M dij=0.0d0 if(E(im).le.WMAX.and.E(im).GE.WMIN)then if(ia.ne.ja)then do 8 ix=1,3 ii=ix+3*(ia-1) ij=ix+3*(ja-1) do 8 jx=1,3 jj=jx+3*(ja-1) ji=jx+3*(ia-1) 8 dij=dij+RAM(k,ii,jj)*S(ii,im)*S(jj,im) 1 +RAM(k,ij,ji)*S(ij,im)*S(ji,im) else do 7 ix=1,3 ii=ix+3*(ia-1) do 7 jx=1,3 jj=jx+3*(ia-1) 7 dij=dij+RAM(k,ii,jj)*S(ii,im)*S(jj,im) endif endif if(dij.lt.0.0d0)then neg=neg+dij else rai=rai+dij endif 6 continue call rec(ma,ia,ja,rai,iram,nram,mram) call rec(ma,ia,ja,neg,iram,nram,mram) else c no separation rai=0.0d0 if(ia.ne.ja)then do 3 ix=1,3 ii=ix+3*(ia-1) ij=ix+3*(ja-1) do 3 jx=1,3 jj=jx+3*(ja-1) ji=jx+3*(ia-1) 3 rai=rai+RAM(k,ii,jj)*SS(ii,jj)+RAM(k,ij,ji)*SS(ij,ji) else do 2 ix=1,3 ii=ix+3*(ia-1) do 2 jx=1,3 jj=jx+3*(ja-1) 2 rai=rai+RAM(k,ii,jj)*SS(ii,jj) endif call rec(ma,ia,ja,rai,iram,nram,mram) endif 1 continue c normalize contributions: call cn(ma,iram) 4 call wr(sk(k),ma,nram,mram,iram) close(9) deallocate(mram,nram,iram) goto 11 END C ========================================== subroutine wr(sk,ma,nram,mram,iram) IMPLICIT none integer*4 ma,nram(ma),mram(ma),ir,i,ic real*8 iram(ma) character*5 sk write(9,905)sk write(6,905)sk 905 format(a5) ir=0 do 5 i=1,ma 5 if(iram(i).ne.0.0d0)ir=ir+1 write(9,*)ir ic=0 do 3 i=1,ma if(iram(i).ne.0.0d0)then ic=ic+1 write(9,601)nram(i),mram(i),iram(i) write(6,602)nram(i),mram(i),iram(i) 601 format(2i7,f12.6) 602 format(2i7,f12.6,$) if(mod(ic,3).eq.0)write(6,*) endif 3 continue if(mod(ic,3).ne.0)write(6,*) return end C ========================================== subroutine cn(ma,iram) IMPLICIT none integer*4 i,ma real*8 iram(ma),s s=0.0d0 do 1 i=1,ma 1 s=s+iram(i)**2 s=dsqrt(s) do 2 i=1,ma 2 iram(i)=iram(i)/s return end C ========================================== subroutine rec(ma,ia,ja,rai,iram,nram,mram) IMPLICIT none integer*4 ma,i,j,ia,ja,nram(ma),mram(ma) real*8 rai,iram(ma) do 3 i=1,ma if(dabs(rai).gt.dabs(iram(i)))then do 31 j=ma,i+1,-1 iram(j)=iram(j-1) nram(j)=nram(j-1) 31 mram(j)=mram(j-1) iram(i)=rai mram(i)=ia nram(i)=ja return endif 3 continue return end C ========================================== subroutine dss(N,M,SS,S,E,WMIN,WMAX) IMPLICIT none integer*4 N,M,i,j,k real*8 SS(N,N),S(N,N),E(N),WMIN,WMAX,ek SS=0.0d0 do 1 k=1,M ek=E(k) if(ek.le.WMAX.and.ek.GE.WMIN)then do 2 i=1,N do 2 j=i,N 2 SS(i,j)=SS(i,j)+S(i,k)*S(j,k) endif 1 continue do 3 i=1,N do 3 j=i+1,N 3 SS(j,i)=SS(i,j) return end C ============================================================== SUBROUTINE READS(N,S,E,M,R,Q) IMPLICIT none integer*4 Q(N/3),i,ix,N,M,j real*8 S(N,N),E(N),R(N) OPEN(4,FILE='F.INP',STATUS='OLD') read(4,*)M,N do 1 i=1,N/3 1 read(4,*)Q(i),(R(ix+3*(i-1)),ix=1,3) read(4,*) DO 2 I=1,N/3 DO 2 j=1,M 2 read(4,*)(s(3*(i-1)+ix,j),ix=1,2),(s(3*(i-1)+ix,j),ix=1,3) read(4,*) READ(4,4000)(E(I),I=1,M) 4000 FORMAT(6F11.6) close(4) RETURN end C SUBROUTINE READTEN(N,NAT,ALPHA,G,A) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,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,*)(ALPHA(IND,I,J),J=1,2),(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,*)(G(IND,I,J),J=1,2),(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,*)(A(IND,I,J,K),K=1,2),(A(IND,I,J,K),K=1,3) CLOSE(2) RETURN END C ============================================================== SUBROUTINE DORAMAN(N,PM,AM,ALPHA,G,A,EXCA,RAM) IMPLICIT none integer*4 IQ,I,JQ,J,N real*8 ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),EXCA,PM(N,3),AM(N,3), 1SAL0,SAL1,SAG0,SAG1,SA1,SPALI,SPALJ,dip,rot,S180,D180,D, 1AMU,BOHR,SpA3,beta2,gpisvejc,betaa,betag,RAM(4,N,N) AMU=1822.0d0 BOHR=0.529177d0 DO 1 IQ=1,N DO 1 JQ=1,N SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SPALI=0.0d0 SPALJ=0.0d0 SA1=0.0d0 DO 3 I=1,3 SPALI=SPALI+ALPHA(IQ,I,I) SPALJ=SPALJ+ALPHA(JQ,I,I) SA1=SA1+ALPHA(IQ,1,I)*(A(JQ,2,3,I)-A(JQ,3,2,I)) 1 +ALPHA(IQ,2,I)*(A(JQ,3,1,I)-A(JQ,1,3,I)) 1 +ALPHA(IQ,3,I)*(A(JQ,1,2,I)-A(JQ,2,1,I)) DO 3 J=1,3 SAL0=SAL0+ALPHA(IQ,I,I)*ALPHA(JQ,J,J) SAL1=SAL1+ALPHA(IQ,I,J)*ALPHA(JQ,I,J) SAG0=SAG0+ALPHA(IQ,I,I)* G(JQ,J,J) 3 SAG1=SAG1+ALPHA(IQ,I,J)* G(JQ,I,J) SA1=SA1/3.0d0 c 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 SpA3=SPALI*SPALJ/9.0d0*(AMU*BOHR**4) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA betaa=SA1*3.0d0/2.0d0 *gpisvejc betag=(3.0d0*SAG1-SAG0)/2.0d0 *gpisvejc dip=PM(IQ,1)*PM(JQ,1)+PM(IQ,2)*PM(JQ,2)+PM(IQ,3)*PM(JQ,3) rot=PM(IQ,1)*AM(JQ,1)+PM(IQ,2)*AM(JQ,2)+PM(IQ,3)*AM(JQ,3) RAM(1,IQ,JQ)=dip RAM(2,IQ,JQ)=rot RAM(3,IQ,JQ)=4.0d0*(45.0d0*SpA3+7.0d0*beta2) 1 RAM(4,IQ,JQ)=96.0d0*(betag+betaa/3.0d0) END SUBROUTINE READDD(N,P,A) C Read dipole derivatives: IMPLICIT none real*8 P(N,3),A(N,3) integer*4 N,NAT,L,J,I open(15,FILE='FILE.TEN',status='old') READ (15,*)NAT DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P((L-1)*3+J,I),I=1,3) DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(3*(L-1)+J,I),I=1,3) close(15) RETURN END C