program uworks c c various tasks associated with behind BO-computations: c c work 1 - open G.OUT and make derivatives of exc coefficients c work 2 - open GD.OUT and read the CP HF (KS) matrices c - needs sos output c work 3 - set up vibrational-electronic Hamiltonian implicit none integer noo,nw0,nl0,nst parameter (noo=100,nl0=100,nw0=100) c icon - control options c iwork - works to do integer*4 icon(noo),work(nw0) c lopt - logical options logical lopt(noo),lzmat,lnorm,lort,lwrt real*8 bohr,ropt(noo) c fn - TDDFT output c fd - derivative output (CP) character*80 fn,fd call readopt(icon,bohr,work,lopt,nl0,noo,nw0,fn,fd,lzmat, 1lnorm,lort,lwrt,nst,ropt) c if(work(1).ne.0)call cij2(fn,lzmat,bohr,lnorm,lort, 1lwrt,nst,lopt) if(work(2).ne.0)call cp(fd,bohr) if(work(3).ne.0)call makeH(icon,lopt) if(work(4).ne.0)call dospec(icon,lopt,ropt) end c ============================================================ subroutine dospec(icon,lopt,ropt) implicit none integer*4 icon(*),ND,no,n,nocc,nvirt,nmo,i,nmo2,n3,ig,ie,ii,ix, 1j,m1,nm,nat,na,nb,nd0,ne1,ne2,nv1,nv2,ib,iv,k,nn parameter (nd0=21) real*8 ropt(*),kt,q,p,ux,uy,uz,mx,my,mz,Cc,ds,rs,w,ef real*8,allocatable::C(:,:),cij0(:),Uo(:),sdm(:),s(:), 1E(:),u(:),uu(:),JdK(:),T(:),dij(:),cijb(:) logical lopt(*),lopen common/important/na,nb,ND,lopen call w36(' ') call w36('dospec') call w36('^^^^^^') open(9,file='ENBO.SCR',form='unformatted') read(9)ND close(9) allocate(E(ND),C(ND,ND)) call rre('ENBO.SCR',E,ND) call rre2('CNBO.SCR',C,ND) write(6,*)'Energies read from ENBO.SCR' write(6,*)'Eigenvectors read from CNBO.SCR' open(9,file='CC0.SCR.TXT',status='old') read(9,*)no,n,nocc,nvirt,nmo allocate(cij0(2*no*n),cijb(2*no*n),dij(2*no*n)) read(9,900)(cij0(i),i=1,2*no*n) read(9,900)(dij(i),i=1,2*no*n) read(9,900)(cijb(i),i=1,2*no*n) 900 format(4d20.14) close(9) write(6,*)'CC0.SCR.TXT read' open(9,file='F.INP',status='old') read(9,*)nm,nat,nat n3=3*nat allocate(s(nm*n3)) do 4 i=1,nat 4 read(9,*) read(9,*) do 5 i=1,nat do 5 j=1,nm ii=3*(i-1)+n3*(j-1) 5 read(9,*)(s(ix+ii),ix=1,2),(s(ix+ii),ix=1,3) close(9) write(6,*)'F.INP read' ne1=icon(1) ne2=icon(2) nv1=icon(3) nv2=icon(4) if(ne1.eq.0)ne1=1 if(nv1.eq.0)nv1=1 if(ne2.eq.0)ne2=n if(nv2.eq.0)nv2=nm nn=ne2-ne1+1 allocate(JdK((nn+1)**2*nd0)) write(6,6001)no,nd0,((no+1)**2*nd0)/1000000 6001 format('dipoles ',i4,i4,i10,'MB') c , J,K .. TDDFT WF, ground or excited: call eldipoles(nmo,JdK,no,n,cij0,nocc,dij,cijb,ne1,ne2) nmo2=nmo*nmo allocate(Uo(nmo2*n3)) c a-coordinate derivative of MO_J = sum(k) U_A_J_k MO_k call readUou(n3,nmo,nocc,nmo2,Uo,icon) no=nocc*nvirt allocate(u(2*n3*no*n),uu(2*n3*no*n)) call readcij2(u,uu,m1) allocate(sdm(N*nmo2*7)) call vz(sdm,N*nmo2*7) call readsdm(sdm,nat,nmo) kT=1.9872041D-3*ropt(1)/627.5d0 q=0.0d0 allocate(T(nd0*(nn+1)*(nv2-nv1+2))) do 1 ig=1,ND ef=(E(ig)-E(1))/kT 1 if(ef.lt.20)q=q+dexp(-ef) write(6,*)'partition ',q open(9,file='U.TAB') write(9,901) 901 format(/,/,'------') do 2 ig=1,ND c , J .. TDDFT WF = Cig_kK call vz(T,nd0*(nn+1)*(nv2-nv1+2)) do 6 J=0,nn c loop over basis functions: TDDFT: ib=0 do 6 K=0,nn c loop over basis functions: vibrational do 6 iv=nv1-1,nv2 ib=ib+1 do 6 i=1,nd0 6 T(i+nd0*(iv+1-nv1)+nd0*(nv2-nv1+2)*J)= 1T(i+nd0*(iv+1-nv1)+nd0*(nv2-nv1+2)*J)+ 2C(ib,ig)*JdK(i+nd0*K+nd0*(nn+1)*J) ef=(E(ig)-E(1))/kT if(ef.lt.20)then p=dexp(-ef)/q else p=0.0d0 endif if(p.gt.ropt(2))then do 3 ie=1,ND if(ie.ne.ig)then ux=0.0d0 uy=0.0d0 uz=0.0d0 mx=0.0d0 my=0.0d0 mz=0.0d0 ib=0 do 7 J=0,nn do 7 iv=nv1-1,nv2 ib=ib+1 Cc=C(ib,ie) ii=nd0*(iv-nv1)+nd0*(nv2-nv1+2)*J ux=ux+T(1+ii)*Cc uy=uy+T(2+ii)*Cc uz=uz+T(3+ii)*Cc mx=mx+T(7+ii)*Cc my=my+T(8+ii)*Cc 7 mz=mz+T(9+ii)*Cc ds=(ux**2+uy**2+uz**2)*p rs=(mx**2+my**2+mz**2)*p if(ie.lt.ig)then ds=-ds rs=-rs endif w=dabs((E(ie)-E(ig))*219474.63d0) write(9,902)w,ds,rs,p,ig,ie 902 format(f15.2,3f10.4,i5,' ->',i6) endif 3 continue endif 2 continue write(9,903) 903 format(60(1h-)) close(9) 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 c ============================================================ subroutine wt(t,n) implicit none integer*4 n,i,ix,jx real*8 t(*) do 1 i=1,n write(6,*) do 1 ix=1,3 1 write(6,600)(t(jx+3*(ix-1)+9*(i-1)),jx=1,3) 600 format(3f8.4) return end c ============================================================ subroutine vz(x,n) implicit none integer*4 n,i real*8 x(*) do 1 i=1,n 1 x(i)=0.0d0 return end c ============================================================ subroutine rdp(x,U,N,M) implicit none integer*4 x,i,j,N,M,io,is,ie,i5,n5 real*8 U(*) io=N*M*(x-1) n5=M/5 if(M.gt.n5*5)n5=n5+1 is=-4 do 1 i5=1,n5 is=is+5 read(9,*) do 1 j=is,N ie=min(is+4,j) read(9,*)U(is+M*(j-1)+io),(U(i+M*(j-1)+io),i=is,ie) do 1 i=is,ie 1 U(j+M*(i-1)+io)=U(i+M*(j-1)+io) return end c ============================================================ subroutine wdur(U,N,M) implicit none integer*4 i,j,N,M,is,ie,i5,n5 real*8 U(*) n5=M/5 if(M.gt.n5*5)n5=n5+1 is=-4 do 1 i5=1,n5 is=is+5 ie=min(is+4,M) read(44,*) do 1 j=1,N 1 read(44,*)U(is+M*(j-1)),(U(i+M*(j-1)),i=is,ie) return end c ============================================================ subroutine rdu(x,U,N,M) implicit none integer*4 x,i,j,N,M,io,is,ie,i5,n5 real*8 U(*) io=N*M*(x-1) n5=M/5 if(M.gt.n5*5)n5=n5+1 is=-4 do 1 i5=1,n5 is=is+5 ie=min(is+4,M) read(9,*) do 1 j=1,N 1 read(9,*)U(is+M*(j-1)+io),(U(i+M*(j-1)+io),i=is,ie) c U(is+M*(j-1)+N*M*(x-1)) c is - expanded MO N .. number of virtual orbitals c j - expanding MO M .. number of occupied orbitals c x - for uR differentiation: 1 -ux c 2 -uy c 3 -uz c 4 -R1x c ... c x -Rlastz return end c ============================================================ subroutine readsdm(s,nat,nmo) implicit none real*8 s(*) integer*4 k,ia,nat,ix,I,N,M,J,L,nmo real*8,allocatable::tt(:,:) allocate(tt(nmo,nmo)) N=3*nat M=N*nmo open(66,file='SX.SCR.TXT',status='old') do 8 k=1,7 read(66,*) do 8 ia=1,nat do 8 ix=1,3 read(66,*) call rmtrio(66,tt,nmo) I=ix+3*(ia-1)+M*nmo*(k-1) do 8 J=1,nmo do 8 L=1,nmo 8 s(I+N*(J-1)+M*(L-1))=tt(J,L) c s(ix+3*(ia-1)+N*(J-1)+M*(L-1)+M*nmo*(k-1)) c ix - nuclear ix c ia - nucleus c M=N*nmo, N=3*nat c k: 1 - overlap c 2 -ux c 3 -uy c 4 -uz close(66) write(6,*)'SX.SCR.TXT read' return end c ============================================================ subroutine rmtrio(io,A,n) implicit none integer*4 n,N1,N3,io,LN,J real*8 A(n,n) N1=1 1 N3=MIN(N1+4,N) read(io,*) DO 130 LN=1,N 130 read(io,*)A(LN,N1),(A(LN,J),J=N1,N3) N1=N1+5 IF(N3.LT.N)GOTO 1 return end c ============================================================ 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 c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ 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 c ============================================================ subroutine readorb(key,e,cij,nmo,nao,io) implicit none integer*4 i,j,i1,j1,i2,nao,ic,io,nmo real*8 e(*),cij(*) logical number character*(*) key character*160 s160 i=1 read(io,*) 377 j=i+4 if(key(8:11).ne.'MOs:'.and.key(7:10).ne.'MOs:')read(io,*) if(j.gt.nmo)j=nmo c c to enable reading of various formats, use s160: read(io,160)s160 160 format(a160) c c add space before negative numbers: i1=22 3775 if(s160(i1:i1).eq.'-'.and.s160(i1-1:i1-1).ne.'E')then do 1 j1=160,i1+1,-1 1 s160(j1:j1)=s160(j1-1:j1-1) s160(i1:i1)=' ' i1=i1+1 endif i1=i1+1 if(i1.lt.160)goto 3775 read(s160(22:160),*)(e(i1),i1=i,j) do 39 i2=1,nao read(io,160)s160 c add space before negative numbers: ic=0 i1=22 3776 if(s160(i1:i1).eq.'-'. 1and.number(s160(i1-1:i1-1)).and.number(s160(i1+1:i1+1)))then ic=ic+1 write(6,*) s160(j1-1:j1-1) do 2 j1=160,i1+1,-1 2 s160(j1:j1)=s160(j1-1:j1-1) s160(i1:i1)=' ' i1=i1+1 endif i1=i1+1 if(i1.lt.160)goto 3776 if(ic.gt.0)then write(6,*)'Corrected:',ic write(6,160)s160 endif 39 read(s160(22:160),*,err=99)(cij(i1+(i2-1)*nmo),i1=i,j) if(j.lt.nmo)then i=i+5 read(io,*) goto 377 endif write(6,*)'MOs read' return 99 write(6,160)s160 write(6,*)i2,i,j call report('input error') end c ============================================================ subroutine wdu(U,N,M,io) implicit none integer*4 i,j,N,M,is,ie,i5,n5,io real*8 U(*) n5=M/5 if(M.gt.n5*5)n5=n5+1 is=-4 do 1 i5=1,n5 is=is+5 ie=min(is+4,M) write(io,440)(j,j=is,ie) 440 format(i18,4i14) do 1 j=1,N 1 write(io,441)j,(U(i+M*(j-1)),i=is,ie) 441 format(i7,5e14.6) return end c ============================================================ function number(s) character*(*) s logical number number=s.eq.'1'.or.s.eq.'2'.or.s.eq.'3'.or.s.eq.'4'.or.s.eq.'5' 1 .or.s.eq.'6'.or.s.eq.'7'.or.s.eq.'8'.or.s.eq.'9'.or.s.eq.'0' return end c ========================================================= subroutine chkcon(icon) integer*4 icon(*) if(icon(1).eq.0)call report('geometry not read') if(icon(2).eq.0)call report('number of occupied orbs missing') if(icon(3).eq.0)call report('number of orbs missing') if(icon(4).eq.0)write(6,*)'Magnetic perturbation missing' if(icon(5).eq.0)write(6,*)'Electric perturbation missing' if(icon(6).eq.0)call report('Nuclear differentiation missing') if(icon(7).eq.0)call report('cij missing') if(icon(8).eq.0)call report('UR missing') if(icon(9).eq.0)call report('P1 missing') if(icon(10).eq.0)call report('S1 missing') if(icon(11).eq.0)write(6,*)'F1 missing' return end c ========================================================= subroutine readopt(icon,bohr,work,lopt,nl0,noo,nw0,fn,fd,lzmat, 1lnorm,lort,lwrt,nst,ropt) implicit none character*(*) fn,fd character*15 key integer*4 nl0,noo,icon(*),work(*),IS,i,iw,nw0,nst,na,nb,ND real*8 bohr,ropt(*) logical lopt(*),lzmat,lnorm,lopen,lort,lwrt common/important/na,nb,ND,lopen do 1 i=1,noo ropt(i)=0.0d0 1 icon(i)=0 do 2 i=1,nw0 2 work(i)=0 do 3 i=1,nl0 3 lopt(i)=.false. bohr=0.529177d0 fn='G.OUT' fd='GD.OUT' lzmat=.true. lnorm=.true. lopen=.false. lort=.true. lwrt=.false. nst=0 ropt(1)=293.15d0 ropt(2)=1.0d-4 OPEN(2,FILE='UWORKS.OPT',STATUS='OLD',IOSTAT=IS) IF(IS.EQ.0)THEN WRITE(6,*)' FILE UWORKS.OPT FOUND ...' 2001 READ(2,1111,end=2002,err=2002)key write(6,1111)key 1111 FORMAT(A15) IF(key(1:4).EQ.'WORK')READ(2,*)iw,work(iw) IF(key(1:4).EQ.'FILE')READ(2,'(a)')fn IF(key(1:4).EQ.'FILD')READ(2,'(a)')fd IF(key(1:4).EQ.'LZMA')READ(2,*)lzmat IF(key(1:4).EQ.'NSTA')READ(2,*)nst IF(key(1:4).EQ.'LNOR')READ(2,*)lnorm IF(key(1:4).EQ.'LORT')READ(2,*)lort IF(key(1:4).EQ.'LWRT')READ(2,*)lwrt c interval of included electronic states NE1 ... NE2: IF(key(1:3).EQ.'NE1') READ(2,*)icon(1) IF(key(1:3).EQ.'NE2') READ(2,*)icon(2) c interval of included vibrational states NV1 ... NV2: IF(key(1:3).EQ.'NV1') READ(2,*)icon(3) IF(key(1:3).EQ.'NV2') READ(2,*)icon(4) IF(key(1:4).EQ.'IUCP')READ(2,*)icon(5) c project translations and rotations from the wavefunction c derivatives: IF(key(1:4).EQ.'PROJ')READ(2,*)lopt(1) IF(key(1:4).EQ.'TEMP')READ(2,*)ropt(1) IF(key(1:4).EQ.'GLIM')READ(2,*)ropt(2) GOTO 2001 2002 CLOSE(2) ENDIF c lopt(1) ... unique lopt(2)=lwrt return end c ========================================================= subroutine cij2(fd,lzmat,bohr,lnorm,lort,lwrt,nst,lopt) implicit none character*(*) fd integer*4 n,nmo,nmo2,nat,nocc,n3,nvirt,I,ierr,ia,ix,IC, 1na,nb,nst,no,ii,nn,l,j,k,n1,n2,m1,m2,nao,id1,id2,jd1,jd2, 1jm,km,ND logical lzmat,lexu,lden,lnorm,lopen,lort,lwrt,lmo,lopt(*) real*8 step,ff,bohr,s1,ad1,ad2,stol real*8,allocatable::r(:,:),r0(:,:),u(:),cmo(:),cmo0(:), 1uu(:),cij0(:),cij(:),e0(:),e(:),dl(:,:),emo0(:),emo(:), 1dv(:,:),eau(:),eau0(:),ev(:),ev0(:),dij(:),r0v(:),r0l(:), 1qv(:,:,:),m(:,:),cijb(:) integer*4,allocatable::moind(:) character*1 xyz(3),sinn(2) data xyz/'x','y','z'/,sinn/'+','-'/ common/important/na,nb,ND,lopen lden=.true. stol=0.5d0 call w36(' ') call w36('cij2') call w36('^^^^') inquire(file='CIJ2.SCR.TXT',exist=lexu) if(lexu)then call w36('CIJ2.SCR.TXT found on disk') return endif write(6,*)' Determining dimensions' call dimensions(n,fd,nmo,nmo2,nat,lzmat,nocc,m1,m2,na,nb) call wdims(6,nat,nmo,nocc,n,m1,m2,na,nb) open(9,file='DIMS.SCR.TXT') call wdims(9,nat,nmo,nocc,n,m1,m2,na,nb) close(9) write(6,*)' Dimensions written to DIMS.SCR.TXT' nvirt=nmo-nocc n3=3*nat no=nocc*nvirt nao=nmo if(n.eq.0.or.nmo.eq.0.or.nat.eq.0)return allocate(r(nat,3),r0(nat,3),moind(nmo),dij(2*no*n), 1cij(2*no*n),cij0(2*no*n),emo(nmo),emo0(nmo),cijb(2*no*n), 2dl(n,3),dv(n,3),m(n,3),e(n),ev(n),eau(n),r0v(n),r0l(n), 3qv(n,3,3),e0(n),ev0(n),eau0(n),cmo(nmo**2),cmo0(nmo**2)) allocate(u(2*n3*no*n),uu(2*n3*no*n)) write(6,6009)dble(4*n3*nocc*nvirt*n*8)/1000000.0d0 6009 format(' for derivatives allocated',f19.3,' MB',/) call vz(u,2*n3*nocc*nvirt*n) call vz(uu,2*n3*nocc*nvirt*n) call vz(dij,2*no*n) call vz(cij,2*no*n) call vz(cijb,2*no*n) call vz(cij0,2*no*n) open(2,file=fd) C c read zero geometry to r0: write(6,*) write(6,*)' Looking for equilibrium geometry ' do 4 ia=1,nat do 4 ix=1,3 4 r(ia,ix)=0.0d0 call chkgeop(r,nat,lzmat,step,I,ff,ierr,ia,ix,r0,0,bohr) c c read in MO coeffs call readcij(2,cmo0,nmo,nao,emo0,lmo,lwrt) c c read-in exc states write(6,*) write(6,*)' Looking for equilibrium states' call readCIEXC(nocc,nvirt,na,nb,n,nst, 1dl,dv,qv,m,r0v,r0l,e0,ev0,eau0,lden,cij0,lopen, 1dij,lnorm,lort,lwrt,cijb) write(6,*)' .... states read' write(6,*) open(9,file='EE.SCR.TXT') write(9,*)n do 20 i=1,n 20 write(9,900)i,eau0(i),e0(i) 900 format(i9,2f15.9) close(9) open(9,file='CC0.SCR.TXT') write(9,901)no,n,nocc,nvirt,nmo 901 format(4i9) write(9,902)(cij0(i),i=1,2*no*n) write(9,902)(dij(i),i=1,2*no*n) write(9,902)(cijb(i),i=1,2*no*n) 902 format(4d20.14) close(9) c N=3*nat c initialize coefficient derivatives c uu(1) = d Ci / d Ra = (ci(+) - ci(-)) / (2 del) c diagonal second are the second part in uu: c uu(2) = d^2 ci / d Ra dRb =(ci(+) + ci(-) - 2 ci(0)) / del^2 i=0 do 5 ia=1,nat do 5 ix=1,3 i=i+1 do 5 l=1,n ii=(l-1)*no nn=n*no+ii n1=ii+2*n*no*(i-1) n2=nn+2*n*no*(i-1) do 5 k=1,nvirt do 5 j=1,nocc c forward and backward uu(j+(k-1)*nocc+n1)=-2.0d0*cij0(j+(k-1)*nocc+ii) 5 uu(j+(k-1)*nocc+n2)=-2.0d0*cij0(j+(k-1)*nocc+nn) c IC=0 step=0.0d0 c 333333333333333333333333333333333333333333333333333333333333 3333 IC=IC+1 c c read geometry and check the diff. step: call chkgeop(r0,nat,lzmat,step,I,ff,ierr,ia,ix,r,1,bohr) if(ierr.eq.99)then write(6,*)'differentiation stoped, ',IC-1,' steps' goto 999 endif if(ff.gt.0.0d0)then l=1 else l=2 endif write(6,6363)IC,step*bohr,I,ia,xyz(ix),sinn(l) 6363 format('Displacement:',I5,' Step(A):',f10.4,' Coord.:',I5, 1' Atom:',I5,a1,1x,a1) if(ix+3*(ia-1).ne.I)then write(6,*)ix,ia,I call report('Inconsistent differentiation') endif c c read the orbitals call readcij(2,cmo,nmo,nao,emo,lmo,lwrt) c read-in exc states call readCIEXC(nocc,nvirt,na,nb,n,nst, 1dl,dv,qv,m,r0v,r0l,e,ev,eau,lden,cij,lopen, 1dij,lnorm,lort,lwrt,cijb) call mapcij(moind,cmo,cmo0,nmo,nao,lmo,emo,emo0) do 51 l=1,n ii=(l-1)*no nn=n*no+ii n1=ii+2*n*no*(i-1) n2=nn+2*n*no*(i-1) c s1=0.0d0 c dominant in equilibrium: id1=0 jd1=0 ad1=0.0d0 c dominant currently: id2=0 jd2=0 ad2=0.0d0 do 511 j=1,nocc jm=moind(j) do 511 k=1,nvirt km=moind(k) if(dabs(cij0(j+(k-1)*nocc+ii)).gt.ad1)then ad1=dabs(cij0(j+(k-1)*nocc+ii)) id1=j jd1=k endif if(dabs(cij(j+(k-1)*nocc+ii)).gt.ad2)then ad2=dabs(cij(j+(k-1)*nocc+ii)) id2=j jd2=k endif 511 s1=s1+cij(jm+(km-1)*nocc+ii)*cij0(j+(k-1)*nocc+ii) 1 -cij(jm+(km-1)*nocc+nn)*cij0(j+(k-1)*nocc+nn) c if(dabs(s1).lt.stol)then write(6,6511)s1,l,id1,jd1,ad1,id2,jd2,ad2 6511 format(' Unclear correspondence',/' =',f7.3,' i:',i5, 1 '(',i4,'->',i4,f7.3,')','(',i4,'->',i4,f7.3,')') endif c c unify phase factor with equilibrium state: if(s1.lt.0.0d0)then do 512 j=1,nocc do 512 k=1,nvirt cij(j+(k-1)*nocc+nn)=-cij(j+(k-1)*nocc+nn) 512 cij(j+(k-1)*nocc+ii)=-cij(j+(k-1)*nocc+ii) endif do 51 j=1,nocc jm=moind(j) do 51 k=1,nvirt km=moind(k) u(j+(k-1)*nocc+n1)=u(j+(k-1)*nocc+n1)+ff*cij(jm+(km-1)*nocc+ii) u(j+(k-1)*nocc+n2)=u(j+(k-1)*nocc+n2)+ff*cij(jm+(km-1)*nocc+nn) uu(j+(k-1)*nocc+n1)=uu(j+(k-1)*nocc+n1)+cij(jm+(km-1)*nocc+ii) 51 uu(j+(k-1)*nocc+n2)=uu(j+(k-1)*nocc+n2)+cij(jm+(km-1)*nocc+nn) c goto 3333 c 333333333333333333333333333333333333333333333333333333333333 c 999 close(2) C if(step.lt.1.0d-9)call report('Diff. step too small!') c i=0 do 15 ia=1,nat do 15 ix=1,3 i=i+1 do 15 l=1,n ii=(l-1)*no nn=n*no+ii n1=ii+2*n*no*(i-1) n2=nn+2*n*no*(i-1) do 15 k=1,nvirt do 15 j=1,nocc u(j+(k-1)*nocc+n1)=u(j+(k-1)*nocc+n1)/(2.0d0*step) u(j+(k-1)*nocc+n2)=u(j+(k-1)*nocc+n2)/(2.0d0*step) uu(j+(k-1)*nocc+n1)=uu(j+(k-1)*nocc+n1)/step**2 15 uu(j+(k-1)*nocc+n2)=uu(j+(k-1)*nocc+n2)/step**2 call writecij(nat,n,no,m1,nocc,nvirt,u,uu,'CIJ2.SCR.TXT',lopt) c write(6,*)'cij2-leaving' RETURN END subroutine dimensions(n,fo,nmo,nmo2,nat,lzmat,nocc,m1,m2,na,nb) c returns number of transitions and atoms implicit none integer*4 n,nat,l,I,q,nmo,nmo2,nocc,na,nb,m1,m2 character*80 s80 character*70 st character*(*) fo logical lzmat,lex real*8 axdum n=0 nat=0 nmo=0 na=0 nb=0 nocc=0 m1=0 m2=0 inquire(file='DIMS.SCR.TXT',exist=lex) if(lex)then write(6,*)' from DIMS.SCR.TXT' open(2,file='DIMS.SCR.TXT') call rdims(2,nat,nmo,nocc,n,m1,m2,na,nb) goto 99 else write(6,*)' from '//fo 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 if(nat.eq.0)then 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 ENDIF if(s80(31:49).eq.'primitive gaussians'.and.nmo.eq.0)then read(s80(1:6),*)nmo endif if(s80(8:22).eq.'alpha electrons')read(s80(1:6),*)na if(s80(33:46).eq.'beta electrons')read(s80(23:31),*)nb if(s80(5:27).eq.'ge of M.O.s used for co')read(s80(38:49),*)m1,m2 if(s80(2:44).eq.'Ground to excited state transition electric')then n=0 read(2,*) 3 read(2,3000,end=1,err=1)st 3000 format(10x,a70) read(st,*,end=1,err=1)axdum,axdum,axdum n=n+1 goto 3 endif if(n.eq.0.or.nmo.eq.0.or.nat.eq.0.or.m1.eq.0.or.m2.eq.0)goto 1 99 close(2) nmo2=nmo**2 if(na.ne.nb)call report('open shell not implemented') nocc=na return end c ============================================================ subroutine chkgeop(r,nat,lzmat,step,ID,asign,ierr,ia,ix,ra,icon, 1bohr) implicit none integer*4 nat,ig98,I,l,ic,ID,ierr,ia,ix,icon,ii real*8 r(nat,3),tol,step,asign,ra(nat,3) character*50 S50 logical lzmat real*8 bohr,x,y,z c ierr=0 ID=0 201 READ(2,1111,END=202)S50 1111 FORMAT(A50) IF((lzmat.and.(s50(19:39).EQ.'Z-Matrix orientation:'.OR. 1 s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(20:37).EQ.'Input orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:')) 1 .OR. 1 ((.not.lzmat).and. 2 (s50(20:40).EQ.'Standard orientation:'.OR. 2 s50(26:46).EQ.'Standard orientation:')))THEN ig98=0 if(s50(26:46).EQ.'Z-Matrix orientation:'.OR. 1 s50(27:44).EQ.'Input orientation:'.OR. 1 s50(26:46).EQ.'Standard orientation:')ig98=1 if(ig98.ne.-1)then DO 41 I=1,4 41 READ(2,*) endif tol=1.0d-4 ic=0 DO 1 l=1,nat if(ig98.eq.0)then READ(2,*)I,I,x,y,z else if(ig98.eq.-1)then READ(2,*)I,x,y,z else READ(2,*)x,x,x,x,y,z endif endif ra(l,1)=x/bohr ra(l,2)=y/bohr ra(l,3)=z/bohr if(icon.eq.1)then do 2 ii=1,3 if(abs(ra(l,ii)-r(l,ii)).gt.tol)then ic=ic+1 ID=3*(l-1)+ii ia=l ix=ii if(ra(l,ii).gt.r(l,ii))then asign=1.0d0 else asign=-1.0d0 endif step=dabs(ra(l,ii)-r(l,ii)) endif 2 continue endif 1 continue if(icon.eq.0)then call w36('Geometry read') else if(ic.ne.1)then write(6,*)ic call report('defective differentiation found') endif endif return ENDIF goto 201 202 call w36('Geometry not found') ierr=99 return end c ============================================================ subroutine readCIEXC(nocc,nvirt,na,nb,n,nst, 1dl,dv,qv,r,r0v,r0l,e,ev,eau,lden,cij,lopen, 1dij,lnorm,lort,lwrt,cijb) implicit none character*80 s80 character*75 st logical lden,lopen,lnorm,lort,lwrt integer*4 I,na,nb,n,nd,k,ie,ni,ic,nocc,nvirt, 1nid,oa,ob,j,nst, 1nib,amax,bmax,nort,no real*8 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,cmax,skj,sum1,sum2,forward,cijb(*) real*8,allocatable::obaj(:),obak(:) real*8 ska call vz(cij,2*nocc*nvirt*n) call vz(dij,2*nocc*nvirt*n) call vz(cijb,2*nocc*nvirt*n) ie=0 lopen=.false. tq=dsqrt(2.0d0) no=nocc*nvirt 1 read(2,2000,end=2,err=2)s80 2000 format(a80) if(s80(8:22).eq.'alpha electrons')then read(s80( 1: 6),*)na read(s80(25:31),*)nb if(lwrt)write(6,2000)s80 endif if(s80(2:44).eq.'Ground to excited state transition electric')then n=0 if(lwrt)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 if(lwrt)write(6,*)n,' dipole length transitions' endif if(s80(2:51).eq. 1'Ground to excited state transition velocity dipole')then if(lwrt)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) if(lwrt)write(6,*)n,' dipole velocity transitions' endif c if(s80(2:44).eq.'Ground to excited state transition magnetic')then if(lwrt)write(6,2000)s80 read(2,*) do 10 i=1,n read(2,3000)st 10 read(st,*)(r(i,j),j=1,3) if(lwrt)write(6,*)n,' magnetic transitions' endif if(s80(2:55).eq. 1'Ground to excited state transition velocity quadrupole')then if(lwrt)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) if(lwrt)write(6,*)n,' quadrupole velocity transitions' endif c if(s80(53:63).eq.'R(velocity)') then if(lwrt)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 if(lwrt)write(6,3000)st if(lwrt)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 if(lwrt)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 if(lwrt)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=0 nib=0 nid=0 sum1=0.0d0 sum2=0.0d0 111 read(2,2000,end=299,err=299)s80 if(s80(10:11).eq.'->')then if(s80(8:8).eq.'A')then c open shell, alpha: lopen=.true. ni=ni+1 read(s80( 1: 7),*)oa read(s80(12:14),*)ob read(s80(18:30),*)forward call setabmax(oa,ob,forward,amax,bmax,cmax) if(oa.gt.nocc)call report('a>nocc') if(ob.lt.nocc)call report('a>nocc') cij(oa+(ob-nocc-1)*nocc+(ie-1)*no)=forward else if(s80(8:8).eq.'B')then c open shell, beta: nid=nid+1 lopen=.true. read(s80( 1: 7),*)oa read(s80(12:14),*)ob read(s80(18:30),*)forward if(oa.gt.nocc)call report('a>nocc') if(ob.lt.nocc)call report('a>nocc') call setabmax(oa,ob, forward,amax,bmax,cmax) dij(oa+(ob-nocc-1)*nocc+(ie-1)*no)=forward else c closed shell: ni=ni+1 read(s80( 1: 8),*)oa read(s80(12:15),*)ob read(s80(18:30),*)forward forward=forward*tq call setabmax(oa,ob,forward,amax,bmax,cmax) sum1=sum1+forward**2 if(oa.gt.nocc)call report('a>nocc') if(ob.lt.nocc)call report('a>nocc') cij(oa+(ob-nocc-1)*nocc+(ie-1)*no)=forward c if(oa.eq.3.and.ob.eq.11.and.ie.eq.1 c 1)write(6,*)forward,forward/tq endif endif goto 111 endif if(s80(10:11).eq.'<-'.and..not.lopen)then nib=nib+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 if(oa.gt.nocc)call report('a>nocc') if(ob.lt.nocc)call report('a>nocc') cback=cback*tq sum2=sum2+cback**2 cijb(oa+(ob-nocc-1)*nocc+(ie-1)*no+n*no)=cback endif endif goto 111 endif 299 continue if(lwrt)then write(6,*)ni,' alpha coefficients' write(6,*)nid,' beta coefficients' write(6,*)nib,' 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,sum1-sum2 1234 format('fnorm, bnorm, cnorm:',3f10.4) endif c seems that fnorm-bnorm=1 if(lnorm)then c renormalize: call normcij(nocc,nvirt,ie,cij,n) if(lwrt)write(6,*)' cij were normalized' endif endif if(nd.eq.n) goto 2 endif goto 1 c 2 continue if(lort)then if(nst.gt.0.and.nst.lt.n)then nort=nst else nort=n endif c check coefs: do 3300 j=1,nort do 3300 ob=1,nvirt do 3300 oa=1,nocc if(dabs(cij(oa+(ob-1)*nocc+(j-1)*no)).gt.10.0d0)then write(6,*)nocc,nvirt,no write(6,*)'state',j,'a',oa,'b',ob, 1 cij(oa+(ob-1)*nocc+(j-1)*no) call report('deffective c') endif if(dabs(cij(oa+(ob-1)*nocc+(j-1)*no+n*no)).gt.10.0d0)then write(6,*)nocc,nvirt,no write(6,*)'state',j,'a',oa,'b',ob, 1 cij(oa+(ob-1)*nocc+(j-1)*no+n*no) call report('deffective c') endif 3300 continue c reorthonormalize: if(lwrt)write(6,*)'orthogonalizing cij for ',nort,' states' if(lopen)call report('lort for open shell not implemented') allocate(obaj(2*no),obak(2*no)) c run twice for better accuracy: do 318 ic=1,2 ska=0.0d0 do 3181 j=1,nort c transcript original cij of state j to obaj: do 324 ob=1,nvirt do 324 oa=1,nocc obaj(oa+nocc*(ob-1) )=cij(oa+(ob-1)*nocc+(j-1)*no ) 324 obaj(oa+nocc*(ob-1)+no)=cij(oa+(ob-1)*nocc+(j-1)*no+n*no) do 319 k=1,j-1 c transcript cij to obak: do 3241 ob=1,nvirt do 3241 oa=1,nocc obak(oa+nocc*(ob-1) )=cij(oa+(ob-1)*nocc+(k-1)*no ) 3241 obak(oa+nocc*(ob-1)+no)=cij(oa+(ob-1)*nocc+(k-1)*no+n*no) c : skj=0.0d0 do 322 ob=1,nvirt do 322 oa=1,nocc 322 skj=skj+obaj(oa+nocc*(ob-1) )*obak(oa+nocc*(ob-1) ) 1 -obaj(oa+nocc*(ob-1)+no)*obak(oa+nocc*(ob-1)+no) c |j'>=|j>-|k>: do 319 ob=1,nvirt do 319 oa=1,nocc obaj(oa+nocc*(ob-1) )=obaj(oa+nocc*(ob-1) ) 1 - skj*obak(oa+nocc*(ob-1) ) 319 obaj(oa+nocc*(ob-1)+no)=obaj(oa+nocc*(ob-1)+no) 1 - skj*obak(oa+nocc*(ob-1)+no) ska=ska+dabs(skj) c c norm obaj call normcij(nocc,nvirt,1,obaj,1) c transcript new obaj to new cij: do 3181 oa=1,nocc do 3181 ob=1,nvirt cij(oa+(ob-1)*nocc+(j-1)*no )=obaj(oa+nocc*(ob-1) ) 3181 cij(oa+(ob-1)*nocc+(j-1)*no+n*no)=obaj(oa+nocc*(ob-1)+no) if(lwrt)write(6,*)ic,ska/dble(nort**2)/2.0d0 318 continue if(lwrt)write(6,*)'double orthogonalization done' endif c write(6,*)'leaving',cij(3+(11-nocc-1)*nocc+(1-1)*no) return end c ============================================================ subroutine readcij(io,c,nmo,nao,e,found,lwrt) implicit none integer*4 io,nmo,nao real*8 c(*),e(*) logical found,lwrt character*50 s50 found=.false. 300 read(io,1111,end=203)S50 1111 format(A50) if(S50(6:35).eq.'Molecular Orbital Coefficients'.or. 1 S50(2:11).eq.'Alpha MOs:') then c eigenvectors listed either as a result of IOP(6/7=1) or IOP(5/33=1) found=.true. if(lwrt)write(6,1111)S50 call readorb(S50,e,c,nmo,nao,2) if(lwrt)write(6,1111)'orbitals read' return endif if(S50(2:20).eq.'Excited states from')return goto 300 203 return end c ============================================================ subroutine cp(fd,bohr) implicit none integer*4 il,NBsUse,Na,i,ia,ix,ic,ii,io1,j,k,jx,Nvirt, 1nd0,nmo2,N,M,nmo,nat,nm,noo,io parameter (nd0=21,noo=100) 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 integer*4 icon(noo) integer*4,allocatable::ty(:),oy(:) real*8,allocatable::r(:),o(:),UB(:),UR(:),P(:),PN(:),Uo(:), 1sdm(:),bp(:),P1(:),U(:),cij(:),e(:),t(:),F1(:),UBT(:), 2A(:),AN(:),PV(:) c U ... dipole + coordinates, total U-matrix c UR ... dipole + coordinate, occ-virt block c UB ... magnetic, occ-virt block c UBT... magnetic, total real*8 bohr logical lpx,lsx,lexu character*80 s80,fd write(6,*) call w36('cp') call w36('^^') inquire(file='Uou.SCR.TXT',exist=lexu) if(lexu)then call w36('Uou.SCR.TXT found on disk') return endif do 27 ic=1,noo 27 icon(ic)=0 il=0 ic=0 Na=0 NBsUse=0 Nvirt=0 ia=0 open(9,file=fd) 1 read(9,900,end=99,err=99)s80 c write(6,900)s80 900 format(a80) il=il+1 if(s80(27:44).eq.'Input orientation:')then write(6,900)s80 icon(1)=1 do 2 i=1,4 2 read(9,*) ia=0 3 read(9,900)s80 if(s80(2:3).eq.'--')goto 4 ia=ia+1 allocate(r(3*ia),ty(ia)) read(s80,*)ty(ia),ty(ia),r(1+3*(ia-1)),(r(ix+3*(ia-1)),ix=1,3) if(ia.gt.1)then do 5 i=1,ia-1 ty(i)=oy(i) do 5 ix=1,3 5 r(ix+3*(i-1))=o(ix+3*(i-1)) deallocate(o,oy) endif allocate(o(3*ia),oy(ia)) do 6 i=1,ia oy(i)=ty(i) do 6 ix=1,3 6 o(ix+3*(i-1))=r(ix+3*(i-1)) deallocate(r,ty) goto 3 4 write(6,*)ia,' atoms' allocate(r(3*ia),ty(ia)) do 11 i=1,ia ty(i)=oy(i) do 11 ix=1,3 11 r(ix+3*(i-1))=o(ix+3*(i-1))/bohr call wgeo(ia,r,ty) endif if(s80(8:22).eq.'alpha electrons')then write(6,900)s80 read(s80(1:6),*)Na icon(2)=1 endif if(s80(2:8).eq.'NBsUse=')then write(6,900)s80 icon(3)=1 read(s80(9:14),*)NBsUse if(ia.eq.0.or.NBsUse.eq.0.or.Na.eq.0) 1 call report('Incomplete reading') Nvirt=NBsUse-Na nmo=NBsUse nmo2=nmo**2 nat=ia N=3*nat write(6,*)Nvirt,Na,ia,NBsUse,nmo2,N allocate(UR(Nvirt*Na*(ia*3+3)),P1(NBsUse**2*(ia*3+3))) allocate(F1(NBsUse**2*(ia*3+3))) allocate(UB(Nvirt*Na*6)) call vz(UB,Nvirt*Na*6) c the total U matrix: allocate(U(nmo2*(N+3)),t(nmo2),UBT(nmo2*3)) call vz(U,nmo2*(N+3)) call vz(UBT,nmo2*3) endif if(s80(1:35).eq.' Molecular Orbital Coefficients'.or. 1 s80(2:11).eq.'Alpha MOs:') then write(6,900)s80 allocate(cij(NBsUse*NBsUse),e(NBsUse)) call readorb(s80,e,cij,NBsUse,NBsUse,9) icon(7)=1 endif if(s80(11:61).eq. 1'Differentiating once with respect to magnetic field')then write(6,900)s80 ic=1 icon(4)=1 endif if(s80(11:66).eq. 1'Differentiating once with respect to nuclear coordinates')then write(6,900)s80 ic=3 icon(6)=1 endif if(s80(11:61).eq. 1'Differentiating once with respect to electric field')then ic=2 write(6,900)s80 icon(5)=1 endif if(s80(2:37).eq.'CPHF results for U (alpha) for IMat=')then read(s80(38:42),*)ix write(6,*)'U ix',ix if(ic.eq.1)call rdu(ix,UB,Nvirt,Na) if(ic.eq.3)then call rdu(ix,UR,Nvirt,Na) icon(8)=1 endif endif if(s80(2:36).eq.'P1 alpha (symm , AO basis) for IC =')then if(ic.eq.3)then read(s80(37:39),*)ix write(6,*)'P ix',ix call rdp(ix,P1,NBsUse,NBsUse) icon(9)=1 endif endif if(s80(2:33).eq.'F1 (MO basis) (alpha) for IMat=')then if(ic.eq.3)then read(s80(33:37),*)ix write(6,*)'F ix',ix call rdp(ix,F1,NBsUse,NBsUse) if(ix.eq.4)icon(11)=1 endif endif goto 1 99 close(9) write(6,*)il,' lines in '//fd M=N*nmo nm=Nvirt*Na write(6,6009)'Number of atoms :',nat write(6,6009)'Number of MOs :',nmo write(6,6009)'Number of alpha electrons :',Na write(6,6009)'Number of virtual MOs :',Nvirt 6009 format(a30,i6) allocate(sdm(N*nmo2*7)) call vz(sdm,N*nmo2*7) inquire(file='SX.SCR.TXT',exist=lsx) if(lsx)then write(6,*)' Reading S-matrix derivatives ...' call readsdm(sdm,nat,nmo) write(6,*)' ... done' icon(10)=1 endif c inquire(file='S1MO.dat',exist=lsx) c if(lsx)then c allocate(sg09((3+N)*n2)) c call rs1(sg09,'S1MO.dat') c do 18 k=1,N+3 c do 18 i=1,nmo c8 sdm(k-3+N*(i-1)+M*(i-1))=sg09(j+((i-1)*i)/2+n2*(k-1))/2.0d0 c icon(10)=1 c endif call chkcon(icon) if(icon(7).eq.1.and.icon(10).eq.1.and.icon(11).eq.1)then open(44,file='UanF.SCR.TXT') write(6,*)'making u from (F-Se)/(e-e)' do 24 k=1,N+3 io=nmo2*(k-1) write(6,*)k,N+3 if(k.le.3)then do 26 i=1,nmo do 26 j=1,nmo 26 if(e(i)-e(j).ne.0.0d0)U(j+nmo*(i-1)+io)= 1 F1(i+nmo*(j-1)+io)/(e(j)-e(i)) else do 23 i=1,nmo do 23 j=1,nmo c if(i.eq.9.and.j.eq.11.and.k.eq.4)then c write(6,*)'i k j ',i,j,k c write(6,*)'F1',F1(i+nmo*(j-1)+io) c write(6,*)'S1', c 1 sdm(k-3+N*(i-1)+M*(j-1))+sdm(k-3+N*(j-1)+M*(i-1)) c write(6,*)'ei',e(i) c write(6,*)'ej',e(j) c stop c endif 23 if(e(i)-e(j).ne.0.0d0)U(j+nmo*(i-1)+io)=(F1(i+nmo*(j-1)+io)- 1 (sdm(k-3+N*(i-1)+M*(j-1))+sdm(k-3+N*(j-1)+M*(i-1)))*e(j)) 1 /(e(j)-e(i)) endif write(6,*)'done',N do 25 i=1,nmo do 25 j=1,nmo 25 t(i+nmo*(j-1))=U(i+nmo*(j-1)+nmo2*(k-1)) write(44,*)k,' - perturbation' 24 call wdu(t,nmo,nmo,44) close(44) write(6,*)' ... U done, written to UanF.SCR.TXT' endif write(6,*)'refining/making u' open(44,file='Uan.SCR.TXT') do 17 k=1,N+3 c ocupied-virtual block was read to UR, replace: do 17 i=1,Na do 17 j=Na+1,NBsUse if(k.lt.4) 1UBT(i+nmo*(j-1)+nmo2*(k-1))=UB(i+Na*(j-Na-1)+nm*(k-1)) 17 U(i+nmo*(j-1)+nmo2*(k-1))=UR(i+Na*(j-Na-1)+nm*(k-1)) write(6,*)' ... occ-virt block replaced' c c ocupied-occupied block = -(1/2) S(1), replace: do 12 k=1,N+3 if(k.gt.3)then do 13 i=1,Na do 13 j=1,Na 13 U(i+nmo*(j-1)+nmo2*(k-1))= 1 -(sdm(k-3+N*(i-1)+M*(j-1))+sdm(k-3+N*(j-1)+M*(i-1)))/2.0d0 endif do 16 i=1,nmo do 16 j=1,nmo 16 t(i+nmo*(j-1))=U(i+nmo*(j-1)+nmo2*(k-1)) write(44,*)k,' - perturbation' 12 call wdu(t,nmo,nmo,44) close(44) write(6,*)' ... occ-occ block replaced' write(6,*)' ... new U done, written in Uan.SCR.TXT' c make our U matrix - nuclear pertirbation only 1...N: write(6,*)'making u~' allocate(Uo(nmo2*N)) open(44,file='Uou.SCR.TXT') do 19 k=1,N do 20 i=1,nmo do 20 j=1,nmo 20 Uo(i+nmo*(j-1)+nmo2*(k-1))=U(i+nmo*(j-1)+nmo2*(k+3-1))+ 1sdm(k+N*(i-1)+M*(j-1)) do 21 i=1,nmo do 21 j=1,nmo 21 t(i+nmo*(j-1))=Uo(i+nmo*(j-1)+nmo2*(k-1)) write(44,*)k,' - perturbation' 19 call wdu(t,nmo,nmo,44) close(44) write(6,*)' ... u~ done, written in Uou.SCR.TXT' inquire(file='PX.MO.SCR.TXT',exist=lpx) if(.not.lpx)call report('*SCR.TXT files needed for lrds') allocate(bp(nd0*nmo2)) write(6,*)' Reading MO integrals ...' call rwm(bp,nmo,nd0) write(6,*)' ... done' c P atomic polar tensor, electronic c PV atomic polar tensor, electronic, gradient c PN atomic polar tensor, nuclear c A atomic axial tensor, electronic c AN atomic axial tensor, nuclear allocate(P(nat*9),PV(9*nat),PN(nat*9),A(nat*9),AN(nat*9)) call vz(P ,9*nat) call vz(PV,9*nat) call vz(PN,9*nat) call vz(A ,9*nat) call vz(AN,9*nat) do 7 i=1,nat AN(1+3*(2-1)+9*(i-1))=-0.25d0*ty(i)*r(3+3*(i-1)) AN(2+3*(1-1)+9*(i-1))= 0.25d0*ty(i)*r(3+3*(i-1)) AN(1+3*(3-1)+9*(i-1))= 0.25d0*ty(i)*r(2+3*(i-1)) AN(3+3*(1-1)+9*(i-1))=-0.25d0*ty(i)*r(2+3*(i-1)) AN(2+3*(3-1)+9*(i-1))=-0.25d0*ty(i)*r(1+3*(i-1)) AN(3+3*(2-1)+9*(i-1))=+0.25d0*ty(i)*r(1+3*(i-1)) PN(1+3*(1-1)+9*(i-1))=ty(i) PN(2+3*(2-1)+9*(i-1))=ty(i) 7 PN(3+3*(3-1)+9*(i-1))=ty(i) c c nucleus: do 8 ia=1,nat c nuclear index: do 8 ix=1,3 I=ix+3*(ia-1) c shift nuclear index by dipole: io1=nmo2*(I+3-1) c dipole index do 8 jx=1,3 ii=jx+3*(I-1) c sum over occupied orbitals: do 8 j=1,Na c probably wrong sdm for ix=3 and jx=3: P( ii)=P( ii)+4.0d0*sdm(I+N*(j-1)+M*(j-1)+M*nmo* jx ) PV(ii)=PV(ii)+4.0d0*sdm(I+N*(j-1)+M*(j-1)+M*nmo* jx ) c s(ix+3*(ia-1)+N*(J-1)+M*(L-1)+M*nmo*(k-1)) c ix - nuclear ix c ia - nucleus c M=N*nmo, N=3*nat c k: 1 - overlap c 2 -ux c 3 -uA c 4 -uz c sum over all orbitals: do 8 k=1,nmo P(ii)=P(ii)+4.0d0*bp(nmo2*(jx-1)+nmo*(j-1)+k) 1*U(j+nmo*(k-1)+io1) 8 A(ii)=A(ii)+UBT(j+nmo*(k-1)+nmo2*(jx-1))*U(j+nmo*(k-1)+io1) c U(is+M*(j-1)+N*M*(x-1)) c is - expanded MO N .. number of virtual orbitals c j - expanding MO M .. number of occupied orbitals c x - for uR differentiation: 1 -ux c 2 -uy c 3 -uz c 4 -R1x c ... c x -Rlastz c d/dR_i_jx d/du_ix c P(ii)=P(ii)+UR(j+Na*(k-1)+io1)*UR(j+Na*(k-1)+io2) c d/dR_i_jx d/dm_ix write(6,*)'P el' call wt(P,nat) write(6,*)'P nu' call wt(PN,nat) do 9 i=1,nat*9 9 P(i)=-P(i)+PN(i) write(6,*)'P to' call wt(P,nat) write(6,*)'A el' call wt(A,nat) write(6,*)'A nu' call wt(AN,nat) do 10 i=1,nat*9 10 A(i)=-A(i)+AN(i) write(6,*)'A to' call wt(A,nat) end subroutine w36(s) character*(*) s write(6,*)s 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 normcij(nocc,nvirt,ie,cij,n) implicit none real*8 anorm,bnorm,cij(*),f integer*4 ie,nocc,nvirt,a,b,n,no,ii,nn no=nocc*nvirt ii=(ie-1)*no nn=ii+n*no anorm=0.0d0 bnorm=0.0d0 do 316 a=1,nocc do 316 b=1,nvirt anorm=anorm+cij(a+(b-1)*nocc+ii)**2 316 bnorm=bnorm+cij(a+(b-1)*nocc+nn)**2 f=1.0d0/dsqrt(anorm-bnorm) do 317 a=1,nocc do 317 b=1,nvirt cij(a+(b-1)*nocc+ii)=cij(a+(b-1)*nocc+ii)*f 317 cij(a+(b-1)*nocc+nn)=cij(a+(b-1)*nocc+nn)*f return end subroutine mapcij(moind,cmo,cmo0,nmo,nao,lmo,emo,emo0) implicit none integer*4 moind(*),nmo,nao,i,j,ii,jbest,few,jsec logical lmo logical,allocatable::used(:) real*8 cmo(*),cmo0(*),abest,ss,asec,mtol,emo(*),emo0(*),etol allocate(used(nmo)) few=4 mtol=0.05d0 etol=0.05d0 if(lmo)then do 2 j=1,nmo 2 used(j)=.false. do 3 i=1,nmo jbest=0 jsec=0 abest=0.0d0 asec=0.0d0 c look at a few orbitals around do 4 j=max(1,i-few),min(nmo,i+few) if(.not.used(j))then if(dabs(emo0(j)-emo(i)).lt.etol)then ss=0.0d0 do 5 ii=1,nao 5 ss=ss+cmo(i+(ii-1)*nmo)*cmo0(j+(ii-1)*nmo) c write(6,*)'<',i,'|',j,'> = ',ss if(dabs(ss).gt.dabs(asec))then if(dabs(ss).lt.dabs(abest))then jsec=j asec=ss else asec=abest jsec=jbest abest=ss jbest=j endif endif endif endif 4 continue c detect possible degeneracy: if(dabs(dabs(abest)-dabs(asec)).lt.mtol.and. 1 dabs(emo(jbest)-emo(jsec)).lt.etol)then write(6,6001)i, abest,jbest,asec,jsec 6001 format(i5,' -> ',f7.3,' x',i5,f7.3,' x',i5) if(i.eq.max(jbest,jsec))then used(jbest)=.true. used(jsec)=.true. endif else if(jbest.ne.i)write(6,*)i,' -> ',jbest,abest used(jbest)=.true. endif 3 moind(i)=jbest else do 1 i=1,nmo 1 moind(i)=i endif return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N DO 100 J = 1, I 100 Z(I,J) = A(I,J) C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END c ============================================================ subroutine makeH(icon,lopt) implicit none integer*4 icon(*),ne1,ne2,nv1,nv2,ND,nm,nat,i,j,n,ie,je, 1iv,jv,IERR,no,nocc,nvirt,n3,nmo,vi,vj,oi,oj,ia,ix,ii,iii, 1jjj,k,m1,n1i,n1j,nmo2,ov,jj real*8 autocm,eau,ecm,eev,eeli,evi,enm,amu,sum1,eaur, 1sum2,sum3,d1,d2,s1,s2,s3,ami real*8,allocatable::H(:,:),C(:,:),E(:),ee(:),ew(:),cij0(:), 1Uo(:),am(:),u(:),uu(:),A(:,:),ka(:),kb1(:),kb2(:) integer*4,allocatable::q(:) integer*4 MENDELEV parameter (MENDELEV=89) logical lopt(*),lex real*4 amas(MENDELEV) data amas/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ amu=0.000548579909d0 ne1=icon(1) ne2=icon(2) nv1=icon(3) nv2=icon(4) IERR=0 c ne1 ...ne2 ... interval of electronic excitations c nv1 .. nv2 ... interval of harmonic states c ND = (ne2-ne1+2) x (nv2-nv1+2) Hamiltonian dimension c - ground state, first excited, etc. ... if(ne1.eq.0)ne1=1 if(nv1.eq.0)nv1=1 autocm=219474.0d0 call w36(' ') call w36('MakeH') call w36('^^^^^') inquire(file='ENBO.SCR',exist=lex) if(lex)then write(6,*) 'ENBO.SCR on the disk, skip this routine' return endif c c read harmonic vibrational energies from F.INP open(9,file='F.INP',status='old') read(9,*)nm,nat,nat allocate(ew(nm),q(nat),am(nat)) do 8 i=1,nat read(9,*)q(i) 8 am(i)=dble(amas(q(i)))/amu do 3 i=1,1+nat*nm+1 3 read(9,*) read(9,*)(ew(i),i=1,nm) do 4 i=1,nm 4 ew(i)=ew(i)/autocm write(6,*)'F.INP read' close(9) n3=3*nat c c read electronic excitation energies open(9,file='EE.SCR.TXT',status='old') read(9,*)n allocate(ee(n)) do 5 i=1,n 5 read(9,*)ee(i),ee(i) write(6,*)'EE.SCR.TXT read' close(9) open(9,file='CC0.SCR.TXT',status='old') c exc state J = sum (k) Cij0_J_K del_K, del_K ... Slater read(9,*)no,n,nocc,nvirt,nmo allocate(cij0(2*no*n)) read(9,902)(cij0(i),i=1,2*no*n) 902 format(4d20.14) close(9) write(6,*)'CC0.SCR.TXT read' nmo2=nmo*nmo allocate(Uo(nmo2*n3)) c a-coordinate derivative of MO_J = sum(k) U_A_J_k MO_k call readUou(n3,nmo,nocc,nmo2,Uo,icon) no=nocc*nvirt allocate(u(2*n3*no*n),uu(2*n3*no*n)) call readcij2(u,uu,m1) if(lopt(1))then write(6,*)'Projecting translations and rotations' allocate(A(n3,n3)) call DOMA(A,n3) call projectcij2(u,uu,m1,A,nat,no,nvirt,nocc,n,lopt) c this would probably be nonsense - coefficients are not c derivatives: c call projectUou(A,n3,nmo,nmo2,Uo) endif if(ne2.eq.0)ne2=n if(nv2.eq.0)nv2=nm ND=(ne2-ne1+2) * (nv2-nv1+2) write(6,601)ne2-ne1+2,nv2-nv1+2,ND 601 format(/,' Dimension of the Hamiltonian:',I6,' x',i6,' =',I9) allocate(H(ND,ND),C(ND,ND),E(ND),ka(no),kb1(no),kb2(no*n3)) i=0 do 1 ie=ne1-1,ne2 c |i> = cij del_j, del .. Slater determinant a->b if(ie.gt.ne1-1)then eeli=ee(ie) else eeli=0.0d0 endif c c pre-calculate some sums: if(ie.gt.ne1-1)then iii=(ie-1)*no do 11 vj=1,nvirt do 11 oj=m1,nocc jj=oj+(vj-1)*nocc s1=0.0d0 s2=0.0d0 do 9 ia=1,nat ami=am(ia) do 9 ix=1,3 ii=ix+3*(ia-1) n1i=iii+2*n*no*(ii-1) s3=0.0d0 do 10 vi=1,nvirt do 10 oi=m1,nocc ov=oi+(vi-1)*nocc s1=s1+cij0(ov+iii)*d2(oi,vi,oj,vj,ia,ix,ia,ix,Uo,nocc,nmo)/ami s2=s2+u(ov+n1i)*d1(oi,vi,oj,vj,ia,ix,Uo,nmo)/ami 10 s3=s3+cij0(ov+iii)*d1(oi,vi,oj,vj,ia,ix,Uo,nmo)/ami 9 kb2(jj+no*(ii-1))=s3 ka(jj)=s1 11 kb1(jj)=s2 else iii=0 endif do 1 iv=nv1-1,nv2 evi=0.0d0 if(iv.gt.nv1-1)evi=evi+ew(iv) i=i+1 j=0 do 1 je=ne1-1,ne2 if(je.gt.ne1-1)then jjj=(je-1)*no else jjj=0 endif do 1 jv=nv1-1,nv2 j=j+1 if(j.eq.i)then H(i,i)=eeli+evi write(6,*)i,H(i,i) else if(j.gt.i)then c The V0 term: sum1=0.0d0 c sum1(l,k) h^2/2Mn C_ie,l C_je,k sum2=0.0d0 c sum2(l) h^2/2Mn d C_ie,l / d R_n d C_je,l / d R_n sum3=0.0d0 c sum3(l,k) h^2/2Mn d C_ie,l / d Rn C_je,k c + h^2/2Mn C_ie,l d C_je,k / d Rn if(ie.eq.ne1-1)then if(je.eq.ne1-1)then c <0||0> do 63 ia=1,nat ami=am(ia) do 63 ix=1,3 ii=ix+3*(ia-1) 63 sum1=sum1+d2(0,0,0,0,ia,ix,ia,ix,Uo,nocc,nmo)/ami else c <0||ej> do 62 ia=1,nat ami=am(ia) do 62 ix=1,3 ii=ix+3*(ia-1) n1j=jjj+2*n*no*(ii-1) do 62 vj=1,nvirt do 62 oj=1,nocc ov=oj+(vj-1)*nocc sum3=sum3-u(ov+n1j)*d1(0,0,oj,vj,ia,ix,Uo,nmo)/ami 62 sum1=sum1+cij0(ov+jjj)*d2(0,0,oj,vj,ia,ix,ia,ix,Uo,nocc,nmo) 1 /ami endif else if(je.eq.ne1-1)then c do 61 ia=1,nat ami=am(ia) do 61 ix=1,3 ii=ix+3*(ia-1) n1i=iii+2*n*no*(ii-1) do 61 vi=1,nvirt do 61 oi=1,nocc ov=oi+(vi-1)*nocc sum3=sum3+u(ov+n1i)*d1(oi,vi,0,0,ia,ix,Uo,nmo)/ami 61 sum1=sum1+cij0(ov+iii)*d2(oi,vi,0,0,ia,ix,ia,ix,Uo,nocc,nmo) 1 /ami else c do 7 ia=1,nat ami=am(ia) do 7 ix=1,3 ii=ix+3*(ia-1) n1i=iii+2*n*no*(ii-1) n1j=jjj+2*n*no*(ii-1) do 7 k=1,nvirt do 7 jj=m1,nocc ov=jj+(k-1)*nocc 7 sum2=sum2+u(ov+n1i)*u(ov+n1j)/ami do 6 vj=1,nvirt do 6 oj=m1,nocc jj=oj+(vj-1)*nocc sum1=sum1+ka(jj) * cij0(jj+jjj) sum3=sum3+kb1(jj)* cij0(jj+jjj) do 6 ii=1,n3 6 sum3=sum3-kb2(jj+no*(ii-1))*u(jj+jjj+2*n*no*(ii-1)) endif endif c write(6,9008)i,j,sum1,sum2,sum3 9008 format(2i3,3g15.4) H(i,j)=(sum1+sum2+sum3)/2.0d0 H(j,i)=H(i,j) endif c j>i endif 1 continue write(6,*)'H made' CALL TRED12(ND,H,C,E,2,IERR) write(6,6001) 6001 format(' Energies: Relative:',/, 1' au au', 2' eV cm-1 nm') do 2 i=1,ND eau=E(i) eaur=E(i)-E(1) eev=eaur*27.211384205943d0 ecm=eev*8065.54476345045d0 if(ecm.ne.0.0d0)then enm=1.0d7/ecm else enm=0.0d0 endif 2 write(6,600)i,eau,eaur,eev,ecm,enm 600 format(i8,f13.6,f12.6,f10.2,f10.1,f9.1) c write energies and eigenvectors call wre('ENBO.SCR',E,ND) call wre2('CNBO.SCR',C,ND) write(6,*)'Energies written in ENBO.SCR' write(6,*)'Eigenvectors written in CNBO.SCR' return end c ============================================================ subroutine wre(f,a,n) implicit none integer*4 n,i character*(*) f real*8 a(*) open(9,file=f,form='unformatted') write(9)n write(9)(a(i),i=1,n) close(9) return end c ============================================================ subroutine wre2(f,a,n) implicit none integer*4 i,j,n character*(*) f real*8 a(n,n) open(9,file=f,form='unformatted') write(9)n write(9)((a(i,j),i=1,n),j=1,n) close(9) return end c ============================================================ subroutine rre(f,a,n) implicit none integer*4 n,i character*(*) f real*8 a(*) open(9,file=f,form='unformatted') read(9)n read(9)(a(i),i=1,n) close(9) return end c ============================================================ subroutine rre2(f,a,n) implicit none integer*4 i,j,n character*(*) f real*8 a(n,n) open(9,file=f,form='unformatted') read(9)n read(9)((a(i,j),i=1,n),j=1,n) close(9) return end c ============================================================ function d1(a,b,c,d,na,nx,Uo,nmo) c < a->b | d c -> d / dRm> c a->b single-excited slater determinant (ground state for a=0) c d / R_n derivative with respect to nuclear coordinate n c uo - our matrix, d MO_j/d Rn = sum(k) U(Rn)jk MO_k implicit none integer*4 a,b,c,d,na,nx,nmo,nmo2,k real*8 d1,Uo(*) nmo2=nmo**2 k=3*(na-1)+nx if(a.eq.0)then if(c.eq.0)then c <0|d0/dR> = 0 d1=0.0d0 else c <0|c->d>= d1m (h) U~_c,h U~_d,h d1=Uo(d+nmo*(c-1)+nmo2*(k-1)) endif else c a->b if(c.eq.0)then c b|0> d1=Uo(b+nmo*(a-1)+nmo2*(k-1)) else c b|c->d> if(a.eq.c)then c b|a->d> if(b.eq.d)then c b|Da->b> = 0 d1=0.0d0 else c b|Da->d> d1=Uo(d+nmo*(b-1)+nmo2*(k-1)) endif else if(b.eq.d)then c b|c->b> d1=Uo(c+nmo*(a-1)+nmo2*(k-1)) else c b|Dc->d> = 0 d1=0.0d0 endif endif endif endif return end c ============================================================ function d2(a,b,c,d,na,nx,ma,mx,Uo,nocc,nmo) c b/dRn | d c -> d / dRm> c a->b single-excited slater determinant (ground state for a=0) c d / R_n derivative with respect to nuclear coordinate n c uo - our matrix, d MO_j/d Rn = sum(k) U(Rn)jk MO_k c c but spin adapted combinations are used --- symmetric in the c spacial wavefunction c implicit none integer*4 a,b,c,d,na,nx,ma,mx,nocc,nmo,nmo2,k,l,i,j real*8 d2,Uo(*),su su=0.0d0 nmo2=nmo**2 k=3*(na-1)+nx l=3*(ma-1)+mx if(a.eq.0)then if(c.eq.0)then c <0|0> = sum (c,d) U~_c,d U~_c,d do 1 i=1,nocc do 1 j=1,nmo 1 su=su+Uo(i+nmo*(j-1)+nmo2*(k-1))*Uo(i+nmo*(j-1)+nmo2*(l-1)) su=2.0d0*su else c <0|c->d>= sum (h) U~_c,h U~_d,h do 2 j=1,nmo 2 su=su+Uo(c+nmo*(j-1)+nmo2*(k-1))*Uo(d+nmo*(j-1)+nmo2*(l-1)) endif else c a->b if(c.eq.0)then c b|0> do 3 j=1,nmo 3 su=su+Uo(a+nmo*(j-1)+nmo2*(k-1))*Uo(b+nmo*(j-1)+nmo2*(l-1)) else c b|c->d> su=Uo(a+nmo*(c-1)+nmo2*(k-1))*Uo(d+nmo*(b-1)+nmo2*(l-1)) 1 +Uo(c+nmo*(a-1)+nmo2*(k-1))*Uo(b+nmo*(d-1)+nmo2*(l-1)) if(b.eq.d)then do 4 j=1,nmo 4 su=su+Uo(a+nmo*(j-1)+nmo2*(k-1))*Uo(c+nmo*(j-1)+nmo2*(l-1)) endif if(a.eq.c)then do 5 j=1,nmo 5 su=su+Uo(b+nmo*(j-1)+nmo2*(k-1))*Uo(d+nmo*(j-1)+nmo2*(l-1)) endif endif endif d2=su return end c ============================================================ subroutine readcij2(u,uu,m1) implicit none real*8 u(*),uu(*) integer*4 nat,n,i,ia,ix,l,ii,no,nocc,nn,n1,n2,m1,j,k,nvirt open(9,file='CIJ2.SCR.TXT') read(9,*)nat,n,no,m1,nocc,nvirt i=0 do 16 ia=1,nat do 16 ix=1,3 i=i+1 do 16 l=1,n read(9,*) ii=(l-1)*no nn=n*no+ii n1=ii+2*n*no*(i-1) n2=nn+2*n*no*(i-1) read(9,*) read(9,6004)((u(j+(k-1)*nocc+n1),k=1,nvirt),j=m1,nocc) 6004 format(5e12.4) read(9,*) read(9,6004)((u(j+(k-1)*nocc+n2),k=1,nvirt),j=m1,nocc) read(9,*) read(9,6004)((uu(j+(k-1)*nocc+n1),k=1,nvirt),j=m1,nocc) read(9,*) 16 read(9,6004)((uu(j+(k-1)*nocc+n2),k=1,nvirt),j=m1,nocc) close(9) write(6,*)'CIJ2.SCR.TXT read' return end c ============================================================ subroutine projectcij2(u,uu,m1,A,nat,no,nvirt,nocc,n,lopt) implicit none real*8 u(*),uu(*),A(3*nat,3*nat),s1,s2,sd1,sd2 integer*4 nat,n,i,l,ii,no,nocc,nn,n1,n2,m1,j,k,nvirt,n3, 1jk,jj logical lopt(*) real*8,allocatable::t11(:),t12(:),td11(:),td12(:) n3=3*nat allocate(t11(n3),t12(n3),td11(n3),td12(n3)) c for each excited state: do 1 l=1,n ii=(l-1)*no nn=n*no+ii c for each coefficient: do 1 k=1,nvirt do 1 j=m1,nocc jk=j+(k-1)*nocc do 2 i=1,n3 s1=0.0d0 s2=0.0d0 sd1=0.0d0 sd2=0.0d0 do 21 jj=1,n3 c forward excitations: n1=ii+2*n*no*(jj-1)+jk c backward excitations: n2=nn+2*n*no*(jj-1)+jk s1=s1+A(i,jj)*u(n1) s2=s2+A(i,jj)*u(n2) sd1=s1+A(i,jj)**2*uu(n1) 21 sd2=s2+A(i,jj)**2*uu(n2) t11(i)=s1 t12(i)=s2 td11(i)=sd1 2 td12(i)=sd2 do 1 i=1,n3 n1=ii+2*n*no*(i-1)+jk n2=nn+2*n*no*(i-1)+jk c first derivatives: u(n1)=t11(i) u(n2)=t12(i) c diagonal second derivatives: uu(n1)=td11(i) 1 uu(n2)=td12(i) write(6,*)'u, uu projected' call writecij(nat,n,no,m1,nocc,nvirt,u,uu,'CIJ2p.SCR.TXT',lopt) return end c ============================================================ subroutine writecij(nat,n,no,m1,nocc,nvirt,u,uu,fn,lopt) implicit none integer*4 nat,no,m1,nocc,nvirt,i,ia,ix,l,n,ii,n1,n2,nn,k,j real*8 u(*),uu(*),sumnorm real*8,allocatable::sumx(:) character*(*) fn character*1 xyz(3) logical lopt(*) data xyz/'x','y','z'/ open(9,file=fn) write(9,6005)nat,n,no,m1,nocc,nvirt 6005 format(6i10) i=0 do 16 ia=1,nat do 16 ix=1,3 i=i+1 do 16 l=1,n write(9,6003)ia,xyz(ix),l 6003 format(' atom ',i5,a1,' state',i9,':') ii=(l-1)*no nn=n*no+ii n1=ii+2*n*no*(i-1) n2=nn+2*n*no*(i-1) write(9,*)'first derivatives, forward, from ',m1,'->',nocc+1,':' write(9,6004)((u(j+(k-1)*nocc+n1),k=1,nvirt),j=m1,nocc) 6004 format(5e12.4) write(9,*)'first derivatives, backward:' write(9,6004)((u(j+(k-1)*nocc+n2),k=1,nvirt),j=m1,nocc) write(9,*)'second derivatives, forward:' write(9,6004)((uu(j+(k-1)*nocc+n1),k=1,nvirt),j=m1,nocc) write(9,*)'second derivatives, backward:' 16 write(9,6004)((uu(j+(k-1)*nocc+n2),k=1,nvirt),j=m1,nocc) write(9,*)'Sum rules' allocate(sumx(no)) do 17 l=1,n ii=(l-1)*no nn=n*no+ii do 17 ix=1,3 write(9,*)'state ',l,'sum ',xyz(ix) call vz(sumx,no) do 18 ia=1,nat i=ix+(ia-1)*3 n1=ii+2*n*no*(i-1) do 18 k=1,nvirt do 18 j=m1,nocc 18 sumx(j+(k-1)*nocc)=sumx(j+(k-1)*nocc)+u(j+(k-1)*nocc+n1) sumnorm=0.0d0 do 19 k=1,nvirt do 19 j=m1,nocc 19 sumnorm=sumnorm+sumx(j+(k-1)*nocc)**2 sumnorm=dsqrt(sumnorm) if(lopt(2))write(6,6006)l,xyz(ix),sumnorm 6006 format(' State:',i6,' sum norm',a1,':',e12.4) 17 write(9,6004)((sumx(j+(k-1)*nocc),k=1,nvirt),j=m1,nocc) close(9) write(6,*)fn//' written' return end c ============================================================ subroutine readUou(n3,nmo,nocc,nmo2,Uo,icon) implicit none integer*4 k,n3,nmo,nmo2,nocc,i,j,icon(*) real*8 Uo(*),uij,uji real*8,allocatable::t(:) allocate(t(nmo2)) c read U~ matrix d MO_j/dR_n = U~(Rn)_jk MO_k: open(44,file='Uou.SCR.TXT',status='old') do 19 k=1,n3 read(44,*) call wdur(t,nmo,nmo) do 19 i=1,nmo do 19 j=1,nmo 19 Uo(i+nmo*(j-1)+nmo2*(k-1))=t(i+nmo*(j-1)) close(44) write(6,*)'U~ read from Uou.SCR.TXT ' c read U-matrix from the numerical differentiation of orbitals: if(icon(5).eq.1)then write(6,*)'Trying to read numerical U~-martix:' open(44,file='Ucp.SCR.TXT',status='old') do 17 k=1,n3 read(44,*) call wdur(t,nmo,nmo) do 17 i=nocc+1,nmo do 17 j=nocc+1,nmo 17 Uo(i+nmo*(j-1)+nmo2*(k-1))=t(i+nmo*(j-1)) close(44) write(6,*)'Virtual-virtual block read from Ucp.SCR.TXT' endif if(icon(5).eq.2)then do 21 k=1,n3 do 21 i=nocc+1,nmo do 21 j=nocc+1,nmo 21 Uo(i+nmo*(j-1)+nmo2*(k-1))=0.0d0 write(6,*)'Virtual-virtual block set to zero' endif c do 18 k=1,n3 do 18 i=1,nocc do 18 j=nocc+1,nmo 18 Uo(j+nmo*(i-1)+nmo2*(k-1))=-Uo(i+nmo*(j-1)+nmo2*(k-1)) write(6,*)'virt-occ made from occ-virt ' do 20 k=1,n3 do 20 i=1,nmo Uo(i+nmo*(i-1)+nmo2*(k-1))=0.0d0 do 20 j=1,i-1 uij=Uo(i+nmo*(j-1)+nmo2*(k-1)) uji=Uo(j+nmo*(i-1)+nmo2*(k-1)) Uo(i+nmo*(j-1)+nmo2*(k-1))=(uij-uji)/2.0d0 20 Uo(j+nmo*(i-1)+nmo2*(k-1))=(uji-uij)/2.0d0 write(6,*)'U anti-symmetrized' open(44,file='Uous.SCR.TXT') do 1 k=1,n3 do 3 i=1,nmo do 3 j=1,nmo 3 t(i+nmo*(j-1))=Uo(i+nmo*(j-1)+nmo2*(k-1)) write(44,*)k,' - perturbation' 1 call wdu(t,nmo,nmo,44) close(44) write(6,*)' ... written in Uous.SCR.TXT' return end c ============================================================ c subroutine projectUou(A,n3,nmo,nmo2,Uo) c implicit none c integer*4 k,kk,n3,nmo,nmo2,i,j,ij c real*8 Uo(*),A(n3,n3),s c real*8,allocatable::t(:) c allocate(t(n3)) c do 1 i=1,nmo c do 1 j=1,nmo c ij=i+nmo*(j-1) c do 3 k=1,n3 c s=0.0d0 c do 2 kk=1,n3 c 2 s=s+A(k,kk)*Uo(ij+nmo2*(kk-1)) c 3 t(k)=s c do 1 k=1,n3 c 1 Uo(ij+nmo2*(k-1)) = t(k) c write(6,*)'U~ projected' c c deallocate(t) c allocate(t(nmo2)) c open(44,file='Uoup.SCR.TXT') c do 19 k=1,n3 c do 21 i=1,nmo c do 21 j=1,nmo c 21 t(i+nmo*(j-1))=Uo(i+nmo*(j-1)+nmo2*(k-1)) c write(44,*)k,' - perturbation' c 19 call wdu(t,nmo,nmo,44) c close(44) c write(6,*)' ... u~ done, written in Uoup.SCR.TXT' c c return c end c ============================================================ SUBROUTINE DOMA(A,N) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N) real*8,allocatable::C(:,:),TEM(:,:) C NAT=N/3 allocate(C(3,NAT),TEM(N,N)) AMACH=1.0D-14 c read coordinates from F.INP open(4,file='F.INP',status='old') read(4,*) DO 11 I=1,NAT 11 READ(4,*)C(1,I),(C(II,I),II=1,3) CLOSE(4) DO 12 I=1,N DO 12 J=1,6 12 A(I,J)=0.0d0 DO 1 I=1,NAT A(3*I-2,1)= 1.00d0 A(3*I-1,2)= 1.00d0 A(3*I ,3)= 1.00d0 A(3*I-2,4)=-C(2,I) A(3*I-1,4)= C(1,I) A(3*I-1,5)=-C(3,I) A(3*I ,5)= C(2,I) A(3*I-2,6)= C(3,I) 1 A(3*I ,6)=-C(1,I) N6=6 CALL ORT(A,N6,NAT) ICOL=1 5 IF(ICOL.GT.N6)GOTO 8 S=0.0d0 DO 6 J=1,N 6 S=S+A(J,ICOL)**2 IF(S.LT.1000.0d0*AMACH)THEN DO 7 J=1,N 7 A(J,ICOL)=A(J,N6) N6=N6-1 GOTO 5 ENDIF ICOL=ICOL+1 GOTO 5 8 CONTINUE WRITE(*,*)N6,' coordinates projected' DO 2 I=1,N DO 2 J=1,N TEM(I,J)=0.0d0 DO 2 II=1,N6 2 TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II) DO 3 I=1,N DO 4 J=1,N 4 A(I,J)=-TEM(I,J) 3 A(I,I)=A(I,I)+1.0d0 RETURN END C ============================================================ SUBROUTINE ORT(A,N6,NAT) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(3*NAT,N6) AMACH=1.0D-14 NAT3=3*NAT C=0.0d0 DO 1 I=1,NAT3 1 C=C+A(I,1)**2 IF(C.GT.AMACH)C=1.0d0/DSQRT(C) DO 2 I=1,NAT3 2 A(I,1)=A(I,1)*C DO 3 IC=2,N6 C Column IC, orthonormalize to previous columns: DO 3 ICP=1,IC-1 SP=0.0d0 DO 5 J=1,NAT3 5 SP=SP+A(J,IC)*A(J,ICP) C Subtract projection to ICP: DO 6 J=1,NAT3 6 A(J,IC)=A(J,IC)-SP*A(J,ICP) C Normalize : SP=0.0d0 DO 7 J=1,NAT3 7 SP=SP+A(J,IC)**2 IF(SP.GT.AMACH)SP=1.0d0/DSQRT(SP) DO 3 J=1,NAT3 3 A(J,IC)=A(J,IC)*SP RETURN END C ============================================================ subroutine rgeo(n,ic,r,ty) implicit none integer*4 ty(*),n,i,ic real*8 r(*),bohr,x,y,z bohr=0.529177d0 open(9,file='FILE.X',status='old') read(9,*) read(9,*)n if(ic.gt.0)then do 1 i=1,n read(9,*)ty(i),x,y,z r(1+3*(i-1))=x/bohr r(2+3*(i-1))=y/bohr 1 r(3+3*(i-1))=z/bohr endif close(9) return end C ============================================================ subroutine wgeo(n,r,ty) implicit none integer*4 ty(*),n,i,ix real*8 r(*),bohr bohr=0.529177d0 open(19,file='FILE.X') write(19,*)'geometry from the CP part' write(19,*)n do 1 i=1,n 1 write(19,900)ty(i),(r(ix+3*(i-1))*bohr,ix=1,3) 900 format(i6,3f12.6,' 0 0 0 0 0 0 0 0.0') close(19) return end C ============================================================ 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 C ============================================================ 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 C ============================================================ 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 C ============================================================ subroutine eldipoles(nmo,JdK,no,n,cij0,nocc,dij0,cijb0,ne1,ne2) c makes dipoles between electronic BO wavefunctions, c , J,K .. TDDFT WF, ground or excited implicit none integer*4 nat,i,j,l,nd0,na,nb,nmo,nmo2,no,n,ND,nocc,id,ie,iep, 1ii,io,iv,k,nn,ne1,ne2 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 integer*4 ,allocatable::ty(:),ni(:),nid(:),nib(:), 1aij(:),bijb(:),aijb(:),dbij(:),daij(:),bij(:) real*8 ,allocatable::r(:),bp(:),bpd(:),cij(:),dij(:), 1cijb(:),pre(:),pred(:),ut(:,:) real*8 un(3),u0(nd0),u0d(nd0),JdK(*),sign,cij0(*), 1djg(nd0),djk(nd0),debye,dij0(*),cijb0(*) logical lopen,lex common/important/na,nb,ND,lopen call w36(' ') call w36('eldipoles') call w36('^^^^^^^^^') inquire(file='JdK.SCR.TXT',exist=lex) if(lex)then open(9,file='JdK.SCR.TXT') read(9,*)nn,nn,nn do 13 i=1,nd0 read(9,*) do 13 j=0,nn read(9,*) 13 read(9,804)(JdK(i+nd0*j+nd0*(nn+1)*k),k=0,nn) close(9) write(6,*)'JdK.SCR.TXT read from disk' write(6,*) return endif debye=2.541765d0 nn=ne2-ne1+1 c read MO integrals : nmo2=nmo**2 allocate(bp(nd0*nmo2),bpd(nd0*nmo2)) call rwm(bp,nmo,nd0) if(lopen)call rwmd(bpd,nmo) do 6 i=0,nn do 6 j=1,nn do 6 id=1,nd0 6 JdK(id+nd0*i+nd0*(nn+1)*j)=0.0d0 c calculate ------------------------------ call rgeo(nat,0,r,ty) allocate(r(3*nat),ty(nat)) call rgeo(nat,1,r,ty) c nuclear dipole moment: do 1 j=1,3 un(j)=0.0d0 do 1 l=1,nat 1 un(j)=un(j)+r(j+3*(l-1))*dble(ty(l)) c c electric dipole, grad,r x grad/2 and (r-rul) x grad/2 etc.: do 2 j=1,nd0 u0(j)=0.0d0 u0d(j)=0.0d0 if(lopen)then do 3 i=1,na 3 u0( j)=u0( j)+bp( i+nmo*(i-1)+nmo2*(j-1)) do 4 i=1,nb 4 u0d(j)=u0d(j)+bpd(i+nmo*(i-1)+nmo2*(j-1)) else c u0 = 2 sum(i,occ) : do 5 i=1,na 5 u0( j)=u0(j) +bp( i+nmo*(i-1)+nmo2*(j-1))*2.0d0 endif 2 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,/) do 7 j=1,nd0 7 JdK(j)=u0(j) c calculate ------------------------------ c rewrite coefficients so that old me can be used: allocate(ni(n),nid(n),nib(n),cij(nmo2*n),dij(nmo2*n), 1cijb(nmo2*n),aij(nmo2*n),bij(nmo2*n),aijb(nmo2*n),bijb(nmo2*n), 1daij(nmo2*n),dbij(nmo2*n)) do 8 ie=ne1,ne2 ni(ie)=no nid(ie)=0 nib(ie)=0 ii=0 do 9 io=1,nocc do 9 iv=nocc+1,nmo ii=ii+1 cij( ii+nmo2*(ie-1))=cij0( io+(iv-nocc-1)*nocc+(ie-1)*no) dij( ii+nmo2*(ie-1))=dij0( io+(iv-nocc-1)*nocc+(ie-1)*no) cijb(ii+nmo2*(ie-1))=cijb0(io+(iv-nocc-1)*nocc+(ie-1)*no) aij( ii+nmo2*(ie-1))=io bij( ii+nmo2*(ie-1))=iv aijb(ii+nmo2*(ie-1))=io bijb(ii+nmo2*(ie-1))=iv daij(ii+nmo2*(ie-1))=io 9 dbij(ii+nmo2*(ie-1))=iv call me(ni,nid,cij,dij,aij,daij,bij,dbij,ie,lopen,djg,bp,bpd,nmo, 1nmo2,3,nib,cijb,aijb,bijb,nd0) do 8 j=1,nd0 c : JdK(j+nd0*(ie-ne1+1))=djg(j) sign=1.0d0 if(j.ge.4.and.j.le.12)sign=-1.0d0 c : 8 JdK(j+nd0*(nn+1)*(ie-ne1+1))=djg(j)*sign c calculate ------------------------------ allocate(ut(nd0,n),pre(nd0*nmo2),pred(nd0*nmo2)) do 10 ie=ne1,ne2 c pre calculate b>: call getieab(pre,pred,lopen,na,nb,nmo,ni,nid,ie,aij,bij, 1cij,daij,dbij,dij,bp,bpd,u0,u0d,nmo2,nd0) c : call getdipole(ut,n,cij,nmo,nmo2,ni,aij,bij, 1dij,nid,daij,dbij,lopen,pre,pred,ie,nd0,nd0) do 11 j=1,nd0 11 JdK(j+nd0*(ie-ne1+1)+nd0*(nn+1)*(ie-ne1+1))=ut(j,ie) do 10 iep=ne1,ne2 if(ie.ne.iep)then c : c djk===> call dodjk(djk,nd0,ni,iep,aij,bij,nmo2,cij,pre,lopen,daij, 1 dbij,pred,nmo,dij,nid) do 12 j=1,nd0 12 JdK(j+nd0*(ie-ne1+1)+nd0*(nn+1)*(iep-ne1+1))=djk(j) endif 10 continue open(9,file='JdK.SCR.TXT') write(9,801)nd0,n,nn 801 format(3i10,' ... nd0 n nn') do 14 i=1,nd0 write(9,802)i 802 format(' perturbation ',i3) do 14 j=0,nn write(9,*)'state ',j 14 write(9,804)(JdK(i+nd0*j+nd0*(nn+1)*k),k=0,nn) 804 format(5g12.6) close(9) write(6,*)'JdK.SCR.TXT written' return end C ============================================================ 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 C ============================================================ subroutine wdims(io,nat,nmo,nocc,n,m1,m2,na,nb) implicit none integer*4 nat,nmo,nocc,n,m1,m2,io,na,nb write(io,6000)nat,nmo,nocc,na,nb,n,m1,m2 6000 format(/,' Number of atoms ',i10,/, 1 ' Number of orbitals ',i10,/, 1 ' Number of occupied orbitals',i10,/, 1 ' Number of alpha electrons ',i10,/, 1 ' Number of beta electrons ',i10,/, 1 ' Number of states ',i10,/, 1 ' Excitations start from MO ',i10,/, 1 ' Excitations end at MO ',i10,/,/) return end C ============================================================ subroutine rdims(io,nat,nmo,nocc,n,m1,m2,na,nb) implicit none integer*4 nat,nmo,nocc,n,m1,m2,io,na,nb read(io,*) read(io,6000)nat read(io,6000)nmo read(io,6000)nocc read(io,6000)na read(io,6000)nb read(io,6000)n read(io,6000)m1 read(io,6000)m2 read(io,*) 6000 format(28x,i10) return end