PROGRAM IRROA IMPLICIT none integer*4 n,i,ix,jx,kx,ii,istart,kk,nm,NW,k,ntt,ic,j, 1jt,kt,ktend,l,lt,ltend,lend,m,iwr,nx,ny,ncub,icub, 1jj,ifi,it,iter,ns,nf,ict,ndim,iw,ntmax,nxc,nyc,nzc, 1ncubc,icubc,IERR,no,nt real*8 r0,dt,df,om,th,fi,ds1,ds2,an1,an2,omold,ds3,tol,sum, 1WMIN,WMAX,DW,GAMMA,sf,W0,TEMP,tempcm,w,tolw,wa,x0(3),v(3), 1dx,dy,u(3),dxc,dyc,dzc,x0c(3),cm,gw character*80,allocatable:: fn(:),ftw(:),fnf(:),ftwas(:) logical,allocatable:: ltr(:) logical lave,luseg,lusea,inversion,debug,lcont,lff,lgauss, 1ldvm,lplane,llog,lcube,ltab,ltdc,lonly,lasym integer*4 ,allocatable::ntw(:),iiz(:) real*8,allocatable:: alpha(:,:,:),gp(:,:,:),a(:,:,:,:), 1alphad(:,:,:,:),gpd(:,:,:,:),ad(:,:,:,:,:),r(:,:),X(:,:), 1alphar(:,:,:),gpr(:,:,:),ar(:,:,:,:),SIT(:,:),Ptr(:,:,:), 1SI(:,:),SIOLD(:,:),wl(:),EXP0(:,:),P0(:,:),P0i(:,:),XP0(:,:), 1tw(:,:),Pt(:,:,:),h(:,:),rr(:,:),d(:),e(:),ee(:),c(:,:), 1alphaas(:,:,:),alphadas(:,:,:,:) n=1 nm=1 ntmax=1 c initial allocation: allocate(r(n,3),fn(n),ftw(n),ntw(n),ltr(n),fnf(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,icub,iwr, 1nxc,nyc,nzc,ncubc,icubc,ndim, ltdc,llog,ltab,lff,ldvm,lgauss, 1luseg,lusea,lcont,lplane,lcube,ltr,lave,inversion,debug, 1lonly,r,alpha,gp,a,alphad,gpd,ad,tw,WMAX,WMIN,W0,GAMMA,WA, 1temp,tol,tolw,x0,v,u,dx,dy,gw,x0c,dxc,dyc,dzc,DW,tempcm,cm, 1fn,ftw,ftwas,fnf,iiz,lasym,alphaas,alphadas) if(ic.eq.0)then write(6,*)n,' groups' deallocate(r,fn,ftw,ntw,ltr,fnf,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),fnf(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,fnf,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),fnf(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,inversion,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 call vz(gpd,n*nm*9) call vz(ad,n*nm*27) call vz(gp,n*9) call vz(a,n*27) write(6,*)' G, A erased' endif allocate(X(ndim,ndim),Ptr(ntmax,ndim,ndim)) call vz(Ptr,ntmax*ndim*ndim) write(6,*)' X,P allocated' c ************* Calculate general distance tensor: ********* call doX(n,r,X) write(6,*)' Distance tensor done' if(iwr.gt.1)call wmtr(X,ndim,ndim,'X-tensor',1,6) c ************* Calculate general polarizability tensor: ********* c (transition polarizabilities) 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: ********* C P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0 c P0:ordinary (non-transition) polarizability allocate(P0(ndim,ndim),P0i(24,24)) call vz(P0,ndim*ndim) do 74 i=1,n call vz(P0i,24*24) 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 742 ii=1,24 do 742 jj=1,24 742 P0(istart-1+ii,istart-1+jj)=P0i(ii,jj) c (ordinary polarizabilities as ntw+1^th component): ntw(i)=ntw(i)+1 tw(i,ntw(i))=0.0d0 do 743 ii=1,24 do 743 jj=1,24 743 Ptr(ntw(i),istart-1+ii,istart-1+jj)=P0i(ii,jj) 74 continue write(6,*)' Zero frequency polarizability added' C P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0P0 c EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0 c (E - X.P0): allocate(EXP0(ndim,ndim),XP0(ndim,ndim)) call mm(ndim,XP0,X,P0) do 75 i=1,ndim do 751 j=1,ndim 751 EXP0(i,j)=-XP0(i,j) 75 EXP0(i,i)=1.0d0+EXP0(i,i) c stop c EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0EXP0 c how many total t-frequencies wi,wi+wj: ntt=0 do 22 ic=1,2 c ic=1 ... only count c ic=2 ... add to list if(ic.eq.2)then allocate(wl(ntt)) ntt=0 endif c look at wi do 17 i=1,n do 17 it=1,ntw(i) w=tw(i,it) c c was it already among wk: do 18 k=1,i ktend=ntw(k) if(k.eq.i)ktend=it-1 do 18 kt=1,ktend 18 if(dabs(w-tw(k,kt)).lt.tolw)goto 17 c w is unique: ntt=ntt+1 if(ic.eq.2)wl(ntt)=w 17 continue c look at wi+wj: if(lff)then do 19 i=1,n do 19 it=1,ntw(i) do 19 j=1,i-1 do 19 jt=1,ntw(j) w=tw(i,it)+tw(j,jt) if(w.gt.WMAX)goto 19 c c was it already among wk: do 20 k=1,n do 20 kt=1,ntw(k) 20 if(dabs(w-tw(k,kt)).lt.tolw)goto 19 c was it already among wk+wl: do 26 k=1,i ktend=ntw(k) if(k.eq.i)ktend=it do 26 kt=1,ktend lend=k-1 if(k.eq.i.and.kt.eq.it)lend=j do 26 l=1,lend ltend=ntw(l) if(k.eq.i.and.kt.eq.it.and.l.eq.j)ltend=jt-1 do 26 lt=1,ltend 26 if(dabs(w-tw(k,kt)-tw(l,lt)).lt.tolw)goto 19 c w is unique ntt=ntt+1 if(ic.eq.2)wl(ntt)=w 19 continue endif 22 continue write(6,6012)ntt 6012 format(i9,' unique frequencies:') write(6,6013)(wl(i),i=1,ntt) 6013 format(8F10.2) write(6,*) c for debug purposes, make total polarizability for W <- 0: c debugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebug if(debug)then allocate 1 (alphar(NW,3,3),gpr(NW,3,3),ar(NW,3,3,3),Pt(NW,ndim,ndim)) do 77 iw=1,NW w=WMIN+dble(iw-1)*DW do 77 i=1,n istart=24*(i-1)+1 do 77 ii=istart,istart+23 do 77 jj=istart,istart+23 sum=0.0d0 do 744 it=1,ntw(i) 744 sum=sum+sf(GAMMA,lgauss,w,tw(i,it))*Ptr(it,ii,jj) 77 Pt(iw,ii,jj)=sum call extr(Pt,alphar,gpr,ar,NW,ndim,r,n) call DORAMAN(NW,alphar,Gpr,Ar,wa,0,WMIN,DW,tempcm) deallocate(alphar,gpr,ar,Pt) c explicit list of T-frequencies: allocate 1 (alphar(ntt,3,3),gpr(ntt,3,3),ar(ntt,3,3,3),Pt(ntt,ndim,ndim)) call vz(Pt,ntt*ndim*ndim) do 24 i=1,n istart=24*(i-1)+1 do 24 it=1,ntw(i) w=tw(i,it) do 25 iw=1,ntt if(dabs(w-wl(iw)).lt.tolw)then do 23 ii=istart,istart+23 do 23 jj=istart,istart+23 23 Pt(iw,ii,jj)=Pt(iw,ii,jj)+Ptr(it,ii,jj) endif 25 continue 24 continue call extr(Pt,alphar,gpr,ar,ntt,ndim,r,n) call DORAMAS(NW,GAMMA,ntt,alphar,Gpr,Ar,wa,WMIN,DW, 1 tempcm,wl,lgauss,ltab) call WRITETTTQ(ntt,1,ntt,alphar,Ar,Gpr,alphar,wl) deallocate(alphar,gpr,ar,Pt) call report('debug stop') endif c debugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebugdebug if(ltdc)then allocate(d(3),e(3)) no=0 do 27 i=1,n call rdf(fnf(i),nt,0,e,d) 27 no=no+nt write(6,*)no,' oscillators' deallocate(d,e) allocate(d(no*3),e(no),rr(no,3)) istart=1 do 28 i=1,n call rdf(fnf(i),nt,1,e(istart),d(3*(istart-1)+1)) do 29 j=1,nt rr(istart+j-1,1)=r(i,1) rr(istart+j-1,2)=r(i,2) 29 rr(istart+j-1,3)=r(i,3) 28 istart=istart+nt call wrdip(no,e,d) allocate(h(no,no),ee(no),c(no,no)) call setv(no,h,e,rr,d) WRITE(6,*)' DIAGONALIZING THE HAMILTONIAN' IERR=0 CALL TRED12(no,h,c,ee,2,IERR) IF(IERR.NE.0)call report('DIAGONALIZATION ERROR') WRITE(6,*)' DIAGONALIZATION SUCCESSFUL' endif if(lcont)then c get redressed polarizabilities, common origin, continuum of freqs write(6,*)'Continuum frequencies' allocate 1 (alphar(NW,3,3),gpr(NW,3,3),ar(NW,3,3,3),Pt(NW,ndim,ndim)) call MDP(n,ndim,x,Ptr,Pt,inversion,NW,WMIN,DW,ntw,ntmax, 1 GAMMA,tw,lgauss) call extr(Pt,alphar,gpr,ar,NW,ndim,r,n) call DORAMAN(NW,alphar,Gpr,Ar,wa,0,WMIN,DW,tempcm) else write(6,*)'Discrete frequencies' c explicit list of T-frequencies: if(ldvm)then write(6,*)'DVM' deallocate(wl) allocate(alphar(m,3,3),gpr(m,3,3),ar(m,3,3,3),wl(m), 1 Pt(m,ndim,ndim)) iw=0 do 101 i=1,n do 101 it=1,ntw(i)-1 iw=iw+1 101 wl(iw)=tw(i,it) call vz(Pt,m*ndim*ndim) call DVM(m,n,ndim,x,Ptr,Pt,ntw,ntmax,EXP0,P0,iwr,wl,gw, 1 XP0) call extrc(Pt,alphar,gpr,ar,m,ndim,r,n,ntw) if(ltdc)then write(6,*)' m: ',m write(6,*)' no: ',no call extd(alphar,gpr,ar,m,c,no) write(6,*)'tdc frequencies / cm-1:' do 102 i=1,no wl(i)=ee(i)*cm 102 write(6,605)i,wl(i) 605 format(i5,f12.4) else no=1 endif call DORAMAS(NW,GAMMA,m,alphar,Gpr,Ar,wa,WMIN,DW, 1 tempcm,wl,lgauss,ltab) call WRITETTTQ(m,1,m,alphar,Ar,Gpr,alphar,wl) else allocate(alphar(ntt,3,3),gpr(ntt,3,3),ar(ntt,3,3,3), 1 Pt(ntt,ndim,ndim)) call vz(Pt,ntt*ndim*ndim) call MDE(ntt,n,ndim,x,Ptr,Pt,ntw,ntmax,tw,wl,tolw,inversion) call extr(Pt,alphar,gpr,ar,ntt,ndim,r,n) call DORAMAS(NW,GAMMA,ntt,alphar,Gpr,Ar,wa,WMIN,DW, 1 tempcm,wl,lgauss,ltab) call WRITETTTQ(ntt,1,ntt,alphar,Ar,Gpr,alphar,wl) endif endif if(lplane)call PLAN(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0,u,v,nx,ny,ncub,icub,dx,dy,r,llog,c,no,ltdc,ntw) if(lcube)call CUBE(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0c,nxc,nyc,nzc,ncubc,icubc,dxc,dyc,dzc,r,llog,c,no,ltdc,ntw) if(lave)then r0=dsqrt(r(1,1)**2+r(1,2)**2+r(1,3)**2) allocate(SIT(2,NW),SI(2,NW),SIOLD(2,NW)) do 15 ii=1,NW do 15 k=1,2 SI(k,ii)=0.0d0 SIT(k,ii)=0.0d0 15 SIOLD(k,ii)=0.0d0 ns=4 omold=0.0d0 write(6,*)'averaging position of the first particle' iter=0 99 iter=iter+1 dt= 4.0d0*atan(1.0d0)/dble(ns) do 151 ii=1,NW 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 r(1,1)=r0*sin(th)*cos(fi) r(1,2)=r0*sin(th)*sin(fi) call doX(n,r,X) call MDP(n,ndim,x,Ptr,Pt,inversion,NW,WMIN,DW,ntw,ntmax, 1 GAMMA,tw,lgauss) call DORAMAN(NW,alphar,Gpr,Ar,wa,1,WMIN,DW,tempcm) om=om+dt*sin(th)*df do 14 jj=1,2 do 14 kk=1,NW 14 SI(jj,kk)=SI(jj,kk)+SIT(jj,kk)*dt*sin(th)*df c errors of Raman and ROA spectra: ds1=0.0d0 ds2=0.0d0 an1=0.0d0 an2=0.0d0 do 16 kk=1,NW 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(NW,alphar,Gpr,Ar,wa,2,WMIN,DW,tempcm) else goto 99 endif else write(6,*)'converged' call DORAMAN(NW,alphar,Gpr,Ar,wa,2,WMIN,DW,tempcm) endif else write(6,*)'averaging switched off' endif 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 call vz(ttt,81) 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 call vz(ttt,27) 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 call vz(tt,9) 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(N,GAMMA,ntt,alpha,G,A,EXCA,WMIN,DW,tempcm, 1wl,lgauss,ltab) implicit none integer*4 N,IQ,I,ntt real*8 alpha(ntt,3,3),G(ntt,3,3),A(ntt,3,3,3),EXCA,ECM,sf, 3roa1,ram1,WMIN,DW,tempcm,const,wl(*),GAMMA,t,doc,ro,YDY,YDX real*8,allocatable::s(:,:) logical lgauss,ltab allocate(s(2,N)) call vz(s,2*N) if(ltab)then open(20,file='ROA.TAB') write(20,200) 200 format('IRROA model',/, 1 ' 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,N t=sf(GAMMA,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,N 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 DORAMAN(N,alpha,G,A,EXCA,ic,WMIN,DW,tempcm) implicit none integer*4 N,IQ,ic real*8 alpha(N,3,3),G(N,3,3),A(N,3,3,3),EXCA,ECM, 3roa1,ram1,WMIN,DW,tempcm,const,doc,ro,YDX,YDY if(ic.eq.0)then OPEN(41,FILE='ram.prn') OPEN(42,FILE='roa.prn') OPEN(43,FILE='doc.prn') else call report('not implemented') endif DO 1 IQ=1,N ECM=WMIN+DW*dble(IQ-1) call rr(ram1,roa1,IQ,alpha,G,A,N,EXCA,doc,ro,YDY,YDX) if(dabs(ECM).lt.0.1d0)ECM=0.1d0 const=1.0d0/ECM/(1.0d0-exp(-ECM/tempcm)) WRITE(41,3000)ECM,ram1*const WRITE(42,3000)ECM,roa1*const 1 WRITE(43,3000)ECM,doc 3000 FORMAT(f9.2,g12.4) close(41) close(42) close(43) RETURN END c ============================================================ 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) call vz(x,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 MDE(ntt,n,ndim,x,Ptr,Pt,ntw,ntmax,tw,wl,tolw,inversion) implicit none integer*4 ndim,i,j,ii,it,jj,jt,kk,n,ntw(*),iw,ntmax,ntt,j1 real*8 x(n*24,n*24),Ptr(ntmax,ndim,ndim),w,tolw,wj, 1Pt(ntt,ndim,ndim),sum,tw(n,ntmax),wi,wl(*) real*8,allocatable::tem(:,:,:),a(:,:),ai(:,:) logical inversion c make X.P, associated transition frequency list is with j-index: allocate(tem(ntmax,ndim,ndim)) call vz(tem,ntmax*ndim*ndim) if(inversion)then allocate(a(ndim,ndim),ai(ndim,ndim)) do 1002 j1=1,n do 1002 it=1,ntw(j1) do 1003 i =1,n do 1003 ii=24*(i-1)+1,24*i do 1003 j =1,n do 1003 jj=24*(j-1)+1,24*j sum=0.0d0 do 10012 kk=24*(j-1)+1,24*j 10012 sum=sum+x(ii,kk)*Ptr(it,kk,jj) a(ii,jj)=-sum 1003 if(ii.eq.jj)a(ii,jj)=a(ii,jj)+1.0d0 call INV(a,ai,ndim,1.0d-25) do 1002 i =1,n do 1002 ii=24*(i-1)+1,24*i do 1002 jj=24*(j1-1)+1,24*j1 1002 tem(it,ii,jj)=ai(ii,jj) deallocate(a,ai) else do 1001 i =1,n do 1001 ii=24*(i-1)+1,24*i do 1001 j =1,n do 1001 jj=24*(j-1)+1,24*j do 1001 it=1,ntw(j) sum=0.0d0 do 10011 kk=24*(j-1)+1,24*j 10011 sum=sum+x(ii,kk)*Ptr(it,kk,jj) if(it.eq.ntw(j).and.ii.eq.jj)sum=sum+1.0d0 1001 tem(it,ii,jj)=sum c perturbation approach, (E-X.P)^-1 ~ (E+X.P) endif write(6,*)' (E-X.P)^(-1) made' c make total polarizability Pt = P(E-X.P)^(-1) for W <- 0: c = Ptr_w w' temi_w'_0= P . (1+XP)= c P.tem: call vz(Pt,ntt*ndim*ndim) do 74 i=1,n do 74 it=1,ntw(i) wi=tw(i,it) do 74 j=1,n do 74 jt=1,ntw(j) wj=tw(j,jt) do 74 iw=1,ntt w=wl(iw) if(dabs(w-wi-wj).lt.tolw)then do 76 ii=24*(i-1)+1,24*i do 76 jj=24*(j-1)+1,24*j c P.tem: sum=0.0d0 do 761 kk=24*(i-1)+1,24*i c actually kk=1,ndim, but Ptr diagonal in i 761 sum=sum+Ptr(it,ii,kk)*tem(jt,kk,jj) 76 Pt(iw,ii,jj)=Pt(iw,ii,jj)+sum endif 74 continue write(6,*)' P(E-X.P)^(-1) made' return end c ============================================================ subroutine PLAN(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0,u,v,nx,ny,ncub,icub,dx,dy,r,llog,c,no,ltdc,ntw) implicit none integer*4 ndim,i,ii,it,jj,n,ntmax,iw, 1nx,ny,ncub,icub,no,ntw(*) real*8 x(n*24,n*24),Ptr(ntmax,ndim,ndim),r(n,3),u(3), 1EXP0(ndim,ndim),P0(ndim,ndim),x0(3),v(3),dx,dy,bohr, 1c(no,no) logical llog,ltdc real*8,allocatable::EPXM(:,:),t1(:,:),t2(:,:),PI(:,:), 1t3(:,:) write(6,600)nx,ny,dx,dy,x0,v,u,ncub,icub 600 format(' Constructing plane ',i5,' x ',i5,/, 1 ' Grid: ',f6.3,' x ',f6.3,' A',/, 1 ' Origin: ',3f12.3,' A',/, 2 ' Point2 ',3f12.3,' A',/, 2 ' Point3 ',3f12.3,' A',/, 3 ' Group ',i5,/, 4 ' Frequency number ',i5,/) bohr=0.529177d0 dx=dx/bohr dy=dy/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-P0.X)-1: call mm(ndim,t1,P0,EPXM) if(icub.eq.0)then c excitation frequency: t3=t1 else c vibrational frequency: c t2=E+(E-P0.X)-1 . P0 X: call mm(ndim,t2,t1,X) do 3 i=1,ndim 3 t2(i,i)=t2(i,i) + 1.0d0 c total polarizability Pt = (E-P0.X)^(-1)P (E+X*(E-XP0)-1.P0) allocate(PI(ndim,ndim)) call vz(PI,ndim*ndim) if(ltdc)then iw=0 do 7 i=1,n do 7 it=1,ntw(i)-1 iw=iw+1 do 7 ii=24*(i-1)+1,24*i do 7 jj=24*(i-1)+1,24*i 7 PI(ii,jj)=PI(ii,jj)+c(icub,iw)*Ptr(it,ii,jj) if(iw.ne.no)call report('iw<>no') else call getpi(PI,ndim,icub,ncub,Ptr,ntmax) endif 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 wrplane(ndim,t3,x0,u,v,n,nx,ny,dx,dy,r,llog) return end c ============================================================ subroutine CUBE(n,ndim,x,Ptr,ntmax,EXP0,P0, 1x0,nx,ny,nz,ncub,icub,dx,dy,dz,r,llog,c,no,ltdc,ntw) implicit none integer*4 ndim,i,ii,it,jj,n,ntmax, 1nx,ny,nz,ncub,icub,no,ntw(*),iw 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, 2c(no,no) logical llog,ltdc 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)) call vz(PI,ndim*ndim) if(ltdc)then iw=0 do 7 i=1,n do 7 it=1,ntw(i)-1 iw=iw+1 do 7 ii=24*(i-1)+1,24*i do 7 jj=24*(i-1)+1,24*i 7 PI(ii,jj)=PI(ii,jj)+c(icub,iw)*Ptr(it,ii,jj) if(iw.ne.no)call report('iw<>no') else call getpi(PI,ndim,icub,ncub,Ptr,ntmax) endif 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 DVM(m,n,ndim,x,Ptr,Pt,ntw,ntmax,EXP0,P0,iwr,wl,gw, 1XP0) implicit none integer*4 ndim,i,j,ii,it,jj,n,ntw(*),iw,ntmax,m,iwr, 1istart real*8 x(n*24,n*24),Ptr(ntmax,ndim,ndim),Pt(m,ndim,ndim), 1EXP0(ndim,ndim),P0(ndim,ndim),wl(*),gw,XP0(ndim,ndim) real*8,allocatable::EPXM(:,:),t1(:,:),t2(:,:),PI(:,:),a(:,:), 1t11(:,:),t22(:,:) if(iwr.eq.1)then allocate(a(24,24)) do 101 i=1,n write(6,*)'group ',i istart=24*(i-1)+1 do 102 ii=1,24 do 102 jj=1,24 102 a(ii,jj)=P0(istart-1+ii,istart-1+jj) 101 call wmtr(a,24,24,'P0 - old',1,6) iw=0 do 103 i=1,n write(6,*)'group ',i istart=24*(i-1)+1 do 103 it=1,ntw(i)-1 iw=iw+1 write(6,6001)it,wl(iw) 6001 format('frequency ',i4,f12.6) do 104 ii=1,24 do 104 jj=1,24 104 a(ii,jj)=Ptr(it,istart-1+ii,istart-1+jj) 103 call wmtr(a,24,24,'PI - old',1,6) endif allocate(EPXM(ndim,ndim)) if(gw.ne.0.0d0)then c make (E-XP0)t [(E-XP0)(E-XP0)t + G^2 XP0 XP0t]-1 call INVG(gw,EXP0,XP0,EPXM,ndim,1.0d-25) else c make (E-XP0)-1 call INV(EXP0,EPXM,ndim,1.0d-25) endif allocate(t1(ndim,ndim)) if(iwr.eq.1)then call wmtr(EXP0,ndim,ndim,'EXP0',1,6) call wmtr(EPXM,ndim,ndim,'EXP0-1',1,6) call vz(t1,ndim*ndim) do 7 i=1,ndim do 7 j=1,ndim do 7 ii=1,ndim 7 t1(i,j)=t1(i,j)+EXP0(i,ii)*EPXM(ii,j) if(iwr.gt.1)call wmtr(t1,ndim,ndim,'EXP0 x EXP0-1',1,6) endif allocate(PI(ndim,ndim),t11(ndim,ndim), 1t2(ndim,ndim),t22(ndim,ndim)) c c P0(E-P0.X)-1 - new effective diagonal polarizabilities: call smooth(EPXM,ndim*ndim,1.0d-100) call smooth(P0,ndim*ndim,1.0d-100) c t1 = P0 .(E-XP0.X)-1 call mm(ndim,t1,P0,EPXM) if(iwr.eq.1)then t2=t1 call pci(n,t2) do 105 i=1,n write(6,*)'group ',i istart=24*(i-1)+1 do 106 ii=1,24 do 106 jj=1,24 106 a(ii,jj)=t2(istart-1+ii,istart-1+jj) 105 call wmtr(a,24,24,'P0 - new',1,6) endif c t2 = P0 .(E-XP0.X)-1 X : call mm(ndim,t2,t1,X) c total polarizability: call vz(Pt,m*ndim*ndim) c t11:first part Pt = (E-P0.X)^(-1)P , for each frequency iw=0 do 41 i=1,n do 41 it=1,ntw(i)-1 if(iwr.eq.1)write(6,*)' group ',i,' freq',it iw=iw+1 call getpi(PI,ndim,it,i,Ptr,ntmax) call smooth(PI,ndim*ndim,1.0d-100) c t11 = PI (E-XP0)-1: call mm(ndim,t11,PI,EPXM) if(iwr.eq.1)then call tk1(n,i,ndim,a,t11) call wmtr(a,24,24,'PI - part1 - new',1,6) endif do 41 ii=1,ndim do 41 jj=1,ndim 41 Pt(iw,ii,jj)=t11(ii,jj) c second part Pt = (E-P0.X)^(-1)P *X*(E-XP0)-1.P0 iw=0 do 4 i=1,n do 4 it=1,ntw(i)-1 if(iwr.eq.1)write(6,*)' group ',i,' freq',it iw=iw+1 call getpi(PI,ndim,it,i,Ptr,ntmax) c t11 = PI (E-XP0)-1: call mm(ndim,t1,PI,EPXM) c t22 = P0 (E-XP0)-1 X PI (E-XP0)-1 : call mm(ndim,t22,t2,t1) if(iwr.eq.1)then call tk1(n,i,ndim,a,t22) call wmtr(a,24,24,'PI - part2 - new',1,6) endif do 4 ii=1,ndim do 4 jj=1,ndim 4 Pt(iw,ii,jj)=Pt(iw,ii,jj)+t22(ii,jj) write(6,*)' (E-XP0)-1 PI (E + X.(E-XP0)^(-1) . P0) made' return end c ============================================================ subroutine MDP(n,ndim,x,Ptr,Pt,inversion,NW,WMIN,DW,ntw,ntmax, 1GAMMA,tw,lgauss) implicit none integer*4 ndim,i,j,ii,it,jj,jt,k,kk,n,NW,ntw(*),iu,iw,ju,ntmax real*8 x(n*24,n*24),Ptr(ntmax,ndim,ndim),w,sf2,GAMMA, 1Pt(NW,ndim,ndim),WMIN,DW,sum,tw(n,ntmax),wi,wj real*8,allocatable::tem(:,:,:),temi(:,:,:),t2(:,:,:),a(:,:), 1ai(:,:) logical inversion,lgauss allocate(tem(ntmax,ndim,ndim),temi(ntmax,ndim,ndim)) c make X.P, associated transition frequency list is with j-index: call vz(tem,ntmax*ndim*ndim) do 1001 i=1,n do 1001 ii=24*(i-1)+1,24*i do 1001 j=1,n do 1001 jj=1,24*(j-1)+1,24*j do 1001 it=1,ntw(j) sum=0.0d0 do 10011 k=1,ndim 10011 sum=sum+x(ii,k)*Ptr(it,k,jj) 1001 tem(it,i,j)=sum write(6,*)' X.P made' c c make (E-X.P)^(-1) = temi if(inversion)then allocate(a(ndim,ndim),ai(ndim,ndim)) do 883 it=1,ntmax do 884 i=1,ndim do 884 j=1,ndim 884 a(i,j)=-tem(it,i,j) do 885 i=1,ndim 885 a(i,i)=1.0d0+a(i,i) call INV(a,ai,ndim,1.0d-25) do 883 i=1,ndim do 883 j=1,ndim 883 temi(it,i,j)=ai(i,j) deallocate(a,ai) else c perturbation approach c (E-X.P)^-1 ~ (E+X.P) do 882 i=1,ndim do 882 j=1,ndim do 882 it=1,ntmax 882 temi(it,i,j)=tem(it,i,j) do 992 i=1,n do 992 ii=24*(i-1)+1,24*i 992 temi(ntw(i),ii,ii)=1.0d0+temi(ntw(i),ii,ii) endif write(6,*)' (E-X.P)^(-1) made' deallocate(tem) c make total polarizability Pt = P(E-X.P)^(-1) for W <- 0: c = Ptr_w w' temi_w'_0 allocate(t2(n*ntmax,ndim,ndim)) do 74 iw=1,NW w=WMIN+dble(iw-1)*DW write(6,*)iw,NW c pre-transform temi to t2: do 75 j=1,n do 75 jj=24*(j-1)+1,24*(j-1)+24 do 75 k=1,n do 75 kk=24*(k-1)+1,24*(k-1)+24 do 75 i=1,n do 75 it=1,ntw(i) wi=tw(i,it) iu=ntmax*(i-1)+it sum=0.0d0 do 742 jt=1,ntw(j) ju=ntmax*(j-1)+jt wj=tw(j,jt) 742 sum=sum+temi(ju,kk,jj)*sf2(GAMMA,lgauss,w,0.0d0,wi,wj) 75 t2(iu,jj,kk)=sum do 74 i=1,n do 74 ii=24*(i-1)+1,24*(i-1)+24 do 74 j=1,n do 74 jj=24*(j-1)+1,24*(j-1)+24 sum=0.0d0 do 741 it=1,ntw(i) wi=tw(i,it) iu=ntmax*(i-1)+it do 741 k=1,n do 741 kk=24*(k-1)+1,24*(k-1)+24 741 sum=sum+Ptr(it,ii,kk)*t2(iu,kk,jj) 74 Pt(iw,ii,jj)=sum write(6,*)' P(E-X.P)^(-1) made' 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 sf2(d,g,a1,a2,wi,wj) implicit none logical g real*8 pi,spi,sf2,d,a1,a2,wi,wj,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((a1-a2-wi-wj)/d)**2 if(g)then if(dd.lt.20.0d0)then sf2=exp(-dd)/d/spi/dsqrt(2.0d0) else sf2=0.0d0 endif else if(dd.lt.1.0d10)then sf2=2.0d0/d/(dd+4.0d0)/pi else sf2=0.0d0 endif endif 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 vz(v,n) implicit none integer*4 n,i real*8 v(*) do 1 i=1,n 1 v(i)=0.0d0 return end c ============================================================ subroutine extd(alphar,gpr,ar,NW,c,no) c a(wj)=c(wj,wi)a(wi) implicit none integer*4 i,j,k,NW,iw,jw,no real*8 alphar(NW,3,3),gpr(NW,3,3),cji, 2ar(NW,3,3,3),c(no,no) real*8,allocatable::alphaf(:,:,:),gpf(:,:,:),af(:,:,:,:) if(no.ne.NW)call report('no <> NW') allocate(alphaf(NW,3,3),gpf(NW,3,3),af(NW,3,3,3)) call vz(alphaf,NW*9) call vz(gpf,NW*9) call vz(af,NW*27) do 1 jw=1,NW do 1 iw=1,NW cji=c(iw,jw) do 1 i=1,3 do 1 j=1,3 alphaf(jw,i,j)=alphaf(jw,i,j)+cji*alphar(iw,i,j) gpf( jw,i,j)=gpf( jw,i,j)+cji*gpr( iw,i,j) do 1 k=1,3 1 af( jw,i,j,k)=af( jw,i,j,k)+cji*ar( iw,i,j,k) alphar=alphaf gpr=gpf ar=af return end c ============================================================ subroutine extr(Pt,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(NW,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 for each frequency iw: do 8 iw=1,NW do 900 j=1,3 do 900 k=1,3 alphar(iw,j,k)=0.0d0 gpr(iw,j,k)=0.0d0 do 900 l=1,3 900 ar(iw,j,k,l)=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(iw,j,k)=alphar(iw,j,k)+Pt(iw,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(iw,ii+3+m,pp+3+j)/2.0d0 gpr(iw,j,k)=gpr(iw,j,k) -Pt(iw,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(iw,ii+kl+12,pp+j) 1+1.5d0 *(rr(l)*Pt(iw,ii+k ,pp+j) 1 +rr(k)*Pt(iw,ii+l ,pp+j)) 1-delta(k,l)*(rr(1)*Pt(iw,ii+1 ,pp+j) 1+ rr(2)*Pt(iw,ii+2 ,pp+j) 1+ rr(3)*Pt(iw,ii+3 ,pp+j)) 8 ar(iw,j,l,k)=ar(iw,j,k,l) write(6,*)' ROA tensors exctracted' 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 vp(a,b,c) real*8 a(*),b(*),c(*) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) 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 wrplane(ndim,P,x0,u,v,n,nx,ny,dx,dy,r,llog) implicit none integer*4 ndim,i,n,jx,kx,ia,ib,nx,ny,ix,istart,nstart,iz,nat real*8 P(ndim,ndim),x0(*),u(*),v(*),a(3),b(3),F(3,3),xa(3),sp, 1da,db,dx,dy,r(n,3),bohr,z,xi(3),c(3) logical llog character*80 fn real*8,allocatable::rr(:,:),XR(:,:) bohr=0.529177d0 do 1 i=1,3 a(i)=v(i)-x0(i) 1 b(i)=u(i)-x0(i) call nv(a) call vp(a,b,c) call nv(c) call vp(c,a,b) 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) open(9,file='PLANE.TXT') da=-dx do 2 ia=1,nx write(6,*)ia,' of ',nx da=da+dx db=-dy do 2 ib=1,ny db=db+dy c take each plane point as n+1 particle: do 3 ix=1,3 3 rr(n+1,ix)=x0(ix)/bohr+da*a(ix)+db*b(ix) call doX(n+1,rr,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 do 5 i=1,n istart=24*(i-1)+1 c F = X . P . F(0) do 5 jx=1,24 5 F(kx,ix)=F(kx,ix) 1+XR(nstart-1+ix,istart-1+jx)*P(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*bohr,db*bohr,z 601 format(2F12.2,G12.4) close(9) 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) 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 rdf(fn,n,ic,e,u) implicit none c read F.INP dipoles c ic=0 ... only find out the number c 1 ... read in character*(*) fn integer n,nat,ic,ia,i,ix real*8 e(*),u(*),cm cm=219474.0d0 open(30,file=fn) read(30,*)n,nat,nat if(ic.eq.1)then do 2 ia=1,2+nat*n+nat 2 read(30,*) read(30,300)(e(i),i=1,n) 300 format(6f11.3) do 3 i=1,n 3 read(30,*)(u(3*(i-1)+ix),ix=1,3) do 5 i=1,n 5 e(i)=e(i)/cm endif close(30) return end c ============================================================ subroutine wrdip(n,e,u) implicit none integer*4 n,i,ix real*8 e(*),u(*),cm cm=219474.0d0 write(6,*)' energy / cm-1 dipole x y z / au' do 1 i=1,n 1 write(6,600)i,e(i)*cm,(u(3*(i-1)+ix),ix=1,3) 600 format(i5,f12.4,3f12.6) return end c ============================================================ subroutine setv(n,v,e,r,u) implicit none integer*4 n,I,J real*8 e(*),u(*),v(n,n),r(n,3),xi,yi,zi,xj,yj,zj, 1uix,uiy,uiz,ujx,ujy,ujz,D2,D5,rij,uij,uir,ujr,xji,yji,zji c setup TDC potential matrix call vz(v,n*n) DO 1 I = 1,n 1 V(I,I)=E(I) c DO 210 I = 1,N xi=R(I,1) yi=R(I,2) zi=R(I,3) uix=u(3*(I-1)+1) uiy=u(3*(I-1)+2) uiz=u(3*(I-1)+3) DO 210 J = I+1,N xj=R(J,1) yj=R(J,2) zj=R(J,3) ujx=u(3*(J-1)+1) ujy=u(3*(J-1)+2) ujz=u(3*(J-1)+3) xji=xj-xi yji=yj-yi zji=zj-zi D2=(xi-xj)**2+(yi-yj)**2+(zi-zj)**2 rij=dsqrt(D2) D5=D2*D2*rij IF(rij.GT.1.0d-4)THEN uij=uix*ujx+uiy*ujy+uiz*ujz uir=uix*xji+uiy*yji+uiz*zji ujr=ujx*xji+ujy*yji+ujz*zji V(I,J)=(uij*D2-3.0d0*uir*ujr)/D5 V(J,I)=V(I,J) ENDIF 210 continue return end c ============================================================ SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(N) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(N),E(N) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END c ============================================================ subroutine INVG(G,EXP0,XP0,M,n,tol) implicit none integer*4 n,i,j,ii real*8 G,tol,XP0(n,n),M(n,n),EXP0(n,n) real*8,allocatable::t1(:,:),t2(:,:),t3(:,:),t4(:,:) allocate(t1(n,n),t2(n,n),t3(n,n),t4(n,n)) c make (E-XP0).(E-XP0)t: do 1 i=1,n do 1 j=1,n t1(i,j)=0.0d0 do 1 ii=1,n 1 t1(i,j)=t1(i,j)+EXP0(i,ii)*EXP0(j,ii) c make G**2.XP0.XP0t: do 2 i=1,n do 2 j=1,n t2(i,j)=0.0d0 do 2 ii=1,n 2 t2(i,j)=t2(i,j)+G**2*XP0(i,ii)*XP0(j,ii) c (E-XP0).(E-XP0)t + G.XP0.XP0t: do 3 i=1,n do 3 j=1,n 3 t3(i,j)=t1(i,j)+t2(i,j) c [(E-XP0).(E-XP0)t + G.XP0.XP0t]-1: call INV(t3,t4,n,tol) c [(E-XP0).(E-XP0)t + G.XP0.XP0t]-1.(E-XP0)t: do 4 i=1,n do 4 j=1,n M(i,j)=0.0d0 do 4 ii=1,n 4 M(i,j)=M(i,j)+EXP0(ii,i)*t4(ii,j) return end c ============================================================ subroutine readopt(ic,n,nm,m,ntw,ntmax,NW,nx,ny,ncub,icub,iwr, 1nxc,nyc,nzc,ncubc,icubc,ndim, ltdc,llog,ltab,lff,ldvm,lgauss, 1luseg,lusea,lcont,lplane,lcube,ltr,lave,inversion,debug, 1lonly,r,alpha,gp,a,alphad,gpd,ad,tw,WMAX,WMIN,W0,GAMMA,WA, 1temp,tol,tolw,x0,v,u,dx,dy,gw,x0c,dxc,dyc,dzc,DW,tempcm,cm, 1fn,ftw,ftwas,fnf,iiz,lasym,alphaas,alphadas) implicit none integer*4 ic,n,m,nm,ntw(*),ntmax,NW,nx,ny,ncub,icub,iwr, 1nxc,nyc,nzc,ncubc,icubc,ndim,i,iiz(*) logical ltdc,llog,ltab,lff,ldvm,lgauss,luseg,lusea,lcont, 1lplane,lcube,ltr(*),lave,inversion,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), 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(*),v(*),u(*),dx,dy,gw,x0c(*), 1dxc,dyc,dzc,DW,tempcm,cm,alphadas(n,nm,3,3) character*80 fn(*),ftw(*),fnf(*),ft,s80,ftwas(*) cm=219474.0d0 c lasym - use asymmetric polarizability c lave average by roattion c lcube - make CUB.CUB with intensity c lplane - make plane file PLANE.TXT with intensity 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 ltdc . transition dipole coupling model c luseg use Gp c lusea use A c inversion use true inversion c debug debug output lasym=.false. lave=.false. lcube=.false. lplane=.false. ldvm=.false. lcont=.false. lff=.false. lgauss=.false. llog=.false. lonly=.false. ltab=.false. ltdc=.false. luseg=.true. lusea=.true. inversion=.false. 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 icub=0 tol=0.01d0 tolw=0.2d0 WMAX=3500.0d0 WMIN=0.0d0 W0=0.0d0 GAMMA=10.0d0 c excitation freqency in A: wa=5320.0d0 TEMP=273.0d0 c grid spacing: dx=0.1d0 dy=0.1d0 c martix inversion damping factor: gw=0.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: call vz(x0,3) c plane points call vz(v,3) call vz(u,3) 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.'INVE')read(9,*)inversion 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.'LTDC')then read(9,*)ltdc do 61 i=1,n 61 read(9,900)fnf(i) 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,*)lplane read(9,*)(x0(i),i=1,3) read(9,*)(v( i),i=1,3) read(9,*)(u( i),i=1,3) read(9,*)nx,ny read(9,*)dx,dy read(9,*)ncub,icub 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,icubc 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=0.6950d0*TEMP 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,inversion,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,inversion,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,inversion 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,/, 2 ' True inversion :',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) call vz(a,n*n) 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 extrc(Pt,alphar,gpr,ar,NW,ndim,r,n,ntw) implicit none integer*4 i,j,k,l,ndim,NW,iw,n,kl,m,ii,jj,ntw(*),it real*8 Pt(NW,ndim,ndim),alphar(NW,3,3),gpr(NW,3,3), 2ar(NW,3,3,3),rr(3),r(n,3),eps,delta,sum,PP(24,24) c call vz(alphar,NW*3*3) call vz(gpr,NW*3*3) call vz(ar,NW*3*3*3) c for each frequency iw: iw=0 c sum over all multipoles on particles i: do 8 i=1,n do 8 it=1,ntw(i)-1 iw=iw+1 c contract last index over particles: call vz(PP,24*24) do 1 ii=1,24 do 1 jj=1,24 do 1 j=1,n 1 PP(ii,jj)=PP(ii,jj)+Pt(iw,24*(i-1)+ii,24*(j-1)+jj) rr(1)=r(i,1) rr(2)=r(i,2) rr(3)=r(i,3) do 8 j=1,3 kl=0 do 8 k=1,3 c from electric moment j = 1..3: c from electric field k = 1..3 alphar(iw,j,k)=alphar(iw,j,k)+PP(j,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)*PP(3+m,3+j)/2.0d0 gpr(iw,j,k)=gpr(iw,j,k) -PP(6+k,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*PP(kl+12,j) 1+1.5d0 *(rr(l)*PP(k ,j) 1 +rr(k)*PP(l ,j)) 1-delta(k,l)*(rr(1)*PP(1 ,j) 1+ rr(2)*PP(2 ,j) 1+ rr(3)*PP(3 ,j)) 8 ar(iw,j,l,k)=ar(iw,j,k,l) write(6,*)' ROA tensors exctracted' return end c ============================================================ subroutine smooth(v,n,t) implicit none real*8 v(*),t integer*4 i,n do 1 i=1,n 1 if(dabs(v(i)).lt.t)v(i)=0.0d0 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 tk1(n,i,ndim,a,b) implicit none integer*4 n,i,ndim,ii,jj,j real*8 a(24,24),b(ndim,ndim) do 108 ii=1,24 do 108 jj=1,24 a(ii,jj)=0.0d0 do 108 j=1,n 108 a(ii,jj)=a(ii,jj)+b(24*(i-1)+ii,24*(j-1)+jj) 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) call vz(PI,ndim*ndim) 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 ============================================================