program hamilton c **************** Program HAMILTON **************** c This program is designed to generate or read in a c list of determinants. You can generate c determinants for a CAS (Complete Active Space) of c orbitals or you can input your own list of c determinants. Next, if you wish, you may read in c the the one- and two-electron MO integrals and c form a Hamiltonian matrix over the determinants. c Finally, if you wish, you may diagonalize the c Hamiltonian matrix constructed over the c determinants generated. c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real(a-h,o-z) logical yesno, cas character*30 dfname, fname, dscrpt character*1 ans include "limits.h" c c Set core needed to: c 1 symmetric array of order MAXORB c 1 array of dimension canonical(MAXORB**4) c 1 symmetric array of order MAXDET c 1 integer vector of length MAXDET c total = (maxorb*(maxorb+1))/2 + c (maxorb*(maxorb+1))/2*(((maxorb*(maxorb+1))/2+1))/2 + c (maxdet*(maxdet+1))/2 + maxdet c c c FORTRAN definitions for standard intrinsics. c c GNU/XLF vs Mac c c and(i,j) iand(i,j) c xor(i,j) ixor(i,j) c rshift(jword,i) ishft(jword,-i) c dimension pcore(mxo2 + mxo4 + mxdet2 + 3*maxdet + maxdet**2 + & 2*maxorb + nbits) dimension n2pow(nbits) dimension iord(maxorb), invord(maxorb) dimension ipntr(10) c itotal = mxo2 + mxo4 + mxdet2 + 3*maxdet + maxdet**2 + 2*maxorb & + nbits + 10 call modinfo('hamilton ') call shwmem('hamilton ',itotal) call goon 10 continue c c Initialize variables c c ... powers of two c do 5 i=1,nbits n2pow(i)=2**(i-1) 5 continue do 6 i = 1, maxorb iord(i) = i invord(i) = i 6 continue write(*,*) write(*,*)' Input the number of orbitals in your system: ' read(*,*)iorb if(iorb.gt.maxorb)then write(*,1010)maxorb goto 10 endif 20 continue write(*,*) write(*,*)' Would you like to generate a determinant ' write(*,*)' list for a CAS (Complete Active Space) ' write(*,*)' of orbitals or would you prefer to input ' write(*,*)' your own list of determinants? < 1,2, or 3 > ' write(*,*) write(*,*)' 1.) CAS ' write(*,*) write(*,*)' 2.) Input determinants ' write(*,*) write(*,*)' 3.) Quit ' write(*,*) read(*,*)igo if(igo.eq.1)then cas = .true. elseif(igo.eq.2)then cas = .false. elseif(igo.eq.3)then call qmexit else write(*,*)' Correct response is < 1,2, or 3 >. ' write(*,*)' Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 20 else call qmexit endif endif nmdets = 0 icl = 0 iva = 0 nelact = 0 ipntr(1) = 1 ipntr(2) = ipntr(1) + maxdet call dodets(pcore(ipntr(1)),n2pow,iord,invord,iorb, & nmdets,icl,iva,nelact,cas) write(*,*) write(*,*)' Would you like to construct a Hamiltonian ' write(*,*)' using this list of determinants? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 25 else call qmexit endif 25 continue iorb2 = iorb*(iorb+1)/2 iorb4 = iorb2*(iorb2+1)/2 nmdets2 = nmdets*(nmdets+1)/2 nmdtsq = nmdets*nmdets ipntr(1) = 1 ipntr(2) = ipntr(1) + nmdets ipntr(3) = ipntr(2) + nmdets2 ipntr(4) = ipntr(3) + iorb2 ipntr(5) = ipntr(4) + iorb4 c c Open the MO integral file. c dscrpt = '1e- and 2e- MO integrals' ids = 24 dfname = 'mo_integrals.dat' idf = 16 fname = ' ' ifn = 30 iflnum = 30 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(.not.yesno)then write(*,*)' You must have the MO integrals available ' write(*,*)' on disk to construct the Hamiltonian. ' call qmexit endif call dohami(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)), & pcore(ipntr(4)),n2pow,iord,invord,iorb,iorb2, & iorb4,nmdets,nmdets2,icl,iva) c c Close 1e- and 2- MO integral file. c close(unit=30) ipntr(1) = 1 ipntr(2) = ipntr(1) + nmdets ipntr(3) = ipntr(2) + nmdets2 ipntr(4) = ipntr(3) + nmdtsq call wthami(pcore(ipntr(2)),pcore(ipntr(3)),nmdets) write(*,*)' Would you like to diagonalize this ' write(*,*)' Hamiltonian? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then ipntr(1) = 1 ipntr(2) = ipntr(1) + nmdets ipntr(3) = ipntr(2) + nmdets2 ipntr(4) = ipntr(3) + nmdtsq ipntr(5) = ipntr(4) + nmdets ipntr(6) = ipntr(5) + nmdets call dodiag(pcore(ipntr(2)),pcore(ipntr(3)),pcore(ipntr(4)), & pcore(ipntr(5)),nmdets) endif call qmexit stop 1010 format(' Maximum number of orbitals is ',I2,'. ',/, & ' Please try again. ') end subroutine dodets(ldets,n2pow,iord,invord,n,nmdets, & icl,iva,nelact,cas) implicit real (a-h,o-z) logical cas character*1 ans include "limits.h" dimension ldets(maxdet) dimension n2pow(nbits) dimension iord(maxorb), invord(maxorb) write(*,*) write(*,*)' Input the desired Ms value of your CI state.' read(*,*)refms write(*,*) write(*,*)' The following couple of questions deal ' write(*,*)' with orbital occupations. ' write(*,*)' You may have to go back to your SCF ' write(*,*)' calculation and look carefully at MO ' write(*,*)' coefficients and orbital energies in ' write(*,*)' order to appropriately determine orbital ' write(*,*)' occupations. ' write(*,*) write(*,*)' Input the number of orbitals that are: ' write(*,*)' doubly occupied in all determinants, ' read(*,*)icl if(icl.gt.0)then write(*,1010)(iord(i),i=1,icl) read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' Input the doubly occupied orbitals: ' read(*,*)(iord(i),i=1,icl) endif endif write(*,*)' Input the number of orbitals that are: ' write(*,*)' partially occupied (active), ' read(*,*)iva write(*,1020)(iord(i),i=icl+1,icl+iva) read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' Input the partially occupied (active) ' write(*,*)' orbitals: ' read(*,*)(iord(i),i=icl+1,icl+iva) endif c c Fill remaining IORD with virtual orbitals. c do 20 i = 1, icl+iva do 10 j = 1, n if(invord(j).eq.iord(i))then invord(j) = 0 goto 11 endif 10 continue 11 continue 20 continue do 40 i = icl+iva+1, n do 30 j = 1, n if(invord(j).ne.0)then iord(i) = invord(j) invord(j) = 0 goto 31 endif 30 continue 31 continue 40 continue c c Calculate INVORD vector. c do 60 i = 1, n do 50 j = 1, n if(iord(j).eq.i)then invord(i) = j goto 51 endif 50 continue 51 continue 60 continue write(*,*)' Input the number of electrons that are ' write(*,*)' in the active space. ' read(*,*)nelact if(cas)then call casgen(refms,ldets,n2pow,nmdets,iva,nelact) & else call usrgen(refms,ldets,n2pow,iord,invord,nmdets,icl, & iva,nelact) endif write(*,1000)nmdets write(*,*)' Would you like to see a list of the ' write(*,*)' determinants generated? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then call wtdets(ldets,n2pow,iord,invord,nmdets,icl,iva,nelact) endif return 1000 format(' The number of determinants generated is: ',I3) 1010 format(' Would you like to specify which orbitals ',/, & ' will be doubly occupied? < y or n > ',/, & ' the default is: ',18(1X,I2)) 1020 format(' Would you like to specify which orbitals ',/, & ' will be partially occupied? < y or n > ',/, & ' the default is: ',18(1X,I2)) end subroutine dohami(ldets,hami,h,g,n2pow,iord,invord,n, & n2,n4,nmdets,nmdets2,icl,iva) implicit real (a-h,o-z) character*1 ans include "limits.h" dimension h(n2), g(n4), hami(nmdets2), ldets(nmdets), n2pow(nbits) dimension iord(maxorb), invord(maxorb) c c Get the 1e- integrals. c call get1es(h,iord,invord,n,n2) c c Get the 2e- integrals. c call get2es(g,iord,invord,n,n4) c c Determine the desired Spin Multiplicity and c occupancies. c icl2 = 2*icl iva2 = 2*iva call formh(h,g,hami,ldets,n2pow,nmdets,nmdets2, & n,n2,n4,icl2,iva2) c c Write out the Hamiltonian matrix. c write(*,*)' Would you like to see the Hamiltonian? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' The Hamiltonian matrix ' call tprntd(hami,nmdets) endif return 1000 format(' The number of determinants generated is: ',I3) end subroutine wthami(hsym,hsqr,n) implicit real (a-h,o-z) logical yesno character*30 dfname, fname, dscrpt dimension hsym(n*(n+1)/2) dimension hsqr(n,n) dscrpt = 'Hamiltonian matrix' ids = 18 dfname = 'hamilton.out' idf = 12 fname = ' ' ifn = 30 iflnum = 32 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then ij = 0 do 20 i = 1,n do 10 j = 1,i ij = ij + 1 hsqr(i,j) = hsym(ij) hsqr(j,i) = hsqr(i,j) 10 continue 20 continue call wtmatd(fname,ifn,dscrpt,ids,hsqr,n,n, & iflnum) c c Close Hamiltonian file. c close(unit=32) endif return end subroutine dodiag(hami,vecs,vals,scr1,n) implicit real (a-h,o-z) logical yesno character*30 dfname, fname, dscrpt character*1 ans dimension hami((n*(n+1))/2) dimension vecs(n,n) dimension vals(n) dimension scr1(n) c c Diagonalize the array. c call diaghe(hami,n,vals,n,vecs,n,scr1) c c Write the eigenvalues to user. c write(*,*)' Would you like to see the eigenvalues ' write(*,*)' of the Hamiltonian? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' the eigenvalues of this matrix are: ' call matout(vals,n,1,n) endif write(*,*)' Would you like to add the nuclear repulsion energy ' write(*,*)' to these eigenvalues? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then c c Get the nuclear repusion. c dscrpt = 'nuclear repulsion' ids = 17 dfname = 'NucRep' idf = 6 fname = ' ' ifn = 30 iflnum = 31 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) c c Read in the nuclear repulsion. c if(yesno)then call rdmatd(fname,ifn,dscrpt,ids,repnuc,1,1,iflnum) c c Close nuclear repulsion file. c close(unit=31) else write(*,*)' What is the nuclear repulsion energy? ' read(*,*)repnuc endif do 20 i = 1,n vals(i) = vals(i) + repnuc 20 continue write(*,*)' The eigenvalues of this matrix plus the ' write(*,*)' nuclear repulsion energy are: ' call matout(vals,n,1,n) endif dscrpt = 'eigenvalues' ids = 11 dfname = 'hami_vals.out' idf = 13 fname = ' ' ifn = 30 iflnum = 34 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call wtmatd(fname,ifn,dscrpt,ids,vals,n,1, & iflnum) c c Close eigenvalues file. c close(unit=34) endif c c Write the eigenvectors to user. c write(*,*)' Would you like to see the eigenvectors ' write(*,*)' of the Hamiltonian? ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then write(*,*)' How many eigenvectors would you like to see? ' read(*,*)nmany write(*,1000)nmany call matout(vecs,n,nmany,n) endif dscrpt = 'eigenvectors' ids = 12 dfname = 'hami_vecs.out' idf = 13 fname = ' ' ifn = 30 iflnum = 33 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call wtmatd(fname,ifn,dscrpt,ids,vecs,n,n, & iflnum) c c Close eigenvectors file. c close(unit=33) endif return 1000 format(' The ',I3,' lowest eigenvectors of this ',/, & ' matrix are: ') end subroutine casgen(refms,ldets,n2pow,nmdets,iva,nelact) c c This subroutine generates a complete active space set of c determinants upon input of spin ms value, number of active c orbitals, and number of electrons. c implicit real(a-h,o-z) include "limits.h" dimension n2pow(nbits), ldets(maxdet) refms2=refms*2.0e0 arefms2=abs(refms2) nextra=nint(arefms2) if(refms2.lt.0.0e0)nextra=-nextra nbeta=(nelact-nextra)/2 nalpha=(nelact+nextra)/2 ibegb=n2pow(nbeta+1)-1 iendb=n2pow(iva+1)-1 ibega=n2pow(nalpha+1)-1 ienda=iendb c c Obtain beta distribution. c do 400 ibeta=ibegb,iendb if(ibtsum(ibeta).ne.nbeta)go to 400 c c 'Zero' beta spin orbital determinant positions. c idetin=0 c c Determine occupancy, ..... etc. c iocct=0 nbsym=1 do 100 nocc=1,iva ickocc=n2pow(nocc) if(and(ibeta,ickocc).eq.0)go to 100 iocct=iocct+1 ipos=2*nocc idetin=idetin+n2pow(ipos) if(iocct.eq.nbeta)go to 101 100 continue 101 continue c c Obtain alpha distribution. c do 300 ialpha=ibega,ienda if(ibtsum(ialpha).ne.nalpha)go to 300 c c 'Zero' alpha spin orbital determinant positions. c jdetin=idetin c c Determine occupancy, ..... etc. c iocct=0 ntsym=nbsym do 200 nocc=1,iva ickocc=n2pow(nocc) if(and(ialpha,ickocc).eq.0)go to 200 iocct=iocct+1 ipos=2*nocc-1 jdetin=jdetin+n2pow(ipos) if(iocct.eq.nalpha)go to 201 200 continue 201 continue call putdet(ldets,nmdets,jdetin) 300 continue 400 continue return end subroutine usrgen(refms,ldets,n2pow,iord,invord,nmdets, & icl,iva,nelact) implicit real(a-h,o-z) logical alpha,beta character*1 ans,ones(10) character*1 line(80) character*80 aline include "limits.h" dimension ldets(maxdet) dimension n2pow(nbits) dimension io(2) dimension iord(maxorb), invord(maxorb) equivalence (line(1),aline) data ones/'1','2','3','4','5','6','7','8','9','0'/ nmdets = 0 write(*,*)' Begin input of determinants: (orbital number ' write(*,*)' followed by spin) e.g. 1A1B2A2B, etc, ' write(*,*)' where A = alpha and B = beta spin. ' 97 continue do 1 i = 1, 80 line(i) = ' ' 1 continue nmdets = nmdets + 1 write(*,1000)nmdets read(*,'(a80)')aline c c Parse line for orbital numbers and As and Bs and c determine LDETS value. c iorbc = 0 do 3 i = 1, 80 alpha = .false. beta = .false. do 2 j = 1, 9 if(line(i).eq.ones(j))then iorbc = iorbc + 1 io(iorbc) = j goto 99 endif 2 continue if(line(i).eq.ones(10))then iorbc = iorbc + 1 io(iorbc) = 0 goto 99 elseif(line(i).eq.'A'.or.line(i).eq.'a')then alpha = .true. goto 98 elseif(line(i).eq.'B'.or.line(i).eq.'b')then beta = .true. goto 98 elseif(line(i).eq.' ')then goto 99 endif write(*,1010)line(i),nmdets write(*,*)' Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 97 else write(*,*)' Exiting program. ' call qmexit endif 98 continue c c Determine orbital number. c if(iorbc.eq.1)then iorb = io(1) elseif(iorbc.eq.2)then iorb = io(1)*10 + io(2) else write(*,1020)nmdets write(*,*)' Retry? < y or n > ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 97 else write(*,*)' Exiting program. ' call qmexit endif endif c c Turn this spin orbital on. c if(alpha)then isp = 1 elseif(beta)then isp = 2 endif isp = isp + 2*(invord(iorb)-icl-1) ldets(nmdets) = ldets(nmdets) + n2pow(isp) iorbc = 0 99 continue 3 continue write(*,*)' Would you like to input another determinant? ' write(*,*)' < y or n> ' read(*,'(a)')ans if(ans.eq.'Y'.or.ans.eq.'y')then goto 97 endif return 1000 format(' Input determinant: ',i3) 1010 format(' Unrecognizable character: ',a1,' in ',/, & ' determinant: ',i3) 1020 format(' Orbital number too large in determinant ',i3) end subroutine formh(h,g,hami,ldets,n2pow,nmdets,nmdets2, & n,n2,n4,icl2,iva2) c c This subroutine sets up the Hamiltonian matrix. c implicit real(a-h,o-z) include "limits.h" dimension h(n2), g(n4), hami(nmdets2), ldets(nmdets), n2pow(nbits) data zilch/1.e-4/ data factor/1.442695040888963e0/ isym2(i,j)=(max0(i,j)*(max0(i,j)-1))/2 + min0(i,j) isym4(i,j,k,l)=(max0(isym2(i,j),isym2(k,l))* & (max0(isym2(i,j),isym2(k,l))-1))/2 + & min0(isym2(i,j),isym2(k,l)) nbit1(num)=log(float(num))*factor+zilch nbit2(num)=log(float(mod(num,2**nbit1(num))))*factor+zilch nspnnc(n)=n+1 nspnwc(n)=n+icl2+1 nsymnc(n)=(nspnnc(n)+1)/2 nsymwc(n)=(nspnwc(n)+1)/2 spinij(i,j)=float(iabs(((-1)**i+(-1)**j)/2)) occ(ndet,nso)=float((and(ndet,n2pow(nso)))/n2pow(nso)) nl = (nmdets*(nmdets+1))/2 call ezero(hami,nl) ij=0 do 60 i=1,nmdets ndeti=ldets(i) do 55 j=1,i ij=ij+1 ndetj=ldets(j) c c Determine particle difference c ndiff=ibtsum(xor(ndeti,ndetj)) c write(*,*)' ndeti, ndetj, ndiff: ',ndeti,ndetj,ndiff if(ndiff.le.4)then if(ndiff.eq.4)then c c Get k+,l+,i,j bit numbers c ijbit=and(ndeti,xor(ndeti,ndetj)) klbit=and(ndetj,xor(ndeti,ndetj)) ibit=nbit1(ijbit) jbit=nbit2(ijbit) kbit=nbit1(klbit) lbit=nbit2(klbit) c write(*,*)' ijbit,klbit,ibit,jbit,kbit,lbit: ', c & ijbit,klbit,ibit,jbit,kbit,lbit c c Get sign. c inum=0 ickocc=n2pow(nspnnc(ibit))-n2pow(nspnnc(jbit)+1) if(ickocc.gt.0)then inum=inum+ibtsum(and(ndeti,ickocc)) endif ickocc=n2pow(nspnnc(kbit))-n2pow(nspnnc(lbit)+1) if(ickocc.gt.0)then ndetm=ndeti-n2pow(nspnnc(jbit))-n2pow(nspnnc(ibit)) inum=inum+ibtsum(and(ndetm,ickocc)) endif sign=(-1.0e0)**inum c c Get coupling elements ... x. c x=g(isym4(nsymwc(ibit),nsymwc(kbit), & nsymwc(jbit),nsymwc(lbit)))* & spinij(nspnnc(ibit),nspnnc(kbit))* & spinij(nspnnc(jbit),nspnnc(lbit))- & g(isym4(nsymwc(ibit),nsymwc(lbit), & nsymwc(jbit),nsymwc(kbit)))* & spinij(nspnnc(ibit),nspnnc(lbit))* & spinij(nspnnc(jbit),nspnnc(kbit)) hami(ij)=sign*x c elseif(ndiff.eq.2)then c c Get j+,i bit numbers. c ibit=nbit1(and(ndeti,xor(ndeti,ndetj))) jbit=nbit1(and(ndetj,xor(ndeti,ndetj))) c write(*,*)' ibit,jbit: ', c & ibit,jbit c c Get sign. c sign=1.0e0 i1=min0(nspnnc(ibit),nspnnc(jbit)) i2=max0(nspnnc(ibit),nspnnc(jbit)) ickocc=n2pow(i2)-n2pow(i1+1) if(ickocc.gt.0)then inum=ibtsum(and(ndeti,ickocc)) sign=(-1.0e0)**inum endif c c Get coupling elements ... x. c ncomm=and(ndeti,ndetj) x=h(isym2(nsymwc(jbit),nsymwc(ibit))) do 5 kbit=0,icl2-1 x=x+g(isym4(nsymwc(ibit),nsymwc(jbit), & nsymnc(kbit),nsymnc(kbit)))* & spinij(nspnnc(ibit),nspnnc(jbit))- & g(isym4(nsymwc(ibit),nsymnc(kbit), & nsymwc(jbit),nsymnc(kbit)))* & spinij(nspnnc(ibit),nspnnc(kbit))* & spinij(nspnnc(jbit),nspnnc(kbit)) 5 continue kind=0 do 10 kbit=icl2,icl2+iva2-1 kind=kind+1 x=x+(g(isym4(nsymwc(ibit),nsymwc(jbit), & nsymnc(kbit),nsymnc(kbit)))* & spinij(nspnnc(ibit),nspnnc(jbit))- & g(isym4(nsymwc(ibit),nsymnc(kbit), & nsymwc(jbit),nsymnc(kbit)))* & spinij(nspnnc(ibit),nspnnc(kbit))* & spinij(nspnnc(jbit),nspnnc(kbit)))* & occ(ncomm,kind) 10 continue hami(ij)=sign*x c elseif(ndiff.eq.0)then c c Get coupling elements ... x. c ncomm=and(ndeti,ndetj) x=0.0e0 do 15 kbit=0,icl2-1 x=x+h(isym2(nsymnc(kbit),nsymnc(kbit))) 15 continue kind=0 do 20 kbit=icl2,icl2+iva2-1 kind=kind+1 x=x+h(isym2(nsymnc(kbit),nsymnc(kbit)))* & occ(ncomm,kind) 20 continue do 30 lbit=0,icl2-1 do 25 kbit=0,lbit x=x+g(isym4(nsymnc(lbit),nsymnc(lbit), & nsymnc(kbit),nsymnc(kbit)))- & g(isym4(nsymnc(lbit),nsymnc(kbit), & nsymnc(kbit),nsymnc(lbit)))* & spinij(nspnnc(lbit),nspnnc(kbit)) 25 continue 30 continue lind=0 do 40 lbit=icl2,icl2+iva2-1 lind=lind+1 do 35 kbit=0,icl2-1 x=x+(g(isym4(nsymnc(lbit),nsymnc(lbit), & nsymnc(kbit),nsymnc(kbit)))- & g(isym4(nsymnc(lbit),nsymnc(kbit), & nsymnc(kbit),nsymnc(lbit)))* & spinij(nspnnc(lbit),nspnnc(kbit)))* & occ(ncomm,lind) 35 continue 40 continue lind=0 do 50 lbit=icl2,icl2+iva2-1 lind=lind+1 kind=0 do 45 kbit=icl2,lbit kind=kind+1 x=x+(g(isym4(nsymnc(lbit),nsymnc(lbit), & nsymnc(kbit),nsymnc(kbit)))- & g(isym4(nsymnc(lbit),nsymnc(kbit), & nsymnc(kbit),nsymnc(lbit)))* & spinij(nspnnc(lbit),nspnnc(kbit)))* & occ(ncomm,lind)*occ(ncomm,kind) 45 continue 50 continue hami(ij)=x c else c hami(ij)=0.0e0 c endif endif 55 continue 60 continue return end integer function ibtsum(iword) implicit real (a-h,o-z) include "limits.h" jword=iword j=0 do 10 i=0,nbits-1 c j=j+and(ishft(jword,-i),1) j=j+and(rshift(jword,i),1) 10 continue ibtsum=j return end subroutine putdet(ldets,nmdets,newdet) implicit real(a-h,o-z) include "limits.h" dimension ldets(maxdet) nmdets = nmdets + 1 if(nmdets.gt.maxdet)then write(*,1000) call qmexit endif ldets(nmdets)=newdet return 1000 format(/,' MAXDET is too small, CASGEN is generating ', & ' too many determinants. ',/, & ' Either increase MAXDET or decrease the ',/, & ' number of active orbitals. ') end subroutine wtdets(ldets,n2pow,iord,invord,nmdets,icl,iva,nelact) implicit real (a-h,o-z) character*1 tdets,ones,alpha,beta,blank,pipe include "limits.h" dimension ldets(maxdet) dimension n2pow(nbits) dimension iord(maxorb), invord(maxorb) dimension tdets(80) dimension ones(10) data alpha/'A'/ data beta/'B'/ data blank/' '/ data pipe/'|'/ data ones/'1','2','3','4','5','6','7','8','9','0'/ write(*,*)' Note, A = alpha spin and B = beta spin. ' write(*,*) c c Put bit packed determinants in character representation. c do 15 k=1,nmdets do 5 j = 1,80 tdets(j) = blank 5 continue tdets(1) = blank tdets(2) = pipe tdets(3) = blank icntr = 3 lrel = 0 do 10 l=1,iva*2,2 lrel = lrel + 1 labs = iord(icl+lrel) if(and(ldets(k),n2pow(l)).ne.0)then icntr = icntr + 1 if(labs.lt.10)then tdets(icntr)=ones(labs) icntr = icntr + 1 elseif(labs.gt.10)then tdets(icntr)=ones(1) icntr = icntr + 1 tdets(icntr)=ones(labs-10) icntr = icntr + 1 else tdets(icntr)=ones(1) icntr = icntr + 1 tdets(icntr)=ones(10) icntr = icntr + 1 endif tdets(icntr)= alpha icntr = icntr + 1 tdets(icntr)= blank endif if(and(ldets(k),n2pow(l+1)).ne.0)then icntr = icntr + 1 if(labs.lt.10)then tdets(icntr)=ones(labs) icntr = icntr + 1 elseif(labs.gt.10)then tdets(icntr)=ones(1) icntr = icntr + 1 tdets(icntr)=ones(labs-10) icntr = icntr + 1 else tdets(icntr)=ones(1) icntr = icntr + 1 tdets(icntr)=ones(10) icntr = icntr + 1 endif tdets(icntr)= beta icntr = icntr + 1 tdets(icntr)= blank endif 10 continue icntr = icntr + 1 tdets(icntr) = pipe icntr = icntr + 1 tdets(icntr) = blank write(*,*)k,(tdets(j),j=1,icntr) 15 continue return end subroutine get1es(h,iord,invord,n,n2) implicit real (a-h,o-z) include "limits.h" dimension h(n2) dimension iord(maxorb), invord(maxorb) isym2(i,j)=(max0(i,j)*(max0(i,j)-1))/2+min0(i,j) c c Read in the one-electron integrals. c rewind(unit=30) do 20 i = 1, n ii = invord(i) do 10 j = 1, i jj = invord(j) read(30,*)xij iijj = isym2(ii,jj) h(iijj) = xij 10 continue 20 continue return end subroutine get2es(g,iord,invord,n,n4) implicit real (a-h,o-z) include "limits.h" dimension g(n4) dimension iord(maxorb), invord(maxorb) isym2(i,j)=(max0(i,j)*(max0(i,j)-1))/2 + min0(i,j) isym4(i,j,k,l)=(max0(isym2(i,j),isym2(k,l))* & (max0(isym2(i,j),isym2(k,l))-1))/2 + & min0(isym2(i,j),isym2(k,l)) c c Read in a canonical list of two-electron integrals c from unit 30 skipping past the one-electron integrals. c c Skip past the one-electron integrals. c rewind(unit=30) il = (n*(n+1))/2 do 10 i = 1,il read(30,*)g(i) 10 continue c c Read in the two-electron integrals. c do 50 i = 1,n ii = invord(i) do 40 k = 1,i kk = invord(k) do 30 j = 1,i jj = invord(j) lend = j if(i.eq.j)lend = k do 20 l = 1,lend ll = invord(l) ijkl = isym4(ii,kk,jj,ll) read(30,*)xijkl g(ijkl) = xijkl 20 continue 30 continue 40 continue 50 continue return end