PROGRAM IRROA2 IMPLICIT none integer*4 n,i,ix,ii,nm,NW,ntt,ic,j,jx,is,l,ls, 1m,iwr,nx,ny,ncub,icen,jfreq,k, 1ndim,ntmax,nxc,nyc,nzc,ncubc,icubc real*8 tol,alpha0(9),gp0(9),a0(27),wnm,wau,wcm, 1WMIN,WMAX,DW,GAMMA,W0,TEMP,tempcm,tolw,wa,x0(3),v(3), 1u(3),dxc,dyc,dzc,x0c(3),cm,gw character*80,allocatable:: fn(:),ftw(:),ftwas(:) logical,allocatable:: ltr(:) logical lave,luseg,lusea,debug,lcont,lff,lgauss, 1ldvm,llog,lcube,ltab,lonly,lasym integer*4 ,allocatable::ntw(:),iiz(:) real*8,allocatable:: alpha(:,:,:),gp(:,:,:),a(:,:,:,:), 1alphad(:,:,:,:),gpd(:,:,:,:),ad(:,:,:,:,:),r(:,:),X(:,:), 1alphar(:,:,:),gpr(:,:,:),ar(:,:,:,:),Ptr(:,:,:), 1wl(:),EXP0(:,:),P0(:,:),XP0(:,:), 1tw(:,:),Pt(:,:),EPXM(:,:),PR(:,:), 1alphaas(:,:,:),alphadas(:,:,:,:), 1PRX(:,:),PRXP(:,:) n=1 nm=1 ntmax=1 c initial allocation: allocate(r(n,3),fn(n),ftw(n),ntw(n),ltr(n),ftwas(n), 1alpha(n,3,3),gp(n,3,3),a(n,3,3,3),alphad(n,nm,3,3),gpd(n,nm,3,3), 1ad(n,nm,3,3,3),tw(n,ntmax),iiz(n), 1alphaas(n,3,3),alphadas(n,nm,3,3)) do 3301 ic=0,2 call readopt(ic,n,nm,m,ntw,ntmax,NW,nx,ny,ncub,iwr, 1nxc,nyc,nzc,ncubc,icubc,ndim, llog,ltab,lff,ldvm,lgauss, 1luseg,lusea,lcont,icen,jfreq,lcube,ltr,lave,debug, 1lonly,r,alpha,gp,a,alphad,gpd,ad,tw,WMAX,WMIN,W0,GAMMA,WA, 1temp,tol,tolw,x0,u,v,gw,x0c,dxc,dyc,dzc,DW,tempcm,cm, 1fn,ftw,ftwas,iiz,lasym,alphaas,alphadas) if(ic.eq.0)then write(6,*)n,' groups' deallocate(r,fn,ftw,ntw,ltr,alpha,gp,a,alphad,gpd,ad,tw,iiz, 1 ftwas,alphadas,alphaas) c temporary allocation: allocate(r(n,3),fn(n),ftw(n),ntw(n),ltr(n),ftwas(n), 1 alpha(n,3,3),gp(n,3,3),a(n,3,3,3),alphad(n,nm,3,3),gpd(n,nm,3,3), 1 ad(n,nm,3,3,3),tw(n,ntmax),iiz(n),alphadas(n,nm,3,3), 1 alphaas(n,3,3)) endif if(ic.eq.1)then deallocate(r,fn,ftw,ntw,ltr,alpha,gp,a,alphad,gpd,ad,tw,iiz, 1 ftwas,alphadas,alphaas) c final allocation: allocate(r(n,3),fn(n),ftw(n),ntw(n),ltr(n),ftwas(n), 1 alpha(n,3,3),gp(n,3,3),a(n,3,3,3),alphad(n,nm,3,3),gpd(n,nm,3,3), 1 ad(n,nm,3,3,3),tw(n,ntmax),iiz(n),alphadas(n,nm,3,3), 1 alphaas(n,3,3)) endif 3301 continue c delete derivatives of desired particles: call iizd(n,nm,iiz,alphad,gpd,ad) call wropt(n,ndim,ntmax,NW,iwr,WMIN,WMAX,DW,GAMMA,W0,tolw,TEMP, 1wa,tol,gw,lave,lusea,luseg,lff,lcont,ldvm,lonly,lasym) c transform G, A to molecular centers call TTT1(alpha,a,gp,1,r,n,alphad,gpd,ad,lusea,luseg, 1nm,ntw,ltr) write(6,*)' Tensors transformed to molecular centers' c eliminate G, A if requested: if(lonly)then gpd=0.0d0 ad=0.0d0 gp=0.0d0 a=0.0d0 write(6,*)' G, A erased' endif allocate(X(ndim,ndim),Ptr(ntmax,ndim,ndim), 1P0(ndim,ndim),EXP0(ndim,ndim), 1XP0(ndim,ndim),PRX(ndim,ndim),PR(ndim,ndim)) Ptr=0.0d0 write(6,*)' X,P allocated' c Calculate general distance tensor: call doX(n,r,X,iwr) c Construct general polarizability tensor c (transition polarizabilities for each group and frequency): call ctp(n,nm,ntmax,ndim,Ptr,alphad,gpd,ad,lonly,ntw) c P0:ordinary (non-transition) polarizability: call ct0(n,ndim,P0,alpha,gp,a,lonly) c EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0 write(6,*)'X:' do 246 k=1,n is=24*(k-1)+1 write(6,*)k do 246 l=k+1,n ls=24*(l-1)+1 246 write(6,600)l,((X(is+ix,ls+jx),ix=0,2),jx=0,2) write(6,*)'P:' do 248 k=1,n is=24*(k-1)+1 248 write(6,600)k,((P0(is+ix,is+jx),ix=0,2),jx=0,2) c (E - X.P0): call mm(ndim,XP0,X,P0) write(6,*)'XP:' do 247 k=1,n is=24*(k-1)+1 do 247 l=k+1,n ls=24*(l-1)+1 247 write(6,600)k,((XP0(is+ix,ls+jx),ix=0,2),jx=0,2) EXP0=-XP0 do 75 i=1,ndim 75 EXP0(i,i)=1.0d0+EXP0(i,i) c make (E-XP0)-1: allocate(EPXM(ndim,ndim)) call INV(EXP0,EPXM,ndim,1.0d-25) write(6,*)'(E-XP0)-1:' do 245 k=1,n is=24*(k-1)+1 245 write(6,600)k,((EPXM(is+ix,is+jx),ix=0,2),jx=0,2) c make redressed polarizability: P~ = P0.(E-XP0)-1: call mm(ndim,PR,P0,EPXM) c make 1 + P~X: call mm(ndim,PRX,PR,X) do 713 ix=1,ndim 713 PRX(ix,ix)=PRX(ix,ix)+1.0d0 c EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0 c c save dressed polarizabilities for excitation frequency: call extr0(PR,alpha0,gp0,a0,ndim,r,n) wnm=wa/10.0d0 wcm=1.0d7/wnm wau=wcm/cm call WRITEPOL(alpha0 ,'POL0.TTT',wau) call WRITEPOL(gp0 ,'GP0.TTT' ,wau) call WRITEPOL(a0 ,'A0.TTT' ,wau) c each group and frequency scatters independently: ntt=0 do 1 i=1,n do 1 j=1,ntw(i) 1 ntt=ntt+1 write(6,*)ntt,' Raman centers' allocate(wl(ntt)) allocate(alphar(ntt,3,3),gpr(ntt,3,3),ar(ntt,3,3,3), 1Pt(ndim,ndim),PRXP(ndim,ndim)) alphar=0.0d0 gpr=0.0d0 ar=0.0d0 c for each frequency/center make polarizabilities: ii=0 do 222 i=1,n do 222 j=1,ntw(i) ii=ii+1 wl(ii)=tw(i,j) write(6,*)'greoup',i,'freq',j c fill transition polarizability P'=Pt for frequency ii: call ftp(i,j,ndim,ntmax,Pt,Ptr) write(6,*)'original' do 242 k=1,n is=24*(k-1)+1 242 write(6,600)k,((Pt(is+ix,is+jx),ix=0,2),jx=0,2) 600 format(i6,/,3(3f10.6,/)) c redress P' polarizability, put it back to Pt: c (1 + P~X)P': call mm(ndim,PRXP,PRX,Pt) write(6,*)'(1 + P~X)Pp:' do 243 k=1,n is=24*(k-1)+1 243 write(6,600)k,((PRXP(is+ix,is+jx),ix=0,2),jx=0,2) c (1 + P~X) P' (1-XP)^-1: call mm(ndim,Pt,PRXP,EPXM) write(6,*)'(1 + P~X) Pp (1-XP)^-1:' do 244 k=1,n is=24*(k-1)+1 244 write(6,600)k,((Pt(is+ix,is+jx),ix=0,2),jx=0,2) c save polarizabilities for this frequency: 222 call extr(Pt,ii,alphar,gpr,ar,ntt,ndim,r,n) c do i=1,ntt c write(6,655)i,((alphar(i,ix,jx),ix=1,3),jx=1,3) c55 format(i6,/,3(3f12.6,/)) c enddo c for each freq, make intensities, add to spectrum, save: call DORAMAS(NW,GAMMA,ntt,alphar,Gpr,Ar,wa,WMIN,DW, 1tempcm,wl,lgauss,ltab) call WRITETTTQ(ntt,1,ntt,alphar,Ar,Gpr,alphar,wl) if(icen.eq.0)call PLA0(n,ndim,x0,u,v,nx,ny,r,llog,EPXM) if(icen.gt.0.and.jfreq.gt.0)then call ftp(icen,jfreq,ndim,ntmax,Pt,Ptr) call mm(ndim,PRXP,PRX,Pt) call mm(ndim,Pt,PRXP,EPXM) call PLAN(n,ndim,x,Pt,PR,x0,u,v,nx,ny,icen,jfreq, 1 r,llog,EPXM,tw,ntmax) endif if(lcube)call CUBE(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0c,nxc,nyc,nzc,ncubc,icubc,dxc,dyc,dzc,r,llog) 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 READTTTQ(nm,n,alphad,Gpd,ad,ftw,tw,ip,ntw,ntmax) IMPLICIT none INTEGER*4 nm,nq,LL,I,J,K,n,ip,ntmax,ntw(*) real*8 alphad(n,nm,3,3),Gpd(n,nm,3,3),Ad(n,nm,3,3,3),tw(n,ntmax) character*80 ftw(*) open(15,file=ftw(ip)) READ(15,*) READ(15,*)nq tw(ip,nq+1)=0.0d0 if(nq.ne.ntw(ip))call report('nq<>ntw') if(nq.eq.0)goto 99 READ(15,*) READ(15,*) DO 1 LL=1,nq READ(15,*)tw(ip,LL),tw(ip,LL) DO 1 I=1,3 1 READ(15,*)alphad(ip,LL,I,1),(alphad(ip,LL,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 LL=1,nq READ(15,*) DO 2 I=1,3 2 READ(15,*)GPD(ip,LL,I,1),(GPD(ip,LL,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 LL=1,nq READ(15,*) DO 3 I=1,3 DO 3 J=1,3 3 READ(15,*)(Ad(ip,LL,I,J,K),K=1,2),(Ad(ip,LL,I,J,K),K=1,3) c first index electric 99 CLOSE(15) RETURN END c ============================================================ SUBROUTINE RDAS(nm,n,alphad,ftw,tw,ip,ntw,ntmax) IMPLICIT none INTEGER*4 nm,nq,LL,I,J,n,ip,ntmax,ntw(*) real*8 alphad(n,nm,3,3),tw(n,ntmax) c Gpd(n,nm,3,3),Ad(n,nm,3,3,3), character*80 ftw(*) open(15,file=ftw(ip)) READ(15,*) READ(15,*)nq tw(ip,nq+1)=0.0d0 if(nq.ne.ntw(ip))call report('nq<>ntw') if(nq.eq.0)goto 99 READ(15,*) READ(15,*) DO 1 LL=1,nq READ(15,*)tw(ip,LL),tw(ip,LL) DO 1 I=1,3 1 READ(15,*)alphad(ip,LL,I,1),(alphad(ip,LL,I,J),J=1,3) 99 CLOSE(15) RETURN END c ============================================================ SUBROUTINE READRX(r,n,fn,im,ic) IMPLICIT none integer*4 Nat,n,i,im,j,ic 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 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 if(ic.eq.1)write(6,600)im,Nat,(r(im,i),i=1,3) 600 format(' center ',i4,', ',i4,' atoms',3F9.4) 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 if(r2.gt.1.0d-5)then 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 else ttt=0.0d0 endif 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) if(rt2.gt.1.0d-5)then 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 else ttt=0.0d0 endif 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) if(rt2.gt.1.0d-5)then 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 else tt=0.0d0 endif return end c ============================================================ 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 c ============================================================ SUBROUTINE DORAMAS(NW,GAMMA,ntt,alpha,G,A,EXCA,WMIN,DW,tempcm, 1wl,lgauss,ltab) implicit none integer*4 NW,IQ,I,ntt real*8 alpha(ntt,3,3),G(ntt,3,3),A(ntt,3,3,3),EXCA,ECM,sf,dd, 3roa1,ram1,WMIN,DW,tempcm,const,wl(*),GAMMA,t,doc,ro,YDY,YDX real*8,allocatable::s(:,:) logical lgauss,ltab allocate(s(2,NW)) s=0.0d0 if(lgauss)then dd=GAMMA/2.0d0/dsqrt(log(2.0d0)) else dd=GAMMA/2.0d0 endif if(ltab)then open(20,file='ROA.TAB') write(20,200) 200 format('IRROA model',/,' freq ram roa ') write(20,3002) 3002 FORMAT(80(1H-)) endif DO 1 IQ=1,ntt ECM=wl(IQ) call rr(ram1,roa1,IQ,alpha,G,A,ntt,EXCA,doc,ro,YDY,YDX) if(ltab)write(20,201) 1IQ,ECM,YDX,YDY,ram1,ram1,ram1,ram1,ro,roa1,doc 201 FORMAT(I5,f11.4,12g12.4) if(ECM.gt.1.0d0)then const=1.0d0/ECM/(1.0d0-exp(-ECM/tempcm)) do 4 i=1,NW t=sf(dd,lgauss,WMIN+DW*dble(i-1),ECM)*const s(1,i)=s(1,i)+t*ram1 4 s(2,i)=s(2,i)+t*roa1 endif 1 continue if(ltab)then write(20,3002) close(20) endif OPEN(41,FILE='rams.prn') OPEN(42,FILE='roas.prn') do 2 i=1,NW WRITE(41,3000)WMIN+DW*dble(i-1),s(1,i) 2 WRITE(42,3000)WMIN+DW*dble(i-1),s(2,i) 3000 FORMAT(f9.2,g12.4) close(41) close(42) RETURN END c ============================================================ subroutine doX1(n,r1,r,X) c field in one place implicit none integer*4 n,j,jstart,ix,jx,kx,lx,ii,jj real*8 rij(3),r1(3),tt(3,3),ttt(3,3,3),tttt(3,3,3,3), 1X(24,n*24),r(n,3),d x=0.0d0 do 6 j=1,n jstart=24*(j-1)+1 do 61 ix=1,3 61 rij(ix)=r1(ix)-r(j,ix) d=rij(1)**2+rij(2)**2+rij(3)**2 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 if(d.gt.1.0d-3)then call mktt(rij,tt) do 62 ix=1,3 do 62 jx=1,3 X(ix, jstart+jx-1 )=tt(ix,jx) X(ix+3,jstart+jx-1+3)=tt(ix,jx) X(ix+6,jstart+jx-1+6)=tt(ix,jx)/137.0d0 62 X(ix+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(ix ,jstart+11+ii)=-ttt(ix,jx,kx) X(12 +ii,jstart+ix- 1)= ttt(ix,jx,kx) X(ix + 3,jstart+17+ii)=-ttt(ix,jx,kx) 65 X(18 +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(ii+12,jstart+11+jj)=-tttt(ix,jx,kx,lx) 64 X(ii+18,jstart+17+jj)=-tttt(ix,jx,kx,lx) endif 6 continue return end c ============================================================ subroutine doX(n,r,X,iwr) implicit none integer*4 i,n,istart,j,jstart,ix,jx,kx,lx,ii,jj,iwr real*8 rij(3),r(n,3),tt(3,3),ttt(3,3,3),tttt(3,3,3,3), 1X(n*24,n*24) x=0.0d0 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 write(6,*)' Distance tensor done' if(iwr.gt.1)call wmtr(X,n*24,n*24,'X-tensor',1,6) return end c ============================================================ subroutine PLA0(n,ndim,x0,u,v,nx,ny,r,llog,PR) c c -1 c F(x) = (1 + X(x)*P*(1-X.P) ) . E0 = X(x) * P~ .F0 c implicit none integer*4 ndim,n,nx,ny real*8 r(n,3),u(3),x0(3),v(3),PR(ndim,ndim) logical llog write(6,600)nx,ny,x0,u,v 600 format(' Constructing plane ',i5,' x ',i5,/, 1 ' Origin: ',3f12.3,' A',/, 2 ' vector1 ',3f12.3,' A',/, 2 ' vector2 ',3f12.3,' A',/, 4 ' Excitation frequency',/) call wrp0(ndim,PR,x0,u,v,n,nx,ny,r,llog,0,0) return end c ============================================================ subroutine PLAN(n,ndim,x,Pp,PR,x0,u,v,nx,ny,icen,ifr,r, 1llog,EPXM,tw,ntmax) c c -1 c F' = X(x) * (1+X.P~) . P' . (1-X.P) . E0 c c -1 -1 -1 c = X(x).[P'.(1-XP) + P . (1-XP) .X.P'(1-XP) ]F0 c c -1 -1 c = X(x).[P'.(1-XP) + P~ .X.P'(1-XP) ]F0 c c c = X(x).[ t1 + P~ .X. t1 ]F0 c = X(x).t1[ t2 ]F0 c = X(x).t3 F0 c c implicit none integer*4 ndim,i,n,nx,ny,icen,ifr,ntmax real*8 x(n*24,n*24),Pp(ndim,ndim),r(n,3),u(3),x0(3),v(3), 1EPXM(ndim,ndim),PR(ndim,ndim),tw(n,ntmax) logical llog real*8,allocatable::t1(:,:),t2(:,:),PI(:,:) write(6,600)nx,ny,x0,u,v,icen,ifr,tw(icen,ifr) 600 format(' Constructing plane ',i5,' x ',i5,/, 1 ' Origin: ',3f12.3,' A',/, 2 ' v1 ',3f12.3,' A',/, 2 ' v2 ',3f12.3,' A',/, 4 ' Centrum ',i5,/, 4 ' Frequency ',i5,f12.2,/) allocate(t1(ndim,ndim),t2(ndim,ndim),PI(ndim,ndim)) c t1 = P'(1-XP)^(-1) call mm(ndim,t1,Pp,EPXM) c t2= 1 + P~ X call mm(ndim,t2,PR,X) do 3 i=1,ndim 3 t2(i,i)=t2(i,i) + 1.0d0 c PI= t1 . t2 = P'(1-XP)^(-1) . (1 + P~ X) call mm(ndim,PI,t1,t2) call wrp0(ndim,PI,x0,u,v,n,nx,ny,r,llog,icen,ifr) return end c ============================================================ subroutine CUBE(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0,nx,ny,nz,ncub,icub,dx,dy,dz,r,llog) implicit none integer*4 ndim,i,n,ntmax, 1nx,ny,nz,ncub,icub real*8 x(n*24,n*24),Ptr(ntmax,ndim,ndim),r(n,3), 1EXP0(ndim,ndim),P0(ndim,ndim),x0(3),dx,dy,dz,bohr logical llog real*8,allocatable::EPXM(:,:),t1(:,:),t2(:,:),PI(:,:), 1t3(:,:) write(6,600)nx,ny,nz,dx,dy,dz,x0,ncub,icub 600 format(/,' Constructing cube ',i6 ,' x',i6, ' x',i6/, 1 ' Grid: ',f6.3,' x',f6.3,' x',f6.3,' A',/, 1 ' Origin: ',3f12.3,' A',/, 3 ' Group ',i5,/, 4 ' Frequency number ',i5,/) bohr=0.529177d0 dx=dx/bohr dy=dy/bohr dz=dz/bohr c make (E-XP0)-1 allocate(EPXM(ndim,ndim)) call INV(EXP0,EPXM,ndim,1.0d-25) c allocate(t1(ndim,ndim),t2(ndim,ndim),t3(ndim,ndim)) c t1 = P0 (E-XP0)-1: call mm(ndim,t1,P0,EPXM) if(icub.eq.0)then c excitation frequency: t3=t1 else c vibrational frequency: c t2 = E+ P0(E-XP0)-1 X: call mm(ndim,t2,t1,X) do 3 i=1,ndim 3 t2(i,i)=t2(i,i) + 1.0d0 c take original polarizability: allocate(PI(ndim,ndim)) PI=0.0d0 call getpi(PI,ndim,icub,ncub,Ptr,ntmax) c total polarizability Pt =(E+P0 (E-XP0)-1 X) Pi (E-P0.X)-1 c t1 = PI (E-XP0)-1: call mm(ndim,t1,PI,EPXM) c t3 = (E + P0.(E-XP0)-1 X ) Pi (E-XP0)-1 call mm(ndim,t3,t2,t1) endif c partially contract last index over particles: call pci(n,t3) call wrcube(ndim,t3,x0,n,nx,ny,nz,dx,dy,dz,r,llog) return end c ============================================================ SUBROUTINE TTT1(alpha,A,G,ICO,r,n,alphad,gd,ad,lusea, 1luseg,nm,ntw,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,ntw(*),iw,nm real*8 alpha(n,3,3),G(n,3,3),A(n,3,3,3),SIGN,SUM,EPS, 1alphad(n,nm,3,3),Gd(n,nm,3,3),Ad(n,nm,3,3,3),r(n,3) logical luseg,lusea,ltr(*) c IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C do 99 iu=1,n if(ltr(iu))then 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)*r(iu,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 iw=1,ntw(iu) sum=0.0d0 DO 11 IG=1,3 DO 11 ID=1,3 11 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*r(iu,IG)*alphad(iu,iw,IA,ID) Gd(iu,iw,IA,IB)=Gd(iu,iw,IA,IB)+SUM c G last index (IB) magnetic 21 IF(ICO.EQ.3.or..not.luseg)Gd(iu,iw,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=r(iu,1)*alpha(iu,IA,1)+ 1r(iu,2)*alpha(iu,IA,2)+r(iu,3)*alpha(iu,IA,3) A(iu,IA,IC,IB)=A(iu,IA,IC,IB)+SIGN* 1(SUM-1.5d0*(r(iu,IB)*alpha(iu,IA,IC)+r(iu,IC)*alpha(iu,IA,IB))) IF(ICO.EQ.3.or..not.lusea)A(iu,IA,IC,IB)=0.0d0 do 31 iw=1,ntw(iu) SUM=0.0d0 IF(IB.EQ.IC)SUM= alphad(iu,iw,IA,1)*r(iu,1) 1 + alphad(iu,iw,IA,2)*r(iu,2) 2 + alphad(iu,iw,IA,3)*r(iu,3) Ad(iu,iw,IA,IC,IB)= Ad(iu,iw,IA,IC,IB)+ 1SIGN*(SUM-1.5d0*( alphad(iu,iw,IA,IC)*r(iu,IB) 2 + alphad(iu,iw,IA,IB)*r(iu,IC))) 31 IF(ICO.EQ.3.or..not.lusea)Ad(iu,iw,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 ============================================================ function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d10)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end c ============================================================ subroutine extr(Pt,iw,alphar,gpr,ar,NW,ndim,r,n) implicit none integer*4 i,j,k,l,ndim,NW,iw,n,kl,m,ii,pp,ip real*8 Pt(ndim,ndim),alphar(NW,3,3),gpr(NW,3,3), 2ar(NW,3,3,3),rr(3),r(n,3),eps,delta,sum c c sum over all multipoles on particles i: do 8 i=1,n ii=24*(i-1) rr(1)=r(i,1) rr(2)=r(i,2) rr(3)=r(i,3) do 8 j=1,3 do 8 ip=1,n pp=24*(ip-1) kl=0 do 8 k=1,3 c from electric moment j = 1..3: c from electric field k = 1..3 on all particles ip: alphar(iw,j,k)=alphar(iw,j,k)+Pt(ii+j,pp+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)*Pt(ii+3+m,pp+3+j)/2.0d0 gpr(iw,j,k)=gpr(iw,j,k) -Pt(ii+6+k,pp+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(iw,j,k,l)=ar(iw,j,k,l) 1 +3.0d0*Pt(ii+kl+12,pp+j) 1+1.5d0 *(rr(l)*Pt(ii+k ,pp+j) 1 +rr(k)*Pt(ii+l ,pp+j)) 1-delta(k,l)*(rr(1)*Pt(ii+1 ,pp+j) 1+ rr(2)*Pt(ii+2 ,pp+j) 1+ rr(3)*Pt(ii+3 ,pp+j)) 8 ar(iw,j,l,k)=ar(iw,j,k,l) return end c ============================================================ subroutine extr0(Pt,alphar,gpr,ar,ndim,r,n) implicit none integer*4 i,j,k,l,ndim,n,kl,m,ii,pp,ip real*8 Pt(ndim,ndim),alphar(9),gpr(9), 2ar(27),rr(3),r(n,3),eps,delta,sum c alphar=0.0d0 gpr=0.0d0 ar=0.0d0 c sum over all multipoles on particles i: do 8 i=1,n ii=24*(i-1) rr(1)=r(i,1) rr(2)=r(i,2) rr(3)=r(i,3) do 8 j=1,3 do 8 ip=1,n pp=24*(ip-1) kl=0 do 8 k=1,3 c from electric moment j = 1..3: c from electric field k = 1..3 on all particles ip: alphar(j+3*(k-1))=alphar(j+3*(k-1))+Pt(ii+j,pp+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)*Pt(ii+3+m,pp+3+j)/2.0d0 gpr(j+3*(k-1))=gpr(j+3*(k-1)) -Pt(ii+6+k,pp+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(j+3*(k-1)+9*(l-1))=ar(j+3*(k-1)+9*(l-1)) 1+3.0d0*Pt(ii+kl+12,pp+j) 1+1.5d0 *(rr(l)*Pt(ii+k ,pp+j) 1 +rr(k)*Pt(ii+l ,pp+j)) 1-delta(k,l)*(rr(1)*Pt(ii+1 ,pp+j) 1+ rr(2)*Pt(ii+2 ,pp+j) 1+ rr(3)*Pt(ii+3 ,pp+j)) 8 ar(j+3*(l-1)+9*(k-1))=ar(j+3*(k-1)+9*(l-1)) return end c ============================================================ subroutine rr(ram1,roa1,IQ,alpha,G,A,N,EXCA,doc,ro,YDY,YDX) implicit none integer*4 IQ,I,J,K,L,N real*8 SAL0,SAL1,SAG0,SAG1,SA1,ram1,roa1,alpha(N,3,3), 1G(N,3,3),A(N,3,3,3),S180,D180,D,S90x,D90X,DX,S90Z,D90Z,DZ, 2YDY,YDX,AMU,BOHR,SpA3,beta2,gpisvejc,EXCA,betaa,betag,EPS, 1doc,ro,down SAL0=0.0d0 SAL1=0.0d0 SAG0=0.0d0 SAG1=0.0d0 SA1=0.0d0 DO 3 I=1,3 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 calculate degree of circularity down= 7.0d0*SAL1+ 1.0d0*SAL0 doc=5.0d0*( 1.0d0*SAL1- 1.0d0*SAL0) if(down.ne.0.0d0)doc=doc/down 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 if(YDX.gt.1.0d-9)then ro=YDY/YDX else ro=0.0d0 endif AMU=1822.0d0 BOHR=0.529177d0 SpA3=dsqrt(SAL0)/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 c 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 c roa2=8.0d0*(3.0d0*betag-betaa) c ram2=12.0d0*beta2 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) return end C ========================================== SUBROUTINE WRITETTTQ(N30,N7,NINT,ALPHA,A,G,alphav,E) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3), 1alphav(N30,3,3),E(*) c open(22,file='TOTAL.Q.TTT') WRITE(22,2003)NINT-N7+1 2003 FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2 ' mode e(cm-1) jx jy jz') DO 1 II=N7,NINT write(22,221)II,E(II) 221 format(i5,f12.2) DO 1 I=1,3 1 WRITE(22,2001)I,(ALPHA(II,I,J),J=1,3) 2001 FORMAT(I3,3F13.7) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' mode e(cm-1) jx(Bx) jy(By) jz(Bz)') DO 2 II=N7,NINT write(22,221)II,E(II) DO 2 I=1,3 2 WRITE(22,2001)I,(G(II,I,J),J=1,3) c last inner index magnetic WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' mode e(cm-1) kx ky kz') DO 3 II=N7,NINT write(22,221)II,E(II) DO 3 I=1,3 DO 3 J=1,3 3 WRITE(22,2002)I,J,(A(II,I,J,K),K=1,3) 2002 FORMAT(2I3,3F13.7) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' mode e(cm-1) vx vy vz') DO 11 II=N7,NINT write(22,221)II,E(II) DO 11 I=1,3 11 WRITE(22,2001)I,(ALPHAV(II,I,J),J=1,3) close(22) 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 wmtr(A,n0,n,st,ic,io) c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n0,n) character*(*) st write(io,*)st write(io,*) N1=1 1 N3=MIN(N1+4,N) WRITE(io,18)(I,I=N1,N3) 18 FORMAT(4X,5I14) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 WRITE(io,17)LN,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(io,*) return end subroutine nv(v) real*8 v(*),a,sp a=dsqrt(sp(v,v)) v(1)=v(1)/a v(2)=v(2)/a v(3)=v(3)/a return end function sp(a,b) real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine wrp0(ndim,M,x0,u,v,n,nx,ny,r,llog,ic,iw) implicit none integer*4 ndim,i,n,jx,kx,ia,ib,nx,ny,ix,istart,iz,nat,iw,ic,j real*8 M(ndim,ndim),x0(*),u(3),v(3),F(3,3),xa(3),sp, 1da,db,r(n,3),bohr,z,xi(3),rr(3),a(3),b(3),uu,vv logical llog character*80 fn character*10 s10,t10 real*8,allocatable::XR(:,:) bohr=0.529177d0 allocate(XR(24,ndim)) vv=dsqrt(sp(v,v)) uu=dsqrt(sp(u,u)) if(ic.eq.0.or.iw.eq.0)then open(9,file='PLANE.TXT') else write(s10,100)ic write(t10,100)iw 100 format(i10) do 1 i=1,len(s10) 1 if(s10(i:i).ne.' ')goto 111 111 do 8 j=1,len(t10) 8 if(t10(j:j).ne.' ')goto 222 222 open(9,file=s10(i:10)//'.'//t10(j:10)//'.PLANE.TXT') endif do 2 ia=1,nx da=dble(ia-1)*uu do 2 ib=1,ny db=dble(ib-1)*vv c take each plane point as n+1 particle: do 3 ix=1,3 3 rr(ix)=(x0(ix)+dble(ia-1)*u(ix)+dble(ib-1)*v(ix))/bohr call doX1(n,rr,r,XR) do 5 ix=1,3 do 5 kx=1,3 c kx ...for initial x,y and z-polarization: F(kx,ix)=0.0d0 if(iw.eq.0.and.ix.eq.kx)F(kx,ix)=1.0d0 do 5 i=1,n istart=24*(i-1)+1 c F = M . F(0) do 5 jx=1,24 c F(x) = X(x) P~ .F0 5 F(kx,ix)=F(kx,ix) 1+XR(ix,istart-1+jx)*M(istart-1+jx,istart-1+kx) z=(F(1,1)**2+F(1,2)**2+F(1,3)**2 1+ F(2,1)**2+F(2,2)**2+F(2,3)**2 1+ F(3,1)**2+F(3,2)**2+F(3,3)**2)/9.0d0 if(llog)then if(z.gt.0.0d0)then z=log(z)/log(10.0d0) else z=0.0d0 endif endif 2 write(9,601)da,db,z 601 format(2F12.2,G12.4) close(9) a=u b=v call nv(a) call nv(b) open(9,file='IRROA.OPT') open(19,file='ATOMS.TXT') read(9,*)n do 6 i=1,n read(9,900)fn 900 format(a80) open(29,file=fn) read(29,*) read(29,*)nat do 61 ia=1,nat read(29,*)iz,(xa(ix),ix=1,3) do 7 ix=1,3 7 xi(ix)=xa(ix)-x0(ix) 61 write(19,190)iz,sp(xi,a),sp(xi,b) 190 format(i3,2f12.3) 6 close(29) close(19) close(9) return end c ============================================================ subroutine wrcube(ndim,P,x0,n,nx,ny,nz,dx,dy,dz,r,llog) implicit none integer*4 ndim,i,n,jx,kx,ix,iy,iz,nx,ny,nz,istart,nstart, 1na,nat,ip,ia real*8 P(ndim,ndim),x0(*),F(3,24),xa(3),B(3,3), 1dx,dy,dz,r(n,3),bohr,x,y,z,zz,bz,ez logical llog character*80 fn real*8,allocatable::rr(:,:),XR(:,:),sz(:),bsz(:),esz(:) bohr=0.529177d0 allocate(rr(n+1,3),XR(24*(n+1),24*(n+1))) nstart=24*(n+1-1)+1 do 4 i=1,n do 4 ix=1,3 4 rr(i,ix)=r(i,ix) c determine number of atoms na=0 open(9,file='IRROA.OPT') read(9,*)n do 8 i=1,n read(9,900)fn 900 format(a80) open(29,file=fn) read(29,*) read(29,*)nat na=na+nat 8 close(29) close(9) open(39,file='CUB.CUB') open(40,file='BCUB.CUB') open(41,file='ECUB.CUB') write(39,901)na,x0(1)/bohr,x0(2)/bohr,x0(3)/bohr, 1nx,dx,0.,0., ny,0.,dy,0., nz,0.,0.,dz write(40,902)na,x0(1)/bohr,x0(2)/bohr,x0(3)/bohr, 1nx,dx,0.,0., ny,0.,dy,0., nz,0.,0.,dz write(41,903)na,x0(1)/bohr,x0(2)/bohr,x0(3)/bohr, 1nx,dx,0.,0., ny,0.,dy,0., nz,0.,0.,dz 901 format(' Electric field density',/, 1' SCF Total Density',/,i5,3f12.6,3(/,i5,3f12.6)) 902 format(' Magnetic field density',/, 1' SCF Total Density',/,i5,3f12.6,3(/,i5,3f12.6)) 903 format(' Gradient field density',/, 1' SCF Total Density',/,i5,3f12.6,3(/,i5,3f12.6)) c write atoms to cube open(9,file='IRROA.OPT') read(9,*)n do 6 i=1,n read(9,900)fn open(29,file=fn) read(29,*) read(29,*)nat do 61 ia=1,nat read(29,*)iz,(xa(ix),ix=1,3) write(39,190)iz,real(iz),(xa(ix)/bohr,ix=1,3) write(40,190)iz,real(iz),(xa(ix)/bohr,ix=1,3) 61 write(41,190)iz,real(iz),(xa(ix)/bohr,ix=1,3) 190 format(i5,4f12.6) 6 close(29) close(9) c write cube allocate(sz(nz),bsz(nz),esz(nz)) x=x0(1)/bohr-dx do 2 ix=1,nx write(6,*)ix,' of ',nx x=x+dx y=x0(2)/bohr-dy do 2 iy=1,ny y=y+dy z=x0(3)/bohr-dz do 21 iz=1,nz z=z+dz c take each plane point as n+1 particle: rr(n+1,1)=x rr(n+1,2)=y rr(n+1,3)=z call doX(n+1,rr,XR,0) do 5 ip=1,24 do 5 kx=1,3 c kx ...for initial x,y and z-polarization: F(kx,ip)=0.0d0 do 5 i=1,n istart=24*(i-1)+1 c F = X . P . F(0) do 5 jx=1,24 5 F(kx,ip)=F(kx,ip) 1+XR(nstart-1+ip,istart-1+jx)*P(istart-1+jx,istart-1+kx) do 9 ip=1,3 do 9 kx=1,3 c kx ...for initial x,y and z-polarization: B(kx,ip)=0.0d0 do 9 i=1,n istart=24*(i-1)+1 c B = X . P . B(0) do 9 jx=1,24 9 B(kx,ip)=B(kx,ip) 1+XR(nstart-1+ip,istart-1+jx)*P(istart-1+jx,istart-1+kx) c electric field: zz=(F(1,1)**2+F(1,2)**2+F(1,3)**2 1+ F(2,1)**2+F(2,2)**2+F(2,3)**2 1+ F(3,1)**2+F(3,2)**2+F(3,3)**2)/9.0d0 if(llog)then if(zz.gt.0.0d0)then zz=log(zz)/log(10.0d0) else zz=0.0d0 endif endif c magnetic field: bz=(B(1,1)**2+B(1,2)**2+B(1,3)**2 1+ B(2,1)**2+B(2,2)**2+B(2,3)**2 1+ B(3,1)**2+B(3,2)**2+B(3,3)**2)/9.0d0 if(llog)then if(bz.gt.0.0d0)then bz=log(bz)/log(10.0d0) else bz=0.0d0 endif endif c gradient field: ez=(F(1,13)**2+F(1,14)**2+F(1,15)**2 1 +F(1,16)**2+F(1,17)**2+F(1,18)**2 1+ F(2,13)**2+F(2,14)**2+F(2,15)**2 1 +F(2,16)**2+F(2,17)**2+F(2,18)**2 1+ F(3,13)**2+F(3,14)**2+F(3,15)**2 1 +F(3,16)**2+F(3,17)**2+F(3,18)**2)/18.0d0 if(llog)then if(ez.gt.0.0d0)then ez=log(ez)/log(10.0d0) else ez=0.0d0 endif endif sz(iz)=zz bsz(iz)=bz 21 esz(iz)=ez write(39,4000)( sz(iz),iz=1,nz) write(40,4000)(bsz(iz),iz=1,nz) 2 write(41,4000)(esz(iz),iz=1,nz) 4000 format(6E13.5) close(39) close(40) close(41) return end c ============================================================ subroutine readopt(ic,n,nm,m,ntw,ntmax,NW,nx,ny,ncub,iwr, 1nxc,nyc,nzc,ncubc,icubc,ndim, llog,ltab,lff,ldvm,lgauss, 1luseg,lusea,lcont,icen,jfreq,lcube,ltr,lave,debug, 1lonly,r,alpha,gp,a,alphad,gpd,ad,tw,WMAX,WMIN,W0,GAMMA,WA, 1temp,tol,tolw,x0,u,v,gw,x0c,dxc,dyc,dzc,DW,tempcm,cm, 1fn,ftw,ftwas,iiz,lasym,alphaas,alphadas) implicit none integer*4 ic,n,m,nm,ntw(*),ntmax,NW,nx,ny,ncub,iwr,icen, 1nxc,nyc,nzc,ncubc,icubc,ndim,i,iiz(*),jfreq logical llog,ltab,lff,ldvm,lgauss,luseg,lusea,lcont, 1lcube,ltr(*),lave,debug,lonly,lasym real*8 r(n,3),alpha(n,3,3),gp(n,3,3),a(n,3,3,3),alphaas(n,3,3),gw, 1alphad(n,nm,3,3),gpd(n,nm,3,3),ad(n,nm,3,3,3),tw(n,ntmax),WMAX, 1WMIN,W0,GAMMA,WA,temp,tol,tolw,x0(3),v(3),u(3),x0c(*), 1dxc,dyc,dzc,DW,tempcm,cm,alphadas(n,nm,3,3) character*80 fn(*),ftw(*),ft,s80,ftwas(*) cm=219474.0d0 c lasym - use asymmetric polarizability c lave average by roattion c lcube - make CUB.CUB with intensity c icen - make plane file PLANE.TXT with intensity c icen =-1 .. nothing c icen > 0 .. required Raaman centrum (molecule) c icen = 0 .. excitation frequency c jfreq .. required Raman frequency c ldvm - the discrete vibration model c lcont ... try the "continuum" frequency model c lff .. try also frequency combinations wi+wj c lgauss use gaussian, not lorentzian peaks c llog - print out logarithm of intensity c lonly only alpha-model c ltab .. write also ROA.TAB, inaddition to .prn files c luseg use Gp c lusea use A c debug debug output gw=0.0d0 icubc=0 lasym=.false. lave=.false. lcube=.false. icen=-1 jfreq=1 ldvm=.false. lcont=.false. lff=.false. lgauss=.false. llog=.false. lonly=.false. ltab=.false. luseg=.true. lusea=.true. debug=.false. c iwr - the amount of output iwr=0 NW=1751 c number of grid points: nx=50 ny=50 c group and frequency for the plane: ncub=1 tol=0.01d0 tolw=0.002d0 WMAX=3500.0d0 WMIN=0.0d0 W0=0.0d0 GAMMA=10.0d0 c excitation freqency in A: WA=5320.0d0 TEMP=300.0d0 open(9,file='IRROA.OPT',status='old') c number of particles/chromophores: read(9,*)n if(ic.eq.0)then close(9) return endif c GEOMETRIES: do 1 i=1,n ltr(i)=.true. read(9,900)fn(i) 900 format(a80) c geometries, get the centers only: 1 CALL READRX(r,n,fn,i,ic) c normal mode polarizability derivatives: c 1) read names and find maximum number of modes: m=0 nm=0 if(ic.eq.1)write(6,*)' group number of transition frequencies:' do 2 i=1,n read(9,900)ftw(i) open(15,file=ftw(i)) READ(15,*) c ntw: number of transition frequencies: READ(15,*)ntw(i) m=m+ntw(i) if(ic.eq.1)write(6,901)i,ntw(i),m 901 format(i5,' ',2i5) close(15) 2 if(ntw(i).gt.nm)nm=ntw(i) ntmax=nm+1 if(ic.eq.1)then write(6,*)m,' of total frequencies',' natmax:',ntmax close(9) return endif c 2) read derivatives/transition polarizabilities/ do 21 i=1,n 21 CALL READTTTQ(nm,n,alphad,Gpd,ad,ftw,tw,i,ntw,ntmax) write(6,*)'Pol. der. read.' if(lasym)then c read asymmetric polarizabilities, c 1) read names and find maximum number of modes: do 201 i=1,n 201 read(9,900)ftwas(i) c 2) read derivatives/transition polarizabilities/ do 211 i=1,n 211 CALL RDAS(nm,n,alphadas,ftwas,tw,i,ntw,ntmax) write(6,*)'Antisym. pol. der. read.' endif c polarizability alpha do 3 i=1,n iiz(i)=0 read(9,900)ft 3 call readpol(alpha,i,ft,n) write(6,*)'Pol. read.' if(lasym)then c asymmetric polarizability alpha do 31 i=1,n iiz(i)=0 read(9,900)ft 31 call readpol(alphaas,i,ft,n) write(6,*)'Antisymmetric polarizability read.' endif c G': do 4 i=1,n read(9,900)ft 4 call readpolg(gp,i,ft,n) write(6,*)'G. read.' c A: do 5 i=1,n read(9,900)ft 5 call readpola(a,i,ft,n) write(6,*)'A. read.' c plane origin in A: x0=0.0d0 c plane vectors v=0.0d0 u=0.0d0 66 read(9,920,end=888,err=888)s80 if(s80(1:4).eq.'AVER')read(9,*)lave if(s80(1:4).eq.'ASYM')read(9,*)lasym if(s80(1:4).eq.'ADER')read(9,*)lusea if(s80(1:4).eq.'GDER')read(9,*)luseg if(s80(1:4).eq.'DAMP')read(9,*)gw if(s80(1:3).eq.'TOL')read(9,*)tol if(s80(1:4).eq.'WMAX')read(9,*)WMAX if(s80(1:4).eq.'WTOL')read(9,*)tolw if(s80(1:4).eq.'WMIN')read(9,*)WMIN if(s80(1:4).eq.'TEMP')read(9,*)TEMP if(s80(1:4).eq.'TOLW')read(9,*)tolw if(s80(1:3).eq.'LFF')read(9,*)lff if(s80(1:4).eq.'CONT')read(9,*)lcont if(s80(1:4).eq.'LGAU')read(9,*)lgauss if(s80(1:4).eq.'WAEX')read(9,*)WA if(s80(1:4).eq.'LDVM')read(9,*)ldvm if(s80(1:4).eq.'DEBU')read(9,*)debug if(s80(1:4).eq.'ONLY')read(9,*)lonly if(s80(1:4).eq.'LLOG')read(9,*)llog if(s80(1:4).eq.'DELE')then read(9,*)i iiz(i)=1 endif if(s80(1:4).eq.'LTAB')read(9,*)ltab if(s80(1:3).eq.'IWR')read(9,*)iwr if(s80(1:4).eq.'PLAN')then read(9,*)icen,jfreq read(9,*)(x0(i),i=1,3) read(9,*)(u( i),i=1,3) read(9,*)(v( i),i=1,3) read(9,*)nx,ny endif if(s80(1:4).eq.'CUBE')then read(9,*)lcube read(9,*)(x0c(i),i=1,3) read(9,*)nxc,dxc read(9,*)nyc,dyc read(9,*)nzc,dzc read(9,*)ncubc endif if(s80(1:4).eq.'TRAN')then do 6 i=1,n 6 read(9,*)ltr(i) endif if(s80(1:2).eq.'W0')read(9,*)W0 if(s80(1:2).eq.'NW')read(9,*)NW if(s80(1:5).eq.'GAMMA')read(9,*)GAMMA goto 66 920 format(a80) 888 close(9) write(6,*)'option and tensors read' DW=(WMAX-WMIN)/dble(NW-1) tempcm=TEMP/1.44d0 ndim=24*n c 1.. 3 electric c 4.. 6 electric' c 7.. 9 mgnetic c 10..12 magnetic c 13..18 quadrupole c 19..24 quadrupole' return end c ============================================================ subroutine wropt(n,ndim,ntmax,NW,iwr,WMIN,WMAX,DW,GAMMA,W0,tolw, 1TEMP,wa,tol,gw,lave,lusea,luseg,lff,lcont,ldvm,lonly, 1lasym) implicit none integer*4 n,ndim,ntmax,NW,iwr real*8 WMIN,WMAX,DW,GAMMA,W0,tolw,TEMP,wa,tol,gw logical lave,lusea,luseg,lff,lcont,ldvm,lonly,lasym write(6,6000)n,ndim,ntmax,NW,iwr,WMIN,WMAX,DW,GAMMA,W0,tolw,TEMP, 1wa,tol,gw,lave,lusea,luseg,lff,lcont,ldvm,lonly,lasym 6000 format(' Number of mols. :',i20,/, 2 ' Total dimension :',i20,/, 2 ' max exc :',i20,/, 2 ' Nw :',i20,/, 2 ' iwr :',i20,/, 2 ' Wmin :',F20.2, ' cm-1',/, 2 ' Wmax :',F20.2, ' cm-1',/, 2 ' DW :',F20.2, ' cm-1',/, 2 ' Gamma :',F20.2, ' cm-1',/, 2 ' Wexc(relative) :',F20.2, ' cm-1',/, 2 ' tolw :',f20.4, ' cm-1',/, 2 ' Temperature :',F20.2, ' K',/, 2 ' Wexc :',f20.2, ' A',/, 2 ' tol :',G20.6,/, 2 ' matrix damping :',G20.6,/, 1 ' Averaging :',l20,/, 1 ' Use A :',l20,/, 2 ' Use G :',l20,/, 2 ' lff :',l20,/, 2 ' continuum model :',l20,/, 2 ' discr. vib. model :',l20,/, 2 ' only alpha :',l20,/, 2 ' assymetric alpha :',l20,/,/) return end c ============================================================ subroutine mm(n,a,b,c) c matrix multiplication A = B x C implicit none integer*4 n,i,j,k real*8 a(n,n),b(n,n),c(n,n) a=0.0d0 do 1 i=1,n do 1 j=1,n do 1 k=1,n 1 a(i,j)=a(i,j) + b(i,k)*c(k,j) return end c ============================================================ subroutine pci(n,t) implicit none integer*4 n,i,j,ii,jj real*8 t(24*n,24*n) do 1 i=1,n do 1 ii=1,24 do 1 jj=1,24 do 1 j=1,n if(j.ne.i)then t(24*(i-1)+ii,24*(i-1)+jj)=t(24*(i-1)+ii,24*(i-1)+jj)+ 1 t(24*(i-1)+ii,24*(j-1)+jj) t(24*(i-1)+ii,24*(j-1)+jj)=0.0d0 endif 1 continue return end c ============================================================ subroutine iizd(n,nm,iiz,alphad,gpd,ad) c delete derivatives of desired particle implicit none integer*4 n,nm,iiz(*),i,iq,ix,jx,kx real*8 alphad(n,nm,3,3),gpd(n,nm,3,3), ad(n,nm,3,3,3) do 1 i=1,n if(iiz(i).eq.1)then write(6,*)'deleting derivatives of particle ',i do 2 iq=1,nm do 2 ix=1,3 do 2 jx=1,3 alphad(i,iq,ix,jx)=0.0d0 gpd(i,iq,ix,jx)=0.0d0 do 2 kx=1,3 2 ad(i,iq,ix,jx,kx)=0.0d0 endif 1 continue return end c ============================================================ subroutine getpi(PI,ndim,it,i,Ptr,ntmax) implicit none integer*4 ndim,it,i,ii,jj,ntmax real*8 PI(ndim,ndim),Ptr(ntmax,ndim,ndim) PI=0.0d0 do 5 ii=24*(i-1)+1,24*i do 5 jj=24*(i-1)+1,24*i 5 PI(ii,jj)=Ptr(it,ii,jj) return end c ============================================================ subroutine ctp(n,nm,ntmax,ndim,Ptr,alphad,gpd,ad,lonly,ntw) implicit none integer*4 i,n,istart,it,ix,ii,jx,ntmax,ndim,nm,kx,ntw(*) real*8 Ptr(ntmax,ndim,ndim),alphad(n,nm,3,3),gpd(n,nm,3,3), 1ad(n,nm,3,3,3) logical lonly do 71 i=1,n istart=24*(i-1)+1 do 71 it=1,ntw(i) do 71 ix=1,3 ii=0 do 71 jx=1,3 Ptr(it,istart+ix-1,istart+jx-1)= alphad(i,it,ix,jx) if(.not.lonly) 1Ptr(it,istart+ix+2,istart+jx+2)= alphad(i,it,ix,jx) Ptr(it,istart+ix-1,istart+jx+8)= gpd(i,it,ix,jx) Ptr(it,istart+jx+8,istart+ix-1)= gpd(i,it,ix,jx) Ptr(it,istart+ix+2,istart+jx+5)= -gpd(i,it,ix,jx) Ptr(it,istart+jx+5,istart+ix+2)= -gpd(i,it,ix,jx) do 71 kx=jx,3 ii=ii+1 Ptr(it,istart+ix- 1,istart+ii+11)=ad(i,it,ix,jx,kx)/3.0d0 Ptr(it,istart+ii+11,istart+ix- 1)=ad(i,it,ix,jx,kx)/3.0d0 Ptr(it,istart+ix+ 2,istart+ii+17)=ad(i,it,ix,jx,kx)/3.0d0 71 Ptr(it,istart+ii+17,istart+ix+ 2)=ad(i,it,ix,jx,kx)/3.0d0 write(6,*)' Transition polarizabilities added' c c E E. B B. EE EE. c ---------------------- c u | a 0 0 G A/3 0 c u. | 0 a -G 0 0 A/3 c m | 0 -Gt 0 0 0 0 c m. | Gt 0 0 0 0 0 c O | A/3 0 0 0 0 0 c O. | 0 A/3 0 0 0 0 c ************* Calculate general polarizability tensor: ********* return end subroutine ct0(n,ndim,P0,alpha,gp,a,lonly) implicit none integer*4 i,n,ix,jx,kx,ndim,istart,ii,jj real*8 P0(ndim,ndim),P0i(24,24),alpha(n,3,3),gp(n,3,3),a(n,3,3,3) logical lonly P0=0.0d0 do 74 i=1,n P0i=0.0d0 do 741 ix=1,3 ii=0 do 741 jx=1,3 P0i(1+ix-1,1+jx-1)= alpha(i,ix,jx) if(.not.lonly) 1P0i(1+ix+2,1+jx+2)= alpha(i,ix,jx) P0i(1+ix-1,1+jx+8)= gp(i,ix,jx) P0i(1+jx+8,1+ix-1)= gp(i,ix,jx) P0i(1+ix+2,1+jx+5)= -gp(i,ix,jx) P0i(1+jx+5,1+ix+2)= -gp(i,ix,jx) do 741 kx=jx,3 ii=ii+1 P0i(1+ix- 1,1+ii+11)=a(i,ix,jx,kx)/3.0d0 P0i(1+ii+11,1+ix- 1)=a(i,ix,jx,kx)/3.0d0 P0i(1+ix+ 2,1+ii+17)=a(i,ix,jx,kx)/3.0d0 741 P0i(1+ii+17,1+ix+ 2)=a(i,ix,jx,kx)/3.0d0 istart=24*(i-1)+1 do 74 ii=1,24 do 74 jj=1,24 74 P0(istart-1+ii,istart-1+jj)=P0i(ii,jj) return end subroutine ftp(i,j,ndim,ntmax,Pt,Ptr) c put transition polarizability of group i and frequency j into Pt implicit none integer*4 i,ix,jx,ndim,ii,jj,istart,ntmax,j real*8 Pt(ndim,ndim),Ptr(ntmax,ndim,ndim) Pt=0.0d0 istart=24*(i-1)+1 do 711 ix=1,24 ii=istart+ix-1 do 711 jx=1,24 jj=istart+jx-1 711 Pt(ii,jj)=Pt(ii,jj)+Ptr(j,ii,jj) end c =============================== SUBROUTINE WRITEPOL(APOL,f,fr) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) real*8 APOL(*) character*(*) f character*1 j j=f(1:1) OPEN(90,FILE=f) if(j.eq.'A')then write(90,701)fr 701 format(' A, dipole-quadrupole polarizability, w = ',f11.6) write(90,902) do 1 IZ=1,3 1 write(90,900)((3.0d0/2.0d0* 1 APOL(IY+(IX-1)*3+(IZ-1)*9),IY=1,IX),IX=1,3) endif if(j.eq.'P')then write(90,702)fr 702 format(' Polarizability, w = ',f11.6) write(90,902) write(90,900)((APOL(IY+(IX-1)*3),IY=1,IX),IX=1,3) write(90,900)(APOL(1)+APOL(5)+APOL(9))/3.0d0 endif if(j.eq.'G')then write(90,704)fr 704 format(' G tensor, w = ',f11.6) write(90,901) write(90,900)((APOL(IY+(IX-1)*3),IY=1,3),IX=1,3) endif 901 format( 1' XX XY XZ YX', 2' YY YZ ZX ZY', 3' ZZ') 902 format( 1' XX XY YY XZ', 2' YZ ZZ') 900 format(9G14.6) close(90) return end