program comparemo c c reads coordinates, basis, orbitals and exc states c from BIG.OUT, BIG.SP c and compare these with those in files c listed in NAMES.LST c IMPLICIT none integer*4 nat0,nt0,isu0,MXB,i,noas,its,nort,j parameter (nat0=96,nt0=50,MXB=40,nort=500,isu0=nort**2) c C nat0: # of atoms, nor: # of orbitals on one atom, C ndf0: # of basis functions, c nt0:# of transitions c isu:# of expaning coefficients for exc states c integer*4 s1(nt0,isu0),s2(nt0,isu0),s1s(nt0,isu0),s2s(nt0,isu0), 1ind(nat0),ITER,it,ii,ia,ix,noa,IERR,nn,multi,multis,inds(nat0), 2ias,indmoa(nort),isn(nt0),isns(nt0),indmob(nort),isnb(nt0), 1s1b(nt0,isu0),s1bs(nt0,isu0),s2b(nt0,isu0),s2bs(nt0,isu0), 1isnsb(nt0) real*8 en(nt0),ws(nt0,isu0),ens(nt0),wss(nt0,isu0),A,XS,XB, 1RTOL,z,zsum,zs,zsums,xv,xvs,cb(3),cs(3),du(nt0,4),dus(nt0,4), 1wsb(nt0,isu0),wssb(nt0,isu0) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,nn,ITER common/atoms/z(nat0),zsum,xv(3,nat0),multi common/atomss/zs(nat0),zsums,xvs(3,nat0),multis c write(6,6767) 6767 format(' Compares electronic transitions',/, 1 ' in BIG and SMALL molecules',/,/, 2 ' Input: BIG.OUT, SMALL.OUT ... RPA outputs',/, 3 ' BIG.SP , SMALL.SP ... SP outputs',/, 4 ' COMPARE.OPT ... options',/, 5 ' CCT.INP ... atom overlap',/, 6 ' Output:COMPARE.LST ... final table',/, 7 ' DIPOLES.LST ... dipole table',/, 7 ' in au small orient',/, 7 ' DIPOLES_ROT.LST ... dipole table',/, 7 ' in au big orient',/, 7 ' COMPARE.OUT ... some controls') C call init write(6,6004) 6004 format(/,' ====== Reading BIG output ======== ',/) open(3,file='COMPARE.OUT') c get R AO MO call init write(6,*)'BIG.SP' call getRAOMO('BIG.SP',noa) c c get exc write(6,*)'BIG.OUT' call getexc('BIG.OUT',it,en,s1,s2,ws,isn,du,isu0,s1b,s2b,wsb,isnb) call writecon call normorb write(6,*)'AOs have been renormalized' c c one electron integrals for starting point call overlap6 call charge(1) call ortog write(6,*)'MOs have been orthonormalized' call charge(0) open(20,file='CCT.INP',status='old') read(20,*) read(20,*) read(20,201)(ind(i),i=1,noa) read(20,201)(ind(i),i=1,noa) 201 format(20i3) close(20) c c center do 7 ix=1,3 cb(ix)=0.0d0 nn=0 do 6 ia=1,noa if(ind(ia).ne.0)then cb(ix)=cb(ix)+xv(ix,ia) nn=nn+1 endif 6 continue 7 cb(ix)=cb(ix)/dble(nn) write(6,6003)cb c record big coordinates nn=0 do 1 ia=1,noa if(ind(ia).ne.0)then nn=nn+1 do 11 ix=1,3 11 XB(nn,ix)=xv(ix,ia)-cb(ix) endif 1 continue c write(6,*)'Big transitions' do 2 ii=1,it if(isnb(ii).eq.0.or.ws(ii,1)**2.gt.wsb(ii,1)**2)then write(6,6100)ii,en(ii),s1(ii,1),s2(ii,1),ws(ii,1) 6100 format(i2,f12.2,i4,'A',1x,i4,'A', f12.6) else write(6,6000)ii,en(ii),s1b(ii,1),s2b(ii,1),wsb(ii,1) 6000 format(i2,f12.2,i4,'B',1x,i4,'B', f12.6) endif 2 continue c C c same for small molecule: c get R AO MO write(6,6005) 6005 format(/,' ====== Reading SMALL output ======== ',/) call getRAOMOS('SMALL.SP',noas) do 8 ias=1,noas 8 inds(ias)=0 do 9 ia=1,noa 9 if(ind(ia).ne.0)inds(ind(ia))=ia write(3,*)'small atom -> big atom' do 10 ias=1,noas 10 if(inds(ias).ne.0)write(3,*)ias,inds(ias) c c get exc call getexc('SMALL.OUT',its,ens,s1s,s2s,wss,isns,dus,isu0,s1bs, 1s2bs,wssb,isnsb) call writecons call normorbs write(6,*)'AOs have been renormalized' c c one electron integrals for starting point call overlap6small call charges(1) call ortogs write(6,*)'MOs have been orthonormalized' call charges(0) c write(6,*)'Small transitions' do 4 ii=1,its if(isnsb(ii).eq.0.or.wss(ii,1)**2.gt.wssb(ii,1)**2)then write(6,6100)ii,ens(ii),s1s(ii,1),s2s(ii,1),wss(ii,1) else write(6,6000)ii,ens(ii),s1bs(ii,1),s2bs(ii,1),wssb(ii,1) endif 4 continue c c center do 15 ix=1,3 cs(ix)=0.0d0 nn=0 do 12 ia=1,noa if(ind(ia).ne.0)then nn=nn+1 cs(ix)=cs(ix)+xvs(ix,ind(ia)) endif 12 continue 15 cs(ix)=cs(ix)/dble(nn) write(6,6003)cs 6003 format('center ',3f10.2) c c record SMALL coordinates nn=0 do 3 ia=1,noa if(ind(ia).ne.0)then nn=nn+1 do 31 ix=1,3 31 XS(nn,ix)=xvs(ix,ind(ia))-cs(ix) endif 3 continue write(3,*)' coordinates BIG SMALL' do 5 ia=1,nn 5 write(3,6009)(XB(ia,ix),ix=1,3),(XS(ia,ix),ix=1,3) 6009 format(6f12.5) do 13 i=1,3 do 14 j=1,3 14 A(i,j)=0.0d0 13 A(i,i)=1.0d0 if(noa.gt.1)call DOMATRIX(IERR) write(6,6001)A 6001 format(' Transformation matrix:',/,3(3F10.4,/)) write(3,*)' coordinates BIG transformed SMALL' do 51 ia=1,nn 51 write(3,6009)(XB(ia,ix),ix=1,3), 1A(1,1)*XS(ia,1)+A(1,2)*XS(ia,2)+A(1,3)*XS(ia,3), 2A(2,1)*XS(ia,1)+A(2,2)*XS(ia,2)+A(2,3)*XS(ia,3), 3A(3,1)*XS(ia,1)+A(3,2)*XS(ia,2)+A(3,3)*XS(ia,3) call rotatesmallorbs call orbitalassignment2(indmoa,indmob,cb,cs) call transitionassignment(it,en,s1,s2,its,ens,s1s,s2s, 1isn,isns,indmoa,du,dus,isu0,s1b,s2b,s1bs,s2bs,isnb,isnsb, 1ws,wsb,wss,wssb) close(3) END subroutine rotatesmallorbs c writes control output implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*2 at character*1 ty parameter (nat0=96,nor=80,nort=500,MXB=40) common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) logical ld5,lf7,lopt,lg9 common/opts/lopt(300),icha equivalence (lopt(4 ),ld5 ),(lopt(8),lf7 ),(lopt(66),lg9) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,nn,ITER real*8 ct(3) c write(6,6000) 6000 format(' Rotating spd small orbitals',/, 1' atom type AO ') iao=1 do 1 i=1,noa do 1 j=1,nop(i) if(ish(i,j+1).gt.ish(i,j))then idel=0 if(ty(i,j).eq.'S')idel=1 if(ty(i,j).eq.'P')idel=3 if(ty(i,j).eq.'L')idel=4 if(ty(i,j).eq.'D'.and.ld5)idel=5 if(ty(i,j).eq.'D'.and..not.ld5)idel=6 if(ty(i,j).eq.'F'.and.lf7)idel=7 if(ty(i,j).eq.'F'.and..not.lf7)idel=10 if(ty(i,j).eq.'G'.and.lg9)idel=9 if(ty(i,j).eq.'G'.and..not.lg9)idel=15 if(ty(i,j).eq.'H')idel=21 if(ty(i,j).eq.'I')idel=28 if(idel.eq.0)call report('uknown shell '//ty(i,j)) if(idel.eq.1)then write(3,6001)i,ty(i,j),iao,iao+idel-1 6001 format(i4,1x,a1,i4,' -',i4,' conserved') elseif(idel.eq.3)then call rp(nmo,ct,cij,iao,A,nort) call rp(nmo,ct,bij,iao,A,nort) write(3,6002)i,ty(i,j),iao,iao+idel-1 6002 format(i4,1x,a1,i4,' -',i4,' rotated') elseif(idel.eq.4)then write(3,6001)i,'S',iao,iao call rp(nmo,ct,cij,iao+1,A,nort) call rp(nmo,ct,bij,iao+1,A,nort) write(3,6002)i,'P',iao+1,iao+idel-1 elseif(idel.eq.6)then call rd(nmo,cij,A,iao,nort) call rd(nmo,bij,A,iao,nort) write(3,6002)i,'D',iao,iao+idel-1 else call rf(idel,nmo,cij,A,iao,nort,ty(i,j)) call rf(idel,nmo,bij,A,iao,nort,ty(i,j)) write(3,6002)i,'FGH',iao,iao+idel-1 endif iao=iao+idel endif 1 continue return end c ============================================================ subroutine orbitalassignment2(indmoa,indmob,cb,cs) implicit none integer*4 nmo,nat,nao,noe,ichm,mtp,ndo,nsa,nsb,nori, 1nop,ish,nat0,nor,nort,i2,i1,io,ia,noa, 2nn,ITER,ip,npi,npii,u0, 2icha,ioo,joo,ipp,jpp,jj,MXB parameter (nat0=96,nor=80,nort=500,MXB=40,u0=28) integer*4 nxi(u0),nyi(u0),nzi(u0), 1indmoa(*),ii,indmob(*) character*1 ty character*2 at real*8 cc,alpha,rcp,cij,bij,e,be,spi,pi,bohr,debye, 1aai(u0,u0),A,XS,XB,RTOL,ai,t,tt,ef,ec,coe,con,aij, 2cci(u0),xi,yi,zi,xij,yij,zij,cb(*),cs(*), 3r2ij,b6(u0,u0),sv(u0,nort),sx,sy,sz,oi,omaxa,omaxb common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) real*8 ccs,alphas,rcps,cijs,bijs,es,bes,suma,sumb, 1aa1(u0,u0),aad(u0,u0),aaf(u0,u0) integer*4 nats,naos,nmos,noes,ichms,mtps,ndos,nsas,nsbs,noas, 1noris,nops,ishs character*1 tys character*2 ats character*3 s3 common/works/ 2ccs(nat0,nor),alphas(nat0,nor), 3rcps(nat0,nor),cijs(nort,nort),bijs(nort,nort), 4es(nort),bes(nort),nats,naos,nmos, 5noes,ichms,mtps,ndos,nsas,nsbs,noas,noris, 6nops(nat0),tys(nat0,nor),ishs(nat0,nor), 7ats(nat0) common/const/pi,bohr,debye logical ld5,lf7,lopt,lg9 common/opts/lopt(300),icha equivalence (lopt(4 ),ld5),(lopt(8),lf7),(lopt(66),lg9) integer*4 i0,i,j,jp,npj,npjj,jo,ja,nxj(u0), 1nyj(u0),nzj(u0),multis,multi,ioj,del,indo(nort),indos(nort), 1ima,imb,imas,imbs real*8 tol,xj,yj,zj,aj,ccj(u0),one,aiii,ajjj, 1z,zs,zsum,zsums,xv,xvs,aaj(u0,u0),sm(nort,nort),tem(nort,nort), 1temb(nort,nort),smb(nort,nort),wa(nort),was(nort), 1wb(nort),wbs(nort),wma,wmb,wmas,wmbs common/atoms/z(nat0),zsum,xv(3,nat0),multi common/atomss/zs(nat0),zsums,xvs(3,nat0),multis COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,nn,ITER del=0 open(77,file='COMPARE.OPT') 773 read(77,771,end=772,err=772)s3 771 format(a3) if(s3.eq.'DEL')read(77,*)del goto 773 772 close(77) if(del.ne.0)write(6,*)' Orbitals explored +/- ',del one=1.0d0 c take all orbitals in big: c make S-matrix AO call inipoly i0=0 call initfd(aa1,aad,aaf) tol=23.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,naos do 10 j=1,u0 10 sv(j,i)=0.0d0 C c =========================================BIG C atom i, orbital io io=0 do 2 ia=1,noa xi=xv(1,ia)-cb(1) yi=xv(2,ia)-cb(2) zi=xv(3,ia)-cb(3) c c primitive functions on atom ia: do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) c c =========================================SMALL C atom j, orbital jo jo=0 do 4 ja=1,noas xj=A(1,1)*(xvs(1,ja)-cs(1))+A(1,2)*(xvs(2,ja)-cs(2)) 1 +A(1,3)*(xvs(3,ja)-cs(3)) yj=A(2,1)*(xvs(1,ja)-cs(1))+A(2,2)*(xvs(2,ja)-cs(2)) 1 +A(2,3)*(xvs(3,ja)-cs(3)) zj=A(3,1)*(xvs(1,ja)-cs(1))+A(3,2)*(xvs(2,ja)-cs(2)) 1 +A(3,3)*(xvs(3,ja)-cs(3)) xij=xj-xi yij=yj-yi zij=zj-zi r2ij=xij*xij+yij*yij+zij*zij c c primitive functions on atom ja: do 5 jp=1,nops(ja) aj=alphas(ja,jp) call nppps(tys(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=0.0d0 if(ec.lt.tol)ef=spi*exp(-ec) c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi coe=con*cci(ioo) sx=oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) sy=oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) sz=oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) 71 b6(ioo,joo)=sx*sy*sz*coe c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo indos(ioj)=ja do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo indos(ioj)=ja 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif c if(ishs(ja,jp+1).gt.ishs(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, write them do 33 ii=1,npii indo(ii+io)=ia do 33 jj=1,naos 33 sm(ii+io,jj)=sv(ii,jj) c io=io+npii c c zero-out temporary buffer: do 9 ii=1,u0 do 9 jj=1,naos 9 sv(ii,jj)=0.0d0 endif 3 continue 2 continue c transform into mos c first index - from big do 901 jj=1,naos do 901 ii=1,nmo suma=0.0d0 sumb=0.0d0 do 903 ioo=1,nao suma=suma+cij(ii,ioo)*sm(ioo,jj) 903 sumb=sumb+bij(ii,ioo)*sm(ioo,jj) tem(ii,jj)=suma 901 temb(ii,jj)=sumb c c second index -from small do 902 jj=1,nmos do 902 ii=1,nmo suma=0.0d0 sumb=0.0d0 do 904 ioo=1,naos suma=suma+cijs(jj,ioo)*tem(ii,ioo) 904 sumb=sumb+bijs(jj,ioo)*temb(ii,ioo) sm(ii,jj)=suma 902 smb(ii,jj)=sumb c write(6,*)' orbital atom ' c do i=1,nao c write(6,*)i,indo(i) c enddo c write(6,*)'small:' c do i=1,naos c write(6,*)i,indos(i) c enddo write(6,6006) 6006 format(10x,'Alpha',40x,'Beta',/, 1' Big MO orbital (at) closest to (at) overlap:') do 1 i1=1,nmo omaxa=dabs(sm(i1,1)) omaxb=dabs(smb(i1,1)) indmoa(i1)=1 indmob(i1)=1 call vz(wa,noa) call vz(wb,noa) if(i1.eq.105)write(6,*)'MO',i1,'ao at cij:' do 14 i=1,nao if(i1.eq.105)write(6,*)i,indo(i),cij(i1,i) wa(indo(i))=wa(indo(i))+cij(i1,i)**2 14 wb(indo(i))=wb(indo(i))+bij(i1,i)**2 ima=1 wma=wa(1) imb=1 wmb=wb(1) do 16 ia=2,noa if(wa(ia).gt.wma)then wma=wa(ia) ima=ia endif if(wb(ia).gt.wmb)then wmb=wb(ia) imb=ia endif 16 continue do 12 i2=2,nmos if(del.eq.0.or.abs(i1-i2).lt.del)then if(omaxa.lt.dabs(sm(i1,i2)))then omaxa=dabs(sm(i1,i2)) indmoa(i1)=i2 endif if(omaxb.lt.dabs(smb(i1,i2)))then omaxb=dabs(smb(i1,i2)) indmob(i1)=i2 endif endif 12 continue c approximate distribution probabilities on atoms: call vz(was,noas) call vz(wbs,noas) if(indmoa(i1).eq.105)write(6,*)'MOs',indmoa(i1),'ao at cij:' do 18 i=1,naos if(indmoa(i1).eq.105)write(6,*)i,indos(i),cijs(indmoa(i1),i) was(indos(i))=was(indos(i))+cijs(indmoa(i1),i)**2 18 wbs(indos(i))=wbs(indos(i))+bijs(indmob(i1),i)**2 imas=1 wmas=was(1) imbs=1 wmbs=wbs(1) do 19 ia=2,noas if(was(ia).gt.wmas)then wmas=was(ia) imas=ia endif if(wbs(ia).gt.wmbs)then wmbs=wbs(ia) imbs=ia endif 19 continue if(wma .gt.1.0d0)wma =1.0d0 if(wmas.gt.1.0d0)wmas=1.0d0 if(wmb .gt.1.0d0)wmb =1.0d0 if(wmbs.gt.1.0d0)wmbs=1.0d0 1 write(6,6005)i1,ima,wma,indmoa(i1),imas,wmas,sm(i1,indmoa(i1)), 1i1,imb,wmb,indmob(i1),imbs,wmbs,smb(i1,indmob(i1)) 6005 format(i10,'(',i3,' ',f5.2,')',i10,'(',i3,' ',f5.2,')', 1f12.6,5x,i10,'(',i3,' ',f5.2,')',i10,'(',i3,' ',f5.2,')',f12.6) return end c ============================================================ subroutine transitionassignment(n,en,s1,s2,ns,ens,s1s,s2s, 1isn,isns,indmo,du,dus,isu0,s1b,s2b,s1bs,s2bs,isnb,isnsb, 1ws,wsb,wss,wssb) implicit none integer*4 nt0,isu0,j,j1,look2,ii,MXB,j1b parameter (nt0=50,MXB=40) integer*4 s1(nt0,isu0),s2(nt0,isu0),s1s(nt0,isu0),s2s(nt0,isu0), 1s1b(nt0,isu0),s2b(nt0,isu0),s1bs(nt0,isu0),s2bs(nt0,isu0), 1i,n,ns,indmo(*),n1,n2,isns(*),isn(*),nn,ITER,isnb(*),isnsb(*) real*8 en(nt0),ens(nt0),du(nt0,4),dus(nt0,4),pi,bohr,debye, 1A,XS,XB,RTOL,urx,ury,urz,sp,ur,alpha,wmax, 1ws(nt0,isu0),wsb(nt0,isu0),wss(nt0,isu0),wssb(nt0,isu0) common/const/pi,bohr,debye COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,nn,ITER character*1 spin write(6,609) 609 format(' BIG SMALL transitions:') open(66,file='COMPARE.LST') open(67,file='DIPOLES.LST') open(68,file='DIPOLES_ROT.LST') write(67,6600) 6600 format(' Ground to excited state Transition', 1' electric dipole moments (Au):',/, 2 ' state X Y', 2' Z Osc.') do 1 i=1,n c just take the biggest weight, no matter of spin: if(isnb(i).eq.0.or.ws(i,1)**2.gt.wsb(i,1)**2)then n1=s1(i,1) n2=s2(i,1) spin='A' else n1=s1b(i,1) n2=s2b(i,1) spin='B' endif c try two excitations to match the n1 -> n2 transition: j1=0 j1b=0 wmax =0.0d0 do 2 j=1,ns do 2 look2=1,2 if(look2.le.isns(j))then c take whichever spin fits best if(indmo(n1).eq.s1s(j,look2).and.indmo(n2).eq.s2s(j,look2))then if(wss(j,look2)**2.gt.wmax)then j1=j j1b=1 wmax=wss(j,look2)**2 endif endif if(indmo(n1).eq.s1bs(j,look2).and.indmo(n2).eq.s2bs(j,look2))then if(wssb(j,look2)**2.gt.wmax)then j1=j j1b=2 wmax=wssb(j,look2)**2 endif endif endif 2 continue if(j1.ne.0)then if(j1b.eq.1) 1 write(6,80)i,en(i),n1,spin,n2,spin,j1,s1s(j1,1),s2s(j1,1), 1 ens(j1),du(i,4),dus(j1,4),look2 if(j1b.eq.2) 1 write(6,80)i,en(i),n1,spin,n2,spin,j1,s1bs(j1,1),s2bs(j1,1), 1 ens(j1),du(i,4),dus(j1,4),look2 80 format(i2,f8.2,i3,a1,'->',i3,a1,i4,i4,' ->',i3,f8.2, 1 ' |u|:',2F10.5,' trial:',i2) write(67,807)i,(dus(j1,ii)/debye,ii=1,3) 807 format(i10,3x,3f12.4) urx=A(1,1)*dus(j1,1)+A(1,2)*dus(j1,2)+A(1,3)*dus(j1,3) ury=A(2,1)*dus(j1,1)+A(2,2)*dus(j1,2)+A(2,3)*dus(j1,3) urz=A(3,1)*dus(j1,1)+A(3,2)*dus(j1,2)+A(3,3)*dus(j1,3) c c take the same phase as in big: sp=urx*du(i,1)+ury*du(i,2)+urz*du(i,3) ur=urx**2+ury**2+urz**2 if(sp.lt.0)then urx=-urx ury=-ury urz=-urz sp=-sp endif if(ur*du(i,4).gt.1.0d-10)then alpha=acos(sp/dsqrt(ur*du(i,4)))*180.0d0/pi else alpha=0.0d0 endif urx=urx/debye ury=ury/debye urz=urz/debye write(68,808)i,urx,ury,urz,ur/debye**2,alpha 808 format(i10,3x,4f12.4,f8.1) write(66,806)i,en(i),n1,n2,j1,s1s(j1,1),s2s(j1,1), 1 ens(j1),du(i,4),dus(j1,4),look2,(du(i,ii),ii=1,3), 2(dus(j1,ii),ii=1,3) 806 format(i2,f8.2,i3,' ->',i3,i4,i4,' ->',i3,f8.2, 1 ' |u|:',2F10.5,' trial:',i2,6f12.6) else write(6,81)i,en(i),n1,spin,n2,spin write(66,81)i,en(i),n1,spin,n2,spin 81 format(i2,f8.2,i3,a1,'->',i3,a1,' no match') endif 1 continue close(66) close(67) close(68) return end c ============================================================== subroutine getRAOMO(ifn,n8) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*(*) ifn character*2 at character*1 ty common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/spdfgh/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) common/atoms/z(nat0),zsum,xv(3,nat0),multi logical lopt common/opts/lopt(300),icha dimension am(nat0),ipse(nat0) call fillvar(nat0,nor,nort,ifn,cc,alpha,rcp,cij,bij, 1e,be,nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,noa,nori,am, 6nop,ty,ish,at,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,nxig,nyig,nzig, 2nxih,nyih,nzih,nxii,nyii,nzii,z,zsum,xv,lopt,icha,multi,ipse) n8=nat return c ============================================================== end function ein(n,a) c int(x**n*exp(-2*a*x**2),x=-inf..inf) implicit none integer*4 n,i real*8 ein,a,b,as,pi,bohr,debye common/const/pi,bohr,debye b=1.0d0 i=1 1 b=b*dble(i) i=i+2 if(i.le.n-1)goto 1 as=dsqrt(a) ein=b*sqrt(2.0d0*pi)/(as*2.0d0)**(n+1) return end SUBROUTINE REPORT(STRING) CHARACTER*(*) STRING WRITE(6,1000) WRITE(3,1000) 1000 FORMAT(80(1H*)) WRITE(6,*)STRING WRITE(3,*)STRING WRITE(6,1000) WRITE(3,1000) WRITE(6,2000) WRITE(3,2000) 2000 FORMAT(' Program terminated') CLOSE(3) STOP END subroutine getexc(ifn,it,en,s1,s2,ws,isn,du,isu0,s1b,s2b,wsb, 1isnb) IMPLICIT none integer*4 nt0,it,isu,isu0,ii,ii1,i,isub character*(*) ifn PARAMETER (nt0=50) real*8 pi,bohr,debye common/const/pi,bohr,debye c Limits: c nt0 -# of transitions on one chromophore c MXB -# of overlaped atoms c npar0 -# of parameters to follow c n40 -# of potential sites c nwl0 -# of fitted transitions REAL*8 en(nt0),ws(nt0,isu0),t,dx,dy,dz,du(nt0,4),wsb(nt0,isu0) integer s1(nt0,isu0),s2(nt0,isu0),nt,isn(*),isnb(*) integer s1b(nt0,isu0),s2b(nt0,isu0) character*80 s80,s10 logical nu c c Input: ifn -- TDDFT G Output c c Read in transition energies and dipoles from gaussian: it=0 write(6,*)ifn open(30,file=ifn,status='old') 10 read(30,300,end=949,err=949)s80 300 format(a80) if(s80(2:44).eq.'Ground to excited state Transition electric'.or. 1 s80(2:44).eq.'Ground to excited state transition electric')then it=0 read(30,*) 31 read(30,310)s10 310 format(a10) if(nu(s10(10:10)))then it=it+1 if(it.gt.nt0)call report('too many transitions') backspace 30 read(30,*)dx,dx,dy,dz c dipole in debye: dx=dx*debye dy=dy*debye dz=dz*debye du(it,1)=dx du(it,2)=dy du(it,3)=dz du(it,4)=dx**2+dy**2+dz**2 goto 31 endif backspace 30 endif if(s80(2:14).eq.'Excited State')then do 203 i=1,79 if(s80(i:i).eq.':')then read(s80(15:i-1),*)it if(it.gt.nt0)call report('too many transitions') endif if(s80(i:i+1).eq.'eV')read(s80(i-10:i-2),*)en(it) 203 if(s80(i:i+1).eq.'nm')read(s80(i- 9:i-2),*)en(it) isu=0 isub=0 1 read(30,300)s80 if(nu(s80(8:8)).or.nu(s80(7:7)))then if(nu(s80(8:8)))then read(s80( 1: 8),*)s1(it,isu) read(s80(12:14),*)s2(it,isu) read(s80(18:len(s80)),*)ws(it,isu) isu=isu+1 if(isu.gt.isu0)call report('isu > isu0') else if(s80(8:8).eq.'A')then read(s80( 1: 7),*)s1(it,isu) read(s80(12:14),*)s2(it,isu) read(s80(18:len(s80)),*)ws(it,isu) isu=isu+1 if(isu.gt.isu0)call report('isu > isu0') else isub=isub+1 if(isub.gt.isu0)call report('isu > isu0') read(s80( 1: 7),*)s1b(it,isub) read(s80(12:14),*)s2b(it,isub) read(s80(18:len(s80)),*)wsb(it,isub) endif endif if(isu.lt.isu0)goto 1 endif c c order contributions: do 2 ii1=1,isu-1 do 2 ii=ii1+1,isu if(dabs(ws(it,ii)).gt.dabs(ws(it,ii1)))then t=ws(it,ii) ws(it,ii)=ws(it,ii1) ws(it,ii1)=t nt=s1(it,ii) s1(it,ii)=s1(it,ii1) s1(it,ii1)=nt nt=s2(it,ii) s2(it,ii)=s2(it,ii1) s2(it,ii1)=nt endif 2 continue do 3 ii1=1,isub-1 do 3 ii=ii1+1,isub if(dabs(wsb(it,ii)).gt.dabs(wsb(it,ii1)))then t=wsb(it,ii) wsb(it,ii)=wsb(it,ii1) wsb(it,ii1)=t nt=s1b(it,ii) s1b(it,ii)=s1b(it,ii1) s1b(it,ii1)=nt nt=s2b(it,ii) s2b(it,ii)=s2b(it,ii1) s2b(it,ii1)=nt endif 3 continue isn(it)=isu isnb(it)=isub endif goto 10 949 close(30) write(6,*)it,' transitions in '//ifn return end function nu(c) character*1 c logical nu nu=c.eq.'1'.or.c.eq.'2'.or.c.eq.'3'.or.c.eq.'4'.or.c.eq.'5'. 1or.c.eq.'6'.or.c.eq.'7'.or.c.eq.'8'.or.c.eq.'9'.or.c.eq.'0' return end SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=40) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=0.0000001d0 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.0.00000000001d0)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=40) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER DIMENSION ANG(3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 R=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 R=R+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=R RETURN END c subroutine ortog c re-ortonormalizes MOs implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*2 at character*1 ty common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/numbers/z0,one,two,three,four,five,half,a10 dimension s(nort,nort),oldvec(nort),oldveb(nort), 1tb(nort),tc(nort) c c read in overlap matrix for AO write(6,*)' ortogonalizing MOs' open(33,file='SAO.SCR',form='unformatted') do 1 i=1,nao read(33)(s(i,j),j=1,i) do 2 j=1,i-1 2 s(j,i)=s(i,j) 1 continue close(33) c c orthonormalize do 3 i=1,nmo c do 6 ii=1,nao oldveb(ii)=bij(i,ii) 6 oldvec(ii)=cij(i,ii) c do 41 jj=1,nao tb(jj)=0.0d0 tc(jj)=0.0d0 do 41 ii=1,nao tb(jj)=tb(jj)+oldveb(ii)*s(ii,jj) 41 tc(jj)=tc(jj)+oldvec(ii)*s(ii,jj) do 4 j=1,i-1 aij=vecpro2(j,cij,nao,tc) aib=vecpro2(j,bij,nao,tb) do 4 ii=1,nao oldveb(ii)=oldveb(ii)-aib*bij(j,ii) 4 oldvec(ii)=oldvec(ii)-aij*cij(j,ii) c an=one/vecnor(oldvec,s,nao) bn=one/vecnor(oldveb,s,nao) write(3,600)i,an,bn 600 format(i5,2f7.4,$) if(mod(i,4).eq.0)write(3,*) do 3 ii=1,nao bij(i,ii)=oldveb(ii)*bn 3 cij(i,ii)=oldvec(ii)*an write(3,*) write(6,*)' MOs have been re-ortonormalized' return end c ============================================================= function vecpro2(j,cij,nao,t) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nort=500) dimension cij(nort,nao),t(*) sum=0.0D0 do 1 jj=1,nao 1 sum=sum+cij(j,jj)*t(jj) vecpro2=sum return end c ============================================================= function vecnor(c,s,n) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nort=500) dimension c(n),s(nort,n) c sum=0.0d0 do 1 ii=1,n do 1 jj=1,n 1 sum=sum+c(ii)*c(jj)*s(ii,jj) vecnor=sqrt(sum) return end C ============ subroutine charge(ic) c calculates charge, for debugging implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*2 at character*1 ty common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/numbers/z0,one,two,three,four,five,half,a10 common/atoms/z(nat0),zsum,xv(3,nat0),multi logical lopt,lnno common/opts/lopt(300),icha equivalence (lopt(57),lnno) dimension s(nort) c qn=z0 do 2 i=1,noa 2 qn=qn+z(i) c qe=0.0d0 c open(33,file='SAO.SCR',form='unformatted') do 1 i=1,nao read(33)(s(j),j=1,i) do 11 k=1,ndo do 11 j=1,i if(j.eq.i)then ff=2.0d0 else ff=4.0d0 endif 11 qe=qe-cij(k,i)*cij(k,j)*s(j)*ff do 12 k=1,nsa do 12 j=1,i if(j.eq.i)then ff=1.0d0 else ff=2.0d0 endif 12 qe=qe-cij(k,i)*cij(k,j)*s(j)*ff do 13 k=1,nsb do 13 j=1,i if(j.eq.i)then ff=1.0d0 else ff=2.0d0 endif 13 qe=qe-bij(k,i)*bij(k,j)*s(j)*ff 1 continue close(33) dt=qn+qe write(6,4000)qe,qn,dt 4000 format(/,' Charge electronic ',F16.8,/, 1 ' nuclear ',F16.8,/, 2 ' [a.u.] ',16(1H-),/, 3 ' total ',F16.8,/) if(abs(dt-icha).gt.0.001d0.and.ic.eq.0) 1call report('Non-zero charge') return end C ============ subroutine overlap6small c the "true" d-orbitals are included here c calculates overlap matrix and other one electron integrals c symmetry properties of overlaps are not used, because of control implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*1 ty parameter (nat0=96,nor=80,nort=500) integer*4 u0,multi parameter (u0=28) common/atomss/q(nat0),zsum,xv(3,nat0),multi character*2 at common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor),at(nat0) common/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half,a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(u0,nort),b6(u0,u0),aai(u0,u0),aaj(u0,u0), 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0) logical ld5,lf7,ldisk character*3 s3 real*8 s(nort,nort) c open(77,file='COMPARE.OPT') 773 read(77,771,end=772,err=772)s3 771 format(a3) if(s3.eq.'LD5')read(77,*)ld5 if(s3.eq.'LF7')read(77,*)lf7 goto 773 772 close(77) call inipoly i0=0 call initfd(aa1,aad,aaf) c inquire(file='SAOS.SCR',exist=ldisk) if(ldisk)return tol=100.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,nao do 10 j=1,u0 10 sv(j,i)=z0 write(3,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C open(33,file='SAOS.SCR',form='unformatted') c C atom i, orbital io io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) c c primitive functions on atom ia: do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppps(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) c C atom j, orbital jo jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij c c primitive functions on atom ja: do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppps(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi 71 b6(ioo,joo)=con*cci(ioo)* 1oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj)* 2oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj)* 3oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, then write them do 33 ii=1,npii iorb=ii+io write(33)(sv(ii,jj),jj=1,iorb) c c check the normalization condition - orbital iorb anorm=one/sqrt(sv(ii,iorb)) if(abs(anorm-one).gt.0.001D0)then write(3,*)anorm,iorb endif 33 continue io=io+npii c c zero-out temporary buffer: do 9 ii=1,u0 do 9 jj=1,nao 9 sv(ii,jj)=z0 endif 3 continue 2 continue close(33) open(33,file='SAOS.SCR',form='unformatted') do 34 i=1,nao 34 read(33)(s(i,j),j=1,i) do 35 i=1,nao do 35 j=1,i-1 35 s(j,i)=s(i,j) close(33) call wmtr(s,nort,nao,'Overlap matrix, AO',i0) write(6,*)' S-matrix calculated by overlap6' return end C ============ subroutine wmtr(A,n0,n,st,ic) c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n0,n) character*(*) st write(3,*)st write(3,*) N1=1 1 N3=MIN(N1+4,N) WRITE(3,18)(I,I=N1,N3) 18 FORMAT(4X,5I14) lnst=n1 if(ic.eq.1)lnst=1 DO 130 LN=lnst,N IEND=MIN(LN,N3) if(ic.eq.1)iend=n3 130 WRITE(3,17)LN,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(3,*) return end C ============ subroutine overlap6 c the "true" d-orbitals are included here c calculates overlap matrix and other one electron integrals c symmetry properties of overlaps are not used, because of control implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*1 ty parameter (nat0=96,nor=80,nort=500) integer*4 u0,multi parameter (u0=28) c n1e ... number of kinds of one-electron integrals common/atoms/q(nat0),zsum,xv(3,nat0),multi character*2 at common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor),at(nat0) common/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half,a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(u0,nort),b6(u0,u0),aai(u0,u0),aaj(u0,u0), 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0) logical ld5,lf7,ldisk character*3 s3 real*8 s(nort,nort) c open(77,file='COMPARE.OPT') 773 read(77,771,end=772,err=772)s3 771 format(a3) if(s3.eq.'LD5')read(77,*)ld5 if(s3.eq.'LF7')read(77,*)lf7 goto 773 772 close(77) call inipoly i0=0 call initfd(aa1,aad,aaf) c inquire(file='SAO.SCR',exist=ldisk) if(ldisk)return tol=100.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,nao do 10 j=1,u0 10 sv(j,i)=z0 write(3,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C open(33,file='SAO.SCR',form='unformatted') c C atom i, orbital io io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) c c primitive functions on atom ia: do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) c C atom j, orbital jo jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij c c primitive functions on atom ja: do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppp(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi 71 b6(ioo,joo)=con*cci(ioo)* 1oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj)* 2oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj)* 3oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, then write them do 33 ii=1,npii iorb=ii+io write(33)(sv(ii,jj),jj=1,iorb) c c check the normalization condition - orbital iorb anorm=one/sqrt(sv(ii,iorb)) if(abs(anorm-one).gt.0.001D0)then write(3,*)anorm,iorb endif 33 continue io=io+npii c c zero-out temporary buffer: do 9 ii=1,u0 do 9 jj=1,nao 9 sv(ii,jj)=z0 endif 3 continue 2 continue close(33) open(33,file='SAO.SCR',form='unformatted') do 34 i=1,nao 34 read(33)(s(i,j),j=1,i) do 35 i=1,nao do 35 j=1,i-1 35 s(j,i)=s(i,j) close(33) call wmtr(s,nort,nao,'Overlap matrix, AO',i0) write(6,*)' S-matrix calculated by overlap6' return end c subroutine inipoly c c initializes the Pascal triangle up to the degree nd implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nd=24) common/poly/pp(nd,nd) one=dble(1) do 1 i=1,nd pp(i,1)=one 1 pp(i,i)=one do 2 i=3,nd do 2 j=2,i-1 2 pp(i,j)=pp(i-1,j-1)+pp(i-1,j) return end subroutine writecon c writes control output implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,bohr,debye character*2 at character*1 ty parameter (nat0=96,nor=80,nort=500) common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/atoms/z(nat0),zsum,xv(3,nat0),multi logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lg9 common/opts/lopt(300),icha equivalence (lopt(1),lsprint),(lopt(2),lpprint),(lopt(3),lqprint), 1(lopt(4 ),ld5 ),(lopt(5),levprint),(lopt(6),lvprint), 1(lopt(7 ),llprint),(lopt(8),lf7 ),(lopt(9),ld2 ), 1(lopt(10),SCALED ),(lopt(11),TTT ),(lopt(12),RAMAN ), 1(lopt(13),XYZ ),(lopt(14),ABINI ),(lopt(15),cc5 ), 1(lopt(16),lan ),(lopt(17),ldo ),(lopt(18),lvv ), 1(lopt(19),lsp ),(lopt(20),lex ),(lopt(21),lscf ), 1(lopt(22),lres ),(lopt(23),lcndo ),(lopt(24),fspec ), 1(lopt(25),lrec ),(lopt(26),lc5 ),(lopt(27),lg94 ), 1(lopt(66),lg9) c mtp=multi write(3,*) 'ROAAI control output' write(3,*) 'MOLECULAR PARAMETERS' write(3,*) '------------------------------------------' write(3,*) 'Number of basis functions : ',nao write(6,*) 'Number of basis functions : ',nao write(3,*) 'Number of electrons : ',noe write(3,*) 'Charge : ',ichm write(3,*) 'Multiplicity : ',mtp write(3,*) 'Number of doubly occupied orbitals : ',ndo write(3,*) 'Number of singly occupied a-orbitals: ',nsa write(3,*) 'Number of singly occupied b-orbitals: ',nsb write(3,*) 'Number of atoms : ',noa write(3,*) 'Sum of atomic charges : ',int(zsum) write(3,*) write(3,*) iao=1 do 1 i=1,noa write(3,4000) i,nop(i) 4000 format(I3,' atom, ',I4,' primitive functions:',/, 1' shell prim alpha expansion P-expansion', 2' AO range') do 1 j=1,nop(i) if(ish(i,j+1).gt.ish(i,j))then idel=0 if(ty(i,j).eq.'S')idel=1 if(ty(i,j).eq.'P')idel=3 if(ty(i,j).eq.'L')idel=4 if(ty(i,j).eq.'D'.and.ld5)idel=5 if(ty(i,j).eq.'D'.and..not.ld5)idel=6 if(ty(i,j).eq.'F'.and.lf7)idel=7 if(ty(i,j).eq.'F'.and..not.lf7)idel=10 if(ty(i,j).eq.'G'.and.lg9)idel=9 if(ty(i,j).eq.'G'.and..not.lg9)idel=15 if(ty(i,j).eq.'H')idel=21 if(ty(i,j).eq.'I')idel=28 if(idel.eq.0)call report('uknown shell '//ty(i,j)) if(ty(i,j).ne.'L')write(3,10001)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),iao,iao+idel-1 if(ty(i,j).eq.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),rcp(i,j),iao,iao+idel-1 iao=iao+idel else if(ty(i,j).ne.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j) if(ty(i,j).eq.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),rcp(i,j) endif 10000 format(1x,i5,a1,I5,3f15.6,I4,' - ',I4) 10001 format(1x,i5,a1,I5,2f15.6,15x,I4,' - ',I4) 1 continue if(nao.gt.nort)call report('Too many basis functions') c if(levprint)then write(3,*)'Eigenvectors' else write(3,*)'Eigenvalues' endif write(3,*)'Alpha' i=1 nori=ndo if(nmo.gt.0)nori=nmo 777 j=i+4 if(j.gt.nori)j=nori write(3,1005)(e(i1),i1=i,j) 1005 format(12x,5f13.6) if(levprint)then write(3,10051)(i1,i1=i,j) 10051 format(12x,5i13) do 3 i2=1,nao 3 write(3,10052)i2,(cij(i1,i2),i1=i,j) 10052 format(i12,5f13.6) write(3,*) endif if(j.lt.nori)then i=i+5 goto 777 endif c write(3,*)'Beta' i=1 7771 j=i+4 if(j.gt.nori)j=nori write(3,1005)(be(i1),i1=i,j) if(levprint)then write(3,10051)(i1,i1=i,j) do 31 i2=1,nao 31 write(3,10052)i2,(bij(i1,i2),i1=i,j) write(3,*) endif if(j.lt.nori)then i=i+5 goto 7771 endif c open(31,file='GEO.X') write(3,*)'Geometry (A): ' write(31,*)'Geometry (A): ' write(31,*)nat do 112 i=1,nat write(31,11021)int(z(i)),(bohr*xv(j,i),j=1,3) 11021 format(i3,3f16.6,' 0 0 0 0 0 0 0 0.0') 112 write(3,1102)at(i),i,(bohr*xv(j,i),j=1,3),z(i) 1102 format(1x,a2,i3,3f16.6,f10.5) close(31) c return end subroutine normorb implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*2 at character*1 ty parameter (nat0=96,nor=80,nort=500) integer*4 u0 parameter (u0=28) common/atoms/q(nat0),zsum,xv(3,nat0),multi common/work/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/spdfgh/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) common/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half,a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(u0,nort),b6(u0,u0),aai(u0,u0),aaj(u0,u0), 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0) logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lg9 common/opts/lopt(300),icha equivalence (lopt(1),lsprint),(lopt(2),lpprint),(lopt(3),lqprint), 1(lopt(4 ),ld5 ),(lopt(5),levprint),(lopt(6),lvprint), 1(lopt(7 ),llprint),(lopt(8),lf7 ),(lopt(9),ld2 ), 1(lopt(10),SCALED ),(lopt(11),TTT ),(lopt(12),RAMAN ), 1(lopt(13),XYZ ),(lopt(14),ABINI ),(lopt(15),cc5 ), 1(lopt(16),lan ),(lopt(17),ldo ),(lopt(18),lvv ), 1(lopt(19),lsp ),(lopt(20),lex ),(lopt(21),lscf ), 1(lopt(22),lres ),(lopt(23),lcndo ),(lopt(24),fspec ), 1(lopt(25),lrec ),(lopt(26),lc5 ),(lopt(27),lg94 ), 1(lopt(66),lg9 ) dimension anorm(nort) c call inipoly i0=0 call initfd(aa1,aad,aaf) c tol=100.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,nao anorm(i)=z0 do 10 j=1,u0 10 sv(j,i)=z0 write(3,4000)nao,noa write(6,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C c first loop, determine norms: io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) if(npi.gt.u0)call report('npi > u0') jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppp(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi coe=con*cci(ioo) sx=oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) sy=oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) sz=oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) 71 b6(ioo,joo)=sx*sy*sz*coe c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif c 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, then write them do 33 ii=1,npii iorb=ii+io 33 anorm(iorb)=sv(ii,iorb) io=io+npii do 9 ii=1,u0 do 9 jj=1,nao 9 sv(ii,jj)=z0 endif 3 continue 2 continue c write(3,*)'atomic orbital norms before normalization:' write(3,3000)(anorm(i),i=1,nao) 3000 format(6f12.6) c c second loop, normalize C atom i, orbital io io=0 fac=1.0d0/sqrt(anorm(1)) do 8 ia=1,noa do 6 ip=1,nop(ia) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) ccs(ia,ip )=ccs(ia,ip )*fac ccp(ia,ip )=ccp(ia,ip )*fac ccd(ia,ip,1)=ccd(ia,ip,1)*fac ccd(ia,ip,2)=ccd(ia,ip,2)*fac ccf(ia,ip,1)=ccf(ia,ip,1)*fac ccf(ia,ip,2)=ccf(ia,ip,2)*fac ccf(ia,ip,3)=ccf(ia,ip,3)*fac do i1=1,15 cci_G(ia,ip,i1)=cci_G(ia,ip,i1)*fac enddo do i1=1,21 cci_H(ia,ip,i1)=cci_H(ia,ip,i1)*fac enddo do i1=1,28 cci_I(ia,ip,i1)=cci_I(ia,ip,i1)*fac enddo if(ish(ia,ip+1).gt.ish(ia,ip))then io=io+npii fac=1.0d0/sqrt(anorm(io+1)) endif 6 continue 8 continue return end function oi(m,ni,nj,ai,aj,xij,t,xi,xj) c c calculates integrals of the type c c infinity c / c | ni m nj -ai*(x-xi)^2 -aj*(x-xj)^2 c i= | (x-xi) x (x-xj) e e dx c | c / c -infinity c c i = sqrt(pi) * oi * exp(-ai*aj*xij^2/(ia+aj)) c c xij=(xj-xi); ai>0; aj>0 c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/numbers/z0,one,two,three,four,five,half,a10 parameter (nd=24) common/poly/pp(nd,nd) c fac=one do 1 i=1,ni+nj+m+1 1 fac=fac*t c sum=z0 ci=aj*xij*t cj=-ai*xij*t cm=(ai*xi+aj*xj)*t do 2 l=0,ni pi=pp(ni+1,l+1) do 3 ll=1,ni-l 3 pi=pi*ci do 2 k=0,nj pj=pi*pp(nj+1,k+1) do 4 kk=1,nj-k 4 pj=pj*cj do 2 n=0,m pm=pj*pp(m+1,n+1) do 5 nn=1,m-n 5 pm=pm*cm 2 sum=sum+pm*ei(n+l+k) oi=fac*sum return end c function ei(N) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c infinity c / 2 c | N -y c i= | y e dy c | c / c -infinity c c i=sqrt(pi)*ei c common/numbers/z0,one,two,three,four,five,half,a10 c nt=N/2 if(nt+nt.ne.N)then ei=z0 return endif ho=-one zl=one do 1 i=1,nt ho=ho+two 1 zl=zl*ho*half ei=zl return end c ============================================================ subroutine nppps(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) implicit none integer*4 u0,nat0,nor,npi,npii,nxi,nyi,nzi, 1ia,ip,i1,i2,nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii parameter (u0=28,nat0=96,nor=80) character*1 shell logical ld5,lf7,lg94 real*8 cci,aai,aa1,aad,aaf,ccs,ccp,ccd,ccf,cci_G, 1cci_H,cci_I dimension cci(*),nxi(*),nyi(*),nzi(*), 1aai(u0,u0),aa1(u0,u0),aad(u0,u0),aaf(u0,u0) common/spdfghs/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) lg94=.true. c nxi(1)=0 nyi(1)=0 nzi(1)=0 do 83 i1=1,u0 do 83 i2=1,u0 83 aai(i1,i2)=aa1(i1,i2) if(shell.eq.'S')then npi=1 npii=1 cci(1)=ccs(ia,ip) return endif if(shell.eq.'P')then npi=3 npii=3 nxi(1)=1 nxi(2)=0 nyi(2)=1 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=1 cci(1)=ccp(ia,ip) cci(2)=ccp(ia,ip) cci(3)=ccp(ia,ip) return endif if(shell.eq.'L')then npi=4 npii=4 nxi(2)=1 nyi(2)=0 nzi(2)=0 nxi(3)=0 nyi(3)=1 nzi(3)=0 nxi(4)=0 nyi(4)=0 nzi(4)=1 cci(1)=ccs(ia,ip) cci(2)=ccp(ia,ip) cci(3)=ccp(ia,ip) cci(4)=ccp(ia,ip) return endif if(shell.eq.'D')then npi=6 npii=6 c order xx,yy,zz,xy,xz,yz nxi(1)=2 nxi(2)=0 nyi(2)=2 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=2 nxi(4)=1 nyi(4)=1 nzi(4)=0 nxi(5)=1 nyi(5)=0 nzi(5)=1 nxi(6)=0 nyi(6)=1 nzi(6)=1 cci(1)=ccd(ia,ip,1) cci(2)=ccd(ia,ip,1) cci(3)=ccd(ia,ip,1) cci(4)=ccd(ia,ip,2) cci(5)=ccd(ia,ip,2) cci(6)=ccd(ia,ip,2) if(ld5)then npii=5 do 78 i1=1,npii do 78 i2=1,npi 78 aai(i1,i2)=aad(i1,i2) endif return endif if(shell.eq.'F')then npi=10 npii=10 c 1 2 3 4 5 6 7 8 9 10 c order xxx,yyy,zzz,xxy,xxz,yyx,yyz,zzx,zzy,xyz (cadpac) c xxx yyy zzz xyy yxx zxx xzz yzz zyy xyz (gaussian) c c use the cadpac order except for cartesian gaussian calculation: nxi(1)=3 nyi(1)=0 nzi(1)=0 nxi(2)=0 nyi(2)=3 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=3 if(lg94.and..not.lf7)then nxi(4)=1 nyi(4)=2 nzi(4)=0 nxi(5)=2 nyi(5)=1 nzi(5)=0 nxi(6)=2 nyi(6)=0 nzi(6)=1 c nxi(4)=1 c nyi(4)=2 c nzi(4)=0 c nxi(5)=2 c nyi(5)=1 c nzi(5)=0 c nxi(6)=2 c nyi(6)=0 c nzi(6)=1 nxi(7)=1 nyi(7)=0 nzi(7)=2 nxi(8)=0 nyi(8)=1 nzi(8)=2 nxi(9)=0 nyi(9)=2 nzi(9)=1 else nxi(4)=2 nyi(4)=1 nzi(4)=0 nxi(5)=2 nyi(5)=0 nzi(5)=1 nxi(6)=1 nyi(6)=2 nzi(6)=0 nxi(7)=0 nyi(7)=2 nzi(7)=1 nxi(8)=1 nyi(8)=0 nzi(8)=2 nxi(9)=0 nyi(9)=1 nzi(9)=2 endif nxi(10)=1 nyi(10)=1 nzi(10)=1 cci(1)=ccf(ia,ip,1) cci(2)=ccf(ia,ip,1) cci(3)=ccf(ia,ip,1) cci(4)=ccf(ia,ip,2) cci(5)=ccf(ia,ip,2) cci(6)=ccf(ia,ip,2) cci(7)=ccf(ia,ip,2) cci(8)=ccf(ia,ip,2) cci(9)=ccf(ia,ip,2) cci(10)=ccf(ia,ip,3) if(lf7)then npii=7 do 79 i1=1,npii do 79 i2=1,npi 79 aai(i1,i2)=aaf(i1,i2) endif return endif if(shell.eq.'G')then npi=15 npii=15 do i1=1,npi nxi(i1)=nxig(i1) nyi(i1)=nyig(i1) nzi(i1)=nzig(i1) cci(i1) =cci_G(ia,ip,i1) enddo return endif if(shell.eq.'H')then npi=21 npii=21 do i1=1,npi nxi(i1)=nxih(i1) nyi(i1)=nyih(i1) nzi(i1)=nzih(i1) cci(i1) =cci_H(ia,ip,i1) enddo return endif if(shell.eq.'I')then npi=28 npii=28 do i1=1,npi nxi(i1)=nxii(i1) nyi(i1)=nyii(i1) nzi(i1)=nzii(i1) cci(i1) =cci_I(ia,ip,i1) enddo return endif write(3,*)'ia=',ia,' ip = ',ip,' shell = ',shell call report(' Unknown shell encountered') end c ============================================================ subroutine nppp(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) implicit none integer*4 u0,nat0,nor,npi,npii,nxi,nyi,nzi, 1ia,ip,i1,i2,nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii parameter (u0=28,nat0=96,nor=80) character*1 shell logical ld5,lf7,lg94 real*8 cci,aai,aa1,aad,aaf,ccs,ccp,ccd,ccf,cci_G, 1cci_H,cci_I dimension cci(*),nxi(*),nyi(*),nzi(*), 1aai(u0,u0),aa1(u0,u0),aad(u0,u0),aaf(u0,u0) common/spdfgh/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) c lg94=.true. nxi(1)=0 nyi(1)=0 nzi(1)=0 do 83 i1=1,u0 do 83 i2=1,u0 83 aai(i1,i2)=aa1(i1,i2) if(shell.eq.'S')then npi=1 npii=1 cci(1)=ccs(ia,ip) return endif if(shell.eq.'P')then npi=3 npii=3 nxi(1)=1 nxi(2)=0 nyi(2)=1 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=1 cci(1)=ccp(ia,ip) cci(2)=ccp(ia,ip) cci(3)=ccp(ia,ip) return endif if(shell.eq.'L')then npi=4 npii=4 nxi(2)=1 nyi(2)=0 nzi(2)=0 nxi(3)=0 nyi(3)=1 nzi(3)=0 nxi(4)=0 nyi(4)=0 nzi(4)=1 cci(1)=ccs(ia,ip) cci(2)=ccp(ia,ip) cci(3)=ccp(ia,ip) cci(4)=ccp(ia,ip) return endif if(shell.eq.'D')then npi=6 npii=6 c order xx,yy,zz,xy,xz,yz nxi(1)=2 nxi(2)=0 nyi(2)=2 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=2 nxi(4)=1 nyi(4)=1 nzi(4)=0 nxi(5)=1 nyi(5)=0 nzi(5)=1 nxi(6)=0 nyi(6)=1 nzi(6)=1 cci(1)=ccd(ia,ip,1) cci(2)=ccd(ia,ip,1) cci(3)=ccd(ia,ip,1) cci(4)=ccd(ia,ip,2) cci(5)=ccd(ia,ip,2) cci(6)=ccd(ia,ip,2) if(ld5)then npii=5 do 78 i1=1,npii do 78 i2=1,npi 78 aai(i1,i2)=aad(i1,i2) endif return endif if(shell.eq.'F')then npi=10 npii=10 c 1 2 3 4 5 6 7 8 9 10 c order xxx,yyy,zzz,xxy,xxz,yyx,yyz,zzx,zzy,xyz (cadpac) c xxx yyy zzz xyy yxx zxx xzz yzz zyy xyz (gaussian) c c use the cadpac order except for cartesian gaussian calculation: nxi(1)=3 nyi(1)=0 nzi(1)=0 nxi(2)=0 nyi(2)=3 nzi(2)=0 nxi(3)=0 nyi(3)=0 nzi(3)=3 if(lg94.and..not.lf7)then nxi(4)=1 nyi(4)=2 nzi(4)=0 nxi(5)=2 nyi(5)=1 nzi(5)=0 nxi(6)=2 nyi(6)=0 nzi(6)=1 nxi(7)=1 nyi(7)=0 nzi(7)=2 nxi(8)=0 nyi(8)=1 nzi(8)=2 nxi(9)=0 nyi(9)=2 nzi(9)=1 else nxi(4)=2 nyi(4)=1 nzi(4)=0 nxi(5)=2 nyi(5)=0 nzi(5)=1 nxi(6)=1 nyi(6)=2 nzi(6)=0 nxi(7)=0 nyi(7)=2 nzi(7)=1 nxi(8)=1 nyi(8)=0 nzi(8)=2 nxi(9)=0 nyi(9)=1 nzi(9)=2 endif nxi(10)=1 nyi(10)=1 nzi(10)=1 cci(1)=ccf(ia,ip,1) cci(2)=ccf(ia,ip,1) cci(3)=ccf(ia,ip,1) cci(4)=ccf(ia,ip,2) cci(5)=ccf(ia,ip,2) cci(6)=ccf(ia,ip,2) cci(7)=ccf(ia,ip,2) cci(8)=ccf(ia,ip,2) cci(9)=ccf(ia,ip,2) cci(10)=ccf(ia,ip,3) if(lf7)then npii=7 do 79 i1=1,npii do 79 i2=1,npi 79 aai(i1,i2)=aaf(i1,i2) endif return endif if(shell.eq.'G')then npi=15 npii=15 do i1=1,npi nxi(i1)=nxig(i1) nyi(i1)=nyig(i1) nzi(i1)=nzig(i1) cci(i1) =cci_G(ia,ip,i1) enddo return endif if(shell.eq.'H')then npi=21 npii=21 do i1=1,npi nxi(i1)=nxih(i1) nyi(i1)=nyih(i1) nzi(i1)=nzih(i1) cci(i1) =cci_H(ia,ip,i1) enddo return endif if(shell.eq.'I')then npi=28 npii=28 do i1=1,npi nxi(i1)=nxii(i1) nyi(i1)=nyii(i1) nzi(i1)=nzii(i1) cci(i1) =cci_I(ia,ip,i1) enddo return endif write(3,*)'ia=',ia,' ip = ',ip,' shell = ',shell call report(' Unknown shell encountered') end c ============================================================ subroutine initfd(aa1,aad,aaf) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) integer*4 u0 parameter (u0=28) common/numbers/z0,one,two,three,four,five,half,a10 dimension aa1(u0,u0),aad(u0,u0),aaf(u0,u0) logical lopt,lturbo common/opts/lopt(300),icha equivalence (lopt(68),lturbo) c do 77 i=1,u0 do 77 j=1,u0 aa1(i,j)=z0 aad(i,j)=z0 aaf(i,j)=z0 77 aa1(i,i)=one c c Special treatment of the "true" D and F orbitals: c c Ds are calculated in the Cadpac/cartesian order, then transformed c into the G94/spherical form c c G94 order Cadpac order Turbomole order c -------------------------------- c 3z^2-R^2 1 xx d0 = (-xx-yy+2zz)/sqrt(12) c xz 2 yy d1a = xz c yz 3 zz d1b = yz c x^2-y^2 4 xy d2a = xy c xy 5 xz d2b = (xx-yy)/2 c - 6 yz c cartesian versus spherical Ds: c 1 2 3 4 5 6 c xx yy zz xy xz yz c c 1 -d1 -d1 2*d1 0 0 0 c 2 0 0 0 0 1 0 c 3 0 0 0 0 0 1 c 4 d5 -d5 0 0 0 0 c 5 0 0 0 1 0 0 c c |normalized z^2-r^2> = d1*(two*|norm.zz>-|norm.xx>-|norm.yy>) d1=half st=sqrt(three) sf=sqrt(five) c |normalized x^2-y^2> = d5*(|norm.xx>-|norm.yy>) d5=st*half aad(1,1)=-d1 aad(1,2)=-d1 aad(1,3)=two*d1 aad(2,5)=one aad(3,6)=one if(lturbo)then aad(5,1)=d5 aad(5,2)=-d5 aad(4,4)=one else aad(4,1)=d5 aad(4,2)=-d5 aad(5,4)=one endif c C F-orbitals c spherical cartesian c -------------------------------- c 5z^3-3R^2z 1 xxx c x(5z^2-R^2) 2 yyy c y(5z^2-R^2) 3 zzz c xyz 4 xxy c z(x^2-y^2) 5 xxz c y(y^2-3x^2) 6 yyx c x(x^2-3y^2) 7 yyz c - 8 zzx c - 9 zzy c - 10 xyx c cartesian versus spherical Fs: c non-normalized functions in the table: c 1 2 3 4 5 6 7 8 9 10 c xxx yyy zzz xxy xxz yyx yyz zzx zzy xyz c1 0 0 2 0 -3 0 -3 0 0 0 c2 -1 0 0 0 0 -1 0 4 0 0 c3 0 -1 0 -1 0 0 0 0 4 0 c4 0 0 0 0 1 0 -1 0 0 0 c5 0 0 0 0 0 0 0 0 0 1 c6 1 0 0 0 0 -3 0 0 0 0 c7 0 -1 0 3 0 0 0 0 0 0 c normalization factors for spherical Fs f1=half/sf f2=half*st/sqrt(a10) f3=half*st/sqrt(a10) f4=st*half f5=one f6=half/sqrt(two) f7=half/sqrt(two) c expansion into normalized cartesians: aaf(1, 3)= f1 * two * sf aaf(1, 5)=-f1 * three aaf(1, 7)=-f1 * three aaf(2, 1)=-f2 * sf aaf(2, 6)=-f2 aaf(2, 8)= f2 * four aaf(3, 2)=-f3 * sf aaf(3, 4)=-f3 aaf(3, 9)= f3 * four aaf(4, 5)= f4 aaf(4, 7)=-f4 aaf(5,10)= f5 aaf(6, 1)= f6 * sf aaf(6, 6)=-f6 * three aaf(7, 2)=-f7 * sf aaf(7, 4)= f7 * three return end c ============================================================ subroutine init implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half,a10 pi=dble(4)*atan(dble(1)) bohr=0.52917705993d0 debye=2.541747758d0 c z0=dble(0) one=dble(1) two=dble(2) three=dble(3) four=dble(4) five=dble(5) half=one/two a10=dble(10) return end c SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS subroutine getRAOMOS(ifn,n8) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*(*) ifn character*2 at character*1 ty common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/spdfghs/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) common/atomss/z(nat0),zsum,xv(3,nat0),multi logical lopt common/opts/lopt(300),icha common/mass/am(nat0),ipse(nat0) call fillvar(nat0,nor,nort,ifn,cc,alpha,rcp,cij,bij, 1e,be,nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,noa,nori,am, 6nop,ty,ish,at,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,nxig,nyig,nzig, 2nxih,nyih,nzih,nxii,nyii,nzii,z,zsum,xv,lopt,icha,multi,ipse) n8=nat return end c ============================================================== subroutine getbasisfck(nao,nat,nop,nor,ish,ld5,lf7, 1 ty,cc,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,alpha,nat0, 1 nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii,rcp, 3 nst,ist,ipps,istam,rpe,rcc,rcl) implicit none integer*4 nao,nat,i,nop(*),k,nor,nat0,ish(nat0,nor), 1ishell,ip,ndao,L,nb,ix,iy,iz,nxig(15),nyig(15),isa, 1nzig(15),nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) real*8 af,ccs(nat0,nor),ccp(nat0,nor),cc(nat0,nor), 1ccd(nat0,nor,2),ccf(nat0,nor,3),cci_G(nat0,nor,15), 1cci_H(nat0,nor,21),cci_I(nat0,nor,28),alpha(nat0,nor), 1a,SNO,pi,bohr,debye,SNP,SND,SNF1,SNF2,SNF3,SNDD, 1rcp(nat0,nor),ein character*2 s2 character*1 ty(nat0,nor) logical ld5,lf7 common/const/pi,bohr,debye integer*4 nst,ist(*),ipps(*),istam(*),ipt real*8 rpe(*),rcc(*),rcl(*) nao=0 do 1 i=1,nat nop(i)=0 do 1 k=1,nor 1 ish(i,k)=32000 c c number of shells on one atom: isa=0 c total number of primitive functions: ipt=0 c loop over shells: write(6,6019) 6019 format('shell primitive primitive atom exp:') do 2 ishell=1,nst isa=isa+1 s2='XX' if(ist(ishell).eq.-1)s2='SP' if(ist(ishell).eq. 0)s2='S ' if(ist(ishell).eq. 1)s2='P ' if(ist(ishell).eq.-2)s2='D ' if(ist(ishell).eq. 2)s2='D ' if(ist(ishell).eq. 3)s2='F ' if(ist(ishell).eq.-3)s2='F ' if(ist(ishell).eq. 4)s2='G ' if(ist(ishell).eq.-4)s2='G ' if(ist(ishell).eq. 5)s2='H ' if(ist(ishell).eq. 6)s2='I ' if(s2.eq.'XX')then write(6,*)ist(ishell) call report('unknown shell type') endif if(s2.eq.'SP')ndao=4 if(s2.eq.'P ')ndao=3 if(s2.eq.'S ')ndao=1 if(s2.eq.'D '.and.ld5)ndao=5 if(s2.eq.'D '.and..not.ld5)ndao=6 if(s2.eq.'F '.and..not.lf7)ndao=10 if(s2.eq.'F '.and.lf7)ndao=7 if(s2.eq.'G ')ndao=15 if(s2.eq.'H ')ndao=21 if(s2.eq.'I ')ndao=28 nao=nao+ndao c shell atom: i=istam(ishell) if(ishell.gt.1.and.i.ne.istam(ishell-1))isa=1 if(i.lt.0.and.i.gt.nat)call report('wrong shell atom assignment') c loop over primitives do 2 ip=1,ipps(ishell) ipt=ipt+1 nop(i)=nop(i)+1 k=nop(i) if(k.gt.nor)call report('Too many primitives on one atom') if(s2.eq.'S ')ty(i,k)='S' if(s2.eq.'SP')ty(i,k)='L' if(s2.eq.'D ')ty(i,k)='D' if(s2.eq.'P ')ty(i,k)='P' if(s2.eq.'F ')ty(i,k)='F' if(s2.eq.'G ')ty(i,k)='G' if(s2.eq.'H ')ty(i,k)='H' if(s2.eq.'I ')ty(i,k)='I' c primitive exponent: af=rpe(ipt) alpha(i,k)=af write(6,6009)ishell,ip,ipt,i,af 6009 format(4i4,f12.6) c contraction coefficients: if(s2.eq.'SP')then ccs(i,k)=rcc(ipt) ccp(i,k)=rcl(ipt) else cc(i,k)=rcc(ipt) endif c normalization factors: a=dble(2)*af/pi SNO=sqrt(sqrt(a*a*a)) a=af a=dble(128)*a*a*a*a*a/(pi*pi*pi) SNP=sqrt(sqrt(a)) a=af a=dble(2048)*a*a*a*a*a*a*a/(pi*pi*pi) SND=sqrt(sqrt(a)) SNDD=SND/sqrt(dble(3)) ccd(i,k,1)=cc(i,k)*SNDD ccd(i,k,2)=cc(i,k)*SND a=dble(2)*af a=sqrt( dble(8)*sqrt(a*a*a*a*a*a*a*a*a)/sqrt(pi*pi*pi) ) SNF1=a/sqrt(dble(15)) SNF2=a/sqrt(dble(3)) SNF3=a ccf(i,k,1)=cc(i,k)*SNF1 ccf(i,k,2)=cc(i,k)*SNF2 ccf(i,k,3)=cc(i,k)*SNF3 a=af L=4 nb=0 do 5 ix=0,L do 5 iy=0,L-ix do 5 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxig(nb)=ix nyig(nb)=iy nzig(nb)=iz cci_G(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 5 continue L=5 nb=0 do 4 ix=0,L do 4 iy=0,L-ix do 4 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxih(nb)=ix nyih(nb)=iy nzih(nb)=iz cci_H(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 4 continue L=6 nb=0 do 6 ix=0,L do 6 iy=0,L-ix do 6 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxii(nb)=ix nyii(nb)=iy nzii(nb)=iz cci_I(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 6 continue if(s2.eq.'SP')then c store exp coef. for control listing: cc(i,k)=ccs(i,k) rcp(i,k)=ccp(i,k) ccs(i,k)=ccs(i,k)*SNO ccp(i,k)=ccp(i,k)*SNP else ccs(i,k)=cc(i,k)*SNO ccp(i,k)=cc(i,k)*SNP endif if(ty(i,k).eq.'D')cc(i,k)=ccd(i,k,1)/SNDD if(ty(i,k).eq.'P')cc(i,k)=ccp(i,k)/SNP 2 ish(i,k)=isa c write(6,*) 'basis loaded from checkpoint' write(6,*) nao,' AOs' return end c ============================================================== subroutine ortogs c re-ortonormalizes MOs implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*2 at character*1 ty common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/numbers/z0,one,two,three,four,five,half,a10 dimension s(nort,nort),oldvec(nort),oldveb(nort), 1tb(nort),tc(nort) c c c read in overlap matrix for AO open(33,file='SAOS.SCR',form='unformatted') do 1 i=1,nao read(33)(s(i,j),j=1,i) do 2 j=1,i-1 2 s(j,i)=s(i,j) 1 continue close(33) c c orthonormalize do 3 i=1,nmo c do 6 ii=1,nao oldveb(ii)=bij(i,ii) 6 oldvec(ii)=cij(i,ii) c do 41 jj=1,nao tb(jj)=0.0d0 tc(jj)=0.0d0 do 41 ii=1,nao tb(jj)=tb(jj)+oldveb(ii)*s(ii,jj) 41 tc(jj)=tc(jj)+oldvec(ii)*s(ii,jj) do 4 j=1,i-1 aij=vecpro2(j,cij,nao,tc) aib=vecpro2(j,bij,nao,tb) do 4 ii=1,nao oldveb(ii)=oldveb(ii)-aib*bij(j,ii) 4 oldvec(ii)=oldvec(ii)-aij*cij(j,ii) c an=one/vecnor(oldvec,s,nao) bn=one/vecnor(oldveb,s,nao) write(3,600)i,an,bn 600 format(i5,2f7.4,$) if(mod(i,4).eq.0)write(3,*) do 3 ii=1,nao bij(i,ii)=oldveb(ii)*bn 3 cij(i,ii)=oldvec(ii)*an write(3,*) write(3,*)' MOs have been re-ortonormalized' return end c ============================================================== subroutine charges(ic) c calculates charge, for debugging implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=500) character*2 at character*1 ty common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/numbers/z0,one,two,three,four,five,half,a10 common/atomss/z(nat0),zsum,xv(3,nat0),multis logical lopt,lnno common/opts/lopt(300),icha equivalence (lopt(57),lnno) dimension s(nort) c qn=z0 do 2 i=1,noa 2 qn=qn+z(i) c qe=0.0d0 c open(33,file='SAOS.SCR',form='unformatted') do 1 i=1,nao read(33)(s(j),j=1,i) do 11 k=1,ndo do 11 j=1,i if(j.eq.i)then ff=2.0d0 else ff=4.0d0 endif 11 qe=qe-cij(k,i)*cij(k,j)*s(j)*ff do 12 k=1,nsa do 12 j=1,i if(j.eq.i)then ff=1.0d0 else ff=2.0d0 endif 12 qe=qe-cij(k,i)*cij(k,j)*s(j)*ff do 13 k=1,nsb do 13 j=1,i if(j.eq.i)then ff=1.0d0 else ff=2.0d0 endif 13 qe=qe-bij(k,i)*bij(k,j)*s(j)*ff 1 continue close(33) dt=qn+qe write(6,4000)qe,qn,dt 4000 format(/,' Charge electronic ',F16.8,/, 1 ' nuclear ',F16.8,/, 2 ' [a.u.] ',16(1H-),/, 3 ' total ',F16.8,/) if(abs(dt-icha).gt.0.001d0.and.ic.eq.0) 1call report('Non-zero charge') return end C ============ c c subroutine writecons c writes control output implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,bohr,debye character*2 at character*1 ty parameter (nat0=96,nor=80,nort=500) common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/atomss/z(nat0),zsum,xv(3,nat0),multis logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lg9 common/opts/lopt(300),icha equivalence (lopt(1),lsprint),(lopt(2),lpprint),(lopt(3),lqprint), 1(lopt(4 ),ld5 ),(lopt(5),levprint),(lopt(6),lvprint), 1(lopt(7 ),llprint),(lopt(8),lf7 ),(lopt(9),ld2 ), 1(lopt(10),SCALED ),(lopt(11),TTT ),(lopt(12),RAMAN ), 1(lopt(13),XYZ ),(lopt(14),ABINI ),(lopt(15),cc5 ), 1(lopt(16),lan ),(lopt(17),ldo ),(lopt(18),lvv ), 1(lopt(19),lsp ),(lopt(20),lex ),(lopt(21),lscf ), 1(lopt(22),lres ),(lopt(23),lcndo ),(lopt(24),fspec ), 1(lopt(25),lrec ),(lopt(26),lc5 ),(lopt(27),lg94 ), 1(lopt(66),lg9) c mtp=multis write(3,*) 'ROAAI control output' write(3,*) 'MOLECULAR PARAMETERS' write(3,*) '------------------------------------------' write(3,*) 'Number of basis functions : ',nao write(6,*) 'Number of basis functions : ',nao write(3,*) 'Number of electrons : ',noe write(3,*) 'Charge : ',ichm write(3,*) 'Multiplicity : ',mtp write(3,*) 'Number of doubly occupied orbitals : ',ndo write(3,*) 'Number of singly occupied a-orbitals: ',nsa write(3,*) 'Number of singly occupied b-orbitals: ',nsb write(3,*) 'Number of atoms : ',noa write(3,*) 'Sum of atomic charges : ',int(zsum) write(3,*) write(3,*) iao=1 do 1 i=1,noa write(3,4000) i,nop(i) 4000 format(I3,' atom, ',I4,' primitive functions:',/, 1' shell prim alpha expansion P-expansion', 2' AO range') do 1 j=1,nop(i) if(ish(i,j+1).gt.ish(i,j))then idel=0 if(ty(i,j).eq.'S')idel=1 if(ty(i,j).eq.'P')idel=3 if(ty(i,j).eq.'L')idel=4 if(ty(i,j).eq.'D'.and.ld5)idel=5 if(ty(i,j).eq.'D'.and..not.ld5)idel=6 if(ty(i,j).eq.'F'.and.lf7)idel=7 if(ty(i,j).eq.'F'.and..not.lf7)idel=10 if(ty(i,j).eq.'G'.and.lg9)idel=9 if(ty(i,j).eq.'G'.and..not.lg9)idel=15 if(ty(i,j).eq.'H')idel=21 if(ty(i,j).eq.'I')idel=28 if(idel.eq.0)call report('uknown shell '//ty(i,j)) if(ty(i,j).ne.'L')write(3,10001)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),iao,iao+idel-1 if(ty(i,j).eq.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),rcp(i,j),iao,iao+idel-1 iao=iao+idel else if(ty(i,j).ne.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j) if(ty(i,j).eq.'L')write(3,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j),rcp(i,j) endif 10000 format(1x,i5,a1,I5,3f15.6,I4,' - ',I4) 10001 format(1x,i5,a1,I5,2f15.6,15x,I4,' - ',I4) 1 continue if(nao.gt.nort)call report('Too many basis functions') c if(levprint)then write(3,*)'Eigenvectors' else write(3,*)'Eigenvalues' endif write(3,*)'Alpha' i=1 nori=ndo if(nmo.gt.0)nori=nmo 777 j=i+4 if(j.gt.nori)j=nori write(3,1005)(e(i1),i1=i,j) 1005 format(12x,5f13.6) if(levprint)then write(3,10051)(i1,i1=i,j) 10051 format(12x,5i13) do 3 i2=1,nao 3 write(3,10052)i2,(cij(i1,i2),i1=i,j) 10052 format(i12,5f13.6) write(3,*) endif if(j.lt.nori)then i=i+5 goto 777 endif c write(3,*)'Beta' i=1 7771 j=i+4 if(j.gt.nori)j=nori write(3,1005)(be(i1),i1=i,j) if(levprint)then write(3,10051)(i1,i1=i,j) do 31 i2=1,nao 31 write(3,10052)i2,(bij(i1,i2),i1=i,j) write(3,*) endif if(j.lt.nori)then i=i+5 goto 7771 endif c open(31,file='GEOS.X') write(3,*)'Geometry (A): ' write(31,*)'Geometry (A): ' write(31,*)nat do 112 i=1,nat write(31,11021)int(z(i)),(bohr*xv(j,i),j=1,3) 11021 format(i3,3f16.6,' 0 0 0 0 0 0 0 0.0') 112 write(3,1102)at(i),i,(bohr*xv(j,i),j=1,3),z(i) 1102 format(1x,a2,i3,3f16.6,f10.5) close(31) c return end subroutine normorbs implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*2 at character*1 ty parameter (nat0=96,nor=80,nort=500) integer*4 u0 parameter (u0=28) common/atomss/q(nat0),zsum,xv(3,nat0),multis common/works/ 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor),cij(nort,nort),bij(nort,nort), 4e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori, 6nop(nat0),ty(nat0,nor),ish(nat0,nor), 7at(nat0) common/spdfghs/ccs(nat0,nor),ccp(nat0,nor),ccd(nat0,nor,2), 1ccf(nat0,nor,3),cci_G(nat0,nor,15),cci_H(nat0,nor,21), 2cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) common/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half,a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(u0,nort),b6(u0,u0),aai(u0,u0),aaj(u0,u0), 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0) logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lnno,lg9 common/opts/lopt(300),icha equivalence (lopt(1),lsprint),(lopt(2),lpprint),(lopt(3),lqprint), 1(lopt(4 ),ld5 ),(lopt(5),levprint),(lopt(6),lvprint), 1(lopt(7 ),llprint),(lopt(8),lf7 ),(lopt(9),ld2 ), 1(lopt(10),SCALED ),(lopt(11),TTT ),(lopt(12),RAMAN ), 1(lopt(13),XYZ ),(lopt(14),ABINI ),(lopt(15),cc5 ), 1(lopt(16),lan ),(lopt(17),ldo ),(lopt(18),lvv ), 1(lopt(19),lsp ),(lopt(20),lex ),(lopt(21),lscf ), 1(lopt(22),lres ),(lopt(23),lcndo ),(lopt(24),fspec ), 1(lopt(25),lrec ),(lopt(26),lc5 ),(lopt(27),lg94 ), 1(lopt(57),lnno ),(lopt(66),lg9) dimension anorm(nort) c call inipoly i0=0 call initfd(aa1,aad,aaf) c tol=100.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,nao anorm(i)=z0 do 10 j=1,u0 10 sv(j,i)=z0 write(3,4000)nao,noa write(6,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C c first loop, determine norms: io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppps(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) if(npi.gt.u0)call report('npi > u0') jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppps(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi coe=con*cci(ioo) sx=oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) sy=oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) sz=oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) 71 b6(ioo,joo)=sx*sy*sz*coe c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif c 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, then write them do 33 ii=1,npii iorb=ii+io 33 anorm(iorb)=sv(ii,iorb) io=io+npii do 9 ii=1,u0 do 9 jj=1,nao 9 sv(ii,jj)=z0 endif 3 continue 2 continue c write(3,*)'atomic orbital norms before normalization:' write(3,3000)(anorm(i),i=1,nao) 3000 format(6f12.6) c c second loop, normalize C atom i, orbital io io=0 fac=1.0d0/sqrt(anorm(1)) do 8 ia=1,noa do 6 ip=1,nop(ia) call nppps(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) ccs(ia,ip )=ccs(ia,ip )*fac ccp(ia,ip )=ccp(ia,ip )*fac ccd(ia,ip,1)=ccd(ia,ip,1)*fac ccd(ia,ip,2)=ccd(ia,ip,2)*fac ccf(ia,ip,1)=ccf(ia,ip,1)*fac ccf(ia,ip,2)=ccf(ia,ip,2)*fac ccf(ia,ip,3)=ccf(ia,ip,3)*fac do i1=1,15 cci_G(ia,ip,i1)=cci_G(ia,ip,i1)*fac enddo do i1=1,21 cci_H(ia,ip,i1)=cci_H(ia,ip,i1)*fac enddo do i1=1,28 cci_I(ia,ip,i1)=cci_I(ia,ip,i1)*fac enddo if(ish(ia,ip+1).gt.ish(ia,ip))then io=io+npii fac=1.0d0/sqrt(anorm(io+1)) endif 6 continue 8 continue return end subroutine rp(nmo,ct,cij,iao,A,nort) implicit none integer*4 nmo,iao,imo,ix,nort,ii real*8 ct(*),cij(nort,nort),A(3,3) do 51 imo=1,nmo do 5 ix=1,3 5 ct(ix)=cij(imo,iao+ix-1) do 51 ix=1,3 cij(imo,iao+ix-1)=0.0d0 do 51 ii=1,3 51 cij(imo,iao+ix-1)=cij(imo,iao+ix-1)+A(ix,ii)*ct(ii) return end subroutine rd(nmo,cij,A,iao,nort) implicit none integer*4 nmo,iao,imo,nort real*8 cij(nort,nort),A(3,3),xx,yy,zz,xy,xz,yx,zx,yz,zy c order xx,yy,zz,xy,xz,yz do 41 imo=1,nmo xx=cij(imo,iao ) yy=cij(imo,iao+1) zz=cij(imo,iao+2) xy=cij(imo,iao+3) yx=cij(imo,iao+3) xz=cij(imo,iao+4) zx=cij(imo,iao+4) yz=cij(imo,iao+5) zy=cij(imo,iao+5) cij(imo,iao)= A(1,1)*A(1,1)*xx+A(1,1)*A(1,2)*xy+A(1,1)*A(1,3)*xz 1 +A(1,2)*A(1,1)*yx+A(1,2)*A(1,2)*yy+A(1,2)*A(1,3)*yz 2 +A(1,3)*A(1,1)*zx+A(1,3)*A(1,2)*zy+A(1,3)*A(1,3)*zz cij(imo,iao+1)=A(2,1)*A(2,1)*xx+A(2,1)*A(2,2)*xy+A(2,1)*A(2,3)*xz 1 +A(2,2)*A(2,1)*yx+A(2,2)*A(2,2)*yy+A(2,2)*A(2,3)*yz 2 +A(2,3)*A(2,1)*zx+A(2,3)*A(2,2)*zy+A(2,3)*A(2,3)*zz cij(imo,iao+2)=A(3,1)*A(3,1)*xx+A(3,1)*A(3,2)*xy+A(3,1)*A(3,3)*xz 1 +A(3,2)*A(3,1)*yx+A(3,2)*A(3,2)*yy+A(3,2)*A(3,3)*yz 2 +A(3,3)*A(3,1)*zx+A(3,3)*A(3,2)*zy+A(3,3)*A(3,3)*zz cij(imo,iao+3)=A(1,1)*A(2,1)*xx+A(1,1)*A(2,2)*xy+A(1,1)*A(2,3)*xz 1 +A(1,2)*A(2,1)*yx+A(1,2)*A(2,2)*yy+A(1,2)*A(2,3)*yz 2 +A(1,3)*A(2,1)*zx+A(1,3)*A(2,2)*zy+A(1,3)*A(2,3)*zz cij(imo,iao+4)=A(1,1)*A(3,1)*xx+A(1,1)*A(3,2)*xy+A(1,1)*A(3,3)*xz 1 +A(1,2)*A(3,1)*yx+A(1,2)*A(3,2)*yy+A(1,2)*A(3,3)*yz 2 +A(1,3)*A(3,1)*zx+A(1,3)*A(3,2)*zy+A(1,3)*A(3,3)*zz 41 cij(imo,iao+5)=A(2,1)*A(3,1)*xx+A(2,1)*A(3,2)*xy+A(2,1)*A(3,3)*xz 1 +A(2,2)*A(3,1)*yx+A(2,2)*A(3,2)*yy+A(2,2)*A(3,3)*yz 2 +A(2,3)*A(3,1)*zx+A(2,3)*A(3,2)*zy+A(2,3)*A(3,3)*zz return end c subroutine rf(idel,nmo,cij,A,iao,nort,shell) c rotation of f,g, and h AO coefficients in MO implicit none integer*4 lmax,u0 character*1 shell parameter (lmax=5,u0=28) integer*4 nmo,iao,imo,nort,jj(lmax),uzer,k(lmax),j1,j2,j3,j4, 1j5,i,idel,isum,ix,iy,iz,npi,npii real*8 cij(nort,nort),A(3,3),findme,cci(21),aai(u0,u0), 1aa1(u0,u0),aad(u0,u0),aaf(u0,u0) real*8,allocatable::cold(:) integer*4 nxi(u0),nyi(u0),nzi(u0) call initfd(aa1,aad,aaf) call nppp(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,.false.,.false.,iao,1) allocate(cold(idel)) do 41 imo=1,nmo do 11 i=1,idel 11 cold(i)=cij(imo,iao+i-1) do 41 i=1,idel ix=nxi(i) iy=nyi(i) iz=nzi(i) isum=ix+iy+iz if(isum.gt.lmax)call report('not implemented') if(isum.eq.0)return cij(imo,iao+i-1)=0.0d0 k(1)=uzer(ix,iy,iz) if(isum.eq.1)then do 1 j1=1,3 jj(1)=j1 1 cij(imo,iao+i-1)=cij(imo,iao+i-1)+A(k(1),j1) 1 *findme(jj,1,nxi,nyi,nzi,cold,idel) goto 41 endif k(2)=uzer(ix,iy,iz) if(isum.eq.2)then do 2 j1=1,3 do 2 j2=1,3 jj(1)=j1 jj(2)=j2 2 cij(imo,iao+i-1)=cij(imo,iao+i-1)+A(k(1),j1)*A(k(2),j2) 1 *findme(jj,2,nxi,nyi,nzi,cold,idel) goto 41 endif k(3)=uzer(ix,iy,iz) if(isum.eq.3)then do 3 j1=1,3 do 3 j2=1,3 do 3 j3=1,3 jj(1)=j1 jj(2)=j2 jj(3)=j3 3 cij(imo,iao+i-1)=cij(imo,iao+i-1)+A(k(1),j1)*A(k(2),j2) 1 *A(k(3),j3)*findme(jj,3,nxi,nyi,nzi,cold,idel) goto 41 endif k(4)=uzer(ix,iy,iz) if(isum.eq.4)then do 4 j1=1,3 do 4 j2=1,3 do 4 j3=1,3 do 4 j4=1,3 jj(1)=j1 jj(2)=j2 jj(3)=j3 jj(4)=j4 4 cij(imo,iao+i-1)=cij(imo,iao+i-1)+A(k(1),j1)*A(k(2),j2) 1 *A(k(3),j3)*A(k(4),j4)*findme(jj,4,nxi,nyi,nzi,cold,idel) goto 41 endif if(isum.eq.5)then k(5)=uzer(ix,iy,iz) do 5 j1=1,3 do 5 j2=1,3 do 5 j3=1,3 do 5 j4=1,3 do 5 j5=1,3 jj(1)=j1 jj(2)=j2 jj(3)=j3 jj(4)=j4 jj(5)=j5 5 cij(imo,iao+i-1)=cij(imo,iao+i-1)+A(k(1),j1)*A(k(2),j2)*A(k(3), 1 j3)*A(k(4),j4)*A(k(5),j5)*findme(jj,5,nxi,nyi,nzi,cold,idel) goto 41 endif 41 continue return end c ============================================================== subroutine getbasis(nao,nat,nop,nor,ish,ld5,lf7, 1 ty,cc,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,alpha,nat0, 1 nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii,rcp) implicit none integer*4 nao,nat,i,nop(*),k,nor,nat0,ish(nat0,nor), 1ishell,ip,ndao,ii,L,nb,ix,iy,iz,nxig(15),nyig(15), 1nzig(15),nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28) real*8 sf,af,ccs(nat0,nor),ccp(nat0,nor),cc(nat0,nor), 1ccd(nat0,nor,2),ccf(nat0,nor,3),cci_G(nat0,nor,15), 1cci_H(nat0,nor,21),cci_I(nat0,nor,28),alpha(nat0,nor), 1a,SNO,pi,bohr,debye,SNP,SND,SNF1,SNF2,SNF3,SNDD, 1rcp(nat0,nor),ein character*2 s2 character*80 s80 character*1 ty(nat0,nor) logical ld5,lf7 common/const/pi,bohr,debye nao=0 do 3 i=1,nat nop(i)=0 do 1 k=1,nor 1 ish(i,k)=32000 ishell=0 read(2,*) 50000 read(2,280)s80 280 format(a80) s2=s80(2:3) read(s80(4:len(s80)),*)ip,sf ishell=ishell+1 ndao=-90 if(s2.eq.'SP')ndao=4 if(s2.eq.'P ')ndao=3 if(s2.eq.'S ')ndao=1 if(s2.eq.'D '.and.ld5)ndao=5 if(s2.eq.'D '.and..not.ld5)ndao=6 if(s2.eq.'F '.and..not.lf7)ndao=10 if(s2.eq.'F '.and.lf7)ndao=7 if(s2.eq.'G ')ndao=15 if(s2.eq.'H ')ndao=21 if(s2.eq.'I ')ndao=28 if(ndao.eq.-90)call report('Unknown shell '//s2//'!') nao=nao+ndao do 2 ii=1,ip nop(i)=nop(i)+1 k=nop(i) if(k.gt.nor)call report('Too many primitives on one atom') if(s2.eq.'S ')ty(i,k)='S' if(s2.eq.'SP')ty(i,k)='L' if(s2.eq.'D ')ty(i,k)='D' if(s2.eq.'P ')ty(i,k)='P' if(s2.eq.'F ')ty(i,k)='F' if(s2.eq.'G ')ty(i,k)='G' if(s2.eq.'H ')ty(i,k)='H' if(s2.eq.'I ')ty(i,k)='I' if(s2.eq.'SP')then read(2,*)af,ccs(i,k),ccp(i,k) c write(6,*)af,ccs(i,k),ccp(i,k) else read(2,*)af,cc(i,k) c write(6,*)af,cc(i,k) endif af=af*sf**2 alpha(i,k)=af ish(i,k)=ishell c SNO normalization factor to S: a=dble(2)*af/pi SNO=sqrt(sqrt(a*a*a)) c c SNP normalization factor for P: a=af a=dble(128)*a*a*a*a*a/(pi*pi*pi) SNP=sqrt(sqrt(a)) c c SND,SNDD normalization factors for D: a=af a=dble(2048)*a*a*a*a*a*a*a/(pi*pi*pi) SND=sqrt(sqrt(a)) SNDD=SND/sqrt(dble(3)) c c SNFx normalization factors for F(fxxx,fxxy,fxyz): a=dble(2)*af a=sqrt( dble(8)*sqrt(a*a*a*a*a*a*a*a*a)/sqrt(pi*pi*pi) ) SNF1=a/sqrt(dble(15)) SNF2=a/sqrt(dble(3)) SNF3=a c c Higher-anglular momentum functions (G,H,I): c G: a=af L=4 nb=0 do 5 ix=0,L do 5 iy=0,L-ix do 5 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxig(nb)=ix nyig(nb)=iy nzig(nb)=iz cci_G(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 5 continue c c H: L=5 nb=0 do 4 ix=0,L do 4 iy=0,L-ix do 4 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxih(nb)=ix nyih(nb)=iy nzih(nb)=iz cci_H(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 4 continue c c I: L=6 nb=0 do 6 ix=0,L do 6 iy=0,L-ix do 6 iz=0,L-ix-iy if(ix+iy+iz.eq.L)then nb=nb+1 nxii(nb)=ix nyii(nb)=iy nzii(nb)=iz cci_I(i,k,nb)=cc(i,k) 1 /dsqrt(ein(2*ix,a)*ein(2*iy,a)*ein(2*iz,a)) endif 6 continue c ccd(i,k,1)=cc(i,k)*SNDD ccd(i,k,2)=cc(i,k)*SND c ccf(i,k,1)=cc(i,k)*SNF1 ccf(i,k,2)=cc(i,k)*SNF2 ccf(i,k,3)=cc(i,k)*SNF3 c c if(s2.eq.'SP')then c store exp coef. for control listing: cc(i,k)=ccs(i,k) rcp(i,k)=ccp(i,k) ccs(i,k)=ccs(i,k)*SNO ccp(i,k)=ccp(i,k)*SNP else ccs(i,k)=cc(i,k)*SNO ccp(i,k)=cc(i,k)*SNP endif if(ty(i,k).eq.'D')cc(i,k)=ccd(i,k,1)/SNDD 2 if(ty(i,k).eq.'P')cc(i,k)=ccp(i,k)/SNP read(2,3336)s2 3336 format(1x,a2) if(s2.ne.'**')then backspace 2 goto 50000 endif 3 continue write(6,*) 'G94 standard basis loaded from output' write(6,*) nao,' AOs' return end c ============================================================ subroutine fillvar(nat0,nor,nort,ifn,cc,alpha,rcp,cij,bij, 1e,be,nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,noa,nori,am, 6nop,ty,ish,at,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,nxig,nyig,nzig, 2nxih,nyih,nzih,nxii,nyii,nzii,z,zsum,xv,lopt,icha,multi,ipse) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*(*) ifn character*2 at,s2,symbols(106) character*1 ty character*80 s80 common/const/pi,bohr,debye dimension cc(nat0,nor),alpha(nat0,nor),rcp(nat0,nor), 1cij(nort,nort),bij(nort,nort),e(nort),be(nort), 5nop(nat0),ty(nat0,nor),ish(nat0,nor),at(nat0),ccs(nat0,nor), 1ccp(nat0,nor),ccd(nat0,nor,2),ccf(nat0,nor,3),cci_G(nat0,nor,15), 1cci_H(nat0,nor,21),cci_I(nat0,nor,28),nxig(15),nyig(15),nzig(15), 2nxih(21),nyih(21),nzih(21),nxii(28),nyii(28),nzii(28),z(nat0), 1xv(3,nat0),am(nat0),ipse(nat0) logical lf7,lzmat,ld5,lconfig,lopt(300) character*80 key dimension amas(89) 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/ character*3 s3 data symbols/'H ', 'He', 2 'Li','Be','B ','C ','N ','O ','F ','Ne', 3 'Na','Mg','Al','Si','P ','S ','Cl','Ar', 4 'K ','Ca', 4 'Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX', 4'XX','XX','XX','XX','XX','XX','XX','XX','XX','XX'/ c gaussian checkpoint: c shell types,number of primitives per shell,shell to atom map, c npe/rpe-primitive exponents c ncc/rcc-contraction coefficients c ncl/rcl-contraction coefficients P in SP integer*4 nst,ist(nort),npps,ipps(nort),nstam,istam(nort), 1npe,ncc,ncl real*8 rpe(nort),rcc(nort),rcl(nort) do i=1,nort rcl(i)=0.0d0 enddo c nst=0 npps=0 nstam=0 npe=0 ncc=0 c ld5=lopt(4 ) lf7=lopt(8) lzmat=lopt(28) lconfig=lopt(65) nmo=0 open(77,file='COMPARE.OPT') 773 read(77,771,end=772,err=772)s3 771 format(a3) if(s3.eq.'LD5')read(77,*)ld5 if(s3.eq.'LF7')read(77,*)lf7 if(s3.eq.'LZM')read(77,*)lzmat if(s3.eq.'NMO')read(77,*)nmo goto 773 772 close(77) do i=1,nat0 do j=1,nor ish(i,j)=0 enddo enddo c try to get some dimensions first open(8,file=ifn) 801 read(8,800,end=8001,err=8001)s80 800 format(a80) if(s80(1:9).eq.' Charge =')then read(s80,878)icha,multi 878 format(9x,i3,15x,i2) write(6,*)' charge read from output ',icha,multi if(multi.gt.1)lconfig=.true. endif if(s80(7:22).eq.' basis functions ')then open(16,file='FMT.SCR') write(16,1009)s80(1:6) 1009 format(a6) rewind 16 read(16,*)nao close(16) goto 8002 endif goto 801 8001 close(8) call report('NMO/NAO not found') 8002 write(6,*)'nao : ',nao if(nmo.ne.0)then write(6,*)'nmo from option file: ',nmo else nmo=nao write(6,*)'nmo set to nao: ',nmo endif nori=nmo close(8) c ieigen=0 igeo=0 ibasis=0 c AO/MO/geometry reading open(2,file=ifn) 111 read(2,1001,end=666)key 1001 format(a80) if(key(1:35).eq.' Molecular Orbital Coefficients'.or. 1 key(2:11).eq.'Alpha MOs:') then write(6,1001)key ieigen=ieigen+1 c eigenvectors listed either as a result of IOP(6/7=1) or IOP(5/33=1) nori=ndo if(nmo.gt.0)nori=nmo if(nao.eq.0)nao=nmo call readorb(key,nori,e,cij,nao,nort) do 377 i1=1,nmo be(i1)=e(i1) do 377 i2=1,nao 377 bij(i1,i2)=cij(i1,i2) read(2,1001)key if(key(2:10).eq.'Beta MOs:')then write(6,1001)key call readorb(key,nori,be,bij,nao,nort) endif endif c if(key(2:10).eq.' SCF done:')write(3,*)key c c GEOMETRY if(key(20:35).eq.'Standard orienta'.or. 1 key(26:41).eq.'Standard orienta')then ig98=0 if(key(26:41).eq.'Standard orienta')ig98=1 write(3,1001)key if(lzmat)then write(3,*)'geometry skipped because lzmat is defined' goto 877 endif goto 776 endif if(key(19:35).eq.'Z-Matrix orientat'.or. 1 key(26:42).eq.'Z-Matrix orientat'.or. 1 key(27:43).eq.'Input orientation'.or. 1 key(20:36).eq.'Input orientation')then ig98=0 if(key(26:42).eq.'Z-Matrix orientat')ig98=1 if(key(27:43).eq.'Input orientation')ig98=1 write(3,1001)key if(.not.lzmat)then write(3,*)'geometry skipped because lzmat is not defined' goto 877 endif goto 776 endif goto 877 776 igeo=igeo+1 do 1031 i=1,4 1031 read(2,*) i=0 1011 i=i+1 read(2,3335)s2 3335 format(a2) if(s2.ne.' -')then backspace 2 if(ig98.eq.0)read(2,*)iza,iza,xv(1,i),xv(2,i),xv(3,i) if(ig98.eq.1)read(2,*)iza,iza,itype,xv(1,i),xv(2,i),xv(3,i) z(i)=dble(iza) am(i)=amas(iza) at(i)=symbols(iza) if(itype.eq.1000)then z(i)=0.0d0 at(i)='Bq' endif goto 1011 endif nat=i-1 write(6,*)'nat = ',nat if(nat.gt.nat0)then write(6,*)'nat0 = ',nat0 call report('too many atoms') endif zsum=0.0d0 do 1012 i=1,nat zsum=zsum+z(i) do 1012 j=1,3 1012 xv(j,i)=xv(j,i)/bohr write(6,*)'G94 geometry read in,',nat,' atoms' write(3,*)'G94 geometry read in,',nat,' atoms' noa=nat noe=nint(zsum)-icha if(lconfig)then write(3,*)' Open shell ' write(6,*)' Open shell ',multi,noe ndo=0 nsa=(multi+noe-1)/2 nsb=noe-nsa else write(3,*)' Closed shell ' ndo=noe/2 nsa=0 nsb=0 endif ichm=icha mtp=1 877 continue c c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB c BASIS GFInput if( key(2:46).eq.'Basis set in the form of general basis input:' 1.or.key(2:46).eq.'AO basis set in the form of general basis inp' 2 )then write(6,*)key write(6,6000)ld5,lf7 write(3,6000)ld5,lf7 6000 format(' LD5 LF7: ',2L2) ibasis=ibasis+1 call getbasis(nao,nat,nop,nor,ish,ld5,lf7, 1 ty,cc,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,alpha,nat0, 1 nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii,rcp) if(nao.gt.nort)call report('too many orbitals') endif c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB if( key(40:50).eq.'Pseudopoten')then call rdps(nat,ipse) do 1 i=1,nat if(ipse(i).ne.0)then write(6,6005)i,ipse(i) 6005 format(' atom ',i4,' atomic charge reduced to',i4) zsum=zsum-z(i)+dble(ipse(i)) z(i)=dble(ipse(i)) endif 1 continue noe=nint(zsum)-icha if(lconfig)then write(3,*)' Open shell ' write(6,*)' Open shell ',multi,noe ndo=0 nsa=(multi+noe-1)/2 nsb=noe-nsa else write(3,*)' Closed shell ' ndo=noe/2 nsa=0 nsb=0 endif ichm=icha mtp=1 endif c GAUSSIAN FORMATTED CHECKPOINT: c FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF if(key(1:22).eq.'Alpha Orbital Energies')then write(6,1001)key write(6,*)nao,nmo read(2,2009)(e(i1),i1=1,nmo) do 10 i1=1,nmo 10 be(i1)=e(i1) endif if(key(1:21).eq.'Alpha MO coefficients')then write(6,1001)key write(6,*)nao,nmo ieigen=ieigen+1 if(nao.eq.0)nao=nmo read(2,2009)((cij(i1,i2),i2=1,nao),i1=1,nmo) 2009 format(5E16.8) write(6,*) 'Eigenvectors read in from verbose checkpoint... ' do 10041 i=1,nmo do 10041 j=1,nao 10041 bij(i,j)=cij(i,j) endif if(key(1:20).eq.'Beta MO coefficients')then write(6,1001)key read(2,2009)((bij(i1,i2),i2=1,nao),i1=1,nmo) endif if(key(1:21).eq.'Beta Orbital Energies')then write(6,1001)key write(6,*)nao,nmo read(2,2009)(be(i1),i1=1,nmo) endif c GEOMETRY if(key(1:15).eq.'Nuclear charges')then read(2,2009)(z(i),i=1,noa) zsum=0.0d0 do 9 i=1,noa zsum=zsum+z(i) 9 at(i)=symbols(int(z(i))) endif if(key(1:17).eq.'Current cartesian')then write(3,1001)key igeo=igeo+1 c in bohrs read(2,2009)((xv(ix,i),ix=1,3),i=1,noa) endif if(key(1:15).eq.'Number of atoms')then read(key(50:61),*)nat if(nat.gt.nat0)call report('too many atoms') noa=nat endif if(key(1:6).eq.'Charge')then read(key(50:61),*)icha ichm=icha endif if(key(1:19).eq.'Number of electrons')read(key(50:61),*)noe if(key(1:25).eq.'Number of alpha electrons')read(key(50:61),*)nsa if(key(1:24).eq.'Number of beta electrons')then read(key(50:61),*)nsb if(nsa.eq.nsb)then lconfig=.false. ndo=nsa nsa=0 nsb=0 write(6,*)' Closed shell ' else write(6,*)' Open shell ',nsa,nsb ndo=0 endif mtp=1 endif if(key(1:25).eq.'Number of basis functions')read(key(50:61),*)noa c if(key(1:20).eq.'Number of contracted')then c nshells endif if(key(1:11).eq.'Shell types')then read(key(50:61),*)nst write(6,1001)key if(nst.gt.nort)call report('too many shell types') read(2,2003)(ist(i),i=1,nst) 2003 format(6i12) endif if(key(1:30).eq.'Number of primitives per shell')then read(key(50:61),*)npps write(6,1001)key if(npps.gt.nort)call report('too many primitives per shell') read(2,2003)(ipps(i),i=1,npps) endif if(key(1:17).eq.'Shell to atom map')then read(key(50:61),*)nstam write(6,1001)key if(nstam.gt.nort)call report('too many shell shell to atom') read(2,2003)(istam(i),i=1,nstam) endif if(key(1:19).eq.'Primitive exponents')then read(key(50:61),*)npe write(6,1001)key if(npe.gt.nort)call report('too many prim exp') read(2,2004)(rpe(i),i=1,npe) 2004 format(5E16.8) endif if(key(1:24).eq.'Contraction coefficients')then read(key(50:61),*)ncc write(6,1001)key if(ncc.gt.nort)call report('too many contract coeff') read(2,2004)(rcc(i),i=1,ncc) read(2,1001)key if(key(1:31).ne.'P(S=P) Contraction coefficients')then ncl=0 if(nst.gt.0.and.npps.gt.0.and.nstam.gt.0.and.npe.gt.0. 1 and.ncc.gt.0)then call getbasisfck(nao,nat,nop,nor,ish,ld5,lf7, 1 ty,cc,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,alpha,nat0, 2 nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii,rcp, 3 nst,ist,ipps,istam,rpe,rcc,rcl) ibasis=ibasis+1 else write(6,*)nst,npps,nstam,npe,ncc,ncl call report('incomplete basis set input') endif endif endif if(key(1:31).eq.'P(S=P) Contraction coefficients')then read(key(50:61),*)ncl if(ncl.gt.nort)call report('too many contract coeff L') read(2,2004)(rcl(i),i=1,ncl) if(nst.gt.0.and.npps.gt.0.and.nstam.gt.0.and.npe.gt.0. 1 and.ncc.gt.0.and.ncl.gt.0)then call getbasisfck(nao,nat,nop,nor,ish,ld5,lf7, 1 ty,cc,ccs,ccp,ccd,ccf,cci_G,cci_H,cci_I,alpha,nat0, 2 nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii,rcp, 3 nst,ist,ipps,istam,rpe,rcc,rcl) ibasis=ibasis+1 else write(6,*)nst,npps,nstam,npe,ncc,ncl call report('incomplete basis set input') endif endif c (end of Gaussian checkpoint) c goto 111 c 666 close(2) noa=nat write(3,4004)ieigen,igeo,ibasis write(6,4004)ieigen,igeo,ibasis 4004 format( 2 ' Eigenvectors read',I4,' times',/, 3 ' Geometry read ',I4,' times',/, 4 ' Basis read ',I4,' times',/) if(ieigen.ne.1.and.igeo.ne.1.and.ibasis.ne.1) 1call report('unbalanced input') return end c ============================================================ subroutine readorb(key,nori,e,cij,nao,cdimension) implicit none integer*4 i,j,i1,j1,i2,nao,nori,ic,cdimension real*8 e(*),cij(cdimension,cdimension) logical number character*(*) key character*160 s160 i=1 read(2,*) 377 j=i+4 if(key(8:11).ne.'MOs:'.and.key(7:10).ne.'MOs:')read(2,*) if(j.gt.nori)j=nori c c to enable reading of various formats, use s160: read(2,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(2,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(3,*)'Corrected:',ic write(3,160)s160 endif 39 read(s160(22:160),*,err=99)(cij(i1,i2),i1=i,j) if(j.lt.nori)then i=i+5 read(2,*) 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 rdps(nat,ipse) c read pseudopotential parameters implicit none character*29 s29 integer*4 i,nat,ipse(*) read(2,*) read(2,*) read(2,*) read(2,*) i=0 1 read(2,200)s29 200 format(a29) if(s29(15:16).eq.' ')goto 1 if(s29(15:16).eq.'==')return i=i+1 if(i.gt.nat)call report('pseudopotential reading error') if(s29(28:29).eq.' ')then ipse(i)=0 else read(s29(20:29),*)ipse(i) endif goto 1 end 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 ============================================================ function uzer(ix,iy,iz) implicit none integer*4 ix,iy,iz,uzer if(ix.gt.0)then ix=ix-1 uzer=1 else if(iy.gt.0)then iy=iy-1 uzer=2 else if(iz.gt.0)then iz=iz-1 uzer=3 else call report('nothing left in uzer') endif endif endif return end c ============================================================ subroutine pretav(x,y,z,jj,m) c 1 1 2 -> x^2 y^1 z^0 c 1 3 2 -> x^1 y^1 z^1 etc. implicit none integer*4 jj(*),m,x,y,z,i x=0 y=0 z=0 do 1 i=1,m if(jj(i).eq.1)x=x+1 if(jj(i).eq.2)y=y+1 1 if(jj(i).eq.3)z=z+1 return end c ============================================================ function findme(jj,m,nxi,nyi,nzi,c,n) implicit none real*8 findme,c(*) integer*4 jj(*),m,nxi(*),nyi(*),nzi(*),x,y,z,n,i call pretav(x,y,z,jj,m) do 1 i=1,n if(x.eq.nxi(i).and.y.eq.nyi(i).and.z.eq.nzi(i))then findme=c(i) return endif 1 continue write(6,*)(jj(i),i=1,m) do 2 i=1,n 2 write(6,*)i,nxi(i),nyi(i),nzi(i) call report('find me could not find me') end c ============================================================ subroutine vz(v,n) implicit none real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end c ============================================================