program fnrgy c **************** Program FENERGY ***************** c This program is designed to read in the LCAO-MO c coefficient matrix, the one- and two-electron AO c integrals (in Dirac <12|12> convention), and the c Fock orbital energies. Upon transformation of the c one- and two-electron integrals from the AO to the c MO basis, the closed shell Hartree - Fock energy c is calculated in two ways. First, the energy is c calculated with the MO integrals, c c Sum(k) 2* + Sum(k,l) (2* - c ) + ZuZv/Ruv. c c Secondly, the energy is calculated with the Fock c orbital energies and one electron energies in the c MO basis, c c Sum(k) (eps(k) + ) + ZuZv/Ruv. c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real (a-h,o-z) logical yesno character*30 dfname, fname, dscrpt include "limits.h" c c Set core needed to: c 2 square arrays of order maxorb c 1 symmetric array of order maxorb c 2 arrays of dimension: c maxorb*maxorb*((maxorb*(maxorb+1))/2) c total = 2*maxorb**2 + (maxorb*(maxorb+1))/2 c 2*maxorb*maxorb*((maxorb*(maxorb+1))/2) c dimension pcore(2*mxosq + mxo2 + 2*mxosq*mxo2) dimension ipntr(10) itotal = 2*mxosq + mxo2 + 2*mxosq*mxo2 + 10 call modinfo('fenergy ') call shwmem('fenergy ',itotal) 10 continue call goon write(*,*) write(*,*)' Input the number of orbitals in your system: ' read(*,*)iorb if(iorb.gt.maxorb)then write(*,1010)maxorb goto 10 endif iorbsq = iorb*iorb iorb2 = (iorb*(iorb+1))/2 ipntr(1) = 1 ipntr(2) = ipntr(1) + iorbsq ipntr(3) = ipntr(2) + iorb2 ipntr(4) = ipntr(3) + iorbsq ipntr(5) = ipntr(4) + iorbsq*iorb2 ipntr(6) = ipntr(5) + iorbsq*iorb2 c c Open the MO-AO coefficient file. c dscrpt = 'MO-AO coefficients' ids = 18 dfname = 'mocoefs.dat' idf = 11 fname = ' ' ifn = 30 iflnum = 31 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call rdmatd(fname,ifn,dscrpt,ids,pcore(ipntr(1)),iorb, & iorb,iflnum) c c Close MO-AO coefficients file. c close(unit=31) else call rdmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb) endif call wtmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb) c c Open the AO integral file. c dscrpt = '1e- and 2e- AO integrals' ids = 24 dfname = 'ao_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 AO integrals available ' write(*,*)' on disk ... this can be done by running ' write(*,*)' program INTEGRAL. ' call qmexit endif call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)), & pcore(ipntr(4)),pcore(ipntr(5)),iorb,iorbsq,iorb2) call qmexit stop 1010 format(' Maximum number of orbitals is ',I2,'. ',/, & ' Please try again. ') end subroutine driver(c,h,p,q,g,n,nsq,n2) implicit real (a-h,o-z) logical yesno character*30 dfname, fname, dscrpt dimension c(n,n), h(n2), p(n,n), q(nsq*n2), g(nsq*n2) n4 = n2*(n2+1)/2 call get1e(h,n,n2) call get2e(g,n,n4) c c Close 1e- and 2- AO integral file. c close(unit=30) c c Transform the one-electron integrals. c write(*,*)' Begin 1e- integral transformation. ' call t1(c,h,p,n,nsq,n2) write(*,*)' End 1e- integral transformation. ' c c Transform the two-electron integrals. c write(*,*)' Begin 2e- integral transformation. ' c c Do the first quarter transformation. c write(*,*)' Begin 1st quarter transformation. ' call q1(c,q,g,n,nsq,n2) write(*,*)' End 1st quarter transformation. ' c c Do the second quarter transformation. c write(*,*)' Begin 2nd quarter transformation. ' call q2(c,q,g,n,nsq,n2) write(*,*)' End 2nd quarter transformation. ' c c Do the third quarter transformation. c write(*,*)' Begin 3rd quarter transformation. ' call q3(c,q,g,n,nsq,n2) write(*,*)' End 3rd quarter transformation. ' c c Do the fourth quarter transformation. c write(*,*)' Begin 4th quarter transformation. ' call q4(c,q,g,n,nsq,n2) write(*,*)' End 4th quarter transformation. ' write(*,*)' End 2e- integral transformation. ' c c Open the Fock orbital energy file. c dscrpt = 'Fock orbital energies' ids = 21 dfname = 'orbital.evals' idf = 13 fname = ' ' ifn = 30 iflnum = 32 call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call rdmatd(fname,ifn,dscrpt,ids,p,n,1,iflnum) c c Close energies file. c close(unit=32) else call rdmatu(dscrpt,ids,p,n,1) endif call wtmatu(dscrpt,ids,p,n,1) 5 continue c c Determine the number of closed shell orbitals. c write(*,*)' Input the number of doubly occupied orbitals. ' read(*,*)ndocc write(*,1000)ndocc if(ndocc.gt.n)then write(*,*)' The number of doubly occupied orbitals cannot ' write(*,*)' exceed the total number of orbitals. Please ' write(*,*)' try again. ' goto 5 endif 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(*,*)' Input the nuclear repulsion energy. ' read(*,*)repnuc endif c c Calculate energy using MO integrals and c using Fock energies. c fe1 = 0.0e0 fe2 = 0.0e0 do 101 k = 1, ndocc kk = (k*(k+1))/2 fe1 = fe1 + h(kk) fe2 = fe2 + p(k,1) + h(kk) 101 continue fe1 = 2.0e0*fe1 do 103 k = 1, ndocc do 102 l = 1, ndocc iklkl = icanon(k,l,k,l) ikllk = icanon(k,l,l,k) fe1 = fe1 + 2.0e0*g(iklkl) - g(ikllk) 102 continue 103 continue fe1 = fe1 + repnuc fe2 = fe2 + repnuc ediff = fe1 - fe2 write(*,1010)fe1,fe2,ediff return 1000 format(' There is ',I2,' doubly occupied orbital(s). ') 1010 format(' The closed shell Fock energy from: ',/, & ' formula: ',//, & ' Sum(k) 2* + Sum(k,l) (2* - ',/, & ' ) + ZuZv/Ruv ',//, & ' is: ',F15.8,//, & ' from formula: ',//, & ' Sum(k) (eps(k) + ) + ZuZv/Ruv ',//, & ' is: ',F15.8,//, & ' the difference is: ',F15.8) end subroutine q1(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) c c Do the first quarter transformation. c nl = nsq*n2 call ezero(q,nl) nsq = n*n iabsym = 0 icx = 0 do 50 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 40 ib = 1,ia ibm1 = ib - 1 ibsq = ibm1*n iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq do 30 ic = 1,ia icm1 = ic - 1 icsq = icm1*n icsym = ic*icm1/2 idend = ic if(ic.eq.ia)idend = ib do 20 id = 1,idend idm1 = id - 1 idsq = idm1*n icdsym = (icsym + id) icdsnn = (icdsym - 1)*nsq icx = icx + 1 xg = g(icx) do 10 l = 1,n issx = iabsnn + icsq + l q(issx) = q(issx) + xg*c(id,l) if(ic.ne.id)then issx = iabsnn + idsq + l q(issx) = q(issx) + xg*c(ic,l) endif if(iabsym.ne.icdsym)then issx = icdsnn + iasq + l q(issx) = q(issx) + xg*c(ib,l) if(ia.ne.ib)then issx = icdsnn + ibsq + l q(issx) = q(issx) + xg*c(ia,l) endif endif 10 continue 20 continue 30 continue 40 continue 50 continue return end subroutine q2(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(g,nl) c c Do the second quarter transformation. c nsq = n*n iabsym = 0 do 150 ia = 1,n do 140 ib = 1,ia iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq do 130 ic = 1,n icm1 = ic - 1 icsq = icm1*n do 120 l = 1,n issx = iabsnn + icsq + l xq = q(issx) do 110 k = l,n km1 = k - 1 ksq = km1*n issx = iabsnn + ksq + l g(issx) = g(issx) + xq*c(ic,k) 110 continue 120 continue 130 continue 140 continue 150 continue return end subroutine q3(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(q,nl) c c Do the third quarter transformation. c nsq = n*n klsym = 0 do 250 k = 1,n km1 = k - 1 ksq = km1*n do 240 l = 1,k klsym = klsym + 1 klsnn = (klsym - 1)*nsq iabsym = 0 do 230 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 220 ib = 1,ia ibm1 = ib - 1 ibsq = ibm1*n iabsym = iabsym + 1 iabsnn = (iabsym - 1)*nsq issx = iabsnn + ksq + l xg = g(issx) do 210 j = 1,n issx = klsnn + iasq + j q(issx) = q(issx) + xg*c(ib,j) if(ia.ne.ib)then issx = klsnn + ibsq + j q(issx) = q(issx) + xg*c(ia,j) endif 210 continue 220 continue 230 continue 240 continue 250 continue return end subroutine q4(c,q,g,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n), q(nsq*n2), g(nsq*n2) nl = nsq*n2 call ezero(g,nl) c c Do the fourth quarter transformation. c nsq = n*n klsym = 0 do 350 k = 1,n do 340 l = 1,k klsym = klsym + 1 klsnn = (klsym - 1)*nsq do 330 ia = 1,n iam1 = ia - 1 iasq = iam1*n do 320 j = 1,n issx = klsnn + iasq + j xq = q(issx) istart = max(j,k) if(l.gt.j)istart = k + 1 do 310 i = istart,n icx = icanon(i,k,j,l) g(icx) = g(icx) + xq*c(ia,i) 310 continue 320 continue 330 continue 340 continue 350 continue return end subroutine t1(c,h,p,n,nsq,n2) implicit real (a-h,o-z) dimension c(n,n),h(n2),p(n,n) c c Transform the one-electron integrals. c nl = nsq call ezero(p,nl) iaib = 0 do 430 ia = 1,n do 420 ib = 1,ia iaib = iaib + 1 do 410 j = 1,n p(ia,j) = p(ia,j) + h(iaib)*c(ib,j) if(ia.ne.ib)then p(ib,j) = p(ib,j) + h(iaib)*c(ia,j) endif 410 continue 420 continue 430 continue nl = n2 call ezero(h,nl) do 530 ia = 1,n ij = 0 do 520 j = 1,n DO 510 I = 1,J ij = ij + 1 h(ij) = h(ij) + p(ia,j)*c(ia,i) 510 continue 520 continue 530 continue return end