program fock c ****************** Program FOCK ****************** c This program is designed to read in the LCAO-MO c coefficient matrix, the one- and two-electron AO c integrals and to form a CLOSED SHELL Fock matrix c (i.e., a Fock matrix for species with all doubly c occupied orbitals). c c Copyright (C) 1992 Jeff Nichols and Jack Simons c c ************************************************** implicit real (a-h,o-z) include "limits.h" logical yesno character*30 dfname, fname, dscrpt c c Set core needed to: c 3 square arrays of order maxorb c 1 array of dimension canonical(maxorb**4) c total = 3*mxosq + mxo4 c dimension pcore(3*mxosq + mxo4) dimension ipntr(10) itotal = 3*mxosq + mxo4 + 10 call modinfo('fock ') call shwmem('fock ',itotal) call goon 10 continue write(*,*) write(*,*)' Input the number of orbitals in your system: ' read(*,*)iorb if(iorb.gt.maxorb)then write(*,1010)maxorb goto 10 endif ipntr(1) = 1 ipntr(2) = ipntr(1) + iorb*iorb ipntr(3) = ipntr(2) + iorb*iorb ipntr(4) = ipntr(3) + iorb*iorb iorb2 = iorb*(iorb+1)/2 iorb4 = iorb2*(iorb2+1)/2 ipntr(5) = ipntr(4) + iorb4 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 integral file. c dscrpt = '1e- and 2e- integrals' ids = 21 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)),iorb,iorb2,iorb4) call qmexit stop 1010 format(' Maximum number of orbitals is ',I2,'. ',/, & ' Please try again. ') end subroutine driver(c,p,f,g,n,n2,n4) implicit real (a-h,o-z) logical yesno character*30 dfname, fname, dscrpt dimension c(n,n),p(n,n),f(n,n),g(n4) c c Get the 1e- integrals. c call get1e(g,n,n2) c c Put the 1e- Hamiltonian directly into the Fock matrix. c ij = 0 do 2 i = 1,n do 1 j = 1,i ij = ij + 1 f(i,j) = g(ij) f(j,i) = f(i,j) 1 continue 2 continue c c Get the 2e- integrals. c call get2e(g,n,n4) c c Close 1e- and 2e- integrals file. c close(unit=30) 5 continue c c Determine the number of closed shell orbitals. c write(*,*)' Please input the number of doubly occupied ' write(*,*)' 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 Calculate the charge bond order matrix P(MR,MS). c do 30 mr = 1,n do 20 ms = 1,n p(mr,ms) = 0.0e0 do 10 k = 1,ndocc p(mr,ms) = p(mr,ms) + c(mr,k)*c(ms,k) 10 continue 20 continue 30 continue c c Write out the charge bond order matrix. c write(*,*)' The charge bond order matrix: ' call matout(p,n,n,n) c c Form the CLOSED SHELL Fock matrix. c do 70 mu = 1,n do 60 nu = 1,n do 50 mr = 1,n do 40 ms = 1,n icx1 = icanon(mu,mr,nu,ms) icx2 = icanon(mu,mr,ms,nu) f(mu,nu) = f(mu,nu) + & p(mr,ms)*(2.0e0*g(icx1) - g(icx2)) 40 continue 50 continue 60 continue 70 continue dscrpt = 'Fock matrix' ids = 11 call wtmatu(dscrpt,ids,f,n,n) dfname = 'fock.out' idf = 8 fname = ' ' ifn = 30 iflnum = 32 call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno) if(yesno)then call wtmatd(fname,ifn,dscrpt,ids,f,n,n, & iflnum) c c Close Fock matrix file. c close(unit=32) endif return 1000 format(' There is ',I2,' doubly occupied orbital(s). ') end