program guvcde c c polarizability from TDDFT by SOS c implicit none character*80 fo,fov,s80,output character*10 s10 integer n,ie,i,j,k,l,np,nmo,ix,jx,ia,nat,grid(3),iep,icpl, 1oa,ob,idum,io,iz,nmo2,na,nb,nd0,iy,ii,nw, 2nt,nst,ifix,iclose,igiao,ntr,nc0,npp,istart,iend,nmms, 2ncmm,n22,nroot,LEXCL,np0,iwr parameter (nd0=21,nc0=4) c nd0 .. number of elements c 1.. 3 el. dipole, c 4.. 6 gradient, c 7.. 9 magnetic dipole, c 10..12 magnetic "GIAO" c 13..18 quadrupole c 19..21 GIAO c nc0 .. number of dipole elements correlated c ncmm ... maximal number of CI coefficients in one state real*8 uux,uuy,uuz,DD,xp,yp,zp,px,py,pz,ddu(3),p0,dw,w, 1p1,p2,p3,de1,de2,de3,rorw,roiw,e00,Gnj(9), 1dg1,dg2,dg3,ucx,ucy,ucz,cijc, 1ax,ay,az,wden,sum,wmin,wmax,gammaau, 1fr,ds1,rs,angle,fi,bohr,djg(nd0),fkj,fkg,ejg,fjg,gammanm, 1rx,ry,rz,dx,dy,dz, 1rsr1,rsv1,rsr2,rsv2,un(3),u0(nd0),dsr,dsv,dkg(nd0),djk(nd0),ekj, 1ekg,rsv,rsr,debye,rsr3,rsv3,pi,dsm,dsq, 1rsa_4,rsa1_4,rsa2_4,rsa3_4, 1rsa_5,rsa1_5,rsa2_5,rsa3_5, 1rsa_6,rsa1_6,rsa2_6,rsa3_6, 1rsa_7,rsa1_7,rsa2_7,rsa3_7, 1degper,degper0,deglim,ecm,u0d(nd0),n3,cc,dst,lambda, 1xx(nc0),xy(nc0),yy(nc0),xa(nc0),ya(nc0),ar(nc0),reg(nc0), 1qxx,qxy,qxz,qyy,qyz,sumt, 1qzz,qi1,qi2,cd2,d2,sumr,sumi,u1kg,u2kg,u3kg,u1kj,u2kj,u3kj, 1qx,qy,qz,hx,hy,hz,com1,com2,com3,org1,org2,org3,org11,org22,org33, 1dif1,dif2,dif3,ecdo,ugx,ugy,ugz,fne,dj1,dj2,dj3,ujx,ujy,ujz, 2brx,bry,brz,qqx,qqy,qqz,ecda, 3org1c,org11c,org2c,org22c,THRC c x atomic coordinates c v atomic weights c iz atomic number logical lden,lzmat,lw,lrds,lpx,lmcd,lcub,ldeg,lopen,lgnm, 1lwrt,lnorm,lgamma,lort,lgau,lwwe,lturbo,ltda,mcdv,LQ1,LDD, 1lgnj,LDE,rgnj,lbck real*8, allocatable::cij(:),rrx(:),orb(:),rry(:),rrz(:),ror(:), 1rix(:),riy(:),riz(:),roi(:),rot(:),cmmg(:),cmme(:), 1rrxw(:),rryw(:),rrzw(:),rixw(:),riyw(:),rizw(:), 1bp(:),dij(:),bpd(:),pre(:),pred(:),dl(:,:),dv(:,:),r(:,:),ut(:,:), 1e(:),ev(:),r0v(:),r0l(:),qv(:,:,:),eau(:),x(:,:),ut6(:,:), 1xau(:,:),wr(:),wi(:),dlrec(:,:),dvrec(:,:),rrec(:,:),cijb(:), 1hd(:),pre6(:),pred6(:),dlg(:),qvrec(:,:,:),wt(:),dijb(:) c cij - alpha orbitals c cijb - alpha orbitals, backward excitations c dij - beta orbitals integer*4, allocatable::ni(:),aij(:),bij(:),ind(:),jnd(:), 1daij(:),dbij(:),nid(:),qq(:),nib(:),aijb(:),bijb(:),daijb(:), 1dbijb(:),nibd(:),immg(:),imme(:) character*1 s1 c MCDV variables: integer*4 nq,nvib c nq -- 3nat-6, dimension c nvib -- number of the actual normal modes real*8 vu0(3),f00,g0(9),mu0(3) real*8,allocatable::VC(:,:),VD(:),VG(:,:),VGP(:,:),ddi(:), 1gnji(:),aai(:) c pi=4.0d0*datan(1.0d0) bohr=0.529177d0 debye=2.541765d0 nmo=0 write(6,6000) 6000 format(/,' This program reads TDDFT Gaussian output,',/, 1 ' produces UVCD abd absorption spectral tables,',/, 1 ' MCD spectrum by SOS computation,',/, 1 ' and SOS polarizabilities.',/) call readoptions(lden,lmcd,lrds,lw,wden,fo,icpl, 1lzmat,ldeg,lwrt,DD,degper0,deglim,lnorm,nst,ifix,nmo,igiao, 2gammaau,wmin,wmax,np,lgamma,lort,gammanm,lgnm,lgau,lwwe,nmms, 3mcdv,nroot,fov,LEXCL,LQ1,np0,iwr,THRC,LDD,lgnj,LDE,rgnj,lbck) allocate(cmmg(nmms),cmme(nmms),immg(nmms),imme(nmms)) c call dimensions(n,fo,nmo,nmo2,nat,lzmat,lturbo,ltda,ncmm) write(6,*)n,' transitions' c remember number of excitations: ntr=n if(n.eq.0)call report('Stop') write(6,*)nmo,' molecular orbitals' if(nmo.eq.0)call report('Stop') write(6,*)nat,' atoms' if(nat.eq.0)call report('Stop') c for cij and ni, allocate also space for deexcitations c allocate for the maximal (ntr) dimension, not to nst: if(ncmm.eq.0)then n22=nmo2 else n22=ncmm endif write(6,*)2*nmo2*n write(6,*)2*n22*n allocate(cij(2*n22*n),aij(n22*n),bij(n22*n),dij(n22*n), 1daij(n22*n),dbij(n22*n),cijb(n22*n),aijb(n22*n),bijb(n22*n), 1daijb(n22*n),dbijb(n22*n),dijb(2*n22*n), 1nib(n),nibd(n),ni(2*n),nid(n),dl(n,3),dv(n,3),r(n,3),e(n),ev(n), 1r0v(n),r0l(n),eau(n),x(3,nat),xau(3,nat),qq(nat),qv(n,3,3), 1ut(3,n),dlrec(n,3),dvrec(n,3),rrec(n,3),qvrec(n,3,3)) if(lturbo)then call readtm(qq,x,bohr,xau,na,nb,nmo2,n,nmo,nst, 1 dl,dv,qv,r,r0v,r0l,e,ev,eau,ni,nid,nat,cij,lopen, 1 aij,bij,ifix,nib,cijb,aijb,bijb,lnorm,lort,lwrt) else call readfile(lzmat,qq,x,bohr,xau,na,nb,nmo2,n,nmo,nst, 1 dl,dv,qv,r,r0v,r0l,e,ev,eau,ni,lden,nid,nat,cij,lopen, 1 aij,bij,dij,daij,dbij,fo,ifix,nib,cijb,aijb,bijb,lnorm,lort,lwrt, 1 dijb,daijb,dbijb,nibd,ltda,lbck) if(ltda)call readtda(na,nb,nmo2,n,nmo,nst,e,ev,eau,ni,lden, 1 nid,cij,lopen,aij,bij,dij,daij,dbij,ifix,nib,lnorm,lort,lwrt, 1 nibd,ncmm) endif c c limit number of excited states: if(nst.gt.0.and.nst.lt.n)then write(6,*)' Number of excited states limited to ',nst n=nst endif c c if desired, arbitrarily perturb degeneracy if(ldeg)then k=0 do 12 i=1,n-1 degper=degper0 do 12 j=i+1,n if(dabs(eau(i)-eau(j)).lt.deglim)then eau(j)=eau(j)+degper ev(j)=eau(j)*27.211384205943d0 ecm=ev(j)*8065.54476345045d0 e(j)=1.0d7/ecm write(6,6003)i,j,e(i),e(j) 6003 format(2i4,f8.2,' -> ',f8.2) degper=degper+degper0 k=k+1 endif 12 continue write(6,6008)k 6008 format('Degeneracy perturbed ',i3,' times.') do 13 i=1,n ev(i)=eau(i)*27.211384205943d0 ecm=ev(i)*8065.54476345045d0 13 e(i)=1.0d7/ecm endif c check the orbital range if(lden)then do 14 i=1,n do 141 j=1,ni(i) if(aij(n22*(i-1)+j).gt.nmo.or.aij(n22*(i-1)+j).lt.1)then write(6,*)' i:',i,' j:',j,' aij:',aij(n22*(i-1)+j),'nmo:',nmo call report('wrong orbital a') endif 141 if(bij(n22*(i-1)+j).gt.nmo.or.bij(n22*(i-1)+j).lt.1) 1 call report('wrong orbital b') do 14 j=1,nid(i) if(daij(n22*(i-1)+j).gt.nmo.or.daij(n22*(i-1)+j).lt.1)then write(6,*)' i:',i,' j:',j,' daij:',daij(n22*(i-1)+j) call report('wrong orbital da') endif if(dbij(n22*(i-1)+j).gt.nmo.or.dbij(n22*(i-1)+j).lt.1)then write(6,*)' i:',i,' j:',j,' dbij:',dbij(n22*(i-1)+j) call report('wrong orbital db') endif 14 continue endif c sign, unit, and Gaussian stupidity correction of dipoles!: do 5 i=1,n do 5 j=1,3 dv(i,j)=-dv(i,j)/eau(i) r( i,j)= r( i,j)/2.0d0 do 5 k=1,3 5 qv(i,j,k)=-qv(i,j,k)/eau(i) c c record spectral tables: c tttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt write(6,*)n,' energies' call titles(fo,'.l.tab',output) call inispec(31,output,'UVCD spectrum in length formalism') write(6,*)output call titles(fo,'.v.tab',output) call inispec(32,output,'UVCD spectrum in velocity formalism') write(6,*)output c do 6 i=1,n call getv(ntr,dl,ds1,r0l,angle,rs,r,i) write(31,3001)i,e(i),ds1,rs,angle 3001 format(i4,f12.2,2f20.10,f8.1) call getv(n,dv,ds1,r0v,angle,rs,r,i) 6 write(32,3001)i,e(i),ds1,rs,angle write(31,3002) write(32,3002) 3002 format(80(1h-)) close(31) close(32) write(6,*) c c computation of polarizabilities from Gaussian transition dipoles: call cpolar('SOS.TTT','spav.prn','spal.prn','spgv.prn','spgl.prn', 1wmax,wmin,np,debye,pi,eau,gammaau,gammanm,lgnm,lgau,dl,dv,r, 2qv,ntr,n) c if(lrds)then c recalculate transition dipole strengths c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd if(.not.lden)call report('lden needed for lrds') inquire(file='PX.MO.SCR.TXT',exist=lpx) if(.not.lpx)call report('*SCR.TXT files needed for lrds') allocate(bp(nd0*nmo2),bpd(nd0*nmo2)) call rwm(bp,nmo,nd0) if(lopen)call rwmd(bpd,nmo) write(6,*) write(6,*)'Transition moment check: r grad m q' call zerstat(xa,ya,xx,yy,xy,ar,nc0) n3=dble(3*n) do 301 ie=1,n c call me( 1 ni,nid,cij,dij,aij,daij,bij,dbij,ie,lopen,djg,bp,bpd,nmo,nmo2, 1 ifix,nib,cijb,aijb,bijb,nd0,n22) do 312 j=1,nd0 312 djg(j)=-djg(j) do 3121 j=4,6 3121 djg(j)=-djg(j)/eau(ie) do 3122 j=13,18 3122 djg(j)=djg(j) ii=0 do 313 j=1,3 call cstat(djg(j ),dl(ie,j),xy(1),xa(1),ya(1),ar(1),xx(1),yy(1)) call cstat(djg(j+3),dv(ie,j),xy(2),xa(2),ya(2),ar(2),xx(2),yy(2)) call cstat(djg(j+6), r(ie,j),xy(3),xa(3),ya(3),ar(3),xx(3),yy(3)) dlrec(ie,j)=djg(j ) dvrec(ie,j)=djg(3+j) rrec( ie,j)=djg(6+j) do 313 k=j,3 ii=ii+1 c 'XX','XY','XZ','YY','YZ','ZZ' call cstat(djg(12+ii),qv(ie,j,k), 1 xy(4),xa(4),ya(4),ar(4),xx(4),yy(4)) qvrec( ie,j,k)=djg(12+ii) 313 qvrec( ie,k,j)=djg(12+ii) 301 if(lwrt)write(6,3006)e(ie),(djg(j),j=1,9),(djg(j),j=13,18), 1 (dl(ie,j),j=1,3),(dv(ie,j),j=1,3),(r(ie,j),j=1,3), 1 ((qv(ie,j,k),k=j,3),j=1,3) 3006 format(f6.1,15f8.4,/,6x,15f8.4) do 107 j=1,3 107 call stn(n*3,ar(j),xa(j),ya(j),reg(j),yy(j),xy(j)) call stn(n*6,ar(4),xa(4),ya(4),reg(4),yy(4),xy(4)) do 1071 j=1,4 write(6,*) if(j.eq.1)write(6,*)'electric:' if(j.eq.2)write(6,*)'gradient:' if(j.eq.3)write(6,*)'magnetic:' if(j.eq.4)write(6,*)'quadrupole:' 1071 write(6,3009)cc(xy(j),n3,xa(j),ya(j),xx(j),yy(j)),reg(j),ar(j) 3009 format(' Correlation coefficient:',f12.6,/, 1 ' Regression :',f12.6,/, 1 ' :',f12.6,/) c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd c c computation of polarizabilities from recalculated dipoles: call cpolar('SOS.r.TTT','spav.r.prn','spal.r.prn','spgv.r.prn', 1'spgl.r.prn',wmax,wmin,np,debye,pi,eau,gammaau,gammanm, 2lgnm,lgau,dlrec,dvrec,rrec,qvrec,ntr,n) c call inispec(31,'ECDL.TAB','UVCD spectrum in length formalism') call inispec(32,'ECDV.TAB','UVCD spectrum in velocity formalism') call zerstat(xa,ya,xx,yy,xy,ar,2) n3=dble(n) do 314 i=1,n call getv(n,dlrec,dsr,r0l,angle,rs,rrec,i) rsr=-(dlrec(i,1)*rrec(i,1)+dlrec(i,2)*rrec(i,2) 1 +dlrec(i,3)*rrec(i,3)) 1 *2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 write(31,3001)i,e(i),dsr,rsr,angle call getv(n,dvrec,dsv,r0v,angle,rs,rrec,i) rsv=-(dvrec(i,1)*rrec(i,1)+dvrec(i,2)*rrec(i,2) 1 +dvrec(i,3)*rrec(i,3)) 1 *2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 call cstat(dsv,dsr,xy(1),xa(1),ya(1),ar(1),xx(1),yy(1)) call cstat(rsv,rsr,xy(2),xa(2),ya(2),ar(2),xx(2),yy(2)) 314 write(32,3001)i,e(i),dsv,rsv,angle write(31,3002) write(32,3002) close(31) close(32) call stn(n,ar(1),xa(1),ya(1),reg(1),yy(1),xy(1)) call stn(n,ar(2),xa(2),ya(2),reg(2),yy(2),xy(2)) write(6,3011)' Grad vs length:', 1 cc(xy(1),n3,xa(1),ya(1),xx(1),yy(1)),' ECD:', 1 cc(xy(2),n3,xa(2),ya(2),xx(2),yy(2)), 1 reg(1) ,' ECD:',reg(2), 1 ar(1) ,' ECD:',ar(2) 3011 format(/,A16,/, 1 ' Correlation coefficients absorption:',f12.6,A5,f12.6,/, 1 ' Regression coefficients absorption:',f12.6,A5,f12.6,/, 1 ' absorption:',f12.6,A5,f12.6,/) if(icpl.ne.0)then c cplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcpl write(6,*)' CPL requested from state ',icpl c nuclear dipole moment: call uun(un,nat,xau,qq) c electric dipole moment: do 4015 j=1,nd0 u0(j)=0.0d0 u0d(j)=0.0d0 if(lopen)then do 4016 i=1,na 4016 u0( j)=u0( j)+bp( i+nmo*(i-1)+nmo2*(j-1)) do 4017 i=1,nb 4017 u0d(j)=u0d(j)+bpd(i+nmo*(i-1)+nmo2*(j-1)) else c u0 = 2 sum(i,occ) : do 4018 i=1,na 4018 u0( j)=u0(j) +bp( i+nmo*(i-1)+nmo2*(j-1))*2.0d0 endif 4015 continue write(6,3101)(un(j),j=1,3),(-u0(j)-u0d(j),j=1,3), 1 (-u0(j)-u0d(j)+un(j),j=1,3) 1 ,((-u0(j)-u0d(j)+un(j))*debye,j=1,3) write(s10,'(i10)')icpl do 4001 istart=1,len(s10) 4001 if(s10(istart:istart).ne.' ')goto 4002 4002 do 4003 iend=len(s10),1,-1 4003 if(s10(iend:iend).ne.' ')goto 4004 4004 call inispec(30,'CPLR.'//s10(istart:iend)//'.TAB', 1 'CPL from state '//s10(istart:iend)//' -length form') call inispec(31,'CPLV.'//s10(istart:iend)//'.TAB', 1 'CPL from state '//s10(istart:iend)//' -velocity form') call inispec(32,'CPLT.'//s10(istart:iend)//'.TAB', 1 'CPL from state '//s10(istart:iend)//' -length+magnetic ') c precalculate b>: allocate (pre(nd0*nmo2),pred(nd0*nmo2)) call getieab(pre,pred,lopen,na,nb,nmo,ni,nid,icpl,aij,bij, 1 cij,daij,dbij,dij,bp,bpd,u0,u0d,nmo2,nd0,n22) c calculate state dipole again and check CI: call getdipole(ut,n,cij,nmo,ni,aij,bij, 1 dij,nid,daij,dbij,lopen,pre,pred,icpl,nd0,3,n22) sum=0.0d0 do 4008 i=1,ni(icpl) 4008 sum=sum+cij(n22*(icpl-1)+i)**2 do 4009 i=1,nid(icpl) 4009 sum=sum+dij(n22*(icpl-1)+i)**2 do 4010 i=1,nib(icpl) 4010 sum=sum-cijb(n22*(icpl-1)+i)**2 write(6,3004)icpl,(un(ia)-ut(ia,icpl),ia=1,3),sum c transitions to ground: c djg=: " singlet (dsqrt(2)) call me( 1 ni,nid,cij,dij,aij,daij,bij,dbij,icpl,lopen,djg,bp,bpd,nmo,nmo2, 1 ifix,nib,cijb,aijb,bijb,nd0,n22) ejg=eau(icpl) do 4011 ia=4,6 4011 djg(ia)=-djg(ia)/ejg c dsv: dsv=djg(3+1)*djg(3+1)+djg(3+2)*djg(3+2)+djg(3+3)*djg(3+3) c dsr: dsr=djg( 1)*djg( 1)+djg( 2)*djg( 2)+djg( 3)*djg( 3) c rsr: rsr=djg(6+1)*djg( 1)+djg(6+2)*djg( 2)+djg(6+3)*djg( 3) c rsv: rsv=djg(6+1)*djg(3+1)+djg(6+2)*djg(3+2)+djg(6+3)*djg(3+3) rsr=rsr*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 rsv=rsv*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 lambda=e(icpl)*10.0d0/bohr dsm= 1 (2.0d0*pi/lambda/eau(icpl))**2*(djg(7)**2+djg(8)**2+djg(9)**2) lambda=1.0d7/(219474.63d0*ejg) write(30,3008)0,lambda,dsr*debye**2,rsr write(31,3008)0,lambda,dsv*debye**2,rsv de1=0.0d0 if(dsr+dsm.ne.0.0d0)de1=rsr/((dsr+dsm)*debye**2) write(32,3008)0,lambda,(dsr+dsm)*debye**2,rsr,de1 c transitions to other states: do 4012 iep=1,n if(iep.ne.icpl)then c djk==> call dodjk(djk,nd0,ni,iep,aij,bij,n22,cij,pre,lopen,daij,dbij, 1 pred,nmo,dij,nid) c transform gradient to dipole: ekj=eau(iep)-ejg if(lgamma)then fkj=ekj/(ekj**2+gammaau**2) else fkj=1.0d0/ekj endif do 4013 ia=4,6 4013 djk(ia)=-djk(ia)*(-fkj) c dsv: dsr=djk(3+1)*djk(3+1)+djk(3+2)*djk(3+2)+djk(3+3)*djk(3+3) c dsr: dsr=djk( 1)*djk( 1)+djk( 2)*djk( 2)+djk( 3)*djk( 3) c rsr: rsr=djk(6+1)*djk( 1)+djk(6+2)*djk( 2)+djk(6+3)*djk( 3) c rsv: rsv=djk(6+1)*djk(3+1)+djk(6+2)*djk(3+2)+djk(6+3)*djk(3+3) rsr=rsr*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 rsv=rsv*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 lambda=-1.0d7/(219474.63d0*ekj)*10.0d0/bohr dsm= 1 (2.0d0*pi/lambda/ekj)**2*(djk(7)**2+djk(8)**2+djk(9)**2) lambda=-1.0d7/(219474.63d0*ekj) write(30,3008)iep,lambda,dsr*debye**2,rsr write(31,3008)iep,lambda,dsv*debye**2,rsv de1=0.0d0 if(dsr+dsm.ne.0.0d0)de1=rsr/((dsr+dsm)*debye**2) write(32,3008)iep,lambda,(dsr+dsm)*debye**2,rsr,de1 endif 4012 continue do 4014 i=30,32 write(i,3002) 4014 close(i) deallocate(pre,pred) endif c cplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcplcpl if(lmcd)then c mcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcd c input parameters for vibrational resolution: if(mcdv)then nq=3*nat-6 allocate(VC(nq,nq),VD(nq),VG(nq,nq),VGP(nq,nq),ddi(9*nat)) allocate(aai(9*nat)) if(lgnj)allocate(gnji(27*nat)) call imcdv(fov,iwr,LDD,nroot,nat,VC,VD,VG,VGP,vu0,ddi, 1 nq,f00,e00,gnji,lgnj,LDE,g0,nvib,mu0,aai) endif c nuclear dipole moment: call uun(un,nat,xau,qq) c c electric dipole, grad,r x grad/2 and (r-rul) x grad/2 etc.: do 310 j=1,nd0 u0(j)=0.0d0 u0d(j)=0.0d0 if(lopen)then do 327 i=1,na 327 u0( j)=u0( j)+bp( i+nmo*(i-1)+nmo2*(j-1)) do 329 i=1,nb 329 u0d(j)=u0d(j)+bpd(i+nmo*(i-1)+nmo2*(j-1)) else c u0 = 2 sum(i,occ) : do 330 i=1,na 330 u0( j)=u0(j) +bp( i+nmo*(i-1)+nmo2*(j-1))*2.0d0 endif 310 continue write(6,3101)(un(j),j=1,3),(-u0(j)-u0d(j),j=1,3), 1 (-u0(j)-u0d(j)+un(j),j=1,3) 1 ,((-u0(j)-u0d(j)+un(j))*debye,j=1,3) 3101 format(/,' Dipole nuclear :',/,3f10.3,/, 1 ' electronic:',/,3f10.3,/, 1 ' total :',/,3f10.3,' debyes:',3f10.3,/) open(28,file='TTENSOR.TAB') write(28,2828) 2828 format(' the MCD tensor, as in',/, 1 ' J. Chem. Theory COmput 2013, 9, 1557, equation (9),',/, 2 ' first index electric',/, 3 ' transition xx',10x,'xy',10x,'xz',10x,'yx',10x,'yy', 4 10x,'yz',10x,'zx',10x,'zy',10x,'zz',/, 4 60(1h-)) open(29,file='MOMENTS.TAB') write(29,*)' electric and effective magnetic moments' call inispec(30,'MCDL.TAB','MCD spectrum in length formalism') call inispec(31,'MCDV.TAB','MCD spectrum in velovity formalism') call inispec(32,'Q.TAB','quadrupole-allowed transitions') call inispec(33,'M.TAB','magnetic-allowed transitions') call inispec(34,'T.TAB','total transition probabilities') call inispec(35,'ECDO.TAB','ECD in LORG formulation') call inispec(36,'ECDA.TAB','ECD in LORGave formulation') if(igiao.eq.6)then call inispec(37,'MCD4.TAB','MCD IGIAO=4') call inispec(38,'MCD5.TAB','MCD IGIAO=5') call inispec(39,'MCD6.TAB','MCD IGIAO=6') call inispec(40,'MCD7.TAB','MCD IGIAO=7') open(41,file='MCD7.TXT') write(41,3014) 3014 format(' MCD7.TXT','MCD IGIAO=7 statistics',/, 1 ' n wavelength(nm)',10x,'Rg',13x,'Re',13x,'R0', 1 11x,'Rtot',/,80(1h-)) else write(s1,'(i1)')igiao call inispec(37,'MCD4.TAB','MCD IGIAO='//s1) endif if(mcdv)call inispec(42,'MCDI.TAB','MCD vibrational') if(mcdv)call inispec(43,'ECDI.TAB','ECD vibrational') write(6,3005) 3005 format(/,' State dipole sum cij**2-cijb**2:') write(6,3004)0,(un(ia)-u0(ia)-u0d(ia),ia=1,3),1.0d0 3004 format(i6,4f10.3) allocate (pre(nd0*nmo2),pred(nd0*nmo2)) c c if LORG (igiao=6), precalculate all and if(igiao.eq.6)then write(6,*)'pre-calculating , x=r,grad' allocate(hd(n*n*6),ut6(6,n),pre6(6*nmo2),pred6(6*nmo2), 1 dlg(6*n)) do 305 ie=1,n c precalculate b>: call getieab(pre6,pred6,lopen,na,nb,nmo,ni,nid,ie,aij,bij, 1 cij,daij,dbij,dij,bp,bpd,u0,u0d,nmo2,6,n22) c diagonal dipole and gradient: call getdipole(ut6,n,cij,nmo,ni,aij,bij, 1 dij,nid,daij,dbij,lopen,pre6,pred6,ie,6,6,n22) c dlg=: " singlet (dsqrt(2)) call me(ni,nid,cij,dij,aij,daij,bij,dbij,ie,lopen,djg,bp, 1 bpd,nmo,nmo2,ifix,nib,cijb,aijb,bijb,6,n22) do 3051 j=1,6 c : dlg(j+6*(ie-1))=djg(j) c : 3051 hd(j+(ie-1+n*(ie-1))*6)=ut6(j,ie) do 305 iep=1,n if(ie.ne.iep)then c djk===> call dodjk(djk,6,ni,iep,aij,bij,n22,cij,pre6,lopen,daij, 1 dbij,pred6,nmo,dij,nid) do 3052 j=1,6 3052 hd(j+(ie-1+n*(iep-1))*6)=djk(j) endif 305 continue deallocate(pre6,pred6) write(6,*)' ... done' endif c call zerstat(xa,ya,xx,yy,xy,ar,3) c 304304304304304304304304304304304304304304304304304304304304304 do 304 ie=1,n ejg=eau(ie) if(lgamma)then fjg=ejg/(ejg**2+gammaau**2) else fjg=1.0d0/ejg endif dsr=0.0d0 dsm=0.0d0 dsq=0.0d0 dst=0.0d0 dsv=0.0d0 rsv1=0.0d0 rsv2=0.0d0 rsr1=0.0d0 rsr2=0.0d0 rsa1_4=0.0d0 rsa2_4=0.0d0 rsa1_5=0.0d0 rsa2_5=0.0d0 com1=0.0d0 com2=0.0d0 com3=0.0d0 org1=0.0d0 org2=0.0d0 org3=0.0d0 org11=0.0d0 org22=0.0d0 org33=0.0d0 c MCD tensor for transition n-> j: c J. Chem. Theory COmput 2013, 9, 1557, equation (9): call vz(Gnj,9) c precalculate b>: call getieab(pre,pred,lopen,na,nb,nmo,ni,nid,ie,aij,bij, 1 cij,daij,dbij,dij,bp,bpd,u0,u0d,nmo2,nd0,n22) c calculate state dipole: call getdipole(ut,n,cij,nmo,ni,aij,bij, 1 dij,nid,daij,dbij,lopen,pre,pred,ie,nd0,3,n22) sum=0.0d0 do 317 i=1,ni(ie) 317 sum=sum+cij(n22*(ie-1)+i)**2 do 326 i=1,nid(ie) 326 sum=sum+dij(n22*(ie-1)+i)**2 do 328 i=1,nib(ie) 328 sum=sum-cijb(n22*(ie-1)+i)**2 write(6,3004)ie,(un(ia)-ut(ia,ie),ia=1,3),sum c djg=: " singlet (dsqrt(2)) call me( 1 ni,nid,cij,dij,aij,daij,bij,dbij,ie,lopen,djg,bp,bpd,nmo,nmo2, 1 ifix,nib,cijb,aijb,bijb,nd0,n22) do 320 ia=4,6 320 djg(ia)=-djg(ia)/eau(ie) c maximum contributions to ground and excited MCD: do 336 i=1,nmms immg(i)=0 imme(i)=0 cmmg(i)=0.0d0 336 cmme(i)=0.0d0 do 307 iep=1,n if(iep.ne.ie)then ekg=eau(iep) ekj=ekg-ejg if(lgamma)then fkj=ekj/(ekj**2+gammaau**2) fkg=ekg/(ekg**2+gammaau**2) else fkj=1.0d0/ekj fkg=1.0d0/ekg endif c dkg=: "" singlet (dsqrt(2)) call me( 1 ni,nid,cij,dij,aij,daij,bij,dbij,iep,lopen,dkg,bp,bpd,nmo,nmo2, 1 ifix,nib,cijb,aijb,bijb,nd0,n22) c djk===> call dodjk(djk,nd0,ni,iep,aij,bij,n22,cij,pre,lopen,daij,dbij, 1 pred,nmo,dij,nid) c transform gradient to dipole: do 319 ia=4,6 djk(ia)=-djk(ia)*(-fkj) 319 dkg(ia)=-dkg(ia)*fkg c ukg: u1kg=djg(2)*djk(3)-djg(3)*djk(2) u2kg=djg(3)*djk(1)-djg(1)*djk(3) u3kg=djg(1)*djk(2)-djg(2)*djk(1) c ukj: u1kj=djg(2)*dkg(3)-djg(3)*dkg(2) u2kj=djg(3)*dkg(1)-djg(1)*dkg(3) u3kj=djg(1)*dkg(2)-djg(2)*dkg(1) c magnetic moment and dipole: p0=(dkg(6+1)*u1kg+dkg(6+2)*u2kg+dkg(6+3)*u3kg)*fkg rsr1=rsr1+p0 com1=com1+p0 p1=(djk(6+1)*u1kj+djk(6+2)*u2kj+djk(6+3)*u3kj)*fkj rsr2=rsr2+p1 com2=com2+p1 if(nroot.eq.ie)then do 15 ix=1,3 do 15 jx=1,3 c G=(J|u|k>/Ekg+/Ekj 15 Gnj(ix+3*(jx-1))=Gnj(ix+3*(jx-1)) 1 +djk(ix)*dkg(6+jx)*fkg+dkg(ix)*djk(6+jx)*fkj endif c c magnetic moment and velocity: rsv1=rsv1+(dkg(6+1)*(djg(2+3)*djk(3+3)-djg(3+3)*djk(2+3)) 1 + dkg(6+2)*(djg(3+3)*djk(1+3)-djg(1+3)*djk(3+3)) 1 + dkg(6+3)*(djg(1+3)*djk(2+3)-djg(2+3)*djk(1+3)))*fkg rsv2=rsv2+(djk(6+1)*(djg(2+3)*dkg(3+3)-djg(3+3)*dkg(2+3)) 1 + djk(6+2)*(djg(3+3)*dkg(1+3)-djg(1+3)*dkg(3+3)) 1 + djk(6+3)*(djg(1+3)*dkg(2+3)-djg(2+3)*dkg(1+3)))*fkj c c GIAO magnetic moment and dipole: if(igiao.ge.4)then c /(Ek-Eg) ~ c common origin terms - p0 c GIAO compensation term: p2=-dkg(9+1)*u1kg-dkg(9+2)*u2kg-dkg(9+3)*u3kg rsa1_4=rsa1_4+p0+p2 rsa1_5=rsa1_5+p0+p2 c c /(Ek-Ej) ~ c common origin terms - p1 c GIAO terms: p2= djk(18+1)*u1kj+djk(18+2)*u2kj+djk(18+3)*u3kj rsa2_4=rsa2_4+p1 rsa2_5=rsa2_5+p1+p2 else rsa1_4=rsa1_4+(dkg(9+1)*u1kg+dkg(9+2)*u2kg+dkg(9+3)*u3kg)*fkg rsa2_4=rsa2_4+(djk(9+1)*u1kj+djk(9+2)*u2kj+djk(9+3)*u3kj)*fkj endif if(igiao.eq.6)then c nonlocal "LORG" like trick: c sum(L,L<>g)x/(El-Eg) etc: call lsums(dx,dy,dz,rx,ry,rz,qx,qy,qz,hx,hy,hz, 1 n,eau,lgamma,gammaau,ie,iep,hd,dlg,na,nb, 2 brx,bry,brz,qqx,qqy,qqz,u0,u0d,ut) c x/(El-Eg)./2: org1c = rx *u1kg+ry *u2kg+rz *u3kg org1 = org1 + org1c c x/(Ek-El)./2: org11c = -brx*u1kg-bry*u2kg-brz*u3kg org11 = org11 + org11c c write(6,6038)'kk',ie,iep,ekg,dkg(9),rz/fkg,u3kg c seems to compensat eperfectly c x/(EK-EL)./2: org2c = qx*u1kj+qy*u2kj+qz*u3kj org2 = org2 +org2c c x/(El-Ej)./2: org22c = -qqx*u1kj-qqy*u2kj-qqz*u3kj org22 = org22 + org22c c write(6,6038)'qq',ie,iep,ekg,djk(9),qz/fkj,u3kj c seems to compensate perfectly if(nroot.eq.ie)then Gnj(1)=Gnj(1)-(djk(1)*(rx-brx)+dkg(1)*(qx-qqx))/2.0d0 Gnj(2)=Gnj(2)-(djk(2)*(rx-brx)+dkg(2)*(qx-qqx))/2.0d0 Gnj(3)=Gnj(3)-(djk(3)*(rx-brx)+dkg(3)*(qx-qqx))/2.0d0 Gnj(4)=Gnj(4)-(djk(1)*(ry-bry)+dkg(1)*(qy-qqy))/2.0d0 Gnj(5)=Gnj(5)-(djk(2)*(ry-bry)+dkg(2)*(qy-qqy))/2.0d0 Gnj(6)=Gnj(6)-(djk(3)*(ry-bry)+dkg(3)*(qy-qqy))/2.0d0 Gnj(7)=Gnj(7)-(djk(1)*(rz-brz)+dkg(1)*(qz-qqz))/2.0d0 Gnj(8)=Gnj(8)-(djk(2)*(rz-brz)+dkg(2)*(qz-qqz))/2.0d0 Gnj(9)=Gnj(9)-(djk(3)*(rz-brz)+dkg(3)*(qz-qqz))/2.0d0 endif endif c store the largest contributions: call dtm(iep,p0-(org1c+org11c)/2.0d0,immg,cmmg,nmms) call dtm(iep,p1-(org2c+org22c)/2.0d0,imme,cmme,nmms) endif c if(iep.ne.ie)then 307 continue c add parts dependent on permanent dipole moment: de1=ut(1,ie)-un(1) de2=ut(2,ie)-un(2) de3=ut(3,ie)-un(3) dg1=u0(1)+u0d(1)-un(1) dg2=u0(2)+u0d(2)-un(2) dg3=u0(3)+u0d(3)-un(3) c -: ddu(1)=de1-dg1 ddu(2)=de2-dg2 ddu(3)=de3-dg3 dif1=djg(2)*ddu(3)-djg(3)*ddu(2) dif2=djg(3)*ddu(1)-djg(1)*ddu(3) dif3=djg(1)*ddu(2)-djg(2)*ddu(1) rsv3=(djg(6+1)*(djg(2+3)*ddu(3)-djg(3+3)*ddu(2)) 1 +djg(6+2)*(djg(3+3)*ddu(1)-djg(1+3)*ddu(3)) 1 +djg(6+3)*(djg(1+3)*ddu(2)-djg(2+3)*ddu(1)))*fjg c .x(-): com3=(djg(6+1)*dif1+djg(6+2)*dif2+djg(6+3)*dif3)*fjg rsr3=com3 c c contribution to the MCD tensor c Gnj=Gnj+(uj-un)/Ejn if(nroot.eq.ie)then do 16 ix=1,3 do 16 iy=1,3 16 Gnj(ix+3*(iy-1))=Gnj(ix+3*(iy-1))+ddu(ix)*djg(6+iy)*fjg endif if(igiao.ge.4)then c GIAO parts: p3=-djg( 9+1)*dif1-djg( 9+2)*dif2-djg( 9+3)*dif3 rsa3_4=com3+p3 rsa3_5=com3+p3 else rsa3_4=(djg(9+1)*(djg(2)*ddu(3)-djg(3)*ddu(2)) 1 + djg(9+2)*(djg(3)*ddu(1)-djg(1)*ddu(3)) 1 + djg(9+3)*(djg(1)*ddu(2)-djg(2)*ddu(1)))*fjg endif if(igiao.eq.6)then call lsums(dx,dy,dz,rx,ry,rz,qx,qy,qz,hx,hy,hz, 1 n,eau,lgamma,gammaau,ie,ie,hd,dlg,na,nb, 2 brx,bry,brz,qqx,qqy,qqz,u0,u0d,ut) c x/(El-Eg).(-)/2/Ne: org3=rx*dif1+ry*dif2+rz*dif3 c x/(Ej-El).(-)/2/Ne: org33=-brx*dif1-bry*dif2-brz*dif3 c write(6,6038)'xx',ie,ie,ejg,djg(9),rz/fjg,dif3 rsa1_6=com1-org1 rsa2_6=com2-org2 rsa3_6=com3-org3 rsa1_7=com1-(org1+org11)/2.0d0 rsa2_7=com2-(org2+org22)/2.0d0 rsa3_7=com3-(org3+org33)/2.0d0 if(nroot.eq.ie)then Gnj(1)=Gnj(1)-ddu(1)*(rx-brx)/2.0d0 Gnj(2)=Gnj(2)-ddu(2)*(rx-brx)/2.0d0 Gnj(3)=Gnj(3)-ddu(3)*(rx-brx)/2.0d0 Gnj(4)=Gnj(4)-ddu(1)*(ry-bry)/2.0d0 Gnj(5)=Gnj(5)-ddu(2)*(ry-bry)/2.0d0 Gnj(6)=Gnj(6)-ddu(3)*(ry-bry)/2.0d0 Gnj(7)=Gnj(7)-ddu(1)*(rz-brz)/2.0d0 Gnj(8)=Gnj(8)-ddu(2)*(rz-brz)/2.0d0 Gnj(9)=Gnj(9)-ddu(3)*(rz-brz)/2.0d0 endif endif c <0|r|1>^2: dsr= djg(1)**2+djg(2)**2+djg(3)**2 lambda=e(ie)*10.0d0/bohr c <0|m|1>^2*4*pi/lambda^2/eau^2: dsm=(2.0d0*pi/lambda/eau(ie))**2*(djg(7)**2+djg(8)**2+djg(9)**2) c <0|q|1>^2*4*pi/lambda^2: qxx=djg(13) qxy=djg(14) qxz=djg(15) qyy=djg(16) qyz=djg(17) qzz=djg(18) qi1=qxx**2+qyy**2+qzz**2+2.0d0*(qxy**2+qxz**2+qyz**2) qi2=(qxx+qyy+qzz)**2 dsq=0.1d0*(2.0d0*pi/lambda)**2*(3.0d0*qi1-qi2) dst=dsr+dsq+dsm dsv= djg(4)**2+djg(5)**2+djg(6)**2 rsv=(rsv1+rsv2+rsv3)/2.0d0 rsr=(rsr1+rsr2+rsr3)/2.0d0 rsa_4=(rsa1_4+rsa2_4+rsa3_4)/2.0d0 rsa_5=(rsa1_5+rsa2_5+rsa3_5)/2.0d0 rsa_6=(rsa1_6+rsa2_6+rsa3_6)/2.0d0 rsa_7=(rsa1_7+rsa2_7+rsa3_7)/2.0d0 if(nroot.eq.ie)then Gnj(1)=Gnj(1)/2.0d0 Gnj(2)=Gnj(2)/2.0d0 Gnj(3)=Gnj(3)/2.0d0 Gnj(4)=Gnj(4)/2.0d0 Gnj(5)=Gnj(5)/2.0d0 Gnj(6)=Gnj(6)/2.0d0 Gnj(7)=Gnj(7)/2.0d0 Gnj(8)=Gnj(8)/2.0d0 Gnj(9)=Gnj(9)/2.0d0 endif call cstat(dsv,dsr,xy(1),xa(1),ya(1),ar(1),xx(1),yy(1)) call cstat(rsv,rsr,xy(2),xa(2),ya(2),ar(2),xx(2),yy(2)) call cstat(rsa_4,rsr,xy(3),xa(3),ya(3),ar(3),xx(3),yy(3)) fne=0.5d0/dble(na+nb)*eau(ie) dg1=-u0(1)-u0d(1) dg2=-u0(2)-u0d(2) dg3=-u0(3)-u0d(3) c (1/2) x /Nel: ugx=(dg2*djg(6)-dg3*djg(5))*fne ugy=(dg3*djg(4)-dg1*djg(6))*fne ugz=(dg1*djg(5)-dg2*djg(4))*fne ecdo=-(djg(1)*(djg(6+1)-ugx) 1 + djg(2)*(djg(6+2)-ugy) 1 + djg(3)*(djg(6+3)-ugz)) 1 *2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 dj1=-ut(1,ie) dj2=-ut(2,ie) dj3=-ut(3,ie) c (1/2) x /Nel: ujx=(dj2*djg(6)-dj3*djg(5))*fne ujy=(dj3*djg(4)-dj1*djg(6))*fne ujz=(dj1*djg(5)-dj2*djg(4))*fne ecda=-(djg(1)*(djg(6+1)-(ugx+ujx)/2.0d0) 1 + djg(2)*(djg(6+2)-(ugy+ujy)/2.0d0) 1 + djg(3)*(djg(6+3)-(ugz+ujz)/2.0d0)) 1 *2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 write(28,2829)ie,((Gnj(iy+3*(ix-1)),iy=1,3),ix=1,3) 2829 format(i4,9f12.6) write(29,3010)ie,e(ie),(djg(ia),ia=1,3),(djg(6+ia),ia=1,3) 3010 format(i4,f12.4,6G15.6) write(30,3008)ie,e(ie),dsr*debye**2,rsr write(31,3008)ie,e(ie),dsv*debye**2,rsv write(32,3008)ie,e(ie),dsq*debye**2 write(33,3008)ie,e(ie),dsm*debye**2 write(34,3008)ie,e(ie),dst*debye**2 write(35,3008)ie,e(ie),dsr*debye**2,ecdo write(36,3008)ie,e(ie),dsr*debye**2,ecda 3008 format(i4,f12.4,2f20.10,f12.5) if(nroot.eq.ie)then call wrgnj(ie,Gnj,djg) if(mcdv)call mcdi(e00,f00,np0,iwr,nq,VC,VD,VGP,LEXCL,vu0, 1 ddi,Gnj,LQ1,LDD,THRC,lwrt,lgnj,gnji,g0,nvib,mu0,aai,rgnj) endif if(igiao.eq.6)then write(37,3008)ie,e(ie),dsr*debye**2,rsa_4 write(38,3008)ie,e(ie),dsr*debye**2,rsa_5 write(39,3008)ie,e(ie),dsr*debye**2,rsa_6 write(40,3008)ie,e(ie),dsr*debye**2,rsa_7 write(41,3012)ie,e(ie),rsa1_7/2.0d0,rsa2_7/2.0d0, 1 rsa3_7/2.0d0,rsa_7 3012 format(i4,f12.4,4f15.8) do 337 i=1,nmms 337 write(41,3013)immg(i),cmmg(i),cmme(i),imme(i) 3013 format(10x,'g ->',i4,f13.8,f15.8,' e ->',i4) else write(37,3008)ie,e(ie),dsr*debye**2,rsa_4 endif 304 continue c 304304304304304304304304304304304304304304304304304304304304304 deallocate(pre,pred) if(igiao.eq.6)then j=41 else j=37 endif do 401 i=28,j write(i,3002) 401 close(i) if(mcdv)then write(42,3002) write(43,3002) close(42) close(43) endif do 109 j=1,3 109 call stn(n,ar(j),xa(j),ya(j),reg(j),yy(j),xy(j)) write(6,3011)' Grad vs length:', 1 cc(xy(1),n3,xa(1),ya(1),xx(1),yy(1)),' MCD:', 1 cc(xy(2),n3,xa(2),ya(2),xx(2),yy(2)), 1 reg(1) ,' MCD:',reg(2), 1 ar(1) ,' MCD:',ar(2) write(6,3011)' GIAO vs length:', 1 1.0d0,' MCD:',cc(xy(3),n3,xa(3),ya(3),xx(3),yy(3)), 1 1.0d0,' MCD:',reg(3), 1 1.0d0,' MCD:',ar(3) endif c mcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcdmcd c deallocate(bp,bpd) endif c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd c c polarization density c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp if(lden)then if(lopen)call report('lden only for closed shells') if(lw)then c initiate atomic weights if(lwwe)then npp=nat*NP else npp=nat endif allocate(wr(npp),wi(npp),wt(npp)) call vz(wi,npp) call vz(wr,npp) call vz(wt,npp) endif inquire(file='orbitals.cub',exist=lcub) if(.not.lcub)call report('orbitals.cub not present') write(6,*) 1 'Polarization density will be calculated with orbitals.cub' open(8,file='orbitals.cub') open(9,file='roi.cub') open(91,file='ror.cub') open(92,file='rot.cub') read(8,2000)s80 2000 format(a80) write(9,2000)s80 write(91,2000)s80 write(92,2000)s80 read(8,*) write(9,901) write(91,901) write(92,901) 901 format(' SCF Total Density') read(8,*)nat,ax,ay,az nat=iabs(nat) write(9,2001)nat,ax,ay,az write(91,2001)nat,ax,ay,az write(92,2001)nat,ax,ay,az 2001 format(i5,3f12.6) read(8,*)grid(1),px write(9 ,2001)grid(1),px,0.0d0,0.0d0 write(91,2001)grid(1),px,0.0d0,0.0d0 write(92,2001)grid(1),px,0.0d0,0.0d0 read(8,*)grid(2),py,py write(9 ,2001)grid(2),0.0d0,py,0.0d0 write(91,2001)grid(2),0.0d0,py,0.0d0 write(92,2001)grid(2),0.0d0,py,0.0d0 read(8,*)grid(3),pz,pz,pz write(9 ,2001)grid(3),0.0d0,0.0d0,pz write(91,2001)grid(3),0.0d0,0.0d0,pz write(92,2001)grid(3),0.0d0,0.0d0,pz do 902 ia=1,nat read(8,2000)s80 read(s80(18:80),*)(x(i,ia),i=1,3) write(9 ,2000)s80 write(91,2000)s80 902 write(92,2000)s80 read(8,905)j,(idum,i=1,j) 905 format(10i5) if(j.ne.nmo)call report('inconsistent number of orbitals') allocate(rrx(grid(3)),orb(grid(3)*nmo),rry(grid(3)),rrz(grid(3)), 1 rix(grid(3)),riy(grid(3)),riz(grid(3)), 1 roi(grid(3)),ror(grid(3)),rot(grid(3))) if(lwwe) 1 allocate(rrxw(grid(3)*np),rryw(grid(3)*np),rrzw(grid(3)*np), 1 rixw(grid(3)*np),riyw(grid(3)*np),rizw(grid(3)*np)) xp=ax-px do 906 ix=1,grid(1) write(6,*)ix,'of',grid(1) xp=xp+px yp=ay-py do 906 iy=1,grid(2) yp=yp+py c do 904 i=1,grid(3) rrx(i)=0.0d0 rry(i)=0.0d0 rrz(i)=0.0d0 rix(i)=0.0d0 riy(i)=0.0d0 904 riz(i)=0.0d0 if(lwwe)then do 9041 i=1,grid(3)*np rrxw(i)=0.0d0 rryw(i)=0.0d0 rrzw(i)=0.0d0 rixw(i)=0.0d0 riyw(i)=0.0d0 9041 rizw(i)=0.0d0 endif read(8,4000)((orb(nmo*(iz-1)+io),io=1,nmo),iz=1,grid(3)) c dw=(wmax-wmin)/dble(np-1) do 907 ie=1,n uux=dv(ie,1)*eau(ie) uuy=dv(ie,2)*eau(ie) uuz=dv(ie,3)*eau(ie) ii=nmo**2*(ie-1) do 907 j=1,ni(ie) cijc=cij(ii+j) ucx=cijc*uux ucy=cijc*uuy ucz=cijc*uuz oa=aij(ii+j) ob=bij(ii+j) c for all frequencies within wmin .. wmax - only |ro|: if(lwwe)then w=wmin-dw do 9071 k=1,np w=w+dw call setv(eau(ie),w,fr,fi,gammaau,gammanm,lgnm,lgau) 9071 call addro(eau(ie),fi,fr,ucx,ucy,ucz,grid(3),grid(3)*(k-1), 1 orb,oa,ob,nmo,rrxw,rryw,rrzw,rixw,riyw,rizw) endif c for requested wden: call setv(eau(ie),wden,fr,fi,gammaau,gammanm,lgnm,lgau) 907 call addro(eau(ie),fi,fr,ucx,ucy,ucz,grid(3),0, 1 orb,oa,ob,nmo,rrx,rry,rrz,rix,riy,riz) c zp=az-pz do 908 iz=1,grid(3) zp=zp+pz ror(iz)=2.0d0*(rrx(iz)**2+rry(iz)**2+rrz(iz)**2) roi(iz)=2.0d0*(rix(iz)**2+riy(iz)**2+riz(iz)**2) rot(iz)=dsqrt(ror(iz)+roi(iz)) roi(iz)=dsqrt(roi(iz)) ror(iz)=dsqrt(ror(iz)) if(lw)then c find closest atom iclose=1 cd2=(xp-x(1,1))**2+(yp-x(2,1))**2+(zp-x(3,1))**2 do 2007 l=2,nat d2 =(xp-x(1,l))**2+(yp-x(2,l))**2+(zp-x(3,l))**2 if(d2.lt.cd2)then cd2=d2 iclose=l endif 2007 continue c increment the real and imaginary weights by the densities: if(lwwe)then c all frequencies: w=wmin-dw do 9081 k=1,np w=w+dw ii=iz+grid(3)*(k-1) rorw=2.0d0*(rrxw(ii)**2+rryw(ii)**2+rrzw(ii)**2) roiw=2.0d0*(rixw(ii)**2+riyw(ii)**2+rizw(ii)**2) ii=iclose+nat*(k-1) wr(ii)=wr(ii)+dsqrt(rorw) wi(ii)=wi(ii)+dsqrt(roiw) 9081 wt(ii)=wt(ii)+dsqrt(rorw+roiw) else c just one frequency: wr(iclose)=wr(iclose)+ror(iz) wi(iclose)=wi(iclose)+roi(iz) wt(iclose)=wt(iclose)+rot(iz) endif endif 908 continue c density only for desired frequency: write(9 ,4000)(roi(iz),iz=1,grid(3)) write(91,4000)(ror(iz),iz=1,grid(3)) 906 write(92,4000)(rot(iz),iz=1,grid(3)) 4000 format(6E13.5) if(lw)then open(50,file='WEIGHTS.TXT') allocate(ind(nat),jnd(nat)) c normalize atomic weights: if(lwwe)then nw=np else nw=1 endif w=wmin-dw do 9082 k=1,nw w=w+dw if(.not.lwwe)w=wden sumr=0.0d0 sumi=0.0d0 sumt=0.0d0 do 2008 i=1,nat ii=i+nat*(k-1) sumi=sumi+wi(ii) sumr=sumr+wr(ii) 2008 sumt=sumt+wt(ii) do 2009 i=1,nat ii=i+nat*(k-1) if(sumi.ne.0.0d0)wi(ii)=wi(ii)/sumi if(sumr.ne.0.0d0)wr(ii)=wr(ii)/sumr 2009 if(sumt.ne.0.0d0)wt(ii)=wt(ii)/sumt write(6,6071)(100.0d0*wr(i+nat*(k-1)),i=1,nat) write(6,6071)(100.0d0*wi(i+nat*(k-1)),i=1,nat) write(6,6071)(100.0d0*wt(i+nat*(k-1)),i=1,nat) 6071 format(12f6.1) c recalculate to get rid of small contributions: write(6 ,28281)k,w write(50,28281)k,w 28281 format(i4,f12.3,' nm, Real weights:') write(6 ,6071)(100.0d0*wr(i+nat*(k-1)),i=1,nat) write(50,6071)(100.0d0*wr(i+nat*(k-1)),i=1,nat) write(6 ,*)'Imaginary weights:' write(50,*)'Imaginary weights:' write(6 ,6071)(100.0d0*wi(i+nat*(k-1)),i=1,nat) write(50,6071)(100.0d0*wi(i+nat*(k-1)),i=1,nat) write(6 ,*)'|ro| weights:' write(50,*)'|ro| weights:' write(6 ,6071)(100.0d0*wt(i+nat*(k-1)),i=1,nat) write(50,6071)(100.0d0*wt(i+nat*(k-1)),i=1,nat) call sort(wr,ind,jnd,nat,k) call trimv(wr,nat,nt,k) write(6 ,*)'Real weights trimmed:' write(50,*)'Real weights trimmed:' write(6 ,6071)(100.0d0*wr(jnd(i)+nat*(k-1)),i=1,nat) write(50,6071)(100.0d0*wr(jnd(i)+nat*(k-1)),i=1,nat) call sort(wi,ind,jnd,nat,k) call trimv(wi,nat,nt,k) write(6 ,*)'Imaginary weights trimmed:' write(50,*)'Imaginary weights trimmed:' write(6 ,6071)(100.0d0*wi(jnd(i)+nat*(k-1)),i=1,nat) write(50,6071)(100.0d0*wi(jnd(i)+nat*(k-1)),i=1,nat) call sort(wt,ind,jnd,nat,k) call trimv(wt,nat,nt,k) write(6 ,*)'|ro| weights trimmed:' write(50,*)'|ro| weights trimmed:' write(6 ,6071)(100.0d0*wt(jnd(i)+nat*(k-1)),i=1,nat) 9082 write(50,6071)(100.0d0*wt(jnd(i)+nat*(k-1)),i=1,nat) close(50) deallocate(ind,jnd) endif deallocate(orb,rrx,rry,rrz,rix,riy,riz,ror,roi) if(lwwe)deallocate(rrxw,rryw,rrzw,rixw,riyw,rizw) close(8) close(9) close(91) write(6,*)' roi.cub ror.cub written' endif c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp stop end subroutine trimv(v,nat,nt,k) implicit none real*8 v(*),t integer*4 nat,i,nt,k,o o=nat*(k-1) t=0.0d0 do 1 nt=1,nat t=t+v(nt+o) 1 if(t.ge.0.5d0)goto 2 2 if(nt.gt.nat)nt=nat do 4 i=1,nt 4 v(i+o)=v(i+o)/t do 3 i=nt+1,nat 3 v(i+o)=0.0d0 return end subroutine sort(v,ind,jnd,nat,k) implicit none real*8 v(*),t integer*4 ind(*),nat,i,j,it,jnd(*),k,o o=nat*(k-1) do 1 i=1,nat jnd(i)=i 1 ind(i)=i do 2 i=1,nat-1 do 2 j=i+1,nat if(v(j+o).gt.v(i+o))then t=v(i+o) v(i+o)=v(j+o) v(j+o)=t it=ind(i) ind(i)=ind(j) ind(j)=it jnd(ind(j))=j jnd(ind(i))=i endif 2 continue return end c subroutine titles(ti,ext,out) character*(*) ti,ext character*80 out do 1 i=1,80 if (ti(i:i).eq.'.'.or.ti(i:i).eq.' ') then k=i-1 goto 2 endif 1 continue 2 out=ti(1:k)//ext return end subroutine setv(eau,w,fr,fi,g,gnm,lgnm,lgau) implicit none real*8 w,jm,fr,fi,g,wau,eau,gnm,a,pi,ln2,ff logical lgnm,lgau c nm to au:= c wau=(1.0d7/w)/219474.63d0 wau=45.5633528d0/w c set bandwidth, either constant in wavelength or energy: if(lgnm)then if(w.gt.gnm)then a=45.5633528d0*gnm/(w-gnm)/(w+gnm) else a=gnm*wau/w endif else a=g endif ff= (eau-wau)*(eau+wau) jm= ff**2+(a*eau)**2 fr= ff / jm if(lgau)then c Gaussian curve for dispersion pi=4.0d0*datan(1.0d0) ln2=log(2.0d0) if(dabs(eau-wau).gt.5.0d0*a)then fi=0.0d0 else fi=dsqrt(pi*ln2)/a/eau**2*exp(-4.0d0*ln2*((eau-wau)/a)**2) endif else c Lorentzian curve G /[(w^2-w0^2)^2 + G^2 w0^2] fi= a / jm endif fr=fr+fr fi=fi+fi return end subroutine getv(n,d,ds1,r0,angle,rs,r,i) implicit none integer*4 i,n real*8 d(n,3),ds1,r0(*),angle,de,dm,ds,r(n,3),rs ds1=(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3))*2.541765d0**2 c ds=2.12689677295232d-30*f(i)/ev(i)/8065.54476345045d0/1.0d-36 rs=r0(i)*1.0d-4 c rs1=(d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3)) c 1*2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36 c c angle between electric and magnetic moment: de=dsqrt(d(i,1)*d(i,1)+d(i,2)*d(i,2)+d(i,3)*d(i,3)) dm=dsqrt(r(i,1)*r(i,1)+r(i,2)*r(i,2)+r(i,3)*r(i,3)) ds= (d(i,1)*r(i,1)+d(i,2)*r(i,2)+d(i,3)*r(i,3)) angle=0.0d0 if(de*dm.gt.0.0d0)angle=acos(ds/de/dm)*180.0d0/4.0d0/atan(1.0d0) return end subroutine report(s) character*(*) s write(6,*)s stop end subroutine rmtr38(A,n0,n,ic) c ic=0 .. reads triangle of a supposedly symmetric matrix c ic=1 .. reads it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n0,n) read(38,*) read(38,*) N1=1 1 N3=MIN(N1+4,N) read(38,*) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=N3 130 read(38,*)J,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 read(38,*) return end subroutine rwm(b,n,nd0) implicit none integer*4 i,j,k,n,n2,nd00,nd0 real*8 b(*) real*8 ,allocatable::p(:,:) parameter (nd00=21) character*2 ss(nd00) character*15 fn data ss/'PX','PY','PZ','VX','VY','VZ','LX','LY','LZ', 1 'DX','DY','DZ','XX','XY','XZ','YY','YZ','ZZ', 1 'EX','EY','EZ'/ c 1 2 3 4 5 6 7 8 9 c 10 11 12 13 14 15 16 17 18 c 19 20 21 logical lex c D .. r x rj c E .. r x ri if(nd0.ne.nd00)call report('nd0 <> nd00') allocate(p(n,n)) n2=n**2 do 1 i=1,nd00 if(i.gt.12.and.i.le.18)then fn='Q'//ss(i)//'.MO.SCR.TXT' else fn=ss(i)//'.MO.SCR.TXT' endif inquire(file=fn,exist=lex) if(lex)then open(38,file=fn) call rmtr38(p,n,n,1) close(38) do 3 j=1,n do 3 k=1,n 3 b(n2*(i-1)+n*(k-1)+j)=p(k,j) write(6,8000)i,fn 8000 format(i4,2x,a15,'read') else write(6,*)ss(i)//'.MO.SCR.TXT not present' do 2 j=1,n do 2 k=1,n 2 b(n2*(i-1)+n*(k-1)+j)=0.0d0 stop endif 1 continue return end subroutine rwmd(b,n) implicit none integer*4 i,j,k,n,n2 real*8 b(*) real*8 ,allocatable::p(:,:) character*3 ss(18) data ss/'PBX','PBY','PBZ','VBX','VBY','VBZ','LBX','LBY','LBZ', 1'WBX','WBY','WBZ','BXX','BXY','BXZ','BYY','BYZ','BZZ'/ logical lex allocate(p(n,n)) n2=n**2 do 1 i=1,18 if(i.gt.12)then inquire(file='Q'//ss(i)//'.MO.SCR.TXT',exist=lex) else inquire(file=ss(i)//'.MO.SCR.TXT',exist=lex) endif if(lex)then if(i.gt.12)then open(38,file='Q'//ss(i)//'.MO.SCR.TXT') else open(38,file=ss(i)//'.MO.SCR.TXT') endif call rmtr38(p,n,n,1) close(38) do 3 j=1,n do 3 k=1,n 3 b(n2*(i-1)+n*(k-1)+j)=p(k,j) else write(6,*)ss(i)//'.MO.SCR.TXT not present' do 2 j=1,n do 2 k=1,n 2 b(n2*(i-1)+n*(k-1)+j)=0.0d0 endif 1 continue return end subroutine getdipole(ut,n,cij,nmo,ni,aij,bij, 1dij,nid,daij,dbij,lopen,pre,pred,ie,nd0,ic,n2) c electronic permanent dipoles of excite state ie implicit none integer*4 ie,j,ni(*),i,oa,ob,aij(*),bij(*),nmo, 1n2,nid(*),daij(*),dbij(*),n,nd0,ic real*8 ut(ic,n),cij(*),dij(*),pre(*),pred(*) logical lopen do 1 j=1,ic c =cie_ab b>: ut(j,ie)=0.0d0 do 2 i=1,ni(ie) oa=aij(n2*(ie-1)+i) ob=bij(n2*(ie-1)+i) 2 ut(j,ie)=ut(j,ie)+cij(n2*(ie-1)+i)*pre(j+(oa-1+(ob-1)*nmo)*nd0) if(lopen)then do 3 i=1,nid(ie) oa=daij(n2*(ie-1)+i) ob=dbij(n2*(ie-1)+i) 3 ut(j,ie)=ut(j,ie)+dij(n2*(ie-1)+i)*pred(j+(oa-1+(ob-1)*nmo)*nd0) endif 1 continue return end subroutine getieab(pre,pred,lopen,na,nb,nmo,ni,nid,ie,aij,bij, 1cij,daij,dbij,dij,bp,bpd,u0,u0d,nmo2,nd0,n22) c b> c note that for singlets states the matrix elements are the same as c for Slater's determinants implicit none logical lopen integer oa,ob,oc,od,na,nb,nmo,j,i,ni(*),nid(*),ie,nmo2, 1aij(*),bij(*),daij(*),dbij(*),nd0,jj,n22 real*8 pre(*),pred(*),cij(*),dij(*),ciji,bp(*),bpd(*), 1u0(*),u0d(*) do 1 oc=1,nd0*nmo*nmo pre(oc)=0.0d0 1 pred(oc)=0.0d0 if(lopen)then do 2 i=1,ni(ie) oa =aij(n22*(ie-1)+i) ob =bij(n22*(ie-1)+i) ciji=cij(n22*(ie-1)+i) do 2 j=1,nd0 jj=nmo2*(j-1) c (ABCD) .. note that A<>B and C<>D: c if(oc.eq.oa) c if(od.eq.ob) : pre(j+(oa-1+(ob-1)*nmo)*nd0)=pre(j+(oa-1+(ob-1)*nmo)*nd0) 1 +ciji*(u0(j)+u0d(j)-bp(oa+nmo*(oa-1)+jj) 1 + bp(ob+nmo*(ob-1)+jj)) c if(od.ne.ob)=: do 3 od=na+1,nmo if(ob.ne.od) 1 pre(j+(oa-1+(od-1)*nmo)*nd0)= 2 pre(j+(oa-1+(od-1)*nmo)*nd0)+ciji*bp(ob+nmo*(od-1)+jj) 3 continue c if(oc.ne.oa) =: do 2 oc=1,na 2 if(oc.ne.oa) 1 pre(j+(oc-1+(ob-1)*nmo)*nd0)= 2 pre(j+(oc-1+(ob-1)*nmo)*nd0)+ciji*bp(oa+nmo*(oc-1)+jj) do 42 i=1,nid(ie) oa =daij(n22*(ie-1)+i) ob =dbij(n22*(ie-1)+i) ciji=dij(n22*(ie-1)+i) c if(oa.ne.oc)then if(ob.eq.od)then: do 4 oc=1,nb if(oc.ne.oa)then do 22 j=1,nd0 22 pred(j+(oc-1+(ob-1)*nmo)*nd0)= 1 pred(j+(oc-1+(ob-1)*nmo)*nd0)+ciji*bpd(oa+nmo*(oc-1)+nmo2*(j-1)) endif 4 continue c if(oa.eq.oc)then if(ob.ne.od)then: do 5 od=nb+1,nmo if(ob.ne.od)then do 32 j=1,nd0 32 pred(j+(oa-1+(od-1)*nmo)*nd0)= 1 pred(j+(oa-1+(od-1)*nmo)*nd0)+ciji*bpd(ob+nmo*(od-1)+nmo2*(j-1)) endif 5 continue c if(oa.eq.oc)then if(ob.eq.od)then: do 42 j=1,nd0 42 pred(j+(oa-1+(ob-1)*nmo)*nd0)=pred(j+(oa-1+(ob-1)*nmo)*nd0) 1 +ciji*(u0(j)+u0d(j)-bpd(oa+nmo*(oa-1)+nmo2*(j-1)) 1 + bpd(ob+nmo*(ob-1)+nmo2*(j-1))) else c closed shell- note that rule is different than for open shell: do 6 i=1,ni(ie) oa =aij(n22*(ie-1)+i) ob =bij(n22*(ie-1)+i) ciji=cij(n22*(ie-1)+i) do 6 j=1,nd0 jj=nmo2*(j-1) c (ABCD) .. note that A<>B and C<>D: c if(oc.eq.oa) c if(od.eq.ob) : pre(j+(oa-1+(ob-1)*nmo)*nd0)=pre(j+(oa-1+(ob-1)*nmo)*nd0) 1 +ciji*(u0(j)+u0d(j)-bp(oa+nmo*(oa-1)+jj) 1 + bp(ob+nmo*(ob-1)+jj)) c if(od.ne.ob)=: do 6 od=na+1,nmo 6 if(ob.ne.od) 1 pre(j+(oa-1+(od-1)*nmo)*nd0)= 2 pre(j+(oa-1+(od-1)*nmo)*nd0)+ciji*bp(ob+nmo*(od-1)+jj) endif return end subroutine readoptions(lden,lmcd,lrds,lw,wden,fo,icpl, 1lzmat,ldeg,lwrt,DD,degper0,deglim,lnorm,nst,ifix,nmo,igiao, 1gammaau,wmin,wmax,np,lgamma,lort,gammanm,lgnm,lgau,lwwe,nmms, 1mcdv,nroot,fov,LEXCL,LQ1,np0,iwr,THRC,LDD,lgnj,LDE,rgnj,lbck) implicit none character*(*) fo,fov logical lden,lmcd,lrds,lw,lzmat,ldeg,lwrt,lnorm,lgamma,lort, 1lgnm,lgau,lwwe,mcdv,LDD,LQ1,lgnj,LDE,rgnj,lbck real*8 DD,degper0,deglim,wden,wmin,wmax,gammaau,gammanm,THRC character*4 key integer*4 nst,ifix,nmo,igiao,np,icpl,nmms,nroot,iwr,np0,LEXCL c Vibrational structure options: c read dipole derivatives from DE.TEN LDE=.false. c read derivatives of the Gnj tensor from GNJD.TEN: lgnj=.false. c read Gnj tensor from GNJR.TEN (debug option): rgnj=.false. c invoke the calculation: mcdv=.false. c vibrational parameter file: fov='GV.OUT' c maximal degree of excitations ("class"): LEXCL=3 c use first dipole derivatives: LQ1=.true. c dimension of the temporary list: np0=10 c optional writing option: iwr=1 c threshold to consider Franck-Condon factors: THRC=1.0d-9 c read dipole derivatives: LDD=.true. c number of th eexcited state nroot=1 c read also backward excitations: lbck=.false. c Other options: degper0=0.0001d0 deglim=0.00001d0 lden=.false. nmms=3 ifix=3 nst=0 lnorm=.false. lort=.true. lgnm=.true. lgamma=.true. lwwe=.false. lgau=.false. lmcd=.true. icpl=0 lrds=.false. lw=.false. wden=532.0d0 fo='G.OUT' lzmat=.true. ldeg=.false. lwrt=.false. DD=2.0d0 nmo=0 igiao=0 gammaau=0.02d0 gammanm=10.0d0 wmin=50.0d0 wmax=500.0d0 np=901 open(9,file='GUVCDE.OPT') 1001 read(9,90,err=109,end=109)key 90 format(a4) c gaussian TDDFT output: if(key.eq.'FILE')read(9,9010)fo 9010 format(a80) c GIAO experimenting: c igiao = 0 (r-Rb) c igiao = 1 (r-R(a+b)) c igiao = 2 (r-R( b))*Dab c igiao = 3 (r-R(a+b))*Dab if(key.eq.'GIAO')read(9,*)igiao c use z-matrix orientation: if(key.eq.'ZMAT')read(9,*)lzmat c remove degeneracy if present: if(key.eq.'LDEG')read(9,*)ldeg c frequency to produce .cub file if(key.eq.'WDEN')read(9,*)wden c produce .cub file: if(key.eq.'LDEN')read(9,*)lden c produce atomic weights: if(key.eq.'LWEI')read(9,*)lw c atomix extent for the weigths in A: if(key.eq.'EXTE')read(9,*)DD c calculate MCD: if(key.eq.'LMCD')read(9,*)lmcd c calculate circular polarized luminiscence from this state if(key.eq.'ICPL')read(9,*)icpl c recalculate dipole strengths: if(key.eq.'LRDS')read(9,*)lrds c degeneracy limit detection: if(key.eq.'DEGL')read(9,*)deglim c degeneracy perturbation: if(key.eq.'DEGP')read(9,*)degper0 c how many E G controbutions for MCD statistics: if(key.eq.'NMMS')read(9,*)nmms c extended writing if(key.eq.'LWRT')read(9,*)lwrt c calculate weights for all NP points between WMIN and WMAX: if(key.eq.'LWEE')read(9,*)lwwe c renormalize cij coefficients if(key.eq.'LNOR')read(9,*)lnorm c renorthonormalize cij coefficients if(key.eq.'LORT')read(9,*)lort c limit number of excited states (put zero to take all): if(key.eq.'NSTA')read(9,*)nst if(key(1:3).eq.'NMO')read(9,*)nmo c fix for backward excitations c 0 .. none c 1 .. add to cij c 2 .. consider separately c 3 .. ad as extra transition c 4 .. ignore, probably same as 0 c 5 .. add complex (not implemented) if(key.eq.'IFIX')read(9,*)ifix c bandwidth in au, for computation of MCD and polarizabilities: if(key.eq.'GAMM')read(9,*)gammaau c bandwidth in nm, for computation of polarizabilities if required c by lgnm: if(key.eq.'GANM')read(9,*)gammanm c whether to use gammaau in MCD: if(key.eq.'LGAM')read(9,*)lgamma c whether to use Gaussian dispersion curves: if(key.eq.'LGAU')read(9,*)lgau c whether to use lgnm: if(key.eq.'LGNM')read(9,*)lgnm c minimal and maximal frequency (nm) c for frequecny-dependent polarizabilities: if(key.eq.'WMIN')read(9,*)wmin if(key.eq.'WMAX')read(9,*)wmax c number of points in the above frequency interval: if(key(1:3).eq.'NPO)')read(9,*)np c make the vibrationally resolved calculation: if(key(1:4).eq.'MCDI')read(9,*)mcdv c number of the excited state for MCDV: if(key(1:4).eq.'NROO')read(9,*)nroot c filename with the parameters for MCDV: if(key(1:4).eq.'FILV')read(9,*)fov c maximum number of excitations: if(key(1:4).eq.'LEXC')read(9,*)LEXCL c use fist dipole derivatives: if(key(1:3).eq.'LQ1')read(9,*)LQ1 c dimension of the expansion space for Franck Condon factors: if(key(1:3).eq.'NP0')read(9,*)np0 c additional writing option for MCDV: if(key(1:3).eq.'IWR')read(9,*)iwr c threshold for vibrational contributions: if(key(1:4).eq.'THRC')read(9,*)THRC c read dipole derivatives: if(key(1:3).eq.'LDD')read(9,*)LDD c read dipole derivatives from DE.TEN: if(key(1:4).eq.'LDEE')read(9,*)LDE c read Gnj tensor derivatives from GNJD.TEN: if(key(1:4).eq.'LGNJ')read(9,*)lgnj c read also backward excitations if(key(1:4).eq.'LBCK')read(9,*)lbck c read Gnj tensor from GNJR.TEN: if(key(1:4).eq.'RGNJ')read(9,*)rgnj goto 1001 109 close(9) write(6,6070)fo ,lzmat,wden,lden,lw,DD,lmcd,lrds,deglim, 1degper0,lwrt,lnorm,lort,nst,ifix,igiao,ldeg,nmo, 1lgnm,lgamma,np,wmin,wmax,gammaau,gammanm,lgau,icpl 6070 format( 1'FILE',/,a80 ,/,'ZMAT ',/,l1 ,/,'WDEN (cm-1)',/,f10.2,/, 1'LDEN',/,l1 ,/,'LWEIGHT',/,l1 ,/,'EXTEND ',/,f10.2,/, 1'LMCD',/,l1 ,/,'LRDS ',/,l1 ,/,'DEGLIM ',/,f10.4,/, 1'DEGP',/,f10.4,/,'LWRT ',/,l1 ,/,'LNORM ',/,l1 ,/, 1'LORT',/,l1 ,/,'NSTA ',/,I10,/,'IFIX ',/,I1 ,/, 1'GIAO',/,i1 ,/,'LDEG ',/,l1 ,/,'NMO ',/,i10 ,/, 1'LGNM',/,l1 ,/,'LGAM ',/,l1 ,/,'NPLGAM' ,/,i10 ,/, 1'WMIN' ,/,f10.2,' nm' ,/,'WMAX',/,f10.2,' nm',/, 1'GAMMAAU',/,f10.5,' hartree',/,'GANM',/,f10.2,' nm',/, 1'LGAUSS',/,l1,/,'ICPL',/,i10) return end subroutine dimensions(n,fo,nmo,nmo2,nat,lzmat,lturbo,ltda,ncmm) c returns number of transitions and atoms implicit none integer*4 n,nat,l,I,q,nmo,nmo2,ncmm,it,nt real*8 axdum character*80 s80 character*70 st character*(*) fo logical lzmat,lturbo,le,ltda n=0 nat=0 ncmm=0 it=0 c Grimme's TDA output - (CI soefficients from cicond) inquire(file='ciss_a.co',exist=ltda) if(ltda)then open(2,file='ciss_a.co') n=0 10 read(2,2000,end=790,err=790)s80 if(s80(1:1).eq.'e'.or.s80(1:2).eq.' ')then n=n+1 c nt: number of expansion coefficients for this state nt=0 11 read(2,2000,end=690,err=690)s80 if(s80(1:1).ne.'e'.and.s80(1:2).ne.' ')then nt=nt+1 goto 11 else backspace 2 endif 690 if(nt.gt.ncmm)then ncmm=nt it=n endif endif goto 10 790 close(2) write(6,*)n,' transitions in ciss_a.co' write(6,*)ncmm,' coefficients for state ',it endif c Turbomole output: inquire(file='control',exist=lturbo) if(lturbo)then write(6,*)'control file from Turbomole found ....' open(2,file='control') nmo=0 100 read(2,2000,end=990,err=990)s80 IF(s80(4:10).EQ.'natoms=')read(s80(11:len(s80)),*)nat IF(s80(4:11).EQ.'nbf(AO)=')read(s80(12:len(s80)),*)nmo if(nat.eq.0.or.nmo.eq.0)goto 100 990 close(2) if(nat.eq.0.or.nmo.eq.0)call report('nat nmo not determined') nmo2=nmo**2 write(6,*)' ..... and read' inquire(file='escf.out',exist=le) if(le)then write(6,*)'escf.out found' open(2,file='escf.out') 200 read(2,2000,end=991,err=991)s80 IF(s80(2:40).EQ.'total number of roots to be determined:')then read(s80(41:len(s80)),*)n goto 991 endif goto 200 991 close(2) if(n.eq.0)call report('n not determined') else call report('escf.out not found') endif return endif c Gaussian output: open(2,file=fo) 1 read(2,2000,end=99,err=99)s80 2000 format(a80) IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(20:37).EQ.'Input orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard orientation:')))THEN write(6,2000)s80 DO 2004 I=1,4 2004 READ(2,*) l=0 2005 READ(2,2000)s80 IF(s80(2:4).NE.'---')THEN l=l+1 BACKSPACE 2 READ(2,*)q,q IF(q.EQ.-1)l=l-1 GOTO 2005 ENDIF nat=l ENDIF if(s80(31:49).eq.'primitive gaussians'.and.nmo.eq.0)then read(s80(1:6),*)nmo endif if(s80(2:44).eq.'Ground to excited state transition electric')then n=0 read(2,*) 3 read(2,3000,end=99,err=99)st 3000 format(10x,a70) read(st,*,end=99,err=99)axdum,axdum,axdum n=n+1 goto 3 endif goto 1 99 close(2) nmo2=nmo**2 return end subroutine readfile(lzmat,qz,x,bohr,xau,na,nb,nmo2,n,nmo,nst, 1dl,dv,qv,r,r0v,r0l,e,ev,eau,ni,lden,nid,nat,cij,lopen, 1aij,bij,dij,daij,dbij,fo,ifix,nib,cijb,aijb,bijb,lnorm,lort,lwrt, 1dijb,daijb,dbijb,nibd,ltda,lbck) implicit none character*80 s80 character*75 st character*(*) fo logical lzmat,lden,lopen,lnorm,lort,lwrt,ltda,lbck integer*4 ig98,I,l,qz(*),na,nb,nmo2,n,nd,k,ie,ni(*),ic,kk, 1nid(*),nat,aij(*),bij(*),dbij(*),daij(*),j,oa,ob,jj,nst, 1ifound,ifix,nib(*),aijb(*),bijb(*),amax,bmax,nmo,nort, 1daijb(*),dbijb(*),nibd(*),imax,ii real*8 x(3,nat),bohr,xau(3,nat),dl(n,3),dv(n,3),qv(n,3,3), 1r(n,3),r0v(*),r0l(*),e(*),ev(*),eau(*),cij(*),dij(*),cback,tq, 1ax,ay,az,cijb(*),cmax,skj,sum1,sum2,dijb(*) real*8,allocatable::obaj(:),obak(:),obajd(:),obakd(:) real*8 ska integer*4,allocatable::ipse(:) nat=0 ie=0 lopen=.false. tq=dsqrt(2.0d0) open(2,file=fo) 1 read(2,2000,end=2,err=2)s80 2000 format(a80) IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(20:37).EQ.'Input orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s80(20:40).EQ.'Standard orientation:'.OR. 2 s80(26:46).EQ.'Standard orientation:')))THEN ig98=0 if(s80(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s80(27:44).EQ.'Input orientation:'.OR. 1 s80(26:46).EQ.'Standard orientation:')ig98=1 write(6,2000)s80 DO 2004 I=1,4 2004 READ(2,*) l=0 2005 READ(2,2000)s80 IF(s80(2:4).NE.'---')THEN l=l+1 BACKSPACE 2 if(ig98.eq.0)then READ(2,*)qz(l),qz(l),(x(i,l),i=1,3) else READ(2,*)qz(l),qz(l),x(1,l),(x(i,l),i=1,3) endif xau(1,l)=x(1,l)/bohr xau(2,l)=x(2,l)/bohr xau(3,l)=x(3,l)/bohr IF(qz(l).EQ.-1)l=l-1 GOTO 2005 ENDIF nat=l ENDIF if(s80(40:50).eq.'Pseudopoten')then allocate(ipse(nat)) call rdps(nat,ipse) do 335 i=1,nat if(ipse(i).ne.0)then write(6,6885)i,ipse(i) 6885 format(' atom ',i4,' atomic charge reduced to',i4) qz(i)=ipse(i) endif 335 continue endif if(s80(8:22).eq.'alpha electrons')then read(s80( 1: 6),*)na read(s80(25:31),*)nb write(6,2000)s80 if(ltda)then write(6,*)'leaving readfile because of tda' close(2) return endif endif if(s80(2:44).eq.'Ground to excited state transition electric')then n=0 write(6,2000)s80 read(2,*) 3 read(2,3000,end=21,err=21)st 3000 format(10x,a75) read(st,*,end=21,err=21)ax,ay,az n=n+1 dl(n,1)=ax dl(n,2)=ay dl(n,3)=az goto 3 21 backspace 2 write(6,*)n,' dipole length transitions' endif if(s80(2:51).eq. 1'Ground to excited state transition velocity dipole')then write(6,2000)s80 read(2,*) do 9 i=1,n read(2,3000)st 9 read(st,*)dv(i,1),dv(i,2),dv(i,3) write(6,*)n,' dipole velocity transitions' endif c if(s80(2:44).eq.'Ground to excited state transition magnetic')then write(6,2000)s80 read(2,*) do 10 i=1,n read(2,3000)st 10 read(st,*)(r(i,j),j=1,3) write(6,*)n,' magnetic transitions' endif if(s80(2:55).eq. 1'Ground to excited state transition velocity quadrupole')then write(6,2000)s80 read(2,*) do 7 i=1,n read(2,3000)st read(st,*)qv(i,1,1),qv(i,2,2),qv(i,3,3),qv(i,1,2),qv(i,1,3), 1 qv(i,2,3) qv(i,3,2)=qv(i,2,3) qv(i,3,1)=qv(i,1,3) 7 qv(i,2,1)=qv(i,1,2) write(6,*)n,' quadrupole velocity transitions' endif c if(s80(53:63).eq.'R(velocity)') then write(6,2000)s80 c Gaussian 98 does not write R(length) nor R(velocity) do 8 i=1,n read(2,3000)st read(st,*,err=99)r0v(i),r0v(i),r0v(i),r0v(i) goto 991 99 write(6,3000)st write(6,*)i stop 991 continue read(2,2000)s80 if(s80(2:6).eq.'Total')then do 81 k=1,4 81 read(2,2000)s80 else backspace 2 endif read(2,2000)s80 if(s80(2:3).eq.'R(')then do 82 k=1,4 82 read(2,2000)s80 else backspace 2 endif 8 continue endif if(s80(54:62).eq.'R(length)') then write(6,2000)s80 do 11 i=1,n read(2,3000)st 11 read(st,*,end=2,err=2)r0l(i),r0l(i),r0l(i),r0l(i) endif c if(s80(2:15).eq.'Excited State '.and.s80(16:16).ne.'s')then write(6,2000)s80 do 203 i=1,79 if(s80(i:i).eq.':')read(s80(15:i-1),*)nd if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)ev(nd) 203 if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)e(nd) c ecm(i)=ev(i)*8065.54476345045d0 eau(nd)=ev(nd)/27.211384205943d0 ie=ie+1 if(ie.ne.nd)call report(' Inconsistent reading') cmax=0.0d0 amax=0 bmax=0 imax=0 if(lden)then ni(ie)=0 nib(ie)=0 nibd(ie)=0 nid(ie)=0 sum1=0.0d0 sum2=0.0d0 111 read(2,2000,end=299,err=299)s80 if(s80(10:11).eq.'->')then if(ni(ie).gt.nmo2)call report('too many coefficients') if(s80(8:8).eq.'A')then c open shell, alpha: lopen=.true. ni(ie)=ni(ie)+1 ii=nmo2*(ie-1)+ni(ie) read(s80( 1: 7),*)aij(ii) if(s80(15:15).eq.'A')then read(s80(12:14),*)bij(ii) else if(s80(16:16).eq.'A')then read(s80(12:15),*)bij(ii) else write(6,2000)s80 call report('Unknown format of excitations') endif endif read(s80(18:30),*)cij(ii) call setabmax(aij(ii),bij(ii),cij(ii), 2 amax,bmax,cmax,ii,imax) sum1=sum1+cij(ii)**2 else if(s80(8:8).eq.'B')then c open shell, beta: nid(ie)=nid(ie)+1 ii=nmo2*(ie-1)+nid(ie) lopen=.true. read(s80( 1: 7),*)daij(ii) if(s80(15:15).eq.'B')then read(s80(12:14),*)dbij(ii) else if(s80(16:16).eq.'B')then read(s80(12:15),*)dbij(ii) else write(6,2000)s80 call report('Unknown format of excitations') endif endif read(s80(18:30),*) dij(ii) call setabmax(daij(ii),dbij(ii),dij(ii), 2 amax,bmax,cmax,ii,imax) sum1=sum1+dij(ii)**2 else c closed shell: ni(ie)=ni(ie)+1 ii=nmo2*(ie-1)+ni(ie) read(s80( 1: 8),*) aij(ii) read(s80(12:15),*) bij(ii) read(s80(18:30),*) cij(ii) cij(ii)=cij(ii)*tq call setabmax(aij(ii),bij(ii),cij(ii), 2 amax,bmax,cmax,ii,imax) sum1=sum1+cij(ii)**2 endif endif goto 111 endif if(s80(10:11).eq.'<-'.and.lbck)then if(s80(8:8).eq.'A')then nib(ie)=nib(ie)+1 read(s80( 1: 7),*)oa if(s80(15:15).eq.'A')then read(s80(12:14),*)ob else if(s80(16:16).eq.'A')then read(s80(12:15),*)ob else write(6,2000)s80 call report('Unknown format of excitations') endif endif read(s80(18:30),*)cback sum2=sum2+cback**2 cijb(nmo2*(ie-1)+nib(ie))=cback aijb(nmo2*(ie-1)+nib(ie))=oa bijb(nmo2*(ie-1)+nib(ie))=ob else if(s80(8:8).eq.'B')then nibd(ie)=nibd(ie)+1 read(s80( 1: 7),*)oa if(s80(15:15).eq.'B')then read(s80(12:14),*)ob else if(s80(16:16).eq.'B')then read(s80(12:15),*)ob else write(6,2000)s80 call report('Unknown format of excitations') endif endif read(s80(18:30),*)cback sum2=sum2+cback**2 dijb( nmo2*(ie-1)+nibd(ie))=cback daijb(nmo2*(ie-1)+nibd(ie))=oa dbijb(nmo2*(ie-1)+nibd(ie))=ob else nib(ie)=nib(ie)+1 read(s80( 1: 8),*)oa read(s80(12:15),*)ob read(s80(18:30),*)cback cback=cback*tq sum2=sum2+cback**2 cijb(nmo2*(ie-1)+nib(ie))=cback aijb(nmo2*(ie-1)+nib(ie))=oa bijb(nmo2*(ie-1)+nib(ie))=ob endif endif c experimental fix for the deexitations, closed shell only: if(ifix.eq.1.and..not.lopen)then ifound=0 do 315 j=1,ni(ie) if(oa.eq.aij(nmo2*(ie-1)+j).and.ob.eq.bij(nmo2*(ie-1)+j)) 1 then ifound=ifound+1 cij(nmo2*(ie-1)+j)=cij(nmo2*(ie-1)+j)+cback endif 315 continue if(ifound.eq.0)then ni(ie)=ni(ie)+1 cij(nmo2*(ie-1)+ni(ie))=cback aij(nmo2*(ie-1)+ni(ie))=oa bij(nmo2*(ie-1)+ni(ie))=ob endif endif c experimental fix,just add to excitations with switched orbitals: if(ifix.eq.3)then if(s80(8:8).ne.'B')then ni(ie)=ni(ie)+1 cij(nmo2*(ie-1)+ni(ie))=cback aij(nmo2*(ie-1)+ni(ie))=ob bij(nmo2*(ie-1)+ni(ie))=oa else nid(ie)=nid(ie)+1 dij(nmo2*(ie-1)+ni(ie))=cback daij(nmo2*(ie-1)+ni(ie))=ob dbij(nmo2*(ie-1)+ni(ie))=oa endif endif goto 111 endif 299 continue if(lwrt)then write(6,*)ni(ie),' alpha coefficients' if(ifix.eq.3)write(6,*)'(with backwards)' write(6,*)nid(ie),' beta coefficients' write(6,*)nib(ie),' backward coefficients' if(s80(8:8).eq.'B')then write(6,6009)ie,amax,bmax,dij(imax),cmax**2*100.0d0 6009 format(i5,': Dominant transition from',i4,' to ',i4,f8.3, 1 ' (',f6.2,' %)') else write(6,6009)ie,amax,bmax,cij(imax),cmax**2*100.0d0 endif write(6,1234)sum1,sum2,sum1-sum2 1234 format('fnorm, bnorm,diff:',3f10.4) endif if(ifix.eq.4.and.lwrt)then nib(ie)=0 write(6,*)'backward ignored' endif c seems that fnorm-bnorm=1 if(lnorm)then c renormalize: call normcijo(ie,cij,dij,ni,nid,nmo2,lopen) write(6,*)' cij were normalized' endif endif if(nd.eq.n) goto 2 endif goto 1 c 2 close(2) if(lort)then if(nst.gt.0.and.nst.lt.n)then nort=nst else nort=n endif c reorthonormalize: write(6,*)'orthogonalizing coefficients for ',nort,' states' allocate(obaj(nmo2),obak(nmo2)) if(lopen)allocate(obajd(nmo2),obakd(nmo2)) write(6,*)'temporary variables allocated' c run twice for better accuracy: do 318 ic=1,2 ska=0.0d0 do 3181 j=1,nort c transcript original cij to obaj: call vz(obaj,nmo2) do 324 jj=1,ni(j) oa=aij(nmo2*(j-1)+jj) ob=bij(nmo2*(j-1)+jj) 324 obaj(ob+nmo*(oa-1))=cij(nmo2*(j-1)+jj) if(lopen)then call vz(obajd,nmo2) do 3241 jj=1,nid(j) oa=daij(nmo2*(j-1)+jj) ob=dbij(nmo2*(j-1)+jj) 3241 obajd(ob+nmo*(oa-1))=dij(nmo2*(j-1)+jj) endif do 334 k=1,j-1 c transcript cij to obak: call vz(obak,nmo2) kk=nmo2*(k-1) do 3242 jj=1,ni(k) oa=aij(kk+jj) ob=bij(kk+jj) 3242 obak(ob+nmo*(oa-1))=cij(kk+jj) if(lopen)then call vz(obakd,nmo2) kk=nmo2*(k-1) do 3243 jj=1,nid(k) oa=daij(kk+jj) ob=dbij(kk+jj) 3243 obakd(ob+nmo*(oa-1))=dij(kk+jj) endif c : skj=0.0d0 do 323 oa=1,na do 323 ob=na+1,nmo c forward excitations: skj=skj+obaj(ob+nmo*(oa-1))*obak(ob+nmo*(oa-1)) c backward excitations: 323 skj=skj+obaj(oa+nmo*(ob-1))*obak(oa+nmo*(ob-1)) if(lopen)then do 331 oa=1,nb do 331 ob=nb+1,nmo skj=skj+obajd(ob+nmo*(oa-1))*obakd(ob+nmo*(oa-1)) 331 skj=skj+obajd(oa+nmo*(ob-1))*obakd(oa+nmo*(ob-1)) endif c |j'>=|j>-|k>: do 319 oa=1,na do 319 ob=na+1,nmo if(oa+nmo*(ob-1).gt.nmo2)write(6,*)'A',oa,ob,nmo2 if(ob+nmo*(oa-1).gt.nmo2)write(6,*)'B',oa,ob,nmo2 obaj(ob+nmo*(oa-1))=obaj(ob+nmo*(oa-1))-skj*obak(ob+nmo*(oa-1)) 319 obaj(oa+nmo*(ob-1))=obaj(oa+nmo*(ob-1))-skj*obak(oa+nmo*(ob-1)) if(lopen)then do 332 oa=1,nb do 332 ob=nb+1,nmo obajd(ob+nmo*(oa-1))= 1 obajd(ob+nmo*(oa-1))-skj*obakd(ob+nmo*(oa-1)) 332 obajd(oa+nmo*(ob-1))= 1 obajd(oa+nmo*(ob-1))-skj*obakd(oa+nmo*(ob-1)) endif 334 ska=ska+dabs(skj) c c transcript new obaj to new cij: if(lopen)then nid(j)=0 do 333 oa=1,nb do 333 ob=nb+1,nmo jj=ob+nmo*(oa-1) if(obajd(jj).ne.0.0d0)then nid(j)=nid(j)+1 dij( nmo2*(j-1)+nid(j))=obajd(jj) daij(nmo2*(j-1)+nid(j))=oa dbij(nmo2*(j-1)+nid(j))=ob endif jj=oa+nmo*(ob-1) if(obajd(jj).ne.0.0d0)then nid(j)=nid(j)+1 dij( nmo2*(j-1)+nid(j))=obajd(jj) daij(nmo2*(j-1)+nid(j))=ob dbij(nmo2*(j-1)+nid(j))=oa endif 333 continue endif ni(j)=0 do 325 oa=1,na do 325 ob=na+1,nmo jj=ob+nmo*(oa-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(nmo2*(j-1)+ni(j))=obaj(jj) aij(nmo2*(j-1)+ni(j))=oa bij(nmo2*(j-1)+ni(j))=ob endif jj=oa+nmo*(ob-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(nmo2*(j-1)+ni(j))=obaj(jj) aij(nmo2*(j-1)+ni(j))=ob bij(nmo2*(j-1)+ni(j))=oa endif 325 continue 3181 call normcijo(j,cij,dij,ni,nid,nmo2,lopen) write(6,*)ic,ska/dble(nort**2)/2.0d0 318 continue write(6,*)'double orthogonalization done' endif return end subroutine normcijo(ie,cij,dij,ni,nid,nmo2,lopen) implicit none real*8 anorm,cij(*),dij(*) integer*4 j,ni(*),ie,nmo2,nid(*) logical lopen anorm=0.0d0 do 315 j=1,ni(ie) 315 anorm=anorm+cij(nmo2*(ie-1)+j)**2 if(lopen)then do 316 j=1,nid(ie) 316 anorm=anorm+dij(nmo2*(ie-1)+j)**2 endif anorm=1.0d0/dsqrt(anorm) do 317 j=1,ni(ie) 317 cij(nmo2*(ie-1)+j)=cij(nmo2*(ie-1)+j)*anorm if(lopen)then do 318 j=1,nid(ie) 318 dij(nmo2*(ie-1)+j)=dij(nmo2*(ie-1)+j)*anorm endif return end subroutine normcij(ie,cij,ni,nmo2) implicit none real*8 anorm,cij(*) integer*4 j,ni(*),ie,nmo2 anorm=0.0d0 do 316 j=1,ni(ie) 316 anorm=anorm+cij(nmo2*(ie-1)+j)**2 anorm=1.0d0/dsqrt(anorm) do 317 j=1,ni(ie) 317 cij(nmo2*(ie-1)+j)=cij(nmo2*(ie-1)+j)*anorm return end subroutine me( 1ni,nid,cij,dij,aij,daij,bij,dbij,ie,lopen,djg,bp,bpd,nmo,nmo2, 1ifix,nib,cijb,aijb,bijb,nd0,n22) c element implicit none integer*4 ni(*),nid(*),j,oa,ob,aij(*),bij(*),daij(*),dbij(*), 1ie,i,nmo,nmo2,nib(*),aijb(*),bijb(*),ifix,nd0,n22 real*8 djg(*),cij(*),dij(*),bp(*),bpd(*),tq,cijb(*) logical lopen tq=dsqrt(2.0d0) do 303 j=1,nd0 djg(j)=0.0d0 do 301 i=1,ni(ie) oa=aij(n22*(ie-1)+i) ob=bij(n22*(ie-1)+i) 301 djg(j)=djg(j)+cij(n22*(ie-1)+i)*bp(ob+nmo*(oa-1)+nmo2*(j-1)) if(ifix.eq.2)then do 302 i=1,nib(ie) oa=aijb(n22*(ie-1)+i) ob=bijb(n22*(ie-1)+i) 302 djg(j)=djg(j)+cijb(n22*(ie-1)+i)*bp(oa+nmo*(ob-1)+nmo2*(j-1)) endif do 312 i=1,nid(ie) oa=daij(n22*(ie-1)+i) ob=dbij(n22*(ie-1)+i) 312 djg(j)=djg(j)+dij(n22*(ie-1)+i)*bpd(ob+nmo*(oa-1)+nmo2*(j-1)) c c for close shells, adjust for singlet-transition: 303 if(.not.lopen)djg(j)=djg(j)*tq return end function cc(xy,n,xa,ya,xx,yy) implicit none real*8 xy,n,xa,ya,xx,yy,cc,d d=(xx-n*xa**2)*(yy-n*ya**2) if(d.eq.0.0d0)then cc=0.0d0 else cc=(xy-n*xa*ya)**2/d endif return end subroutine cstat(v,r,xy,xa,ya,ar,xx,yy) implicit none real*8 v,r,xy,xa,ya,ar,xx,yy xy=xy+v*r xx=xx+v*v yy=yy+r*r xa=xa+v ya=ya+r if(dabs(r).gt.1.0d-4)then ar=ar+v/r else ar=ar+1.0d0 endif return end subroutine stn(n,ar,xa,ya,reg,yy,xy) implicit none integer*4 n real*8 ar,xa,ya,reg,yy,xy,n3 n3=dble(n) ar=ar/n3 xa=xa/n3 ya=ya/n3 reg=1.0d0 if(yy.ne.0.0d0)reg=xy/yy return end subroutine zerstat(xa,ya,xx,yy,xy,ar,n) implicit none real*8 xa(*),ya(*),xx(*),yy(*),xy(*),ar(*) integer*4 i,n do 107 i=1,n xa(i)=0.0d0 ya(i)=0.0d0 xx(i)=0.0d0 yy(i)=0.0d0 xy(i)=0.0d0 107 ar(i)=0.0d0 return end subroutine setabmax(a,b,c,amax,bmax,cmax,n,i) implicit none integer*4 a,b,amax,bmax,n,i real*8 c,cmax if(dabs(c).gt.cmax)then cmax=dabs(c) amax=a bmax=b i=n endif return end subroutine dodjk(djk,nd0,ni,iep,aij,bij,n22,cij,pre,lopen, 1daij,dbij,pred,nmo,dij,nid) implicit none real*8 djk(*),cij(*),pre(*),pred(*),dij(*) logical lopen integer*4 j,nd0,ip,ni(*),iep,oc,od,aij(*),bij(*),n22, 1daij(*),dbij(*),nid(*),nmo do 309 j=1,nd0 djk(j)=0.0d0 do 308 ip=1,ni(iep) oc=aij(n22*(iep-1)+ip) od=bij(n22*(iep-1)+ip) 308 djk(j)=djk(j)+cij(n22*(iep-1)+ip)*pre(j+(oc-1+(od-1)*nmo)*nd0) if(lopen)then do 3071 ip=1,nid(iep) oc=daij(n22*(iep-1)+ip) od=dbij(n22*(iep-1)+ip) 3071 djk(j)=djk(j) 1 +dij(n22*(iep-1)+ip)*pred(j+(oc-1+(od-1)*nmo)*nd0) endif 309 continue return end subroutine lsums(dx,dy,dz,rx,ry,rz,qx,qy,qz,hx,hy,hz, 1n,eau,lgamma,gammaau,ie,iep,dd,dlg,na,nb, 1rrx,rry,rrz,qqx,qqy,qqz,u0,u0d,ut) implicit none real*8 dx,dy,dz,rx,ry,rz,qx,qy,qz,hx,hy,hz,elg,ekl, 1eau(*),gammaau,dd(*),dlg(*),flg,fkl,fne, 2rrx,rry,rrz,qqx,qqy,qqz,ejl,fjl,ejg,u0(*),u0d(*), 3ut(3,n),dg1,dg2,dg3,dk1,dk2,dk3, 4ekg,fkg,fkj,ekj logical lgamma integer*4 n,iel,kl,iep,ie,jl,lg,lk,na,nb,kg,jk dx=0.0d0 dy=0.0d0 dz=0.0d0 rx=0.0d0 ry=0.0d0 rz=0.0d0 qx=0.0d0 qy=0.0d0 qz=0.0d0 c permanent el moment (), ground state: dg1=u0(1)+u0d(1) dg2=u0(2)+u0d(2) dg3=u0(3)+u0d(3) c permanent el moment (), state k: dk1=ut(1,ie) dk2=ut(2,ie) dk3=ut(3,ie) kg=6*(iep-1) ekg=eau(iep) ekj=eau(iep)-eau(ie) if(lgamma)then fkg=ekg/(ekg**2+gammaau**2) fkj=ekj/(ekj**2+gammaau**2) else fkg=1.0d0/ekg fkj=1.0d0/ekj endif kg=6*(iep-1) c x/(Ek-Eg)/Nel: rrx=(dlg(5+kg)*dg3-dlg(6+kg)*dg2)*fkg rry=(dlg(6+kg)*dg1-dlg(4+kg)*dg3)*fkg rrz=(dlg(4+kg)*dg2-dlg(5+kg)*dg1)*fkg jk=(ie-1+n*(iep-1))*6 c x/(EK-EJ): qqx=(dd(5+jk)*dk3-dd(6+jk)*dk2)*fkj qqy=(dd(6+jk)*dk1-dd(4+jk)*dk3)*fkj qqz=(dd(4+jk)*dk2-dd(5+jk)*dk1)*fkj hx=0.0d0 hy=0.0d0 hz=0.0d0 fne=0.5d0/dble(na+nb) do 308 iel=1,n elg=eau(iel) ejg=eau(ie ) ekl=ekg-elg ejl=ejg-elg if(lgamma)then flg=elg/(elg**2+gammaau**2) fkl=ekl/(ekl**2+gammaau**2) else if(elg.eq.0.0d0)then call report('elg division by zero') write(6,*)iel,iep,elg,ekg,ekl stop endif fkl=1.0d0/ekl flg=1.0d0/elg endif kl=(iep-1+n*(iel-1))*6 lk=(iel-1+n*(iep-1))*6 lg=6*(iel-1) c sum(L<>n)x/(El-Eg)/Nel: rx=rx+(dd(2+kl)*dlg(6+lg)-dd(3+kl)*dlg(5+lg))*flg ry=ry+(dd(3+kl)*dlg(4+lg)-dd(1+kl)*dlg(6+lg))*flg rz=rz+(dd(1+kl)*dlg(5+lg)-dd(2+kl)*dlg(4+lg))*flg c x/(Ek-Eg)/Nel+ c sum(L<>k,L<>g)x/(Ek-EL)/Nel: if(iel.ne.iep)then rrx=rrx+(dd(5+kl)*dlg(3+lg)-dd(6+kl)*dlg(2+lg))*fkl rry=rry+(dd(6+kl)*dlg(1+lg)-dd(4+kl)*dlg(3+lg))*fkl rrz=rrz+(dd(4+kl)*dlg(2+lg)-dd(5+kl)*dlg(1+lg))*fkl endif c sum(L)x): dx=dx+dd(2+kl)*dlg(3+lg)-dd(3+kl)*dlg(2+lg) dy=dy+dd(3+kl)*dlg(1+lg)-dd(1+kl)*dlg(3+lg) dz=dz+dd(1+kl)*dlg(2+lg)-dd(2+kl)*dlg(1+lg) jl=(ie-1+n*(iel-1))*6 if(iel.ne.iep)then if(lgamma)then fkl=ekl/(ekl**2+gammaau**2) fjl=ejl/(ejl**2+gammaau**2) else if(ekl.eq.0.0d0)then write(6,*)iel,iep,elg,ekg,ekl call report('ekl division by zero') endif fkl=1.0d0/ekl fjl=1.0d0/ejl endif c sum(L)x/(EK-EL): qx=qx+(dd(2+jl)*dd(6+lk)-dd(3+jl)*dd(5+lk))*fkl qy=qy+(dd(3+jl)*dd(4+lk)-dd(1+jl)*dd(6+lk))*fkl qz=qz+(dd(1+jl)*dd(5+lk)-dd(2+jl)*dd(4+lk))*fkl c x/(EK-EJ)+ c sum(L<>j,L<>k)x/(EL-EJ): if(iel.ne.ie.and.iel.ne.iep)then qqx=qqx+(dd(5+jl)*dd(3+lk)-dd(6+jl)*dd(2+lk))*fjl qqy=qqy+(dd(6+jl)*dd(1+lk)-dd(4+jl)*dd(3+lk))*fjl qqz=qqz+(dd(4+jl)*dd(2+lk)-dd(5+jl)*dd(1+lk))*fjl endif c sum(L)x): hx=hx+dd(2+jl)*dd(3+lk)-dd(3+jl)*dd(2+lk) hy=hy+dd(3+jl)*dd(1+lk)-dd(1+jl)*dd(3+lk) hz=hz+dd(1+jl)*dd(2+lk)-dd(2+jl)*dd(1+lk) endif 308 continue dx=dx*fne dy=dy*fne dz=dz*fne rx=rx*fne ry=ry*fne rz=rz*fne qx=qx*fne qy=qy*fne qz=qz*fne rrx=rrx*fne rry=rry*fne rrz=rrz*fne qqx=qqx*fne qqy=qqy*fne qqz=qqz*fne hx=hx*fne hy=hy*fne hz=hz*fne return end c ============================================================ subroutine inispec(i,f,s) implicit none integer*4 i character*(*)f,s open(i,file=f) write(i,*)' '//s write(i,3000) 3000 format(' n wavelength (nm) dipole strength (D^2)', 1' rotatory strength (cgs/10**-36)',/,80(1h-)) return end c ============================================================ subroutine vz(v,n) implicit none integer*4 n,i real*8 v(*) do 1 i=1,n 1 v(i)=0.0d0 return end c ============================================================ subroutine cpolar(st,sav,sal,sgv,sgl,wmax,wmin,np,debye,pi,eau, 1gammaau,gammanm,lgnm,lgau,dl,dv,r,qv,ntr,n) implicit none integer*4 np,i,j,k,l,ie,ntr,n real*8 wmax,wmin,dw,cab,debye,pi,ccd,w,wau,al(3,3),ali(3,3), 1alit(3,3),allv(3,3),alilv(3,3),eau(*),gammaau,gammanm,t, 1dl(ntr,3),dv(ntr,3),r(ntr,3),qua,qv(ntr,3,3),fr,fi character*(*) st,sav,sal,sgv,sgl logical lgnm,lgau dw=(wmax-wmin)/dble(np-1) c constants making tensor traces epsilon c (factor of two missing for cab - do not know why): cab=debye**2 *108.7d0/pi c ccd=2.0d0*2.541765d0*9.2740154d-21*1.0d-18/1.0d-36/pi ccd=debye**2*4.0d0*108.7d0/pi/137.5d0 open(44,file=st) c velocity-based trace: open(9,file=sav) c dipole-based trace: open(91,file=sal) write(44,446)np,wmin,wmax 446 format(i6,2f12.2) write(44,443) 443 format(' Dynamic polarizabilities') c alpha start ----------------------------------------------------- write(44,444) 444 format(10x,' Alpha ( f ) Alpha ( g) length-length ', 1'and length-velocity : ', 1/,'l(nm) Re ullx Re ully Re ullz Im ullx Im ully', 1' Im ullz Re ulvx Re ulvy Re ulvz Im ulvx Im ulvy', 1' Im ulvz') w=wmin-dw do 102 k=1,np w=w+dw wau=(1.0d7/w)/219474.63d0 do 101 i=1,3 do 101 j=1,3 al(i,j)=0.0d0 ali(i,j)=0.0d0 alit(i,j)=0.0d0 allv(i,j)=0.0d0 alilv(i,j)=0.0d0 do 101 ie=1,n call setv(eau(ie),w,fr,fi,gammaau,gammanm,lgnm,lgau) c a_ab = 2 sum wjn (1/(wjn**2-w**2)) Re ua ub c length-length: t=dl(ie,i)*dl(ie,j)*eau(ie) al(i,j) =al(i,j) +t*fr ali(i,j)=ali(i,j)+t*fi c length-velocity: t=dl(ie,i)*dv(ie,j)*eau(ie) allv(i,j) =allv(i,j) +t*fr alilv(i,j)=alilv(i,j)+t*fi c velocity-velocity (for spectra calculation only): 101 alit(i,j)=alit(i,j)+dv(ie,i)*dv(ie,j)*fi*eau(ie) write(9,69)w, 1wau**2*cab*(alit(1,1)+alit(2,2)+alit(3,3)), 1 al( 1,1)+al( 2,2)+al( 3,3) write(91,69)w, 1wau**2*cab*(ali(1,1)+ali(2,2)+ali(3,3)), 1 al(1,1)+al( 2,2)+al( 3,3) 69 format(f7.1,2f10.2) write(44,441)w,(al( 1,j),j=1,3),(ali( 1,j),j=1,3), 1 (allv(1,j),j=1,3),(alilv(1,j),j=1,3) write(44,442) (al( 2,j),j=1,3),(ali( 2,j),j=1,3), 1 (allv(2,j),j=1,3),(alilv(2,j),j=1,3) 102 write(44,442) (al( 3,j),j=1,3),(ali( 3,j),j=1,3), 1 (allv(3,j),j=1,3),(alilv(3,j),j=1,3) 441 format(f8.2,12f12.4) 442 format(8x ,12f12.4) close(9) write(6,*)sav//' (velocity) written' close(91) write(6,*)sal//' (length ) written' c alpha end ----------------------------------------------------- c Gp start ----------------------------------------------------- write(44,445) 445 format(10x,' Gp/w ( f) Gp/w ( g )', 1/,'l(nm) ux uy uz') open(9 ,file=sgv) open(91,file=sgl) w=wmin-dw do 103 k=1,np w=w+dw wau=(1.0d7/w)/219474.63d0 do 104 i=1,3 do 104 j=1,3 al( i,j)=0.0d0 ali( i,j)=0.0d0 do 104 ie=1,n call setv(eau(ie),w,fr,fi,gammaau,gammanm,lgnm,lgau) c length: t=dl(ie,i)*r(ie,j) al( i,j) =al(i,j) +t*fr ali(i,j)=ali(i,j) +t*fi c velocity: 104 alit(i,j)=alit(i,j) +dv(ie,i)*r(ie,j)*fi write(91,69)w,( -ali(1,1)-ali( 2,2)-ali( 3,3))*wau**3*ccd, 1 al(1,1)+al( 2,2)+al( 3,3) write(9 ,69)w,(-alit(1,1)-alit(2,2)-alit(3,3))*wau**3*ccd, 1 al( 1,1)+al( 2,2)+al( 3,3) write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 103 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(9) close(91) write(6,*)sgv//' written' write(6,*)sgl//' written' c Gp end ----------------------------------------------------- c A start ----------------------------------------------------- write(44,447) 447 format(10x,' A ( f) A ( g) ', 1/,'l(nm) ux uy uz') w=wmin-dw do 105 k=1,np w=w+dw do 105 l=1,3 do 106 i=1,3 do 106 j=1,3 al( i,j)=0.0d0 ali(i,j)=0.0d0 do 106 ie=1,n if(i.eq.j)then qua=(3.0d0*qv(ie,i,j)-qv(ie,1,1)-qv(ie,2,2)-qv(ie,3,3))/2.0d0 else qua=1.5d0*qv(ie,i,j) endif call setv(eau(ie),w,fr,fi,gammaau,gammanm,lgnm,lgau) al(i,j) =al( i,j) +dl(ie,l)*qua*fr*eau(ie) 106 ali(i,j)=ali(i,j) +dl(ie,l)*qua*fi*eau(ie) if(l.eq.1)then write(44,441)w,(al(1,j),j=1,3),(ali(1,j),j=1,3) else write(44,442) (al(1,j),j=1,3),(ali(1,j),j=1,3) endif write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 105 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) close(44) c A end ----------------------------------------------------- write(6,*)st//' written' return end c ============================================================ subroutine addro(e,fi,fr,ucx,ucy,ucz,ng,jc, 1orb,oa,ob,nmo,rrxw,rryw,rrzw,rixw,riyw,rizw) implicit none real*8 e,fi,fr,ucx,ucy,ucz,ucxi,ucyi,uczi,ucxr,ucyr,uczr, 1orb(*),fab,rrxw(*),rryw(*),rrzw(*),rixw(*),riyw(*),rizw(*) integer*4 iz,ng,jc,oa,ob,nmo,ii fi=fi*e fr=fr*e ucxi=ucx*fi ucyi=ucy*fi uczi=ucz*fi ucxr=ucx*fr ucyr=ucy*fr uczr=ucz*fr do 9071 iz=1,ng fab=orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) ii=iz+jc rrxw(ii)=rrxw(ii)+ucxr*fab rryw(ii)=rryw(ii)+ucyr*fab rrzw(ii)=rrzw(ii)+uczr*fab rixw(ii)=rixw(ii)+ucxi*fab riyw(ii)=riyw(ii)+ucyi*fab 9071 rizw(ii)=rizw(ii)+uczi*fab return end c ============================================================ subroutine readtm(qz,x,bohr,xau,na,nb,nmo2,n,nmo,nst, 1dl,dv,qv,r,r0v,r0l,e,ev,eau,ni,nid,nat,cij,lopen, 1aij,bij,ifix,nib,cijb,aijb,bijb,lnorm,lort,lwrt) implicit none character*80 s80 logical lopen,lnorm,lort,lwrt integer*4 I,l,qz(*),na,nb,nmo2,n,k,ie,ni(*),ic,kk, 1nid(*),nat,aij(*),bij(*),j,oa,ob,jj,nst,ii,icc,it, 1ifix,nib(*),aijb(*),bijb(*),amax,bmax,nmo,nort,imax, 1ichek(100),nex,nf0,i1,i2,i3 parameter (nf0=64) integer*4 ifrom(nf0),ito(nf0),nff,iac,ns,nmosy(nf0), 1noccsy(nf0),ofrom(nf0),oto(nf0),syt(nf0,nf0),svec(nf0) logical lff character*50 ffn(nf0) real*8 x(3,nat),xau(3,nat),dl(n,3),dv(n,3),qv(n,3,3), 1r(n,3),r0v(*),r0l(*),e(*),ev(*),eau(*),cij(*), 1cijb(*),cmax,skj,sum1,sum2,bohr real*8,allocatable::obaj(:),obak(:) real*8 ska do 3 i=1,100 3 ichek(i)=0 write(6,*)'re-reading escf.out' n=0 open(2,file='escf.out') 2001 read(2,2000,end=23,err=23)s80 2000 format(a80) if(s80(15:32).EQ.'atomic coordinates')then do 1001 l=1,nat read(2,2000,end=23,err=23)s80 read(s80(1:43),*)(xau(i,l),i=1,3) x(1,l)=xau(1,l)*bohr x(2,l)=xau(2,l)*bohr x(3,l)=xau(3,l)*bohr 1001 read(s80(58:60),*)qz(l) ichek(1)=1 endif if(s80(2:30).EQ.'number of occupied orbitals :')then read(s80(31:len(s80)),*)na nb=na ichek(2)=1 endif if(s80(2:24).EQ.'Excitation energy: ')then n=n+1 read(s80(25:len(s80)),*)eau(n) ichek(3)=1 endif if(s80(2:24).EQ.'Excitation energy / eV:')then read(s80(25:len(s80)),*)ev(n) ichek(4)=1 endif if(s80(2:24).EQ.'Excitation energy / nm:')then read(s80(25:len(s80)),*)e(n) ichek(5)=1 endif if(s80(2:39).eq.'Electric transition dipole moment (vel')then read(2,*) read(2,2009)dv(n,1) read(2,2009)dv(n,2) read(2,2009)dv(n,3) 2009 format(5x,f16.6) ichek(6)=1 endif if(s80(2:39).eq.'Electric transition dipole moment (len')then read(2,*) read(2,2009)dl(n,1) read(2,2009)dl(n,2) read(2,2009)dl(n,3) ichek(7)=1 endif if(s80(2:34).eq.'Magnetic transition dipole moment')then read(2,*) read(2,2009)r(n,1) read(2,2009)r(n,2) read(2,2009)r(n,3) ichek(8)=1 endif if(s80(2:31).eq.'Electric quadrupole transition')then read(2,*) read(2,2009)qv(n,1,1) read(2,2009)qv(n,2,2) read(2,2009)qv(n,3,3) read(2,2009)qv(n,1,2) read(2,2009)qv(n,1,3) read(2,2009)qv(n,2,3) qv(n,3,2)=qv(n,2,3) qv(n,3,1)=qv(n,1,3) qv(n,2,1)=qv(n,1,2) ichek(9)=1 endif if(s80(2:19).eq.'Rotatory strength:')then read(2,2000,end=23,err=23)s80 read(2,2000,end=23,err=23)s80 read(s80(39:len(s80)),*)r0v(n) read(2,2000,end=23,err=23)s80 read(2,2000,end=23,err=23)s80 read(2,2000,end=23,err=23)s80 read(2,2000,end=23,err=23)s80 read(s80(39:len(s80)),*)r0l(n) ichek(10)=1 endif if(s80(11:40).eq.'IRREP tensor space dimension')then nff=0 iac=0 read(2,*) 944 read(2,2000)s80 if(s80(11:13).ne.' ')then nff=nff+1 if(nff.gt.nf0)call report('nff.gt.nf0') read(s80(20:30),*)ii ifrom(nff)=iac+1 ito(nff)=iac+ii write(6,*)'symmetry',nff,' from state ', 1 ifrom(nff),' to ',ito(nff) iac=iac+ii goto 944 endif ichek(11)=1 endif if(s80(4:13).eq.'irrep mo')then nff=0 iac=0 945 read(2,2000)s80 if(s80(3:8).ne.' ')then nff=nff+1 if(nff.gt.nf0)call report('nff.gt.nf0') read(s80(8:len(s80)),*)ii,noccsy(nff) nmosy(nff)=ii ofrom(nff)=iac+1 oto(nff)=iac+ii write(6,*)'symmetry',nff,' from MO ',iac+1,' to ',oto(nff), 1 ' occ:',noccsy(nff) iac=iac+ii goto 945 endif ichek(12)=1 endif goto 2001 23 close(2) write(6,*)n,' excitations from escf.out' if(ichek( 1).eq.0)call report(' Atomic coordinates missing') if(ichek( 2).eq.0)call report(' Nocc missing') if(ichek( 3).eq.0)call report(' E missing') if(ichek( 4).eq.0)call report(' E-eV missing') if(ichek( 5).eq.0)call report(' E-nm missing') if(ichek( 6).eq.0)call report(' Dipole-v missing') if(ichek( 7).eq.0)call report(' Dipole-l missing') if(ichek( 8).eq.0)call report(' Magnet missing') if(ichek( 9).eq.0)call report(' Quadrupole missing') if(ichek(10).eq.0)call report(' Rotatory strength missing') if(ichek(11).eq.0)call report(' Number of symmetries missing') if(ichek(12).eq.0)call report(' Orbital symmetries missing') if(nff.ne.1)then open(2,file='sym',status='old') read(2,*)ie if(ie.ne.nff)call report('ie<>nff') do 556 i=1,nff 556 read(2,*)(syt(i,j),j=1,nff) write(6,*)'sym with symmetry table found' else syt(1,1)=1 endif ie=0 lopen=.false. inquire(file='F.LST',exist=lff) if(lff)then open(8,file='F.LST') do 800 i=1,nff 800 read(8,8000)ffn(i) 8000 format(a50) else if(nff.ne.1)call report('F.LST missing') ffn(1)='sing_a' ifrom(1)=1 ito(1)=nmo endif it=0 do 555 ic=1,nff write(6,*)'reading '//ffn(ic) open(2,file=ffn(ic),status='old') 1 read(2,2000,end=2,err=2)s80 if(s80(2:27).eq.'current subspace dimension') 1read(s80(28:36),*),ns if(s80(13:24).eq.'eigenvalue =')then write(6,2000)s80 read(s80(1:9),*)ie it=it+1 write(6,*)' total eigenvalue ',it ii=0 c occupied symmetries do 557 i1=1,nff c virtual symmetries do 557 i2=1,nff c does is give symmetry ic?: do 558 i3=1,nff 558 svec(i3)=syt(i1,i3)*syt(i2,i3) if(ie.eq.1)write(6,*)i1,i2,(svec(i3),i3=1,nff) icc=0 do 559 i3=1,nff do 5591 jj=1,nff 5591 if(syt(i3,jj).ne.svec(jj))goto 559 icc=i3 559 continue if(icc.eq.0)call report('symmetry calc. error') if(ie.eq.1)write(6,*)i1,'x',i2,'corresponds to ',icc,' want',ic if(icc.eq.ic)then do 300 i=ofrom(i1),ofrom(i1)+noccsy(i1)-1 do 300 j=ofrom(i2)+noccsy(i2),ofrom(i2)+nmosy(i2)-1 ii=ii+1 aij( nmo2*(it-1)+ii)=i aijb(nmo2*(it-1)+ii)=i bijb(nmo2*(it-1)+ii)=j if(j.gt.nmo.or.i.gt.nmo)then write(6,*)i1,i2,i,j,ofrom(i2),noccsy(i2),nmosy(i2) call report('overflow2') endif 300 bij( nmo2*(it-1)+ii)=j endif 557 continue if(ie.eq.1)write(6,*)ii,' states determined from orbitals' ni(it)=ns if(ii.ne.ns)then write(6,*)ic,ie,it,ni(ie),ns,noccsy(ic),nmosy(ic) call report('incosistent numbering') endif nib(it)=ni(it) nid(it)=0 cmax=0.0d0 amax=0 bmax=0 imax=0 sum1=0.0d0 sum2=0.0d0 read(2,900)(cij(nmo2*(it-1)+ii),ii=1,ni(it)), 1 (cijb(nmo2*(it-1)+ii),ii=1,nib(it)) 900 format(4d20.14) nex=0 do 301 ii=1,ni(it) call setabmax(aij(nmo2*(it-1)+ii),bij(nmo2*(it-1)+ii), 1 cij(nmo2*(it-1)+ii),amax,bmax,cmax,ii,imax) sum1=sum1+cij(nmo2*(it-1)+ii)**2 sum2=sum2+cijb(nmo2*(it-1)+ii)**2 c experimental fix for the deexitations, closed shell only: if(ifix.eq.1) 1 cij(nmo2*(it-1)+ii)=cij(nmo2*(it-1)+ii)+cijb(nmo2*(it-1)+ii) c experimental fix,just add to excitations with switched orbitals: if(ifix.eq.3)then nex=nex+1 cij(nmo2*(it-1)+ni(it)+nex)=cijb(nmo2*(it-1)+ii) aij(nmo2*(it-1)+ni(it)+nex)=bijb(nmo2*(it-1)+ii) bij(nmo2*(it-1)+ni(it)+nex)=aijb(nmo2*(it-1)+ii) nib(it)=0 endif 301 continue ni(it)=ni(it)+nex if(lwrt)then write(6,*)ni(it),' alpha coefficitnts' if(ifix.eq.3)write(6,*)'(with backwards)' write(6,*)nid(it),' beta coefficitnts' write(6,*)nib(it),' backward coefficitnts' write(6,6009)amax,bmax,cij(nmo2*(it-1)+imax),cmax**2*100.0d0 6009 format(' Dominant transition from',i4,' to ',i4,f8.3, 1 ' (',f6.2,' %)') if(.not.lopen)write(6,1234)sum1,sum2 1234 format('fnorm, bnorm:',2f10.4) endif if(ifix.eq.4.and.lwrt)then nib(it)=0 write(6,*)'backward ignored' endif c seems that fnorm-bnorm=1 endif goto 1 c 2 close(2) 555 continue if(lnorm)then c renormalize: call normcij(it,cij,ni,nmo2) write(6,*)' cij were normalized' endif if(lort)then if(nst.gt.0.and.nst.lt.n)then nort=nst else nort=n endif if(lwrt)then do i=1,n write(6,6000)i,e(i) 6000 format(' state ',i4,f12.5) do ii=1,ni(i) j=nmo2*(i-1) if(dabs(cij(j+ii)).gt.0.1d0) 1 write(6,6001)aij(j+ii),bij(j+ii),cij(j+ii) 6001 format(2i4,f12.6) enddo enddo endif c reorthonormalize: write(6,*)'orthogonalizing coefficients for ',nort,' states' if(lopen)call report('lort for open shell not implemented') allocate(obaj(nmo2),obak(nmo2)) c run twice for better accuracy: do 318 ic=1,2 ska=0.0d0 do 3181 j=1,nort c transcript original cij to obaj: do 320 jj=1,nmo2 320 obaj(jj)=0.0d0 do 324 jj=1,ni(j) oa=aij(nmo2*(j-1)+jj) ob=bij(nmo2*(j-1)+jj) if(ob+nmo*(oa-1).gt.nmo2)then write(6,*)oa,ob,nmo2 call report('overflow') endif 324 obaj(ob+nmo*(oa-1))=cij(nmo2*(j-1)+jj) do 319 k=1,j-1 c transcript cij to obak: do 321 kk=1,nmo2 321 obak(kk)=0.0d0 kk=nmo2*(k-1) do 3241 jj=1,ni(k) oa=aij(kk+jj) ob=bij(kk+jj) 3241 obak(ob+nmo*(oa-1))=cij(kk+jj) c : skj=0.0d0 do 322 oa=1,na do 322 ob=na+1,nmo c forward excitations: skj=skj+obaj(ob+nmo*(oa-1))*obak(ob+nmo*(oa-1)) c backward excitations: 322 skj=skj+obaj(oa+nmo*(ob-1))*obak(oa+nmo*(ob-1)) c |j'>=|j>-|k>: do 319 oa=1,na do 319 ob=na+1,nmo obaj(ob+nmo*(oa-1))=obaj(ob+nmo*(oa-1))-skj*obak(ob+nmo*(oa-1)) 319 obaj(oa+nmo*(ob-1))=obaj(oa+nmo*(ob-1))-skj*obak(oa+nmo*(ob-1)) ska=ska+dabs(skj) c c transcript new obaj to new cij: ni(j)=0 do 325 oa=1,na do 325 ob=na+1,nmo jj=ob+nmo*(oa-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(nmo2*(j-1)+ni(j))=obaj(jj) aij(nmo2*(j-1)+ni(j))=oa bij(nmo2*(j-1)+ni(j))=ob endif jj=oa+nmo*(ob-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(nmo2*(j-1)+ni(j))=obaj(jj) aij(nmo2*(j-1)+ni(j))=ob bij(nmo2*(j-1)+ni(j))=oa endif 325 continue 3181 call normcij(j,cij,ni,nmo2) write(6,*)ic,ska/dble(nort**2)/2.0d0 318 continue write(6,*)'double orthogonalization done' deallocate(obaj,obak) endif return end subroutine uun(un,nat,xau,qq) implicit none real*8 un(*),xau(3,nat) integer*4 nat,j,l,qq(*) do 311 j=1,3 un(j)=0.0d0 do 311 l=1,nat 311 un(j)=un(j)+xau(j,l)*qq(l) return end subroutine rdps(nat,ipse) c read pseudopotential parameters implicit none character*29 s29 integer*4 i,nat,ipse(*) read(2,*) read(2,*) read(2,*) read(2,*) i=0 1 read(2,200)s29 200 format(a29) if(s29(15:16).eq.' ')goto 1 if(s29(15:16).eq.'==')return i=i+1 if(i.gt.nat)call report('pseudopotential reading error') if(s29(28:29).eq.' ')then ipse(i)=0 else read(s29(20:29),*)ipse(i) endif goto 1 end subroutine dtm(e,a,i,c,m) real*8 a,c(*) integer*4 e,i(*),m do 1 j=1,m if(dabs(a).gt.dabs(c(j)))then do 2 k=m,j+1,-1 i(k)=i(k-1) 2 c(k)=c(k-1) i(j)=e c(j)=a return endif 1 continue return end subroutine readtda(na,nb,nmo2,n,nmo,nst,e,ev,eau,ni,lden,nid,cij, 1lopen,aij,bij,dij,daij,dbij,ifix,nib,lnorm,lort,lwrt,nibd,ncmm) implicit none character*80 s80 logical lden,lopen,lnorm,lort,lwrt integer*4 na,nb,nmo2,n,nd,k,ie,ni(*),ic,kk,nid(*),aij(*),bij(*), 1dbij(*),daij(*),j,oa,ob,jj,nst,ifix,nib(*),amax,bmax,nmo,nort, 1nibd(*),ii,n22,ncmm,imax real*8 e(*),ev(*),eau(*),cij(*),dij(*),cmax,skj,sum1,sum2,ecm,ska real*8,allocatable::obaj(:),obak(:),obajd(:),obakd(:) if(ncmm.ne.0)then n22=ncmm else n22=nmo2 endif ie=0 lopen=.false. open(2,file='ciss_a.co') 1 read(2,2000,end=2,err=2)s80 write(6,2000)s80 2000 format(a80) c if(s80(1:1).eq.'e'.or.s80(1:2).eq.' ')then ie=ie+1 read(s80(3:len(s80)),*)nd,eau(nd) ev(nd)=eau(nd)*27.211384205943d0 ecm=ev(nd)*8065.54476345045d0 e(nd)=1.0d7/ecm if(ie.ne.nd)call report(' Inconsistent reading') cmax=0.0d0 amax=0 bmax=0 imax=0 if(lden)then ni(ie)=0 nib(ie)=0 nibd(ie)=0 nid(ie)=0 sum1=0.0d0 sum2=0.0d0 111 read(2,2000,end=299,err=299)s80 if(s80(1:1).ne.'e'.and.s80(1:2).ne.' ')then if(ni(ie).gt.n22)call report('too many coefficients') c closed shell: ni(ie)=ni(ie)+1 ii=n22*(ie-1)+ni(ie) read(s80,*)cij(ii),aij(ii),bij(ii) call setabmax(aij(ii),bij(ii),cij(ii),amax,bmax,cmax, 1 ii,imax) sum1=sum1+cij(ii)**2 goto 111 else backspace 2 endif 299 continue if(lwrt)then write(6,*)ni(ie),' alpha coefficients' if(ifix.eq.3)write(6,*)'(with backwards)' write(6,*)nid(ie),' beta coefficients' write(6,*)nib(ie),' backward coefficients' write(6,6009)amax,bmax,cij(imax),cmax**2*100.0d0 6009 format(' Dominant transition from',i4,' to ',i4,2f8.3) write(6,1234)sum1,sum2,sum1-sum2 1234 format('fnorm, bnorm,diff:',3f10.4) endif if(lnorm)then c renormalize: call normcijo(ie,cij,dij,ni,nid,n22,lopen) write(6,*)' cij were normalized' endif endif if(nd.eq.n) goto 2 endif goto 1 c 2 close(2) if(lort)then if(nmo2.ne.n22)call report('cannot re-orthogonalize' 1 //'within limited coefficient space') if(nst.gt.0.and.nst.lt.n)then nort=nst else nort=n endif c reorthonormalize: write(6,*)'orthogonalizing coefficients for ',nort,' states' allocate(obaj(nmo2),obak(nmo2)) if(lopen)allocate(obajd(nmo2),obakd(nmo2)) write(6,*)'temporary variables allocated' c run twice for better accuracy: do 318 ic=1,2 ska=0.0d0 do 3181 j=1,nort c transcript original cij to obaj: call vz(obaj,nmo2) do 324 jj=1,ni(j) oa=aij(n22*(j-1)+jj) ob=bij(n22*(j-1)+jj) 324 obaj(ob+nmo*(oa-1))=cij(n22*(j-1)+jj) if(lopen)then call vz(obajd,nmo2) do 3241 jj=1,nid(j) oa=daij(n22*(j-1)+jj) ob=dbij(n22*(j-1)+jj) 3241 obajd(ob+nmo*(oa-1))=dij(n22*(j-1)+jj) endif do 334 k=1,j-1 c transcript cij to obak: call vz(obak,nmo2) kk=n22*(k-1) do 3242 jj=1,ni(k) oa=aij(kk+jj) ob=bij(kk+jj) 3242 obak(ob+nmo*(oa-1))=cij(kk+jj) if(lopen)then call vz(obakd,nmo2) kk=n22*(k-1) do 3243 jj=1,nid(k) oa=daij(kk+jj) ob=dbij(kk+jj) 3243 obakd(ob+nmo*(oa-1))=dij(kk+jj) endif c : skj=0.0d0 do 323 oa=1,na do 323 ob=na+1,nmo c forward excitations: skj=skj+obaj(ob+nmo*(oa-1))*obak(ob+nmo*(oa-1)) c backward excitations: 323 skj=skj+obaj(oa+nmo*(ob-1))*obak(oa+nmo*(ob-1)) if(lopen)then do 331 oa=1,nb do 331 ob=nb+1,nmo skj=skj+obajd(ob+nmo*(oa-1))*obakd(ob+nmo*(oa-1)) 331 skj=skj+obajd(oa+nmo*(ob-1))*obakd(oa+nmo*(ob-1)) endif c |j'>=|j>-|k>: do 319 oa=1,na do 319 ob=na+1,nmo if(oa+nmo*(ob-1).gt.nmo2)write(6,*)'A',oa,ob,nmo2 if(ob+nmo*(oa-1).gt.nmo2)write(6,*)'B',oa,ob,nmo2 obaj(ob+nmo*(oa-1))=obaj(ob+nmo*(oa-1))-skj*obak(ob+nmo*(oa-1)) 319 obaj(oa+nmo*(ob-1))=obaj(oa+nmo*(ob-1))-skj*obak(oa+nmo*(ob-1)) if(lopen)then do 332 oa=1,nb do 332 ob=nb+1,nmo obajd(ob+nmo*(oa-1))= 1 obajd(ob+nmo*(oa-1))-skj*obakd(ob+nmo*(oa-1)) 332 obajd(oa+nmo*(ob-1))= 1 obajd(oa+nmo*(ob-1))-skj*obakd(oa+nmo*(ob-1)) endif 334 ska=ska+dabs(skj) c c transcript new obaj to new cij: if(lopen)then nid(j)=0 do 333 oa=1,nb do 333 ob=nb+1,nmo jj=ob+nmo*(oa-1) if(obajd(jj).ne.0.0d0)then nid(j)=nid(j)+1 dij( n22*(j-1)+nid(j))=obajd(jj) daij(n22*(j-1)+nid(j))=oa dbij(n22*(j-1)+nid(j))=ob endif jj=oa+nmo*(ob-1) if(obajd(jj).ne.0.0d0)then nid(j)=nid(j)+1 dij( n22*(j-1)+nid(j))=obajd(jj) daij(n22*(j-1)+nid(j))=ob dbij(n22*(j-1)+nid(j))=oa endif 333 continue endif ni(j)=0 do 325 oa=1,na do 325 ob=na+1,nmo jj=ob+nmo*(oa-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(n22*(j-1)+ni(j))=obaj(jj) aij(n22*(j-1)+ni(j))=oa bij(n22*(j-1)+ni(j))=ob endif jj=oa+nmo*(ob-1) if(obaj(jj).ne.0.0d0)then ni(j)=ni(j)+1 cij(n22*(j-1)+ni(j))=obaj(jj) aij(n22*(j-1)+ni(j))=ob bij(n22*(j-1)+ni(j))=oa endif 325 continue 3181 call normcijo(j,cij,dij,ni,nid,n22,lopen) write(6,*)ic,ska/dble(nort**2)/2.0d0 318 continue write(6,*)'double orthogonalization done' endif return end c ============================================================ subroutine imcdv(fn,iwr,LDD,nroot,nat,C,D,G,GP,u0,ddi,N,fac,e00, 1gnji,lgnj,LDE,g0,NQ,m0,aai) implicit none character*(*) fn character*80 s80 integer*4 iwr,nroot,ip(10),ix,iy,iz,ia,iw,N,nat,IERR,NQ logical LDD,lgnj,LDE real*8 C(N,N),D(*),u0(*),e00,G(N,N),GP(N,N),ddi(*),fac, 1sp1,sp2,sp,proc,gnji(*),g0(*),m0(*),aai(*) real*16 GETDET real*8,allocatable::J(:,:),JT(:,:),K(:),DT(:), 1A(:,:),B(:),E(:,:),F(:,:),FI(:,:),W(:),WP(:),dd(:),T(:,:),TP(:,:), 1TU(:),TV(:),EE(:,:),gnjx(:),CT(:,:),GT(:,:),GPT(:,:),aa(:) NQ=N if(LDD)then allocate(dd(9*nat),aa(9*nat)) if(LDE)then call rde(dd,nat,u0,aa,m0) else call rdd(fn,dd,nat,nroot,u0,iwr,m0,aa) endif call trafo(dd,ddi,nat,NQ) call trafo(aa,aai,nat,NQ) deallocate(dd,aa) call wrda(ddi,aai,NQ,nroot) endif if(lgnj)then allocate(gnjx(27*nat)) open(88,file='GNJD.TEN') do 1 ia=1,nat do 1 iz=1,3 read(88,*) do 1 ix=1,3 1 read(88,*)(gnjx(ix+3*(iy-1)+9*(iz-1)+27*(ia-1)),iy=1,3) c ix-electric,iy-magnetic,iz-atomic read(88,*) do 5 ix=1,3 5 read(88,*)(g0(ix+3*(iy-1)),iy=1,3) close(88) call trafog(gnjx,gnji,nat,NQ) deallocate(gnjx) call wrgn(gnji,NQ,nroot) endif allocate(JT(NQ,NQ),K(NQ),A(NQ,NQ),B(NQ),CT(NQ,NQ),J(NQ,NQ), 1FI(NQ,NQ),F(NQ,NQ),EE(NQ,2*NQ),E(NQ,NQ),GT(NQ,NQ),GPT(NQ,NQ), 1T(NQ,NQ),TP(NQ,NQ),TU(NQ),TV(NQ),DT(NQ)) call zm(A,NQ) call zv(B,NQ) call zm(CT,NQ) call zm(C,N) call zv(D,N) call zv(DT,NQ) call zm(E,NQ) call zm(J,NQ) call zm(JT,NQ) call zv(K,NQ) allocate(W(N),WP(N)) call zv(W,N) call zv(WP,N) e00=0.0d0 do 3 ix=1,10 3 ip(ix)=0 iw=0 open(9,file=fn,status='old') 2 read(9,900,end=88,err=88)s80 900 format(a80) c explore alternatives of the fucking Gaussian output: if(s80(4:20).eq.'Duschinsky Matrix')call rdm(9,NQ,JT,ip(1)) if(s80(2:18).eq.'Duschinsky matrix')then read(9,*) read(9,*) read(9,*) call rdm(9,NQ,JT,ip(1)) endif if(s80(2:24).eq.'Final Duschinsky matrix')then read(9,*) read(9,*) read(9,*) call rdm(9,NQ,JT,ip(1)) endif if(s80(4:15).eq.'Shift Vector')call rdv(9,NQ,K,ip(2)) if(s80(2:13).eq.'Shift Vector')call rdv(9,NQ,K,ip(2)) if(s80(4:11).eq.'A Matrix')call rdm(9,NQ,A,ip(3)) if(s80(2:9).eq.'A Matrix' )call rdm(9,NQ,A,ip(3)) if(s80(4:11).eq.'B Vector')call rdv(9,NQ,B,ip(4)) if(s80(2:9).eq.'B Vector' )call rdv(9,NQ,B,ip(4)) if(s80(4:11).eq.'C Matrix')call rdm(9,NQ,CT,ip(5)) if(s80(2:9).eq.'C Matrix' )call rdm(9,NQ,CT,ip(5)) if(s80(4:11).eq.'D Vector')call rdv(9,NQ,DT,ip(6)) if(s80(2:9).eq.'D Vector' )call rdv(9,NQ,DT,ip(6)) if(s80(4:11).eq.'E Matrix')call rdm(9,NQ,E,ip(7)) if(s80(2:9).eq.'E Matrix' )call rdm(9,NQ,E,ip(7)) if(s80(2:30).eq.'Energy of the 0-0 transition:') 1read(s80(31:42),*)e00 if(s80(2:22).eq.'Harmonic frequencies ')call rw(9,N,W,WP,iw) goto 2 88 close(9) 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.ne.2)call report('Frequencies not found') c reduce number of frequencies according to the c excited state vibrations: do 7 ix=NQ,N-1 do 7 iy=1,N-ix+NQ-1 WP(iy)=WP(iy+1) 7 W(iy)=W(iy+1) call mz(G,N) call mz(GP,N) call mz(GT,NQ) call mz(GPT,NQ) do 4 ix=1,NQ GT(ix,ix)=W(ix) GPT(ix,ix)=WP(ix) G(ix,ix)=W(ix) 4 GP(ix,ix)=WP(ix) do 6 ix=1,NQ D(ix)=DT(ix) do 6 iy=1,NQ 6 C(ix,iy)=CT(ix,iy) call mt(J,JT,NQ) call mm(T,GT,J,NQ) call mm(TP,JT,T,NQ) call ms(F,TP,GPT,NQ) call INV(NQ,F,FI,NQ,EE,IERR) if(IERR.ne.0)call report('Inversion error') call mm(T,GT,GPT,NQ) c fac=dsqrt(2.0d0**NQ)*dbleq(sqrt(sqrt(abs(GETDET(T,NQ)))) fac=dsqrt(2.0d0**NQ)*real(sqrt(sqrt(abs(GETDET(T,NQ)))) 1*sqrt(abs(GETDET(J,NQ)/GETDET(F,NQ))),8) call mv(TV,GT,K,NQ) sp1=sp(K,TV,NQ) call mv(TU,GT,K,NQ) call mv(TV,JT,TU,NQ) call mv(TU,FI,TV,NQ) call mv(TV,J,TU,NQ) call mv(TU,GT,TV,NQ) sp2=sp(K,TU,NQ) fac=fac*exp(-0.50d0*(sp1-sp2)) proc=fac*fac*100.0d0 write(6,300)fac,proc 300 format('<0|0*> = ',f20.13,' (',f6.2,'%)') return end subroutine mcdi(e00,fac,np0,iwr,N,C,D,GP,LEXCL,u0,ddi,Gnj, 1LQ1,LDD,THRC,HW,lgnj,gnji,g0,NQ,m0,aai,rgnj) c e00 ... energy of the 0-0 transition c fac = <0|0*> c dimension of the expansion space c iwr .. extra writing for debug c Gnj ... "MCD tensor" for this transition c gnji ... MCD tensor derivatives implicit none integer*4 nt,N,LEXCL,ix,kk,np0,NQ1,NExc,iex,I,icen,II,ic, 1ik,ip,iwr,nuk,NQ real*8 pp,FC,e00,fac,ut(3),u0(*),ddi(*),Gnj(*),C(N,N),D(N), 1GP(N,N),proc,ds,enm,EJ,CM,THRC,pm,qk,rs,gnji(*),gt(9),g0(*), 1m0(*),aai(*),mt(3),mcd,vt(3) logical LQ1,LDD,HW,lgnj,rgnj integer,allocatable::si(:),sit(:) CM=219470.0d0 proc=0.0d0 nt=0 c <0|0*>: nt=nt+1 proc=proc+100.0d0*fac**2 open(70,file='MCDI.MOM.TXT') allocate(si(LEXCL+1),sit(LEXCL+1)) do 11 kk=1,LEXCL+1 11 si(kk)=0 if(rgnj)call rdgnj(Gnj,u0,m0) c transition el and magn dipole <0|u|*>=u0<0|*>+du/dQk <0|Qk|*>,etc do 12 ix=1,3 mt(ix)=m0(ix)*fac 12 ut(ix)=u0(ix)*fac if(lgnj)then c transition tensor <0|G|*>=G0<0|*>+dG/dQk <0|Qk|*> do 1 ix=1,9 1 Gt(ix)=g0(ix)*fac else do 101 ix=1,9 101 Gt(ix)=Gnj(ix)*fac endif if(LQ1)then do 10 kk=1,NQ si(1)=kk c <0|Q*k|0*>=sqrt(h/(2*wk)) <0|1*k> pp=FC(fac,si,1,np0,LEXCL,iwr,C,D,N)/dsqrt(2.0d0*GP(kk,kk)) if(LDD)then if(lgnj)then do 2 ix=1,9 2 Gt(ix)=Gt(ix)+gnji(ix+9*(kk-1))*pp endif do 13 ix=1,3 mt(ix)=mt(ix)+ aai(ix+3*(kk-1))*pp 13 ut(ix)=ut(ix)+ ddi(ix+3*(kk-1))*pp endif 10 write(6,1600)kk,pp 1600 format('<0|Q_',i3,'|0*> = ',g13.4) endif vt(1)=Gt(2+3*(3-1))-Gt(3+3*(2-1)) vt(2)=Gt(3+3*(1-1))-Gt(1+3*(3-1)) vt(3)=Gt(1+3*(2-1))-Gt(2+3*(1-1)) ds =ut(1)*ut(1)+ut(2)*ut(2)+ut(3)*ut(3) rs =ut(1)*mt(1)+ut(2)*mt(2)+ut(3)*mt(3) mcd=ut(1)*vt(1)+ut(2)*vt(2)+ut(3)*vt(3) c cm-1 to nm: enm=1.0d7/e00 write(42,550)nt,enm,ds,mcd,proc write(43,550)nt,enm,ds,rs ,proc 550 format(i8,f10.2,2(' ',g13.4),' ',f10.2,'%',10i3) write(70,700)nt,enm,(ut(I),I=1,3),(mt(I),I=1,3),(vt(I),I=1,3) 700 format(i8,f10.2,9g13.4) c states Nexc excited NQ1=LEXCL do 30000 Nexc=1,LEXCL c distribute Nexc excitations upon centers 1..NQ: c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=1 50000 do 11000 I=1,Nexc c check that the state is allowed: icen=1 do 11000 II=I+1,Nexc if(si(II).eq.si(I))icen=icen+1 11000 if(icen.gt.NQ1)goto 12000 c EJ=0.0d0 DO 1002 II=1,Nexc 1002 EJ=EJ+GP(si(II),si(II)) if(iwr.gt.1)write(6,604)ii,(si(iex),iex=1,NExc) 604 format(' state ',i2,':',10i3) c <0|*> pp=FC(fac,si,Nexc,np0,LEXCL,iwr,C,D,N) proc=proc+pp*pp*100.0d0 if(HW)then write(6,599) 599 format('<0|',$) do 107 I=1,Nexc 107 write(6,601)si(I) 601 format(i3,$) write(6,598) 598 format('>: ',$) write(6,602)pp,proc,EJ*CM 602 format(g12.4,' (',f10.3,'%), E= ',F10.2) endif c <0|u|*>=u0<0|*>+du/dQk<0|Qk|*> c <0|Qk|*>=sqrt(h/(2wk))(sqrt(vk)<0|*-1k>+sqrt(vk+1)<0|*+1k>) do 14 ix=1,3 ut(ix)=u0(ix)*pp 14 mt(ix)=m0(ix)*pp if(lgnj)then do 3 ix=1,9 3 Gt(ix)=g0(ix)*pp endif if(LQ1)then if(LDD)then do 15 kk=1,NQ nuk=0 do 16 ic=1,Nexc 16 if(si(ic).eq.kk)nuk=nuk+1 pm=0.0d0 if(nuk.gt.0)then ip=0 ik=0 do 17 ic=1,Nexc if(si(ic).eq.kk.and.ik.eq.0)then ik=ik+1 else ip=ip+1 sit(ip)=si(ic) endif 17 continue pm=dsqrt(dble(nuk))*FC(fac,sit,Nexc-1,np0,LEXCL,iwr,C,D,N) endif do 18 ic=1,Nexc 18 sit(ic)=si(ic) sit(Nexc+1)=kk pp=dsqrt(dble(nuk+1))*FC(fac,sit,Nexc+1,np0,LEXCL,iwr,C,D,N) qk=(pm+pp)/dsqrt(2.0d0*GP(kk,kk)) if(lgnj)then do 4 ix=1,9 4 Gt(ix)=Gt(ix)+gnji(ix+9*(kk-1))*qk endif do 15 ix=1,3 ut(ix)=ut(ix)+ddi(ix+3*(kk-1))*qk 15 mt(ix)=mt(ix)+aai(ix+3*(kk-1))*qk endif endif vt(1)=Gt(2+3*(3-1))-Gt(3+3*(2-1)) vt(2)=Gt(3+3*(1-1))-Gt(1+3*(3-1)) vt(3)=Gt(1+3*(2-1))-Gt(2+3*(1-1)) ds =ut(1)*ut(1)+ut(2)*ut(2)+ut(3)*ut(3) rs =ut(1)*mt(1)+ut(2)*mt(2)+ut(3)*mt(3) mcd=ut(1)*vt(1)+ut(2)*vt(2)+ut(3)*vt(3) if(dabs(ds).gt.THRC)then nt=nt+1 enm=1.0d7/(e00+EJ*CM) write(42,550)nt,enm,ds,mcd,proc,(si(i),i=1,Nexc) write(43,550)nt,enm,ds,rs ,proc,(si(i),i=1,Nexc) write(70,700)nt,enm,(ut(I),I=1,3),(mt(I),I=1,3),(vt(I),I=1,3) endif c c find index to be changed 12000 do 80000 ic=Nexc,1,-1 80000 if(si(ic).lt.NQ)goto 90000 goto 30000 90000 do 10000 i=ic+1,Nexc 10000 si(i)=si(ic)+1 si(ic)=si(ic)+1 c goto 50000 30000 continue return end subroutine mz(A,N) implicit none integer*4 n,i,j real*8 A(N,N) do 3 i=1,N do 3 j=1,N 3 A(i,j)=0.0d0 return end SUBROUTINE INV(nr,a,ai,n,e,IERR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension ai(nr,nr),a(nr,nr),e(nr,2*nr) C TOL=1.0d-10 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 mm(A,B,C,N) c matrix multiplication A = B x C implicit none integer*4 N,i,j,k real*8 A(N,N),B(N,N),C(N,N) do 1 i=1,N do 1 j=1,N A(i,j)=0.0d0 do 1 k=1,N 1 A(i,j)=A(i,j)+B(i,k)*C(k,j) return end subroutine mv(A,B,C,N) c multiplication of matrix B by vector C, A = B . C implicit none integer*4 N,i,k real*8 A(*),B(N,N),C(*) do 1 i=1,N A(i)=0.0d0 do 1 k=1,N 1 A(i)=A(i)+B(i,k)*C(k) return end function sp(A,C,N) c scalar product of two vectors implicit none integer*4 N,k real*8 A(*),C(*),sp sp=0.0d0 do 1 k=1,N 1 sp=sp+A(k)*C(k) return end subroutine ms(A,B,C,N) c matrix sum A = B + C implicit none integer*4 N,i,j real*8 A(N,N),B(N,N),C(N,N) do 1 i=1,N do 1 j=1,N 1 A(i,j)=B(i,j)+C(i,j) return end subroutine mt(AT,A,N) c matrix transposition AT = A^T implicit none integer*4 N,i,j real*8 A(N,N),AT(N,N) do 1 i=1,N do 1 j=1,N 1 AT(i,j)=A(j,i) return end subroutine rw(o,N,a,b,i) implicit none integer*4 o,N,i,ns,ne,j real*8 a(*),b(*),cmtoau character*80 s cmtoau=219474.63d0 if(i.ge.2)return i=i+1 ns=1 1 read(o,9)s 9 format(a80) if(s(2:15).eq.'Frequencies --')then ne=min(N,ns+2) if(i.eq.1)read(s(16:len(s)),*)(a(j),j=ns,ne) if(i.eq.2)read(s(16:len(s)),*)(b(j),j=ns,ne) c convert to Hartree: do 2 j=ns,ne if(i.eq.1)a(j)=a(j)/cmtoau 2 if(i.eq.2)b(j)=b(j)/cmtoau ns=min(ns+3,N) if(ne.eq.N)return endif goto 1 end SUBROUTINE rdm(io,N,M,ic) IMPLICIT none integer*4 N,N1,N3,io,LN,J,ic real*8 M(N,N) ic=ic+1 read(io,*) N1=1 1 N3=min(N1+4,N) 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.N)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(9,*)V(LN),V(LN) return end FUNCTION GETDET(A,N) c determinant of a matrix IMPLICIT none integer*4 N,I,J,K logical DETEXISTS real*8 A(N,N) real*16 M,TEMP,L,GETDET real*16, allocatable::ELEM(:,:) allocate(ELEM(N,N)) DO 1 I=1,N DO 1 J=1,N c ELEM(I,J)=QEXTD(A(I,J)) 1 ELEM(I,J)=real(A(I,J),16) DETEXISTS=.TRUE. L=1.0q0 DO 2 K=1,N-1 IF(ABS(ELEM(K,K)).LE.1.0q-20)THEN DETEXISTS=.FALSE. DO 3 I=K+1,N IF(ELEM(I,K).NE.0.0q0)THEN DO 4 J=1,N TEMP=ELEM(I,J) ELEM(I,J)=ELEM(K,J) 4 ELEM(K,J)=TEMP DETEXISTS=.TRUE. L=-L EXIT ENDIF 3 CONTINUE IF (DETEXISTS.EQV..FALSE.)THEN GETDET = 0.0q0 RETURN ENDIF ENDIF DO 2 J=K+1,N M=ELEM(J,K)/ELEM(K,K) DO 2 I=K+1,N 2 ELEM(J,I)=ELEM(J,I)-M*ELEM(K,I) GETDET =L DO 5 I=1,N 5 GETDET=GETDET*ELEM(I,I) RETURN END function jeje(np,nx,it,NN,ee,np0,LEXCL) implicit none integer*4 np,nx,it(*),NN(*),np0,LEXCL,ee(np0,LEXCL),je,jeje,jj,ie je=0 do 104 jj=1,np if(nx.ne.NN(jj))goto 104 do 1041 ie=1,NN(jj) 1041 if(it(ie).ne.ee(jj,ie))goto 104 je=jj goto 1042 104 continue 1042 jeje=je return end subroutine digest(np,ie,np0,LEXCL,p,NN,it,iexc,T) implicit none integer*4 je,np,np0,LEXCL,NN(*),it(*),ie(np0,LEXCL),jj,iexc, 1jeje real*8 T,p(*) c does it exist already within 1 ... np?: je=jeje(np,iexc,it,NN,ie,np0,LEXCL) if(je.ne.0)then c the term already exist in the expansion as je^th - just add it p(je)=p(je)+T else c term not found - add as new term np=np+1 if(np.gt.np0)call report('Too many terms') p(np)=T NN(np)=iexc do 1 jj=1,iexc 1 ie(np,jj)=it(jj) endif return end function FC(fac,si,Nexc,np0,LEXCL,iwr,C,D,N) c FC = Franc Condon factor for <0|si> c fac = <0|0*> c np0 ... dimension of working buffer implicit none integer*4 si(*),Nexc,np,iex,ii,nu,jj,ip,N,ic,jold, 1nuj,jc,kk,iwr,LEXCL,np0,iq,jq real*8 FC,fac,D(*),C(N,N),pini real*8,allocatable::p(:) integer*4,allocatable::NN(:),ie(:,:),it(:) np=1 allocate(p(np0),NN(np0),ie(np0,LEXCL),it(np0)) p(np)=1.0d0 NN(np)=NExc c if(NExc.eq.2.and.si(1).eq.1.and.si(2).eq.12)then c write(6,*)'D1',D(1),' D12',D(12),'C1 12',C(1,12) c iwr=2 c else c iwr=1 c endif if(iwr.gt.1)write(6,604)(si(iex),iex=1,NExc) 604 format(/,'<0|e>, e = ',10i3,' requested:') do 102 iex=1,NN(np) 102 ie(np,iex)=si(iex) c expand reccurent formula into a sum and shrink back to <0|0> 777 do 101 ii=1,np c term ii - reduce excitation one by one: pini=p(ii) if(iwr.gt.1)write(6,607)ii,p(ii),(ie(ii,iex),iex=1,NN(ii)) 607 format(' term ',i2,':',f10.6,10i3) do 101 iex=1,NN(ii) iq=ie(ii,iex) if(iq.gt.0)then c <0| si_nu> reduce to <0| i-1>, <0|i-2>, <0|i-1,j-1>(j<>i) nu=0 do 1032 jj=1,NN(ii) 1032 if(ie(ii,jj).eq.iq)nu=nu+1 c write term <0|v'-1_iex> into string it: ip=0 ic=0 do 1031 jj=1,NN(ii) if(iq.eq.ie(ii,jj).and.ic.lt.1)then ic=ic+1 else ip=ip+1 it(ip)=ie(ii,jj) endif 1031 continue c call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-1, 1 pini*D(iq)/dsqrt(dble(2*nu))) if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-1) 606 format('digest ',10i3) if(nu.gt.1)then c write term <0|v'-2_iex> into string it: ic=0 ip=0 do 1033 jj=1,NN(ii) if(ie(ii,jj).eq.iq.and.ic.lt.2)then ic=ic+1 else ip=ip+1 it(ip)=ie(ii,jj) endif 1033 continue call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-2, 1 pini*dsqrt(dble(nu-1)/dble(nu))*C(iq,iq)) if(iwr.gt.1)write(6,606)(it(jj),jj=1,NN(ii)-2) endif jold=0 do 106 jj=1,NN(ii) jq=ie(ii,jj) if(jq.ne.iq.and.jq.ne.jold)then nuj=0 do 1034 kk=1,NN(ii) 1034 if(jq.eq.ie(ii,kk))nuj=nuj+1 c write term <0|v'-1_iex-1_j, j<>iex> into it: ic=0 jc=0 ip=0 do 1035 kk=1,NN(ii) if(ie(ii,kk).eq.iq.and.ic.lt.1)then ic=ic+1 else if(ie(ii,kk).eq.jq.and.jc.lt.1)then jc=jc+1 else ip=ip+1 it(ip)=ie(ii,kk) endif endif 1035 continue call digest(np,ie,np0,LEXCL,p,NN,it,NN(ii)-2, 1 pini*dsqrt(dble(nuj)/dble(nu))*C(iq,jq)) if(iwr.gt.1)write(6,606)(it(kk),kk=1,NN(ii)-2) endif 106 jold=jq if(iwr.gt.1)then write(6,*)'np now ',np do 108 jj=1,np 108 write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj)) endif c eliminate the old term and start over do 105 jj=ii,np-1 p(jj)=p(jj+1) NN(jj)=NN(jj+1) do 105 kk=1,NN(jj) 105 ie(jj,kk)=ie(jj+1,kk) np=np-1 if(iwr.gt.1)then write(6,*)'np after elimination ',np do 107 jj=1,np 107 write(6,605)jj,NN(jj),p(jj),(ie(jj,kk),kk=1,NN(jj)) 605 format(i3,i2,g9.3,10i3) endif goto 777 endif 101 continue if(np.ne.1)call report('np <> 1') FC=p(1)*fac return end subroutine rdd(f,d,nat,nroot,u0,iwr,m0,a) c read transition dipole derivatives from excited state freq calc implicit none character*(*) f real*8 d(*),u(3),u0(*),step,sign,stepau,m0(*),a(*),m(3) integer*4 nat,ln,ia,xa,i,nroot,xar,iar,ibr,ib,ic,ii,ix,iwr character*80 s80 open(9,file=f,status='old') ln=0 iar=0 xar=0 ibr=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)-3 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 2 endif endif goto 1 88 close(88) call report('TD freq not found') 2 ia=1 xa=1 ib=0 ii=0 step=0.001d0 do 7 ix=1,9*nat a(ix)=0.0d0 7 d(ix)=0.0d0 4 read(9,900,end=78,err=78)s80 if(s80(2:14).eq.'Nuclear step=')then write(6,*)s80(1:23) read(s80(15:23),*)step endif if(s80(2:24).eq.'Re-enter D2Numr: IAtom=')then read(s80(25:27),*)iar read(s80(34:34),*)xar read(s80(42:43),*)ibr if(ia.ne.iar)call report('ia <> iar') if(xa.ne.xar)call report('xa <> xar') if(ib.ne.ibr)call report('ib <> ibr') endif if(s80(2:65).eq.'Ground to excited state transition electric dipol 1e moments (Au):')then if(ii.eq.0)then write(6,*)'Zero point' do 5 i=1,nroot 5 read(9,*) read(9,*)u0(1),(u0(ix),ix=1,3) else ib=ib+1 if(ib.gt.2)then ib=1 xa=xa+1 if(xa.gt.3)then xa=1 ia=ia+1 endif endif if(iwr.gt.1)write(6,*)' d ',ia,xa,ib,iar,xar,ibr do 6 i=1,nroot 6 read(9,*) read(9,*)u(1),(u(ix),ix=1,3) if(ib.eq.1)then sign=1.0d0 else sign=-1.0d0 endif do 8 ix=1,3 8 d(ix+3*(xa-1)+9*(ia-1))=d(ix+3*(xa-1)+9*(ia-1))+sign*u(ix) endif endif if(s80(2:65).eq.'Ground to excited state transition magnetic dipol 1e moments (Au):')then if(ii.eq.0)then write(6,*)'Zero point magnetic' do 51 i=1,nroot 51 read(9,*) read(9,*)m0(1),(m0(ix),ix=1,3) else if(iwr.gt.1)write(6,*)' m ',ia,xa,ib,iar,xar,ibr do 61 i=1,nroot 61 read(9,*) read(9,*)m(1),(m(ix),ix=1,3) do 81 ix=1,3 81 a(ix+3*(xa-1)+9*(ia-1))=a(ix+3*(xa-1)+9*(ia-1))+sign*m(ix) endif ii=ii+1 if(ia.eq.nat.and.ix.eq.3.and.ib.eq.2)goto 78 endif goto 4 78 close(9) stepau=step/0.529177d0 do 9 ix=1,9*nat d(ix)=d(ix)/(stepau+stepau) 9 a(ix)=a(ix)/(stepau+stepau) 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) 400 format(i5,i2,3g12.4) close(40) write(6,*)'DEG.TEN written' return end c ============================================================== subroutine rde(d,nat,u0,a,m0) c read transition dipole derivatives from file implicit none real*8 d(*),u0(*),a(*),m0(*) integer*4 nat,ia,xa,ii,ix open(40,file='DE.TEN') read(40,*) do 10 ia=1,nat do 10 xa=1,3 ii=3*(xa-1)+9*(ia-1) 10 read(40,*)(d(ix+ii),ix=1,2),(d(ix+ii),ix=1,3) read(40,*) read(40,*)u0(1),u0(1),u0(1),u0(2),u0(3) read(40,*) do 11 ia=1,nat do 11 xa=1,3 ii=3*(xa-1)+9*(ia-1) 11 read(40,*)(a(ix+ii),ix=1,2),(a(ix+ii),ix=1,3) read(40,*) read(40,*)m0(1),m0(1),m0(1),m0(2),m0(3) close(40) return end c ============================================================== SUBROUTINE READS(N3,S,E,NQ) IMPLICIT none integer*4 N3,NQ,N,nat3,nat,i,J,ix,iz real*8 S(N3,N3),E(*),CM,SFAC CM=219470.0d0 SFAC=0.0234280d0 C 1 a.u.= 2.1947E5 cm^-1 open(4,file='F.INP') read(4,*)NQ,nat3,nat N=NAT3 do 1 i=1,NAT 1 read(4,*) read(4,*) DO 2 I=1,NAT DO 2 J=1,NQ 2 read(4,*)(s(3*(i-1)+ix,J),ix=1,2), 1(s(3*(i-1)+ix,J),ix=1,3) read(4,*) READ(4,4000)(E(I),I=1,NQ) 4000 FORMAT(6F11.6) close(4) c c delete zero and negative modes: iz=0 66 do 6 i=1,NQ if(E(i).lt.0.1d0)then do 7 j=i,NQ-1 E(j)=E(j+1) do 7 ix=1,nat3 7 s(ix,j)=s(ix,j+1) do 8 ix=1,nat3 8 s(ix,NQ)=0.0d0 E(NQ)=0.0d0 iz=iz+1 NQ=NQ-1 goto 66 endif 6 continue if(iz.gt.0)write(6,*)iz,' zero and negative modes deleted' write(6,*)NQ,' vibrational modes considered' DO 3 I=1,NQ 3 E(I)=E(I)/CM DO 5 I=1,N DO 5 J=1,N 5 S(I,J)=S(I,J)*SFAC c now SMS=1 if masses are in au RETURN end c ============================================================== subroutine trafog(dd,ddi,nat,NQ) implicit none integer*4 nat,NQ,ix,iq,ia,xa real*8 dd(*),ddi(*) real*8,allocatable::e(:),s(:,:) allocate(e(3*nat),s(3*nat,3*nat)) call READS(3*nat,s,e,NQ) do 1 ix=1,9 do 1 iq=1,NQ ddi(ix+9*(iq-1))=0.0d0 do 1 ia=1,nat do 1 xa=1,3 c s-opposite order of normal modes: 1 ddi(ix+9*(iq-1))=ddi(ix+9*(iq-1))+dd(ix+9*(xa-1)+27*(ia-1)) 1*s(xa+3*(ia-1),NQ-iq+1) write(6,*)'Gnj transformed to normal modes' return end c ============================================================== subroutine trafo(dd,ddi,nat,NQ) implicit none integer*4 nat,NQ,ix,iq,ia,xa real*8 dd(*),ddi(*) real*8,allocatable::e(:),s(:,:) allocate(e(3*nat),s(3*nat,3*nat)) call READS(3*nat,s,e,NQ) do 1 ix=1,3 do 1 iq=1,NQ ddi(ix+3*(iq-1))=0.0d0 do 2 ia=1,nat do 2 xa=1,3 c s-opposite order of normal modes: 2 ddi(ix+3*(iq-1))=ddi(ix+3*(iq-1))+dd(ix+3*(xa-1)+9*(ia-1)) 1*s(xa+3*(ia-1),NQ-iq+1) 1 continue write(6,*)'dipoles transformed to normal modes' return end c ============================================================== subroutine zm(M,N) implicit none integer*4 N,i,j real*8 M(N,N) do 1 i=1,N do 1 j=1,N 1 M(i,j)=0.0d0 return end c ============================================================== subroutine zv(J,N) implicit none real*8 J(*) integer*4 N,i do 1 i=1,N 1 J(i)=0.0d0 return end c ============================================================== subroutine wrgnj(i,g,d) implicit none integer*4 i,ix,iy real*8 g(*),d(*) open(43,file='GNJ.TEN') write(43,431)i 431 format(i5,' transition, Gnj tensor:') do 1 ix=1,3 1 write(43,4343)(g(ix+3*(iy-1)),iy=1,3) c ix index - electric 4343 format(3G14.6) write(43,432)i 432 format(i5,' transition electric dipole:') write(43,4343)(d(ix),ix=1,3) write(43,433)i 433 format(i5,' transition magnetic dipole:') write(43,4343)(d(ix),ix=7,9) close(43) return end c ============================================================== subroutine rdgnj(g,d,m) implicit none integer*4 ix,iy real*8 g(*),d(*),m(*) open(43,file='GNJR.TEN') read(43,*) do 1 ix=1,3 1 read(43,*)(g(ix+3*(iy-1)),iy=1,3) read(43,*) read(43,*)(d(ix),ix=1,3) read(43,*) read(43,*)(m(ix),ix=1,3) close(43) write(6,*)' Gnj0 d0 m0 read from GNJR.TEN' return end c ============================================================== subroutine wrgn(gnji,NQ,nroot) implicit none integer*4 NQ,i,ix,iy,nroot real*8 gnji(*) open(9,file='GNJDI.TEN') write(9,900)nroot 900 format(i3,' transition, derivatives of MCD tensor') do 1 i=1,NQ write(9,901)i 901 format(i5,' mode') do 1 ix=1,3 1 write(9,902)(gnji(ix+3*(iy-1)+9*(i-1)),iy=1,3) 902 format(3F12.6) close(9) end c ============================================================== subroutine wrda(ddi,aai,NQ,nroot) implicit none integer*4 NQ,i,ix,nroot real*8 ddi(*),aai(*) open(9,file='DEI.TEN') write(9,900)nroot 900 format(i3,' transition, derivatives of electric dipole',/, 1' mode ux uy uz') do 1 i=1,NQ 1 write(9,901)i,(ddi(ix+3*(i-1)),ix=1,3) 901 format(i5,3F12.6) write(9,902)nroot 902 format(i3,' transition, derivatives of magnetic dipole',/, 1' mode mx my mz') do 2 i=1,NQ 2 write(9,901)i,(aai(ix+3*(i-1)),ix=1,3) close(9) end