PROGRAM NEW6D 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),f,IERR, 1jjs,iwr,nlim,block,nga c parameter (ns0=19,maxcen=10) parameter (ns0=19) c INTEGER*4 nc,fc,ncen(maxcen),fcen(maxcen),fexc(maxcen), c 1nexc(maxcen) real*8 u0(3),m0(3),v0(3),q0(6),alr(3,3),ali(3,3),gtr(3,3), 1gti(3,3),gcr(3,3),gci(3,3),atr(3,3,3),ati(3,3,3),uu(81), 1acr(3,3,3),aci(3,3,3),wexc,we,Gau,wmin,wmax,pis, 1aasr(3,3),aasi(3,3),gsr(3,3),gsi(3,3),gcsr(3,3),gcsi(3,3), 1asr(3,3,3),asi(3,3,3),acsr(3,3,3),acsi(3,3,3),iint(81), 1jint(81),av2,avn,f0,g0,f1,g1,CM,unorm,ddinorm,dtol,wz,wz0 real*8,allocatable::ALPHA(:,:,:),A(:,:,:,:),G(:,:,:), 1dd(:),vv(:),aa(:),qq(:),ddi(:),vvi(:),rri(:),qqi(:), 1Gi(:,:,:),ALPHAi(:,:,:),Ai(:,:,:,:),gradi(:), 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(:,:),ff(:,:),gg(:,:) LOGICAL asym(5),luseA,luseG,ldo(ns0),lstat,slist,lvert,eattt, 1lalt character*80,allocatable::sx(:) C WRITE(6,6666) 6666 FORMAT(/,/, A' * * * * * * * * * * * * * * * * * * * * * *',/, 1' * * ',/, 2' * RESONANCE RAMAN OPTICAL ACTIVITY * ',/, 3' * * ',/, 4' * * * * * * * * * * * * * * * * * * * * * *',/, 8' ',/, A'INPUT FILES: NEW6D.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,asym,pis, 1lstat,slist,iwr,lvert,eattt,nlim,dtol,lalt,block) 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 allocate(ALPHA(N,3,3),A(N,3,3,3),G(N,3,3),ALPHAi(N,3,3), 1Ai(N,3,3,3),Gi(N,3,3),AAi(N,3,3,3),GGi(N,3,3),AAr(N,3,3,3), 1GGr(N,3,3),grad(N),gradi(N),sg(N,N),se(N,N),Am(N,N),B(N),C(N,N), 1D(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)) c zero-out ~ tensors: call t0(N,ALPHA,ALPHAi,G,Gi,GGr,GGi,A,Ai,AAr,AAi) jj=jjs(sx(1)) call rduschs(N,NQ,sg,sx(1)(1:jj)//'.DUSCH','Ground S-Matrix') if(iwr.ne.0)call control(N,NQ,sg) call rduschw(N,W,sx(1)(1:jj)//'.DUSCH','Ground w') write(6,*)'frequencies read' if(lstat)then c read static polarizabilities: CALL READTEN(N,NAT,ALPHA,G,A,0) CALL READTEN(N,NAT,ALPHAi,Gi,Ai,3) C transform them into normal modes CALL TRTEN(sg,3*NAT,ALPHA,A,G) CALL TRTEN(sg,3*NAT,ALPHAi,Ai,Gi) c write derivatives in normal modes: CALL WTQi(N,NQ,ALPHA,A,G,ALPHAi,Ai,Gi,ALPHA,W,WMIN,WMAX) c put them in ~ tensors: call wtt(N,NQ,G,Gi,GGr,GGi,A,Ai,AAr,AAi) write(6,*)'Static tensors read' else write(6,*)'Static tensors skipped' endif write(6,608)asym 608 format(' alpha G GC A AC:',5l2) write(6,617)wexc*CM,Gau*CM 617 format(' wexc/cm-1:',f10.0,' G/cm-1:',f10.0,/, 1' state alpha norm we/cm-1 wexc/cm-1 gf0 gf1 D Dd: ') c cycle over files/electronic excited states: c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee do 1000 i=1,ni c read ground and excited state s-matrix: jj=jjs(sx(i)) call rduschs(N,NQ,sg,sx(i)(1:jj)//'.DUSCH','Ground S-Matrix') call rduschs(N,NQ,se,sx(i)(1:jj)//'.DUSCH','Excited S-Matri') write(6,609)i,N,NQ 609 format(i6,'. state: N, NQ:',2i6) c read the Dushchinsky matrices iw=0 ip=0 allocate(JT(NQ,NQ),K(NQ)) 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 wzz(NQ,N,W,WP,wz,wz0) write(6,607)wz*CM,wz0*CM 607 format(' ZPE ground, excited:',2f12.2,' cm-1') call chkj(NQ,JT,block,W,WP) allocate(J(NQ,NQ),jtj(NQ,NQ),jtji(NQ,NQ), 1U(NQ,NQ),V(NQ,NQ),X(NQ,NQ),Y(NQ,NQ),ki(NQ),sj(NQ,NQ), 1ff(81,NQ),gg(81,NQ)) 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') write(6,604)1,1,jtj(1,1),jtji(1,1),1,NQ,jtj(1,NQ),jtji(1,NQ), 1NQ,NQ,jtj(NQ,NQ),jtji(NQ,NQ) 604 format(' Inversion info, i j M Mi:',/, 13(2i4,2e12.4,/)) 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,we) c read moment derivatives: call rdd(sx(i),nat,NROOT,u0,v0,m0,q0,dd,vv,aa,qq,grad) c erase derivatives if required: if(i.gt.nlim)then write(6,*)' state ',i,' derivatives erased' dd=0.0d0 vv=0.0d0 aa=0.0d0 qq=0.0d0 endif c erase some atom derivatives if required: if(eattt)call ede(nat,dd,vv,aa,qq,grad) write(6,67)' Dipole ',u0 write(6,67)' Velocity ',v0 write(6,67)' Magnetic ',m0 write(6,67)' Quadrupole',q0 67 format(A11,6f12.6) 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:') call wrge('DEGS.TEN',nat,dd,aa,qq,vv,u0,m0,v0,q0) 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 uu,ff,gg: call sc(N,NQ,u0,m0,q0,ddi,rri,qqi,uu,ff,gg) unorm=dsqrt(u0(1)**2+u0(2)**2+u0(3)**2) ddinorm=0.0d0 do 201 f=1,3*NQ 201 ddinorm=ddinorm+ddi(f)**2 ddinorm=dsqrt(ddinorm) 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,ki,U,V,X,Y,W,WP,sj,wz,block) c |0>: c nc=0 c nexc=0 c ncen=0 c make freqeuncy functions: call mkf0(we,wexc,Gau,f0,g0,f1,g1,lalt) write(6,557)we,f0,g0,f1,g1 557 format(' we:',f12.6,/,' f0:',e12.4,' g0:',e12.4, 1' f1:',e12.4,' g1:',e12.4) c final vibrational states (fundamental modes only now): avn=0.0d0 do 1 f=block+1,NQ c fundamental only: call mkijf(N,NQ,f,W,uu,ff,gg,iint,jint,X,Y,U,V,ki,block) c construct i-> f transition polarizabilities: call ctp(we,wexc,iint,jint,asym,aasr,aasi,gsr,gsi, 1gcsr,gcsi,asr,asi,acsr,acsi,f0,g0,f1,g1) if(f.eq.block+1)then write(6,*)'control list for the first mode:' write(6,*)'assr:' write(6,699)aasr write(6,*)'I:' call wnn(9,iint,'norm:') write(6,699)(iint(jj),jj=1,9) write(6,*)'J:' call wnn(9,jint,'norm:') write(6,699)(jint(jj),jj=1,9) 699 format(3e13.3) endif c add them to ~ tensor f-derivatives by pis factor call adss(f,N,W,wexc,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1gcsr,gcsi,GGr,GGi,asr,asi,A,Ai,acsr,acsi,AAr,AAi,pis,av2,iwr) 1 avn=avn+av2 write(6,606)i,avn/dble(NQ),we*CM,dsqrt(g0**2+f0**2), 1dsqrt(g1**2+f1**2),unorm,ddinorm 606 format(' ##',i6,e12.4,f10.0,4e12.4) 1000 deallocate(J,JT,K,jtj,jtji,U,V,X,Y,ki,sj,ff,gg) c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee call initab(ldo) c do 2 i=block+1,NQ c rewrite into to smaller arrays: call trsten(N,i,ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi,alr, 1ali,gtr,gti,gcr,gci,atr,ati,acr,aci) c get intensities: 2 call wrram3(wexc,0.0d0,W(i),i,alr,ali,gtr,gti,gcr,gci,atr, 1ati,acr,aci,luseA,luseG,ldo) call closetab(ldo) 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 wrram3(w,ei,ef,nt,alr,ali,gtr,gti,gtcr,gtci,atr, 1ati,atcr,atci,luseA,luseG,ldo) implicit none integer*4 nt,ns0,ni0, 1ia,b,id,e,ii,is parameter (ns0=19,ni0=13) c ns0 : number of experimental setups c ni0 : number of invariants real*8 ef,ei,CM,AMU,BOHR,ECM,w, 1YDY,YDX,gpisvejc,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 ECM=(ef-ei)*CM 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) YDX=0.0d0 YDY=0.0d0 WRITE(50+is,3001)nt,ECM,YDX,YDY,a(1),a(2) 3001 FORMAT(I7,f9.2,2f3.0,g12.4,' 0 0 0 0',2g12.4) endif 3 continue WRITE(44,4305)nt,ECM,inv 4305 FORMAT(i6,f10.2,13e11.4) 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, 1ntr,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) 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 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:') 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,NAT,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.gt.N/3)call report('too many 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,*,err=9)(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,*)NAT,' atoms, tensors read in from FILE.TTT' RETURN 9 WRITE(6,*)'Could not read tensors, set to zero' P=0.0d0 G=0.0d0 A=0.0d0 RETURN END C ============================================================== SUBROUTINE IOPAR(w,wmin,wmax,luseA,luseG,ldo,Gau,asym, 1pis,lstat,slist,iwr,lvert,eattt,nlim,dtol,lalt, 1block) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) parameter (ns0=19) integer*4 block LOGICAL lex,luseA,luseG,ldo(ns0),asym(5),lstat,slist,lvert,eattt, 1lalt 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 number of frozen modes block=6 c alternative form of complex frequency function lalt=.false. c maximum excitation with derivatives considered: nlim=1000 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 asym=.true. Gammacm=500.0d0 ldo=.false. wmin=0.0d0 wmax=1.0d90 c excitation light wavelength in nm: EXCNM=532.0d0 luseA=.true. c use G'-tensor derivatives: luseG=.true. INQUIRE(FILE='NEW6D.OPT',exist=lex) if(lex)then OPEN(2,FILE='NEW6D.OPT') WRITE(6,*)' FILE NEW6D.OPT FOUND ...' else goto 99 endif 2001 READ(2,1111)STRING 1111 FORMAT(A15) IF(STRING(1:4).EQ.'BLOC')READ(2,*)block IF(STRING(1:4).EQ.'LALT')READ(2,*)lalt IF(STRING(1:4).EQ.'DTOL')READ(2,*)dtol IF(STRING(1:4).EQ.'NLIM')READ(2,*)nlim 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.'ASYM')READ(2,*)ias,asym(ias) 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 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 6009 format(' l_exc ',f5.1,' nm, Gamma:',f10.1,' cm-1') w=1.0d7/EXCNM/CM Gau=Gammacm/CM 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 adss(f,N,W,wexc,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1gcsr,gcsi,Gcr,Gci,asr,asi,A,Ai,acsr,acsi,Acr,Aci,pis,av2,iwr) implicit none integer*4 N,i,j,f,l,iwr real*8 aasr(3,3),aasi(3,3),ALPHA(N,3,3),W(N),wexc,clight, 1ALPHAi(N,3,3),gsr(3,3),gsi(3,3),gcsr(3,3),gcsi(3,3),h,o, 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),pis,av1,av2 clight=137.03599d0 h=1.0d0/dsqrt(2.0d0*w(f)) o=h*clight/wexc av1=h*(ALPHA(f,1,1)+ALPHA(f,2,2)+ALPHA(f,3,3)) av2= (aasr( 1,1)+aasr( 2,2)+ aasr( 3,3)) do 1 i=1,3 do 1 j=1,3 ALPHA( f,i,j)=h* ALPHA( f,i,j)+pis*aasr(i,j) ALPHAi(f,i,j)=h* ALPHAi(f,i,j)+pis*aasi(i,j) G( f,i,j)=h*wexc*G( f,i,j)+pis*gsr( i,j) Gi( f,i,j)=h*wexc*Gi( f,i,j)+pis*gsi( i,j) Gcr( f,i,j)=h*wexc*Gcr( f,i,j)+pis*gcsr(i,j) Gci( f,i,j)=h*wexc*Gci( f,i,j)+pis*gcsi(i,j) do 1 l=1,3 A( f,i,j,l)=o*A( f,i,j,l)+pis*asr( i,j,l) Ai( f,i,j,l)=o*Ai( f,i,j,l)+pis*asi( i,j,l) Acr( f,i,j,l)=o*Acr( f,i,j,l)+pis*acsr(i,j,l) 1 Aci( f,i,j,l)=o*Aci( f,i,j,l)+pis*acsi(i,j,l) if(iwr.ge.2)write(6,600)f,av1,av2,pis 600 format(i3,' stat',e12.4,': vib',e12.4,' pis=',e10.2) 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 close(44) return end c ============================================================ subroutine initab(ldo) implicit none integer*4 i,j,ns0 parameter(ns0=19) character*11 st(ns0) character*6 iss(13) 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 '/ data iss/' a2', 'bs(a)2', 'ba(a)2', ' aG', 'bs(G)2', 1 'ba(G)2', 'bs(A)2', 'ba(A)2', ' aG', 'bs(G)2', 2 'ba(G)2', 'bs(A)2', 'ba(A)2'/ 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,*) OPEN(44,FILE='INV.TXT') write(44,4000) 4000 format(' n w',$) do 4 i=1,13 4 write(44,4001)iss(i) 4001 format(5x,a6,$) write(44,*) 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 trsten(N,i,ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi,alr, 1ali,gtr,gti,gcr,gci,atr,ati,acr,aci) implicit none integer*4 N,d,b,c,i real*8 alr(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),ALPHA(N,3,3), 1A(N,3,3,3),G(N,3,3),ALPHAi(N,3,3),Ai(N,3,3,3),Gi(N,3,3), 1AAr(N,3,3,3),GGr(N,3,3),AAi(N,3,3,3),GGi(N,3,3) do 2 d=1,3 do 2 b=1,3 alr(d,b)= ALPHA( i,d,b) ali(d,b)= ALPHAi(i,d,b) gtr(d,b)= G ( i,d,b) gti(d,b)= Gi ( i,d,b) gcr(d,b)= GGr( i,d,b) gci(d,b)= GGi( i,d,b) do 2 c=1,3 atr(d,b,c)= A( i,d,b,c) ati(d,b,c)= Ai( i,d,b,c) acr(d,b,c)= AAr(i,d,b,c) 2 aci(d,b,c)= AAi(i,d,b,c) 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.3)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 wtt(N,NQ,G,Gi,Gcr,Gci,Ar,Ai,Acr,Aci) c fill script tensors with FFR ones implicit none integer*4 a,b,c,N,NQ,i real*8 G(N,3,3),Gi(N,3,3),Gcr(N,3,3),Gci(N,3,3), 1Ar(N,3,3,3),Ai(N,3,3,3),Acr(N,3,3,3),Aci(N,3,3,3) 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') G=-G Gi=-Gi do 1 i=1,NQ do 1 a=1,3 do 1 b=1,3 Gcr(i,a,b)=-G( i,b,a) Gci(i,a,b)=-Gi(i,b,a) do 1 c=1,3 Acr(i,a,b,c)=Ar(i,a,b,c) 1 Aci(i,a,b,c)=Ai(i,a,b,c) 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 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) 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 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 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(we,w,Gau,f0,g0,f1,g1,lalt) implicit none real*8 we,w,Gau,f0,g0,f1,g1 logical lalt f0=(we-w)*(we+w) / (((we-w)*(we+w))**2+(2.0d0*we*Gau)**2) g0=-2.0d0*we*Gau / (((we-w)*(we+w))**2+(2.0d0*we*Gau)**2) if(lalt)then f1=1.0d0/((we-w)**2+Gau**2) g1=0.0d0 else f1=((we-w)**2-Gau**2) / ((we-w)**2+Gau**2)**2 g1=-2.0d0*(we-w)*Gau / ((we-w)**2+Gau**2)**2 endif return end subroutine ctp(we,w,iint,jint,asym,aasr,aasi,gsr,gsi, 1gcsr,gcsi,asr,asi,acsr,acsi,f0,g0,f1,g1) implicit none integer*4 i,j,a,b,c real*8 we,w,f0,g0,f1,g1,iint(81),jint(81), 1aasr(3,3),aasi(3,3),gsr(3,3),gsi(3,3),gcsr(3,3), 1gcsi(3,3),asr(3,3,3),asi(3,3,3),acsr(3,3,3),acsi(3,3,3) logical asym(5) aasr=0.0d0 aasi=0.0d0 gsr=0.0d0 gsi=0.0d0 gcsr=0.0d0 gcsi=0.0d0 asr=0.0d0 asi=0.0d0 acsr=0.0d0 acsi=0.0d0 c make asymmetric alpha i=0 j=0 do 4 a=1,3 do 4 b=1,3 i=i+1 if(asym(1))then aasr(a,b)=2.0d0*we*iint(i)*f0-jint(i)*f1 aasi(a,b)=2.0d0*we*iint(i)*g0-jint(i)*g1 endif if(asym(2))then c Re = - Im * Im gsr(a,b)=-2.0d0*w*iint(i+9)*g0+jint(i+9)*g1 c Im = Im * Re gsi(a,b)= 2.0d0*w*iint(i+9)*f0-jint(i+9)*f1 endif if(asym(3))then gcsr(a,b)=-2.0d0*w*iint(i+18)*g0+jint(i+18)*g1 gcsi(a,b)= 2.0d0*w*iint(i+18)*f0-jint(i+18)*f1 endif do 4 c=1,3 j=j+1 if(asym(4))then asr(a,b,c)=2.0d0*we*iint(j+27)*f0-jint(j+27)*f1 asi(a,b,c)=2.0d0*we*iint(j+27)*g0-jint(j+27)*g1 endif if(asym(4))then acsr(a,b,c)=2.0d0*we*iint(j+54)*f0-jint(j+54)*f1 acsi(a,b,c)=2.0d0*we*iint(j+54)*g0-jint(j+54)*g1 endif 4 continue return end subroutine mkijf(N,NQ,f,W,uu,ff,gg,iint,jint,X,Y,U,V,ki,b) c fundamental transitions only 0->1f implicit none integer*4 N,NQ,f,ii,j,b real*8 W(N),U(NQ,NQ),V(NQ,NQ),X(NQ,NQ),Y(NQ,NQ),ki(NQ), 1uu(81),ff(81,NQ),gg(81,NQ),iint(81),jint(81),s do 1 ii=1,81 1 iint(ii)=ff(ii,f)/dsqrt(2.0d0*W(f)) s=0.0d0 do 2 j=b+1,NQ 2 s=s-2.0d0*V(j,f)*ki(j) do 21 ii=1,81 21 jint(ii)=uu(ii)*s s=0.0d0 do 3 j=b+1,NQ 3 s=s+Y(j,j) do 31 ii=1,81 31 jint(ii)=jint(ii)+iint(ii)*s do 4 j=b+1,NQ do 4 ii=1,81 4 jint(ii)=jint(ii)+ff(ii,j)*dsqrt(2.0d0/W(j))*X(j,f) do 5 j=b+1,NQ do 5 ii=1,81 5 jint(ii)=jint(ii)+gg(ii,j)*dsqrt(2.0d0/W(j))*2.0d0*U(f,j) 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 subroutine wzz(NQ,N,W,WP,wz,wz0) integer*4 NQ,N,i real*8 W(N),WP(N),wz,wz0 wz=0.0d0 wz0=0.0d0 do 1 i=1,NQ wz0=wz0+0.5d0*W(i) 1 wz=wz+0.5d0*WP(i) 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=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=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 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