PROGRAM redchk WRITE(6,*)'Basis redundancy check, in G.OUT' call init CALL IOPAR CALL ROAAI END subroutine init implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,clau,spt,bohr,debye common/numbers/z0,one,two,three,four,five,six,half,qua, 1a8,a10 pi=dble(4)*datan(dble(1)) spt=sqrt(dble(2)/pi) bohr=0.52917705993d0 debye=2.541747758d0 clau=137.037d0 c z0=dble(0) one=dble(1) two=dble(2) three=dble(3) four=dble(4) five=dble(5) six=dble(6) half=one/two qua=half*half a8=dble(8) a10=dble(10) return end SUBROUTINE IOPAR IMPLICIT REAL*8(A-H,O-Z) implicit integer*4 (i-n) parameter (nat0=7,nor=80,nort=1610) character*1 ok,ty character*2 at common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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,lzmat,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,ltt,lsubst,lpara,lcazfac,lkorb,lcazcub, 1lcazexc,lbasred,lcazdip,lklim,ldamp common/kelvin/T,ltt 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) common/ropts/ropt(30) 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) common/opts/lopt(300),icha,n1 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(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) c c other occupied options: 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 1 i=1,300 iopt(i)=0 1 lopt(i)=.false. 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. lrdox2=.true. lhdisk=.true. ldipole=.true. ioe=0 nmo=0 ncnmr=10 inormo=2 nprc=3 kmin=0 jmin=0 jmax=0 kmax=0 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 T=298.0d0 ltt=.false. idk=4 c OPEN(2,FILE='ROA.OPT',STATUS='OLD',IOSTAT=IS) IF(IS.EQ.0)THEN WRITE(6,*)' FILE ROA.OPT FOUND ...' 2001 READ(2,1111)key 1111 FORMAT(A15) IF(key(1:4).EQ.'N1N2')READ(2,*)N1 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')READ(2,*)WEXC 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.'JMIN')read(2,*)JMIN if(key(1:4).eq.'KMIN')read(2,*)KMIN if(key(1:4).eq.'NKNJ')read(2,*)NK,NJ if(key(1:4).eq.'KMAX')read(2,*)KMAX if(key(1:4).eq.'JMAX')read(2,*)JMAX 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.'LZM')read(2,*)LZMAT 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: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: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 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:6).eq.'CONFIG')then read(2,*)lconfig read(2,*)multi,ndog,nsi write(6,*)'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(6,*)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(6,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(6,*)'General electron configurations defined' multi=2*ispin+1 write(6,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 c ltt - non-zero temperature 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.'CAZLTT' )READ(2,*)ltt 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:6).EQ.'TEMPER' )READ(2,*)T 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:5).EQ.'NKLIM')READ(2,*)nklim 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(6,1111)contr call report('Unknown option selected') endif GOTO 2001 2002 CLOSE(2) ENDIF c if(nmo.eq.0)then write(6,*)'trying to get nmo from nao' 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 if(nmo.gt.nort)call report('too many orbitals') 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 SUBROUTINE ROAAI IMPLICIT REAL*8 (A-H,O-Z) implicit integer*4 (i-n) parameter (nat0=7,nor=80,nort=1610) 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/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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) 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 common/iopts/iopt(300) common/opts/lopt(300),icha,n1 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(iopt(1),ncnmr ),(iopt(2),inormo ),(iopt(7),ioe ), 1(lopt(65),lconfig),(iopt(11),multi ),(iopt(12),ndog ), 1(iopt(13),nsi ),(lopt(73),lroax ),(lopt(76),lnormao), 1(lopt(77),lfuzzy ),(lopt(78),lgradc ),(iopt(23),nk ), 1(iopt(24),nj ),(lopt(81),laolo ),(lopt(82),ldipole), 1(lopt(83),laol3 ),(lopt(84),lcaz ),(lopt(92),lbasred ), 1(lopt(94),lcazdip) common/rundiff/step,idif,igeo,ibasis,ieigen character*80 ifn common/filenames/ifn common/sobs/sbe(nort*2),ndogc,nsic,nsorb,isos(2*nort), 1isomo(2*nort) parameter (nkkjj0=nort*(nort+1)) logical lnknj common/kkjj/ckkjj(nkkjj0),enknj(nort),lnknj c C idif=0 lnknj=.false. do 1 i=1,nat0 do 1 j=1,nor 1 ish(i,j)=32000 c open(2,file='G.OUT') ieigen=0 igeo=0 ibasis=0 call roainp close(2) write(6,*)nat call writecon if(lnormao)call normorb call overlap6 call aoredcheck RETURN END subroutine roainp c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) c parameter (nat0=7,nor=80,nort=1610) character*2 at,s2,symbols(106) character*1 ok,ty common/const/pi,clau,spt,bohr,debye common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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/step,idif,igeo,ibasis,ieigen common/atoms/z(nat0),zsum,xv(3,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,n1 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*50 key 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 c c c Browse at the G94 or C5 output: ic=0 111 ic=ic+1 read(2,1001,end=666)key 1001 format(a49) 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(6,1001)key if(lzmat)then write(6,*)'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(6,1001)key if(.not.lzmat)then write(6,*)'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,*)ia,iza,xv(1,i),xv(2,i),xv(3,i) if(ig98.eq.1)read(2,*)ia,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 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' noa=nat noe=int(zsum)-icha if(lconfig)then write(6,*)' Open shell ',multi,noe ndo=0 nsa=(multi+noe-1)/2 nsb=noe-nsa else 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 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,*)idum1,idum2 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 ')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' cif(lbasred)then write(6,*)'Basis check only' return cendif if(nao.gt.nort)call report('too many orbitals') endif c BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB c c noa=nat c c if(ibasis.gt.0)idif=ieigen-1 if(idif.eq.ndiff.and.ibasis.gt.0)then 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',/) return endif goto 111 666 close(2) call report('end of the ab initio output file') 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,clau,spt,bohr,debye common/const/pi,clau,spt,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 aoredcheck c check if redundancies are present in AO basis implicit none integer*4 nat0,nor,nort,nat,nao,nmo,noe,ichm,mtp, 1ndo,nsa,nsb,noa,nori,nop,ish,i,j,k,l,ielim,ir,in,jn,kn character*2 at character*1 ok,ty parameter (nat0=7,nor=80,nort=1610) real*8 escf,cc,alpha,rcp,e,be, 1cot,s,sp,norm,tol,ropt,sp1,norm1 common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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 dimension cot(nort*nort),s(nort*nort),ir(nort) common/ropts/ropt(30) c tol=ropt(17) write(6,4000)nao,noa,tol 4000 format('Redundancy check ',I4,' orbitals read in',I4,' atoms',/, 1' Tolerance ',g10.3) c c read S-matrix open(33,file='SAO.SCR',form='unformatted') do 34 i=1,nao 34 read(33)(s((i-1)*nao+j),j=1,i) do 35 i=1,nao do 35 j=1,i-1 35 s((j-1)*nao+i)=s((i-1)*nao+j) close(33) c c make trial MOs: do 6 i=1,nao*nao 6 cot(i)=0.0d0 do 5 i=1,nao ir(i)=0 5 cot((i-1)*nao+i)=1.0d0 C ielim=0 c orthonormalize do 2 i=1,nao in=(i-1)*nao c subtract all non-redundant previous: do 1 j=1,i-1 jn=(j-1)*nao if(ir(j).eq.0)then sp=0.0d0 do 3 k=1,i kn=(k-1)*nao sp1=0.0d0 do 31 l=1,i 31 sp1=sp1+cot(jn+l)*s(kn+l) 3 sp=sp+sp1*cot(in+k) do 11 k=1,i 11 cot(in+k)=cot(in+k)-sp*cot(jn+k) endif 1 continue norm=0.0d0 do 4 k=1,i kn=(k-1)*nao norm1=0.0d0 do 41 l=1,i 41 norm1=norm1+cot(in+l)*s(kn+l) 4 norm=norm+norm1*cot(in+k) if(norm.gt.tol)then write(6,6000)i,norm 6000 format(' Orbital ',i4,' rest norm ',f12.8) else ir(i)=1 ielim=ielim+1 write(6,6001)i,norm 6001 format(' Orbital ',i4,' rest norm ',f12.8,' *** redundant **') endif sp=dsqrt(norm) do 2 k=1,i 2 cot(in+k)=cot(in+k)/sp write(6,6003)ielim,nao-ielim 6003 format(i6,' orbitals can be eliminated,',i5,' left') call report('Basis check only') end subroutine normorb implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*2 at character*1 ok,ty parameter (nat0=7,nor=80,nort=1610) integer*4 u0 parameter (u0=28) common/atoms/q(nat0),zsum,xv(3,nat0) common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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,clau,spt,bohr,debye common/numbers/z0,one,two,three,four,five,six,half,qua, 1a8,a10 dimension nxi(u0),nyi(u0),nzi(u0), 1nxj(u0),nyj(u0),nzj(u0),cci(u0),ccj(u0), 3sv(u0,nort),b6(u0,u0),aai(u0,u0),aaj(u0,u0), 4aad(u0,u0),aaf(u0,u0),aa1(u0,u0) logical lsprint,lpprint,lqprint,ld5,levprint,lvprint,llprint,lf7, 1ld2,SCALED,TTT,RAMAN,XYZ,ABINI,CC5,LAN,LDO,lvv,lsp,lex, 1lscf,lres,lcndo,fspec,lrec,lc5,lg94, 9lopt,lnno common/opts/lopt(300),icha,n1 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) 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(6,4000)nao,noa 4000 format(I4,' orbitals read in',I4,' atoms') C c first loop, determine norms: io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) if(npi.gt.u0)call report('npi > u0') jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppp(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) do 71 ioo=1,npi coe=con*cci(ioo) sx=oi(i0,nxi(ioo),nxj(joo),ai,aj,xij,t,xi,xj) sy=oi(i0,nyi(ioo),nyj(joo),ai,aj,yij,t,yi,yj) sz=oi(i0,nzi(ioo),nzj(joo),ai,aj,zij,t,zi,zj) 71 b6(ioo,joo)=sx*sy*sz*coe c if(npi.ne.npii.or.npjj.ne.npj)then do 74 ioo=1,npii do 74 joo=1,npjj ioj=joo+jo do 74 ipp=1,npi aiii=aai(ioo,ipp) do 74 jpp=1,npj ajjj=aaj(joo,jpp)*aiii 74 sv(ioo,ioj)=sv(ioo,ioj)+ajjj*b6(ipp,jpp) else do 7 ioo=1,npi do 7 joo=1,npj ioj=joo+jo 7 sv(ioo,ioj)=sv(ioo,ioj)+b6(ioo,joo) endif c 777 if(ish(ja,jp+1).gt.ish(ja,jp)) jo=jo+npjj c 5 continue 4 continue C if(ish(ia,ip+1).gt.ish(ia,ip))then c if last primitive of the orbitals io+1 ... io+npii, then write them do 33 ii=1,npii iorb=ii+io 33 anorm(iorb)=sv(ii,iorb) io=io+npii do 9 ii=1,u0 do 9 jj=1,nao 9 sv(ii,jj)=z0 endif 3 continue 2 continue c write(6,*)'atomic orbital norms before normalization:' write(6,3000)(anorm(i),i=1,nao) 3000 format(6f12.6) c c second loop, normalize C atom i, orbital io io=0 fac=1.0d0/sqrt(anorm(1)) do 8 ia=1,noa do 6 ip=1,nop(ia) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) ccs(ia,ip )=ccs(ia,ip )*fac ccp(ia,ip )=ccp(ia,ip )*fac ccd(ia,ip,1)=ccd(ia,ip,1)*fac ccd(ia,ip,2)=ccd(ia,ip,2)*fac ccf(ia,ip,1)=ccf(ia,ip,1)*fac ccf(ia,ip,2)=ccf(ia,ip,2)*fac ccf(ia,ip,3)=ccf(ia,ip,3)*fac do i1=1,15 cci_G(ia,ip,i1)=cci_G(ia,ip,i1)*fac enddo do i1=1,21 cci_H(ia,ip,i1)=cci_H(ia,ip,i1)*fac enddo do i1=1,28 cci_I(ia,ip,i1)=cci_I(ia,ip,i1)*fac enddo if(ish(ia,ip+1).gt.ish(ia,ip))then io=io+npii fac=1.0d0/sqrt(anorm(io+1)) endif 6 continue 8 continue return end c ============================================================ subroutine writecon c writes control output implicit real*8 (a-h,o-z) implicit integer*4 (i-n) common/const/pi,clau,spt,bohr,debye character*2 at character*1 ok,ty parameter (nat0=7,nor=80,nort=1610) common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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) 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 common/opts/lopt(300),icha,n1 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 ) c mtp=multi write(6,*) 'ROAAI control output' write(6,*) 'MOLECULAR PARAMETERS' write(6,*) '------------------------------------------' write(6,*) 'Number of basis functions : ',nao write(6,*) 'Number of electrons : ',noe write(6,*) 'Charge : ',ichm write(6,*) 'Multiplicity : ',mtp write(6,*) 'Number of doubly occupied orbitals : ',ndo write(6,*) 'Number of singly occupied a-orbitals: ',nsa write(6,*) 'Number of singly occupied b-orbitals: ',nsb write(6,*) 'Number of atoms : ',noa write(6,*) 'Sum of atomic charges : ',int(zsum) write(6,*) write(6,*) iao=1 do 1 i=1,noa write(6,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')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(6,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(6,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(6,10000)ish(i,j), 1 ty(i,j),j,alpha(i,j),cc(i,j) if(ty(i,j).eq.'L')write(6,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 c open(31,file='GEO.X') write(6,*)'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(6,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 ===================================== 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,six,half,qua, 1a8,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,six,half,qua, 1a8,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 REPORT(STRING) CHARACTER*(*) STRING WRITE(6,1000) 1000 FORMAT(80(1H*)) WRITE(6,*)STRING WRITE(6,1000) WRITE(6,2000) 2000 FORMAT(' Program terminated') CLOSE(3) STOP END C ============================================================ subroutine wmtr(A,n0,n,st,ic) c c ic=0 .. writes triangle of a supposedly symmetric matrix c ic=1 .. writes it all c IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) dimension A(n0,n) character*(*) st write(6,*)st write(6,*) N1=1 1 N3=MIN(N1+4,N) WRITE(6,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(6,17)LN,(A(LN,J),J=N1,iend) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) write(6,*) return end subroutine nppp(shell,npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) implicit none integer*4 u0,nat0,nor,npi,npii,nxi,nyi,nzi,icha,n1, 1ia,ip,i1,i2,nxig,nyig,nzig,nxih,nyih,nzih,nxii,nyii,nzii parameter (u0=28,nat0=7,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 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,n1 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 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(6,*)'ia=',ia,' ip = ',ip,' shell = ',shell call report(' Unknown shell encountered') 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=7,nor=80,nort=1610,n1e=16) integer*4 u0 parameter (u0=28) c n1e ... number of kinds of one-electron integrals common/atoms/q(nat0),zsum,xv(3,nat0) common/work/escf, 2cc(nat0,nor),alpha(nat0,nor), 3rcp(nat0,nor), 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,clau,spt,bohr,debye common/numbers/z0,one,two,three,four,five,six,half,qua, 1a8,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) 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 common/opts/lopt(300),icha,n1 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 ) c call inipoly i0=0 ione=1 itwo=2 call initfd(aa1,aad,aaf) c inquire(file='V.SCR',exist=ldisk) if(ldisk)return c tol=100.00d0 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(6,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') c C atom i, orbital io io=0 do 2 ia=1,noa xi=xv(1,ia) yi=xv(2,ia) zi=xv(3,ia) c c primitive functions on atom ia: do 3 ip=1,nop(ia) ai=alpha(ia,ip) call nppp(ty(ia,ip),npi,npii,cci,nxi,nyi,nzi, 1aai,aa1,aad,aaf,ld5,lf7,ia,ip) c C atom j, orbital jo jo=0 do 4 ja=1,noa xj=xv(1,ja) yj=xv(2,ja) zj=xv(3,ja) xij=xj-xi yij=yj-yi zij=zj-zi x2ij=xij*xij y2ij=yij*yij z2ij=zij*zij r2ij=x2ij+y2ij+z2ij c c primitive functions on atom ja: do 5 jp=1,nop(ja) aj=alpha(ja,jp) call nppp(ty(ja,jp),npj,npjj,ccj,nxj,nyj,nzj, 1aaj,aa1,aad,aaf,ld5,lf7,ja,jp) aij=ai+aj tt=one/aij t=sqrt(tt) ec=ai*aj*r2ij*tt ef=z0 if(ec.lt.tol)then ef=spi*exp(-ec) else goto 777 endif c C Loop over orbitals to which the primitives ip, jp contribute do 71 joo=1,npj con=ef*ccj(joo) 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 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 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(6,*)anorm,iorb if(.not.lnno)call report('Non-normalized aos') endif if(lsprint)write(6,*)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) 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) 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)then open(33,file='P.SCR',form='unformatted') do 36 i=1,nao read(33)(s(i,j),j=1,nao) read(33) 36 read(33) close(33) call wmtr(s,nort,nao,'Electric Dipole Overlaps, AO, PX',ione) open(33,file='P.SCR',form='unformatted') do 38 i=1,nao read(33) read(33)(s(i,j),j=1,nao) 38 read(33) close(33) call wmtr(s,nort,nao,'Electric Dipole Overlaps, AO, PY',ione) open(33,file='P.SCR',form='unformatted') do 40 i=1,nao read(33) read(33) 40 read(33)(s(i,j),j=1,nao) close(33) call wmtr(s,nort,nao,'Electric Dipole Overlaps, AO, PZ',ione) endif 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 open(33,file='V.SCR',form='unformatted') do 47 i=1,nao read(33)(s(i,j),j=1,nao) read(33) 47 read(33) call wmtr(s,nort,nao,'Dipole velocities, AO, DX',ione) rewind 33 do 48 i=1,nao read(33) read(33)(s(i,j),j=1,nao) 48 read(33) call wmtr(s,nort,nao,'Dipole velocities, AO, DY',ione) rewind 33 do 49 i=1,nao read(33) read(33) 49 read(33)(s(i,j),j=1,nao) call wmtr(s,nort,nao,'Dipole velocities, AO, DZ',ione) close(33) endif c if(llprint)then open(33,file='L.SCR',form='unformatted') do 50 i=1,nao read(33)(s(i,j),j=1,nao) read(33) 50 read(33) call wmtr(s,nort,nao,'Angular moment, AO, LX',ione) rewind 33 do 51 i=1,nao read(33) read(33)(s(i,j),j=1,nao) 51 read(33) call wmtr(s,nort,nao,'Angular moment, AO, LY',ione) rewind 33 do 52 i=1,nao read(33) read(33) 52 read(33)(s(i,j),j=1,nao) call wmtr(s,nort,nao,'Angular moment, AO, LZ',ione) close(33) endif 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 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,six,half,qua, 1a8,a10 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