PROGRAM fof IMPLICIT none real*8 alpha write(6,600) 600 format(' Decomposition of Gaussian spherical harmonics',/, 1 ' into MOs',/,/,' Input alpha: ',$) read(5,*)alpha open(3,file='FOF.OUT') call init CALL IOPAR CALL ROAAIc call whatdiff(3,alpha) close(3) end subroutine whatdiff(l,a0) implicit none real*8 a0,tr,w1,wJ integer*4 l,m,MCart,i,j1,J,im,imend,MPure,i1,i2 real*8,allocatable::Coef2P(:,:),Coef2C(:,:),cao(:,:) real*8 fac,f0 common/ffac/f0,fac(6) integer*4 nat0,nor,nort parameter (nat0=190,nor=80,nort=700) real*8 cc,alpha,rcp,cij,bij,e,be integer*4 nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,noa,nori,nop, 1ish character*1 ty,ok character*2 at common/work/cc(nat0,nor),alpha(nat0,nor),rcp(nat0,nor), 1cij(nort,nort),bij(nort,nort),e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori,nop(nat0),ty(nat0,nor), 1ish(nat0,nor),at(nat0),ok integer*4 ii,ia,ila,iz real*8 z,zsum,xv common/atoms/z(nat0),zsum,xv(3,nat0),iz(nat0) real*8,allocatable::t(:,:),smoa(:,:),smob(:,:),smoar(:,:), 1smobr(:,:),smoai(:,:),smobi(:,:),list(:),listb(:) logical lopt,lsprint integer*4 icha common/opts/lopt(300),icha equivalence (lopt(1),lsprint) integer*4 n10 parameter (n10=10) real*8 cmax(n10) integer*4 imax(n10) character*1 ab(n10) if(l.gt.6)call report('maximum l overflow') c find lanthanide ila=0 ii=0 do 1 ia=1,nat if(iz(ia).ge.57.and.iz(ia).le.70)then ii=ii+1 ila=ia endif 1 continue if(ii.eq.0)call report('no lanthanide found') if(ii.gt.1)then write(6,*)(iz(ia),ia=1,nat) call report('more than one lanthanide found') endif i=0 fac(i) = 1.0d0 Do 10 i = 1, (2*6) 10 fac(i) = dble(i)*fac(i-1) c number of cartesian functions: MCart = ((l+1)*(l+2))/2 c number of spherical harmonics: MPure=2*l+1 c calculate the transformation coefficients pure-cart: allocate(Coef2P(MCart,MPure),Coef2C(MCart,MPure)) Call C2PGn1(L,MCart,MPure,Coef2P,Coef2C) c calculate overlap integrals > allocate(cao(MCart,nao)) call ccao(cao,MCart,nao,l,ila,a0) c calculate overlap integrals > allocate(t(Mpure,nao),list(nao)) if(lsprint) write(3,*)'Spherical overlaps:' ia=0 do 2 m =0,l if(m.eq.0)then imend=1 else imend=2 endif do 2 im=1,imend c im=1 (|lm>+|l-m>/sqrt( 2) c im=2 (|lm>-|l-m>/sqrt(-2) ia=ia+1 if(lsprint)then if(m.eq.0)then write(3,6002)l,m 6002 format(' (<',i1,i3,'|:') else if(im.eq.1)then write(3,6001)l,m,'+',l,-m 6001 format(' (<',i1,i3,'| ',a1,' <',i1,i3,'|)/sqrt(2):') else write(3,6001)l,m,'-',l,-m endif endif endif do 21 j=1,nao t(ia,j)=0.0d0 do 22 ii=1,MCart 22 t(ia,j)=t(ia,j)+Coef2P(ii,ia)*cao(ii,j) 21 list(j)=t(ia,j) if(lsprint) call nlist(3,list,nao) 2 continue c calculate overlap integrals > allocate(smoa(Mpure,nao),smob(Mpure,nao),listb(nao)) if(lsprint)write(3,*)'Spherical-real|MO>:' do 3 ia=1,MPure do 31 J=1,nao smoa(ia,J)=0.0d0 smob(ia,J)=0.0d0 do 311 i=1,nao smoa(ia,J)=smoa(ia,J)+t(ia,i)*cij(J,i) 311 smob(ia,J)=smob(ia,J)+t(ia,i)*bij(J,i) list( J)=smoa(ia,J) 31 listb(J)=smob(ia,J) if(lsprint)write(3,*)ia,' alpha:' if(lsprint) call nlist(3,list ,nsa) if(lsprint)write(3,*)ia,' beta :' 3 if(lsprint) call nlist(3,listb,nsb) allocate(smoar(Mpure,nao),smobr(Mpure,nao)) allocate(smoai(Mpure,nao),smobi(Mpure,nao)) tr=dsqrt(2.0d0) ia=0 do 4 m =0,l if(m.eq.0)then ia=ia+1 do 41 J=1,nao smoar(m+l+1,J)=smoa(ia,J) smobr(m+l+1,J)=smob(ia,J) smoai(m+l+1,J)=0.0d0 41 smobi(m+l+1,J)=0.0d0 else ia=ia+1 do 42 J=1,nao c |l m> smoar( m+l+1,J)= smoa(ia ,J)/tr smobr( m+l+1,J)= smob(ia ,J)/tr smoai( m+l+1,J)= smoa(ia+1,J)/tr smobi( m+l+1,J)= smob(ia+1,J)/tr c |l-m> smoar(-m+l+1,J)= smoa(ia ,J)/tr smobr(-m+l+1,J)= smob(ia ,J)/tr smoai(-m+l+1,J)=-smoa(ia+1,J)/tr 42 smobi(-m+l+1,J)=-smob(ia+1,J)/tr ia=ia+1 endif 4 continue if(lsprint)then write(3,*)'Spherical|alpha-MO>:' do 6 m=-l,l write(3,3001)l,m,'real' 3001 format(' <',i1,i3,'|',a4,'>:') do 61 J=1,nsa 61 list(J)=smoar(m+l+1,J) call nlist(3,list,nsa) write(3,3001)l,m,'imag' do 62 J=1,nsa 62 list(J)=smoai(m+l+1,J) 6 call nlist(3,list,nsa) write(3,*)'Spherical|beta-MO>:' do 7 m=-l,l write(3,3001)l,m,'real' do 71 J=1,nsb 71 list(J)=smobi(m+l+1,J) call nlist(3,list,nsb) write(3,3001)l,m,'imag' do 72 J=1,nsb 72 list(J)=smobi(m+l+1,J) 7 call nlist(3,list,nsb) endif do 5 m =-l,l do 500 i1=1,n10 imax(i1)=0 cmax(i1)=0.0d0 ab(i1)='X' do 51 J=1,nsa do 53 i2=1,i1-1 53 if(imax(i2).eq.J.and.ab(i2).eq.'A')goto 51 wJ=smoar(m+l+1,J)**2+smoai(m+l+1,J)**2 if(wJ.gt.cmax(i1))then cmax(i1)=wJ imax(i1)=J ab(i1)='A' endif 51 continue do 52 J=1,nsb do 54 i2=1,i1-1 54 if(imax(i2).eq.J.and.ab(i2).eq.'B')goto 52 wJ=smobr(m+l+1,J)**2+smobi(m+l+1,J)**2 if(wJ.gt.cmax(i1))then cmax(i1)=wJ imax(i1)=J ab(i1)='B' endif 52 continue 500 continue 5 write(6,609)l,m,imax(1),ab(1),cmax(1)*100.0d0, 1imax(2),ab(2),cmax(2)*100.0d0,imax(3),ab(3),cmax(3)*100.0d0 609 format(' |LM> = |',i1,i3,'> J_MO = ',3(i5,a1,f10.2,' %')) return end SUBROUTINE ROAAIc IMPLICIT none integer*4 NATF,nmos nmos=0 call ROAAI(NATF,nmos) return end c SUBROUTINE ROAAI(NATF,nmos) c c ab initio (SOS) calculations c IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) parameter (nat0=190,nor=80,nort=700) c C nat0: # of atoms, nor: # of orbitals on one atom, C ndf0: # of basis functions, c character*2 at character*1 ok,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),ok logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,lzmat,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94,lnmr,lapt,lader,lcnmr, 9lopt,lhell,lsuper,lgrad,lssup,lsupinv,lplot,lcroa,lgvir, 1lci,li3,lconfig,lroax,lnormao,lfuzzy,lgradc,laolo,ldipole, 1laol3,lcaz,lbasred,lcazdip,lord,lboys,lsrx,lorban 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(28),lzmat ),(lopt(29),lnmr ),(lopt(30),lapt ), 1(lopt(31),lader ),(lopt(41),lhell ),(lopt(43),lsuper), 1(lopt(45),lgrad ),(lopt(46),lssup ),(lopt(47),lsupinv), 1(lopt(48),lplot ),(lopt(52),lcroa ),(lopt(54),lgvir ), 1(lopt(56),lcnmr ),(lopt(59),lci ),(lopt(54),li3 ), 1(lopt(83),laol3 ),(lopt(84),lcaz ),(lopt(92),lbasred ), 1(lopt(65),lconfig),(lopt(76),lnormao),(lopt(73),lroax ), 1(lopt(81),laolo ),(lopt(82),ldipole),(lopt(77),lfuzzy ), 1(lopt(78),lgradc ),(lopt(94),lcazdip),(lopt(97),lord), 1(lopt(103),lboys),(lopt(104),lsrx),(lopt(106),lorban) common/rundiff/idif,igeo,ibasis,ieigen character*80 ifn common/filenames/ifn c C idif=0 nmo=nmos do 1 i=1,nat0 do 1 j=1,nor 1 ish(i,j)=32000 c open(2,file=ifn) c c read in basic geometry, basis set, and wavefunction NDIFF=0 ieigen=0 igeo=0 ibasis=0 call roainp(NDIFF) natf=nat write(6,*)nat c c control output call writecon(NDIFF) c if(lnormao)then call normorb write(3,*)'AOs have been renormalized' write(6,*)'AOs have been renormalized' endif c c one electron integrals for starting point call overlap6 return END subroutine roainp(ndiff) c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c parameter (nat0=190,nor=80,nort=700,MENDELEV=89) character*2 at,s2,symbols(MENDELEV) character*1 ok,ty common/const/pi,bohr 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),ok 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/rundiff/idif,igeo,ibasis,ieigen common/atoms/z(nat0),zsum,xv(3,nat0),iz(nat0) logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,lzmat,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lconfig,lturbo,lbasred 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(28),lzmat ),(lopt(65),lconfig),(lopt(68),lturbo), 1(lopt(92),lbasred ) common/iopts/iopt(300) equivalence (iopt(1),ncnmr),(iopt(2),inormo ),(iopt(3),nek1 ), 1 (iopt(4),nek2 ),(iopt(5),nej1 ),(iopt(6),nej2 ), 1 (iopt(7),ioe ),(iopt(8),nstart ),(iopt(9),nend ), 1 (iopt(10),nprc),(iopt(11),multi ),(iopt(12),ndog ), 1 (iopt(13),nsi ),(iopt(14),nconfig),(iopt(15),nproc ), 1 (iopt(16),neci),(iopt(17),idiago ),(iopt(18),igradc), 1 (iopt(19),nev ),(iopt(20),mmit ),(iopt(21),iham ), 1 (iopt(22),ncut),(iopt(23),nk ),(iopt(24),nj ) character*80 key integer*4 ipse(nat0) data symbols/'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne', 3'Na','Mg','Al','Si','P ','S ','Cl','Ar', 4'K ','Ca','Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 5'Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', 5 'In','Sn','Sb','Te','I ','Xe', 6'Cs','Ba','La', 6 'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho', 6 'Er','Tm','Yb','Lu', 6'Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg', 6 'Tl','Pb','Bi','Po','At','Rn', 7'Fr','Ra','Ac'/ 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) integer*4 n10 parameter (n10=10) real*8 cmax(n10) integer*4 imax(n10) c do i=1,nort rcl(i)=0.0d0 enddo nst=0 npps=0 nstam=0 npe=0 ncc=0 if(lturbo)then call readtm return endif c c Browse at the G94 or C5 output: ic=0 111 ic=ic+1 read(2,1001,end=666)key 1001 format(a80) c { 1234567890123456789012345678901234567890} C c c GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG c GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG c GAUSSIAN 94 OUTPUT CATCHING: if(LG94)then c EIGENVALUES (PRINTED WITH THE OPTION IOP(6/7=1)) if(key(1:35).eq.' Molecular Orbital Coefficients'.or. 1 key(2:11).eq.'Alpha MOs:'.and..not.lex) 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(lex)goto 877 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) 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 ndog=ndo nsi=nsa-nsb 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(lbasred)then write(6,*)'Basis check only' return endif 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 ndog=ndo nsi=nsa-nsb ichm=icha mtp=1 endif c c BASIS GFPrint c if you want to use this, uncomment next and comment following line: c if(ibasis.eq.0)THEN if(1.eq.0)THEN if((key(2:16).eq.'Standard basis:'.or.key(6:14).eq.'Exponent='. 1 or.key(2:14).eq.'General basis').and.ndiff.eq.0)then if(key(6:14).eq.'Exponent=')then 2222 read(2,1001)key if(key(2:10).ne.'*********')goto 2222 backspace 2 else write(6,*)key endif c c check whether general basis follows read(2,1001)key 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 backspace 2 goto 111 endif c ibasis=ibasis+1 read(2,3336)s2 3336 format(1x,a2) c if(s2.eq.'**')then do 4 i=1,5 4 read(2,*) nao=0 c naolast=0 do 7 i=1,nat nop(i)=0 do 6 k=1,nor 6 ish(i,k)=32000 ishell=0 read(2,*) 3906 read(2,3901,err=7,iostat=ioss)nao,sf 3901 format(47x,I3,10x,F9.2) if(ioss.ne.0)goto 7 if(lopt(298))write(3,3901)nao,sf ishell=ishell+1 ndao=nao-naolast naolast=nao 3904 nop(i)=nop(i)+1 k=nop(i) if(k.gt.nor)call report('Too many primitives on one atom') read(2,3336)s2 if(s2.eq.'**'.or.s2.eq.'*-')goto 3903 backspace 2 read(2,3902,err=3903,iostat=ioss)af,ccs(i,k),ccp(i,k),cc(i,k) 1 ,dum 3902 format(71x,5g12.6) if(abs(af).lt.0.000001d0)goto 3903 if(ioss.ne.0)goto 3903 if(lopt(298))write(3,3902)af,ccs(i,k),ccp(i,k),cc(i,k),dum af=af*sf**2 alpha(i,k)=af ish(i,k)=ishell c 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 ccd(i,k,1)=cc(i,k)*SNDD ccd(i,k,2)=cc(i,k)*SND c ccf(i,k,1)=dum*SNF1 ccf(i,k,2)=dum*SNF2 ccf(i,k,3)=dum*SNF3 c c store exp coef. for control listing: cc(i,k)=ccs(i,k) rcp(i,k)=ccp(i,k) c ccs(i,k)=ccs(i,k)*SNO ccp(i,k)=ccp(i,k)*SNP if(ndao.eq.1)ty(i,k)='S' if(ndao.eq.3)ty(i,k)='P' if(ndao.eq.4)ty(i,k)='L' if(ndao.eq.5)ty(i,k)='D' if(ndao.eq.5)ld5=.true. if(ndao.eq.6)ty(i,k)='D' if(ndao.eq.7)ty(i,k)='F' if(ndao.eq.7)lf7=.true. if(ndao.eq.10)ty(i,k)='F' 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 if(ty(i,k).eq.'F')cc(i,k)=dum goto 3904 3903 nop(i)=nop(i)-1 backspace 2 read(2,3336)s2 if(s2.eq.'*-'.or.s2.eq.'**')goto 7 backspace 2 goto 3906 7 continue write(6,*) 'G94 standard basis loaded from output' write(6,*) nao,' AOs' if(nao.gt.nort)call report('too many orbitals') c c else c call report('set the iop(5/33=1) for the calculation !') call report('set the GFPRINT option for the calculation !') c endif c endif ENDIF c noa=nat c endif c (end of Gaussian output catching) c c GAUSSIAN FORMATTED CHECKPOINT: c FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF if(key(1:22).eq.'Alpha Orbital Energies')then write(6,1001)key if(nmo.eq.0)nmo=nao write(6,*)nao,nmo read(2,2009)(e(i1),i1=1,nmo) write(6,*)nao,nmo do 10 i1=1,nmo 10 be(i1)=e(i1) write(6,*)nao,nmo endif if(key(1:21).eq.'Alpha MO coefficients')then write(6,1001)key write(3,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) 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 write(3,1001)key read(2,2009)((bij(i1,i2),i2=1,nao),i1=1,nmo) 99 write(6,*)'Jakej AO chces vedet? ' read(5,*)j if(j.gt.0)then do 910 i1=1,n10 cmax(i1)=0.0d0 imax(i1)=0 do 910 i=1,nsa do 901 i2=1,i1-1 901 if(imax(i2).eq.i)goto 910 if(cij(i,j)**2.gt.cmax(i1))then cmax(i1)=cij(i,j)**2 imax(i1)=i endif 910 continue write(6,*)'Occupied alpha MOS:' write(6,609)imax(1),cij(imax(1),j), 1 imax(2),cij(imax(2),j),imax(3),cij(imax(3),j) 609 format(3(i4,f10.4)) do 810 i1=1,n10 cmax(i1)=0.0d0 imax(i1)=0 do 810 i=1,nsb do 801 i2=1,i1-1 801 if(imax(i2).eq.i)goto 810 if(bij(i,j)**2.gt.cmax(i1))then cmax(i1)=bij(i,j)**2 imax(i1)=i endif 810 continue write(6,*)'Occupied beta MOS:' write(6,609)imax(1),bij(imax(1),j), 1 imax(2),bij(imax(2),j),imax(3),bij(imax(3),j) goto 99 endif 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:14).eq.'Atomic numbers')then read(2,2109)(iz(i),i=1,noa) 2109 format(6i12) do 91 i=1,noa zsum=zsum+z(i) 91 at(i)=symbols(iz(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 ndog=ndo nsi=nsa-nsb 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 if(ibasis.gt.0)idif=ieigen-1 c if(idif.eq.ndiff.and.ibasis.gt.0)return goto 111 666 close(2) write(3,4004)idif,ieigen,igeo,ibasis write(6,4004)idif,ieigen,igeo,ibasis 4004 format(I4,' point of the differentiation:',/, 1 '----', ' -----------------------------',/, 2 ' Eigenvectors read',I4,' times',/, 3 ' Geometry read ',I4,' times',/, 4 ' Basis read ',I4,' times',/) if(igeo.eq.0)call report('geometry not found') if(ibasis.eq.0)call report('basis not found') if(ieigen.eq.0)call report('eigenvectors not found') return 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 common/const/pi,bohr 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 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,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 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 ')then if(ld5)then ndao=5 else ndao=6 endif endif if(lf7)then if(s2.eq.'F ')ndao=7 if(s2.eq.'G ')ndao=9 if(s2.eq.'H ')ndao=11 if(s2.eq.'I ')ndao=13 else if(s2.eq.'F ')ndao=10 if(s2.eq.'G ')ndao=15 if(s2.eq.'H ')ndao=21 if(s2.eq.'I ')ndao=28 endif 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 ============================================================ function number(s) character*(*) s logical number number=s.eq.'1'.or.s.eq.'2'.or.s.eq.'3'.or.s.eq.'4'.or.s.eq.'5' 1 .or.s.eq.'6'.or.s.eq.'7'.or.s.eq.'8'.or.s.eq.'9'.or.s.eq.'0' return end c ============================================================ subroutine 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*50 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(6,*)'Corrected:',ic write(6,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 call w36('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 ============================================================== 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,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 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.-5)s2='H ' if(ist(ishell).eq. 6)s2='I ' 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 ')then if(ld5)then ndao=5 else ndao=6 endif endif if(lf7)then if(s2.eq.'F ')ndao=7 if(s2.eq.'G ')ndao=9 if(s2.eq.'H ')ndao=11 if(s2.eq.'I ')ndao=13 else if(s2.eq.'F ')ndao=10 if(s2.eq.'G ')ndao=15 if(s2.eq.'H ')ndao=21 if(s2.eq.'I ')ndao=28 endif 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 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*2 at character*1 ok,ty parameter (nat0=190,nor=80,nort=700,n1e=28) c n1e: c 1 c 2 c 3 c 4 c 5 c 6 c 7 c 8 c 9 c 10 c 11 c 12 c 13 c 14 c 15 c 16 c 17 (1/2) (r-rj) x grad W c 18 (1/2) (r-rj) x grad c 19 (1/2) (r-rj) x grad c 20 (1/2) (r-(ri+rj)/2) x grad T c 21 (1/2) (r-(ri+rj)/2) x grad c 22 (1/2) (r-(ri+rj)/2) x grad c 23 (1/2) rj x r / 2 D c 24 (1/2) rj x r / 2 c 25 (1/2) rj x r / 2 c 26 (1/2) ri x r / 2 E c 27 (1/2) ri x r / 2 c 28 (1/2) ri x r / 2 integer*4 u0 parameter (u0=28) c n1e ... number of kinds of one-electron integrals common/atoms/q(nat0),zsum,xv(3,nat0),iz(nat0) 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),ok common/const/pi,bohr common/numbers/z0,one,two,three,four,five,half,a10 dimension s(nort,nort),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), 4aag(u0,u0),aah(u0,u0),aaii(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,ldisk,leat 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(102),leat ) c write(6,*)'Overlap6 started' call inipoly i0=0 ione=1 itwo=2 call initfd(aa1,aad,aaf,aag,aah,aaii) c inquire(file='V.SCR',exist=ldisk) if(ldisk)then write(6,*)'V.SCR found on disk, skip calculation' return endif 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='W.SCR',form='unformatted') open(39,file='T.SCR',form='unformatted') open(40,file='D.SCR',form='unformatted') open(41,file='E.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,aag,aah,aaii,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 xx=(xi+xj)*half yy=(yi+yj)*half zz=(zi+zj)*half 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,aag,aah,aaii,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) anjx=dble(nxj(joo)) anjy=dble(nyj(joo)) anjz=dble(nzj(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) 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 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 = (1/2) r x grad 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 angular momentum (1/2) (r - rj) x grad b6(17,ioo,joo)=((oy-yj*sy)*dz-(oz-zj*sz)*dy)*sx*half b6(18,ioo,joo)=((oz-zj*sz)*dx-(ox-xj*sx)*dz)*sy*half b6(19,ioo,joo)=((ox-xj*sx)*dy-(oy-yj*sy)*dx)*sz*half c angular momentum (1/2) (r - (ri+rj)/2) x grad b6(20,ioo,joo)=((oy-yy*sy)*dz-(oz-zz*sz)*dy)*sx*half b6(21,ioo,joo)=((oz-zz*sz)*dx-(ox-xx*sx)*dz)*sy*half b6(22,ioo,joo)=((ox-xx*sx)*dy-(oy-yy*sy)*dx)*sz*half c GIAO like r x rj /2 b6(23,ioo,joo)=(oy*zj*sz-oz*yj*sy)*sx*half b6(24,ioo,joo)=(oz*xj*sx-ox*zj*sz)*sy*half b6(25,ioo,joo)=(ox*yj*sy-oy*xj*sx)*sz*half c GIAO like r x ri /2 b6(26,ioo,joo)=(oy*zi*sz-oz*yi*sy)*sx*half b6(27,ioo,joo)=(oz*xi*sx-ox*zi*sz)*sy*half b6(28,ioo,joo)=(ox*yi*sy-oy*xi*sx)*sz*half 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 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) do 335 kk=17,19 335 write(38)(sv(kk,ii,jj),jj=1,nao) do 336 kk=20,22 336 write(39)(sv(kk,ii,jj),jj=1,nao) do 337 kk=23,25 337 write(40)(sv(kk,ii,jj),jj=1,nao) do 338 kk=26,28 338 write(41)(sv(kk,ii,jj),jj=1,nao) 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 do 11 kk=33,41 11 close(kk) c if(lsprint)then 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) endif c if(lpprint)call wff('P.SCR','Dipole, AO, P',nao,s,nort) c if(lqprint)then open(33,file='Q.SCR',form='unformatted') do 41 i=1,nao read(33)(s(i,j),j=1,nao) read(33) read(33) read(33) read(33) 41 read(33) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, XX',ione) rewind 33 do 42 i=1,nao read(33) read(33)(s(i,j),j=1,nao) read(33) read(33) read(33) 42 read(33) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, XY',ione) rewind 33 do 43 i=1,nao read(33) read(33) read(33)(s(i,j),j=1,nao) read(33) read(33) 43 read(33) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, XZ',ione) rewind 33 do 44 i=1,nao read(33) read(33) read(33) read(33)(s(i,j),j=1,nao) read(33) 44 read(33) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, YY',ione) rewind 33 do 45 i=1,nao read(33) read(33) read(33) read(33) read(33)(s(i,j),j=1,nao) 45 read(33) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, YZ',ione) rewind 33 do 46 i=1,nao read(33) read(33) read(33) read(33) read(33) 46 read(33)(s(i,j),j=1,nao) call wmtr(s,nort,nao,'Quadrupole Overlaps, AO, ZZ',ione) close(33) endif c if(lvprint)then call wff('V.SCR','Dipole velocities, AO, M',nao,s,nort) call wff('W.SCR','Angular GIAO w, AO, M',nao,s,nort) call wff('T.SCR','Angular GIAO t, AO, M',nao,s,nort) call wff('D.SCR','Angular GIAO d, AO, M',nao,s,nort) call wff('E.SCR','Angular GIAO e, AO, M',nao,s,nort) endif c if(llprint)call wff('L.SCR','Angular moment, AO, L',nao,s,nort) 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 wff(s1,s2,nao,s,nort) character*(*) s1,s2 integer*4 nao,nort,i,j real*8 s(nort,nort) open(33,file=s1,form='unformatted') do 53 i=1,nao read(33)(s(i,j),j=1,nao) read(33) 53 read(33) call wmtr(s,nort,nao,s2//'X',1) rewind 33 do 54 i=1,nao read(33) read(33)(s(i,j),j=1,nao) 54 read(33) call wmtr(s,nort,nao,s2//'Y',1) rewind 33 do 55 i=1,nao read(33) read(33) 55 read(33)(s(i,j),j=1,nao) call wmtr(s,nort,nao,s2//'Z',1) close(33) 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 subroutine nppp(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,aag,aah,aaii,ld5,lf7,ia,ip) implicit none integer*4 u0,nat0,nor,npi,npii,nxi,nyi,nzi,icha, 1ia,ip,i1,i2,nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii parameter (u0=28,nat0=190,nor=80) character*1 shell logical ld5,lf7,lopt,lg94 real*8 cci,aai,aa1,aad,aaf,ccs,ccp,ccd,ccf,cci_G, 1cci_H,cci_I,aag,aah,aaii dimension cci(*),nxi(*),nyi(*),nzi(*), 1aai(u0,u0),aa1(u0,u0),aad(u0,u0),aaf(u0,u0), 1aag(u0,u0),aah(u0,u0),aaii(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 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 xxy xxz xzz yzz yyz 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(lf7)then npii=9 do 80 i1=1,npii do 80 i2=1,npi 80 aai(i1,i2)=aag(i1,i2) 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 if(lf7)then npii=11 do 81 i1=1,npii do 81 i2=1,npi 81 aai(i1,i2)=aah(i1,i2) endif 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 if(lf7)then npii=13 do 82 i1=1,npii do 82 i2=1,npi 82 aai(i1,i2)=aaii(i1,i2) endif return endif write(3,*)'ia=',ia,' ip = ',ip,' shell = ',shell call report(' Unknown shell encountered') end SUBROUTINE REPORT(STRING) CHARACTER*(*) STRING WRITE(6,*)STRING STOP END SUBROUTINE IOPAR IMPLICIT REAL*8(A-H,O-Z) implicit integer*4 (i-n) parameter (nat0=190,nor=80,nort=700) character*1 ok,ty,c 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), 7at(nat0),ok real*8 ktol,kxmax,kymax,kzmax logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,zmat,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94,lopt,lnmr,lapt,lader,lgiao, 1l333,nosym,tight,lspace,loose,lcomb,lpro,lrdox2,lhell,l66,lsuper, 1lall,lgrad,lssup,lsupinv,lplot,lrul,lpolar,ldoroa,lcroa,ldisk, 1lgvir,luvcd,lcnmr,lnno,lsppl,lci,li3,lmp2,lmp1,lmpdisk,lmpn4, 1lconfig,gconfig,lcall,lturbo,lonly,lhdisk,lciss,ltriplets, 1lroax,ldox31,lnormao,lfuzzy,lgradc,lrgr,lorg,laolo,ldipole, 1laol3,lcaz,lkcaz,lkprn,lsubst,lpara,lcazfac,lkorb,lcazcub, 1lcazexc,lbasred,lcazdip,lklim,ldamp,lord,lloc,ltxt,leat,lboys, 1lsxr,lab,lorban common/iopts/iopt(300) equivalence (iopt(1),ncnmr),(iopt(2),inormo ),(iopt(3),nek1 ), 1 (iopt(4),nek2 ),(iopt(5),nej1 ),(iopt(6),nej2 ), 1 (iopt(7),ioe ),(iopt(8),nstart ),(iopt(9),nend ), 1 (iopt(10),nprc),(iopt(11),multi ),(iopt(12),ndog ), 1 (iopt(13),nsi ),(iopt(14),nconfig),(iopt(15),nproc ), 1 (iopt(16),neci),(iopt(17),idiago ),(iopt(18),igradc), 1 (iopt(19),nev ),(iopt(20),mmit ),(iopt(21),iham ), 1 (iopt(22),ncut),(iopt(23),nk ),(iopt(24),nj ), 1 (iopt(25),nxmax),(iopt(26),nymax ),(iopt(27),nzmax ), 1 (iopt(28),idk) ,(iopt(29),nklim),(iopt(30),ilebed), 1 (iopt(31),idiff),(iopt(32),ibiter) common/ropts/ropt(300) equivalence (ropt( 1),wexc),(ropt( 2),rfuz), 1 (ropt( 3),ox ),(ropt( 4),oy ),(ropt( 5),oz), 1 (ropt( 6),alx ),(ropt( 7),aly ),(ropt( 8),alz), 1 (ropt( 9),ktol),(ropt(10),akmax), 1 (ropt(11),kxmax),(ropt(12),kymax ),(ropt(13),kzmax ), 1 (ropt(14),px ),(ropt(15),py ),(ropt(16),pz), 1 (ropt(17),redorb),(ropt(18),rklim),(ropt(19),eklim), 1 (ropt(20),pr),(ropt(21),deklim),(ropt(22),btol) 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( 28),zmat ),(lopt( 29),lnmr ),(lopt( 30),lapt ), 1(lopt( 31),lader ),(lopt( 32),lgiao ),(lopt( 33),l333 ), 1(lopt( 34),nosym ),(lopt( 35),tight ),(lopt( 36),lspace ), 1(lopt( 37),loose ),(lopt( 38),lcomb ),(lopt( 39),lpro ), 1(lopt( 40),lrdox2 ),(lopt( 41),lhell ),(lopt( 42),l66 ), 1(lopt( 43),lsuper ),(lopt( 44),lall ),(lopt( 45),lgrad ), 1(lopt( 46),lssup ),(lopt( 47),lsupinv ),(lopt( 48),lplot ), 1(lopt( 49),lrul ),(lopt( 50),lpolar ),(lopt( 51),ldoroa ), 1(lopt( 52),lcroa ),(lopt( 53),ldisk ),(lopt( 54),lgvir ), 1(lopt( 55),luvcd ),(lopt( 56),lcnmr ),(lopt( 57),lnno ), 1(lopt( 58),lsppl ),(lopt( 59),lci ),(lopt( 60),li3 ), 1(lopt( 61),lmp2 ),(lopt( 62),lmp1 ),(lopt( 63),lmpdisk ), 1(lopt( 64),lmpn4 ),(lopt( 65),lconfig ),(lopt( 66),gconfig ), 1(lopt( 67),lcall ),(lopt( 68),lturbo ),(lopt( 69),lonly ), 1(lopt( 70),lhdisk ),(lopt( 71),lciss ),(lopt( 72),ltriplets), 1(lopt( 73),lroax ),(lopt( 74),ldox31 ),(lopt( 75),lmcon ), 1(lopt( 76),lnormao),(lopt( 77),lfuzzy ),(lopt( 78),lgradc ), 1(lopt( 79),lrgr ),(lopt( 80),lorg ),(lopt( 81),laolo ), 1(lopt( 82),ldipole),(lopt( 83),laol3 ),(lopt( 84),lcaz ), 1(lopt( 85),lkcaz ),(lopt( 86),lkprn ),(lopt( 87),lkorb ), 1(lopt( 88),lsubst ),(lopt( 89),lpara ),(lopt( 90),lcazfac ), 1(lopt( 91),lcazcub),(lopt( 92),lbasred ),(lopt( 93),lcazexc ), 1(lopt( 94),lcazdip),(lopt( 95),lklim ),(lopt( 96),ldamp ), 1(lopt( 97),lord ),(lopt( 99),lcmx ),(lopt(100),lloc ), 1(lopt(101),ltxt ),(lopt(102),leat ),(lopt(103),lboys ), 1(lopt(104),lsxr ),(lopt(105),lab ),(lopt(106),lorban ) c c lboys - make localized MOs (Boys' scheme) c c other occupied options: c 300 - read g94 output/MO coeffs. in a fixed format defined in roainp.f c 299 - read g94 output/MO coeffs. in a g94 default format c 298 - also occupied orbitals for NMR coupling (dcoupl.f) c 297 - try the gaussian fit/diamagn coupling c CHARACTER*15 key,contr character*80 ifn,s80 common/filenames/ifn dimension icf(nort*2) c do 11 i=1,30 11 ropt(i)=0.0d0 do 3 i=1,300 iopt(i)=0 3 lopt(i)=.false. idiff=22 btol=1.0d-10 ibiter=1000 redorb=1.0d-6 ifn='G.OUT' rfuz=0.01d0 igradc=7 idiago=1 neci=2 icha=0 lrul=.true. lciss=.true. LAN=.true. ldisk=.true. lc5=.true. lg94=.true. ldoroa=.true. lloc=.true. lrdox2=.true. lhdisk=.true. ldipole=.true. ioe=0 nmo=0 ncnmr=10 inormo=2 nprc=3 lkcaz=.true. nxmax=10 nymax=10 nzmax=10 px=1.0d0 py=1.0d0 pz=1.0d0 alx=10.0d0 aly=10.0d0 alz=10.0d0 kxmax=30.0d0 kymax=30.0d0 kzmax=30.0d0 rklim=40.0d0 idk=4 wexc=1.0d0/589.0d-7/219470.0d0 c OPEN(2,FILE='SOS.OPT',STATUS='OLD',IOSTAT=IS) IF(IS.EQ.0)THEN WRITE(3,*)' FILE SOS.OPT FOUND ...' 2001 READ(2,1111)key write(6,1111)key 1111 FORMAT(A15) IF(key(1:4).EQ.'RFUZ')READ(2,*)RFUZ IF(key(1:5).EQ.'FUZZY')READ(2,*)LFUZZY IF(key(1:5).EQ.'SUBST')READ(2,*)LSUBST IF(key(1:4).EQ.'WEXC')then read(2,200)s80 200 format(a80) do 1 istart=1,len(s80) c=s80(istart:istart) 1 if( c.ne.'0'.and.c.ne.'1'.and.c.ne.'2'.and.c.ne.'3'. 1 and.c.ne.'4'.and.c.ne.'5'.and.c.ne.'6'.and.c.ne.'7'. 1 and.c.ne.'8'.and.c.ne.'9'.and.c.ne.'.'.and.c.ne.' '. 1 and.c.ne.'e'.and.c.ne.'d'.and.c.ne.'E'.and.c.ne.'D')goto 2 2 READ(s80(1:istart-1),*)WEXC if(s80(istart:istart+1).eq.'nm')WEXC=1.0d7/219470.0d0/WEXC if(s80(istart:istart+1).eq.'Nm')WEXC=1.0d7/219470.0d0/WEXC if(s80(istart:istart+1).eq.'NM')WEXC=1.0d7/219470.0d0/WEXC if(s80(istart:istart+1).eq.'nM')WEXC=1.0d7/219470.0d0/WEXC if(s80(istart:istart+1).eq.'cm')WEXC=WEXC/219470.0d0 if(s80(istart:istart+1).eq.'Cm')WEXC=WEXC/219470.0d0 if(s80(istart:istart+1).eq.'CM')WEXC=WEXC/219470.0d0 if(s80(istart:istart+1).eq.'cM')WEXC=WEXC/219470.0d0 endif IF(key(1:4).EQ.'BOYS')READ(2,*)lboys IF(key(1:4).EQ.'BTOL')READ(2,*)btol IF(key(1:4).EQ.'BITE')READ(2,*)ibiter IF(key(1:7).EQ.'DIPOLES')READ(2,*)LDIPOLE IF(key(1:6).EQ.'SCALED') READ(2,*)SCALED IF(key(1:5).EQ.'RAMAN')READ(2,*)RAMAN IF(key(1:5).EQ.'ABINI')READ(2,*)ABINI IF(key(1:6).EQ.'NORMAO')READ(2,*)LNORMAO IF(key(1:4).EQ.'NCUT')READ(2,*)NCUT IF(key(1:3).EQ.'TTT') READ(2,*)TTT IF(key(1:5).EQ.'NPROC') READ(2,*)nproc IF(key(1:3).EQ.'CC5') READ(2,*)CC5 IF(key(1:3).EQ.'XYZ') READ(2,*)XYZ if(key(1:3).eq.'NMO')read(2,*)nmo if(key(1:3).eq.'FIL')read(2,'(a)')ifn if(key(1:3).eq.'SPR')read(2,*)LSPRINT if(key(1:4).eq.'LEAT')read(2,*)leat if(key(1:4).eq.'NKNJ')read(2,*)NK,NJ if(key(1:3).eq.'PPR')read(2,*)LPPRINT if(key(1:4).eq.'MCON')read(2,*)LMCON if(key(1:4).eq.'IHAM')read(2,*)IHAM if(key(1:3).eq.'QPR')read(2,*)LQPRINT if(key(1:3).eq.'LD5')read(2,*)LD5 if(key(1:3).eq.'LEV')read(2,*)levprint if(key(1:3).eq.'VPR')read(2,*)lvprint if(key(1:3).eq.'LF7')read(2,*)lf7 if(key(1:3).eq.'LPR')read(2,*)llprint if(key(1:3).eq.'LAN')read(2,*)LAN if(key(1:3).eq.'LDO')read(2,*)LDO if(key(1:3).eq.'LD2')read(2,*)LD2 if(key(1:3).eq.'LSP')read(2,*)LSP if(key(1:3).eq.'LSC')read(2,*)LSCF if(key(1:3).eq.'RES')read(2,*)LRES if(key(1:3).eq.'REC')read(2,*)LREC if(key(1:3).eq.'LVV')read(2,*)LVV if(key(1:3).eq.'LEX')read(2,*)LEX if(key(1:3).eq.'IOE')read(2,*)IOE if(key(1:3).eq.'CND')read(2,*)LCNDO if(key(1:3).eq.'SPE')read(2,*)FSPEC if(key(1:3).eq.'LC5')read(2,*)LC5 if(key(1:3).eq.'LG9')read(2,*)LG94 if(key(1:3).eq.'NMR')read(2,*)LNMR if(key(1:3).eq.'APT')read(2,*)LAPT if(key(1:3).eq.'ADE')read(2,*)LADER if(key(1:3).eq.'GIA')read(2,*)LGIAO if(key(1:3).eq.'LOO')read(2,*)LOOSE if(key(1:3).eq.'333')read(2,*)L333 if(key(1:5).eq.'NOSYM')read(2,*)NOSYM if(key(1:5).eq.'TIGHT')read(2,*)TIGHT if(key(1:5).eq.'COMBI')read(2,*)lcomb if(key(1:5).eq.'PROJE')read(2,*)lpro if(key(1:5).eq.'LRDOX')read(2,*)LRDOX2 if(key(1:5).eq.'DOX31')read(2,*)LDOX31 if(key(1:4).eq.'HELL')read(2,*)LHELL if(key(1:4).eq.'CO66')read(2,*)L66 if(key(1:4).eq.'CHAR')read(2,*)icha if(key(1:5).eq.'SUPER')read(2,*)LSUPER if(key(1:4).eq.'LALL')read(2,*)LALL if(key(1:5).eq.'LGRAD')read(2,*)LGRAD if(key(1:4).eq.'LSXR')read(2,*)LSXR if(key(1:5).eq.'GRADC')read(2,*)LGRADC if(key(1:4).eq.'SSUP')read(2,*)LSSUP if(key(1:3).eq.'VCD')read(2,*)LSUPINV if(key(1:4).eq.'PLOT')read(2,*)LPLOT if(key(1:8).eq.'ROARULES')read(2,*)LRUL if(key(1:5).eq.'POLAR')read(2,*)LPOLAR if(key(1:5).eq.'DOROA')read(2,*)LDOROA if(key(1:5).eq.'ORBAN')read(2,*)LORBAN if(key(1:4).eq.'ROAX')read(2,*)LROAX if(key(1:4).eq.'CROA')read(2,*)LCROA if(key(1:4).eq.'DISK')read(2,*)LDISK if(key(1:4).eq.'UVCD')read(2,*)LUVCD c for open shell, make also integrals between alpha and beta MOs: if(key(1:3).eq.'LAB')read(2,*)lab if(key(1:4).eq.'CNMR')read(2,*)LCNMR if(key(1:5).eq.'NCNMR')read(2,*)NCNMR if(key(1:6).eq.'INORMO')read(2,*)INORMO if(key(1:3).eq.'NNO')read(2,*)LNNO if(key(1:4).eq.'SPPL')read(2,*)LSPPL if(key(1:2).eq.'CI')read(2,*)LCI if(key(1:3).eq.'LI3')read(2,*)LI3 if(key(1:4).eq.'NEK1')read(2,*)NEK1 if(key(1:4).eq.'NEK2')read(2,*)NEK2 if(key(1:4).eq.'NEJ1')read(2,*)NEJ1 if(key(1:4).eq.'NEJ2')read(2,*)NEJ2 if(key(1:3).eq.'MP2')read(2,*)LMP2 if(key(1:3).eq.'MP1')read(2,*)LMP1 if(key(1:6).eq.'MPDISK')read(2,*)LMPDISK if(key(1:4).eq.'MPN4')read(2,*)LMPN4 if(key(1:4).eq.'NEND')read(2,*)NEND if(key(1:6).eq.'NSTART')read(2,*)NSTART if(key(1:5).eq.'TURBO')read(2,*)LTURBO if(key(1:4).eq.'CALL')read(2,*)lcall if(key(1:6).eq.'LHDISK')read(2,*)lhdisk if(key(1:4).eq.'ONLY')read(2,*)lonly if(key(1:4).eq.'NECI')read(2,*)neci if(key(1:4).eq.'IDIA')read(2,*)idiago if(key(1:4).eq.'IGRA')read(2,*)igradc if(key(1:3).eq.'NEV')read(2,*)nev if(key(1:4).eq.'MMIT')read(2,*)mmit if(key(1:6).eq.'BASRED')read(2,*)lbasred if(key(1:6).eq.'BASLIM')read(2,*)redorb c check basis redundancy if(key(1:4).eq.'LORG')read(2,*)lorg if(key(1:4).eq.'AOLO')read(2,*)laolo if(key(1:4).eq.'AOL3')read(2,*)laol3 if(key(1:3).eq.'RGR')read(2,*)lrgr if(key(1:5).eq.'LCISS')read(2,*)lciss if(key(1:5).eq.'TRIPL')read(2,*)ltriplets if(key(1:5).eq.'LZMAT')read(2,*)zmat if(key(1:4).eq.'LCMX')read(2,*)lcmx if(key(1:6).eq.'CONFIG')then read(2,*)lconfig read(2,*)multi,ndog,nsi write(3,*)'electron configuration defined' endif if(key(1:7).eq.'GCONFIG')then read(2,*)gconfig open(33,file='CONFIG.SCR',form='unformatted') read(2,*)nconfig write(33)nconfig write(3,*)nconfig,' configurations' ispin=0 ndog=0 nsi=0 do ii=1,nconfig read(2,5555)(icf(i),i=1,2*nmo) write(33)(icf(i),i=1,2*nmo) write(3,5555)(icf(i),i=1,2*nmo) 5555 format(80I1) if(ii.eq.1)then do i=1,2*nmo if((i/2)*2.eq.i)then c beta orbitals if(icf(i).eq.1)ispin=ispin-1 if(icf(i-1).eq.1.and.icf(i).eq.1)ndog=ndog+1 if(icf(i-1)+icf(i).eq.1)nsi=nsi+1 else c alpha orbitals if(icf(i).eq.1)ispin=ispin+1 endif enddo endif enddo write(3,*)'General electron configurations defined' multi=2*ispin+1 write(3,3111)ispin,multi,ndog,nsi 3111 format(' ground state:',/, 1 ' spin ',I4,/, 1 ' multiplicity ',I4,/, 1 ' number of doubly occupied MOs ',I4,/, 1 ' number of singly occupied MOs ',I4) close(33) endif c lcaz-do Cazimir calculation, lkcaz - n and k type input, c nxmax,nymax,nzmax -maximum photon quantum number c ox,oy,oz - box origin c alx,aly,alz - box dimensions c kxmax,kymax,kzmax -maximum photon wave vector c px,py,pz - k-increment for integration c lkprn - print k-contributions in a projection c idk - k-integration path IF(key(1:7).EQ.'CAZIMIR')READ(2,*)lcaz IF(key(1:7).EQ.'KCAZIMI')READ(2,*)lkcaz IF(key(1:6).EQ.'CAZIN3' )READ(2,*)nxmax,nymax,nzmax IF(key(1:6).EQ.'CAZIO3' )READ(2,*)ox,oy,oz IF(key(1:6).EQ.'CAZIA3' )READ(2,*)alx,aly,alz IF(key(1:6).EQ.'CAZIK3' )READ(2,*)kxmax,kymax,kzmax IF(key(1:6).EQ.'CAZIP3' )READ(2,*)px,py,pz IF(key(1:6).EQ.'CAZIPR' )READ(2,*)pr IF(key(1:6).EQ.'CAZLEB' )READ(2,*)ilebed IF(key(1:7).EQ.'CAZLORB' )READ(2,*)lkorb IF(key(1:7).EQ.'CAZPARA' )READ(2,*)lpara IF(key(1:6).EQ.'CAZFAC' )READ(2,*)lcazfac IF(key(1:6).EQ.'CAZCUB' )READ(2,*)lcazcub IF(key(1:6).EQ.'CAZEXC' )READ(2,*)lcazexc IF(key(1:6).EQ.'CAZDIP' )READ(2,*)lcazdip IF(key(1:5).EQ.'AKMAX' )READ(2,*)akmax IF(key(1:4).EQ.'KTOL' )READ(2,*)ktol IF(key(1:5).EQ.'LKPRN' )READ(2,*)lkprn IF(key(1:6).EQ.'CAZIDK' )READ(2,*)idk IF(key(1:5).EQ.'RKLIM')READ(2,*)rklim IF(key(1:5).EQ.'EKLIM')READ(2,*)eklim IF(key(1:5).EQ.'DKLIM')READ(2,*)deklim IF(key(1:5).EQ.'LKLIM')READ(2,*)lklim IF(key(1:5).EQ.'LDAMP')READ(2,*)ldamp IF(key(1:4).EQ.'LORD')READ(2,*)lord IF(key(1:4).EQ.'LLOC')READ(2,*)lloc IF(key(1:5).EQ.'IDIFF')READ(2,*)idiff IF(key(1:5).EQ.'NKLIM')READ(2,*)nklim IF(key(1:4).EQ.'LTXT')READ(2,*)ltxt if(key(1:3).eq.'###')then read(2,*)i read(2,*)lopt(i) endif IF(key(1:3).EQ.'END')GOTO 2002 IF(key(1:6).EQ.'REMARK')GOTO 2001 backspace 2 read(2,1111)contr if(contr.eq.key)then write(3,1111)contr write(6,1111)contr call report('Unknown option selected') endif GOTO 2001 2002 CLOSE(2) ENDIF c if(nmo.eq.0)then call w36('trying to get nmo from nao') open(8,file=ifn) 801 read(8,800,end=8001,err=8001)s80 800 format(a80) c checkpoint format: if(s80(1:6).eq.'Charge')read(s80(50:61),*)icha if(s80(1:12).eq.'Multiplicity')then read(s80(50:61),*)multi write(6,*)'Multiplicity ',multi if(multi.gt.1)lconfig=.true. endif if(s80(1:22).eq.'Alpha Orbital Energies')then read(s80(50:61),*)nmo goto 8002 endif 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 write(3,*)' charge read from output ',icha,multi if(multi.gt.1)lconfig=.true. endif if(s80(7:22).eq.' basis functions ')then read(s80(1:6),*)nmo goto 8002 endif goto 801 8001 close(8) call report('NMO/NAO not found') 8002 write(6,*)'nmo set to ',nmo write(3,*)'nmo set to ',nmo close(8) endif if(ldo.and.lgiao)call report('Conflicting options ldo/lgiao') if(lmp2.and.(nstart.eq.0.or.nend.eq.0)) 1call report('MP2 orbital range must be defined') RETURN END c ============================================================ subroutine w36(s) character*(*) s write(6,*)s return end 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 subroutine initfd(aa1,aad,aaf,aag,aah,aai) 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) dimension aag(u0,u0),aah(u0,u0),aai(u0,u0) logical lopt,lturbo common/opts/lopt(300),icha equivalence (lopt(68),lturbo) real*8,allocatable::Coef2P(:,:),Coef2C(:,:) real*8 fac,f0 common/ffac/f0,fac(6) c do 77 i=1,u0 do 77 j=1,u0 aa1(i,j)=z0 aad(i,j)=z0 aaf(i,j)=z0 aag(i,j)=z0 aah(i,j)=z0 aai(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 xyz 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-1391.72064776868 -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 c g,h,i-functions (f for control): i=0 fac(i) = 1.0d0 Do 10 i = 1, (2*6) 10 fac(i) = dble(i)*fac(i-1) DO 1 L=1,6 MCart = ((L+1)*(L+2))/2 MPure = 2*L+1 allocate(Coef2P(MCart,MPure),Coef2C(MCart,MPure)) Call C2PGn1(L,MCart,MPure,Coef2P,Coef2C) ia=0 do 2 m =0,l if(m.eq.0)then imend=1 else imend=2 endif do 2 im=1,imend c im=1 (|lm>+|l-m>/sqrt( 2) c im=2 (|lm>-|l-m>/sqrt(-2) ia=ia+1 ix=0 do 2 Lx=0,L do 2 Ly=0,(L-Lx) c Lz=L-Lx-Ly ix=ix+1 if(L.eq.4)aag(ia,ix)=Coef2P(ix,ia) if(L.eq.5)aah(ia,ix)=Coef2P(ix,ia) 2 if(L.eq.6)aai(ia,ix)=Coef2P(ix,ia) 1 deallocate(Coef2P,Coef2C) return end c ================================================== function c2pcff(L,M,Lx,Ly) implicit real*8 (a-h,o-z) c c Generate coefficients for conversion from normalized cartesian c gaussians to pure spherical harmonics c c v(L,M) = Sum c2pcff(L,M,Lx,Ly,fac) v(Lx,Ly,Lz) c (Lx=0,L;Ly=0,L-Lx) c c L,M for spherical harmonic c Lx,Ly,Lz=L-Lx-Ly for cartesian c fac(n) = n! c c H. B. Schlegel May 1994 c Save Zero, Two Data Zero/0.D0/, Two/2.D0/ real*8 fac,f0 common/ffac/f0,fac(6) C Lz = L - Lx - Ly Ma = IAbs(M) c2pcff = Zero If(Ma.gt.L.or.Lz.lt.0) Return a = Sqrt(fac(2*Lx)*fac(2*Ly)*fac(2*Lz)*fac(L)/ $ (fac(2*L)*fac(Lx)*fac(Ly)*fac(Lz))) b = Sqrt(fac(L-Ma)/fac(L+Ma))/((2**L)*fac(L)) ilim = (L-Ma)/2 j = (L-Ma-Lz)/2 g = Zero If((2*j).eq.(L-Ma-Lz)) then do 20 i = 0, ilim c = Zero If(j.ge.0.and.j.le.i) then c = (fac(L)/(fac(i)*fac(L-i)))* $ (fac(2*L-2*i)*((-1)**i)/fac(L-Ma-2*i))* $ (fac(i)/(fac(j)*fac(i-j))) do 10 k = 0, j d = Zero If((Lx-2*k).ge.0.and.(Lx-2*k).le.Ma) then d = (fac(j)/(fac(k)*fac(j-k)))* $ (fac(Ma)/(fac(Lx-2*k)*fac(Ma+2*k-Lx))) e = Zero kmLx2 = k - Lx/2 If(M.eq.0.and.mod(Lx,2).eq.0) e=(-1)**kmLx2 kTMaLx = ((2*k+Ma-Lx)/2) If(M.gt.0.and.mod(IAbs(Ma-Lx),2).eq.0) $ e = Sqrt(Two)*(-1)**kTMaLx If(M.lt.0.and.mod(IAbs(Ma-Lx),2).eq.1) $ e = Sqrt(Two)*(-1)**kTMaLx g = g + c*d*e endIf 10 Continue endIf 20 Continue endIf c2pcff = a*b*g Return End c ================================================== Subroutine C2PGn1(L,MCart,MPure,Coef2P,Coef2C) Implicit Real*8(A-H,O-Z) Dimension Coef2P(MCart,MPure), Coef2C(MCart,MPure) real*8 fac,f0 common/ffac/f0,fac(6) c c M=0 case c do i=1,MCart do j=1,Mpure Coef2C(i,j)=0.0d0 enddo enddo lxyz = 0 do 10 Lx = 0, L do 10 Ly = 0, (L-Lx) Lz = L-Lx-Ly lxyz = lxyz+1 coef2p(lxyz,1) = c2pcff(L,0,Lx,Ly) Do 10 M = 1, L coef2p(lxyz,2*M) = c2pcff(L,M,Lx,Ly) 10 coef2p(lxyz,2*M+1) = c2pcff(L,-M,Lx,Ly) c c compute the overlap and form coef2c c i = 0 do 230 Lx1 = 0, L do 230 Ly1 = 0, (L-Lx1) Lz1 = L - Lx1 - Ly1 i = i + 1 j = 0 a1 = Sqrt(fac(Lx1)*fac(Ly1)*fac(Lz1)/ $ (fac(2*Lx1)*fac(2*Ly1)*fac(2*Lz1))) do 220 Lx2 = 0, L do 220 Ly2 = 0, (L-Lx2) Lz2 = L - Lx2 - Ly2 j = j + 1 a2 = Sqrt(fac(Lx2)*fac(Ly2)*fac(Lz2)/ $ (fac(2*Lx2)*fac(2*Ly2)*fac(2*Lz2))) Lx = Lx1 + Lx2 Ly = Ly1 + Ly2 Lz = Lz1 + Lz2 if(mod(Lx,2).eq.0.and.mod(Ly,2).eq.0.and.mod(Lz,2).eq.0) $ then s = a1*a2*fac(Lx)*fac(Ly)*fac(Lz)/ $ (fac(Lx/2)*fac(Ly/2)*fac(Lz/2)) do 210 k = 1, mpure 210 coef2c(i,k) = coef2c(i,k)+s*coef2p(j,k) endIf 220 Continue 230 Continue Return End subroutine readtm c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c parameter (nat0=190,nor=80,nort=700) character*2 at character*1 ok,ty character*80 s80 common/const/pi,bohr 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),ok 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/rundiff/idif,igeo,ibasis,ieigen common/atoms/z(nat0),zsum,xv(3,nat0),iz(nat0) logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,lzmat,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lconfig,lqq,lsymaos 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(28),lzmat ),(lopt(65),lconfig) dimension it(nor) character*2 tmsy(nat0) character*4 statesym real*8,allocatable::ds(:) c call w36('Turbomole files coord,basis and mos') c id5=6 if7=10 if(ld5)id5=5 if(lf7)if7=7 do 10000 i=1,nort do 10000 j=1,nort bij(i,j)=0.0d0 10000 cij(i,j)=0.0d0 c inquire(file='coord',exist=lqq) if(.not.lqq)call report('coord does not exist') igeo=igeo+1 open(4,file='coord') 1 read(4,4010)ok 4010 format(a1) if(ok.eq.'$')goto 400 goto 1 400 i=0 33 read(4,4010)ok if(ok.eq.'$')goto 401 backspace 4 i=i+1 if(i.gt.nat0)call report('too many atoms') read(4,4001)(xv(ix,i),ix=1,3),ok 4001 format(f20.14,2f22.14,6x,a1) z(i)=0.0d0 if(ok.eq.'h')z(i)=1 if(ok.eq.'c')z(i)=6 if(ok.eq.'n')z(i)=7 if(ok.eq.'o')z(i)=8 tmsy(i)(1:1)=ok if(z(i).lt.1)call report('unknown atom') goto 33 401 close(4) noa=i write(6,*)noa,' atoms' c inquire(file='mos',exist=lqq) if(.not.lqq)call report('mos does not exist') statesym='1234' nao=0 ieigen=ieigen+1 open(4,file='mos') naot=0 2 read(4,4009)s80 4009 format(a80) if(s80(16:25).eq.'eigenvalue')goto 299 goto 2 299 i=0 3 i=i+1 if(i.gt.nort)call report('too many eigenvalues') if(s80(7:9).ne.statesym)then naot=naot+nao endif statesym=s80(7:9) read(s80,4002)e(i),nao if(nao.gt.nort)call report('too many orbitals') if(naot+nao.gt.nort)call report('too many tot. orbitals') 4002 format(26x,d20.14,9x,i3) read(4,4003)(cij(i,naot+j),j=1,nao) 4003 format(4d20.14) be(i)=e(i) do 50 j=1,nao 50 bij(i,j+naot)=cij(i,j+naot) read(4,4009)s80 if(s80(1:4).ne.'$end') goto 3 close(4) nori=i nmo=nori nao=nao+naot write(6,*)nmo,' MOs read in, ',nao,' AOs' c inquire(file='basis',exist=lqq) if(.not.lqq)call report('basis does not exist') ibasis=ibasis+1 open(4,file='basis') nao=0 isht=0 c do 4 j=1,noa at(j)=' ' at(j)(1:1)=tmsy(j)(1:1) nop(j)=0 rewind 4 5 read(4,4010)ok if(ok.eq.tmsy(j)(1:1))goto 555 goto 5 555 k=0 7 read(4,4010)ok if(ok.eq.'*')goto 778 goto 7 778 read(4,4005)ipr,ok isht=isht+1 4005 format(i4,2x,a1) if(ok.eq.'s')nao=nao+1 if(ok.eq.'p')nao=nao+3 if(ok.eq.'d')nao=nao+id5 if(ok.eq.'f')nao=nao+if7 do 6 i=1,ipr k=k+1 if(k.gt.nor)call report('Too many primitives on one atom') if(ok.eq.'s')ty(j,k)='S' if(ok.eq.'p')ty(j,k)='P' if(ok.eq.'d')ty(j,k)='D' if(ok.eq.'f')ty(j,k)='F' c c temporarily record angular momentum for arrangement: if(ok.eq.'s')it(k)=0 if(ok.eq.'p')it(k)=1 if(ok.eq.'d')it(k)=2 if(ok.eq.'f')it(k)=3 c ish(j,k)=isht read(4,*)alpha(j,k),cc(j,k) c aa=alpha(j,k) c c SNO normalization factor to S: SNO=sqrt(sqrt( 8.0d0*aa**3/pi**3)) c c SNP normalization factor for P: SNP=sqrt(sqrt( 128.0d0*aa**5/pi**3)) c c SND normalization factor for D: SND=sqrt(sqrt( 2048.0d0*aa**7/pi**3)) SNDD=SND/sqrt(dble(3)) c c SNFx normalization factors for F(fxxx,fxxy,fxyz): SNF=sqrt(sqrt(32786.0d0*aa**9/pi**3)) SNF1=SNF/sqrt(dble(15)) SNF2=SNF/sqrt(dble(3)) SNF3=SNF c ccs(j,k)=cc(j,k)*SNO ccp(j,k)=cc(j,k)*SNP ccd(j,k,1)=cc(j,k)*SNDD ccd(j,k,2)=cc(j,k)*SND ccf(j,k,1)=cc(j,k)*SNF1 ccf(j,k,2)=cc(j,k)*SNF2 ccf(j,k,3)=cc(j,k)*SNF3 6 continue read(4,4010)ok if(ok.ne.'*')then backspace 4 goto 778 endif nop(j)=k write(8,*)k,' primitives for atom ',j c c reorder according to the angular momentum 403 do 402 k=1,nop(j)-1 if(it(k).gt.it(k+1))then call swi( it( k ), it( k+1)) call swi( ish(j,k ), ish(j,k+1)) call swh( ty(j,k ), ty(j,k+1)) call sw8(alpha(j,k ),alpha(j,k+1)) call sw8( cc(j,k ), cc(j,k+1)) call sw8( ccs(j,k ), ccs(j,k+1)) call sw8( ccp(j,k ), ccp(j,k+1)) call sw8( ccd(j,k,1), ccd(j,k+1,1)) call sw8( ccd(j,k,2), ccd(j,k+1,2)) call sw8( ccf(j,k,1), ccf(j,k+1,1)) call sw8( ccf(j,k,2), ccf(j,k+1,2)) call sw8( ccf(j,k,3), ccf(j,k+1,3)) goto 403 endif 402 continue c c adjpust reordered ish ist=ish(j,1) do 405 k=1,nop(j) if(ish(j,k).lt.ist)ist=ish(j,k) 405 it(k)=ish(j,k) c ish(j,1)=ist do 404 k=2,nop(j) if(it(k-1).ne.it(k))ist=ist+1 404 ish(j,k)=ist c 4 continue close(4) write(6,*)nao,' AOs' c zsum=0.0d0 nat=noa do 1000 iat=1,noa 1000 zsum=zsum+z(iat) noe=nint(zsum)-icha ndo=noe/2 ichm=0 mtp=1 nsa=noe-2*ndo nsb=0 c if(ibasis.gt.0)idif=ieigen-1 write(3,4004)idif,ieigen,igeo,ibasis write(6,4004)idif,ieigen,igeo,ibasis 4004 format(I4,' point of the differentiation:',/, 1 '----', ' -----------------------------',/, 2 ' Eigenvectors read',I4,' times',/, 3 ' Geometry read ',I4,' times',/, 4 ' Basis read ',I4,' times',/) inquire(file='symaos',exist=lsymaos) if(lsymaos)then write(6,*)'symmetrized AO coefficients found in symaos' allocate(ds(nao*nao)) call vz(ds,0.0d0,nao*nao) call writecon(0) call readsymaos('symaos',nao,ds) call trafosymmo(cij,ds,nao,nmo,nort) call trafosymmo(bij,ds,nao,nmo,nort) write(6,*)'cij transformed into non-symmetric AOs' write(3,*)'cij transformed into non-symmetric AOs' endif close(2) return end c subroutine sw8(a,b) real*8 a,b,c c=a a=b b=c return end c subroutine swi(a,b) integer*4 a,b,c c=a a=b b=c return end c subroutine swh(a,b) character*1 a,b,c c=a a=b b=c return end subroutine readsymaos(f,nao,ds) implicit none integer*4 nao,icum,icur,atom,shell,nort,ii,j,i parameter (nort=700) character*(*) f character*80 s80 real*8 ds(*) character*3 aot,typ integer*4 iat,isha common/aotypes/aot(nort),iat(nort),isha(nort) write(6,*)f open(9,file=f) 1 read(9,900,err=99,end=99)s80 write(6,900)s80 900 format(a80) 31 read(s80(1:35),*,err=88,end=88)icum,icur 3 read(s80(35:44),*,end=88,err=88)atom read(s80(48:52),*,end=88,err=88)shell typ=s80(54:56) write(6,*)icum,icur,atom,shell,typ j=0 do 2 ii=1,nao 2 if(iat(ii).eq.atom.and.isha(ii).eq.shell. 1and.aot(ii).eq.typ)j=ii write(6,*)'j=',j if(j.eq.0)then write(6,900)s80 write(6,*)'atom ',atom write(6,*)'shell ',shell write(6,*)'type ',typ call report('AO not found') endif read(s80(59:71),*)ds(icum+nao*(j-1)) read(9,900,end=99,err=99)s80 write(6,900)s80 if(s80(1:10).eq.' ') goto 3 goto 31 88 goto 1 99 close(9) if(icum.ne.nao)call report('nao does not agree') do i=1,nao write(6,*)'SAO',i,':' do j=1,nao if(dabs(ds(i+(j-1)*nao)).gt.0.0d0)write(6,89)j,ds(i+(j-1)*nao) 89 format(i8,f12.6) enddo enddo return end subroutine trafosymmo(cij,ds,nao,nmo,nort) implicit none integer*4 nao,nmo,nort,i,j,is real*8 cij(nort,nort),ds(*) real*8,allocatable:: tij(:,:) allocate(tij(nort,nort)) do 1 i=1,nmo do 1 j=1,nao tij(i,j)=0.0d0 do 1 is=1,nao 1 tij(i,j)=tij(i,j)+cij(i,is)*ds(is+nao*(j-1)) do 2 i=1,nmo do 2 j=1,nao 2 cij(i,j)=tij(i,j) return end subroutine normorb implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*2 at character*1 ok,ty parameter (nat0=190,nor=80,nort=700) integer*4 u0 parameter (u0=28) common/atoms/q(nat0),zsum,xv(3,nat0),iz(nat0) 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),ok 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 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), 4aag(u0,u0),aah(u0,u0),aaii(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 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 ) dimension anorm(nort) c call inipoly i0=0 call initfd(aa1,aad,aaf,aag,aah,aaii) 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,aag,aah,aaii,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,aag,aah,aaii,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(6,*)'atomic orbital norms before normalization:' write(3,3000)(anorm(i),i=1,nao) write(6,3000)anorm(1) write(6,3000)anorm(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,aag,aah,aaii,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 c ============================================================ subroutine writecon(ndiff) c writes control output implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,bohr character*2 at character*1 ok,ty parameter (nat0=190,nor=80,nort=700) 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),ok common/atoms/z(nat0),zsum,xv(3,nat0),iz(nat0) 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,lturbo equivalence (lopt(68),lturbo) 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 ) common/iopts/iopt(300) equivalence (iopt(1),ncnmr),(iopt(2),inormo ),(iopt(3),nek1 ), 1 (iopt(4),nek2 ),(iopt(5),nej1 ),(iopt(6),nej2 ), 1 (iopt(7),ioe ),(iopt(8),nstart ),(iopt(9),nend ), 1 (iopt(10),nprc),(iopt(11),multi ),(iopt(12),ndog ), 1 (iopt(13),nsi ),(iopt(14),nconfig),(iopt(15),nproc ), 1 (iopt(16),neci),(iopt(17),idiago ),(iopt(18),igradc), 1 (iopt(19),nev ),(iopt(20),mmit ),(iopt(21),iham ), 1 (iopt(22),ncut),(iopt(23),nk ),(iopt(24),nj ) character*3 aot integer*4 iat,isha common/aotypes/aot(nort),iat(nort),isha(nort) parameter (nshtypes=49) character*3 shtypes(nshtypes) data shtypes/'s ','px ','py ','pz ','d0 ','d1a','d1b', 1 'd2a','d2b', 1'f 0','f+1','f-1','f+2','f-2','f+3','f-3', 2'g 0','g+1','g-1','g+2','g-2','g+3','g-3','g+4','g-4', 2'h 0','h+1','h-1','h+2','h-2','h+3','h-3','h+4','h-4','h+5','h-5', 2'i 0','i+1','i-1','i+2','i-2','i+3','i-3','i+4','i-4','i+5','i-5', 2'i+6','i-6'/ character*1 spd(7) data spd/'s','p','d','f','g','h','i'/ c mtp=multi if(ndiff.eq.0)then 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')then idel=1 aot(iao)='s ' endif if(ty(i,j).eq.'P')then idel=3 aot(iao )='px ' aot(iao+1)='py ' aot(iao+2)='pz ' endif if(ty(i,j).eq.'L')then idel=4 aot(iao )='s ' aot(iao+1)='px ' aot(iao+2)='py ' aot(iao+3)='pz ' endif if(ty(i,j).eq.'D'.and.ld5)then idel=5 aot(iao )='d0 ' aot(iao+1)='d1a' aot(iao+2)='d1b' if(lturbo)then aot(iao+3)='d2a' aot(iao+4)='d2b' else aot(iao+3)='d2b' aot(iao+4)='d2a' endif endif if(ty(i,j).eq.'D'.and..not.ld5)then idel=6 aot(iao )='xx ' aot(iao+1)='yy ' aot(iao+2)='zz ' aot(iao+3)='xy ' aot(iao+4)='xz ' aot(iao+5)='yz ' endif if(ty(i,j).eq.'F'.and.lf7)then idel=7 aot(iao )='f 0' aot(iao+ 1)='f+1' aot(iao+ 2)='f-1' aot(iao+ 3)='f+2' aot(iao+ 4)='f-2' aot(iao+ 5)='f+3' aot(iao+ 6)='f-3' endif if(ty(i,j).eq.'F'.and..not.lf7)then L=3 idel=10 aot(iao )='xxx' aot(iao+ 1)='yyy' aot(iao+ 2)='zzz' aot(iao+ 3)='xxy' aot(iao+ 4)='xxz' aot(iao+ 5)='yyx' aot(iao+ 6)='yyz' aot(iao+ 7)='zzx' aot(iao+ 8)='zzy' aot(iao+ 9)='xyz' endif if(ty(i,j).eq.'G'.or.ty(i,j).eq.'H'.or.ty(i,j).eq.'I')then if(ty(i,j).eq.'G')L=4 if(ty(i,j).eq.'H')L=5 if(ty(i,j).eq.'I')L=6 idel=0 if(lf7)then do 6 m=0,L ipend=2 if(m.eq.0)ipend=1 do 6 ip=1,ipend aot(iao+idel)(1:1)=spd(L+1) if(ip.eq.1)then aot(iao+idel)(2:2)='+' else aot(iao+idel)(2:2)='-' endif if(ipend.eq.1)aot(iao+idel)(2:2)=' ' write(aot(iao+idel)(3:3),'(i1)')m 6 idel=idel+1 else do 7 ix=0,L do 7 iy=0,L-ix iz=L-ix-iy write(aot(iao+idel)(1:1),'(i1)')ix write(aot(iao+idel)(2:2),'(i1)')iy write(aot(iao+idel)(3:3),'(i1)')iz 7 idel=idel+1 endif endif do 200 k=0,idel-1 200 iat(iao+k)=i 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 endif c Turbomole shell numbering - restarts each atom, separate for each shell do ii=1,nshtypes iaold=0 isn=0 do i=1,nao if(iat(i).ne.iaold)isn=0 if(aot(i).eq.shtypes(ii))then isn=isn+1 isha(i)=isn endif iaold=iat(i) enddo enddo if(nao.gt.nort)call report('Too many basis functions') write(3,*)' AO type atom shell' do 5 i=1,nao 5 write(3,300)i,aot(i),iat(i),isha(i) 300 format(i5,1x,a3,1x,2i4) 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 c ===================================== subroutine vz(a,z,n) real*8 a(*),z integer*4 i,n do 1 i=1,n 1 a(i)=z return end subroutine init implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,bohr common/numbers/z0,one,two,three,four,five,half,a10 pi=dble(4)*datan(dble(1)) bohr=0.52917705993d0 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 subroutine ccao(cao,MCart,naoc,l,ila,ai) c calculates overlap matrix c fi = N x^lx y^ly z^lz epx(-alpha r^2) implicit none integer*4 nat0,nor,nort,ila,l,MCart,naoc,iz parameter (nat0=190,nor=80,nort=700) real*8 q,zsum,xv,cao(MCArt,naoc),ai,ein,xj,yj,zj common/atoms/q(nat0),zsum,xv(3,nat0),iz(nat0) real*8 cc,alpha,rcp,cij,bij,e,be integer*4 nat,nao,nmo,noe,ichm,mtp,ndo,nsa,nsb,noa,nori,nop, 1ish character*2 at character*1 ok,ty common/work/cc(nat0,nor),alpha(nat0,nor),rcp(nat0,nor), 1cij(nort,nort),bij(nort,nort),e(nort),be(nort),nat,nao,nmo, 5noe,ichm,mtp,ndo,nsa,nsb,noa,nori,nop(nat0),ty(nat0,nor), 1ish(nat0,nor),at(nat0),ok real*8 pi,bohr common/const/pi,bohr real*8 z0,one,two,three,four,five,half,a10 common/numbers/z0,one,two,three,four,five,half,a10 integer*4 u0 parameter (u0=28) integer*4 nxj(u0),nyj(u0),nzj(u0) real*8 ccj(u0),sv(nort),b6(u0), 1aaj(u0,u0),aad(u0,u0),aaf(u0,u0),aa1(u0,u0), 4aag(u0,u0),aah(u0,u0),aaii(u0,u0) logical ld5,lf7,lopt,lsprint integer*4 icha common/opts/lopt(300),icha equivalence (lopt(1),lsprint),(lopt(4),ld5),(lopt(8),lf7) integer*4 i0,lxyz,Lx,Ly,Lz,i,ioj,ja,jj,jo,joo,jp,jpp,npj,npjj real*8 tol,spi,xi,yi,zi,No,sai,xij,yij,zij,aj,tt,ec, 1ef,x2ij,y2ij,z2ij,aij,r2ij,t,oi,sx,sy,sz real*8 fac,f0 common/ffac/f0,fac(6) c call inipoly i0=0 call initfd(aa1,aad,aaf,aag,aah,aaii) tol=100.0d0 spi=sqrt(pi*pi*pi) sai=sqrt(ai*ai*ai) call vz(sv,z0,nao) if(lsprint)write(3,*)' Cartesian overlaps:' xi=xv(1,ila) yi=xv(2,ila) zi=xv(3,ila) lxyz = 0 do 2 Lx = 0, L do 2 Ly = 0, (L-Lx) Lz = L-Lx-Ly lxyz = lxyz+1 if(lsprint)write(3,6001)Lx,Ly,Lz,ai 6001 format('< x^',i1,' y^',i1,' z^',i1,' exp(-',f5.2,'r^2)|:') c No x^Lx y^ly z^lz exp(-ai r^2) No=fac(2*Lx)*fac(2*Ly)*fac(2*Lz)*spi No=No/(2**(2*l)*fac(Lx)*fac(Ly)*fac(Lz)*ai**l*sai) No=1.0d0/dsqrt(No) write(6,*)No,1.0d0/dsqrt(ein(2*Lx,ai)*ein(2*Ly,ai)*ein(2*Lz,ai)) 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,aag,aah,aaii,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 sx=oi(i0,Lx,nxj(joo),ai,aj,xij,t,xi,xj) sy=oi(i0,Ly,nyj(joo),ai,aj,yij,t,yi,yj) sz=oi(i0,Lz,nzj(joo),ai,aj,zij,t,zi,zj) 71 b6(joo)=sx*sy*sz*ef*ccj(joo)*No c if(npjj.ne.npj)then do 74 joo=1,npjj ioj=joo+jo do 74 jpp=1,npj 74 sv(ioj)=sv(ioj)+aaj(joo,jpp)*b6(jpp) else do 7 joo=1,npj ioj=joo+jo 7 sv(ioj)=sv(ioj)+b6(joo) endif 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue c jp 4 continue c ja if(lsprint) call nlist(3,sv,nao) C do 33 jj=1,nao 33 cao(lxyz,jj)=sv(jj) c c zero-out temporary buffer: call vz(sv,z0,nao) 2 continue c lx,ly,lz return end subroutine nlist(io,l,n) implicit none integer*4 n0 parameter (n0=10) integer*4 io,n,il(n0),nt,i,istart,iend real*8 l(*) istart=1 9 iend=min(istart+n0-1,n) nt=iend-istart+1 do 1 i=1,nt 1 il(i)=istart+i-1 write(io,601)(il(i),i=1,nt) 601 format(10i8) write(io,602)(l(i),i=istart,iend) 602 format(10f8.4) if(iend.lt.n)then istart=iend+1 goto 9 endif return end