1C> \file eij.F
2C> Coupling coefficients eij
3C>
4C> \ingroup selci
5C> @{
6C>
7C> \brief Coupling coefficients involving single excitations
8C>
9C> Calculate coupling coefficients [1,2]
10C> \f{eqnarray*}{
11C>    e(u,v) &=& <Iu|E_{ij}|Jv>
12C> \f}
13C> where u, v label the spin functions
14C> and the orbital occupation of I is specified by ns, iocc
15C> and indbar. Only for i.ne.j.
16C>
17C> \f{eqnarray*}{
18C>    <Iu|Eij|Jv> &=& \sum_q <Iu|E_{ia}|Sq><Sq|E_{aj}|Jv>
19C> \f}
20C>
21C> possible cases are:
22C>
23C> iocc(i) | iocc(j) | case
24C> --------|---------|-----
25C> 1       |  0      | 1
26C> 1       |  1      | 2
27C> 3       |  0      | 3
28C> 3       |  1      | 4
29C>
30C> ### References ###
31C>
32C> [1] P.J. Knowles, H.-J. Werner, "An efficient method for the
33C>     evaluation of coupling coefficients in configuration interaction
34C>     calculations", Chem. Phys. Lett. <b>145</b> (1988) 514-522, doi:
35C>     <a href="https://doi.org/10.1016/0009-2614(88)87412-8">
36C>     10.1016/0009-2614(88)87412-8</a>.
37C>
38C> [2] P.E.M. Siegbahn, "A new direct CI method for large CI expansions
39C>     in a small orbital space", Chem. Phys. Lett. <b>109</b> (1984)
40C>     417-423, doi:
41C>     <a href="https://doi.org/10.1016/0009-2614(84)80336-X">
42C>     10.1016/0009-2614(84)80336-X</a>.
43C>
44      subroutine selci_eij(e,i,j,ns,indbar,iocc,w1,w2)
45*
46* $Id$
47*
48#include "implicit.fh"
49#include "errquit.fh"
50#include "ciinfo.fh"
51#include "ceij.fh"
52c
53      dimension e(*),indbar(*),iocc(*),
54     $     w1(nfmax,nfmax,nsmax),w2(nfmax2,nfmax,nsmax-1)
55c
56c     e(u,v) = <Iu|Eij|Jv>, where u,v label the spin functions
57c     and the orbital occupation of I is specified by ns, iocc
58c     and indbar. Only for i.ne.j.
59c
60c     <Iu|Eij|Jv> = sum(q) <Iu|Eia|Sq><Sq|Eaj|Jv>
61c
62c     possible cases are:
63c
64c     iocc(i)  iocc(j)  case
65c     1         0       1
66c     1         1       2
67c     3         0       3
68c     3         1       4
69c
70      if (i.eq.j) call errquit('eij, i=j',i, UNKNOWN_ERR)
71c
72      numf = nf(ns)
73      icase = iocc(i)+iocc(j)
74      goto (10,20,30,40) icase
75      call errquit('eij: fell thru goto',icase, UNKNOWN_ERR)
76c
77c     ni=1, nj=0 : e(u,v) = sum(q) w1(u,q,ib)*w1(v,q,jb)
78c
79 10   if (ns.le.nseij) then
80         nsdiff = nseij - ns
81         ib = indbar(i) + nsdiff
82         jb = indbar(j) + nsdiff
83         if (numf.eq.1) then
84            e(1) = case1(1,1,ib,jb)
85         else
86            ie = 0
87cvd$ shortloop
88            do 11 is = 1,numf
89cvd$ shortloop
90               do 12 js = 1,numf
91                  e(ie+js) = case1(js,is,ib,jb)
92 12            continue
93               ie = ie + numf
94 11         continue
95         endif
96      else
97         nsdiff = nsmax - ns
98         ib = indbar(i) + nsdiff
99         jb = indbar(j) + nsdiff
100         if (j.gt.i) jb = jb - 1
101c     call selci_mxma(w1(1,1,ib),1,nfmax,w1(1,1,jb),nfmax,1,e,1,numf,
102c     $     numf,numf,numf)
103         call selci_axbt(w1(1,1,ib),nfmax, w1(1,1,jb),nfmax, e,numf,
104     $        numf, numf, numf)
105      endif
106      return
107c
108c     ni=1, nj=1 : e(u,v) = sum(q) w1(u,q,ib)*w2(v,q,jb)
109c
110 20   numf2 = nf(ns-2)
111      if (ns.le.nseij) then
112         nsdiff = nseij - ns
113         ib = indbar(i) + nsdiff
114         jb = indbar(j) + nsdiff
115         ie = 0
116cvd$  shortloop
117         do 22 js = 1,numf2
118cvd$  shortloop
119            do 23 is = 1,numf
120               e(ie+is) = case2(is,js,ib,jb)
121 23         continue
122            ie = ie + numf
123 22      continue
124      else
125         nsdiff = nsmax - ns
126         ib = indbar(i) + nsdiff
127         jb = indbar(j) + nsdiff
128         if (j.gt.i) jb = jb - 1
129c     call selci_mxma(w1(1,1,ib),1,nfmax,w2(1,1,jb),nfmax2,1,e,1,numf,
130c     $     numf,numf,numf2)
131         call selci_axbt(w1(1,1,ib),nfmax, w2(1,1,jb),nfmax2, e,numf,
132     $        numf, numf, numf2)
133      endif
134      return
135c
136c     ni=2, nj=0 : e(u,v) = sum(q) w2(u,q,ib)*w1(v,q,jb)
137c
138 30   numf2 = nf(ns+2)
139      if (ns+2.le.nseij) then
140         nsdiff = nseij - ns - 2
141         ib = indbar(i) + nsdiff
142         jb = indbar(j) + nsdiff
143         if (i.gt.j) then
144            ib = ib + 1
145         else
146            jb = jb + 1
147         endif
148         ie = 1
149cvd$  shortloop
150         do 32 is = 1,numf
151            iee = ie
152cvd$  shortloop
153            do 33 js = 1,numf2
154               e(iee) = case2(js,is,jb,ib)
155               iee = iee + numf
156 33         continue
157            ie = ie + 1
158 32      continue
159      else
160c
161         nsdiff = nsmax - ns - 2
162         ib = indbar(i) + nsdiff
163         jb = indbar(j) + nsdiff
164         if (j.gt.i) jb = jb + 1
165c     call selci_mxma(w2(1,1,ib),1,nfmax2,w1(1,1,jb),nfmax,1,e,1,numf,
166c     $     numf,numf2,numf2)
167         call selci_axbt(w2(1,1,ib),nfmax2, w1(1,1,jb),nfmax, e,numf,
168     $        numf, numf2, numf2)
169      endif
170      return
171c
172c     ni=2, nj=1 : e(u,v) = sum(q) w2(u,q,ib)*w2(v,q,jb)
173c
174 40   continue
175      if (ns.le.ns4eij) then
176         nsdiff = ns4eij - ns
177         ib = indbar(i) + nsdiff
178         jb = indbar(j) + nsdiff
179         if (numf.eq.1) then
180            e(1) = case4(1,1,ib,jb)
181         else
182            ie = 0
183cvd$  shortloop
184            do 42 is = 1,numf
185cvd$  shortloop
186               do 43 js =1,numf
187                  e(ie+js) = case4(js,is,ib,jb)
188 43            continue
189               ie = ie + numf
190 42         continue
191         endif
192      else
193         numf2 = nf(ns+2)
194         nsdiff = nsmax - ns - 2
195         ib = indbar(i) + nsdiff
196         jb = indbar(j) + nsdiff
197         if (j.gt.i) jb = jb + 1
198c     call selci_mxma(w2(1,1,ib),1,nfmax2,w2(1,1,jb),nfmax2,1,e,1,numf,
199c     $     numf,numf2,numf)
200         call selci_axbt(w2(1,1,ib),nfmax2, w2(1,1,jb),nfmax2, e,numf,
201     $        numf, numf2, numf)
202      endif
203      return
204c
205      end
206C>
207C> @}
208