1C> \file makeh.F
2C> The Hamiltonian generator
3C>
4C> \ingroup selci
5C> @{
6C>
7C> \brief Compute the Hamiltonian matrix
8C>
9C> This routine computes the Hamiltonian matrix using a rather
10C> naive approach.
11C>
12C> For each orbital configuration I as many interactions
13C> as possible are accumulated in hbuf. When this is full
14C> or we start processing another I value this buffer is flushed
15C> into the I/O buffer proper.
16C>
17      subroutine selci_makeh(h,g,int12,int34,w1,w2,ioconf,indxci,hd,
18     $     work1, work2, work3, f)
19*
20* $Id$
21*
22#include "implicit.fh"
23#include "errquit.fh"
24#include "ciinfo.fh"
25#include "cbuff.fh"
26#include "global.fh"
27      dimension h(*),g(*),int12(*),int34(*),
28     &     w1(nfmax,nfmax,nsmax),w2(nfmax2,nfmax,nsmax-1),
29     &     ioconf(nintpo,noconf),indxci(noconf),hd(nci)
30c
31      parameter (lenhbf = 20000)
32      dimension iocc(255),indbar(255), lists(255), listd(255),
33     $     hbuf(lenhbf), jbuf(2,lenhbf)
34c
35      dimension work2(nfmax*nfmax),
36     $     work1(nfmax*nfmax),work3(nfmax*nfmax),
37     $     f(norbs*(norbs+1)/2),ind(4)
38c
39      integer selci_iodiff
40c
41c     naively compute the hamiltonian matrix
42c
43c     For each orbital configuration I as many interactions
44c     as possible are accumulated in hbuf. When this is full
45c     or we start processing another I value this buffer is flushed
46c     into the I/O buffer proper
47c
48c     no. of elements in hbuf
49      ninhb = 0
50      ninjb = 0
51c     no. of elements in I/O buffer
52      nrinb = 0
53      niinb = 0
54      nrec = 0
55      nval = 0
56      me = ga_nodeid()
57      nproc = ga_nnodes()
58      icount = -1
59      call dfill(nci, 0.0d0, hd, 1)
60c
61c     loop through I occupancies
62c
63      do 10 iconf = 1,noconf
64c
65c     get required information on I
66c
67         call selci_upkcon(norbs, iocc, ioconf(1,iconf), nintpo, nbitpi)
68         call selci_mkindb(norbs, iocc, indbar, listd, lists, ns, nd)
69         call selci_makef(f, h, g, int12, int34, iocc, listd, lists,
70     $        ns, nd, .false.)
71         ibase = indxci(iconf)
72c
73         do 20 jconf = iconf,noconf
74            icount = icount + 1
75            if (mod(icount, nproc) .ne. me) goto 20
76c
77c     locate interacting J occupancies
78c     all this can be vectorised at a later date
79c
80            iexcit = selci_iodiff(ioconf(1,iconf), ioconf(1,jconf),
81     $           nintpo)
82            if (iexcit .gt. 4) goto 20
83c
84c     jconf interacts with iconf
85c
86            if (iexcit.eq.0) then
87               call selci_makehd(work1,.false.,
88     $              ns,nd,indbar,iocc,lists,listd,w1,w2,
89     $              work2,work3,f,h,
90     $              g,int12,int34,numf)
91c    zero the upper half and diagonal of work1 to simplify putinb
92               do 33 id = 1,numf
93                  idid = (id-1)*numf
94                  hd(ibase+id) = work1(idid+id)
95                  do 34 jd = id,numf
96                     work1(idid+jd) = 0.0d0
97 34               continue
98 33            continue
99               numf2 = numf
100            else  if (iexcit.eq.2) then
101               call selci_getij(i,j,ioconf(1,iconf),ioconf(1,jconf),
102     $              nintpo,nbitpi,iocc)
103               call selci_makehs(work1,i,j,ns,indbar,iocc,lists,
104     $              w1,w2,work2,work3,f,g,int12,int34,numf,numf2)
105            else if(iexcit.eq.4) then
106               call selci_gtijkl(ind,ioconf(1,iconf),ioconf(1,jconf),
107     $              nintpo,nbitpi,iocc,icase)
108               call selci_makhdb(work1,ind,icase,ns,indbar,iocc,
109     $              w1,w2,work2,work3,g,int12,int34,numf,numf2)
110            else
111               call errquit('strange excitation value ',iexcit,
112     &       INPUT_ERR)
113            endif
114            if (numf*numf2+ninhb.gt.lenhbf)
115     $           call selci_putinb(iflham,indxci(iconf),numf,hbuf,ninhb,
116     $           jbuf,ninjb)
117            call dcopy(numf*numf2,work1,1,hbuf(ninhb+1),1)
118            ninhb = ninhb + numf*numf2
119            ninjb = ninjb + 1
120            jbuf(1,ninjb) = indxci(jconf)
121            jbuf(2,ninjb) = numf2
122 20      continue
123         call selci_putinb(iflham,indxci(iconf),numf,hbuf,ninhb,jbuf,
124     $        ninjb)
125 10   continue
126c
127      ninjb = 0
128      call selci_putinb(iflham,-1,0,hbuf,0,jbuf,ninjb)
129c     if running in parallel need to get all the diags together
130      call ga_dgop (99, hd, nci, '+')
131      if (me .eq. 0) then
132        rewind iflhdg
133        call selci_swrite(iflhdg,hd,nci)
134        close(iflhdg, status='keep')
135      endif
136      call ga_igop(991, nrec, 1, '+')
137      call ga_igop(992, nval, 1, '+')
138      if (me .eq. 0) then
139         write(6,99) nrec,nval
140 99      format(/' total no. of hamiltonian records  ',i9/
141     $        ' total no. of hamiltonian elements ',i9/)
142      endif
143      close(iflham, status='keep')
144c
145      end
146C>
147C> @}
148