PROGRAM NEW6E c experimental - asymmetric polarizabilities from transition dipole c derivatives IMPLICIT none INTEGER*4 ns0,N,NAT,NQ,NQE,NROOT,i parameter (ns0=19) 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), 1acr(3,3,3),aci(3,3,3),w,we,Gau,wmin,wmax,pis,c,fsk,qmin,qmax real*8,allocatable::S(:,:),ALPHA(:,:,:),A(:,:,:,:),G(:,:,:),E(:), 1dd(:),vv(:),aa(:),qq(:),ddi(:),vvi(:),rri(:),qqi(:), 1Gi(:,:,:),ALPHAi(:,:,:),Ai(:,:,:,:),aasr(:,:,:),aasi(:,:,:), 1gsr(:,:,:),gsi(:,:,:),gcsr(:,:,:),gcsi(:,:,:),gradi(:), 1asr(:,:,:,:),asi(:,:,:,:),acsr(:,:,:,:),acsi(:,:,:,:), 1GGi(:,:,:),AAi(:,:,:,:),GGr(:,:,:),AAr(:,:,:,:),grad(:),K(:), 1se(:,:),ee(:),KP(:),x(:) LOGICAL asym(5),luseA,luseG,ldo(ns0),ltd C WRITE(6,6666) 6666 FORMAT(/,/, A' * * * * * * * * * * * * * * * * * * * * * *',/, 1' * * ',/, 2' * RAMAN OPTICAL ACTIVITY * ',/, 3' * * ',/, 4' * * * * * * * * * * * * * * * * * * * * * *',/, 8' ',/, 8' Petr Bour 1997-2023 ',/,/,/, A'INPUT FILES: NEW6E.OPT, FILE.TTT, F.INP, TD.OUT',/, B'OUTPUT : ROA.TAB',/) CALL IOPAR(w,wmin,wmax,luseA,luseG,NROOT,ldo,Gau,ltd,asym,pis,c, 1fsk,qmin,qmax) c read transformed tensors from FILE.Q.TTT: OPEN(2,FILE='FILE.TTT',STATUS='OLD') READ(2,*) READ(2,*)NAT close(2) N=3*NAT allocate(S(N,N),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),E(N),aasr(N,3,3),aasi(N,3,3),x(N), 1AAi(N,3,3,3),GGi(N,3,3),AAr(N,3,3,3),GGr(N,3,3),K(N),KP(N), 1gsr(N,3,3),gsi(N,3,3),gcsr(N,3,3),gcsi(N,3,3),grad(N),gradi(N), 1asr(N,3,3,3),asi(N,3,3,3),acsr(N,3,3,3),acsi(N,3,3,3), 1se(N,N),ee(N)) CALL READTEN(N,NAT,ALPHA,G,A,0) CALL READTEN(N,NAT,ALPHAi,Gi,Ai,3) c read S-matrix in atomic units: CALL READS(N,NQ,S,E,'F.INP') C transform the tensors into normal modes CALL TRTEN(S,3*NAT,ALPHA,A,G) CALL TRTEN(S,3*NAT,ALPHAi,Ai,Gi) WRITE(6,*)'(Complex part)' c write derivatives in normal modes: CALL WTQi(N,NQ,ALPHA,A,G,ALPHAi,Ai,Gi,ALPHA,E,WMIN,WMAX) c fill script tensors with FFR ones: call wtt(N,NQ,G,Gi,GGr,GGi,A,Ai,AAr,AAi) if(ltd)then allocate(dd(3*N),aa(3*N),vv(3*N),qq(6*N), 1 ddi(3*N),rri(3*N),vvi(3*N),qqi(6*N)) c check that we have the desired state: call chkroot('TD.OUT',NROOT) c read moments: call r00('TD.OUT',NROOT,u0,m0,v0,q0,we) c read moment derivatives: call rdd('TD.OUT',nat,NROOT,u0,v0,m0,q0,dd,vv,aa,qq,grad) C transform into normal modes call trafN(3,dd,ddi,nat,NQ,s) call trafN(3,aa,rri,nat,NQ,s) call trafN(6,qq,qqi,nat,NQ,s) call trafN(3,vv,vvi,nat,NQ,s) c excited state normal modes c read S-matrix for the excited state in atomic units: CALL READS(N,NQE,se,ee,'e0/F.INP') call trafN(1,grad,gradi,nat,NQE,se) c excited normal mode shifts: call wmse(N,NQE,gradi,ee,we,KP,qmin,qmax,fsk,se) c mode shifts in ground normal modes: call gms(N,NQ,NQE,s,se,K,KP) c ground mode shifts: call wms1(N,NQ,E,we,K,s,x) deallocate(dd,aa,vv,qq,grad) c make asymmetric alpha 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 write(6,*)asym if(asym(1))call mkaa( NQ,N,u0,ddi,aasr,aasi,Gau,we,w,E,gradi) if(asym(2))call mkaG( NQ,N,u0,m0,ddi,rri, gsr, gsi,Gau,we,w,E, 1 gradi) if(asym(3))call mkaGC(NQ,N,u0,m0,ddi,rri,gcsr,gcsi,Gau,we,w,E, 1 gradi) if(asym(4))call mka( NQ,N,u0,q0,ddi,qqi, asr, asi,Gau,we,w,E, 1 gradi) if(asym(5))call mkac( NQ,N,u0,q0,ddi,qqi,acsr,acsi,Gau,we,w,E, 1 gradi) c add it to the FFR parts call adss(N,NQ,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1 gcsr,gcsi,GGr,GGi,asr,asi,A,Ai,acsr,acsi,AAr,AAi,pis) endif call initab(ldo) c do 1 i=1,NQ call trsten(N,i,ALPHA,ALPHAi,G,Gi,A,Ai,GGr,GGi,AAr,AAi,alr, 1ali,gtr,gti,gcr,gci,atr,ati,acr,aci) 1 call wrram3(w,0.0d0,E(i),i,alr,ali,gtr,gti,gcr,gci,atr, 1ati,acr,aci,luseA,luseG,ldo) call closetab(ldo) if(.not.luseA)write(6,*)'A tensor deleted' if(.not.luseG)write(6,*)'G tensor deleted' WRITE(6,*)' PROGRAM FINISHED OK' 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,13f11.3) 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) write(6,*)'Analytical transition derivatives read' 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,*)'E' c dipole derivatives from gaussian: open(40,file='DEG.TEN') 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,*)'DEG.TEN written' return end C ========================================== SUBROUTINE READX(N,X,Z) IMPLICIT none integer*4 N,i,ix,Z(N/3) real*8 X(N) OPEN(4,FILE='F.INP',STATUS='OLD') read(4,*) do 1 i=1,N/3 1 read(4,*)Z(i),(X(ix+3*(i-1)),ix=1,3) close(4) return end C ========================================== SUBROUTINE READS(N3,nint,S,E,f) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION S(N3,N3),E(*) character*(*) f CM=219470.0d0 SFAC=0.0234280d0 C 1 a.u.= 2.1947E5 cm^-1 OPEN(4,FILE=f,STATUS='OLD') read(4,*)nint,nat3,nat N=NAT3 do 1 i=1,NAT 1 read(4,*) read(4,*) DO 2 I=1,NAT DO 2 J=1,NINT 2 read(4,*)(s(3*(i-1)+ix,j),ix=1,2), 1(s(3*(i-1)+ix,j),ix=1,3) read(4,*)nint READ(4,4000)(E(I),I=1,NINT) 4000 FORMAT(6F11.6) CLOSE(4) DO 3 I=1,NINT 3 E(I)=E(I)/CM DO 5 I=1,N DO 5 J=1,N 5 S(I,J)=S(I,J)*SFAC NAT=N/3 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,NROOT,ldo,Gau,ltd,asym, 1pis,c,fsk,qmin,qmax) IMPLICIT none INTEGER*4 NROOT,ns0,i,ias,ie REAL*8 w,wmin,wmax,Gau,pis,c,fsk,CM,Gammacm,EXCNM,qmin,qmax parameter (ns0=19) LOGICAL lex,luseA,luseG,ldo(ns0),ltd,asym(5) 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 c minimum frequency of involved modes (cm-1): qmin=50.0d0 c maximum frequency of involved modes (cm-1): qmax=9999.0d0 c mode shift scaling factor: fsk=1.0d0 pis=1.0d0 c=0.7d0 CM=219470.0d0 asym=.true. ltd=.true. Gammacm=500.0d0 ldo=.false. NROOT=1 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='NEW6E.OPT',exist=lex) if(lex)then OPEN(2,FILE='NEW6E.OPT') WRITE(6,*)' FILE NEW6E.OPT FOUND ...' else goto 99 endif 2001 READ(2,1111)STRING 1111 FORMAT(A15) IF(STRING(1:4).EQ.'CFAC')READ(2,*)c IF(STRING(1:3).EQ.'FSK') READ(2,*)fsk IF(STRING(1:4).EQ.'QMIN')READ(2,*)qmin IF(STRING(1:4).EQ.'QMAX')READ(2,*)qmax IF(STRING(1:4).EQ.'TDDF')READ(2,*)ltd 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.'ROOT')READ(2,*)NROOT IF(STRING(1:4).EQ.'ASYM')READ(2,*)ias,asym(ias) IF(STRING(1:4).EQ.'PISV')READ(2,*)pis IF(STRING(1:5).EQ.'NROOT')READ(2,*)NROOT 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 mkaG(NQ,N,u0,m0,ddi,aai,gr,gi,G,we,w,e,gradi) c ua mb implicit none integer*4 NQ,N,iq,a,b real*8 m0(3),ddi(3*N),gr(N,3,3),G,gi(N,3,3),we,w, 1fr,fi,aai(3*N),as,u0(3),a1,a2,a3,a4,u(3),e(*),gradi(*) c du0a/dQ<1|Q|0> m0b c G = ----------------- --> ------------------ (we - w - iG) c 2 2 c we - w + iG (we - w) + G fr= (we-w)/((we-w)**2+G**2) fi= -G/((we-w)**2+G**2) do 1 iq=1,NQ call seta(a1,a2,a3,a4,e(iq),gradi(iq)) call setu(u,ddi,iq) do 1 a=1,3 do 1 b=1,3 c supplied is m0 = <0|m|e> = - as=u0(a)*(-m0(b) )*a1+u(a)*(-m0(b) )*a2 1+ u0(a)*(-aai(b+3*(iq-1)))*a3+u(a)*(-aai(b+3*(iq-1)))*a4 gi(iq,a,b)= fr*as 1 gr(iq,a,b)= -fi*as return end subroutine mkaGC(NQ,N,u0,m0,ddi,aai,gr,gi,G,we,w,e,gradi) c ma ub implicit none integer*4 NQ,N,iq,a,b real*8 u0(3),aai(3*N),gr(N,3,3),G,gi(N,3,3),we,w, 1fr,fi,ddi(3*N),as,m0(3),a1,a2,a3,a4,u(3),e(*),gradi(*) c dm0a/dQ<1|Q|0> u0b c GC = ----------------- --> ------------------ (we - w - iG) c 2 2 c we - w + iG (we - w) + G fr= (we-w)/((we-w)**2+G**2) fi= -G/((we-w)**2+G**2) do 1 iq=1,NQ call seta(a1,a2,a3,a4,e(iq),gradi(iq)) call setu(u,ddi,iq) do 1 a=1,3 do 1 b=1,3 c supplied is dm/dQ = d<0|m|e>/dQ as=m0(a)* u0(b)*a1+aai(a+3*(iq-1))*u0(b)*a2 1+ m0(a)* u (b)*a2+aai(a+3*(iq-1))*u (b)*a4 gi(iq,a,b)= fr*as 1 gr(iq,a,b)= fi*as return end subroutine mkac(NQ,N,u0,q0,ddi,qqi,ar,ai,G,we,w,e,gradi) c qbc ua implicit none integer*4 NQ,N,iq,a,b,c real*8 u0(3),qqi(6*N),ar(N,3,3,3),G,ai(N,3,3,3),we,w,q(3,3), 1fr,fi,qd(3,3),as,ddi(3*N),q0(6),a1,a2,a3,a4,u(3),e(*),gradi(*) c dq0bc/dQ<1|Q|0> u0a c A = ----------------- --> ------------------- (we - w - iG) c 2 2 c we - w + iG (we - w) + G 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) do 1 iq=1,NQ call seta(a1,a2,a3,a4,e(iq),gradi(iq)) call setu(u,ddi,iq) 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) fr= (we-w)/((we-w)**2+G**2) fi= -G/((we-w)**2+G**2) do 1 a=1,3 do 1 b=1,3 do 1 c=1,3 as=q(b,c)* u0(a)*a1+qd(b,c)*u0(a)*a2 1+ q(b,c)* u (a)*a2+qd(b,c)*u (a)*a4 ar(iq,a,b,c)= fr*as 1 ai(iq,a,b,c)= fi*as return end subroutine mka(NQ,N,u0,q0,ddi,qqi,ar,ai,G,we,w,e,gradi) c ua qbc implicit none integer*4 NQ,N,iq,a,b,c real*8 q0(6),ddi(3*N),ar(N,3,3,3),G,ai(N,3,3,3),we,w,q(3,3), 1fr,fi,as,qqi(6*N),qd(3,3),u0(3),a1,a2,a3,a4,u(3),e(*),gradi(*) 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) fr= (we-w)/((we-w)**2+G**2) fi= -G/((we-w)**2+G**2) do 1 iq=1,NQ call seta(a1,a2,a3,a4,e(iq),gradi(iq)) call setu(u,ddi,iq) 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) do 1 a=1,3 do 1 b=1,3 do 1 c=1,3 as=u0(a)* q (b,c)*a1+u(a)*q (b,c)*a2 1+ u0(a)* qd(b,c)*a2+u(a)*qd(b,c)*a4 ar(iq,a,b,c)= fr*as 1 ai(iq,a,b,c)= fi*as return end subroutine mkaa(NQ,N,u0,ddi,aasr,aasi,G,we,w,e,gradi) implicit none integer*4 NQ,N,iq,a,b real*8 u0(3),ddi(3*N),aasr(N,3,3),G,aasi(N,3,3),we,w,as, 1fr,fi,gradi(N),e(N),a1,a2,a3,a4,u(3) c du0a du0b c ---- u0b(1-c)+u0a ----(1+c) c a b ddQ dQ c a = ----------------- -> ---------------------------(we-w-iG)<1|Q|0> c 2 2 c we - w + iG 2[(we - w) + G ] fr= (we-w)/((we-w)**2+G**2) fi= -G/((we-w)**2+G**2) do 1 iq=1,NQ call seta(a1,a2,a3,a4,e(iq),gradi(iq)) call setu(u,ddi,iq) do 1 a=1,3 do 1 b=1,3 as=u0(a)* u0(b)*a1+u(a)*u0(b)*a2 1+ u0(a)* u (b)*a2+u(a)*u (b)*a4 aasr(iq,a,b)=fr*as 1 aasi(iq,a,b)=fi*as write(6,*)'we',we*219470.0d0 write(6,*)'w ',w *219470.0d0 write(6,*)'G ',G *219470.0d0 write(6,*)'fr fi ',fr,fi 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 seta(a1,a2,a3,a4,w,g) implicit none real*8 q0,y0,a1,a2,a3,a4,ff,w,g c 2 2 c w 2 dE 2 d E 2 c E = -- Q' g'= -- = w Q', s = --- = w c 2 dQ' dQ'^2 c c g' = g(0') + s K', K' = g'/s c Q',K' ... excited electronic state modes c Q, K ... ground electronic state modes c K = -K', Q = -Q', g'=dE/dQ'=-g=-dE/dQ, c g is actually supplied as force = -g c K=q0: q0=-g/w**2 y0=dsqrt(w)*q0 a1=y0**3/2.0d0/dsqrt(2.0d0) a2=(y0**4+y0**2+4.0d0)/4.0d0/dsqrt(2.0d0*w) a3=y0*(y0**4-2.0d0*y0**2+4.0d0)/4.0d0/dsqrt(2.0d0*w) a4=y0*(y0**4-y0**2+6.0d0)/8.0d0/dsqrt(2.0d0*w) ff=exp(-y0**2/4.0d0) a1=a1*ff a2=a2*ff a3=a3*ff a4=a4*ff return end subroutine adss(N,NQ,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1gcsr,gcsi,Gcr,Gci,asr,asi,A,Ai,acsr,acsi,Acr,Aci,pis) implicit none integer*4 N,i,j,k,l,NQ real*8 aasr(N,3,3),aasi(N,3,3),ALPHA(N,3,3), 1ALPHAi(N,3,3),gsr(N,3,3),gsi(N,3,3),gcsr(N,3,3),gcsi(N,3,3), 1G(N,3,3),Gi(N,3,3),asr(N,3,3,3),asi(N,3,3,3),acsr(N,3,3,3), 1Gcr(N,3,3),Gci(N,3,3),acr(N,3,3,3),aci(N,3,3,3), 1acsi(N,3,3,3),A(N,3,3,3),Ai(N,3,3,3),pis do 1 k=1,NQ do 1 i=1,3 do 1 j=1,3 ALPHA( k,i,j)=ALPHA( k,i,j)+pis*aasr(k,i,j) ALPHAi(k,i,j)=ALPHAi(k,i,j)+pis*aasi(k,i,j) G( k,i,j)=G( k,i,j)+pis*gsr( k,i,j) Gi( k,i,j)=Gi( k,i,j)+pis*gsi( k,i,j) Gcr( k,i,j)=Gcr( k,i,j)+pis*gcsr(k,i,j) Gci( k,i,j)=Gci( k,i,j)+pis*gcsi(k,i,j) do 1 l=1,3 A( k,i,j,l)=A( k,i,j,l)+pis*asr( k,i,j,l) Ai( k,i,j,l)=Ai( k,i,j,l)+pis*asi( k,i,j,l) Acr( k,i,j,l)=Acr( k,i,j,l)+pis*acsr(k,i,j,l) 1 Aci( k,i,j,l)=Aci( k,i,j,l)+pis*acsi(k,i,j,l) 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 do 2 j=len(st),1,-1 2 if(st(i)(j:j).ne.' ')goto 3 3 open(50+i,file=st(i)(1:j)//'.tab') write(6,600)st(i) 600 format(a11) write(50+i,3000)st(i) 3000 format(20x,a11,/,' n wavelength (cm) Raman ROA',/,80(1h-)) endif 1 continue 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 chkroot(f,NROOT) implicit none integer*4 NROOT,found,nd,i character*(*) f character*80 s c 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)nd,found,NROOT 600 format(i6,' excited states',/, 1i6,' root found',/ 1i6,' root required') if(found.ne.NROOT)call report('Not consistent') 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 wmse(N,NQ,gradi,E,we,K,qmin,qmax,fsk,se) implicit none integer*4 N,NQ,i,j real*8 E(N),gradi(N),CM,we,fsk,wep,K(N),qmin,qmax, 1ecm,se(N,N) CM=219474.63d0 write(6,600) 600 format(' Excited mode frequencies and shifts') wep=we j=0 do 1 i=1,NQ ecm=E(i)*CM K(i)=0.0d0 if(ecm.gt.qmin.and.ecm.lt.qmax)then j=j+1 c excited energy = E0 + w^2/2 Q^2, grad = w^2 Q -> Q = grad/w^2, desired c shift K = -Q = -grad/w^2, in gradi the force = -grad is actually c supplied, so that K = gradi/w^2 K(i)=gradi(i)/E(i)**2*fsk wep=wep-(E(i)*K(i))**2/2.0d0 write(6,601)ecm,K(i) 601 format(2f10.2,$) if(mod(j,4).eq.0)write(6,*) endif 1 continue write(6,602)we*CM,wep*CM 602 format(/,' Vertical and interpolated energy (cm-1):',2f12.2) write(6,603)j,qmin,qmax,fsk 603 format(i6,' modes within',f10.1,' -',f10.1,' cm-1 involved,',/, 1' shift scaled by',f10.4) call shgeo(N,NQ,se,K,'FIP.X') return end subroutine wms1(N,NQ,E,we,K,s,x) implicit none integer*4 N,NQ,i real*8 E(N),CM,we,wep,K(N),s(N,N),x(N) CM=219474.63d0 write(6,600) 600 format(' Energy ground mode shift x') wep=we do 1 i=1,NQ x(i)=E(i)*K(i)**2/2 wep=wep-(E(i)*K(i))**2/2.0d0 write(6,601)E(i)*CM,K(i),x(i) 1 if(mod(i,2).eq.0)write(6,*) 601 format(3f10.2,$) write(6,602)we*CM,wep*CM 602 format(/,' Vertical and interpolated energy (cm-1):',2f12.2) call shgeo(N,NQ,s,K,'FI.X') end subroutine shgeo(N,NQ,s,K,f) implicit none integer*4 N,NQ,i,j real*8 s(N,N),K(N),BOHR real*8,allocatable::X0(:),X(:) integer*4,allocatable::Z(:) character*(*) f BOHR=0.529177d0 allocate(X0(N),X(N),Z(N/3)) call READX(N,X0,Z) do 2 i=1,N X(i)=X0(i) do 2 j=1,NQ 2 X(i)=X(i)+s(i,j)*K(j)*BOHR call wrx(Z,X,N,f) write(6,603) 603 format(' Interpolated geometry in ',$) do 1 i=1,len(f) 1 write(6,604)f(i:i) 604 format(a1,$) write(6,*) return end subroutine wrx(Z,X,N,s) implicit none integer*4 N,i,ix,Z(N/3) real*8 X(N) character*(*) s open(9,file=s) write(9,*)'interpolated geometry' write(9,*)N/3 do 1 i=1,N/3 1 write(9,90)Z(i),(X(ix+3*(i-1)),ix=1,3) 90 format(i6,3f12.6) close(9) return end subroutine mws(N,s,ss,Z,NQ) implicit none integer*4 MENDELEV,N,Z(N),NQ,i,ia,ix,j,k real*8 SFAC,amu,ma,s(N,N),ss(N,N) real*8,allocatable::e(:,:) parameter (Mendelev=89) real*4 amas(MENDELEV) data amas/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,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/ SFAC=0.0234280d0 amu=1.0d0/SFAC**2 do 1 ia=1,N/3 ma=dsqrt(dble(amas(Z(ia)))*amu) do 1 ix=1,3 i=ix+3*(ia-1) do 1 j=1,NQ 1 ss(i,j)=s(i,j)*ma allocate(e(N,N)) do 2 i=1,NQ do 2 j=1,NQ e(i,j)=0.0d0 do 2 k=1,N 2 e(i,j)=e(i,j)+ss(k,i)*ss(k,j) write(6,601) 601 format(' s-test:') write(6,600)e(1 ,1),e(1 ,2),e(1 ,3),e(1 ,NQ-1),e (1,NQ) write(6,600)e(2 ,1),e(2 ,2),e(2 ,3),e(2 ,NQ-1),e (2,NQ) write(6,600)e(3 ,1),e(3 ,2),e(3 ,3),e(3 ,NQ-1),e (3,NQ) write(6,*)' ...' write(6,600)e(NQ-1,1),e(NQ-1,2),e(NQ-1,3),e(NQ-1,NQ-1),e(NQ-1,NQ) write(6,600)e(NQ ,1),e(NQ ,2),e(NQ ,3),e(NQ ,NQ-1),e(NQ ,NQ) 600 format(3f10.6,' ... ',2f10.6) return end subroutine gms(N,NQ,NQE,s,se,K,KP) c transform excite mode shifts into ground ones implicit none integer*4 N,NQ,NQE,i,j,a real*8 s(N,N),se(N,N),K(N),KP(N) real*8 ,allocatable::X(:),ss(:,:),sse(:,:) integer*4 ,allocatable::Z(:) allocate(X(N),Z(N/3),ss(N,N),sse(N,N)) c read atomic numbers call READX(N,X,Z) c ss=mass-weighted s write(6,602)' ground s: ' 602 format(a11) call mws(N,s,ss,Z,NQ) c sse=mass-weighted se write(6,602)' excited s:' call mws(N,se,sse,Z,NQE) do 1 i=1,NQ K(i)=0.0d0 do 1 j=1,N do 1 a=1,NQE 1 K(i)=K(i)+ss(j,i)*sse(j,a)*KP(a) return end