PROGRAM SEROA IMPLICIT none integer*4 n,i,nat,ix,jx,kx,ii,istart,k,kk,l, 1j,jj,natm,ifi,it,iter,ns,natr,nf,ict integer*4 ,allocatable::n1(:),n2(:) real*8,allocatable:: alpha(:,:,:),gp(:,:,:),a(:,:,:,:), 1alphad(:,:,:),gpd(:,:,:),ad(:,:,:,:),r(:,:),P(:,:),X(:,:), 1alphar(:,:,:),gpr(:,:,:),ar(:,:,:,:),SIT(:,:), 1pd(:,:,:),temq(:,:),S(:,:),E(:),SI(:,:),SIOLD(:,:), 1tem(:,:),temi(:,:),pt(:,:),temk(:,:),teml(:,:),temp(:,:), 1alpharm(:,:,:),gprm(:,:,:),arm(:,:,:,:) real*8 r0,dt,df,om,th,fi,ds1,ds2,an1,an2,bohr,omold,ds3,tol character*80,allocatable:: fn(:) character*80 ft,s80 logical lex,lave,ldeb,luseg,lusea,lusep1,lusep2 character*10 s10 logical,allocatable:: ltr(:) c bohr=0.529177d0 ldeb=.false. luseg=.true. lusea=.true. lusep1=.true. lusep2=.true. tol=0.01d0 inquire(file='SEROA.OPT',exist=lex) if(.not.lex)call report('SEROA.OPT not found') c iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii open(9,file='SEROA.OPT') read(9,*)n write(6,*)n nat=0 allocate(r(n,3),fn(n),n1(n),n2(n)) c geometries,get centers do 1 i=1,n read(9,900)fn(i) write(6,*)i,fn(i) 900 format(a80) 1 CALL READRX(r,n,fn,i,nat,n1,n2) write(6,*)nat,'atoms' allocate(alphad(3*nat,3,3),gpd(3*nat,3,3),ad(3*nat,3,3,3), 1alpha(n,3,3),gp(n,3,3),a(n,3,3,3),ltr(n)) c polarizability derivatives do 2 i=1,n ltr(i)=.false. read(9,900)ft write(6,*)i,ft 2 CALL READTTT(nat,n,i,alphad,Gpd,ad,ft,fn,n1,n2) c polarizability alpha do 3 i=1,n read(9,900)ft 3 call readpol(alpha,i,ft,n) c G': do 4 i=1,n read(9,900)ft 4 call readpolg(gp,i,ft,n) c A: do 5 i=1,n read(9,900)ft 5 call readpola(a,i,ft,n) lave=.false. 66 read(9,920,end=888,err=888)s80 if(s80(1:4).eq.'AVER')read(9,*)lave if(s80(1:4).eq.'DEBU')read(9,*)ldeb if(s80(1:4).eq.'ADER')read(9,*)lusea if(s80(1:4).eq.'GDER')read(9,*)luseg if(s80(1:4).eq.'P1TE')read(9,*)lusep1 if(s80(1:4).eq.'P2TE')read(9,*)lusep2 if(s80(1:3).eq.'TOL')read(9,*)tol if(s80(1:4).eq.'TRAN')then do i=1,n read(9,*)ltr(i) enddo endif goto 66 920 format(a80) 888 close(9) c iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii write(6,*)'tensors read' c transform G, A to small molecule center call TTT1(alpha,a,gp,1,r,n,alphad,gpd,ad,nat,n1,n2,lusea,luseg, 1ltr) allocate(X(n*24,n*24),P(n*24,n*24),tem(n*24,n*24), 1temi(n*24,n*24),pt(n*24,n*24),temk(n*24,n*24), 1teml(n*24,n*24),temp(n*24,n*24),temq(n*24,n*24)) do 63 i=1,n*24 do 63 j=1,n*24 tem(i,j)=0.0d0 P(i,j)=0.0d0 63 X(i,j)=0.0d0 c ************* Calculate general distance tensor: ********* call doX(n,r,X) c ************* Calculate general polarizability tensor: ********* do 7 i=1,n istart=24*(i-1)+1 do 7 ix=1,3 ii=0 do 7 jx=1,3 P(istart+ix-1,istart+jx-1)=alpha(i,ix,jx) P(istart+ix+2,istart+jx+2)=alpha(i,ix,jx) P(istart+ix-1,istart+jx+8)=gp(i,ix,jx) P(istart+jx+8,istart+ix-1)=gp(i,ix,jx) P(istart+ix+2,istart+jx+5)=-gp(i,ix,jx) P(istart+jx+5,istart+ix+2)=-gp(i,ix,jx) do 7 kx=jx,3 ii=ii+1 P(istart+ix- 1,istart+ii+11)=a(i,ix,jx,kx)/3.0d0 P(istart+ii+11,istart+ix- 1)=a(i,ix,jx,kx)/3.0d0 P(istart+ix+ 2,istart+ii+17)=a(i,ix,jx,kx)/3.0d0 7 P(istart+ii+17,istart+ix+ 2)=a(i,ix,jx,kx)/3.0d0 c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 c ************* Calculate general polarizability derivative tensor: * allocate(pd(3*nat,24,24)) do 71 l=1,3*nat do 71 ix=1,3 ii=0 do 71 jx=1,3 pd(l,ix ,jx )=alphad(l,ix,jx) pd(l,ix+3,jx+3)=alphad(l,ix,jx) pd(l,ix ,jx+9)=gpd(l,ix,jx) pd(l,jx+9,ix )=gpd(l,ix,jx) pd(l,ix+3,jx+6)=-gpd(l,ix,jx) pd(l,jx+6,ix+3)=-gpd(l,ix,jx) do 71 kx=jx,3 ii=ii+1 pd(l,ix ,ii+12)=ad(l,ix,jx,kx)/3.0d0 pd(l,ii+12,ix )=ad(l,ix,jx,kx)/3.0d0 pd(l,ix+ 3,ii+18)=ad(l,ix,jx,kx)/3.0d0 71 pd(l,ii+18,ix+ 3)=ad(l,ix,jx,kx)/3.0d0 c get redressed polarizability derivatives, common origin allocate(alphar(3*nat,3,3),gpr(3*nat,3,3),ar(3*nat,3,3,3)) call MDP(x,p,n,tem,temi,temk,pt,alphar,gpr,ar,nat,n1,n2, 1pd,teml,temp,temq,r,lusep1,lusep2) call WRITETTT(NAT,alphaR,AR,GPR,'FILE.TTT') do 10 i=1, n if(fn(i)(1:4).ne.'ball')then natm=n2(i)-n1(i)+1 allocate(alpharm(3*natm,3,3),gprm(3*natm,3,3),arm(3*natm,3,3,3)) call rewrit(alpharm,gprm,arm,natm,NAT,alphar,gpr,ar,i,n1,n2) write(s10,100)i 100 format(i10) do 12 k=1,len(s10) 12 if(s10(k:k).ne.' ')goto 13 13 call WRITETTT(natm,alphaRm,ARm,GPRm,s10(k:10)//'r.ttt') if(i.eq.n)then allocate(S(3*natm,3*natm),E(3*natm)) call READS(S,E,3*natm,natr) if(natr.ne.natm)call report('Wrong atom numbers') call TRTEN(S,3*natm,alpharm,Arm,Gprm) allocate(SI(2,3*natm),SIOLD(2,3*natm),SIT(2,3*natm)) c Raman and ROA intensities: do 15 ii=1,3*natm do 15 k=1,2 SI(k,ii)=0.0d0 SIT(k,ii)=0.0d0 15 SIOLD(k,ii)=0.0d0 omold=0.0d0 call DORAMAN(3*natm,alpharm,Gprm,Arm,5320.0d0,E,SIT,0) ns=4 r0= 1 dsqrt((r(1,1)-r(2,1))**2+(r(1,2)-r(2,2))**2+(r(1,3)-r(2,3))**2) write(6,6001)r0*bohr,ns,tol 6001 format(' R12 = ',f18.2,' A',/, 1 ' Ns = ',i18,/, 1 ' tol = ',f18.3) if(lave)then write(6,*)'averaging' iter=0 99 iter=iter+1 dt= 4.0d0*atan(1.0d0)/dble(ns) do 151 ii=1,3*natm SI(1,ii)=0.0d0 151 SI(2,ii)=0.0d0 om=0.0d0 th=-dt/2.0d0 ict=0 do 14 it=1,ns th=th+dt r(1,3)=r0*cos(th) df=dt/sin(th) nf=max(1,nint(2.0d0*4.0d0*atan(1.0d0)/df)) df= 2.0d0*4.0d0*atan(1.0d0)/dble(nf) fi=-df/2.0d0 do 14 ifi=1,nf ict=ict+1 fi=fi+df if(ldeb)write(6,*)ict,it,th,ifi,fi r(1,1)=r0*sin(th)*cos(fi) r(1,2)=r0*sin(th)*sin(fi) call doX(n,r,X) call MDP(x,p,n,tem,temi,temk,pt,alphar,gpr,ar,nat,n1,n2, 1 pd,teml,temp,temq,r,lusep1,lusep2) call rewrit(alpharm,gprm,arm,natm,NAT,alphar,gpr,ar,i,n1,n2) call TRTEN(S,3*natm,alpharm,Arm,Gprm) call DORAMAN(3*natm,alpharm,Gprm,Arm,5320.0d0,E,SIT,1) if(ldeb) 1 call DORAMAN(3*natm,alpharm,Gprm,Arm,5320.0d0,E,SIT,ict+2) om=om+dt*sin(th)*df do 14 jj=1,2 do 14 kk=1,3*natm 14 SI(jj,kk)=SI(jj,kk)+SIT(jj,kk)*dt*sin(th)*df if(ldeb)call report('debug stop') c errors of Raman and ROA spectra: ds1=0.0d0 ds2=0.0d0 an1=0.0d0 an2=0.0d0 do 16 kk=1,3*natm SI(1,kk)=SI(1,kk)/om SI(2,kk)=SI(2,kk)/om an1=an1+dabs(SI(1,kk)) an2=an2+dabs(SI(2,kk)) ds1=ds1+dabs(SI(1,kk)-SIOLD(1,kk)) ds2=ds2+dabs(SI(2,kk)-SIOLD(2,kk)) SIOLD(1,kk)=SI(1,kk) 16 SIOLD(2,kk)=SI(2,kk) ds3=dabs(om-omold)/om omold=om if(an1.ne.0.0d0)ds1=ds1/an1 if(an2.ne.0.0d0)ds2=ds2/an2 write(6,6002)iter,ns,ds1,ds2,ds3,om/4.0d0/4.0d0/atan(1.0d0) 6002 format(2i5,4f15.4) if(ds1.gt.tol.or.ds2.gt.tol)then ns=ns*2 if(ds3.lt.1.0d-4)then write(6,*)'not converged' call DORAMAN(3*natm,alpharm,Gprm,Arm,5320.0d0,E,SI,2) else goto 99 endif else write(6,*)'converged' call DORAMAN(3*natm,alpharm,Gprm,Arm,5320.0d0,E,SI,2) endif else write(6,*)'averaging switched off' endif endif deallocate(alpharm,gprm,arm) endif 10 continue end c ============================================================ SUBROUTINE readpol(a,im,fn,n) implicit none integer*4 im,n,i,j character*(*)fn real*8 a(n,3,3) open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(im,i,j),i=1,j),j=1,3) close(90) do 1 j=1,3 do 1 i=1,j 1 a(im,j,i)=a(im,i,j) return end c ============================================================ SUBROUTINE readpolg(a,im,fn,n) implicit none integer*4 im,n,i,j character*(*)fn real*8 a(n,3,3) open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(im,j,i),i=1,3),j=1,3) close(90) return end c ============================================================ SUBROUTINE readpola(a,im,fn,n) implicit none integer*4 im,n,i,j,k character*(*)fn real*8 a(n,3,3,3) open(90,file=fn) read(90,*) read(90,*) do 1 k=1,3 read(90,*)((a(im,k,i,j),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(im,k,j,i)=a(im,k,i,j) c first index electric close(90) return end c ============================================================ SUBROUTINE READTTT(nat,n,im,alphad,gpd,Ad,ft,fn,n1,n2) IMPLICIT none INTEGER*4 nat,n,im,ix,jx,kx,n1(*),natc,i,j,k,l,ll,ld,n2(*), 1IIND real*8 alphad(3*nat,3,3),Gpd(3*nat,3,3),Ad(3*nat,3,3,3) character*80 ft,fn(n) if(fn(im)(1:4).eq.'ball')then do 11 I=1,3 IIND=3*(n1(im)-1)+I do 11 ix=1,3 do 11 jx=1,3 alphad(IIND,ix,jx)=0.0d0 gpd(IIND,ix,jx)=0.0d0 do 11 kx=1,3 11 ad(IIND,ix,jx,kx)=0.0d0 write(6,*)im,' - derivatives zeroed out' else open(15,file=ft) READ(15,*) READ(15,*)NATC if(NATC.NE.n2(im)-n1(im)+1) 1 call report('Wrong number of atoms in .TTT file') READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 LL=1,NATC L=n1(im)+LL-1 DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)LD,LD,(alphad(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 LL=1,NATC L=n1(im)+LL-1 DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)LD,LD,(GPD(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 LL=1,NATC L=n1(im)+LL-1 DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)LD,LD,(Ad(IIND,I,J,K),K=1,3) c first index electric CLOSE(15) write(6,*)im,' - derivatives read' endif RETURN END c ============================================================ SUBROUTINE READRX(r,n,fn,im,Nattot,n1,n2) IMPLICIT none integer*4 Nat,n,i,im,j,Nattot,n1(*),n2(*) real*8 r(n,3),q,x,y,z,bohr CHARACTER*(*) fn(n) bohr=0.529177d0 OPEN(22,FILE=fn(im)) READ(22,*) read(22,*)Nat n1(im)=Nattot+1 Nattot=Nattot+Nat n2(im)=Nattot q=0.0d0 r(im,1)=0.0d0 r(im,2)=0.0d0 r(im,3)=0.0d0 do 1 i=1,Nat READ(22,*)j,x,y,z x=x/bohr y=y/bohr z=z/bohr q=q+dble(j) r(im,1)=r(im,1)+x*dble(j) r(im,2)=r(im,2)+y*dble(j) 1 r(im,3)=r(im,3)+z*dble(j) r(im,1)=r(im,1)/q r(im,2)=r(im,2)/q r(im,3)=r(im,3)/q close(22) RETURN END c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ subroutine mktttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3,3),r2,r4,r5,r9,delta integer*4 ix,jx,kx,lx r2=rt(1)**2+rt(2)**2+rt(3)**2 r4=r2*r2 r5=r4*dsqrt(r2) r9=r5*r4 do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 do 613 lx=1,3 613 ttt(ix,jx,kx,lx)=3.0d0*(35.0d0*rt(ix)*rt(jx)*rt(kx)*rt(lx) 1-5.0d0*r2*(delta(ix,kx)*rt(jx)*rt(lx) 1 +delta(ix,lx)*rt(jx)*rt(kx) 1 +delta(ix,jx)*rt(kx)*rt(lx) 1 +delta(kx,jx)*rt(lx)*rt(ix) 1 +delta(kx,lx)*rt(ix)*rt(jx) 1 +delta(jx,lx)*rt(ix)*rt(kx)) 1+r4*(delta(ix,lx)*delta(kx,jx) 1 +delta(kx,lx)*delta(ix,jx) 1 +delta(jx,lx)*delta(ix,kx)))/r9 return end c ============================================================ function delta(i,j) implicit none real*8 delta integer*4 i,j if(i.eq.j)then delta=1.0d0 else delta=0.0d0 endif return end c ============================================================ subroutine mkttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3),rt2,rt5 integer*4 ix,jx,kx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 ttt(ix,jx,kx)=-15.0d0*rt(ix)*rt(jx)*rt(kx)/rt5/rt2 if(ix.eq.jx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(kx)/rt5 if(ix.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(jx)/rt5 613 if(jx.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(ix)/rt5 return end c ============================================================ subroutine mktt(rt,tt) implicit none real*8 rt(3),tt(3,3),rt2,rt5 integer*4 ix,jx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 611 ix=1,3 do 611 jx=1,3 611 tt(ix,jx)=3.0d0*rt(ix)*rt(jx)/rt5 do 610 ix=1,3 610 tt(ix,ix)=tt(ix,ix)-rt2/rt5 return end c ============================================================ subroutine mulmg(a,b,c,n) implicit none integer*4 n,i,j,k real*8 a(n,n),b(n,n),c(n,n) do 1 i=1,n do 1 j=1,n c(i,j)=0.0d0 do 1 k=1,n 1 c(i,j)=c(i,j)+a(i,k)*b(k,j) return end c ============================================================ SUBROUTINE INV(A,AI,N,TOL) IMPLICIT none integer*4 N,oo,ii,jj,iw,io,kk,i2,j2 REAL*8 TOL, A(N,N),AI(N,N),w real*8, allocatable ::E(:,:) allocate(E(N,2*N)) DO 1 ii=1,N DO 1 jj=1,N e(ii,jj)=a(ii,jj) E(II,JJ+N)=0.0D0 1 if (ii.EQ.jj)e(ii,jj+N)=1.0D0 DO 2 ii=1,N-1 iw=ii if (DABS(e(iw,iw)).LT.TOL) then DO 3 io=iw+1,N oo=IO if (io.GT.N)oo=io-N 3 if (DABS(e(oo,iw)).GE.TOL) goto 11 call report('Inverse cannot be done') 11 CONTINUE DO 4 kk=1, 2*N w=e(iw,kk) e(iw,kk)=e(oo,kk) 4 e(oo,kk)=w ENDIF DO 2 jj=ii+1,N DO 6 kk=ii+1, 2*N 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii) 2 e(jj,ii)=0.0D0 DO 7 ii=1, N-1 i2=N-ii+1 DO 7 jj=1,i2-1 j2=i2-jj DO 9 kk=1, N 9 e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2) 7 e(j2,i2)=0.0D0 DO 10 ii=1,N DO 10 jj=1,N 10 AI(ii,jj)=e(ii,jj+N)/e(ii,ii) deallocate(E) RETURN END SUBROUTINE WRITETTT(NAT,alpha,A,G,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION alpha(3*NAT,3,3),G(3*NAT,3,3),A(3*NAT,3,3,3) character*(*)s OPEN(2,FILE=s) WRITE(2,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 WRITE(2,2001)L,IX,(alpha(IIND,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F18.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx jy jz') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 WRITE(2,2001)L,IX,(G(IIND,I,J),J=1,3) WRITE(2,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(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3) CLOSE(2) WRITE(6,*)s//' written ' RETURN END FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K 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 READS(S,E,N,nat) implicit none integer*4 N,nint,nat,kd,i,j,ix real*8 S(N,N),E(N),CM,SFAC OPEN(4,FILE='F.INP') read(4,*)nint,nat,nat do 1 i=1,NAT 1 read(4,*) read(4,*) do 4 i=1,3*NAT E(i)=0.0d0 do 4 j=1,3*NAT 4 S(i,j)=0.0d0 DO 2 I=1,NAT DO 2 J=1,NINT 2 read(4,*)kd,kd,(s(3*(i-1)+ix,j),ix=1,3) read(4,*) READ(4,4000)(E(I),I=1,nint) 4000 FORMAT(6F11.6) close(4) CM=219470.0d0 DO 3 I=1,nint 3 E(I)=E(I)/CM SFAC=0.0234280d0 C 1 a.u.= 2.1947E5 cm^-1 DO 5 I=1,N DO 5 J=1,N 5 S(I,J)=S(I,J)*SFAC RETURN end SUBROUTINE TRTEN(S,N,alpha,A,G) IMPLICIT none integer*4 N,i,j,jq,IX,K,IQ real*8 S(N,N),alpha(N,3,3),A(N,3,3,3),G(N,3,3),alphaS,GS,AS real*8,allocatable::V(:) allocate(V(N)) 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) 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) 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) 12 A(IQ,I,J,K)=AS deallocate(V) RETURN END SUBROUTINE DORAMAN(N,alpha,G,A,EXCA,E,Sit,ic) c ic = 0 ... calculate and write ROA.TAB c 1 ... fill Sit c 2 ... only write Sit to ROA.R.TAB c >2 ... calculate and write as (ic-2).tab implicit none integer*4 N,IQ,I,J,K,L,ic,istart real*8 alpha(N,3,3),G(N,3,3),A(N,3,3,3),E(*),EXCA,CM,ECM,SAL0, 1SAL1,SAG1,SA1,SPAL,S180,D,D180,S90X,D90X,DX,S90Z,SAG0, 2D90Z,DZ,YDY,YDX,P1,AMU,BOHR,SpA3,beta2,gpisvejc ,betaa,betag,roa3, 3roa2,ram2,CID2,roa1,ram1,CID1,Sit(2,N),EPS character*10 s10 logical lw CM=219470.0d0 if(ic.eq.2)then OPEN(4,FILE='ROA.R.TAB') WRITE(4,3304) write(4,3002) DO 2 IQ=1,N ECM=E(IQ)*CM IF(ABS(ECM).GT.20.0d0)THEN YDX=0.0d0 YDY=0.0d0 CID1=0.0d0 CID2=0.0d0 DX=0.0d0 P1=0.0d0 roa3=0.0d0 ram1=sit(1,iq) roa1=sit(2,iq) WRITE(4,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3 endif 2 continue WRITE(4,3002) CLOSE(4) endif if(ic.eq.0)OPEN(4,FILE='ROA.TAB') if(ic.gt.2)then write(s10,100)ic-2 100 format(i10) do 101 istart=1,len(s10) 101 if(s10(istart:istart).ne.' ')goto 102 102 open(4,file=s10(istart:len(s10))//'.tab') endif lw=ic.eq.0.or.ic.ge.2 if(lw)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(lw)write(4,3002) 3002 FORMAT(80(1H-)) c DO 1 IQ=1,N ECM=E(IQ)*CM IF(ABS(ECM).GT.20.0d0)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 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 alphag=SAG0/9.0d0*gpisvejc 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(lw)then WRITE(4,3000)IQ,ECM,YDX,YDY,ram1,CID2,DX,CID1,P1,roa1,roa3 3000 FORMAT(I5,f9.2,6g12.3,f6.3,2g12.4) else sit(1,iq)=ram1 sit(2,iq)=roa1 endif ENDIF 1 CONTINUE if(lw)then WRITE(4,3002) CLOSE(4) endif RETURN END subroutine doX(n,r,X) implicit none integer*4 i,n,istart,j,jstart,ix,jx,kx,lx,ii,jj real*8 rij(3),r(n,3),tt(3,3),ttt(3,3,3),tttt(3,3,3,3), 1X(n*24,n*24) do 6 i=1,n istart=24*(i-1)+1 do 6 j=1,n jstart=24*(j-1)+1 if(i.ne.j)then do 61 ix=1,3 61 rij(ix)=r(i,ix)-r(j,ix) c T 0 0 0 -U 0 T=grad(1/rij), U=grad(T), G=grad(U), W=T/c c 0 T 0 0 0 -U c 0 0 W 0 0 0 c 0 0 0 W 0 0 c U 0 0 0 -G 0 c 0 U 0 0 0 -G call mktt(rij,tt) do 62 ix=1,3 do 62 jx=1,3 X(istart+ix-1, jstart+jx-1 )=tt(ix,jx) X(istart+ix-1+3,jstart+jx-1+3)=tt(ix,jx) X(istart+ix-1+6,jstart+jx-1+6)=tt(ix,jx)/137.0d0 62 X(istart+ix-1+9,jstart+jx-1+9)=tt(ix,jx)/137.0d0 call mkttt(rij,ttt) do 65 ix=1,3 ii=0 do 65 jx=1,3 do 65 kx=jx,3 ii=ii+1 X(istart+ix-1,jstart+11+ii)=-ttt(ix,jx,kx) X(istart+11+ii,jstart+ix-1)= ttt(ix,jx,kx) X(istart+ix+2,jstart+17+ii)=-ttt(ix,jx,kx) 65 X(istart+17+ii,jstart+ix+2)= ttt(ix,jx,kx) call mktttt(rij,tttt) ii=0 do 64 ix=1,3 do 64 jx=ix,3 ii=ii+1 jj=0 do 64 kx=1,3 do 64 lx=kx,3 jj=jj+1 X(istart+ii+11,jstart+11+jj)=-tttt(ix,jx,kx,lx) 64 X(istart+ii+17,jstart+17+jj)=-tttt(ix,jx,kx,lx) endif 6 continue return end c ============================================================ subroutine MDP(x,p,n,tem,temi,temk,pt,alphar,gpr,ar,nat,n1,n2, 1pd,teml,temp,temq,r,lusep1,lusep2) implicit none integer*4 n,i,j,nat,k,l,n1(*),n2(*),ix,ii,jj,kk,kkk,jjj, 1im,ixx,m,kl,ia,ik real*8 x(24*n,24*n),p(24*n,24*n),tem(24*n,24*n),temi(24*n,24*n), 1temk(24*n,24*n),pt(24*n,24*n),alphar(3*nat,3,3),gpr(3*nat,3,3), 2ar(3*nat,3,3,3),pd(3*nat,24,24),teml(24*n,24*n),temp(24*n,24*n), 3temq(24*n,24*n),r(n,3),sum,eps,rr(3),delta real*8, allocatable::ptd(:,:) logical lusep1,lusep2 c c make (E-X.P) call mulmg(x,p,tem,n*24) do 881 i=1,n*24 do 881 j=1,n*24 881 tem(i,j)=-tem(i,j) do 991 i=1,n*24 991 tem(i,i)=1.0d0+tem(i,i) c make (E-X.P)^(-1) = temi call INV(tem,temi,24*n,1.0d-25) c make total polarizability Pt = P(E-X.P)^(-1) call mulmg(p,temi,pt,n*24) c make P(E-X.P)^(-1)X = pt . X =temk call mulmg(pt,x,temk,24*n) do 900 i=1,3*nat do 900 j=1,3 do 900 k=1,3 alphar(i,j,k)=0.0d0 gpr(i,j,k)=0.0d0 do 900 l=1,3 900 ar(i,j,k,l)=0.0d0 c c calculate Pt derivatives, local origins c = P'(E-X.P)^(-1)+P(E-X.P)^(-1).X.P'(E-X.P)^(-1) allocate(ptd(24*n,24*n)) c c loop over the atomic coordinates c LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL do 8 i=1,n do 8 ia=n1(i),n2(i) do 8 ix=1,3 ii=3*(ia-1)+ix c c P': do 81 j=1,n do 81 jj=1,24 jjj=jj+24*(j-1) do 81 k=1,n do 81 kk=1,24 kkk=kk+24*(k-1) if(j.eq.i.and.k.eq.i)then temp(jjj,kkk)=pd(ii,jj,kk) else temp(jjj,kkk)=0.0d0 endif 81 continue c c make P'(E-X.P)^(-1) =teml call mulmg(temp,temi,teml,24*n) c make P(E-XP)^(-1)XP'(E-X.P)^(-1)= temk.teml=temq call mulmg(temk,teml,temq,24*n) if(.not.lusep1)call tz(teml,24*n) if(.not.lusep2)call tz(temq,24*n) do 83 j=1,24*n do 83 k=1,24*n 83 ptd(j,k)=teml(j,k)+temq(j,k) do 8 im=1,n do 84 ixx=1,3 84 rr(ixx)=r(im,ixx) do 8 j=1,3 do 8 ik=1,n kl=0 do 8 k=1,3 c from electric moment 1..3: c from electric field 1..3: alphar(ii,j,k)=alphar(ii,j,k)+ptd(24*(im-1)+j,24*(ik-1)+k) c from magnetic moment 7-9: c from electric field derivative 4..6: sum=0.0d0 do 85 l=1,3 do 85 m=1,3 85 sum=sum+eps(k,l,m)*rr(l)*ptd(24*(im-1)+3+m,24*(ik-1)+j+3)/2.0d0 gpr(ii,j,k)=gpr(ii,j,k)-ptd(24*(im-1)+6+k,24*(ik-1)+3+j)-sum do 8 l=k,3 kl=kl+1 c from quadrupole 13-18 c from electric field 1..3: c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 ar(ii,j,k,l)=ar(ii,j,k,l)+ptd(24*(im-1)+kl+12,24*(ik-1)+j)*3.0d0 1+1.5d0 *(rr(l)*ptd(24*(im-1)+k,24*(ik-1)+j) 1 +rr(k)*ptd(24*(im-1)+l,24*(ik-1)+j)) 1-delta(k,l)*(rr(1)*ptd(24*(im-1)+1,24*(ik-1)+j) 1+ rr(2)*ptd(24*(im-1)+2,24*(ik-1)+j) 1+ rr(3)*ptd(24*(im-1)+3,24*(ik-1)+j)) 8 ar(ii,j,l,k)=ar(ii,j,k,l) deallocate(ptd) return end c ============================================================ subroutine rewrit(alpharm,gprm,arm,natm,NAT,alphar,gpr,ar,i,n1,n2) implicit none integer*4 natm,n1(*),n2(*),ia,ix,jx,kx,i,NAT,iax,ii,jj real*8 alpharm(3*natm,3,3),gprm(3*natm,3,3),arm(3*natm,3,3,3) real*8 alphar(3*nat,3,3),gpr(3*nat,3,3),ar(3*nat,3,3,3) do 11 ia=n1(i),n2(i) do 11 iax=1,3 ii=3*(ia-1)+iax jj=3*(ia-n1(i))+iax do 11 ix=1,3 do 11 jx=1,3 alpharm(jj,ix,jx)=alphar(ii,ix,jx) gprm(jj,ix,jx)=gpr(ii,ix,jx) do 11 kx=1,3 11 arm(jj,ix,jx,kx)=ar(ii,ix,jx,kx) return end c ============================================================ SUBROUTINE TTT1(alpha,A,G,ICO,r,n,alphad,gd,ad,nat,n1,n2,lusea, 1luseg,ltr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C c supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,n,iu,IA,IB,IG,ID,IC,ix,n1(*),n2(*),iat,iax,ii,nat real*8 alpha(n,3,3),G(n,3,3),A(n,3,3,3),c(3),SIGN,SUM,EPS real*8 alphad(3*nat,3,3),Gd(3*nat,3,3),Ad(3*nat,3,3,3) real*8 r(n,3) logical luseg,lusea,ltr(*) do 99 iu=1,n if(ltr(iu))then do 52 ix=1,3 52 c(ix)=r(iu,ix) C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*c(IG)*alpha(iu,IA,ID) G(iu,IA,IB)=G(iu,IA,IB)+SUM IF(ICO.EQ.3.or..not.luseg)G(iu,IA,IB)=0.0d0 do 21 iat=n1(iu),n2(iu) do 21 iax=1,3 ii=3*(iat-1)+iax sum=0.0d0 DO 11 IG=1,3 DO 11 ID=1,3 11 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*c(IG)*alphad(ii,IA,ID) Gd(ii,IA,IB)=Gd(ii,IA,IB)+SUM c G last index (IB) magnetic 21 IF(ICO.EQ.3.or..not.luseg)Gd(ii,IA,IB)=0.0d0 2 CONTINUE C c A .. first index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1c(1)*alpha(iu,IA,1)+c(2)*alpha(iu,IA,2)+c(3)*alpha(iu,IA,3) A(iu,IA,IC,IB)=A(iu,IA,IC,IB)+ 1SIGN*(SUM-1.5d0*(c(IB)*alpha(iu,IA,IC)+c(IC)*alpha(iu,IA,IB))) IF(ICO.EQ.3.or..not.lusea)A(iu,IA,IC,IB)=0.0d0 do 31 iat=n1(iu),n2(iu) do 31 iax=1,3 ii=3*(iat-1)+iax SUM=0.0d0 IF(IB.EQ.IC)SUM= 1c(1)*alphad(ii,IA,1)+c(2)*alphad(ii,IA,2)+c(3)*alphad(ii,IA,3) Ad(ii,IA,IC,IB)=Ad(ii,IA,IC,IB)+ 1SIGN*(SUM-1.5d0*(c(IB)*alphad(ii,IA,IC)+c(IC)*alphad(ii,IA,IB))) 31 IF(ICO.EQ.3.or..not.lusea)Ad(ii,IA,IC,IB)=0.0d0 3 CONTINUE write(6,600)iu 600 format(' molecule ',i4,' transformed') else write(6,601)iu 601 format(' molecule ',i4,' not transformed') endif 99 continue RETURN END c ============================================================ subroutine tz(t,n) implicit none integer*4 n,i,j real*8 t(n,n) do 1 i=1,n do 1 j=1,n 1 t(i,j)=0.0d0 return end