PROGRAM NEW6U c experimental - asymmetric polarizabilities from transition dipole c derivatives IMPLICIT none INTEGER*4 ns0,N,NAT,NQ,NROOT,i,jj,ni,nl,iw,ip(10),IERR,NA, 1jjs,iwr,mmaxi,LEXCI,iwf,nb0,i2,i8,iex,ic,nproc, 1si,ii,mfc,iii,nc,mmi,nsi,tp,block,nga,nb,i3,i4,ib parameter (ns0=19,nb=10000) real*8 u0(3),m0(3),v0(3),q0(6),wz,wexc,we1,Gau,wmin,wmax,pis,dw, 1CM,dtol,kelvin,ee,ktlim,kt,qred,qp,wrax,enq,ei,bo, 1pf(nb),uu1(81) real*8,allocatable::ALPHA(:,:,:),A(:,:,:,:),G(:,:,:), 1dd(:),vv(:),aa(:),qq(:),ddi(:),vvi(:),rri(:),qqi(:), 1Gi(:,:,:),ALPHAi(:,:,:),Ai(:,:,:,:),gradi(:),we(:), 1GGi(:,:,:),AAi(:,:,:,:),GGr(:,:,:),AAr(:,:,:,:),grad(:), 1JT(:,:),J(:,:),K(:),Am(:,:),B(:),C(:,:),D(:),Em(:,:),W(:),WP(:), 1sg(:,:),se(:,:),jtj(:,:),jtji(:,:),X(:,:,:),Y(:,:,:),U(:,:,:), 1V(:,:,:),ki(:,:),sj(:,:),uu(:,:),ff(:,:,:),gg(:,:,:),Gw(:,:), 1X1(:,:),Y1(:,:),U1(:,:),V1(:,:),K1(:),f0(:),g0(:),f1(:),g1(:), 1roa(:,:,:),droa(:,:,:),ff1(:,:),gg1(:,:) integer*4,allocatable::ncen(:),nexc(:),wmaxi(:),sii(:), 1ntc(:),nte(:),ncen_(:,:),nexc_(:,:),nc_(:) LOGICAL luseA,luseG,ldo(ns0),lstat,slist,lvert,eattt, 1lt,limag character*80,allocatable::sx(:) C WRITE(6,6666) 6666 FORMAT(/,/, A' * * * * * * * * * * * * * * * * * * * * * *',/, 1' * * ',/, 2' * RESONANCE RAMAN OPTICAL ACTIVITY * ',/, 3' * * ',/, 4' * * * * * * * * * * * * * * * * * * * * * *',/, 8' ',/, A'INPUT FILES: NEW6U.OPT, FILE.TTT ',/, A' TD.OUT and TD.OUT.DUSCH, or list of files',/, B'OUTPUT : .TAB files with intensities',/) CALL IOPAR(wexc,wmin,wmax,luseA,luseG,ldo,Gau,pis, 1lstat,slist,iwr,lvert,eattt,dtol,kelvin, 1mmaxi,LEXCI,ktlim,wrax,NA,block,dw,kt,nproc,limag) CM=219474.63d0 if(slist)then ni=nl('SLIST.TXT') allocate(sx(ni)) call fl(ni,'SLIST.TXT',sx) else ni=1 allocate(sx(ni)) sx(1)='TD.OUT' endif jj=jjs(sx(1)) NAT=nga(sx(1)(1:jj)) N=3*NAT write(6,6090)NAT 6090 format(i6,' atoms') allocate(ALPHA(N,3,3),ALPHAi(N,3,3),G(N,3,3),Gi(N,3,3), 1A(N,3,3,3),Ai(N,3,3,3),GGr(N,3,3),GGi(N,3,3),AAr(N,3,3,3), 1AAi(N,3,3,3),roa(2,ns0,NA),droa(2,ns0,NA), 1grad(N),gradi(N),sg(N,N),Gw(N,N),se(N,N),Am(N,N),B(N), 1C(N,N),D(N),Em(N,N),W(N),WP(N),dd(3*N),aa(3*N),vv(3*N),qq(6*N), 1ddi(3*N),rri(3*N),vvi(3*N),qqi(6*N),we(ni),f0(ni),g0(ni),f1(ni), 1g1(ni)) roa=0.0d0 c ground state S-matrix and frequencies: call t0(N,ALPHA,ALPHAi,G,Gi,GGr,GGi,A,Ai,AAr,AAi) Gw=0.0d0 call rduschs(N,NQ,sg,sx(1)(1:jj)//'.DUSCH','Ground S-Matrix') call rduschw(N,W,sx(1)(1:jj)//'.DUSCH','Ground w') write(6,6091)NQ 6091 format(i6,' normal modes') do 1 ii=1,NQ 1 Gw(ii,ii)=W(ii) if(iwr.ne.0)call control(N,NQ,sg) allocate(J(NQ,NQ),JT(NQ,NQ),K(NQ),jtj(NQ,NQ),jtji(NQ,NQ), 1U1(NQ,NQ),V1(NQ,NQ),X1(NQ,NQ),Y1(NQ,NQ),K1(NQ),sj(NQ,NQ), 1U(ni,NQ,NQ),V(ni,NQ,NQ),X(ni,NQ,NQ),Y(ni,NQ,NQ),ki(ni,NQ), 1uu(ni,81),ff(ni,81,NQ),gg(ni,81,NQ),ff1(81,NQ),gg1(81,NQ)) c static polarizability if(lstat)call rdstat(NQ,N,ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi, 1sg,W,wmin,wmax,limag) c cycle over files/electronic excited states: c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee do 1000 i=1,ni write(6,609)i 609 format(i6,'. electronic state:') jj=jjs(sx(i)) c read S and the Dushchinsky matrices: call rduschs(N,NQ,se,sx(i)(1:jj)//'.DUSCH','Excited S-Matri') call rdusch(NQ,N,JT,K,Am,B,C,D,Em,W,WP,iw,ip, 1sx(i)(1:jj)//'.DUSCH') if(ip(1).eq.0)call report('Duschinsky Matrix not found') if(ip(2).eq.0)call report('Shift Vector not found') if(ip(3).eq.0)call report('A Matrix not found') if(ip(4).eq.0)call report('B Vector not found') if(ip(5).eq.0)call report('C Matrix not found') if(ip(6).eq.0)call report('D Vector not found') if(ip(7).eq.0)call report('E Matrix not found') if(iw.lt.2)call report('Frequencies not found') call chkj(NQ,JT,block,W,WP) call mt(J,JT,NQ) c jtj = Jt . J call mm(jtj,JT,J,NQ) c jtji = (Jt . J)^-1 call INV(jtj,jtji,NQ,IERR) if(IERR.ne.0)call report('Inversion failed') c sj = j = (Jt . J)^-1 . Jt call mm(sj,jtji,JT,NQ) c find which excited state is calculated call chkroot(sx(i),NROOT) c read moments: call r00(sx(i),NROOT,u0,m0,v0,q0,we1) c frequency multipliers: call mkf0(i,ni,we1,we,wexc,Gau,f0,g0,f1,g1) c read moment derivatives: call rdd(sx(i),nat,NROOT,u0,v0,m0,q0,dd,vv,aa,qq,grad) call wrge('DEGS.TEN',nat,dd,aa,qq,vv,u0,m0,v0,q0) write(6,*)'Analytical transition derivatives, current norms:' call wnn(9*nat,dd,'NS_d:') call wnn(9*nat,vv,'NS_v:') call wnn(9*nat,aa,'NS_m:') call wnn(18*nat,qq,'NS_q:') c erase some atom derivatives if required: if(eattt)call ede(nat,dd,vv,aa,qq,grad) C transform into ground state normal modes call trafN(3,dd,ddi,nat,NQ,sg) call trafN(3,aa,rri,nat,NQ,sg) call trafN(6,qq,qqi,nat,NQ,sg) call trafN(3,vv,vvi,nat,NQ,sg) call trafN(1,grad,gradi,nat,NQ,sg) c store transition dipoles and derivatives in uu1,ff1,gg1: call sc(N,NQ,u0,m0,q0,ddi,rri,qqi,uu1,ff1,gg1) c store in larger arrays, divide ff1 by dsqrt(2w): call sl(i,ni,NQ,N,uu1,ff1,gg1,uu,ff,gg,W) c vertical approximation, make new shift vectors: if(lvert)call mksv(NQ,N,K,gradi,WP,J,W,iwr,block) c make matrices call m3(N,NQ,K1,U1,V1,X1,Y1,W,WP,sj,wz,block) c store in larger arrays: 1000 call s3(ni,i,NQ,K1,U1,V1,X1,Y1,ki,U,V,X,Y) c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee c partition function exact, for information: qp=1.0d0 do 2 ii=block+1,NQ 2 qp=qp/(1.0d0-exp(-W(ii)*CM/kt)) write(6,648)qp 648 format(' Partition function, exact:',e12.4) c states mmi excited allocate(sii(LEXCI)) ei=0.0d0 bo=exp(-ei*CM/kt) qred=bo nb0=1 do 3 mmi=1,LEXCI c distribute mmi excitations upon centers b+1..NQ: c initial indices - all to the first center: ei=0.0d0 do 41 iex=1,mmi ei=ei+W(block+1) 41 sii(iex)=block+1 5 if(ei*CM.lt.20.0d0*kt)then qred=qred+exp(-ei*CM/kt) nb0=nb0+1 endif c find index to be changed do 8 ic=mmi,1,-1 8 if(sii(ic).lt.NQ)goto 9 goto 3 9 do 10 jj=ic+1,mmi ei=ei-W(sii(jj))+W(sii(ic)+1) 10 sii(jj)=sii(ic)+1 ei=ei-W(sii(ic))+W(sii(ic)+1) sii(ic)=sii(ic)+1 goto 5 3 continue write(6,650)' LEXCI = ',LEXCI,qred,nb0 call omp_set_num_threads(nproc) iwf=LEXCI+1 nb0=0 qred=0.0d0 do 55 i2=1,2 c i2: count states only, i2: make integrals tp=0 si=0 i8=0 do 2011 mmi=0,mmaxi ii=max(1,mmi) allocate(ncen(ii),nexc(ii),wmaxi(ii),ntc(ii),nte(ii), 1nc_(nb),ncen_(ii,nb),nexc_(ii,nb)) ib=0 wmaxi=LEXCI+1 nsi=iwf**mmi do 201 iii=1,mfc(NQ-block,mmi) call getdis(iii,ncen,NQ-block,mmi,ii) Do 201 jj = 1, nsi Call FCA2NQ(mmi,jj,wmaxi,lt,nc,nexc) if(lt)then ei=enq(nc,N,ncen,Gw,nexc,block) if(ei*CM.lt.10.0d0*kt)then bo=exp(-ei*CM/kt) if(bo.gt.ktlim)then si=si+1 if(i2.eq.2)then if((100*si)/nb0.gt.tp)then write(6,646)tp 646 format(i4,'%',$) tp=(100*si)/nb0 i8=i8+1 if(mod(i8,12).eq.0)write(6,*) endif ib=ib+1 nc_(ib)=mmi do 301 i3=1,mmi ncen_(i3,ib)=ncen(i3) 301 nexc_(i3,ib)=nexc(i3) pf(ib)=bo/qred if(ib.eq.nb)then C$OMP Parallel do Default(Shared) C$OMP+ Private(i3,i4,ntc,nte,droa) do 299 i3=1,ib do 3011 i4=1,nc_(i3) ntc(i4)=ncen_(i4,i3) 3011 nte(i4)=nexc_(i4,i3) droa=0.0d0 call mkijf(ni,N,NA,NQ,W,nc_(i3),ntc,nte,uu,ff,gg, 1 X,Y,U,V,ki,wmin,dw,pf(i3),droa,pis,lstat, 1 ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi, l luseA,luseG,ldo,we,wexc,f0,g0,f1,g1,block,iwr) C$OMP Critical roa=roa+droa C$OMP End Critical 299 continue ib=0 endif else qred=qred+bo if(iwr.gt.0)call wssr(' From: ',si,ei*CM,nc,ncen,nexc) endif endif endif endif 201 continue if(ib.gt.0)then C$OMP Parallel do Default(Shared) C$OMP+ Private(i3,i4,ntc,nte,droa) do 601 i3=1,ib do 6011 i4=1,nc_(i3) ntc(i4)=ncen_(i4,i3) 6011 nte(i4)=nexc_(i4,i3) droa=0.0d0 call mkijf(ni,N,NA,NQ,W,nc_(i3),ntc,nte,uu,ff,gg, 1 X,Y,U,V,ki,wmin,dw,pf(i3),droa,pis,lstat, 1 ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi, l luseA,luseG,ldo,we,wexc,f0,g0,f1,g1,block,iwr) C$OMP Critical roa=roa+droa C$OMP End Critical 601 continue ib=0 endif 2011 deallocate(ncen,wmaxi,nexc,ntc,nte,nc_,ncen_,nexc_) nb0=si if(i2.eq.1)write(6,650)' mmaxi = ',mmaxi,qred,si 650 format(' Partition function for',a9,I6,':',e12.4,',',i9,' states') 55 continue write(6,*) c call initab(ldo) ee=wmin-dw do 56 i=1,NA ee=ee+dw do 56 jj=1,ns0 if(ldo(jj))then if(dabs(roa(1,jj,i)).gt.1.0d-12)WRITE(50+jj,3001)i,ee*CM, 1 0.0d0,0.0d0,roa(1,jj,i),roa(2,jj,i) 3001 FORMAT(I7,f9.2,2f3.0,g12.4,' 0 0 0 0',2g12.4) endif 56 continue call closetab(ldo) c END c ============================================================== function nl(f) c number of lines in a file implicit none character*(*) f integer*4 nl,n open(9,file=f) n=0 1 read(9,*,end=99,err=99) n=n+1 goto 1 99 close(9) nl=n return end c ============================================================== subroutine rdd(f,nat,NROOT,u0,v0,m0,q0,d,v,a,q,grad) c read transition dipole derivatives from excited state freq calc implicit none integer*4 nat,ia,xa,i,NROOT,ii,ix,nd,ntr,j,ir real*8 d(9*nat),u0(*),m0(3),a(9*nat), 1v(9*nat),v0(3),q(18*nat),q0(6),ev,eau,grad(3*nat) character*(*) f character*80 s80 real*8,allocatable::v16(:,:) c write(6,*) write(6,*)' Rdd: reading '//f write(6,*) c total number of roots (exc. states) calculated: ntr=0 open(9,file=f,status='old') 1 read(9,900,end=88,err=88)s80 900 format(a80) if(s80(38:59).eq.'Forces (Hartrees/Bohr)')then read(9,*) read(9,*) do 2 i=1,nat 2 read(9,*)(grad(ix+3*(i-1)),ix=1,2),(grad(ix+3*(i-1)),ix=1,3) endif if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then do 20 i=1,79 if(s80(i:i) .eq.':') read(s80(15:i-1),*)nd 20 if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev if(nd.eq.NROOT)then eau=ev/27.211384205943d0 endif endif if(s80(2:31).eq.'Electronic transition elements')then c read NROOT moments allocate(v16(16,NROOT)) call ru(9,v16,16,NROOT) do 21 ix=1,3 u0(ix)=v16(ix+ 1,NROOT) v0(ix)=v16(ix+ 4,NROOT) 21 m0(ix)=v16(ix+ 7,NROOT) q0(1 )=v16(11,NROOT) q0(2 )=v16(14,NROOT) q0(3 )=v16(15,NROOT) q0(4 )=v16(12,NROOT) q0(5 )=v16(16,NROOT) q0(6 )=v16(13,NROOT) deallocate(v16) write(6,*)'Analytical transitions read' endif if(s80(2:34).eq.'Electronic Transition Derivatives')then allocate(v16(16,3 + 3*nat)) call ru(9,v16,16,3 + 3*nat) do 22 ia=1,nat do 22 xa=1,3 ii=xa+3*(ia-1) do 221 ix=1,3 d(ix +3*(ii-1))=v16(ix+ 1,3+ii) v(ix +3*(ii-1))=v16(ix+ 4,3+ii) 221 a(ix +3*(ii-1))=v16(ix+ 7,3+ii) q(1+6*(ii-1))=v16(11,3+ii) q(2+6*(ii-1))=v16(14,3+ii) q(3+6*(ii-1))=v16(15,3+ii) q(4+6*(ii-1))=v16(12,3+ii) q(5+6*(ii-1))=v16(16,3+ii) 22 q(6+6*(ii-1))=v16(13,3+ii) deallocate(v16) write(6,*)'Analytical transition derivatives read, norms:' call wnn(9*nat,d,'N_d:') call wnn(9*nat,v,'N_v:') call wnn(9*nat,a,'N_m:') call wnn(18*nat,q,'N_q:') endif c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA c checkpoint style c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC if(s80(1:18).eq.'ETran state values')then if(ntr.eq.0)call report('nstates not found') allocate(v16(1,16*ntr+48+16*3*nat)) read(9,*)v16 ii=0 do 237 ir=1,ntr ii=ii+1 do 231 ix=1,3 ii=ii+1 231 if(ir.eq.NROOT)u0(ix)=v16(1,ii) do 232 ix=1,3 ii=ii+1 232 if(ir.eq.NROOT)v0(ix)=v16(1,ii) do 233 ix=1,3 ii=ii+1 233 if(ir.eq.NROOT)m0(ix)=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(1 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(4 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(6 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(2 )=v16(1,ii) ii=ii+1 if(ir.eq.NROOT)q0(3 )=v16(1,ii) ii=ii+1 237 if(ir.eq.NROOT)q0(5 )=v16(1,ii) ii=ii+48 do 238 j=1,3*nat ii=ii+1 do 234 ix=1,3 ii=ii+1 234 d(ix+3*(j-1))=v16(1,ii) do 235 ix=1,3 ii=ii+1 235 v(ix+3*(j-1))=v16(1,ii) do 236 ix=1,3 ii=ii+1 236 a(ix+3*(j-1))=v16(1,ii) ii=ii+1 q(1+6*(j-1))=v16(1,ii) ii=ii+1 q(4+6*(j-1))=v16(1,ii) ii=ii+1 q(6+6*(j-1))=v16(1,ii) ii=ii+1 q(2+6*(j-1))=v16(1,ii) ii=ii+1 q(3+6*(j-1))=v16(1,ii) ii=ii+1 238 q(5+6*(j-1))=v16(1,ii) deallocate(v16) write(6,*)'Analytical transition derivatives read from chk' endif c CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC goto 1 88 close(9) if(nd.eq.0)call report('Energy not found') v = -v/eau q = -q/eau a = a/2.0d0 v0 = -v0/eau q0 = -q0/eau m0 = m0/2.0d0 call wrge('DEG.TEN',nat,d,a,q,v,u0,m0,v0,q0) return end C ========================================== subroutine wrge(f,nat,d,a,q,v,u0,m0,v0,q0) implicit none integer*4 nat,ia,xa,ix real*8 d(9*nat),v(9*nat),a(9*nat),q(18*nat), 1u0(3),v0(3),m0(3),q0(6) character*(*) f open(40,file=f) write(40,*)'atom coord dipx dipy dipz' do 10 ia=1,nat do 10 xa=1,3 10 write(40,400)ia,xa,(d(ix+3*(xa-1)+9*(ia-1)),ix=1,3) write(40,*)'equilibrium transition electric dipole:' write(40,400)0,0,(u0(ix),ix=1,3) write(40,*)'atom coordmdipx mdipy mdipz' do 101 ia=1,nat do 101 xa=1,3 101 write(40,400)ia,xa,(a(ix+3*(xa-1)+9*(ia-1)),ix=1,3) write(40,*)'equilibrium transition magnetic dipole:' write(40,400)0,0,(m0(ix),ix=1,3) write(40,*)'atom coord Qv xx xy xz yy yz zz:' do 102 ia=1,nat do 102 xa=1,3 102 write(40,400)ia,xa,(q(ix+6*(xa-1)+18*(ia-1)),ix=1,6) write(40,*)'equilibrium transition quadrupole velocity:' write(40,400)0,0,(q0(ix),ix=1,6) write(40,*)'atom coord grad x y z:' do 103 ia=1,nat do 103 xa=1,3 103 write(40,400)ia,xa,(v(ix+3*(xa-1)+9*(ia-1)),ix=1,3) write(40,*)'equilibrium transition velocity:' write(40,400)0,0,(v0(ix),ix=1,3) 400 format(i5,i2,6g12.4) close(40) write(6,*)f//' written' return end C ========================================== subroutine wnn(N,v,f) integer*4 N,i real*8 v(N),a character*4 f a=0.0d0 do 1 i=1,N 1 a=a+v(i)**2 a=dsqrt(a) write(6,600)f,a 600 format(a4,e12.4) return end C ========================================== SUBROUTINE READTEN(N,P,G,A,ic) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N,3,3),G(N,3,3),A(N,3,3,3) OPEN(2,FILE='FILE.TTT',STATUS='OLD') READ(2,*) READ(2,*)NAT if(NAT.ne.N/3)call report('inconsistent number of atoms') READ(2,*) READ(2,*) DO 1 I=1,3 READ(2,*) DO 1 L=1,NAT DO 1 IX=1,3 M=3*(L-1)+IX 1 READ(2,*)(P(M,I,J),J=1,2),(P(M,I,J),J=1,ic),(P(M,I,J),J=1,3) READ(2,*) READ(2,*) DO 2 I=1,3 READ(2,*) DO 2 L=1,NAT DO 2 IX=1,3 M=3*(L-1)+IX c last index magnetic, inner 2 READ(2,*)(G(M,I,J),J=1,2),(G(M,I,J),J=1,ic),(G(M,I,J),J=1,3) READ(2,*) READ(2,*) DO 3 I=1,3 DO 3 J=1,3 READ(2,*) DO 3 L=1,NAT DO 3 IX=1,3 M=3*(L-1)+IX 3 READ(2,*)(A(M,I,J,K),K=1,2),(A(M,I,J,K),K=1,ic), 1(A(M,I,J,K),K=1,3) CLOSE(2) WRITE(6,*)' tensors read in from FILE.TTT' RETURN END C ============================================================== SUBROUTINE IOPAR(w,wmin,wmax,luseA,luseG,ldo,Gau, 1pis,lstat,slist,iwr,lvert,eattt,dtol,kelvin, 1mmaxi,LEXCI,ktlim,wrax,NA,block,dw,kt,nproc,limag) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) parameter (ns0=19) integer*4 block LOGICAL lex,luseA,luseG,ldo(ns0),lstat,slist,lvert,eattt,limag real*8 ktlim,kelvin,kt CHARACTER*15 STRING character*11 st(ns0) data st/'ICP_0 ','ICPx_90 ','ICPz_90 ','ICPs_90 ', 1 'ICPu_90 ','ICP_180 ','SCP_0 ','SCPx_90 ', 1 'SCPz_90 ','SCPs_90 ','SCPu_90 ','SCP_180 ', 1 'DCPI_0 ','DCPI_90 ','DCPI_180 ','DCPII_0 ', 1 'DCPII_90 ','DCPII_180 ','FCTAB '/ c imaginary static tensors: limag=.false. c number of processors nproc=1 c number of modes to block: block=6 c maximal initial and final excitation classess,excitations: mmaxi=3 LEXCI=3 c number of frequency points for alphas: NA=1000 c maximal Raman frequency cm-1: wrax=4000.0d0 c limit for exp(-E/kT): ktlim=0.20d0 c temperature kelvin=300.0d0 c tolerance for geometry: dtol=0.2d0 iwr=0 c output amount: iwr=0 c add complex static polarizability: lstat=.false. c SLIST.TXT of excited state outputs slist=.false. c vertical approximation: lvert=.true. c consider EA.LST: eattt=.false. pis=1.0d0 CM=219470.0d0 Gammacm=500.0d0 ldo=.false. c spectral limits, in cm-1: wmin=10.0d0 wmax=4000.0d0 c excitation light wavelength in nm: EXCNM=532.0d0 luseA=.true. c use G'-tensor derivatives: luseG=.true. INQUIRE(FILE='NEW6U.OPT',exist=lex) if(lex)then OPEN(2,FILE='NEW6U.OPT') WRITE(6,*)' FILE NEW6U.OPT FOUND ...' else goto 99 endif 2001 READ(2,1111)STRING 1111 FORMAT(A15) IF(STRING(1:4).EQ.'IMAG')READ(2,*)limag IF(STRING(1:4).EQ.'BLOC')READ(2,*)block IF(STRING(1:4).EQ.'NPRO')READ(2,*)nproc IF(STRING(1:4).EQ.'DTOL')READ(2,*)dtol IF(STRING(1:4).EQ.'NAAL')READ(2,*)NA IF(STRING(1:3).EQ.'IWR')READ(2,*)iwr IF(STRING(1:5).EQ.'EATTT')READ(2,*)eattt IF(STRING(1:4).EQ.'STAT')READ(2,*)lstat IF(STRING(1:4).EQ.'SLIS')READ(2,*)slist IF(STRING(1:4).EQ.'WMIN')READ(2,*)wmin IF(STRING(1:4).EQ.'WMAX')READ(2,*)wmax IF(STRING(1:5).EQ.'EXCNM')READ(2,*)EXCNM IF(STRING(1:4).EQ.'USEA')READ(2,*)luseA IF(STRING(1:4).EQ.'USEG')READ(2,*)luseG IF(STRING(1:4).EQ.'PISV')READ(2,*)pis IF(STRING(1:4).EQ.'VERT')READ(2,*)lvert IF(STRING(1:5).EQ.'GAMMA')READ(2,*)Gammacm IF(STRING(1:5).EQ.'KELVI')READ(2,*)kelvin IF(STRING(1:5).EQ.'MMAXI')READ(2,*)mmaxi IF(STRING(1:5).EQ.'LEXCI')READ(2,*)lexci IF(STRING(1:5).EQ.'KTLIM')READ(2,*)ktlim IF(STRING(1:4).EQ.'WRAX')READ(2,*)wrax do 4 ie=1,len(STRING) 4 if(STRING(ie:ie).eq.' ')goto 3 3 ie=ie-1 do 5 i=1,ns0 5 if(STRING(1:ie).eq.st(i)(1:ie))read(2,*)ldo(i) IF(STRING(1:3).EQ.'END')GOTO 2002 GOTO 2001 2002 CLOSE(2) 99 write(6,6009)EXCNM,Gammacm,wmin,wmax 6009 format(' l_exc ',f5.1,' nm, Gamma:',f10.1,' cm-1',/, 1' wmin:',f10.1,' cm-1, wmax:',f10.1,' cm-1 (converted to au now)') w=1.0d7/EXCNM/CM Gau=Gammacm/CM wmin=wmin/CM wmax=wmax/CM dw=(wmax-wmin)/dble(NA-1) write(6,599)NA,dw*cm,w*CM,Gau*CM 599 format(i6,' alpha points, dw:',f9.2,' cm-1',/, 1' wexc/cm-1:',f10.0,' G/cm-1:',f10.0,/) kt=kelvin*0.6950356d0 RETURN END SUBROUTINE TRTEN(S,N,ALPHA,A,G) IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) DIMENSION S(N,N),A(N,3,3,3),ALPHA(N,3,3),G(N,3,3) real*8,allocatable::V(:) allocate(V(N)) DO 1 I=1,3 DO 1 J=1,3 DO 3 JQ=1,N 3 V(JQ)=ALPHA(JQ,I,J) DO 1 IQ=1,N ALPHAS=0.0d0 DO 2 IX=1,N 2 ALPHAS=ALPHAS+S(IX,IQ)*V(IX) 1 ALPHA(IQ,I,J)=ALPHAS DO 11 I=1,3 DO 11 J=1,3 DO 31 JQ=1,N 31 V(JQ)=G(JQ,I,J) DO 11 IQ=1,N GS=0.0d0 DO 21 IX=1,N 21 GS=GS+S(IX,IQ)*V(IX) 11 G(IQ,I,J)=GS DO 12 I=1,3 DO 12 J=1,3 DO 12 K=1,3 DO 32 JQ=1,N 32 V(JQ)=A(JQ,I,J,K) DO 12 IQ=1,N AS=0.0D0 DO 22 IX=1,N 22 AS=AS+S(IX,IQ)*V(IX) 12 A(IQ,I,J,K)=AS RETURN END subroutine report(s) character*(*) s write(6,*)s stop end subroutine r00(f,NROOT,u0,m0,v0,q0,eau) c read transition dipoles from excited state freq calc implicit none integer*4 ln,i,NROOT,ic,ix,nd,i1,i2,i3,i4,i5 real*8 u0(*),m0(3),v0(3),q0(6),ev,eau character*(*) f character*80 s80 real*8,allocatable::v16(:,:) c write(6,*) write(6,*)' R00: reading '//f write(6,*) open(9,file=f,status='old') ln=0 i1=0 i2=0 i3=0 i4=0 i5=0 c where it starts: 1 read(9,900,end=88,err=88)s80 900 format(a80) ln=ln+1 if(s80(1:2).eq.' #')then ic=0 do 3 i=1,len(s80)-4 if(s80(i:i+4).eq.' freq'.or.s80(i:i+4).eq.' freq')ic=ic+1 3 if(s80(i:i+2).eq.' td' .or.s80(i:i+2).eq.' TD' )ic=ic+1 if(ic.eq.2)then write(6,*)'TD starts at line ',ln goto 4 endif endif goto 1 88 close(9) call report('TD freq not found') 4 read(9,900,end=78,err=78)s80 if(s80(2:65).eq.'Ground to excited state transition electric dipol 1e moments (Au):'.and.i1.eq.0)then write(6,900)s80 do 5 i=1,NROOT 5 read(9,*) read(9,*)u0(1),(u0(ix),ix=1,3) i1=i1+1 endif if(s80(2:65).eq.'Ground to excited state transition velocity dipol 1e moments (Au):'.and.i2.eq.0)then write(6,900)s80 do 11 i=1,NROOT 11 read(9,*) read(9,*)v0(1),(v0(ix),ix=1,3) i2=i2+1 endif if(s80(2:65).eq.'Ground to excited state transition magnetic dipol 1e moments (Au):'.and.i3.eq.0)then write(6,900)s80 do 14 i=1,NROOT 14 read(9,*) read(9,*)m0(1),(m0(ix),ix=1,3) i3=i3+1 endif if(s80(2:69).eq.'Ground to excited state transition velocity quadr 1upole moments (Au):'.and.i4.eq.0)then write(6,900)s80 do 17 i=1,NROOT 17 read(9,*) read(9,*)q0(1),q0(1),q0(4),q0(6),q0(2),q0(3),q0(5) i4=i4+1 endif if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s'. 1and.i5.eq.0)then do 20 i=1,79 if(s80(i:i) .eq.':') read(s80(15:i-1),*)nd 20 if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev if(nd.eq.NROOT)then write(6,900)s80 eau=ev/27.211384205943d0 i5=i5+1 endif endif c analytical frequency output: c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA if(s80(2:31).eq.'Electronic transition elements')then c read NROOT moments write(6,900)s80 allocate(v16(16,NROOT)) call ru(9,v16,16,NROOT) do 21 ix=1,3 u0(ix )=v16(ix+ 1,NROOT) v0(ix )=v16(ix+ 4,NROOT) 21 m0(ix )=v16(ix+ 7,NROOT) q0(1)=v16(11,NROOT) q0(2)=v16(14,NROOT) q0(3)=v16(15,NROOT) q0(4)=v16(12,NROOT) q0(5)=v16(16,NROOT) q0(6)=v16(13,NROOT) c our order 1 2 3 4 5 6 c our order xx xy xz yy yz zz c 1 2 3 4 5 6 c 11 12 13 14 15 16 c gaussian # xx yy zz xy xz yz deallocate(v16) write(6,*)'Analytical transitions read' endif c AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA goto 4 78 close(9) if(nd.eq.0)call report('Energy not found') v0 = -v0 / eau q0 = -q0 / eau m0 = m0 / 2.0d0 c looks like that we have this situation: c Gaussian gives <0|r|e>, , <0|-r x grad|e> c and c To transform it to what we want we can use c hbar^2 c <0|r|e> = ------- <0|grad|e> c m Ee0 c and c hbar^2 c <0|rr|e> = ------- <0|r grad + grad r|e> c m E0 c Eeo=Ee-E0 c so that (in atomic units hbar=m=1) c <0|r|e>(in velocity) = - / Ee0 c Im <0|m|e>= Im <0|-r x grad/2|e> = <0|-r x grad|e> / 2 c <0|rr|e> = - / Ee0 c return end subroutine setu(u,ddi,iq) implicit none integer*4 iq,a real*8 u(3),ddi(*) do 1 a=1,3 1 u(a)=ddi(a+3*(iq-1)) return end subroutine ads1(wexc,w,f,N,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1gcsr,gcsi,Gcr,Gci,asr,asi,A,Ai,acsr,acsi,Acr,Aci,pis,iwr) implicit none integer*4 N,i,j,f,l,iwr real*8 aasr(3,3),aasi(3,3),ALPHA(N,3,3),wexc,w,clight,h,o,pis, 1ALPHAi(N,3,3),gsr(3,3),gsi(3,3),gcsr(3,3),gcsi(3,3), 1G(N,3,3),Gi(N,3,3),asr(3,3,3),asi(3,3,3),acsr(3,3,3), 1Gcr(N,3,3),Gci(N,3,3),acr(N,3,3,3),aci(N,3,3,3), 1acsi(3,3,3),A(N,3,3,3),Ai(N,3,3,3) clight=137.03599d0 h=1.0d0/dsqrt(2.0d0*w) o=h*clight/wexc do 1 i=1,3 do 1 j=1,3 aasr(i,j)=aasr(i,j)+pis*h* ALPHA( f,i,j) aasi(i,j)=aasi(i,j)+pis*h* ALPHAi(f,i,j) gsr( i,j)=gsr( i,j)+pis*h*wexc*G( f,i,j) gsi( i,j)=gsi( i,j)+pis*h*wexc*Gi( f,i,j) gcsr(i,j)=gcsr(i,j)+pis*h*wexc*Gcr( f,i,j) gcsi(i,j)=gcsi(i,j)+pis*h*wexc*Gci( f,i,j) do 1 l=1,3 asr( i,j,l)=asr( i,j,l)+pis*o*A( f,i,j,l) asi( i,j,l)=asi( i,j,l)+pis*o*Ai( f,i,j,l) acsr(i,j,l)=acsr(i,j,l)+pis*o*Acr( f,i,j,l) 1 acsi(i,j,l)=acsi(i,j,l)+pis*o*Aci( f,i,j,l) if(iwr.gt.0)then write(6,600)f,aasr(1,1)+aasr(2,2)+aasr(3,3), 1 h*(ALPHA(f,1,1)+ALPHA(f,2,2)+ALPHA(f,3,3)),pis 600 format(' transition',i6,' Sp av / Sp as:',e11.3,' /', 1 e11.3,' Pis:',e10.2) endif return end c ============================================================ subroutine closetab(ldo) implicit none integer*4 i,ns0 parameter(ns0=19) logical ldo(ns0) do 1 i=1,ns0 if(ldo(i))then write(50+i,3000) 3000 format(60(1h-)) close(50+i) endif 1 continue return end c ============================================================ subroutine initab(ldo) implicit none integer*4 i,j,ns0 parameter(ns0=19) character*11 st(ns0) logical ldo(ns0) data st/'ICP_0 ','ICPx_90 ','ICPz_90 ','ICPs_90 ', 1 'ICPu_90 ','ICP_180 ','SCP_0 ','SCPx_90 ', 1 'SCPz_90 ','SCPs_90 ','SCPu_90 ','SCP_180 ', 1 'DCPI_0 ','DCPI_90 ','DCPI_180 ','DCPII_0 ', 1 'DCPII_90 ','DCPII_180 ','FCTAB '/ do 1 i=1,ns0 if(ldo(i))then write(6,600)st(i) 600 format(a11,$) do 2 j=len(st(i)),1,-1 2 if(st(i)(j:j).ne.' ')goto 3 3 open(50+i,file=st(i)(1:j)//'.tab') write(50+i,3000)st(i) 3000 format(20x,a11,/,' n wavelength (cm) Raman ROA',/,80(1h-)) endif 1 continue write(6,*) return end c ============================================================== subroutine trafN(N,dd,ddi,nat,NQ,s) implicit none integer*4 N,nat,NQ,ix,iq,ia,xa,iqi,ii real*8 dd(*),ddi(*),s(3*nat,3*nat) do 1 ix=1,N do 1 iq=1,NQ iqi=ix+N*(iq-1) ddi(iqi)=0.0d0 do 2 ia=1,nat do 2 xa=1,3 ii=xa+3*(ia-1) c s-opposite order of normal modes: 2 ddi(iqi)=ddi(iqi)+dd(ix+N*(ii-1))*s(ii,iq) 1 continue return end c ============================================================ subroutine calcep(e,a) c eap_ab = eps_abc A_d,cd implicit none real*8 e(3,3),a(3,3,3) e(1,1)=0.0d0 e(2,2)=0.0d0 e(3,3)=0.0d0 e(1,2)= a(1,3,1)+a(2,3,2)+a(3,3,3) e(1,3)=-a(1,2,1)-a(2,2,2)-a(3,2,3) e(2,3)= a(1,1,1)+a(2,1,2)+a(3,1,3) e(2,1)=-e(1,2) e(3,1)=-e(1,3) e(3,2)=-e(2,3) return end subroutine dsa(tr,ti,tsr,tsi,tar,tai) implicit none real*8 tr(3,3),ti(3,3),tsr(3,3),tsi(3,3),tar(3,3),tai(3,3) integer*4 a,b do 1 a=1,3 do 1 b=1,3 tsr(a,b)=0.5d0*(tr(a,b)+tr(b,a)) tsi(a,b)=0.5d0*(ti(a,b)+ti(b,a)) tar(a,b)=0.5d0*(tr(a,b)-tr(b,a)) 1 tai(a,b)=0.5d0*(ti(a,b)-ti(b,a)) return end subroutine abab(tr,ti,ar,ai,br,bi) c t = a_ab b*_ab implicit none integer*4 a real*8 tr,ti,ar(3,3),ai(3,3),br(3,3),bi(3,3) tr=0.0d0 ti=0.0d0 do 1 a=1,3 tr=tr+ar(a,1)*br(a,1)+ai(a,1)*bi(a,1) 1 +ar(a,2)*br(a,2)+ai(a,2)*bi(a,2) 1 +ar(a,3)*br(a,3)+ai(a,3)*bi(a,3) 1 ti=ti-ar(a,1)*bi(a,1)+ai(a,1)*br(a,1) 1 -ar(a,2)*bi(a,2)+ai(a,2)*br(a,2) 1 -ar(a,3)*bi(a,3)+ai(a,3)*br(a,3) return end subroutine ru(io,A,n,m) c matrix n x m implicit none REAL*8 A(n,m) INTEGER*4 io,n,m,N1,N3,LN,J N1=1 1 N3=MIN(N1+4,m) read(io,*) DO 130 LN=1,N 130 read(io,*)A(LN,N1),(A(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.m)GOTO 1 return end subroutine t0(N,ALPHA,ALPHAi,G,Gi,GGr,GGi,A,Ai,AAr,AAi) implicit none integer*4 N real*8 ALPHA(N,3,3),ALPHAi(N,3,3),G(N,3,3),Gi(N,3,3),A(N,3,3,3), 1GGr(N,3,3),GGi(N,3,3),AAr(N,3,3,3),Ai(N,3,3,3),AAi(N,3,3,3) ALPHA=0.0d0 ALPHAi=0.0d0 G =0.0d0 Gi =0.0d0 A=0.0d0 Ai=0.0d0 GGr=0.0d0 GGi =0.0d0 AAr=0.0d0 AAi=0.0d0 return end subroutine chkroot(f,NROOT) implicit none integer*4 NROOT,found,nd,i character*(*) f character*80 s nd=0 found=0 open(9,file=f,status='old') 1 read(9,900,end=88,err=88)s 900 format(a80) if(s(2:15).eq.'Excited State '.and.s(16:16).ne.'s')then do 20 i=1,79 20 if(s(i:i) .eq.':') read(s(15:i-1),*)nd endif if(s(2:29).eq.'This state for optimization ')found=nd goto 1 88 close(9) write(6,600)found,nd 600 format(i6,' of',i6,' excited states in ',$) do 21 i=1,len(f) 21 if(f(i:i).ne.' ')write(6,601)f(i:i) 601 format(a1,$) write(6,*) NROOT=found return end SUBROUTINE WTQi(N,NINT,ALPHA,A,G,ALPHAi,Ai,Gi,av,E,WMIN,WMAX) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) integer*4 o DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3), 1ALPHAi(N,3,3),Gi(N,3,3),Ai(N,3,3,3),av(N,3,3),E(*) c CM=219470.0d0 im=0 DO 4 II=1,NINT 4 if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)im=im+1 open(22,file='FILE.Q.TTT') do 1 o=1,4 if(o.eq.1)WRITE(22,2001)im 2001 FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2 ' mode e(cm-1) jx jy jz') if(o.eq.2)WRITE(22,2002) 2002 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' mode e(cm-1) jx(Bx) jy(By) jz(Bz)') if(o.eq.3)WRITE(22,2003) 2003 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' mode e(cm-1) kx ky kz') if(o.eq.4)WRITE(22,2004) 2004 FORMAT('The velocity form of alpha',/, 1' mode e(cm-1) vx vy vz') im=0 DO 1 II=1,NINT if(E(II)*CM.GE.WMIN.AND.E(II).LE.WMAX)then im=im+1 write(22,221)im,E(II)*CM 221 format(i5,f12.2) if(o.eq.1)then DO 11 I=1,3 11 WRITE(22,2005)I,(ALPHA(II,I,J),J=1,3),(ALPHAi(II,I,J),J=1,3) 2005 FORMAT(I3,6F13.7) endif if(o.eq.2)then DO 21 I=1,3 21 WRITE(22,2005)I,(G(II,I,J),J=1,3),(Gi(II,I,J),J=1,3) c last inner index magnetic endif if(o.eq.3)then DO 31 I=1,3 DO 31 J=1,3 31 WRITE(22,2006)I,J,(A(II,I,J,K),K=1,3),(Ai(II,I,J,K),K=1,3) 2006 FORMAT(2I3,6F13.7) endif if(o.eq.4)then DO 51 I=1,3 51 WRITE(22,2005)I,(av(II,I,J),J=1,3),(av(II,I,J),J=1,3) endif endif 1 continue close(22) RETURN END subroutine tq(q,q0) implicit none real*8 q(3,3),q0(6) q(1,1)=q0(1) q(1,2)=q0(2) q(1,3)=q0(3) q(2,2)=q0(4) q(2,3)=q0(5) q(3,3)=q0(6) q(2,1)=q(1,2) q(3,1)=q(1,3) q(3,2)=q(2,3) return end subroutine tqi(qd,qqi,iq) implicit none integer*4 iq real*8 qd(3,3),qqi(*) qd(1,1)=qqi(1+6*(iq-1)) qd(1,2)=qqi(2+6*(iq-1)) qd(1,3)=qqi(3+6*(iq-1)) qd(2,2)=qqi(4+6*(iq-1)) qd(2,3)=qqi(5+6*(iq-1)) qd(3,3)=qqi(6+6*(iq-1)) qd(2,1)=qd(1,2) qd(3,1)=qd(1,3) qd(3,2)=qd(2,3) return end subroutine fl(ni,s,sx) implicit none integer*4 i,ni character*80 sx(ni) character*(*) s open(9,file=s,status='old') do 1 i=1,ni 1 read(9,90)sx(i) 90 format(a80) close(9) return end subroutine rdusch(NQ,N,JT,K,A,B,C,D,E,W,WP,iw,ip,f) IMPLICIT none character*(*) f integer*4 NQ,N,ip(10),iw real*8 JT(NQ,NQ),K(NQ),A(N,N),B(N),C(N,N),D(N),E(N,N),W(N),WP(N) character*80 s80 iw=0 ip=0 JT=0.0d0 K=0.0d0 A=0.0d0 B=0.0d0 C=0.0d0 D=0.0d0 E=0.0d0 W=0.0d0 WP=0.0d0 open(20,file=f,status='old') read(20,*)NQ write(6,*)' number of modes :',NQ write(6,*)' matrix dimension:', N 2 read(20,900,end=88,err=88)s80 900 format(a80) if(s80(2:18).eq.'Duschinsky Matrix') 1 call rdm(20,NQ,NQ,JT,ip(1),s80) if(s80(2:13).eq.'Shift Vector') call rdv(20,NQ,K ,ip(2)) if(s80(2: 9).eq.'A Matrix' ) call rdm(20,NQ,N ,A ,ip(3),s80) if(s80(2: 9).eq.'B Vector' ) call rdv(20,NQ,B , ip(4)) if(s80(2: 9).eq.'C Matrix' ) call rdm(20,NQ,N ,C ,ip(5),s80) if(s80(2: 9).eq.'D Vector' ) call rdv(20,NQ,D , ip(6)) if(s80(2: 9).eq.'E Matrix' ) call rdm(20,NQ,N ,E ,ip(7),s80) if(s80(2: 9).eq.'Ground w' ) call rdw(20,NQ,W ,iw) if(s80(2:10).eq.'Excited w') call rdw(20,NQ,WP,iw) goto 2 88 close(20) write(6,*)f//' read for Duschinsky parameters' return end subroutine rduschw(N,e,f,s) c read ground or excited vibrational frequencies c s='Ground w' or s ='Excited w' IMPLICIT none character*(*) f,s integer*4 NQ,N,iw real*8 e(N) character*80 s80 open(20,file=f,status='old') read(20,*)NQ iw=0 2 read(20,900,end=88,err=88)s80 900 format(a80) if(s80(2: 9).eq.s)then call rdw(20,NQ,e,iw) goto 88 endif if(s80(2:10).eq.s)then call rdw(20,NQ,e,iw) goto 88 endif goto 2 88 close(20) if(iw.eq.0)call report(s//' not found in '//f) write(6,600) 600 format(' Frequecies read') return end subroutine mt(AT,A,NQ) c matrix transposition AT = A^T implicit none integer*4 NQ,i,j real*8 A(NQ,NQ),AT(NQ,NQ) do 1 i=1,NQ do 1 j=1,NQ 1 AT(i,j)=A(j,i) return end subroutine mm(A,B,C,NQ) c matrix multiplication A = B x C implicit none integer*4 NQ,i,j,k real*8 A(NQ,NQ),B(NQ,NQ),C(NQ,NQ) do 1 i=1,NQ do 1 j=1,NQ A(i,j)=0.0d0 do 1 k=1,NQ 1 A(i,j)=A(i,j)+B(i,k)*C(k,j) return end subroutine INV(a,ai,n,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(n,n),a(n,n) real*8,allocatable::e(:,:) C TOL=1.0d-10 allocate(e(n,2*n)) 10000 IERR=0 DO 1 ii=1,n DO 1 jj=1,n e(ii,jj)=a(ii,jj) 1 e(ii,jj+n)=0.0D0 do 13 ii=1,n 13 e(ii,ii+n)=1.0D0 c DO 2 ii=1,n-1 if (ABS(e(ii,ii)).LE.TOL) then DO 3 io=ii+1,n 3 if (ABS(e(io,ii)).GT.TOL) goto 11 IERR=1 write(6,*)ii write(6,*)'tol = ',tol tol=tol*0.50d0 if(tol.gt.1.0d-20)goto 10000 RETURN c 11 CONTINUE DO 4 kk=1,2*n w=e(ii,kk) e(ii,kk)=e(io,kk) 4 e(io,kk)=w ENDIF eii=e(ii,ii) DO 5 jj=ii+1,n e1=e(jj,ii)/eii DO 6 kk=ii+1, 2*n 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e1 5 e(jj,ii)=0.0D0 2 CONTINUE c DO 7 i2=n,2,-1 eii=e(i2,i2) DO 7 j2=i2-1,1,-1 e1=e(j2,i2)/eii DO 9 kk=1, n 9 e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1 7 e(j2,i2)=0.0d0 c DO 10 ii=1,n ei=1.0d0/e(ii,ii) DO 12 jj=1,n 12 ai(ii,jj)=e(ii,jj+n)*ei 10 CONTINUE c RETURN END subroutine s3(ni,i,NQ,K1,U1,V1,X1,Y1,ki,U,V,X,Y) implicit none integer*4 ni,i,NQ,a,b real*8 K1(NQ),U1(NQ,NQ),V1(NQ,NQ),X1(NQ,NQ),Y1(NQ,NQ), 1ki(ni,NQ),U(ni,NQ,NQ),V(ni,NQ,NQ),X(ni,NQ,NQ),Y(ni,NQ,NQ) do 1 a=1,NQ ki(i,a)=K1(a) do 1 b=1,NQ U(i,a,b)=U1(a,b) V(i,a,b)=V1(a,b) X(i,a,b)=X1(a,b) 1 Y(i,a,b)=Y1(a,b) return end subroutine m3(N,NQ,ki,U,V,X,Y,W,WP,sj,wz,block) implicit none integer*4 N,NQ,i,j,k,block real*8 ki(NQ),wj,wk,WP(N),U(NQ,NQ),V(NQ,NQ),X(NQ,NQ),Y(NQ,NQ), 1sj(NQ,NQ),W(N),wz real*8,allocatable::su(:,:),sv(:,:) U=0.0d0 V=0.0d0 X=0.0d0 Y=0.0d0 ki=0.0d0 allocate(su(NQ,NQ),sv(NQ,NQ)) wz=0.0d0 do 1 j=block+1,NQ wz=wz+WP(j)/2.0d0 do 1 k=block+1,NQ sv(j,k)=0.0d0 su(j,k)=0.0d0 do 1 i=block+1,NQ su(j,k)=su(j,k)+ sj(i,k)*sj(i,j) 1 sv(j,k)=sv(j,k)+WP(i)**2*sj(i,k)*sj(i,j) do 2 j=block+1,NQ wj=W(j) ki(j)=dsqrt(2.0d0*wj) do 2 k=block+1,NQ wk=W(k) U(j,k)=0.25d0*su(j,k)*dsqrt(wj*wk) V(j,k)=0.25d0*sv(j,k)/dsqrt(wj*wk) X(j,k)=V(j,k)-U(j,k) 2 Y(j,k)=V(j,k)+U(j,k) return end subroutine sl(i,ni,NQ,N,uu1,ff1,gg1,uu,ff,gg,W) implicit none integer*4 ni,NQ,N,i,a,iq real*8 uu1(81),uu(ni,81),ff1(81,NQ),ff(ni,81,NQ),gg1(81,NQ), 1gg(ni,81,NQ),W(N) do 1 a=1,81 uu(i,a)=uu1(a) do 1 iq=1,NQ ff(i,a,iq)=ff1(a,iq)/dsqrt(2.0d0*W(iq)) 1 gg(i,a,iq)=gg1(a,iq) return end subroutine sc(N,NQ,u0,m0,q0,ddi,rri,qqi,uu,ff,gg) c ii ji index: c 1.. 9 alpha_ab c 10..18 G_ab c 19..27 GC_ab c 27..54 A_abc c 25..81 AC_abc implicit none integer*4 N,NQ,i,j,a,b,c,iq real*8 uu(81),q(3,3),u0(3),q0(6),m0(3),ddi(3*N),rri(3*N),qqi(6*N), 1qd(3,3),ud(3),md(3),ff(81,NQ),gg(81,NQ) c constants: call tq(q,q0) i=0 j=0 do 1 a=1,3 do 1 b=1,3 i=i+1 uu(i )=u0(a)* u0(b) uu(i+ 9)=u0(a)*(-m0(b)) uu(i+18)=m0(a)* u0(b) do 1 c=1,3 j=j+1 uu(j+27)=u0(a)* q(b,c) 1 uu(j+54)=u0(a)* q(b,c) c derivatives, ff-symmetric,gg-assymetric: do 2 iq=1,NQ call tqi(qd,qqi,iq) call setu(ud,ddi,iq) call setu(md,rri,iq) j=0 i=0 do 2 a=1,3 do 2 b=1,3 i=i+1 ff(i ,iq)=ud(a)* u0(b) + u0(a)* ud(b) ff(i+ 9,iq)=ud(a)*(-m0(b)) +u0(a)*(-md(b)) ff(i+18,iq)=md(a)* u0(b) +m0(a)* ud(b) gg(i ,iq)= u0(a)* ud(b) gg(i+ 9,iq)= u0(a)*(-md(b)) gg(i+18,iq)= m0(a)* ud(b) do 2 c=1,3 j=j+1 ff(j+27,iq)=ud(a) *q(b,c) + u0(a)*qd(b,c) ff(j+54,iq)=qd(b,c)*u0(a) + q(b,c)*ud(a) gg(j+27,iq)= u0(a)*qd(b,c) 2 gg(j+54,iq)= q(b,c)*ud(a) return end subroutine mkf0(i,ni,we1,we,w,Gau,f0,g0,f1,g1) implicit none integer*4 i,ni real*8 we1,we(ni),w,Gau,f0(ni),g0(ni),f1(ni),g1(ni),down we(i)=we1 down=((we1-w)*(we1+w))**2+(2.0d0*we1*Gau)**2 f0(i)=(we1-w)*(we1+w)/down g0(i)=-2.0d0*we1*Gau/down down=((we1-w)**2+Gau**2)**2 f1(i)=1.0d0/((we1-w)**2+Gau**2) g1(i)=0.0d0 return end subroutine ztp(lr,li,gr,gi,gcr,gci,ar,ai,acr,aci) implicit none real*8 lr(3,3),li(3,3),gr(3,3),gi(3,3),gcr(3,3), 1gci(3,3),ar(3,3,3),ai(3,3,3),acr(3,3,3),aci(3,3,3) lr=0.0d0 li=0.0d0 gr=0.0d0 gi=0.0d0 gcr=0.0d0 gci=0.0d0 ar=0.0d0 ai=0.0d0 acr=0.0d0 aci=0.0d0 return end subroutine ctp(we,w,lr,li,gr,gi,gcr,gci,ar,ai,acr,aci, 1f0,g0,f1,g1,iint,jint) implicit none integer*4 i,j,a,b,c real*8 we,w,f0,g0,f1,g1,iint(81), 1jint(81),lr(3,3),li(3,3),gr(3,3),gi(3,3),gcr(3,3), 1gci(3,3),ar(3,3,3),ai(3,3,3),acr(3,3,3),aci(3,3,3),tw data tw/2.0d0/ i=0 j=0 do 4 a=1,3 do 4 b=1,3 i=i+1 lr(a,b)=lr(a,b)+ tw*we*iint(i)*f0-jint(i)*f1 li(a,b)=li(a,b)+ tw*we*iint(i)*g0-jint(i)*g1 c Re = - Im * Im gr(a,b) =gr(a,b) -tw*w*iint(i+ 9)*g0+jint(i+ 9)*g1 c Im = Im * Re gi(a,b) =gi(a,b) +tw*w*iint(i+ 9)*f0-jint(i+ 9)*f1 gcr(a,b) =gcr(a,b)-tw*w*iint(i+18)*g0+jint(i+18)*g1 gci(a,b) =gci(a,b)+tw*w*iint(i+18)*f0-jint(i+18)*f1 do 4 c=1,3 j=j+1 ar(a,b,c) =ar( a,b,c)+tw*we*iint(j+27)*f0-jint(j+27)*f1 ai(a,b,c) =ai( a,b,c)+tw*we*iint(j+27)*g0-jint(j+27)*g1 acr(a,b,c)=acr(a,b,c)+tw*we*iint(j+54)*f0-jint(j+54)*f1 4 aci(a,b,c)=aci(a,b,c)+tw*we*iint(j+54)*g0-jint(j+54)*g1 return end function plus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) c A |nkp> = a+k |n>, A=sqrt(plus) implicit none integer*4 nkp,nc,nkpcen(*),nkpexc(*),ncen(*),nexc(*),k,a,plus nkp=nc do 2 a=1,nc nkpexc(a)=nexc(a) 2 nkpcen(a)=ncen(a) do 1 a=1,nc if(k.eq.ncen(a))then nkpexc(a)=nexc(a)+1 plus=nkpexc(a) return endif 1 continue plus=1 nkp=nc+1 nkpexc(nkp)=1 nkpcen(nkp)=k return end function plusn(k,nc,ncen,nexc) c A |nkp> = a+k |n>, A=sqrt(plusn) implicit none integer*4 nc,ncen(*),nexc(*),k,a,plusn do 1 a=1,nc if(k.eq.ncen(a))then plusn=nexc(a)+1 return endif 1 continue plusn=1 return end function minus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) c A |nkp> = a-k |n>, A=sqrt(minus) implicit none integer*4 nkp,nc,nkpcen(*),nkpexc(*),ncen(*), 1nexc(*),k,a,minus,b nkp=nc do 3 a=1,nc nkpexc(a)=nexc(a) 3 nkpcen(a)=ncen(a) do 1 a=1,nc if(k.eq.ncen(a))then nkpexc(a)=nexc(a)-1 minus=nexc(a) if(nkpexc(a).eq.0)then do 2 b=a,nkp-1 nkpexc(b)=nkpexc(b+1) 2 nkpcen(b)=nkpcen(b+1) nkp=nkp-1 endif return endif 1 continue minus=0 return end function minusn(k,nc,ncen,nexc) c A |nkp> = a-k |n>, A=sqrt(minusn) implicit none integer*4 nc,ncen(*),nexc(*),k,a,minusn do 1 a=1,nc if(k.eq.ncen(a))then minusn=nexc(a) return endif 1 continue minusn=0 return end subroutine mkijf(ni,N,NA,NQ,W,nc,ncen,nexc,uu,ff,gg, 1X,Y,U,V,ki,wmin,dw,pf,roa,pis,lstat, 1ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi, lluseA,luseG,ldo,we,wexc,f0,g0,f1,g1,b,iwr) c should be fast version implicit none integer*4 ni,N,NQ,maxcen,nn,ii,j,k,akp,akm,app,amm,c,ie, 1ap,NA,minusn,plusn,ns0,nik,akpp,akmm,b,anm,iwr,at parameter (ns0=19) parameter (maxcen=10) INTEGER*4 nc,ncen(*),nexc(*), 1nkp,nkpcen(maxcen),nkpexc(maxcen), 1nkpp,nkppcen(maxcen),nkppexc(maxcen), 1nt,tcen(maxcen),texc(maxcen), 1nkm,nkmcen(maxcen),nkmexc(maxcen), 1nkmm,nkmmcen(maxcen),nkmmexc(maxcen), 1nmm,mmcen(maxcen),mmexc(maxcen), 1nm,nmcen(maxcen),nmexc(maxcen), 1npp,ppcen(maxcen),ppexc(maxcen), 1plus,minus,mpp real*8 W(N),U(ni,NQ,NQ),V(ni,NQ,NQ),X(ni,NQ,NQ),Y(ni,NQ,NQ), 1uu(ni,81),ff(ni,81,NQ),gg(ni,81,NQ),wn,tw,t,tj,ki(ni,NQ), 1wmin,dw,pf,wk,wj,roa(2,ns0,NA), 1pis,anik,tk,dakp, 1ALPHA(N,3,3),ALPHAi(N,3,3),G(N,3,3),Gi(N,3,3), 1A(N,3,3,3),Ai(N,3,3,3),GGr(N,3,3),GGi(N,3,3),AAr(N,3,3,3), 1AAi(N,3,3,3),CM,iint(81),jint(81), 1alr(3,3),ali(3,3),gtr(3,3),gti(3,3),gcr(3,3),gci(3,3), 1atr(3,3,3),ati(3,3,3),acr(3,3,3),aci(3,3,3), 1we(ni),wexc,f0(ni),g0(ni),f1(ni),g1(ni),z0,two logical lstat,luseA,luseG,ldo(ns0) real*8,allocatable::sk(:,:),sy(:),syk(:,:),sxk(:,:) data z0,two,CM/0.0d0,2.0d0,219470.0d0/ c=b+1 allocate(sk(ni,NQ),sy(ni),sxk(ni,81),syk(ni,81)) sk=z0 sy=z0 do 25 ie=1,ni do 25 k=c,NQ nik=minusn(k,nc,ncen,nexc) sy(ie)=sy(ie)+Y(ie,k,k)*dble(nik+nik+1) do 25 j=c,NQ 25 sk(ie,k)=sk(ie,k)-ki(ie,j)*V(ie,j,k) sk=sk+sk do 2 k=c,NQ sxk=z0 syk=z0 do 24 j=c,NQ if(j.ne.k)then t=dsqrt(two/w(j)) tj=t*dble(minusn(j,nc,ncen,nexc)) t=tj+t do 23 ie=1,ni do 23 ii=1,81 sxk(ie,ii)=sxk(ie,ii)+t *FF(ie,ii,j)*X(ie,j,k) 23 syk(ie,ii)=syk(ie,ii)+tj*FF(ie,ii,j)*Y(ie,j,k) endif 24 continue wk=w(k) tw=dsqrt(wk+wk) akm=minus(nkm,nkmcen,nkmexc,k,nc,ncen,nexc) akmm=minus(nkmm,nkmmcen,nkmmexc,k,nkm,nkmcen,nkmexc)*akm nik=plusn(k,nkm,nkmcen,nkmexc) anik=dble(nik) akp=plus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) akpp=plus(nkpp,nkppcen,nkppexc,k,nkp,nkpcen,nkpexc)*akp dakp=dsqrt(dble(akp)) c a+k |i> ap=int((wk-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then c zero transition polarizabilities: call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) c add all electronic states: do 261 ie=1,ni do 3 ii=1,81 tk=z0 do 26 j=c,NQ t=dsqrt(8.0d0/W(j)) 26 tk=tk+t*U(ie,j,k)*gg(ie,ii,j) iint(ii)=ff(ie,ii,k)*dakp 3 jint(ii)=(uu(ie,ii)*sk(ie,k)+tk+tw*ff(ie,ii,k)*(sy(ie) 1 +X(ie,k,k)*(anik+two))/tw+sxk(ie,ii)+syk(ie,ii) )*dakp c add to transition polarizabilities 261 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) c add the static derivatives by the pis factor: if(lstat)call ads1(wexc,wk,k,NQ,alr,ali,ALPHA,ALPHAi,gtr,gti,G, 1 Gi,gcr,gci,GGr,GGi,atr,ati,A,Ai,acr,aci,AAr,AAi,pis,iwr) c save intensities (and forget this transition) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)call wssr(' to : ',k,wk*CM,nkp,nkpcen,nkpexc) endif c a+k a+k |i>: ap=int((wk+wk-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=dsqrt(dble(akpp)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 262 ie=1,ni do 34 ii=1,81 34 jint(ii)=(sk(ie,k)*ff(ie,ii,k)+uu(ie,ii)*X(ie,k,k))*t 262 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1) 1 call wssr(' to : ',k,(wk+wk)*CM,nkpp,nkppcen,nkppexc) endif c a+k a+k a+k |i>: ap=int((wk+wk+wk-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=dsqrt(dble(plusn(k,nkpp,nkppcen,nkppexc)*akpp)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 263 ie=1,ni do 36 ii=1,81 36 jint(ii)=ff(ie,ii,k)*X(ie,k,k)*t 263 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,k,nkpp,nkppcen,nkppexc)*akpp call wssr(' to : ',k,3.0d0*wk*CM,nt,tcen,texc) endif endif do 2 j=c,NQ if(j.ne.k)then wj=w(j) app=plus(npp,ppcen,ppexc,k,nkp,nkpcen,nkpexc)*akp amm=minus(nmm,mmcen,mmexc,j,nkm,nkmcen,nkmexc)*akm c a+j a-k |i> j<>k : ap=int((wj-wk-wmin+dw)/dw) if(akm.ne.0.and.ap.ge.1.and.ap.le.NA)then call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) t=dsqrt(dble(plusn(j,nkm,nkmcen,nkmexc)*akm)) do 264 ie=1,ni do 35 ii=1,81 35 jint(ii)=(sk(ie,k)*ff(ie,ii,j)+sk(ie,j)*ff(ie,ii,k) 1 +two*uu(ie,ii)*Y(ie,j,k))*t 264 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,j,nkm,nkmcen,nkmexc) call wssr(' to : ',k,(wj-wk)*CM,nt,tcen,texc) endif endif c a+j a+k |i> j>k : ap=int((wj+wk-wmin+dw)/dw) if(j.gt.k.and.ap.ge.1.and.ap.le.NA)then t=dsqrt(dble(app)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 265 ie=1,ni do 33 ii=1,81 33 jint(ii)=(sk(ie,k)*ff(ie,ii,j)+sk(ie,j)*ff(ie,ii,k) 1 +two*uu(ie,ii)*X(ie,j,k))*t 265 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then call wssr(' to : ',k,(wk+wj)*CM,npp,ppcen,ppexc) endif endif c a+j a+k a+k |i> j<>k: ap=int((wj+wk+wk-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=dsqrt(dble(plusn(j,nkpp,nkppcen,nkppexc)*akpp)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 266 ie=1,ni do 37 ii=1,81 37 jint(ii)= 1 (ff(ie,ii,j)*X(ie,k,k)+two*ff(ie,ii,k)*X(ie,j,k))*t 266 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,j,nkpp,nkppcen,nkppexc) call wssr(' to : ',k,(wk+wk+wj)*CM,nt,tcen,texc) endif endif c a+j a-k a-k |i> j<>k: ap=int((wj-wk-wk-wmin+dw)/dw) if(akmm.ne.0.and.ap.ge.1.and.ap.le.NA)then t=dsqrt(dble(plusn(j,nkmm,nkmmcen,nkmmexc)*akmm)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 267 ie=1,ni do 38 ii=1,81 38 jint(ii)= 1 (ff(ie,ii,j)*X(ie,k,k)+two*ff(ie,ii,k)*Y(ie,j,k))*t 267 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,j,nkmm,nkmmcen,nkmmexc) call wssr(' to : ',k,(-wk+wj-wk)*CM,nt,tcen,texc) endif endif c a-j a+k a+k |i> j<>k: ap=int((-wj+wk+wk-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then mpp=minusn(j,nkpp,nkppcen,nkppexc)*akpp if(mpp.ne.0)then t=dsqrt(dble(mpp)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 268 ie=1,ni do 39 ii=1,81 39 jint(ii)= 1 (ff(ie,ii,j)*X(ie,k,k)+two*ff(ie,ii,k)*Y(ie,j,k))*t 268 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=minus(nt,tcen,texc,j,nkpp,nkppcen,nkppexc) call wssr(' to : ',k,(wk-wj+wk)*CM,nt,tcen,texc) endif endif endif endif c j<>k 2 continue c separate triple sums with special filtering c a+n a+j a+k |i> nn > j > k: do 52 k=c,NQ wk=W(k) if(int((wk-wmin+dw)/dw).le.NA)then akp=plus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) do 53 j=k+1,NQ wj=W(j) if(int((wj+wk-wmin+dw)/dw).le.NA)then app=plus(npp,ppcen,ppexc,k,nkp,nkpcen,nkpexc)*akp do 54 nn=j+1,NQ wn=w(nn) ap=int((wn+wk+wj-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=2.0d0*dsqrt(dble(plusn(nn,npp,ppcen,ppexc)*app)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 269 ie=1,ni do 40 ii=1,81 40 jint(ii)=(ff(ie,ii,nn)*X(ie,j,k)+ff(ie,ii,j)*X(ie,nn,k)+ 1 ff(ie,ii,k)*X(ie,nn,j))*t 269 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,nn,npp,ppcen,ppexc) call wssr(' to : ',k,(wk+wj+wn)*CM,nt,tcen,texc) endif endif 54 continue endif 53 continue endif 52 continue c a+n a-j a-k |i> nn <> j, nn <> k, j > k: do 62 k=c,NQ akm=minus(nkm,nkmcen,nkmexc,k,nc,ncen,nexc) if(akm.ne.0)then wk=W(k) do 63 j=k+1,NQ amm=minus(nmm,mmcen,mmexc,j,nkm,nkmcen,nkmexc)*akm if(amm.ne.0)then wj=W(j) do 64 nn=c,NQ if(nn.ne.j.and.nn.ne.k)then wn=w(nn) ap=int((wn-wk-wj-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=2.0d0*dsqrt(dble(plusn(nn,nmm,mmcen,mmexc)*amm)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 270 ie=1,ni do 41 ii=1,81 41 jint(ii)=(ff(ie,ii,nn)*X(ie,j,k)+ff(ie,ii,j)*Y(ie,nn,k) 1 +ff(ie,ii,k)*Y(ie,nn,j))*t 270 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,nn,nmm,mmcen,mmexc) call wssr(' to : ',k,(-wk-wj+wn)*CM,nt,tcen,texc) endif endif endif 64 continue endif 63 continue endif 62 continue c a-n a+j a+k |i> nn <> j, nn <> k, j > k: do 72 nn=c,NQ anm=minus(nm,nmcen,nmexc,nn,nc,ncen,nexc) if(anm.ne.0)then wn=w(nn) do 73 k=c,NQ if(k.ne.nn)then akp=plus(nkp,nkpcen,nkpexc,k,nm,nmcen,nmexc)*anm wk=W(k) do 74 j=k+1,NQ if(j.ne.nn)then wj=W(j) ap=int((-wn+wk+wj-wmin+dw)/dw) if(ap.ge.1.and.ap.le.NA)then t=2.0d0*dsqrt(dble(plusn(j,nkp,nkpcen,nkpexc)*akp)) call ztp(alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci) do 271 ie=1,ni do 42 ii=1,81 42 jint(ii)=(ff(ie,ii,nn)*X(ie,j,k)+ff(ie,ii,j)*Y(ie,nn,k) 1 +ff(ie,ii,k)*Y(ie,nn,k))*t 271 call ctp(we(ie),wexc,alr,ali,gtr,gti,gcr,gci,atr,ati,acr,aci, 1 f0(ie),g0(ie),f1(ie),g1(ie),iint,jint) call wrram1(iwr,NA,ap,wexc,alr,ali,gtr,gti,gcr,gci,atr, 1 ati,acr,aci,luseA,luseG,ldo,roa,pf) if(iwr.gt.1)then at=plus(nt,tcen,texc,j,nkp,nkpcen,nkpexc) call wssr(' to : ',k,(wk+wj-wn)*CM,nt,tcen,texc) endif endif endif 74 continue endif 73 continue endif 72 continue return end subroutine rds(io,NQ,N,M,ic) IMPLICIT none integer*4 NQ,N,N1,N3,io,LN,J,ic real*8 M(N,N) ic=ic+1 read(io,*) N1=1 1 N3=min(N1+4,NQ) read(io,*) DO 130 LN=1,N 130 READ(io,*)M(LN,N1),(M(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.NQ)GOTO 1 return end subroutine rdv(io,N,V,ic) IMPLICIT none integer*4 N,io,LN,ic real*8 V(*) ic=ic+1 read(io,*) read(io,*) DO 130 LN=1,N 130 READ(io,*)V(LN),V(LN) return end subroutine rdw(io,N,e,iw) IMPLICIT none integer*4 i,io,iw,N real*8 e(*),CM read(io,*) read(io,100)(e(i),i=1,N) 100 format(6f11.2) CM=219474.63d0 do 1 i=1,N 1 e(i)=e(i)/CM iw=iw+1 return end subroutine rdm(io,NQ,N,M,ic,s) IMPLICIT none integer*4 NQ,N,N1,N3,io,LN,J,ic real*8 M(N,N) character*(*) s write(6,*)s(1:len(s)/2) ic=ic+1 read(io,*) N1=1 1 N3=min(N1+4,NQ) read(io,*) DO 130 LN=1,NQ 130 READ(io,*)M(LN,N1),(M(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.NQ)GOTO 1 return end subroutine rduschs(N,NQ,s,f,ff) IMPLICIT none character*(*) f,ff integer*4 NQ,N,ic,I,J real*8 s(N,N),SFAC character*80 s80 SFAC=0.0234280d0 open(20,file=f,status='old') ic=0 read(20,*)NQ 2 read(20,900,end=88,err=88)s80 900 format(a80) if(s80(2:16).eq.ff)then call rds(20,NQ,N,s,ic) DO 5 I=1,N DO 5 J=1,NQ 5 s(I,J)=s(I,J)*SFAC close(20) write(6,*)f write(6,*)'S read, ',NQ,' modes' return endif goto 2 88 close(20) write(6,*)f if(ic.eq.0)call report('S not found') end function jjs(s) implicit none integer*4 jjs,jj character*(*) s do 1001 jj=len(s),1,-1 1001 if(s(jj:jj).ne.' ')goto 1002 1002 jjs=jj return end subroutine control(N,NQ,s) integer*4 N,NQ,i,j,k,ity,ii,ix,nat real*8 s(N,N) real*8 ,allocatable::t(:,:),ss(:,:),m(:) real amas(89) data amas/1.007825,4.002603, 2 6.941, 9.012, 10.810,12.000,14.003074,15.9949,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ allocate(t(N,N),ss(N,N),m(N/3)) write(6,*)N,NQ open(9,file='ground/FILE.X',status='old') read(9,*) read(9,*)nat do 3 i=1,nat read(9,*)ity 3 m(i)=dble(amas(ity)*1822.89) close(9) do 4 i=1,nat do 4 ix=1,3 ii=ix+3*(i-1) do 4 j=1,NQ 4 ss(ii,j)=s(ii,j)*sqrt(m(i)) do 1 i=1,NQ do 1 j=1,NQ t(i,j)=0.0d0 do 1 k=1,N 1 t(i,j)=t(i,j)+ss(k,i)*ss(k,j) call wpart(N,t,NQ,8) do 2 i=1,N do 2 j=1,N t(i,j)=0.0d0 do 2 k=1,NQ 2 t(i,j)=t(i,j)+ss(i,k)*ss(j,k) call wpart(N,t,N,8) return end subroutine wpart(Ni,t,N,it) integer*4 N,it,j,i,Ni real*8 t(Ni,Ni) write(6,601)(i,i=1,it),(N-it+i,i=1,it) 601 format(16i7) do 5 j=1,it 5 write(6,600)(t(j,i),i=1,it),(t(j,N-it+i),i=1,it) do 6 j=1,it 6 write(6,600)(t(N-it+j,i),i=1,it),(t(N-it+j,N-it+i),i=1,it) 600 format(16f7.4) return end subroutine mksv(NQ,N,K,gradi,WP,J,W,iwr,b) implicit none integer*4 NQ,N,i,p,iwr,b real*8 K(NQ),WP(N),gradi(N),J(NQ,NQ),W(N),CM,a CM=219474.63d0 K=0.0d0 do 1 i=b+1,NQ do 1 p=b+1,NQ 1 K(i)=K(i)+J(i,p)*gradi(p)/WP(p)**2 a=0.0d0 do 3 i=b+1,NQ 3 a=a+K(i)**2 a=dsqrt(a) if(iwr.ne.0)then write(6,600) 600 format(' K-shifts:') do 2 i=b+1,NQ write(6,601)W(i)*CM,K(i) 601 format(f9.1,e11.3,$) 2 if(mod(i,4).eq.0)write(6,*) write(6,*) endif write(6,602)a 602 format(' K-norm:',e11.3) return end subroutine ede(nat,d,v,a,q,g) implicit none integer*4 nat,nae,i,ii,il,j,ix,jx real*8 d(9*nat),a(9*nat),v(9*nat),q(18*nat),g(3*nat) integer*4,allocatable::ial(:) allocate(ial(nat)) open(9,file='EA.LST',status='old') read(9,*)nae if(abs(nae).gt.nat)call report('too many atoms') read(40,*)(ial(i),i=1,abs(nae)) write(6,*)'Deleting derivatives according to EA.LST:' if(nae.gt.0)then write(6,*)nae,' atoms deleted' do 1 il=1,nae i=ial(il) do 1 ix=1,3 g(ix+3*(i-1))=0.0d0 ii=ix+3*(i-1) do 1 jx=1,3 d(jx +3*(ii-1))=0.0d0 v(jx +3*(ii-1))=0.0d0 a(jx +3*(ii-1))=0.0d0 q(jx +6*(ii-1))=0.0d0 1 q(jx+3+6*(ii-1))=0.0d0 else do 2 i=1,nat do 22 j=1,abs(nae) 22 if(i.eq.ial(j))goto 2 do 3 ix=1,3 g(ix+3*(i-1))=0.0d0 ii=ix+3*(i-1) do 3 jx=1,3 d(jx +3*(ii-1))=0.0d0 v(jx +3*(ii-1))=0.0d0 a(jx +3*(ii-1))=0.0d0 q(jx +6*(ii-1))=0.0d0 3 q(jx+3+6*(ii-1))=0.0d0 2 continue write(6,*)abs(nae),' atoms conserved' endif return end function mfc(N,m) c N! c mfc = -------- c m!(N-m)! implicit none integer*4 N,m,mf,u,mfc,i mf=1 do 1 i=1,m 1 mf=mf*i u=1 do 2 i=N-m+1,N 2 u=u*i mfc=u/mf return end subroutine getdis(i,sm,N,m,ns) implicit none integer*4 i,ns,sm(ns),N,m,mfc,j,jj,mr,ip,is,Na,js,jt integer*4,allocatable::um(:) c ip: number of placed excitations c is: when the nex series starts c mr: rest excitations that need to be placed c js:at which center we start actually allocate(um(N)) sm=0 ip=0 is=1 Na=N js=1 99 ip=ip+1 mr=m-ip if(mr.ge.0)then do 1 j=1,Na-mr um(j)=0 do 1 jj=1,j 1 um(j)=um(j)+mfc(Na-jj,mr) sm(ip)=js jt=1 do 2 j=2,Na-mr if(i.gt.is+um(j-1)-1.and.i.le.is+um(j)-1)then is=is+um(j-1) sm(ip)=js+j-1 jt=j endif 2 continue js=js+jt Na=N-js+1 endif if(mr.ge.1)goto 99 return end Subroutine FCA2NQ(NOsc,Idx,NQm,InCl,LNQ,NQ) Implicit Integer(A-H,O-Z) C C Franck-Condon storage Array to Number of Quanta C Get the vector of quantum numbers NQ from the index in the storage C array C C Input: C NOsc : Max. num. of excited oscillators (=class) C Idx : Index in the storage array C NQm : (*) Maximum number reachable by each quantum in NQ C C Output: C InCl : True if the state is in class NOsc C LNQF : Length of NQ actually used C NQ : (LNQ) Vector of the quantum numbers C C Dimensions Integer NOsc, LNQ C Input Integer NQm(*), Idx C Output Integer NQ(*) Logical InCl C Local Integer i, j, N, NQi, NVS C LNQ = 0 C C Case |0>: nothing to do If(Idx.le.1) then C -- FUNDAMENTAL STATE |0> -- else if(Idx.le.NQm(1)) then C -- ONLY THE 1ST MODE IS EXCITED -- NQ(1) = Idx - 1 LNQ = 1 else C -- GENERAL CASE -- NVS = NQm(1)*NQm(2) i = 2 100 If(NVS.lt.Idx) then i = i + 1 NVS = NVS * NQm(i) Goto 100 endIf N = Idx LNQ = i 200 NVS = NVS/NQm(i) NQi = Int(N/NVS) N = N - NQi*NVS If(N.eq.0) then C All prev modes (j') return end function nga(f) c number of atoms in gaussian output implicit none integer*4 nga,l,I,KA character*(*) f CHARACTER*21 z,s CHARACTER*50 FN CHARACTER*18 t data z/'Z-Matrix orientation:'/ data t/'Input orientation:'/ data s/'Standard orientation:'/ l=0 OPEN(2,FILE=f) 1 READ(2,2000,END=1000)FN 2000 FORMAT(A50) IF(FN(19:39).EQ.z.OR.FN(26:46).EQ.z.OR. 1 FN(20:37).EQ.t.OR.FN(27:44).EQ.t.OR. 1 FN(20:40).EQ.s.OR.FN(26:46).EQ.s)THEN DO 4 I=1,4 4 READ(2,*) 5 READ(2,2000)FN IF(FN(2:4).NE.'---')THEN l=l+1 BACKSPACE 2 READ(2,*)KA,KA IF(KA.EQ.-1)l=l-1 GOTO 5 ENDIF goto 1000 ENDIF goto 1 1000 close(2) nga=l END subroutine rdstat(NQ,N,ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi, 1sg,W,WMIN,WMAX,limag) implicit none integer*4 N,NQ real*8 ALPHA(N,3,3),ALPHAi(N,3,3),G(N,3,3),Gi(N,3,3), 1A(N,3,3,3),Ai(N,3,3,3),GGr(N,3,3),GGi(N,3,3),AAr(N,3,3,3), 1AAi(N,3,3,3),sg(N,N),W(N),WMIN,WMAX real*8,allocatable::lr_(:,:,:),li_(:,:,:),gr_(:,:,:),gi_(:,:,:), 1ar_(:,:,:,:),ai_(:,:,:,:) logical limag allocate(lr_(N,3,3),li_(N,3,3),gr_(N,3,3), 1gi_(N,3,3),ar_(N,3,3,3),ai_(N,3,3,3)) c read real and imaginary Cartesian tensro derivatives: CALL READTEN(N,lr_,gr_,ar_,0) li_=0.0d0 gi_=0.0d0 ai_=0.0d0 if(limag)CALL READTEN(N,li_,gi_,ai_,3) C transform them into normal modes CALL TRTEN(sg,N,lr_,ar_,gr_) CALL TRTEN(sg,N,li_,ai_,gi_) c write derivatives in normal modes: CALL WTQi(N,NQ,lr_,ar_,gr_,li_,ai_,gi_,lr_,W,WMIN,WMAX) c put them in generic ~ tensors: c Gaussians seems to give G in opposite sign: c should be c G~ = G - iG' c Re G~ = Im (G') c Im G~ = -Re (G') ALPHA = lr_ ALPHAi= li_ G = -gr_ Gi = -gi_ GGr = gr_ GGi = gi_ A = ar_ Ai = ai_ AAr = ar_ AAi = ai_ return end subroutine chkj(NQ,JT,block,W,WP) implicit none integer*4 NQ,block,i,f,j real*8 JT(NQ,NQ),sp,sf,W(*),WP(*),CM,an real*8,allocatable:: v(:,:) allocate(v(NQ,NQ)) v=JT CM=219474.63d0 write(6,602) 602 format('Checking linear dependence, rest norms:') do 1 i=1,NQ do 2 f=1,i-1 sp=0.0d0 sf=0.0d0 do 3 j=1,NQ sp=sp+v(i,j)*v(f,j) 3 sf=sf+v(f,j)*v(f,j) do 2 j=1,NQ 2 v(i,j)=v(i,j)-sp*v(f,j)/sf an=0.0d0 do 4 j=1,NQ 4 an=an+v(i,j)**2 an=dsqrt(an) if(i.le.block)then write(6,556)i,W(i)*CM,an,' blocked' 556 format(i6,f10.2,f10.4,A8,$) else write(6,556)i,W(i)*CM,an,' ' endif if(mod(i,2).eq.0)write(6,*) 1 continue write(6,*) do 5 i=block+1,NQ do 5 j=i+1,NQ 5 if(dabs(CM*(W( i)-W( j))).lt.0.001d0)write(6,601)'Ground ', 1i,j,W( i)*CM,W( j)*CM do 6 i=block+1,NQ do 6 j=i+1,NQ 6 if(dabs(CM*(WP(i)-WP(j))).lt.0.001d0)write(6,601)'Excited', 1i,j,WP(i)*CM,WP(j)*CM 601 format(1x,A7,' state, similar frequencies',2i4,2f12.2) return end subroutine wrram1(iwr,NA,ap,w,alr,ali,gtr,gti,gtcr,gtci,atr, 1ati,atcr,atci,luseA,luseG,ldo,roa,p) implicit none integer*4 ns0,ni0,NA,ap,ia,b,id,e,ii,is,iwr parameter (ns0=19,ni0=13) c ns0 : number of experimental setups c ni0 : number of invariants real*8 CM,AMU,BOHR,w,roa(2,ns0,NA),p, 1gpisvejc,EXCA,tr,ti,a(2), 1inv(ni0),co(ni0+2,ns0),aaar,aaai,gtaar,gtaai,gtcaar,gtcaai, 1alr(3,3),ali(3,3),gtr(3,3),gti(3,3),gtcr(3,3),gtci(3,3), 1atr(3,3,3),ati(3,3,3),atcr(3,3,3),atci(3,3,3), 1alsr(3,3),alsi(3,3),gtsr(3,3),gtsi(3,3),gtcsr(3,3),gtcsi(3,3), 1alar(3,3),alai(3,3),gtar(3,3),gtai(3,3),gtcar(3,3),gtcai(3,3), 1ear(3,3),eapr(3,3),easr(3,3),eapsr(3,3),eaar(3,3),eapar(3,3), 1eai(3,3),eapi(3,3),easi(3,3),eapsi(3,3),eaai(3,3),eapai(3,3), 1eacr(3,3),eapcr(3,3),eacsr(3,3),eapcsr(3,3),eacar(3,3), 1eapcar(3,3), 1eaci(3,3),eapci(3,3),eacsi(3,3),eapcsi(3,3),eacai(3,3),eapcai(3,3) logical lusea,luseg,ldo(ns0) c prefix RAM, prefix ROA c a2 bs(a)2 ba(a)2 aG bs(G)2 ba(G)2 bs(A)2 c ba(A)2 aG bs(G)2 ba(G)2 bs(A)2 ba(A)2 c script: ^ ^ ^ ^ ^ c c co: Raman prefactor, ROA prefactor,3 Raman invariants, c 10 ROA invariants, all for ns0=18 spectral kinds data co/ 4.0d0, 8.0d0, ! 1 1 2 ICP(0o): 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0,-45.0d0, 5.0d0, -5.0d0, -3.0d0, 1.0d0, 1 4.0d0, 2.0d0, ! 2 3 4 ICPx(90o): 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 8.0d0, 4.0d0, ! 2 3 4 ICPx(90o): 1 0.0d0, 3.0d0, 5.0d0, 0.0d0, 3.0d0, 5.0d0, -1.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 13.3333d0,6.666d0, ! 4 7 8 ICP*(90o): (magic) 1 9.0d0, 2.0d0, 2.0d0, 9.0d0, 2.0d0, 2.0d0, 0.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 4.0d0, 4.0d0, ! 5 9 10 ICPu(90o): 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 1.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 4.0d0, 8.0d0, ! 6 11 12 ICP(180o): 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, 7.0d0, 5.0d0, 1.0d0, 1 -1.0d0, 45.0d0, -5.0d0, 5.0d0, 3.0d0, -1.0d0, 1 4.0d0, 8.0d0, ! 7 13 14 SCP(0o): 1 45.0d0, 7.0d0, 5.0d0, 45.0d0, -5.0d0, 5.0d0, -3.0d0, 1 -1.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, 1 2.0d0, 4.0d0, ! 8 15 16 SCPx(90o): 1 45.0d0, 7.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, 1 4.0d0, 8.0d0, ! 9 17 18 SCPz(90o): 1 0.0d0, 3.0d0, 5.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0, 0.0d0, -3.0d0, -5.0d0, -1.0d0, -1.0d0, 1 13.3333d0,6.666d0, ! 10 19 20 SCP*(90o): 1 9.0d0, 2.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0, -9.0d0, -2.0d0, -2.0d0, 0.0d0, 0.0d0, 1 4.0d0, 4.0d0, ! 11 21 22 SCPu(90o): 1 45.0d0, 13.0d0, 15.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0, 1 4.0d0, 8.0d0, ! 12 23 24 SCP(180o): 1 45.0d0, 7.0d0, 5.0d0,-45.0d0, 5.0d0, -5.0d0, 3.0d0, 1 1.0d0,-45.0d0, -7.0d0, -5.0d0, 1.0d0, 1.0d0, 1 4.0d0, 8.0d0, ! 13 25 26 DCPI(0o): 1 45.0d0, 1.0d0, 5.0d0, 45.0d0, 1.0d0, 5.0d0, -1.0d0, 1 -1.0d0,-45.0d0, -1.0d0, -5.0d0, -1.0d0, +1.0d0, 1 2.0d0, 2.0d0, ! 14 27 28 DCPI(90o): 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 +1.0d0,-45.0d0,-13.0d0,-15.0d0, -1.0d0, -1.0d0, 1 24.0d0, 16.0d0, ! 15 29 30 DCPI(180o): 1 0.0d0, 1.0d0, 0.0d0, 0.0d0, 3.0d0, 0.0d0, 1.0d0, 1 0.0d0, 0.0d0, -3.0d0, 0.0d0, 1.0d0, 0.0d0, 1 24.0d0, 16.0d0, ! 16 31 32 DCPII(0o): 1 0.0d0, 1.0d0, 0.0d0, 0.0d0, 3.0d0, 0.0d0, 1.0d0, 1 0.0d0, 0.0d0, 3.0d0, 0.0d0, -1.0d0, 0.0d0, 1 2.0d0, 2.0d0, ! 17 33 34 DCPII(90o): 1 45.0d0, 13.0d0, 15.0d0, 45.0d0, 13.0d0, 15.0d0, -1.0d0, 1 1.0d0, 45.0d0, 13.0d0, 15.0d0, 1.0d0, 1.0d0, 1 4.0d0, 8.0d0, ! 18 35 36 DCPII(180o): 1 45.0d0, 1.0d0, 5.0d0, 45.0d0, 1.0d0, 5.0d0, -1.0d0, 1 -1.0d0, 45.0d0, 1.0d0, 5.0d0, 1.0d0, -1.0d0, 1 1.0d0, 0.0d0, ! 19 special 1 1.0d0, .0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 1 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/ c Prefix Raman ROA c a2 bs(a)2 ba(a)2 aG bs(G)2 ba(G)2 bs(A)2 c ba(A)2 aG bs(G)2 ba(G)2 bs(A)2 ba(A)2 c script script script script script c c invariants: c 1 a2 =(1/9)Re (as_aa as*_bb) c 2 bs(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb) c 3 ba(a)2 =(1/2)Re(3as_ab as*_ab-as_aa as*_bb) c 4 aG =(1/9)Im(as_aa Gs*_bb) c 5 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb) c 6 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab) c 7 bs(A)2 =(w/2)Im(i as_ab eas*_ab) c ea_ab=eps_adg A_d,gb c 8 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab c eap_ab=eps_abg A_d,gd c script tensors: c 9 aG =(1/9)Im(as_aa Gs*_bb) c 10 bs(G)2 =(1/2)Im(3as_ab Gs*_ab-as_aa Gs*_bb) c 11 ba(G)2 =(3/2)Im(3aa_ab Ga*_ab) c 12 bs(A)2 =(w/2)Im(i as_ab eas*_ab) c 13 ba(A)2 =(w/2)Im(i aa_ab eaa*_ab)+i aa_ab eap_ab CM=219474.63d0 AMU=1822.0d0 BOHR=0.529177d0 EXCA=(1.0d7/(w*CM))*10.0d0 gpisvejc=(AMU*BOHR**5)*1.0d4*2.0d0*4.0d0*atan(1.0d0)/EXCA c clight=137.03599d0 c eap_ab = eps_abc A_d,cd call calcep(eapr ,atr ) call calcep(eapi ,ati ) call calcep(eapci,atci) call calcep(eapcr,atcr) c c ea_ab = eps_adc A_d,cb ear =0.0d0 eai =0.0d0 eacr=0.0d0 eaci=0.0d0 do 2 ia=1,3 id=ia+1 if(id.gt.3)id=1 e=id+1 if(e.gt.3)e=1 do 2 b=1,3 ear( ia,b)=atr (id,e,b)-atr (e,id,b) eai( ia,b)=ati (id,e,b)-ati (e,id,b) eacr(ia,b)=atcr(id,e,b)-atcr(e,id,b) 2 eaci(ia,b)=atci(id,e,b)-atci(e,id,b) c symmetric and antisymmetric tensor combinations call dsa(alr , ali, alsr, alsi, alar, alai) call dsa(gtr , gti, gtsr, gtsi, gtar, gtai) call dsa(gtcr,gtci,gtcsr,gtcsi,gtcar,gtcai) call dsa(ear ,eai ,easr ,easi ,eaar ,eaai ) call dsa(eapr ,eapi ,eapsr ,eapsi ,eapar ,eapai ) call dsa(eacr ,eaci ,eacsr ,eacsi ,eacar ,eacai ) call dsa(eapcr,eapci,eapcsr,eapcsi,eapcar,eapcai) inv=0.0d0 c a2: (1/9) Re (als_aa als*_bb) aaar=alsr(1,1)+alsr(2,2)+alsr(3,3) aaai=alsi(1,1)+alsi(2,2)+alsi(3,3) inv(1)=(aaar*aaar+aaai*aaai)/9.0d0*(AMU*BOHR**4) c beta_s(alpha)2:(1/2)Re(3als_ab als*_ab-als_aa als*_bb) call abab(tr,ti,alsr,alsi,alsr,alsi) inv(2)=1.5d0*tr*(AMU*BOHR**4)-4.5d0*inv(1) c beta_a(alpha)2:(3/2)Re(als_ab als*_ab) call abab(tr,ti,alar,alai,alar,alai) inv(3)=1.5d0*tr*(AMU*BOHR**4) if(luseg)then c aG: (1/9) Im(als_aa gs*_bb) gtaar=gtsr(1,1)+gtsr(2,2)+gtsr(3,3) gtaai=gtsi(1,1)+gtsi(2,2)+gtsi(3,3) inv(4)=(aaai*gtaar-aaar*gtaai)/9.0d0*gpisvejc c beta_s(G)2:(1/2)Im(3als_ab Gs*_ab-als_aa Gs*_bb) call abab(tr,ti,alsr,alsi,gtsr,gtsi) inv(5)=1.5d0*ti*gpisvejc-4.5d0*inv(4) c beta_a(G)2: (3/2)Im(ala_ab Ga*_ab) call abab(tr,ti,alar,alai,gtar,gtai) inv(6)=1.5d0*ti*gpisvejc c aGc: gtcaar=gtcsr(1,1)+gtcsr(2,2)+gtcsr(3,3) gtcaai=gtcsi(1,1)+gtcsi(2,2)+gtcsi(3,3) inv(9)=(aaai*gtcaar-aaar*gtcaai)/9.0d0*gpisvejc c beta_s(Gc)2: call abab(tr,ti,alsr,alsi,gtcsr,gtcsi) inv(10)=1.5d0*ti*gpisvejc-4.5d0*inv(9) c beta_a(Gc)2: (3/2)Im(ala_ab Gca*_ab) call abab(tr,ti,alar,alai,gtcar,gtcai) inv(11)=1.5d0*ti*gpisvejc endif if(lusea)then c betasA2: (w/2)Im(i als_ab e_agd As*g,db) call abab(tr,ti,alsr,alsi,easr,easi) c inv(7)=0.5d0*w*tr/clight*gpisvejc c gaussian units: inv(7)=0.5d0*tr*gpisvejc c betaaA2: (w/2)Im(i[ala_ab[e_agd Aa*g,db+ e_abg Aa*d,gd) call abab(tr,ti,alar,alai,eaar,eaai) inv(8)=tr call abab(tr,ti,alar,alai,eapar,eapai) inv(8)=0.5d0*(inv(8)+tr)*gpisvejc c betasAc2: (w/2)Im(i als_ab e_agd Acs*g,db) call abab(tr,ti,alsr,alsi,eacsr,eacsi) inv(12)=0.5d0*tr*gpisvejc c betaaAc2: call abab(tr,ti,alar,alai,eacar,eacai) inv(13)=tr call abab(tr,ti,alar,alai,eapcar,eapcai) inv(13)=0.5d0*(inv(13)+tr)*gpisvejc endif c spectral intensities: do 3 is=1,ns0 if(ldo(is))then c Raman: a(1)= 1 co(1,is)*(co(3,is)*inv(1)+co(4,is)*inv(2)+co(5,is)*inv(3)) c ROA: a(2)=0.0d0 do 4 ii=4,ni0 4 a(2)=a(2)+inv(ii)*co(2+ii,is) a(2)=co(2,is)*a(2) roa(1,is,ap)=roa(1,is,ap)+a(1)*p roa(2,is,ap)=roa(2,is,ap)+a(2)*p if(iwr.gt.1)then write(6,600)p,is,a(1),a(2) 600 format(' p:',e10.2,i3,' ram:',e12.4,' roa:',e12.4) endif endif 3 continue return end c ==============================================================