PROGRAM NEW6T 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,NN,N3, 1jjs,iwr,nlim,mmaxi,mmaxf,maxt,LEXCI,LEXCF,iwf,nb0,nkt,ng,nif, 1si,mmf,il,iif,nsf,ii,LNQf,mfc,kk,iii,kf,LNQi,mmi,nsi,maxcen c parameter (ns0=19,maxcen=10) parameter (ns0=19) parameter (maxcen=10) INTEGER*4 nc,fc,ncen(maxcen),fcen(maxcen),fexc(maxcen), 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),wz, 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,f0,g0,f1,g1,CM,dtol,kelvin, 1ktlim,kt,qred,qp,qpar,emax,wrax,enq,ei,bo,ef 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(:,:),uu(:),ff(:,:),gg(:,:),Gw(:,:),eib(:), 1efs(:),bf(:),efb(:),wi(:),wf(:) logical,allocatable::doi(:),dof(:),doif(:) integer*4,allocatable::wmi(:),wmf(:),idi(:),idf(:),smf(:),smi(:), 1wmaxf(:),NQsf(:),NQsi(:),wmaxi(:) LOGICAL asym(5),luseA,luseG,ldo(ns0),lstat,slist,lvert,eattt, 1lt character*80,allocatable::sx(:) C WRITE(6,6666) 6666 FORMAT(/,/, A' * * * * * * * * * * * * * * * * * * * * * *',/, 1' * * ',/, 2' * RESONANCE RAMAN OPTICAL ACTIVITY * ',/, 3' * * ',/, 4' * * * * * * * * * * * * * * * * * * * * * *',/, 8' ',/, A'INPUT FILES: NEW6T.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,kelvin, 1mmaxi,mmaxf,LEXCI,LEXCF,ktlim,wrax) CM=219474.63d0 kt=kelvin*0.6950356d0 emax=(wrax-kt*log(ktlim))/CM OPEN(2,FILE='FILE.TTT',STATUS='OLD') READ(2,*) READ(2,*)NAT close(2) N=3*NAT allocate(grad(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)) 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 allocate(dd(3*N),aa(3*N),vv(3*N),qq(6*N), 1ddi(3*N),rri(3*N),vvi(3*N),qqi(6*N)) write(6,608)asym 608 format(' alpha G GC A AC:',5l2) write(6,607)wexc*CM,Gau*CM 607 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)) if(i.eq.1)then call rduschs(N,NQ,sg,sx(i)(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') Gw=0.0d0 write(6,*)'frequencies read' do 4 f=1,NQ 4 Gw(f,f)=W(f) endif 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(J(NQ,NQ),JT(NQ,NQ),K(NQ),jtj(NQ,NQ),jtji(NQ,NQ), 1U(NQ,NQ),V(NQ,NQ),X(NQ,NQ),Y(NQ,NQ),ki(NQ),sj(NQ,NQ), 1uu(81),ff(81,NQ),gg(81,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 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,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 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 uu,ff,gg: call sc(N,NQ,u0,m0,q0,ddi,rri,qqi,uu,ff,gg) c vertical approximation, make new shift vectors: if(lvert)call mksv(NQ,N,K,gradi,WP,J,W,iwr) c make matrices call m3(N,NQ,ki,U,V,X,Y,W,WP,sj,wz) c partition function c IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII if(i.eq.1)then qp=1.0d0 do 3 f=1,NQ 3 qp=qp/(1.0d0-exp(-W(f)/0.00000316678d0/kelvin)) c make excitations write(6,*) c maximal excitations for initial/final: allocate(wmi(NQ),wmf(NQ)) maxt=max(mmaxf,mmaxi) iwf=max(LEXCF+1,LEXCI+1) c debug: wmi=1 wmf=1 qred=1.0d0 c allocate(bf(100000)) c goto 8989 c number of states, partition function: call nkti(nb0,nkt,ktlim,maxt,iwf,NQ,N,Gw,CM,kt,0,qpar, 1emax,mmaxi,mmaxf,LEXCI,LEXCF,qred,wmi,wmf,nif) c nb0 number of initial states c nkt number of final states c ktlim for exp(-E/kT) c number of states, probabilities: write(6,611)qp,qpar 611 format(' Partition functions analytical numerical:',2e12.4) call nkti(nb0,nkt,ktlim,maxt,iwf,NQ,N,Gw,CM,kt,1,qpar, 1emax,mmaxi,mmaxf,LEXCI,LEXCF,qred,wmi,wmf,nif) c number of states with reduced excitations: call nktr(nb0,nkt,ng,ktlim,maxt,NQ,N,Gw,CM,kt, 1emax,mmaxi,mmaxf,LEXCI,LEXCF,wmi,nif) c nb0: number of states with exp(-E/kT): nc=mmi do 401 ii=1,nc ncen(ii)=smi(ii) 401 nexc(ii)=NQSi(ii) c fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff kf=0 do 2021 mmf=0,mmaxf il=max(1,mmf) allocate(smf(il),NQsf(il),wmaxf(il)) wmaxf=iwf nsf=iwf**mmf do 202 iif=1,mfc(NQ,mmf) call getdis(iif,smf,NQ,mmf,il) c call as3(nsf,mmf,wmaxf,smf,wmf) Do 2023 kk = 1, nsf Call FCA2NQ(mmf,kk,wmaxf,lt,LNQf,NQsf) ef=enq(LNQf,N,smf,Gw,NQsf) if(lt)then if(ef.lt.emax.and.mmf.le.mmaxf.and.LNQf.ne.0. 1 and.ef.gt.ei+1.0d-5)then do 301 ii=1,LNQf 301 if(NQSf(ii).gt.LEXCF)goto 3011 kf=kf+1 if(iwr.gt.0)call wssr(' to: ',kf,ef*CM,LNQf,smf,NQsf) f=f+1 wi(f)=ei wf(f)=ef c |f>: fc=mmf do 402 ii=1,fc fcen(ii)=smf(ii) 402 fexc(ii)=NQSf(ii) iint=0.0d0 jint=0.0d0 call mkijf(N,NQ,W,fc,fcen,fexc,nc,ncen,nexc,uu,ff,gg, 1 iint,jint,X,Y,U,V,ki) c construct i-> f transition polarizabilities: call ctp(we,wexc,iint,jint,asym,aasr,aasi,gsr,gsi, 1 gcsr,gcsi,asr,asi,acsr,acsi,f0,g0,f1,g1) c add them to ~ tensor f-derivatives by pis factor call adss(f,NN,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1 gcsr,gcsi,GGr,GGi,asr,asi,A,Ai,acsr,acsi,AAr,AAi,pis,av2) 3011 continue endif endif 2023 continue 202 continue 2021 deallocate(smf,wmaxf,NQsf) endif endif 201 continue 2011 deallocate(smi,wmaxi,NQsi) c ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff 3012 continue 1000 deallocate(J,JT,K,jtj,jtji,U,V,X,Y,ki,sj,uu,ff,gg) c eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee c N3=f write(6,*)N3,' transitions with positive energy (Stokes)' call initab(ldo) do 2 i=1,N3 c rewrite into to smaller arrays: call trsten(NN,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,wi(i),wf(i),i,alr,ali,gtr,gti,gcr,gci,atr, 1ati,acr,aci,luseA,luseG,ldo) 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 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, 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,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,kelvin, 1mmaxi,mmaxf,LEXCI,LEXCF,ktlim,wrax) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) parameter (ns0=19) LOGICAL lex,luseA,luseG,ldo(ns0),asym(5),lstat,slist,lvert,eattt real*8 ktlim,kelvin 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 maximal initial and final excitation classess,excitations: mmaxi=3 mmaxf=1 LEXCI=3 LEXCF=1 c maximal Raman frequency cm-1: wrax=4000.0d0 c limit for exp(-E/kT): ktlim=0.20d0 c temperature kelvin=300.0d0 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='NEW6T.OPT',exist=lex) if(lex)then OPEN(2,FILE='NEW6T.OPT') WRITE(6,*)' FILE NEW6T.OPT FOUND ...' else goto 99 endif 2001 READ(2,1111)STRING 1111 FORMAT(A15) 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 IF(STRING(1:5).EQ.'KELVI')READ(2,*)kelvin IF(STRING(1:5).EQ.'MMAXI')READ(2,*)mmaxi IF(STRING(1:5).EQ.'MMAXF')READ(2,*)mmaxf IF(STRING(1:5).EQ.'LEXCI')READ(2,*)lexci IF(STRING(1:5).EQ.'LEXCF')READ(2,*)lexcf 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 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,aasr,aasi,ALPHA,ALPHAi,gsr,gsi,G,Gi, 1gcsr,gcsi,Gcr,Gci,asr,asi,A,Ai,acsr,acsi,Acr,Aci,pis,av2) implicit none integer*4 N,i,j,f,l real*8 aasr(3,3),aasi(3,3),ALPHA(N,3,3), 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),pis,av1,av2 av1=0.0d0 av2=0.0d0 do 1 i=1,3 do 1 j=1,3 av1=av1+dabs(ALPHA(f,i,j)) av2=av2+dabs(aasr(i,j)) ALPHA( f,i,j)=ALPHA( f,i,j)+pis*aasr(i,j) ALPHAi(f,i,j)=ALPHAi(f,i,j)+pis*aasi(i,j) G( f,i,j)=G( f,i,j)+pis*gsr( i,j) Gi( f,i,j)=Gi( f,i,j)+pis*gsi( i,j) Gcr( f,i,j)=Gcr( f,i,j)+pis*gcsr(i,j) Gci( f,i,j)=Gci( f,i,j)+pis*gcsi(i,j) do 1 l=1,3 A( f,i,j,l)=A( f,i,j,l)+pis*asr( i,j,l) Ai( f,i,j,l)=Ai( f,i,j,l)+pis*asi( i,j,l) Acr( f,i,j,l)=Acr( f,i,j,l)+pis*acsr(i,j,l) 1 Aci( f,i,j,l)=Aci( f,i,j,l)+pis*acsi(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 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) implicit none integer*4 N,NQ,i,j,k 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(:,:) allocate(su(NQ,NQ),sv(NQ,NQ)) wz=0.0d0 do 1 j=1,NQ wz=wz+WP(j)/2.0d0 do 1 k=1,NQ sv(j,k)=0.0d0 su(j,k)=0.0d0 do 1 i=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=1,NQ wj=W(j) ki(j)=dsqrt(2.0d0*wj) do 2 k=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 function l0(fc,fcen,fexc,nc,ncen,nexc) c fc,nc:number of centers c fcen:list of centers/oscillators c fexc:list of excitations c l0 = c implicit none integer*4 fc,fcen(*),fexc(*),nc,ncen(*),nexc(*),i,j real*8 l0 c : l0=0.0d0 if(nc.ne.fc)return do 4 i=1,nc do 5 j=1,fc 5 if(fcen(j).eq.ncen(i).and.fexc(j).eq.nexc(i))goto 4 return 4 continue l0=1.0d0 return end function l1(fc,fcen,fexc,p,n,nc,ncen,nexc) c fc,nc:number of centers c fcen:list of centers/oscillators c fexc:list of excitations c c p = 1: l1: c p =-1: l1: implicit none integer*4 maxcen parameter (maxcen=10) integer*4 fc,fcen(maxcen),fexc(maxcen),nc,ncen(maxcen), 1nexc(maxcen),n,p,tc,texc(maxcen),tcen(maxcen),i,j,dn,nold real*8 l1 logical le dn=0 do 2 i=1,fc 2 dn=dn+fexc(i) do 3 i=1,nc 3 dn=dn-nexc(i) l1=0.0d0 if(abs(dn).eq.1)then c a+-|n>: call el(tc,tcen,texc,p,n,nc,ncen,nexc,nold,le) if(le)return c : if(tc.ne.fc)return do 4 i=1,tc do 5 j=1,fc 5 if(fcen(j).eq.tcen(i).and.fexc(j).eq.texc(i))goto 4 return 4 continue l1=dsqrt(dble(nold+(p+1)/2)) endif return end function l2(fc,fcen,fexc,nj,nk,j,k,nc,ncen,nexc) c c fc,nc:number of centers c fcen:list of centers/oscillators c fexc:list of excitations implicit none integer*4 maxcen parameter (maxcen=10) integer*4 fc,fcen(maxcen),fexc(maxcen),nc,ncen(maxcen), 1nexc(maxcen),nj,nk,j,k,dn,tc,texc(maxcen),tcen(maxcen), 1uc,uexc(maxcen),ucen(maxcen),i,jj,nold,mold real*8 l2 logical le dn=0 do 2 i=1,fc 2 dn=dn+fexc(i) do 3 i=1,nc 3 dn=dn-nexc(i) l2=0.0d0 if(abs(dn).eq.0.or.abs(dn).eq.2)then c ak+-|n>: call el(tc,tcen,texc,nk,k,nc,ncen,nexc,nold,le) if(le)return c aj+- ak+-|n> = aj+- |t> call el(uc,ucen,uexc,nj,j,tc,tcen,texc,mold,le) if(le)return c : if(uc.ne.fc)return do 4 i=1,uc do 5 jj=1,fc 5 if(fcen(jj).eq.ucen(i).and.fexc(jj).eq.uexc(i))goto 4 return 4 continue l2=dsqrt(dble((nold+(nk+1)/2)*(mold+(nj+1)/2))) endif return end function l3(fc,fcen,fexc,nj,nk,nl,j,k,l,nc,ncen,nexc) c implicit none integer*4 maxcen parameter (maxcen=10) integer*4 fc,fcen(maxcen),fexc(maxcen),nc,ncen(maxcen), 1nexc(maxcen),nj,nk,nl,i,jj,j,k,l,dn,tc,texc(maxcen),tcen(maxcen), 1uc,uexc(maxcen),ucen(maxcen),t1c,t1exc(maxcen),t1cen(maxcen), 1lold,nold,mold real*8 l3 logical le dn=0 do 2 i=1,fc 2 dn=dn+fexc(i) do 3 i=1,nc 3 dn=dn-nexc(i) l3=0.0d0 if(abs(dn).eq.1.or.abs(dn).eq.3)then c al+-|n>: call el(t1c,t1cen,t1exc,nl,l,nc,ncen,nexc,lold,le) if(le)return c ak+- al+-|n>: call el(tc,tcen,texc,nk,k,t1c,t1cen,t1exc,nold,le) if(le)return c aj+- ak+- al+i|n>: call el(uc,ucen,uexc,nj,j,tc,tcen,texc,mold,le) if(le)return c : if(uc.ne.fc)return do 6 i=1,uc do 5 jj=1,fc 5 if(fcen(jj).eq.ucen(i).and.fexc(jj).eq.uexc(i))goto 6 return 6 continue l3=dsqrt(dble((nold+(nk+1)/2)*(mold+(nj+1)/2)*(lold+(nl+1)/2))) endif return end subroutine el(tc,tcen,texc,nl,l,nc,ncen,nexc,nold,le) c |t>= a+-l |n> implicit none integer*4 tc,maxcen,nold,i,nl,l,nc,j parameter (maxcen=10) integer*4 texc(maxcen),tcen(maxcen),nexc(maxcen),ncen(maxcen) logical le le=.false. tc=nc tcen=ncen texc=nexc nold=0 do 7 i=1,nc if(ncen(i).eq.l)then c center l is already excited: nold=nexc(i) texc(i)=texc(i)+nl if(texc(i).eq.0)then do 71 j=i,tc-1 texc(i)=texc(i+1) 71 tcen(i)=tcen(i+1) tc=tc-1 return endif endif 7 continue c center k is not yet excited: if(nl.eq.-1)then c state is annihilated: le=.true. else tc=tc+1 tcen(tc)=l texc(tc)=1 endif return end subroutine mkf0(we,w,Gau,f0,g0,f1,g1) implicit none real*8 we,w,Gau,f0,g0,f1,g1,down down=(we**2-w**2-Gau**2)**2+(2.0d0*we*Gau)**2 f0=(we**2-w**2-Gau**2)/down g0=-2.0d0*we*Gau/down down=((we-w)**2+Gau**2)**2 f1=((we-w)**2-Gau**2)/down g1=-2.0d0*(we-w)*Gau/down 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 mkij(N,NQ,W,fc,fcen,fexc,nc,ncen,nexc,uu,ff,gg, 1iint,jint,X,Y,U,V,ki) implicit none integer*4 N,NQ,maxcen,nn,ii,j,k parameter (maxcen=10) INTEGER*4 nc,fc,ncen(maxcen),fcen(maxcen),fexc(maxcen), 1nexc(maxcen) 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),wn,tw,t, 1l1,l2,l3,t1,t2 do 2 nn=1,NQ wn=W(nn) tw=dsqrt(2.0d0*wn) t=(l1(fc,fcen,fexc, 1,nn,nc,ncen,nexc) 1+ l1(fc,fcen,fexc,-1,nn,nc,ncen,nexc))/tw do 2 ii=1,81 2 iint(ii)=iint(ii)+ff(ii,nn)*t do 3 j=1,NQ do 3 k=1,NQ t1=(l2(fc,fcen,fexc, 1, 1,j,k,nc,ncen,nexc) 1 +l2(fc,fcen,fexc,-1,-1,j,k,nc,ncen,nexc))*X(j,k) 1 +(l2(fc,fcen,fexc,-1, 1,j,k,nc,ncen,nexc) 1 +l2(fc,fcen,fexc, 1,-1,j,k,nc,ncen,nexc))*Y(j,k) 1 -2.0d0*((l1(fc,fcen,fexc, 1,k,nc,ncen,nexc) 1 +l1(fc,fcen,fexc,-1,k,nc,ncen,nexc))*ki(j))*V(j,k) t2=U(j,k)*(l1(fc,fcen,fexc, 1,j,nc,ncen,nexc) 1 -l1(fc,fcen,fexc,-1,j,nc,ncen,nexc)) 1*2.0d0*dsqrt(2.0d0/W(k)) do 31 ii=1,81 31 jint(ii)=jint(ii)+uu(ii)*t1+gg(ii,k)*t2 do 3 nn=1,NQ t =(l3(fc,fcen,fexc, 1, 1, 1,nn,j,k,nc,ncen,nexc) 1 +l3(fc,fcen,fexc, 1,-1,-1,nn,j,k,nc,ncen,nexc))*X(j,k) 1 +(l3(fc,fcen,fexc, 1,-1, 1,nn,j,k,nc,ncen,nexc) 1 +l3(fc,fcen,fexc, 1, 1,-1,nn,j,k,nc,ncen,nexc))*Y(j,k) 1 -2.0d0*((l2(fc,fcen,fexc, 1, 1,nn,k,nc,ncen,nexc) 1 +l2(fc,fcen,fexc, 1,-1,nn,k,nc,ncen,nexc))*ki(j))*V(j,k) 1 +(l3(fc,fcen,fexc,-1, 1, 1,nn,j,k,nc,ncen,nexc) 1 +l3(fc,fcen,fexc,-1,-1,-1,nn,j,k,nc,ncen,nexc))*X(j,k) 1 +(l3(fc,fcen,fexc,-1,-1, 1,nn,j,k,nc,ncen,nexc) 1 +l3(fc,fcen,fexc,-1, 1,-1,nn,j,k,nc,ncen,nexc))*Y(j,k) 1 -2.0d0*((l2(fc,fcen,fexc,-1, 1,nn,k,nc,ncen,nexc) 1 +l2(fc,fcen,fexc,-1,-1,nn,k,nc,ncen,nexc))*ki(j))*V(j,k) t=t/dsqrt(2.0d0*W(nn)) do 3 ii=1,81 3 jint(ii)=jint(ii)+ff(ii,nn)*t return end function plus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) c A |nkp> = a+k |n>, A=sqrt(plus) implicit none integer*4 maxcen parameter (maxcen=10) integer*4 nkp,nc,nkpcen(maxcen),nkpexc(maxcen),ncen(maxcen), 1nexc(maxcen),k,a,plus nkp=nc nkpexc=nexc nkpcen=ncen 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 minus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) c A |nkp> = a-k |n>, A=sqrt(minus) implicit none integer*4 maxcen parameter (maxcen=10) integer*4 nkp,nc,nkpcen(maxcen),nkpexc(maxcen),ncen(maxcen), 1nexc(maxcen),k,a,minus nkp=nc nkpexc=nexc nkpcen=ncen do 1 a=1,nc if(k.eq.ncen(a))then nkpexc(a)=nexc(a)-1 minus=nexc(a) return endif 1 continue minus=0 return end subroutine mkijf(N,NQ,W,fc,fcen,fexc,nc,ncen,nexc,uu,ff,gg, 1iint,jint,X,Y,U,V,ki) c should be fast version implicit none integer*4 N,NQ,maxcen,nn,ii,j,k,akp,akm,app,amm,apm,amp,anp, 1anm parameter (maxcen=10) INTEGER*4 nc,ncen(maxcen),nexc(maxcen),fc,fcen(maxcen), 1fexc(maxcen),nkp,nkpcen(maxcen),nkpexc(maxcen), 1nkm,nkmcen(maxcen),nkmexc(maxcen), 1njmkm,jmkmcen(maxcen),jmkmexc(maxcen), 1njpkm,jpkmcen(maxcen),jpkmexc(maxcen), 1njmkp,jmkpcen(maxcen),jmkpexc(maxcen), 1njpkp,jpkpcen(maxcen),jpkpexc(maxcen), 1nap,napcen(maxcen),napexc(maxcen), 1nam,namcen(maxcen),namexc(maxcen), 1plus,minus real*8 W(N),U(NQ,NQ),V(NQ,NQ),X(NQ,NQ),Y(NQ,NQ),ki(NQ),l1, 1uu(81),ff(81,NQ),gg(81,NQ),iint(81),jint(81),wn,tw,t,l0, 1t1,t2,t3,kp,km do 2 nn=1,NQ wn=W(nn) tw=dsqrt(2.0d0*wn) t=(l1(fc,fcen,fexc, 1,nn,nc,ncen,nexc) 1+ l1(fc,fcen,fexc,-1,nn,nc,ncen,nexc))/tw do 2 ii=1,81 2 iint(ii)=iint(ii)+ff(ii,nn)*t do 3 k=1,NQ c a+k |i> akp=plus(nkp,nkpcen,nkpexc,k,nc,ncen,nexc) c a-k |i> akm=minus(nkm,nkmcen,nkmexc,k,nc,ncen,nexc) if(akp.ne.0.or.akm.ne.0)then kp=l0(fc,fcen,fexc,nkp,nkpcen,nkpexc)*dsqrt(dble(akp)) 1 +l0(fc,fcen,fexc,nkm,nkmcen,nkmexc)*dsqrt(dble(akm)) km=l0(fc,fcen,fexc,nkp,nkpcen,nkpexc)*dsqrt(dble(akp)) 1 -l0(fc,fcen,fexc,nkm,nkmcen,nkmexc)*dsqrt(dble(akm)) t1=0.0d0 t2=0.0d0 do 9 j=1,NQ c a+j a+k |i>,etc. app=plus (njpkp,jpkpcen,jpkpexc,j,nkp,nkpcen,nkpexc)*akp apm=plus (njpkm,jpkmcen,jpkmexc,j,nkm,nkmcen,nkmexc)*akm amp=minus(njmkp,jmkpcen,jmkpexc,j,nkp,nkpcen,nkpexc)*akp amm=minus(njmkm,jmkmcen,jmkmexc,j,nkm,nkmcen,nkmexc)*akm if(app.ne.0)t2=t2+ 1 l0(fc,fcen,fexc,njpkp,jpkpcen,jpkpexc)*dsqrt(dble(app))*X(j,k) if(amm.ne.0)t2=t2+ 1 l0(fc,fcen,fexc,njmkm,jmkmcen,jmkmexc)*dsqrt(dble(amm))*X(j,k) if(amp.ne.0)t2=t2+ 1 l0(fc,fcen,fexc,njmkp,jmkpcen,jmkpexc)*dsqrt(dble(amp))*Y(j,k) if(apm.ne.0)t2=t2+ 1 l0(fc,fcen,fexc,njpkm,jpkmcen,jpkmexc)*dsqrt(dble(apm))*Y(j,k) t1=t1+ki(j)*V(j,k) t=U(j,k)*km*2.0d0*dsqrt(2.0d0/W(j)) do 9 ii=1,81 9 jint(ii)=jint(ii)+gg(ii,j)*t c j t1=t2-2.0d0*t1*kp do 32 ii=1,81 32 jint(ii)=jint(ii)+uu(ii)*t1 endif c akm akp <> 0 3 continue c k do 6 nn=1,NQ anm=minus(nam,namcen,namexc,nn,fc,fcen,fexc) anp=plus( nap,napcen,napexc,nn,fc,fcen,fexc) if(anm.ne.0.or.anp.ne.0)then t3=0.0d0 do 7 k=1,NQ akp=plus( nkp,nkpcen,nkpexc,k,nc,ncen,nexc) akm=minus(nkm,nkmcen,nkmexc,k,nc,ncen,nexc) if(akp.ne.0.or.akm.ne.0)then do 4 j=1,NQ app=plus (njpkp,jpkpcen,jpkpexc,j,nkp,nkpcen,nkpexc)*akp apm=plus (njpkm,jpkmcen,jpkmexc,j,nkm,nkmcen,nkmexc)*akm amp=minus(njmkp,jmkpcen,jmkpexc,j,nkp,nkpcen,nkpexc)*akp amm=minus(njmkm,jmkmcen,jmkmexc,j,nkm,nkmcen,nkmexc)*akm if(app+apm+amp+amm.ne.0)then t3=t3+ 1(l0(nam,namcen,namexc,njpkp,jpkpcen,jpkpexc)*dsqrt(dble(app*anm)) 1+l0(nam,namcen,namexc,njmkm,jmkmcen,jmkmexc)*dsqrt(dble(amm*anm)) 1+l0(nap,napcen,napexc,njpkp,jpkpcen,jpkpexc)*dsqrt(dble(app*anp)) 1+l0(nap,napcen,napexc,njmkm,jmkmcen,jmkmexc)*dsqrt(dble(amm*anp))) 1*X(j,k)+ 1(l0(nam,namcen,namexc,njmkp,jmkpcen,jmkpexc)*dsqrt(dble(amp*anm)) 1+l0(nam,namcen,namexc,njpkm,jpkmcen,jpkmexc)*dsqrt(dble(apm*anm)) 1+l0(nap,napcen,napexc,njmkp,jmkpcen,jmkpexc)*dsqrt(dble(amp*anp)) 1+l0(nap,napcen,napexc,njpkm,jpkmcen,jpkmexc)*dsqrt(dble(apm*anp))) 1*Y(j,k)-2.0d0*ki(j)*V(j,k)* 1(l0(nam,namcen,namexc,nkp,nkpcen,nkpexc)*dsqrt(dble(akp*anm)) 1+l0(nam,namcen,namexc,nkm,nkmcen,nkmexc)*dsqrt(dble(akm*anm)) 1+l0(nap,napcen,napexc,nkp,nkpcen,nkpexc)*dsqrt(dble(akp*anp)) 1+l0(nap,napcen,napexc,nkm,nkmcen,nkmexc)*dsqrt(dble(akm*anp))) endif 4 continue endif 7 continue t3=t3/dsqrt(2.0d0*W(nn)) do 8 ii=1,81 8 jint(ii)=jint(ii)+ff(ii,nn)*t3 endif 6 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) implicit none integer*4 NQ,N,i,p,iwr real*8 K(NQ),WP(N),gradi(N),J(NQ,NQ),W(N),CM,a CM=219474.63d0 K=0.0d0 do 1 i=1,NQ do 1 p=1,NQ 1 K(i)=K(i)+J(i,p)*gradi(p)/WP(p)**2 a=0.0d0 do 3 i=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=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 nkti(ki,kf,ktlim,mx,w,NQ,N,G,CM,kt,ic,qpar, 1emax,mmaxi,mmaxf,LEXCI,LEXCF,qred, 1wmi,wmf,nif) c ki: number of initial states with exp(-E/kT): 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') endif 3011 continue endif if(b.gt.ktlim.and.mi.le.mmaxi)then do 302 j=1,LNQi 302 if(NQSi(j).gt.LEXCI)goto 3012 ki=ki+1 ini=.true. bf(ki)=b/qred prec=prec+100.0d0*b/qpar if(ic.eq.0)then write(6,601)ki,e*CM,b,b/qpar*100.0d0,b/qred*100.0d0 601 format(i6,f12.2,f11.4,f9.4,f8.3,' |',$) if(LNQi.eq.0)write(6,603)0 603 format(i4,' ',$) do 6 j=1,LNQi 6 write(6,600)smi(j),NQsi(j) write(6,605) endif 3012 continue endif if(ini.or.inf)nif=nif+1 endif 3 continue 22 continue 7 deallocate(smi,NQsi,wmaxi) if(ic.eq.0)write(6,400)ki,kg,prec 400 format(/,i6,' initial states of',i11,' tried,', 1f9.4,' % population',/) if(ic.eq.1)write(6,404)kf,kg,nif 404 format(/,i6,' final states of',i11,/, 1i6,' states initial or final',/) return end subroutine as3(nsf,mmf,wmaxf,smf,wmf) integer*4 nsf,mmf,wmaxf(*),smf(*),kk,wmf(*) nsf=1 do 681 kk=1,mmf wmaxf(kk)=wmf(smf(kk)) 681 nsf=nsf*wmaxf(kk) return end subroutine wssr(s,kf,ecm,LNQi,smi,NQsi) integer*4 kf,LNQi,smi(*),NQsi(*),j real*8 ecm character*7 s write(6,6)s 6 format(a7,$) write(6,6011)kf,ecm 6011 format(i6,f12.2,' |',$) do 8 j=1,LNQi 8 write(6,600)smi(j),NQsi(j) 600 format(i4,'^',i2,$) write(6,605) 605 format('>') return end