program guvcde c c polarizability from TDDFT by SOS c implicit none character*80 fo,s80,output integer n,ie,i,j,k,l,np,nmo,ix,ia,nat,grid(3),iep, 1oa,ob,idum,io,iz,nmo2,na,nb,nd0,iy, 2nt,nst,ifix,iclose,igiao parameter (nd0=21) 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 real*8 uux,uuy,uuz,DD,xp,yp,zp,px,py,pz,ddu(3),p0,t, 1uuxi,uuyi,uuzi,uuxr,uuyr,uuzr,p1,p2,p3,de1,de2,de3, 1ucxi,ucyi,uczi,ucxr,ucyr,uczr,dg1,dg2,dg3, 1al(3,3),ax,ay,az,qua,wden,sum,wmin,wmax,dw,ali(3,3),gammaau, 1w,fr,ds1,rs,angle,fi,bohr,djg(nd0),fkj,fkg,ejg,fjg,gammanm, 1rx,ry,rz,dx,dy,dz,alit(3,3), 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(3),xy(3),yy(3),xa(3),ya(3),ar(3),reg(3),qxx,qxy,qxz,qyy,qyz, 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, 3ttl(9),cab,ccd,wau 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 real*8, allocatable::cij(:),rrx(:),orb(:),rry(:),rrz(:),ror(:), 1rix(:),riy(:),riz(:),roi(:), 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(:) 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(:) character*1 s1 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, 1lzmat,ldeg,lwrt,DD,degper0,deglim,lnorm,nst,ifix,nmo,igiao, 2gammaau,wmin,wmax,np,lgamma,lort,gammanm,lgnm,lgau) c call dimensions(n,fo,nmo,nmo2,nat,lzmat) write(6,*)n,' transitions' 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 allocate(cij(2*nmo2*n),aij(nmo2*n),dij(nmo2*n), 1cijb(nmo2*n),aijb(nmo2*n),bijb(nmo2*n),nib(n), 1bij(nmo2*n),ni(2*n),daij(nmo2*n),dbij(nmo2*n),nid(n), 2dl(n,3),dv(n,3),r(n,3),e(n),ev(n),r0v(n),r0l(n),eau(n), 3x(3,nat),xau(3,nat),qq(nat),qv(n,3,3),ut(3,n), 4dlrec(n,3),dvrec(n,3),rrec(n,3)) call readfile(lzmat,qq,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) 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(nmo2*(i-1)+j).gt.nmo.or.aij(nmo2*(i-1)+j).lt.1)then write(6,*)' i:',i,' j:',j,' aij:',aij(nmo2*(i-1)+j),'nmo:',nmo call report('wrong orbital a') endif 141 if(bij(nmo2*(i-1)+j).gt.nmo.or.bij(nmo2*(i-1)+j).lt.1) 1 call report('wrong orbital b') do 14 j=1,nid(i) if(daij(nmo2*(i-1)+j).gt.nmo.or.daij(nmo2*(i-1)+j).lt.1)then write(6,*)' i:',i,' j:',j,' dij:',daij(nmo2*(i-1)+j) call report('wrong orbital da') endif 14 if(dbij(nmo2*(i-1)+j).gt.nmo.or.dbij(nmo2*(i-1)+j).lt.1) 1 call report('wrong orbital db') 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(n,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 tttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt c c computation of polarizabilities: 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='SOS.TTT') c velocity-based trace: open(9,file='spav.prn') c dipole-based trace: open(91,file='spal.prn') write(44,446)np,wmin,wmax 446 format(i6,2f12.2) write(44,443) 443 format(' Dynamic polarizabilities') write(44,444) 444 format(10x,' Alpha ( f ) Alpha ( g) ', 1/,'l(nm) ux uy uz') 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 do 101 ie=1,n call setv(eau(ie),w,fr,fi,gammaau,gammanm,lgnm,0,lgau) c a_ab = 2 sum wjn (1/(wjn**2-w**2)) Re ua ub c length: t=dl(ie,i)*dl(ie,j) al(i,j) =al(i,j) +t*fr ali(i,j)=ali(i,j)+t*fi c velocity: t=dv(ie,i)*dv(ie,j) 101 alit(i,j)=alit(i,j)+t*fi 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) write(44,442) (al(2,j),j=1,3),(ali(2,j),j=1,3) 102 write(44,442) (al(3,j),j=1,3),(ali(3,j),j=1,3) 441 format(f8.2,6f10.2) 442 format(8x ,6f10.2) close(9) write(6,*)' spav.prn (velocity) written' close(91) write(6,*)' spal.prn (length ) written' write(44,445) 445 format(10x,' Gp/w ( f) Gp/w ( g )', 1/,'l(nm) ux uy uz') open(9 ,file='spgv.prn') open(91,file='spgl.prn') 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,1,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,*)' spgv.prn written' write(6,*)' spgl.prn written' 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,0,lgau) al(i,j) =al( i,j) +dl(ie,l)*qua*fr 106 ali(i,j)=ali(i,j) +dl(ie,l)*qua*fi 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) write(6,*)' SOS.TTT written' c c recalculate dipole strengths c dddddddddddddddddddddddddddddddddddddddddddddddddddddddddddd if(lrds)then 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:' call zerstat(xa,ya,xx,yy,xy,ar,3) 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) do 312 j=1,nd0 312 djg(j)=-djg(j) djg(4)=-djg(4)/eau(ie) djg(5)=-djg(5)/eau(ie) djg(6)=-djg(6)/eau(ie) 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) 313 rrec( ie,j)=djg(6+j) 301 write(6,3006)e(ie),(djg(j),j=1,9), 1 (dl(ie,j),j=1,3),(dv(ie,j),j=1,3),(r(ie,j),j=1,3) 3006 format(f6.1,9f8.4,/,6x,9f8.4) do 107 j=1,3 107 call stn(n*3,ar(j),xa(j),ya(j),reg(j),yy(j),xy(j)) write(6,3009)(cc(xy(j),n3,xa(j),ya(j),xx(j),yy(j)),j=1,3), 1 (reg(j),j=1,3),(ar(j),j=1,3) 3009 format(/,' Correlation coefficients:',/, 1 ' electric:',f12.6,' gradient:',f12.6,' magnetic:',f12.6,/, 1 ' Regression:',/, 1 ' ',f12.6,' ',f12.6,' ',f12.6,/, 1 ' :',/, 1 ' ',f12.6,' ',f12.6,' ',f12.6,/) 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(lmcd)then c nuclear dipole moment: do 311 j=1,3 un(j)=0.0d0 do 311 l=1,nat 311 un(j)=un(j)+xau(j,l)*qq(l) 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') else write(s1,'(i1)')igiao call inispec(37,'MCD4.TAB','MCD IGIAO='//s1) endif 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 allocate(hd(n*n*6),ut6(6,n),pre6(6*nmo2),pred6(6*nmo2), 1 dlg(6*n)) write(6,*)'pre-calculating , x=r,grad' 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) c diagonal dipole and gradient: call getdipole(ut6,n,cij,nmo,nmo2,ni,aij,bij, 1 dij,nid,daij,dbij,lopen,pre6,pred6,ie,6,6) 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) do 3051 j=1,6 c : dlg(j+6*(ie-1))=djg(j) c : 3051 hd(j+(ie+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,nmo2,cij,pre6,lopen,daij, 1 dbij,pred6,nmo,dij,nid) do 3052 j=1,6 3052 hd(j+(ie+n*(iep-1))*6)=djk(j) endif 305 continue deallocate(pre6,pred6) 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 the whole tensor J. Chem. Theory COmput 2013, 9, 1557, equation (9): call vz(ttl,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) c calculate state dipole: call getdipole(ut,n,cij,nmo,nmo2,ni,aij,bij, 1 dij,nid,daij,dbij,lopen,pre,pred,ie,nd0,3) sum=0.0d0 do 317 i=1,ni(ie) 317 sum=sum+cij(nmo2*(ie-1)+i)**2 do 326 i=1,nid(ie) 326 sum=sum+dij(nmo2*(ie-1)+i)**2 do 328 i=1,nib(ie) 328 sum=sum-cijb(nmo2*(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) do 320 ia=4,6 320 djg(ia)=-djg(ia)/eau(ie) 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) c djk===> call dodjk(djk,nd0,ni,iep,aij,bij,nmo2,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 c c contribution to the MCD tensor c /ekg + /ekj do 1 ix=1,3 do 1 iy=1,3 1 ttl(iy+3*(ix-1))=ttl(iy+3*(ix-1))+ 1 djk(ix)*dkg(6+iy)*fkg+dkg(ix)*djk(6+iy)*fkj 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 compansation 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: org1 =org1 +rx *u1kg+ry *u2kg+rz *u3kg c x/(Ek-El)./2: org11=org11-brx*u1kg-bry*u2kg-brz*u3kg c write(6,6038)'kk',ie,iep,ekg,dkg(9),rz/fkg,u3kg c seems to compensat eperfectly c x/(EK-EL)./2: org2=org2+qx*u1kj+qy*u2kj+qz*u3kj c x/(El-Ej)./2: org22=org22-qqx*u1kj-qqy*u2kj-qqz*u3kj c write(6,6038)'qq',ie,iep,ekg,djk(9),qz/fkj,u3kj c seems to compensate perfectly endif 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 (-)/ejg do 2 ix=1,3 do 2 iy=1,3 2 ttl(iy+3*(ix-1))=ttl(iy+3*(ix-1))+ddu(ix)*djg(6+iy)*fjg 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 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 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,((ttl(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,f8.1) 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 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=40 else j=37 endif do 401 i=28,j write(i,3002) 401 close(i) 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 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 allocate(wr(nat),wi(nat)) do 2006 i=1,nat wi(i)=0.0d0 2006 wr(i)=0.0d0 endif inquire(file='orbitals.cub',exist=lcub) if(.not.lcub)call report('orbitals.cub not present') open(8,file='orbitals.cub') open(9,file='roi.cub') open(91,file='ror.cub') read(8,2000)s80 2000 format(a80) write(9,2000)s80 write(91,2000)s80 read(8,*) write(9,901) write(91,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 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 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 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 do 902 ia=1,nat read(8,2000)s80 read(s80(18:80),*)(x(i,ia),i=1,3) write(9 ,2000)s80 902 write(91,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))) xp=ax-px do 906 ix=1,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 read(8,4000)((orb(nmo*(iz-1)+io),io=1,nmo),iz=1,grid(3)) c do 907 ie=1,n call setv(eau(ie),wden,fr,fi,gammaau,gammanm,lgnm,0,lgau) uux=dv(ie,1)*eau(ie) uuy=dv(ie,2)*eau(ie) uuz=dv(ie,3)*eau(ie) uuxi=uux*fi uuyi=uuy*fi uuzi=uuz*fi uuxr=uux*fr uuyr=uuy*fr uuzr=uuz*fr do 907 j=1,ni(ie) ucxi=cij(nmo**2*(ie-1)+j)*uuxi ucyi=cij(nmo**2*(ie-1)+j)*uuyi uczi=cij(nmo**2*(ie-1)+j)*uuzi ucxr=cij(nmo**2*(ie-1)+j)*uuxr ucyr=cij(nmo**2*(ie-1)+j)*uuyr uczr=cij(nmo**2*(ie-1)+j)*uuzr oa=aij(nmo**2*(ie-1)+j) ob=bij(nmo**2*(ie-1)+j) do 907 iz=1,grid(3) rrx(iz)=rrx(iz)+ucxr*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) rry(iz)=rry(iz)+ucyr*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) rrz(iz)=rrz(iz)+uczr*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) rix(iz)=rix(iz)+ucxi*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) riy(iz)=riy(iz)+ucyi*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) 907 riz(iz)=riz(iz)+uczi*orb(nmo*(iz-1)+oa)*orb(nmo*(iz-1)+ob) c zp=az-pz do 908 iz=1,grid(3) zp=zp+pz ror(iz)=dsqrt(2.0d0*(rrx(iz)**2+rry(iz)**2+rrz(iz)**2)) roi(iz)=dsqrt(2.0d0*(rix(iz)**2+riy(iz)**2+riz(iz)**2)) 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: wr(iclose)=wr(iclose)+ror(iz) wi(iclose)=wi(iclose)+roi(iz) endif 908 continue write(9 ,4000)(roi(iz),iz=1,grid(3)) 906 write(91,4000)(roi(iz),iz=1,grid(3)) 4000 format(6E13.5) if(lw)then c normalize atomic weights: sumr=0.0d0 sumi=0.0d0 do 2008 i=1,nat sumi=sumi+wi(i) 2008 sumr=sumr+wr(i) do 2009 i=1,nat if(sumi.ne.0.0d0)wi(i)=wi(i)/sumi 2009 if(sumr.ne.0.0d0)wr(i)=wr(i)/sumr write(6,6071)(100.0d0*wr(i),i=1,nat) write(6,6071)(100.0d0*wi(i),i=1,nat) 6071 format(12f6.1) c recalculate to get rid of small contributions: allocate(ind(nat),jnd(nat)) open(50,file='WEIGHTS.TXT') write(6 ,*)'Real weights:' write(50,*)'Real weights:' write(6 ,6071)(100.0d0*wr(i),i=1,nat) write(50,6071)(100.0d0*wr(i),i=1,nat) write(6 ,*)'Imaginary weights:' write(50,*)'Imaginary weights:' write(6 ,6071)(100.0d0*wi(i),i=1,nat) write(50,6071)(100.0d0*wi(i),i=1,nat) call sort(wr,ind,jnd,nat) call trimv(wr,nat,nt) write(6 ,*)'Real weights trimmed:' write(50,*)'Real weights trimmed:' write(6 ,6071)(100.0d0*wr(jnd(i)),i=1,nat) write(50,6071)(100.0d0*wr(jnd(i)),i=1,nat) call sort(wi,ind,jnd,nat) call trimv(wi,nat,nt) write(6 ,*)'Imaginary weights trimmed:' write(50,*)'Imaginary weights trimmed:' write(6 ,6071)(100.0d0*wi(jnd(i)),i=1,nat) write(50,6071)(100.0d0*wi(jnd(i)),i=1,nat) close(50) deallocate(ind,jnd) endif deallocate(orb,rrx,rry,rrz,rix,riy,riz,ror,roi) close(8) close(9) write(6,*)' roi.cub ror.cub written' endif c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp stop end subroutine trimv(v,nat,nt) implicit none real*8 v(*),t integer*4 nat,i,nt t=0.0d0 do 1 nt=1,nat t=t+v(nt) 1 if(t.ge.0.5d0)goto 2 2 if(nt.gt.nat)nt=nat do 4 i=1,nt 4 v(i)=v(i)/t do 3 i=nt+1,nat 3 v(i)=0.0d0 return end subroutine sort(v,ind,jnd,nat) implicit none real*8 v(*),t integer*4 ind(*),nat,i,j,it,jnd(*) 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).gt.v(i))then t=v(i) v(i)=v(j) v(j)=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,ic,lgau) implicit none real*8 w,jm,fr,fi,g,wau,eau,gnm,a,pi,ln2 integer*4 ic logical lgnm,lgau c nm to au:= wau=(1.0d7/w)/219474.63d0 c set bandwidth, either constant in wavelength or energy: if(lgnm)then if(w.gt.gnm)then a=1.0d7/219474.63d0*gnm/(w**2-gnm**2) else a=gnm*wau/w endif else a=g endif jm=((eau+wau)*(eau-wau))**2+(a*eau)**2 fr= (eau**2-wau**2) / 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 if(ic.eq.0)then fr= fr*eau fi= fi*eau 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 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,n2 ,ni,aij,bij, 1dij,nid,daij,dbij,lopen,pre,pred,ie,nd0,ic) 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) 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 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 do 2 i=1,ni(ie) oa=aij(nmo2*(ie-1)+i) ob=bij(nmo2*(ie-1)+i) ciji=cij(nmo2*(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) 1pre(j+(oa-1+(od-1)*nmo)*nd0)= 2pre(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 if(oc.ne.oa) 1pre(j+(oc-1+(ob-1)*nmo)*nd0)= 2pre(j+(oc-1+(ob-1)*nmo)*nd0)+ciji*bp(oa+nmo*(oc-1)+jj) 2 continue if(lopen)then do 42 i=1,nid(ie) oa=daij(nmo2*(ie-1)+i) ob=dbij(nmo2*(ie-1)+i) ciji=dij(nmo2*(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))) endif return end subroutine readoptions(lden,lmcd,lrds,lw,wden,fo, 1lzmat,ldeg,lwrt,DD,degper0,deglim,lnorm,nst,ifix,nmo,igiao, 1gammaau,wmin,wmax,np,lgamma,lort,gammanm,lgnm,lgau) implicit none character*(*) fo logical lden,lmcd,lrds,lw,lzmat,ldeg,lwrt,lnorm,lgamma,lort, 1lgnm,lgau real*8 DD,degper0,deglim,wden,wmin,wmax,gammaau,gammanm character*4 key integer*4 nst,ifix,nmo,igiao,np degper0=0.0001d0 deglim=0.00001d0 lden=.false. ifix=3 nst=0 lnorm=.false. lort=.true. lgnm=.true. lgamma=.false. lgau=.false. lmcd=.true. 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 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 extended writing if(key.eq.'LWRT')read(9,*)lwrt 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:2).eq.'NP')read(9,*)np 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 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) return end subroutine dimensions(n,fo,nmo,nmo2,nat,lzmat) c returns number of transitions and atoms implicit none integer*4 n,nat,l,I,q,nmo,nmo2 real*8 axdum character*80 s80 character*70 st character*(*) fo logical lzmat n=0 nat=0 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) implicit none character*80 s80 character*75 st character*(*) fo logical lzmat,lden,lopen,lnorm,lort,lwrt 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 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 real*8,allocatable::obaj(:),obak(:) real*8 ska 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(8:22).eq.'alpha electrons')then read(s80( 1: 6),*)na read(s80(25:31),*)nb write(6,2000)s80 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 if(lden)then ni(ie)=0 nib(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 read(s80( 1: 7),*)aij(nmo2*(ie-1)+ni(ie)) read(s80(12:14),*)bij(nmo2*(ie-1)+ni(ie)) read(s80(18:30),*)cij(nmo2*(ie-1)+ni(ie)) call setabmax(aij(nmo2*(ie-1)+ni(ie)), 1 bij(nmo2*(ie-1)+ni(ie)),cij(nmo2*(ie-1)+ni(ie)), 2 amax,bmax,cmax) else if(s80(8:8).eq.'B')then c open shell, beta: nid(ie)=nid(ie)+1 lopen=.true. read(s80( 1: 7),*)daij(nmo2*(ie-1)+nid(ie)) read(s80(12:14),*)dbij(nmo2*(ie-1)+nid(ie)) read(s80(18:30),*) dij(nmo2*(ie-1)+nid(ie)) call setabmax(daij(nmo2*(ie-1)+ni(ie)), 1 dbij(nmo2*(ie-1)+ni(ie)),dij(nmo2*(ie-1)+ni(ie)), 2 amax,bmax,cmax) else c closed shell: ni(ie)=ni(ie)+1 read(s80( 1: 8),*) aij(nmo2*(ie-1)+ni(ie)) read(s80(12:15),*) bij(nmo2*(ie-1)+ni(ie)) read(s80(18:30),*) cij(nmo2*(ie-1)+ni(ie)) cij(nmo2*(ie-1)+ni(ie))=cij(nmo2*(ie-1)+ni(ie))*tq call setabmax(aij(nmo2*(ie-1)+ni(ie)), 1 bij(nmo2*(ie-1)+ni(ie)),cij(nmo2*(ie-1)+ni(ie)), 2 amax,bmax,cmax) sum1=sum1+cij(nmo2*(ie-1)+ni(ie))**2 endif endif goto 111 endif if(s80(10:11).eq.'<-'.and..not.lopen)then nib(ie)=nib(ie)+1 if(s80(8:8).eq.'A')then else if(s80(8:8).eq.'B')then else 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 c experimental fix for the deexitations, closed shell only: if(ifix.eq.1)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 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 endif 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' write(6,6009)amax,bmax,cmax 6009 format(' Dominant transition from',i4,' to ',i4,f6.3) if(.not.lopen)write(6,1234)sum1,sum2 1234 format('fnorm, bnorm:',2f10.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 normcij(ie,cij,ni,nmo2) 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' 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) 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 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) c element implicit none integer*4 ni(*),nid(*),j,oa,ob,aij(*),bij(*),daij(*),dbij(*), 1ie,i,nmo,nmo2,nib(*),aijb(*),bijb(*),ifix,nd0 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(nmo2*(ie-1)+i) ob=bij(nmo2*(ie-1)+i) 301 djg(j)=djg(j)+cij(nmo2*(ie-1)+i)*bp(ob+nmo*(oa-1)+nmo2*(j-1)) if(ifix.eq.2)then do 302 i=1,nib(ie) oa=aijb(nmo2*(ie-1)+i) ob=bijb(nmo2*(ie-1)+i) 302 djg(j)=djg(j)+cijb(nmo2*(ie-1)+i)*bp(oa+nmo*(ob-1)+nmo2*(j-1)) endif do 312 i=1,nid(ie) oa=daij(nmo2*(ie-1)+i) ob=dbij(nmo2*(ie-1)+i) 312 djg(j)=djg(j)+dij(nmo2*(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) implicit none integer*4 a,b,amax,bmax real*8 c,cmax if(dabs(c).gt.cmax)then cmax=dabs(c) amax=a bmax=b endif return end subroutine dodjk(djk,nd0,ni,iep,aij,bij,nmo2,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(*),nmo2, 1daij(*),dbij(*),nid(*),nmo do 309 j=1,nd0 djk(j)=0.0d0 do 308 ip=1,ni(iep) oc=aij(nmo2*(iep-1)+ip) od=bij(nmo2*(iep-1)+ip) 308 djk(j)=djk(j)+cij(nmo2*(iep-1)+ip)*pre(j+(oc-1+(od-1)*nmo)*nd0) if(lopen)then do 3071 ip=1,nid(iep) oc=daij(nmo2*(iep-1)+ip) od=dbij(nmo2*(iep-1)+ip) 3071 djk(j)=djk(j) 1 +dij(nmo2*(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+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+n*(iel-1))*6 lk=(iel+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+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 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 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