program mvcd c magnetic vibrational circular dichroism implicit none real*8 CM,srt,af1,af2,unj(3),unk(3),mjk(3),mnkj(3),ds,b1,b2,AUD, 1AAf,Ej,AF,CLIGHT,del,af3,b3 integer*4 i,k,J,NM,ix,N,N7,nat,IA character*20 s20 real*8,allocatable::QP(:,:),QM(:,:),aptaq(:,:),aataq(:,:),S(:,:), 1E(:),aatapt(:,:) srt=dsqrt(2.0d0) AUD=2.541732d0 CM=219470.0d0 CLIGHT=137.053d0 c energy "uncertainty", cm-1: if(iargc().ne.0)then call getarg(1,s20) read(s20,*)del del=del/CM else del=1.0d-2/CM endif write(6,6001) 6001 format(/,' MVCD from anharmonic parameters',/,/, 1 ' Input: TENAQ.TXT.SCR',/, 1 ' TENQ.TXT.SCR',/, 1 ' F.INP',/,/, 1 ' Output: MVCD.TAB',/) open(9,file='F.INP',status='old') read(9,*)NM,N close(9) N7=N-NM+1 allocate(QP(3,N),QM(3,N),aptaq(2*N**2,3),aataq(N**2,3),S(N,N), 1E(N)) c read vibrational frequencies: call READS(S,N,E) c read dipole derivatives: call rdten1q(QP,QM,N,N) c read second dipole derivatives: call rdtenq(aptaq,aataq,N,N) c AAT derivatives through APT model: nat=N/3 allocate(aatapt(N**2,3)) do 10 I=N7,N do 10 J=N7,N aatapt(I+(J-1)*N,1)=0.0d0 aatapt(I+(J-1)*N,2)=0.0d0 aatapt(I+(J-1)*N,3)=0.0d0 do 10 IA=1,nat aatapt(I+(J-1)*N,1)= 10.5d0*(S(3*(IA-1)+2,J)*QP(3,I)-S(3*(IA-1)+3,J)*QP(2,I)) aatapt(I+(J-1)*N,2)= 10.5d0*(S(3*(IA-1)+3,J)*QP(1,I)-S(3*(IA-1)+1,J)*QP(3,I)) 10 aatapt(I+(J-1)*N,3)= 10.5d0*(S(3*(IA-1)+1,J)*QP(2,I)-S(3*(IA-1)+2,J)*QP(1,I)) open(87,file='compar.txt') do 9 I=N7,N do 9 J=N7,N do 9 ix=1,3 9 write(87,800)I,J,ix,E(I),E(J), 1aatapt(I+(J-1)*N,ix),aataq(I+(J-1)*N,ix) 800 format(2I5,i2,2F10.2,2G12.4) close(87) open(87,file='MVCD.TAB') open(88,file='MVCD-1.TAB') open(89,file='MVCD-2.TAB') open(90,file='MVCD-3.TAB') do 5 i=87,90 5 write(i,8787)del*CM 8787 format(' MVCD',/,' transition w D DE: ',g12.4,' cm-1'/,60(1h-)) do 1 j=N,N7,-1 do 2 ix=1,3 c electric dipole =P/sqrt(2) 2 unj(ix)=QP(ix,j)/srt ds=(unj(1)**2+unj(2)**2+unj(3)**2)*AUD**2 b1=0.0d0 b2=0.0d0 b3=0.0d0 Ej=E(j) do 3 k=N7,N do 8 ix=1,3 c electric dipole =P/sqrt(2) 8 unk(ix)=QP(ix,k)/srt if(k.ne.j)then do 4 ix=1,3 c magnetic dipole =M(2)/2 mjk(ix)=aataq(K+(J-1)*N,ix)-aataq(J+(K-1)*N,ix) c magnetic dipole =0 4 mnkj(ix)=aataq(K+(J-1)*N,ix)+aataq(J+(K-1)*N,ix) c -1/(-Ej - Ek) af1=AAf(j,k,E(j)+E(k),0.0d0,del) b1=b1 1 +(mnkj(1)*(unj(2)*unk(3)-unj(3)*unk(2)) 1 + mnkj(2)*(unj(3)*unk(1)-unj(1)*unk(3)) 1 + mnkj(3)*(unj(1)*unk(2)-unj(2)*unk(1)))*af1/del**2 c 1/(Ej - Ek) af2=AAf(j,k,E(j),E(k),del) b2=b2 1 +(mjk( 1)*(unj(2)*unk(3)-unj(3)*unk(2)) 1 + mjk( 2)*(unj(3)*unk(1)-unj(1)*unk(3)) 1 + mjk( 3)*(unj(1)*unk(2)-unj(2)*unk(1)))*af2/del**2 Ej=Ej+af2 else c -1/(-Ej - Ej) af3=AAf(j,k,E(j)+E(j),0.0d0,del) do 7 ix=1,3 c magnetic dipole <2|M(2)PjQj|0>=M(2)PjQj 7 mnkj(ix)=aataq(J+(J-1)*N,ix) b3=b3 1 +(mnkj(1)*(unj(2)*unk(3)-unj(3)*unk(2)) 1 + mnkj(2)*(unj(3)*unk(1)-unj(1)*unk(3)) 1 + mnkj(3)*(unj(1)*unk(2)-unj(2)*unk(1)))*af3/del**2 endif 3 continue c magn. moment = MSsqrt(w) x 2/c (all au): AF=2.0d0/CLIGHT c Deff = 2 Bterm B -> Reff = Bterm B / 2: AF=AF/2.0D0 c make mB/Ejk dimensionless, normalize to one Tesla c B(Tesla)=B(au)*2.5xe5 AF=AF/2.5D5 c from the <0|PQ|1> magnetic element: AF=AF/2.0D0 c D^2 as for dipole strength: AF=AF*AUD**2 c pisvejc: AF=AF*CLIGHT write(87,8789)j,Ej*CM,ds,(b1+b2+b3)*AF,b1*AF,b2*AF,b3*AF write(88,8789)j,Ej*CM,ds,b1*AF write(89,8789)j,Ej*CM,ds,b2*AF 1 write(90,8789)j,Ej*CM,ds,b3*AF 8789 format(i5,f14.4,5g12.4) do 6 i=87,90 write(i,8788) 8788 format(60(1h-)) 6 close(i) write(6,*)'Done' end c ============================================================ subroutine rdten1q(QP,QM,N,NT3) implicit none integer*4 N,NT3,i,ix,ndif,j,iq real*8 QP(3,NT3),QM(3,NT3) do 1 i=1,N do 1 ix=1,3 QP(ix,i)=0.0d0 1 QM(ix,i)=0.0d0 open(9,file='TENQ.TXT.SCR') read(9,*) read(9,*) read(9,*)ndif DO 16 I=ndif,1,-1 do 16 j=1,3 read(9,*)iq,iq,ix,QP(ix,iq),QP(ix,iq),QP(ix,iq) 16 continue read(9,*) DO 17 I=ndif,1,-1 do 17 j=1,3 17 read(9,*)iq,iq,ix,QM(ix,iq) close(9) return end c ============================================================ subroutine rdtenq(aptaq,aataq,N,NT3) implicit none integer*4 N,NT3,i,ix,ndif,J,k,iq,jq real*8 aptaq(2*NT3**2,3),aataq(NT3**2,3) do 2 i=1,N**2 do 2 ix=1,3 aptaq(i,ix)=0.0d0 aptaq(N**2+i,ix)=0.0d0 2 aataq(i,ix)=0.0d0 open(9,file='TENAQ.TXT.SCR',status='old') read(9,*)ndif read(9,*) read(9,*) DO 18 I=ndif,1,-1 DO 18 J=I,1,-1 do 18 k=1,3 18 read(9,*)iq,jq,ix,aptaq(iq+(jq-1)*N,ix),aptaq(jq+(iq-1)*N,ix), 1aptaq(N**2+iq+(jq-1)*N,ix) read(9,*) read(9,*) DO 1 I=ndif,1,-1 DO 1 J=I,1,-1 do 1 k=1,3 1 read(9,*)iq,jq,ix,aataq(iq+(jq-1)*N,ix),aataq(jq+(iq-1)*N,ix) close(9) return end c ============================================================ function AAf(j,k,enp,en,v) c for non-deg case goes to v^2/(enp-en) else to +/- |V| implicit none real*8 AAf,en,enp,v integer*4 j,k if(en.gt.enp)then AAf=-0.5d0*(enp-en+dsqrt((en-enp)**2+4.0d0*v**2)) else if(en.lt.enp)then AAf=-0.5d0*(enp-en-dsqrt((en-enp)**2+4.0d0*v**2)) else if(k.gt.j)then AAf=-0.5d0*(enp-en-dsqrt((en-enp)**2+4.0d0*v**2)) else AAf=-0.5d0*(enp-en+dsqrt((en-enp)**2+4.0d0*v**2)) endif endif endif return end c ============================================================ SUBROUTINE READS(S,N,E) implicit none integer*4 N,NM,nat,I,IA,K,iq real*8 S(N,N),E(*),cm,SFACR cm=219470.0d0 SFACR=0.02342179D0 OPEN(7,FILE='F.INP') READ(7,*)NM,N,nat DO 2 I=1,N E(I)=0.0d0 DO 2 IA=1,N 2 S(IA,I)=0.0d0 DO 1 I=1,nat 1 READ(7,*) READ(7,*) DO 3 IA=1,NAT DO 3 I=1,NM 3 READ(7,*)(S(3*(IA-1)+K,N-I+1),K=1,2),(S(3*(IA-1)+K,N-I+1),K=1,3) READ(7,*) READ(7,*)(E(N-I+1),I=1,NM) CLOSE(7) do 9 iq=1,nm E(N-iq+1)=E(N-iq+1)/cm do 9 i=1,nat do 9 k=1,3 9 S(3*(i-1)+k,N-iq+1)=S(3*(i-1)+k,N-iq+1)*SFACR RETURN END