1C> \ingroup selci
2C> @{
3C>
4      subroutine selci_eirerj(e,i,ir,j,ns,indbar,iocc,w1,w2,work)
5*
6* $Id$
7*
8#include "implicit.fh"
9#include "errquit.fh"
10#include "ciinfo.fh"
11      dimension e(*),indbar(*),iocc(*),work(*),
12     $     w1(nfmax,nfmax,nsmax),w2(nfmax2,nfmax,nsmax-1)
13c
14c     e(u,v) = <Iu|EirErj|Iv>, where u,v label the spin functions
15c     and the orbital occupation of I is specified by ns, iocc
16c     and indbar. Only for i.ne.j.ne.r, and r singly occupied.
17c
18c     e(u,v) = sum(q) <Iu|Eir|Sq><Sq|Erj|Iv>
19c
20c     <Iu|Eir|Sq> = sum(p) <Iu|Eia|Tp><Tp|Ear|Sq>
21c
22c     possible cases are:
23c
24c     iocc(i)  iocc(j)  case
25c       1         0       1
26c       1         1       2
27c       3         0       3
28c       3         1       4
29c
30      if (iocc(ir).ne.1 .or. ir.eq.i .or. ir.eq.j .or. i.eq.j)
31     $     call errquit('eirerj: wrong i,r,j',ir, INPUT_ERR)
32c
33      icase = iocc(i)+iocc(j)
34      goto (10,20,30,40) icase
35      call errquit('eirerj: fell thru goto',icase, INPUT_ERR)
36c
37 10   continue
38c
39c     ni=1, nj=0 : e(u,v) = w1(u,p,ib)*w2(q,p,rb)*w2(q,s,rb)*w1(v,s,jb)
40c                         ir    ->    ra    ->   rr    ->   ra   ->   jr
41c
42      numf = nf(ns)
43      numf2 = nf(ns-2)
44      ib = indbar(i) + (nsmax-ns)
45      irb = indbar(ir) + (nsmax-ns)
46      if (ir.gt.i) irb = irb - 1
47      jb = indbar(j) + (nsmax-ns)
48      if (j.gt.i) jb = jb - 1
49*      call selci_mxma(w1(1,1,ib),1,nfmax,w2(1,1,irb),nfmax2,1,
50*     $     e,1,numf,numf,numf,numf2)
51      call selci_axbt(w1(1,1,ib),nfmax,w2(1,1,irb),nfmax2,
52     $     e,numf,numf,numf,numf2)
53*      call selci_mxma(e,1,numf,w2(1,1,irb),1,nfmax2,
54*     $     work,1,numf,numf,numf2,numf)
55      call selci_axb(e,numf,w2(1,1,irb),nfmax2,
56     $     work,numf,numf,numf2,numf)
57*      call selci_mxma(work,1,numf,w1(1,1,jb),nfmax,1,
58*     $     e,1,numf,numf,numf,numf)
59      call selci_axbt(work,numf,w1(1,1,jb),nfmax,
60     $     e,numf,numf,numf,numf)
61      return
62c
63 20   continue
64c
65c     ni=1, nj=1 : e(u,v) = w1(u,p,ib)*w2(q,p,rb)*w2(q,s,rb)*w2(v,s,jb)
66c                        ijr    ->   jra    ->  rrj   ->   jra  ->   jjr
67c     same as 10 except for last mxma
68      numf = nf(ns)
69      numf2 = nf(ns-2)
70      ib = indbar(i) + (nsmax-ns)
71      irb = indbar(ir) + (nsmax-ns)
72      if (ir.gt.i) irb = irb - 1
73      jb = indbar(j) + (nsmax-ns)
74      if (j.gt.i) jb = jb - 1
75*      call selci_mxma(w1(1,1,ib),1,nfmax,w2(1,1,irb),nfmax2,1,
76*     $     e,1,numf,numf,numf,numf2)
77      call selci_axbt(w1(1,1,ib),nfmax,w2(1,1,irb),nfmax2,
78     $     e,numf,numf,numf,numf2)
79*      call selci_mxma(e,1,numf,w2(1,1,irb),1,nfmax2,
80*     $     work,1,numf,numf,numf2,numf)
81      call selci_axb(e,numf,w2(1,1,irb),nfmax2,
82     $     work,numf,numf,numf2,numf)
83*      call selci_mxma(work,1,numf,w2(1,1,jb),nfmax2,1,
84*     $     e,1,numf,numf,numf,numf2)
85      call selci_axbt(work,numf,w2(1,1,jb),nfmax2,
86     $     e,numf,numf,numf,numf2)
87      return
88c
89 30   continue
90c
91c     ni=2, nj=0 : e(u,v) = w2(u,p,ib)*w2(q,p,rb)*w2(q,s,rb)*w1(v,s,jb)
92c                        iir    ->   ira    ->  rri   ->   ira  ->   ijr
93c
94      numf = nf(ns)
95      numf2 = nf(ns+2)
96      ib = indbar(i) + (nsmax-ns-2)
97      irb = indbar(ir) + (nsmax-ns-2)
98      if (ir.gt.i) irb = irb + 1
99      jb = indbar(j) + (nsmax-ns-2)
100      if (j.gt.i) jb = jb + 1
101*      call selci_mxma(w2(1,1,ib),1,nfmax2,w2(1,1,irb),nfmax2,1,
102*     $     e,1,numf,numf,numf2,numf)
103      call selci_axbt(w2(1,1,ib),nfmax2,w2(1,1,irb),nfmax2,
104     $     e,numf,numf,numf2,numf)
105*      call selci_mxma(e,1,numf,w2(1,1,irb),1,nfmax2,
106*     $     work,1,numf,numf,numf,numf2)
107      call selci_axb(e,numf,w2(1,1,irb),nfmax2,
108     $     work,numf,numf,numf,numf2)
109*      call selci_mxma(work,1,numf,w1(1,1,jb),nfmax,1,
110*     $     e,1,numf,numf,numf2,numf2)
111      call selci_axbt(work,numf,w1(1,1,jb),nfmax,
112     $     e,numf,numf,numf2,numf2)
113      return
114c
115 40   continue
116c
117c     ni=2, nj=1 : e(u,v) = w2(u,p,ib)*w2(q,p,rb)*w2(q,s,rb)*w2(v,s,jb)
118c                        iijr   ->   ijra   ->  rrij  ->   ijra  ->   jjir
119c     same as 30 except for last mxma
120      numf = nf(ns)
121      numf2 = nf(ns+2)
122      ib = indbar(i) + (nsmax-ns-2)
123      irb = indbar(ir) + (nsmax-ns-2)
124      if (ir.gt.i) irb = irb + 1
125      jb = indbar(j) + (nsmax-ns-2)
126      if (j.gt.i) jb = jb + 1
127*      call selci_mxma(w2(1,1,ib),1,nfmax2,w2(1,1,irb),nfmax2,1,
128*     $     e,1,numf,numf,numf2,numf)
129      call selci_axbt(w2(1,1,ib),nfmax2,w2(1,1,irb),nfmax2,
130     $     e,numf,numf,numf2,numf)
131*      call selci_mxma(e,1,numf,w2(1,1,irb),1,nfmax2,
132*     $     work,1,numf,numf,numf,numf2)
133      call selci_axb(e,numf,w2(1,1,irb),nfmax2,
134     $     work,numf,numf,numf,numf2)
135*      call selci_mxma(work,1,numf,w2(1,1,jb),nfmax2,1,
136*     $     e,1,numf,numf,numf2,numf)
137      call selci_axbt(work,numf,w2(1,1,jb),nfmax2,
138     $     e,numf,numf,numf2,numf)
139      return
140c
141      end
142C>
143C> @}
144