program findrydberg c c reads coordinates, basis, orbitals and exc states c from BIG.OUT, BIG.SP c and estimates c IMPLICIT none integer*4 nat0,nt0,isu0,nort,imo,jj,nor,u0,nb0 parameter (nat0=96,nt0=160,isu0=10000,nort=800,nor=80,u0=28, 1nb0=200) 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),it,ii,ia,ix,noa,multi,kk, 1isn(nt0),icha,npl,nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,nori, 2nop,ish,i,io,ia1,ia2,ir1,ir2 real*8 en(nt0),ws(nt0,isu0),qtot,enev(nt0),zq,zsum,xv,cb(3), 1du(nt0,4),qxx,qyy,qzz,rxx,ryy,rzz,rxxg,ryyg,rzzg, 2r2,nxx,nxy,nxz,nyy,nyz,nzz,pi,bohr,debye,xxt,xyt,xzt,yyt,ani, 2yzt,zzt,fp,qxy,qxz,qyz,rxyg,rxzg,ryzg,bxx,byy,bzz,bxy,bxz,byz, 3l2,b2,l22,a2,z0,x,y,z,psioc,psivi,f,r,forbxp4,t, 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0),cc,alpha,rcp,cij,bij,e,be character*1 ty character*2 at,rty character*10 oname common/atoms/zq(nat0),zsum,xv(3,nat0),multi common/qq/qxx(nort,nort),qyy(nort,nort),qzz(nort,nort), 1qxy(nort,nort),qxz(nort,nort),qyz(nort,nort),l2(nort,nort) common/qb/bxx(nort,nort),byy(nort,nort),bzz(nort,nort), 1bxy(nort,nort),bxz(nort,nort),byz(nort,nort),b2(nort,nort) logical lopt common/opts/lopt(300),icha,npl(10) common/const/pi,bohr,debye 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 cut(3,nb0,nort) common/cuts/cut c write(6,6767) 6767 format(' Electronic transition analysis',/, 2 ' Input: BIG.OUT ... RPA outputs',/, 3 ' BIG.SP ... SP outputs',/, 3 ' RYDBERG.OPT ... options ',/, 6 ' Output: ',/, 7 ' ',/, 7 ' FINDRYDBERG.OUT ... some controls') C open(3,file='FINDRYDBERG.OUT') call init z0=0.0d0 call initfd(aa1,aad,aaf) call getRAOMO('BIG.SP') call getexc('BIG.OUT',it,en,s1,s2,ws,isn,du,enev) call writecon if(lopt(1))then call normorb write(6,*)'AOs have been renormalized' endif c c one electron integrals for starting point call overlap6 call charge(1) if(lopt(1))then call ortog write(6,*)'MOs have been orthonormalized' call charge(0) endif call mu call thetarr nxx=0.0d0 nxy=0.0d0 nxz=0.0d0 nyy=0.0d0 nyz=0.0d0 nzz=0.0d0 qtot=0.0d0 do 1 ix=1,3 1 cb(ix)=0.0d0 do 2 ia=1,noa qtot=qtot+zq(ia) nxx =nxx +zq(ia)*xv(1,ia)*xv(1,ia) nxy =nxy +zq(ia)*xv(1,ia)*xv(2,ia) nxz =nxz +zq(ia)*xv(1,ia)*xv(3,ia) nyy =nyy +zq(ia)*xv(2,ia)*xv(2,ia) nyz =nyz +zq(ia)*xv(2,ia)*xv(3,ia) nzz =nzz +zq(ia)*xv(3,ia)*xv(3,ia) do 2 ix=1,3 2 cb(ix)=cb(ix)+xv(ix,ia)*zq(ia) do 3 ix=1,3 3 cb(ix)=cb(ix)/qtot c rxxg=0.0d0 rxyg=0.0d0 rxzg=0.0d0 ryyg=0.0d0 ryzg=0.0d0 rzzg=0.0d0 l22=0.0d0 do 51 imo=1,nsa rxxg=rxxg-qxx(imo,imo) rxyg=rxyg-qxy(imo,imo) rxzg=rxzg-qxz(imo,imo) ryyg=ryyg-qyy(imo,imo) ryzg=ryzg-qyz(imo,imo) rzzg=rzzg-qzz(imo,imo) 51 l22=l22-l2(imo,imo) do 52 imo=1,nsb rxxg=rxxg-bxx(imo,imo) rxyg=rxyg-bxy(imo,imo) rxzg=rxzg-bxz(imo,imo) ryyg=ryyg-byy(imo,imo) ryzg=ryzg-byz(imo,imo) rzzg=rzzg-bzz(imo,imo) 52 l22=l22-b2(imo,imo) do 5 imo=1,ndo rxxg=rxxg-qxx(imo,imo)*2.0d0 rxyg=rxyg-qxy(imo,imo)*2.0d0 rxzg=rxzg-qxz(imo,imo)*2.0d0 ryyg=ryyg-qyy(imo,imo)*2.0d0 ryzg=ryzg-qyz(imo,imo)*2.0d0 rzzg=rzzg-qzz(imo,imo)*2.0d0 5 l22=l22-l2(imo,imo)*2.0d0 xxt=rxxg+nxx xyt=rxyg+nxy xzt=rxzg+nxz yyt=ryyg+nyy yzt=ryzg+nyz zzt=rzzg+nzz fp=debye*bohr l22=(dsqrt(1.0d0+4.0d0*dabs(l22))-1)/2.0d0 write(3,6005)cb,nxx,nxy,nxz,nyy,nyz,nzz, 1rxxg,rxyg,rxzg,ryyg,ryzg,rzzg, 1xxt,xyt,xzt,yyt,yzt,zzt, 1xxt*fp,xyt*fp,xzt*fp,yyt*fp,yzt*fp,zzt*fp,l22 write(6,6005)cb,nxx,nxy,nxz,nyy,nyz,nzz, 1rxxg,rxyg,rxzg,ryyg,ryzg,rzzg, 1xxt,xyt,xzt,yyt,yzt,zzt, 1xxt*fp,xyt*fp,xzt*fp,yyt*fp,yzt*fp,zzt*fp,l22 6005 format('Ground state:',/, 1'Molecular center ',3f12.6,/, 1'Quadrupoles:',/, 1'Nuclear :',6f10.2,/, 1'Electronic :',6f10.2,/, 1'Total :',6f10.2,' (au)',/, 1'Total :',6f10.2,' (D*A)',/, 1'l :',1f10.2,' (au)',/) write(6,6004) write(3,6004) 6004 format('Electronic Transitions/exc states:',/, 1'number energy i -> j cij (main)) dxx dyy dzz dr2 l2') do 4 ii=1,it rxx=0.0d0 ryy=0.0d0 rzz=0.0d0 a2=0.0d0 do 6 jj=1,isn(ii) do 6 kk=1,isn(ii) if(s1(ii,jj).eq.s1(ii,kk))then rxx=rxx+ws(ii,jj)*ws(ii,kk)*qxx(s2(ii,jj),s2(ii,kk)) ryy=ryy+ws(ii,jj)*ws(ii,kk)*qyy(s2(ii,jj),s2(ii,kk)) rzz=rzz+ws(ii,jj)*ws(ii,kk)*qzz(s2(ii,jj),s2(ii,kk)) a2=a2 +ws(ii,jj)*ws(ii,kk)*l2(s2(ii,jj),s2(ii,kk)) endif if(s2(ii,jj).eq.s2(ii,kk))then rxx=rxx+ws(ii,jj)*ws(ii,kk)*qxx(s1(ii,jj),s1(ii,kk)) ryy=ryy+ws(ii,jj)*ws(ii,kk)*qyy(s1(ii,jj),s1(ii,kk)) rzz=rzz+ws(ii,jj)*ws(ii,kk)*qzz(s1(ii,jj),s1(ii,kk)) a2=a2 +ws(ii,jj)*ws(ii,kk)*l2(s1(ii,jj),s1(ii,kk)) endif 6 continue r2=rxx+ryy+rzz r=dsqrt(r2) ani=(rxx+ryy-rzz)/(rxx+ryy+rzz)*3.0d0 c density plots: if(ii.ge.npl(5).and.ii.le.npl(6))then write(oname,2006)ii 2006 format(i10) do 2003 io=1,3 if(io.eq.1)open(45,file=oname//'sdx.prn') if(io.eq.2)open(45,file=oname//'sdy.prn') if(io.eq.3)open(45,file=oname//'sdz.prn') do 31 i=1,npl(3)+1 x=z0 y=z0 z=z0 t=(-1.0d0+dble(2*(i-1))/dble(npl(3)))*r*dble(npl(4)) if(io.eq.1)x=t if(io.eq.2)y=t if(io.eq.3)z=t f=0.0d0 c if(npl(7).gt.0)then c do 62 jj=1,nsa c62 f=f+forbxp4(x,y,z,s2(ii,jj),1,aa1,aad,aaf)**2 c do 63 jj=1,nsb c63 f=f+forbxp4(x,y,z,s2(ii,jj),2,aa1,aad,aaf)**2 c do 64 jj=1,ndo c64 f=f+forbxp4(x,y,z,s2(ii,jj),1,aa1,aad,aaf)**2*2.0d0 c endif c do 61 jj=1,isn(ii) c psioc=forbxp4(x,y,z,s1(ii,jj),1,aa1,aad,aaf) c psivi=forbxp4(x,y,z,s2(ii,jj),1,aa1,aad,aaf) c if(npl(7).gt.0)f=f-ws(ii,jj)**2*psioc**2 c61 f=f+ws(ii,jj)**2*psivi**2 c c density change only do 61 jj=1,isn(ii) ia1=s1(ii,jj) ir1=s2(ii,jj) do 61 kk=1,isn(ii) ia2=s1(ii,kk) ir2=s2(ii,kk) if(ir1.eq.ir2)then if(ia1.eq.ia2)then c f=f+ws(ii,jj)*ws(ii,kk)*(cut(io,i,ir1)**2-cut(io,i,ia1)**2) c neglect original orbitals: f=f+ws(ii,jj)*ws(ii,kk)*cut(io,i,ir1)**2 else c f=f+ws(ii,jj)*ws(ii,kk)*cut(io,i,ia1)*cut(io,i,ia2) endif endif 61 if(ia1.eq.ia2.and.ir1.ne.ir2) 1 f=f+ws(ii,jj)*ws(ii,kk)*cut(io,i,ir1)*cut(io,i,ir2) 31 write(45,4545)t*bohr,f*t**2 4545 format(f7.2,f15.10) 2003 close(45) endif rty=' ' if(abs(ani-1.0d0).lt.0.1d0)rty=' S' c a2=(dsqrt(1.0d0+4.0d0*a2)-1.0d0)/2.0d0 write(3,6000)ii,enev(ii),s1(ii,1),s2(ii,1),ws(ii,1), 1rxx,ryy,rzz,r2,rty,a2 6000 format(i3,f8.2,i3,'->',i3,f5.2,3f9.1,1x,f10.1,A2,f8.1) 4 write(6,6000)ii,en(ii),s1(ii,1),s2(ii,1),ws(ii,1), 1rxx,ryy,rzz,r2,rty,a2 close(3) stop END c ============================================================ subroutine getRAOMO(ifn) c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c parameter (nat0=96,nor=80,nort=800) character*(*) ifn character*2 at,s2,symbols(106) character*1 ty character*160 s160 character*80 s80 common/const/pi,bohr,debye 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 lf7,lzmat,lg94,ld5,lconfig,lopt,lg9,lnorm,lst common/opts/lopt(300),icha,npl(10) equivalence (lopt(1),lnorm),(lopt(2),lst), 1(lopt(4 ),ld5 ),(lopt(8),lf7 ),(lopt(27),lg94 ), 1(lopt(28),lzmat ),(lopt(65),lconfig),(lopt(66),lg9) character*50 key 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 lzmat=.true. lst=.false. lnorm=.true. ld5=.true. lf7=.true. lg9=.true. do 8989 i=1,10 8989 npl(i)=0 npl(3)=200 c plot orbital cuts for mos within npl(1)..npl(2) c plot density cuts for mos within npl(5)..npl(6) c with npl(3) points c with npl(4)*ave r range c if npl(7)<>0 total density open(77,file='RYDBERG.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.'LG9')read(77,*)lg9 if(s3.eq.'LZM')read(77,*)lzmat if(s3.eq.'NOR')read(77,*)lnorm if(s3.eq.'STR')read(77,*)lst if(s3.eq.'NPL')read(77,*)i,npl(i) 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,*)nmo close(16) goto 8002 endif goto 801 8001 close(8) call report('NMO/NAO not found') 8002 write(6,*)'nmo set to ',nmo nao=nmo nori=nmo close(8) c if(lst)write(6,*)'fixed format' ieigen=0 igeo=0 ibasis=0 c AO/MO/geometry reading open(2,file=ifn) 111 read(2,1001,end=666)key 1001 format(a49) c checkpoint format: if(key(1:22).eq.'Alpha Orbital Energies')then write(6,1001)key read(2,*)(e(i),i=1,nmo) do i=1,nmo be(i)=e(i) enddo endif if(key(1:21).eq.'Alpha MO coefficients')then write(6,1001)key ieigen=ieigen+1 read(2,*)((cij(i,j),j=1,nmo),i=1,nao) do i=1,nao do j=1,nmo bij(j,i)=cij(j,i) enddo enddo endif if(key.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) i=1 read(2,*) 377 j=i+4 if(key(8:11).ne.'MOs:')read(2,*) if(j.gt.nmo)j=nmo c c to enable reading of various formats: if(lst)then c fixed format read(2,1160)(e(i1),i1=i,j) 1160 format(21x,5f10.4) do 88 i2=1,nao 88 read(2,1160)(cij(i1,i2),i1=i,j) else read(2,1600)s160 1600 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 j1=160,i1+1,-1 s160(j1:j1)=s160(j1-1:j1-1) enddo 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) c do 391 i2=1,nao read(2,1600)s160 read(s160(1:159),*,err=9888,end=9888)idum,(cij(i1,i2),i1=i,j) goto 391 9888 write(6,1600)s160 write(6,*)i,j,i2,nao,nmo call report('MO error reading') 391 continue endif if(j.lt.nmo)then i=i+5 read(2,*) goto 377 endif write(6,*) 'Eigenvectors read in from verbose output... ' do 1004 i=1,nmo be(i)=e(i) do 1004 j=1,nao 1004 bij(i,j)=cij(i,j) read(2,1001,end=666)key c beta-spin orbitals if(key(2:11).eq.'Beta MOs:') then write(6,1001)key i=1 read(2,*) 3771 j=i+4 if(j.gt.nmo)j=nmo c c write to scratch to enable reading of various formats: if(lst)then c fixed format read(2,1160)(be(i1),i1=i,j) do 881 i2=1,nao 881 read(2,1160)(bij(i1,i2),i1=i,j) else read(2,'(a)')s160 i1=22 3776 if(s160(i1:i1).eq.'-'.and.s160(i1-1:i1-1).ne.'E')then do j1=160,i1+1,-1 s160(j1:j1)=s160(j1-1:j1-1) enddo s160(i1:i1)=' ' i1=i1+1 endif i1=i1+1 if(i1.lt.160)goto 3776 read(s160(22:160),*)(be(i1),i1=i,j) do 3911 i2=1,nao 3911 read(2,*)idum,(bij(i1,i2),i1=i,j) endif if(j.lt.nmo)then i=i+5 read(2,*) goto 3771 endif write(6,*) 'Beta eigenvectors read in from verbose output... ' endif endif c c ENERGY 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)then read(2,*)idum,iza,xv(1,i),xv(2,i),xv(3,i) else read(2,*)idum,iza,itype,xv(1,i),xv(2,i),xv(3,i) endif z(i)=dble(iza) at(i)=symbols(iza) if(itype.eq.1000)then z(i)=0.0d0 at(i)='Bq' endif goto 1011 endif nat=i-1 if(nat.gt.nat0)then write(6,*)'nat = ',nat 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=int(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 cBasis set in the form of general basis input: c2345678901234567890123456789012345678901234567890 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 ibasis=ibasis+1 nao=0 do i=1,nat nop(i)=0 do k=1,nor ish(i,k)=32000 enddo ishell=0 read(2,*) 50000 read(2,3907)s2,ip,sf c write(6,3907)s2,ip,sf 3907 format(1x,a2,i4,f5.2) 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 '.and.lg9)ndao=9 if(s2.eq.'G '.and..not.lg9)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 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 ix=0,L do iy=0,L-ix do 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 enddo enddo enddo c H: L=5 nb=0 do ix=0,L do iy=0,L-ix do 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 enddo enddo enddo c I: L=6 nb=0 do ix=0,L do iy=0,L-ix do 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 enddo enddo enddo 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 if(ty(i,k).eq.'P')cc(i,k)=ccp(i,k)/SNP c if(ty(i,k).eq.'F')cc(i,k)=dum enddo read(2,3336)s2 3336 format(1x,a2) if(s2.ne.'**')then backspace 2 goto 50000 endif enddo write(6,*) 'G94 standard basis loaded from output' write(6,*) nao,' AOs' if(nao.gt.nort)call report('too many orbitals') endif c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB 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') end c ============================================================== 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,enev) IMPLICIT none integer*4 nt0,it,isu,isu0,ii,ii1,idumm character*(*) ifn PARAMETER (isu0=10000,nt0=160) real*8 pi,bohr,debye common/const/pi,bohr,debye c Limits: c nt0 -# of transitions on one chromophore 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),enev(*) integer s1(nt0,isu0),s2(nt0,isu0),nt,isn(*) 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')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,*)idumm,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 read(s80(15:18),*)it if(it.gt.nt0)call report('too many transitions') read(s80(35:43),*)enev(it) read(s80(49:55),*)en(it) isu=0 1 read(30,300)s80 if(nu(s80(7:7)))then isu=isu+1 if(isu.gt.isu0)call report('isu > isu0') read(s80,301)s1(it,isu),s2(it,isu),ws(it,isu) 301 format(i7,4x,i3,1x,f15.5) goto 1 endif if(nu(s80(8:8)))then isu=isu+1 if(isu.gt.isu0)call report('isu > isu0') read(s80,341)s1(it,isu),s2(it,isu),ws(it,isu) 341 format(i8,3x,i3,1x,f15.5) goto 1 endif isn(it)=isu 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 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 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=800) 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, 1a10 dimension s(nort,nort),oldvec(nort),oldveb(nort) c c read in overlap matrix for AO 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 write(6,*)' MO norms:' do 3 i=1,nmo c do 6 ii=1,nao oldveb(ii)=bij(i,ii) 6 oldvec(ii)=cij(i,ii) c do 4 j=1,i-1 aij=vecpro(i,j,cij,s,nao) aib=vecpro(i,j,bij,s,nao) 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(6,609)i,vecnor(oldvec,s,nao),vecnor(oldveb,s,nao) 609 format(i4,2f10.4) do 3 ii=1,nao bij(i,ii)=oldveb(ii)*bn 3 cij(i,ii)=oldvec(ii)*an write(3,*)' MOs have been re-ortonormalized' write(6,*)' MOs have been re-ortonormalized' return end c ============================================================= function vecpro(i,j,cij,s,nao) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nort=800) dimension cij(nort,nao),s(nort,nao) c sum=0.0D0 do 1 ii=1,nao d1=cij(i,ii) do 1 jj=1,nao 1 sum=sum+d1*cij(j,jj)*s(ii,jj) vecpro=sum return end c ============================================================= function vecnor(c,s,n) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nort=800) 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 subroutine mu c calculates dipole moment implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=96,nor=80,nort=800) 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, 1a10 common/atoms/z(nat0),zsum,xv(3,nat0),multi dimension pix(nort,nort),piy(nort,nort),piz(nort,nort) amux=z0 amuy=z0 amuz=z0 open(33,file='P.SCR',form='unformatted') do 1 i=1,nao read(33)(pix(i,j),j=1,i) read(33)(piy(i,j),j=1,i) read(33)(piz(i,j),j=1,i) do 1 j=1,i pix(j,i)=pix(i,j) piy(j,i)=piy(i,j) 1 piz(j,i)=piz(i,j) close(33) call trfo(pix,nmo,nao,cij) call trfo(piy,nmo,nao,cij) call trfo(piz,nmo,nao,cij) do 111 imo=1,ndo amux=amux-pix(imo,imo)*2.0d0 amuy=amuy-piy(imo,imo)*2.0d0 111 amuz=amuz-piz(imo,imo)*2.0d0 ax=z0 ay=z0 az=z0 do 2 i=1,noa ax=ax+xv(1,i)*z(i) ay=ay+xv(2,i)*z(i) 2 az=az+xv(3,i)*z(i) dx=amux+ax dy=amuy+ay dz=amuz+az dd=sqrt(dx*dx+dy*dy+dz*dz) write(3,*) write(6,*)' Dipole moment in atomic units:' write(6,4000)amux,amuy,amuz,ax,ay,az,dx,dy,dz,dd write(3,*)' Dipole moment in atomic units:' write(3,4000)amux,amuy,amuz,ax,ay,az,dx,dy,dz,dd 4000 format(/,' ',15x,'X',15x,'Y',15x,'Z',/, 1 ' electronic ',3F16.8,/, 1 ' nuclear ',3F16.8,/, 2 ' ',48(1H-),/, 3 ' total ',3F16.8,/, 4 ' |u| ', F16.8,/) write(3,*)' Dipole moment in Debyes:' write(6,*)' Dipole moment in Debyes:' debye=2.541765d0 write(3,4000)amux*debye,amuy*debye,amuz*debye, 1ax*debye,ay*debye,az*debye,dx*debye,dy*debye,dz*debye,dd*debye write(6,4000)amux*debye,amuy*debye,amuz*debye, 1ax*debye,ay*debye,az*debye,dx*debye,dy*debye,dz*debye,dd*debye c 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=800) 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, 1a10 common/atoms/z(nat0),zsum,xv(3,nat0),multi logical lopt,lnno common/opts/lopt(300),icha,npl(10) 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(3,4000)qe,qn,dt 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 overlap6 c 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 none real*8 q,zsum,xv,cc,alpha,rcp,cij,bij,e,be,pi,bohr,debye, 1z0,one,two,three,four,five,half,a10,b6,sv,aai,aaj,aad,aaf,aa1, 1cci,ccj,ai,aiii,aij,aj,ajjj,an2x,an2y,an2z,anorm,anx,any,anz, 2coe,con,dx,dxx,dy,dyy,dz,dzz,ec,ef,oi,osxm,osxp,osym,osyp, 3oszp,oszm,ox,oxdx,oxx,oy,oydy,oyy,oz,ozdz,ozz,r2ij, 4sx,sxm,sxmm,sy,sym,symm,sz,szm,szmm,szp,szpp,sxp,sxpp,syp,sypp, 5t,tol,tt,x2ij,yi,z2ij,zi,zij,zj,anjx,anjy,anjz,spi,xi,xij, 6y2ij,yij,yj,xj integer*4 nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,nori, 1nop,multi,nxi,nyi,nzi,nxj,nyj,nzj,ione,itwo,i,i0,ia,icha,npl, 2ioj,ioo,iorb,ip,ipp,ish,j,ja,jo,joo,jp,k,kk,n1e,nat0, 3nor,nort,npi,npii,npj,npjj,ii,io,jj,jpp,noa character*2 at character*1 ty parameter (nat0=96,nor=80,nort=800,n1e=17) integer*4 u0 parameter (u0=28) c n1e ... number of kinds of one-electron integrals 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/const/pi,bohr,debye common/numbers/z0,one,two,three,four,five,half, 1a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(n1e,u0,nort),b6(n1e,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,npl(10) 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 ) c call inipoly i0=0 ione=1 itwo=2 call initfd(aa1,aad,aaf) c tol=100.0d0 spi=sqrt(pi*pi*pi) do 10 i=1,nao do 10 j=1,u0 do 10 k=1,n1e 10 sv(k,j,i)=z0 write(3,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C open(33,file='SAO.SCR',form='unformatted') open(34,file='P.SCR',form='unformatted') open(35,file='Q.SCR',form='unformatted') open(36,file='V.SCR',form='unformatted') open(37,file='L.SCR',form='unformatted') open(38,file='L2.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,lg9) 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,lg9) 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) anjx=dble(nxj(joo)) anjy=dble(nyj(joo)) anjz=dble(nzj(joo)) anx=anjx*(anjx-1.0d0) any=anjy*(anjy-1.0d0) anz=anjz*(anjz-1.0d0) an2x=2.0d0*(2.0d0*anjx+1.0d0) an2y=2.0d0*(2.0d0*anjy+1.0d0) an2z=2.0d0*(2.0d0*anjz+1.0d0) 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) sxm=z0 sym=z0 szm=z0 if(nxj(joo).gt.0) 1sxm=oi(i0,nxi(ioo),nxj(joo)-1,ai,aj,xij,t,xi,xj) sxp=oi(i0,nxi(ioo),nxj(joo)+1,ai,aj,xij,t,xi,xj) if(nyj(joo).gt.0) 1sym=oi(i0,nyi(ioo),nyj(joo)-1,ai,aj,yij,t,yi,yj) syp=oi(i0,nyi(ioo),nyj(joo)+1,ai,aj,yij,t,yi,yj) if(nzj(joo).gt.0) 1szm=oi(i0,nzi(ioo),nzj(joo)-1,ai,aj,zij,t,zi,zj) szp=oi(i0,nzi(ioo),nzj(joo)+1,ai,aj,zij,t,zi,zj) ox =oi(ione,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) oy =oi(ione,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) oz =oi(ione,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) oxx=oi(itwo,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) oyy=oi(itwo,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) ozz=oi(itwo,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) dx=anjx*sxm-two*aj*sxp dy=anjy*sym-two*aj*syp dz=anjz*szm-two*aj*szp osxm=0.0d0 osym=0.0d0 oszm=0.0d0 if(nxj(joo).gt.0) 1osxm=oi(ione,nxi(ioo),nxj(joo)-1,ai,aj,xij,t,xi,xj) if(nyj(joo).gt.0) 1osym=oi(ione,nyi(ioo),nyj(joo)-1,ai,aj,yij,t,yi,yj) if(nzj(joo).gt.0) 1oszm=oi(ione,nzi(ioo),nzj(joo)-1,ai,aj,zij,t,zi,zj) osxp=oi(ione,nxi(ioo),nxj(joo)+1,ai,aj,xij,t,xi,xj) osyp=oi(ione,nyi(ioo),nyj(joo)+1,ai,aj,yij,t,yi,yj) oszp=oi(ione,nzi(ioo),nzj(joo)+1,ai,aj,zij,t,zi,zj) oxdx=anjx*osxm-two*aj*osxp oydy=anjy*osym-two*aj*osyp ozdz=anjz*oszm-two*aj*oszp sxmm=0.0d0 symm=0.0d0 szmm=0.0d0 if(nxj(joo).gt.1) 1sxmm=oi(i0,nxi(ioo),nxj(joo)-2,ai,aj,xij,t,xi,xj) if(nyj(joo).gt.1) 1symm=oi(i0,nyi(ioo),nyj(joo)-2,ai,aj,yij,t,yi,yj) if(nzj(joo).gt.1) 1szmm=oi(i0,nzi(ioo),nzj(joo)-2,ai,aj,zij,t,zi,zj) sxpp=oi(i0,nxi(ioo),nxj(joo)+2,ai,aj,xij,t,xi,xj) sypp=oi(i0,nyi(ioo),nyj(joo)+2,ai,aj,yij,t,yi,yj) szpp=oi(i0,nzi(ioo),nzj(joo)+2,ai,aj,zij,t,zi,zj) dxx=anx*sxmm-an2x*aj*sx+4.0d0*aj*aj*sxpp dyy=any*symm-an2y*aj*sy+4.0d0*aj*aj*sypp dzz=anz*szmm-an2z*aj*sz+4.0d0*aj*aj*szpp c overlap b6(1,ioo,joo)=sx*sy*sz c electric dipole b6(2,ioo,joo)=ox*sy*sz b6(3,ioo,joo)=sx*oy*sz b6(4,ioo,joo)=sx*sy*oz c electric quadrupole xx,xy,xz,yy,yz,zz: b6(5,ioo,joo)=oxx*sy*sz b6(6,ioo,joo)=ox*oy*sz b6(7,ioo,joo)=ox*sy*oz b6(8,ioo,joo)=oyy*sz*sx b6(9,ioo,joo)=sx*oy*oz b6(10,ioo,joo)=ozz*sx*sy c gradient/velocity d/dx, d/dy, d/dz b6(11,ioo,joo)=dx*sy*sz b6(12,ioo,joo)=dy*sz*sx b6(13,ioo,joo)=dz*sx*sy c angular momentum lx,ly,lz b6(14,ioo,joo)=(oy*dz-oz*dy)*sx*half b6(15,ioo,joo)=(oz*dx-ox*dz)*sy*half b6(16,ioo,joo)=(ox*dy-oy*dx)*sz*half c square of angular momentum b6(17,ioo,joo)= 1 -oxx*dyy*sz -oxx*dzz*sy 1-oyy*dxx*sz -oyy*dzz*sx 1-ozz*dxx*sy -ozz*dyy*sx 2+2.0d0*(oxdx*sy*sz+oydy*sx*sz+ozdz*sx*sy) 1 +oxdx*oydy*sz +oxdx*ozdz*sy 1+oydy*oxdx*sz +oydy*ozdz*sx 1+ozdz*oxdx*sy +ozdz*oydy*sx do 71 k=1,n1e 71 b6(k,ioo,joo)=b6(k,ioo,joo)*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 do 74 k=1,n1e 74 sv(k,ioo,ioj)=sv(k,ioo,ioj)+ajjj*b6(k,ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo do 7 k=1,n1e 7 sv(k,ioo,ioj)=sv(k,ioo,ioj)+b6(k,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 write(33)(sv(1,ii,jj),jj=1,iorb) c c check the normalization condition - orbital iorb anorm=one/sqrt(sv(1,ii,iorb)) if(abs(anorm-one).gt.0.001D0)then write(3,*)anorm,iorb if(.not.lnno)call report('Non-normalized aos') endif if(lsprint)write(3,*)iorb,anorm do 331 kk=2,4 331 write(34)(sv(kk,ii,jj),jj=1,nao) do 332 kk=5,10 332 write(35)(sv(kk,ii,jj),jj=1,nao) do 333 kk=11,13 333 write(36)(sv(kk,ii,jj),jj=1,nao) do 334 kk=14,16 334 write(37)(sv(kk,ii,jj),jj=1,nao) write(38)(sv(17,ii,jj),jj=1,nao) write(6,*)io+ii,io+ii,sv(17,ii,io+ii) 33 continue io=io+npii c c zero-out temporary buffer: do 9 ii=1,u0 do 9 jj=1,nao do 9 kk=1,n1e 9 sv(kk,ii,jj)=z0 endif 3 continue 2 continue close(33) close(34) close(35) close(36) close(37) close(38) c c write(6,*)' One electron integrals 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=800) 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,npl(10) 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 trfo(qxx,nmo,nao,cij) c trafo to MO implicit none integer*4 nmo,nao,nort,imo,j,jmo,i parameter (nort=800) real*8 qxx(nort,nort),cij(nort,nort),tem(nort,nort),sum do 2 imo=1,nmo do 2 j=1,nao sum=0.0d0 do 3 i=1,nao 3 sum=sum+cij(imo,i)*qxx(i,j) 2 tem(imo,j)=sum do 4 imo=1,nmo do 4 jmo=1,nmo sum=0.0d0 do 5 j=1,nao 5 sum=sum+cij(jmo,j)*tem(imo,j) 4 qxx(imo,jmo)=sum return end subroutine thetarr implicit real*8 (a-h,o-z) implicit integer*4 (i-n) integer*4 u0,nb0 parameter (nb0=200,nat0=96,nor=80,nort=800,u0=28) 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) real*8 l2,aad(u0,u0),aaf(u0,u0),aa1(u0,u0) common/qq/qxx(nort,nort),qyy(nort,nort),qzz(nort,nort), 1qxy(nort,nort),qxz(nort,nort),qyz(nort,nort),l2(nort,nort) common/qb/bxx(nort,nort),byy(nort,nort),bzz(nort,nort), 1bxy(nort,nort),bxz(nort,nort),byz(nort,nort),b2(nort,nort) logical lconfig,lopt equivalence (lconfig,lopt(65)) common/opts/lopt(300),icha,npl(10) character*10 oname real*8 cut(3,nb0,nort) common/cuts/cut c bohr=0.529177d0 call initfd(aa1,aad,aaf) c do 14 is=1,2 if(.not.lconfig.and.is.eq.2)goto 14 c if(is.eq.1)then c read AO integrals, skip off-diagonals open(33,file='Q.SCR',form='unformatted') open(38,file='L2.SCR',form='unformatted') do 1 i=1,nao read(33)(qxx(i,j),j=1,nao) read(33)(qxy(i,j),j=1,nao) read(33)(qxz(i,j),j=1,nao) read(33)(qyy(i,j),j=1,nao) read(33)(qyz(i,j),j=1,nao) read(33)(qzz(i,j),j=1,nao) 1 read(38)(l2(i,j) ,j=1,nao) close(33) close(38) write(6,*)'Transforming alpha-orbitals' call trfo(qxx,nmo,nao,cij) call trfo(qxy,nmo,nao,cij) call trfo(qxz,nmo,nao,cij) call trfo(qyy,nmo,nao,cij) call trfo(qyz,nmo,nao,cij) call trfo(qzz,nmo,nao,cij) call trfo(l2 ,nmo,nao,cij) else open(33,file='Q.SCR',form='unformatted') open(38,file='L2.SCR',form='unformatted') do 6 i=1,nao read(33)(bxx(i,j),j=1,nao) read(33)(bxy(i,j),j=1,nao) read(33)(bxz(i,j),j=1,nao) read(33)(byy(i,j),j=1,nao) read(33)(byz(i,j),j=1,nao) read(33)(bzz(i,j),j=1,nao) 6 read(38)(b2(i,j) ,j=1,nao) close(33) close(38) write(6,*)'Transforming beta-orbitals' call trfo(bxx,nmo,nao,bij) call trfo(bxy,nmo,nao,bij) call trfo(bxz,nmo,nao,bij) call trfo(byy,nmo,nao,bij) call trfo(byz,nmo,nao,bij) call trfo(bzz,nmo,nao,bij) call trfo(b2 ,nmo,nao,bij) endif write(6,6001) write(3,6001) 6001 format(' Orbital (au) sqrt()(au) sqrt()(A) xx yy zz l',/, 1 ' ---------------------------------------------------') do 14 imo=1,nmo if(is.eq.1)then rxx=qxx(imo,imo) ryy=qyy(imo,imo) rzz=qzz(imo,imo) a2=l2(imo,imo) else rxx=bxx(imo,imo) ryy=byy(imo,imo) rzz=bzz(imo,imo) a2=b2(imo,imo) endif r2=rxx+ryy+rzz a2=(dsqrt(1.0d0+4.0d0*a2)-1.0d0)/2.0d0 r=dsqrt(r2) c c investigate cuts to find nodes: nx=0 ny=0 nz=0 z0=0.0d0 do 312 i=1,npl(3)+1 do 312 io=1,3 312 cut(io,i,imo)=z0 write(oname,2006)imo 2006 format(i10) do 2 io=1,3 if(imo.ge.npl(1).and.imo.le.npl(2))then if(is.eq.1)then if(io.eq.1)open(45,file=oname//'x.prn') if(io.eq.2)open(45,file=oname//'y.prn') if(io.eq.3)open(45,file=oname//'z.prn') if(io.eq.1)open(46,file=oname//'dx.prn') if(io.eq.2)open(46,file=oname//'dy.prn') if(io.eq.3)open(46,file=oname//'dz.prn') else if(io.eq.1)open(45,file=oname//'bx.prn') if(io.eq.2)open(45,file=oname//'by.prn') if(io.eq.3)open(45,file=oname//'bz.prn') if(io.eq.1)open(46,file=oname//'bdx.prn') if(io.eq.2)open(46,file=oname//'bdy.prn') if(io.eq.3)open(46,file=oname//'bdz.prn') endif endif if(npl(3)+1.gt.nb0)call report('too many plotting points') do 311 i=1,npl(3)+1 t=(-1.0d0+dble(2*(i-1))/dble(npl(3)))*r*dble(npl(4)) if(io.eq.1)cut(1,i,imo)=forbxp4(t,z0,z0,imo,1,aa1,aad,aaf) if(io.eq.2)cut(2,i,imo)=forbxp4(z0,t,z0,imo,1,aa1,aad,aaf) 311 if(io.eq.3)cut(3,i,imo)=forbxp4(z0,z0,t,imo,1,aa1,aad,aaf) do 31 i=1,npl(3)+1 if(i.lt.npl(3))then f =cut(io,i ,imo) fp =cut(io,i+1,imo) fpp=cut(io,i+2,imo) if(((f.lt.z0.and.fp.ge.z0).or.(f.gt.z0.and.fp.le.z0)). 1 and.(fp.ne.z0.or.fpp.ne.z0))then if(io.eq.1)nx=nx+1 if(io.eq.2)ny=ny+1 if(io.eq.3)nz=nz+1 endif endif if(imo.ge.npl(1).and.imo.le.npl(2))write(46,4545)t*bohr,f**2 31 if(imo.ge.npl(1).and.imo.le.npl(2))write(45,4545)t*bohr,f 4545 format(f7.2,f15.10) if(imo.ge.npl(1).and.imo.le.npl(2))close(46) 2 if(imo.ge.npl(1).and.imo.le.npl(2))close(45) write(3,6002)imo,r2,r,r*bohr,rxx,ryy,rzz,a2,nx,ny,nz 6002 format(i5,6f9.1,f5.1,1x,3i3) write(6,6002)imo,r2,r,r*bohr,rxx,ryy,rzz,a2,nx,ny,nz 14 continue 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=800) 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, 1a10 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,npl(10) 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 nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip,lg9) 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,lg9) 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(6,*)'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,lg9) 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, 1a10 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, 1a10 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 subroutine nppp(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip,lg9) implicit none integer*4 u0,nat0,nor,npi,npii,nxi,nyi,nzi,icha,npl, 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,lopt,lg94,lg9 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) common/opts/lopt(300),icha,npl(10) equivalence (lopt(27),lg94) 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 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 if(lg9)then npii=9 npi=9 endif 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 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, 1a10 dimension aa1(u0,u0),aad(u0,u0),aaf(u0,u0) 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 c -------------------------------- c 3z^2-R^2 1 xx c xz 2 yy c yz 3 zz c x^2-y^2 4 xy c xy 5 xz 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 aad(4,1)=d5 aad(4,2)=-d5 aad(5,4)=one 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 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, 1a10 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 =============================================== function forbxp4(rx,ry,rz,i,ispin,aa1,aad,aaf) c c orbital i at the point rx,ry,rz,without initfd c implicit none integer*4 nat0,nor,nort,u0,multi,icha,npl parameter (nat0=96,nor=80,nort=800,u0=28) character*2 at character*1 ty real*8 q,zsum,xv common/atoms/q(nat0),zsum,xv(3,nat0),multi logical ld5,lf7,lopt,special,lg9 common/opts/lopt(300),icha,npl(10) real*8 cc,alpha,rcp,cij,bij,e,be integer*4 nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,nori, 1nop,ish,noa 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 z0,one,two,three,four,five,half,a10 common/numbers/z0,one,two,three,four,five,half, 1a10 equivalence (lopt(4),ld5),(lopt(8),lf7),(lopt(66),lg9 ) real*8 aai(u0,u0),cci(u0),axs(4),ays(4),azs(4), 1aa1(u0,u0),aad(u0,u0),aaf(u0,u0) integer*4 nxi(u0),nyi(u0),nzi(u0),io,ia,ip,ipx,npi,npii,i,ioo, 1ispin real*8 psii,forbxp4,x,y,z,ar2,a,ecp,ef,rx,ry,rz,b6 c psii=z0 c c MO psi(i)=sum over AOs io=0 do 2 ia=1,noa x=rx-xv(1,ia) y=ry-xv(2,ia) z=rz-xv(3,ia) ar2=x*x+y*y+z*z axs(1)=one axs(2)=x axs(3)=x*x axs(4)=x*x*x ays(1)=one ays(2)=y ays(3)=y*y ays(4)=y*y*y azs(1)=one azs(2)=z azs(3)=z*z azs(4)=z*z*z c c sum over primitives do 3 ip=1,nop(ia) a=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip,lg9) ecp=a*ar2 if(ecp.gt.12.0D0)goto 3 ef=exp(-ecp) special=npi.ne.npii c do 4 ipx=1,npi b6=cci(ipx)*ef*axs(nxi(ipx)+1)*ays(nyi(ipx)+1)*azs(nzi(ipx)+1) if(ispin.eq.1)then if(special)then do 74 ioo=1,npii 74 psii=psii+aai(ioo,ipx)*b6*cij(i,io+ioo) else psii=psii+b6*cij(i,io+ipx) endif else if(special)then do 741 ioo=1,npii 741 psii=psii+aai(ioo,ipx)*b6*bij(i,io+ioo) else psii=psii+b6*bij(i,io+ipx) endif endif 4 continue c 3 if(ish(ia,ip+1).gt.ish(ia,ip))io=io+npii 2 continue forbxp4=psii return end