1c $Id$
2*
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Computes 3 center 2-electron integrals
7C>
8C> Computes 3 center 2-electron integrals of the following kind:
9C> \f{eqnarray*}{
10C> ({\mu}|{\nu}{\lambda}) = \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\frac{1}{r_{12}}
11C> g_{\nu}(X_{\nu},r_{2})g_{\lambda}(X_{\lambda},r_{2})dr_{1}dr_{2}
12C> \f}
13C>
14c:tex-% this is part of the API Standard Integral routines.
15c:tex-\subsection{int\_2e3c}
16c:tex-this routine computes the 3 center 2 electron integrals:
17c:tex-\begin{eqnarray*}
18c:tex-({\mu}|{\nu}{\lambda}) = \int_{-\infty}^{\infty} g_{\mu}(X_{\mu},r_{1})\frac{1}{r_{12}}
19c:tex-g_{\nu}(X_{\nu},r_{2})g_{\lambda}(X_{\lambda},r_{2})dr_{1}dr_{2}
20c:tex-\end{eqnarray*}
21c:tex-
22c:tex-{\it Syntax:}
23c:tex-\begin{verbatim}
24      subroutine int_2e3c(brain, ish, ketin, jsh, ksh,
25     &       lscr, scr, leri, eri)
26c:tex-\end{verbatim}
27      implicit none
28c
29c basic api routine to generate a block of 3 center two electron integrals
30c eri = <bra_g(ish)|ket_g(jsh).ket_g(ksh)>
31c
32#include "apiP.fh"
33#include "errquit.fh"
34#include "bas.fh"
35#include "nwc_const.fh"
36#include "basP.fh"
37#include "basdeclsP.fh"
38#include "geomP.fh"
39#include "geobasmapP.fh"
40#include "mafdecls.fh"
41#include "bas_exndcf_dec.fh"
42#include "bas_ibs_dec.fh"
43#include "int_nbf.fh"
44#include "stdio.fh"
45#include "rel_nwc.fh"
46#include "util.fh"
47c
48c::external subroutines used
49c errquit
50c::functions
51      logical cando_nw
52      logical cando_sp
53      logical cando_sim
54      logical int_chk_sh
55      logical int_chk_init
56      external cando_nw
57      external cando_sp
58      external cando_sim
59      external int_chk_sh
60      external int_chk_init
61*----------------------------#define USE_TEXAS
62cedo#if defined(USE_TEXAS_BROKE)
63      logical cando_txs
64      external cando_txs
65cedo#endif
66c::functions
67      integer int_nint_cart
68      external int_nint_cart
69c:: passed
70c:tex-\begin{verbatim}
71      integer brain !< [Input] bra basis set handle
72      integer ish   !< [Input] shell/contraction index
73      integer ketin !< [Input] ket basis set handle
74      integer jsh   !< [Input] shell/contraction index
75      integer ksh   !< [Input] shell/contraction index
76      integer lscr  !< [Input] length of scratch array
77      double precision scr(lscr) !< [Scratch] array
78      integer leri  !< [Input] length of integral array
79      double precision eri(leri) !< [Output] 2e3c integrals
80c:tex-\end{verbatim}
81c:: local
82      logical shells_ok
83      integer bra, ket
84      integer p_geom, cd_geom, ucont
85      integer Lp, p_prim, p_gen, p_iexp, p_icfp, p_cent
86      integer Lc, c_prim, c_gen, c_iexp, c_icfp, c_cent
87      integer Ld, d_prim, d_gen, d_iexp, d_icfp, d_cent
88      logical status_nw, status_sp, status_sim
89cedo#if defined(USE_TEXAS_BROKE)
90      logical status_gen
91      logical status_txs
92      integer num_quart, dummy_lab(2)
93      double precision roff(3), q4
94      integer txs_i,txs_j,txs_k,txs_d
95      logical dum_log
96cedo#endif
97      integer nint
98      logical OFALSE
99      logical any_spherical
100      logical used_txs
101c.rel-dmd
102      logical status_rel, ket_rel, j_rel, k_rel
103      integer sbas, abas, kets, c_icfps, d_icfps
104c
105      integer WarnP
106      save WarnP
107      data WarnP /0/
108c
109#include "bas_exndcf_sfn.fh"
110#include "bas_ibs_sfn.fh"
111c
112      OFALSE = .false.
113      used_txs=.false.
114c
115c check if ERI is big enough - TLW
116c
117      nint = int_nint_cart(brain,ish,ketin,jsh,ketin,ksh,ketin,0)
118      if (nint.gt.leri) then
119        write(luout,*) 'nint = ',nint
120        write(luout,*) 'leri     = ',leri
121        call errquit('int_2e3c: nint>leri error',911, INT_ERR)
122      endiF
123c
124c check initialization
125c
126      if (.not.int_chk_init('int_2e3c'))
127     &       call errquit('int_2e3c: int_init was not called' ,0,
128     &       INT_ERR)
129c
130c     check input shell ranges
131c
132      shells_ok = int_chk_sh(brain,ish)
133      shells_ok = shells_ok .and. int_chk_sh(ketin,jsh)
134      shells_ok = shells_ok .and. int_chk_sh(ketin,ksh)
135      if (.not.shells_ok)
136     &       call errquit('int_2e3c: invalid contraction/shell',0,
137     &       BASIS_ERR)
138cedo#ifndef USE_TEXAS_BROKE
139      call int_nogencont_check(brain,'int_2e3c:bra')
140      call int_nogencont_check(ketin,'int_2e3c:ket')
141chvd  I currently find no evidence that the code cannot handle SP-shells
142c     call int_nospshell_check(brain,'int_2e3c:bra')
143c     call int_nospshell_check(ketin,'int_2e3c:ket')
144cedo#endif
145c
146c     define center information required
147c
148      bra = brain + BASIS_HANDLE_OFFSET
149      ket = ketin + BASIS_HANDLE_OFFSET
150      kets = ket
151      p_geom  = ibs_geom(bra)
152      cd_geom = ibs_geom(ket)
153c
154c  check if spherical
155c
156      any_spherical = bas_spherical(bra).or.bas_spherical(ket)
157c
158      if (p_geom.ne.cd_geom.and.WarnP.eq.0) then
159        write(luout,*)
160     &      'int_2e3c: WARNING: possible geometry inconsistency'
161        write(luout,*)'bra geometry handle:',p_geom
162        write(luout,*)'ket geometry handle:',cd_geom
163        WarnP = 1
164      endif
165c
166      p_cent  = (sf_ibs_cn2ce(ish,bra))
167      c_cent  = (sf_ibs_cn2ce(jsh,ket))
168      d_cent  = (sf_ibs_cn2ce(ksh,ket))
169c
170      ucont   = (sf_ibs_cn2ucn(ish,bra))
171      Lp      = infbs_cont(CONT_TYPE ,ucont,bra)
172      p_gen   = infbs_cont(CONT_NGEN ,ucont,bra)
173c
174      ucont   = (sf_ibs_cn2ucn(jsh,ket))
175      Lc      = infbs_cont(CONT_TYPE ,ucont,ket)
176      c_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
177c
178      ucont   = (sf_ibs_cn2ucn(ksh,ket))
179      Ld      = infbs_cont(CONT_TYPE ,ucont,ket)
180      d_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
181
182*
183c
184c set integral code status
185c
186      status_nw = cando_nw(brain,ish,0).and.cando_nw(ketin,jsh,ksh)
187      status_sp = cando_sp(brain,ish,0).and.cando_sp(ketin,jsh,ksh)
188      status_sim = cando_sim(brain,ish,0).and.cando_sim(ketin,jsh,ksh)
189cedo#if defined(USE_TEXAS_BROKE)
190      status_txs = cando_txs(brain,ish,0).and.cando_txs(ketin,jsh,ksh)
191      status_gen = (max(p_gen,c_gen,d_gen)) .gt. 1  ! if general contraction texas is only option
192
193*                                (p|ff) or greater do texas
194      if (.not.status_gen)
195     &    status_txs = status_txs .and. ((abs(Lp)+abs(Lc)+abs(Ld)).ge.7)
196cedo#endif
197*
198      status_rel = dyall_mod_dir .and. .not.nesc_1e_approx
199     &    .and. (ketin .eq. ao_bsh)
200      k_rel=.false.
201      j_rel=.false.
202      if (status_rel) then
203c
204c     get basis set handles; relativistic integral option valid
205c     if bra or ket are the ao basis and bra and ket have both
206c     functions relativistic
207c
208        ket_rel = .false.
209        sbas = sc_bsh + BASIS_HANDLE_OFFSET
210        abas = ao_bsh + BASIS_HANDLE_OFFSET
211        kets = sbas
212        ket_rel = ket .eq. abas
213        if (ket_rel) then
214          ucont = sf_ibs_cn2ucn(jsh,ket)
215          j_rel = infbs_cont(CONT_RELLS ,ucont,ket) .ne. 0
216          ucont = sf_ibs_cn2ucn(ksh,ket)
217          k_rel = infbs_cont(CONT_RELLS ,ucont,ket) .ne. 0
218          ket_rel = ket_rel .and. j_rel .and. k_rel
219        end if
220        status_rel = status_rel .and. ket_rel
221      end if
222c
223      if (status_sp .and. .not.status_rel) then
224        call genr70(
225     &         brain,ish,coords(1,p_cent,p_geom),
226     &                 0,coords(1,p_cent,p_geom),
227     &         ketin,jsh,coords(1,c_cent,cd_geom),
228     &               ksh,coords(1,d_cent,cd_geom),
229     &         eri)
230cedo#if defined(USE_TEXAS_BROKE)
231      else if (status_txs .and. .not.status_rel) then
232        num_quart = 1
233        dummy_lab(1) = 0
234        dummy_lab(2) = 0
235        roff(1) = 0.0d00
236        roff(2) = 0.0d00
237        roff(3) = 0.0d00
238        txs_i = ish
239        txs_j = jsh
240        txs_k = ksh
241        txs_d = 0
242        dum_log = OFALSE
243        q4 = 1.0d00
244        nint = 0
245        call texas_hf2_m(
246     &        brain,txs_i,txs_d,
247     &        ketin,txs_j,txs_k,num_quart,
248     &        q4,OFALSE,roff,roff,roff,roff,OFALSE,
249     &        eri, leri, dummy_lab, dummy_lab, dummy_lab, dummy_lab,
250     &        nint, OFALSE, dum_log, scr, lscr, 0.0d0,'scfd_int')
251        used_txs=.true.
252cedo#endif
253      else if(status_nw) then
254c
255        ucont   = (sf_ibs_cn2ucn(ish,bra))
256        Lp      = infbs_cont(CONT_TYPE ,ucont,bra)
257        p_prim  = infbs_cont(CONT_NPRIM,ucont,bra)
258        p_gen   = infbs_cont(CONT_NGEN ,ucont,bra)
259        p_iexp  = infbs_cont(CONT_IEXP ,ucont,bra)
260        p_icfp  = infbs_cont(CONT_ICFP ,ucont,bra)
261c
262        ucont   = (sf_ibs_cn2ucn(jsh,ket))
263        Lc      = infbs_cont(CONT_TYPE ,ucont,ket)
264        c_prim  = infbs_cont(CONT_NPRIM,ucont,ket)
265        c_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
266        c_iexp  = infbs_cont(CONT_IEXP ,ucont,ket)
267        c_icfp  = infbs_cont(CONT_ICFP ,ucont,ket)
268        if (j_rel) ucont = ao_to_ls(ucont)
269        c_icfps = infbs_cont(CONT_ICFP ,ucont,kets)
270c
271        ucont   = (sf_ibs_cn2ucn(ksh,ket))
272        Ld      = infbs_cont(CONT_TYPE ,ucont,ket)
273        d_prim  = infbs_cont(CONT_NPRIM,ucont,ket)
274        d_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
275        d_iexp  = infbs_cont(CONT_IEXP ,ucont,ket)
276        d_icfp  = infbs_cont(CONT_ICFP ,ucont,ket)
277        if (k_rel) ucont = ao_to_ls(ucont)
278        d_icfps = infbs_cont(CONT_ICFP ,ucont,kets)
279c
280        if (status_rel) then
281          call rel_2e4c_sf (
282     &        coords(1,p_cent,p_geom), dbl_mb(mb_exndcf(p_iexp,bra)),
283     &        dbl_mb(mb_exndcf(p_icfp,bra)),
284     &        dbl_mb(mb_exndcf(p_icfp,bra)),p_prim,p_gen,Lp,p_cent,
285     &        coords(1,p_cent,p_geom),DCexp,DCcoeff,
286     &        DCcoeff,1,1,0,p_cent,
287     &        coords(1,c_cent,cd_geom), dbl_mb(mb_exndcf(c_iexp,ket)),
288     &        dbl_mb(mb_exndcf(c_icfp,ket)),
289     &        dbl_mb(mb_exndcf(c_icfps,kets)), c_prim,c_gen,Lc,c_cent,
290     &        coords(1,d_cent,cd_geom), dbl_mb(mb_exndcf(d_iexp,ket)),
291     &        dbl_mb(mb_exndcf(d_icfp,ket)),
292     &        dbl_mb(mb_exndcf(d_icfps,kets)), d_prim,d_gen,Ld,d_cent,
293c...................... canAB   canCD   canPQ   DryRun
294     &        eri,leri,.false.,.false.,.false.,.false.,scr,lscr,
295c............  bra_rel                    do_ssss
296     &        .false.,ket_rel,ss_one_cent,.false.,rel_dbg)
297        else
298          call hf2(
299     &        coords(1,p_cent,p_geom), dbl_mb(mb_exndcf(p_iexp,bra)),
300     &        dbl_mb(mb_exndcf(p_icfp,bra)), p_prim, p_gen, Lp,
301     &        coords(1,p_cent,p_geom), DCexp,
302     &        DCcoeff, 1, 1, 0,
303     &        coords(1,c_cent,cd_geom), dbl_mb(mb_exndcf(c_iexp,ket)),
304     &        dbl_mb(mb_exndcf(c_icfp,ket)), c_prim, c_gen, Lc,
305     &        coords(1,d_cent,cd_geom), dbl_mb(mb_exndcf(d_iexp,ket)),
306     &        dbl_mb(mb_exndcf(d_icfp,ket)), d_prim,d_gen,Ld,
307c......................... canAB    canCD    canPQ
308     &        eri, leri, OFALSE, OFALSE, OFALSE,
309c............ dryrun
310     &        OFALSE, scr, lscr)
311        end if
312c
313      elseif(status_sim .and. .not.status_rel)  then
314        call nwcsim_hf2_3c(
315     &        bra,ish,
316     &        ket,jsh,ksh,
317     &        nint, eri, leri, scr, lscr)
318      else
319        write(luout,*)'int_2e3c: could not do nw integrals'
320        write(luout,*)' brain :',brain
321        write(luout,*)' ketin :',ketin
322        write(luout,*)' ish   :',ish
323        write(luout,*)' jsh   :',jsh
324        write(luout,*)' ksh   :',ksh
325        call errquit('int_2e3c: fatal error ',0, INT_ERR)
326      endif
327      if (any_spherical.and.(.not.used_txs))then
328c ... reset general contractions for sp shells to 1 since they are handled
329c     as a block of 4. Since int_nbf_* arrays are set to the appropriate size.
330          if (Lp.eq.-1) p_gen = 1
331          if (Lc.eq.-1) c_gen = 1
332          if (Ld.eq.-1) d_gen = 1
333          call spcart_3ctran(eri,scr,lscr,
334     &        int_nbf_x(Lp),int_nbf_s(Lp),Lp,p_gen,bas_spherical(bra),
335     &        int_nbf_x(Lc),int_nbf_s(Lc),Lc,c_gen,bas_spherical(ket),
336     &        int_nbf_x(Ld),int_nbf_s(Ld),Ld,d_gen,bas_spherical(ket),
337     &        OFALSE)
338        endif
339      end
340C> @}
341