1c $Id$
2*
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Calculate the 4-center 2-electron integrals between GIAO basis functions
7C>
8C> The GIAO basis functions are important in the treatment of magnetic
9C> properties. These properties are not automatically gauge invariant but
10C> by multiplying the usual Atomic Orbitals (AO) with a magnetic field dependent
11C> factor Gauge Invariant Atomic Orbitals (GIAO) can be defined. Using these
12C> basis functions the properties can be cast in an gauge invariant form (see
13C> [1,2]).
14C>
15C> In the presence of a magnetic field a vector potential is generated given by
16C> \f{eqnarray*}{
17C>   \vec{A}_A = \vec{B} \otimes \vec{R}_A
18C> \f}
19C> where \f$ \vec{B} \f$ is the magnetic field, \f$ \vec{R}_A \f$ is the
20C> nuclear position relative to the origin of the systems,
21C> \f$ \vec{A} \f$ the generated vector potential.
22C> A GIAO is now given by
23C> \f{eqnarray*}{
24C>    g^{GIAO}_\mu(X_\mu,r) = g_\mu(X_\mu,r)e^{-i \vec{A}_A\cdot\vec{r}}
25C> \f}
26C> The integrals computed here can now be represented as
27C> \f{eqnarray*}{
28C>    ({\mu}{\rho}|{\nu}{\lambda}) = \int_{-\infty}^{\infty}
29C>    g^{GIAO}_{\mu}(X_{\mu},r_{1})g^{GIAO}_{\rho}(X_{\rho},r_{1})\frac{1}{r_{12}}
30C>    g^{GIAO}_{\nu}(X_{\nu},r_{2})g^{GIAO}_{\lambda}(X_{\lambda},r_{2})dr_{1}dr_{2}
31C> \f}
32C>
33C> [1] M. Dupuis,
34C>     "New integral transforms for molecular properties and application
35C>     to a massively parallel GIAO-SCF implementation",
36C>     Comp. Phys. Comm. <b>134</b>, 150-166 (2001), DOI:
37C>     <a href="https://doi.org/10.1016/S0010-4655(00)00195-8">
38C>     10.1016/S0010-4655(00)00195-8</a>.
39C>
40C> [2] K. Wolinski, J.F. Hinton, P. Pulay,
41C>     "Efficient implementation of the gauge-independent atomic orbital
42C>      method for NMR chemical shift calculations",
43C>     J. Am. Chem. Soc. <b>112</b>, 8251-8260 (1990), DOI:
44C>     <a href="https://doi.org/10.1021/ja00179a005">10.1021/ja00179a005</a>.
45C>
46      subroutine int_giao_2e(brain, ish, jsh, ketin, ksh, lsh,
47     &       lscr, scr, leri, eri)
48      implicit none
49c
50c basic api routine to generate a block of two electron integrals
51c eri = <bra_g(ish).bra_g(jsh) | ket_g(ksh).ket_g(lsh)>
52c
53#include "bas.fh"
54#include "errquit.fh"
55#include "nwc_const.fh"
56#include "apiP.fh"
57#include "basP.fh"
58#include "basdeclsP.fh"
59#include "geomP.fh"
60#include "geobasmapP.fh"
61#include "mafdecls.fh"
62#include "bas_exndcf_dec.fh"
63#include "bas_ibs_dec.fh"
64#include "int_nbf.fh"
65#include "stdio.fh"
66#include "rel_nwc.fh"
67#include "util.fh"
68#include "global.fh"
69      common/testdata/timing(20),irepeat
70      double precision timing
71      integer irepeat
72c
73c::external subroutines used
74c errquit
75c::functions
76      logical cando_nw
77      logical cando_sp
78      logical cando_txs
79      external cando_nw
80      external cando_sp
81      external cando_txs
82c:: passed
83      integer brain !< [Input] bra basis set handle
84      integer ish   !< [Input] shell/contraction index
85      integer jsh   !< [Input] shell/contraction index
86      integer ketin !< [Input] ket basis set handle
87      integer ksh   !< [Input] shell/contraction index
88      integer lsh   !< [Input] shell/contraction index
89      integer lscr  !< [Input] length of scratch array
90      double precision scr(lscr) !< [Scratch] array
91      integer leri  !< [In|Output] length of integral array
92      double precision eri(*) !< [Output] 2e4c integrals
93c:: local
94      integer bra, ket
95      integer ab_geom, cd_geom, ucont, ityp
96      integer La, a_prim, a_gen, a_iexp, a_icfp, a_cent, a_icfps
97      integer Lb, b_prim, b_gen, b_iexp, b_icfp, b_cent, b_icfps
98      integer Lc, c_prim, c_gen, c_iexp, c_icfp, c_cent, c_icfps
99      integer Ld, d_prim, d_gen, d_iexp, d_icfp, d_cent, d_icfps
100c
101      double precision roff(3)
102      double precision q4
103      integer nint
104      logical dum_log, do_bra, do_ket
105      logical status_sp, status_nw, status_txs, status_gen
106      integer texas_ang_limit
107      logical ieqj, keql
108      integer sbas, abas, bras, kets
109c
110      logical any_spherical
111      integer a_nbf,   b_nbf,   c_nbf,   d_nbf
112      integer a_nbf_s, b_nbf_s, c_nbf_s, d_nbf_s
113      integer ab_gen, ab_cmp, cd_gen, cd_cmp,i,j
114c
115#include "bas_exndcf_sfn.fh"
116#include "bas_ibs_sfn.fh"
117c
118c     timing(10)=timing(10)-util_wallsec()
119      bra = brain + BASIS_HANDLE_OFFSET
120      ket = ketin + BASIS_HANDLE_OFFSET
121      bras = bra
122      kets = ket
123      ab_geom = ibs_geom(bra)
124      cd_geom = ibs_geom(ket)
125      a_cent  = (sf_ibs_cn2ce(ish,bra))
126      b_cent  = (sf_ibs_cn2ce(jsh,bra))
127      c_cent  = (sf_ibs_cn2ce(ksh,ket))
128      d_cent  = (sf_ibs_cn2ce(lsh,ket))
129c
130c
131      any_spherical = bas_spherical(bra).or.bas_spherical(ket)
132c
133      status_sp = cando_sp(brain,ish,jsh).and.cando_sp(ketin,ksh,lsh)
134      if (.not.status_sp) then
135c
136        ieqj = ish.eq.jsh
137        keql = ksh.eq.lsh
138c
139        ucont   = sf_ibs_cn2ucn(ish,bra)
140        La      = infbs_cont(CONT_TYPE ,ucont,bra)
141        a_prim  = infbs_cont(CONT_NPRIM,ucont,bra)
142        a_gen   = infbs_cont(CONT_NGEN ,ucont,bra)
143        a_iexp  = infbs_cont(CONT_IEXP ,ucont,bra)
144        a_icfp  = infbs_cont(CONT_ICFP ,ucont,bra)
145        a_icfps = infbs_cont(CONT_ICFP ,ucont,bras)
146c
147        ucont   = sf_ibs_cn2ucn(jsh,bra)
148        Lb      = infbs_cont(CONT_TYPE ,ucont,bra)
149        b_prim  = infbs_cont(CONT_NPRIM,ucont,bra)
150        b_gen   = infbs_cont(CONT_NGEN ,ucont,bra)
151        b_iexp  = infbs_cont(CONT_IEXP ,ucont,bra)
152        b_icfp  = infbs_cont(CONT_ICFP ,ucont,bra)
153        b_icfps = infbs_cont(CONT_ICFP ,ucont,bras)
154c
155        ucont   = sf_ibs_cn2ucn(ksh,ket)
156        Lc      = infbs_cont(CONT_TYPE ,ucont,ket)
157        c_prim  = infbs_cont(CONT_NPRIM,ucont,ket)
158        c_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
159        c_iexp  = infbs_cont(CONT_IEXP ,ucont,ket)
160        c_icfp  = infbs_cont(CONT_ICFP ,ucont,ket)
161        c_icfps = infbs_cont(CONT_ICFP ,ucont,kets)
162c
163        ucont   = sf_ibs_cn2ucn(lsh,ket)
164        Ld      = infbs_cont(CONT_TYPE ,ucont,ket)
165        d_prim  = infbs_cont(CONT_NPRIM,ucont,ket)
166        d_gen   = infbs_cont(CONT_NGEN ,ucont,ket)
167        d_iexp  = infbs_cont(CONT_IEXP ,ucont,ket)
168        d_icfp  = infbs_cont(CONT_ICFP ,ucont,ket)
169        d_icfps = infbs_cont(CONT_ICFP ,ucont,kets)
170c
171        a_nbf = int_nbf_x(La)
172        b_nbf = int_nbf_x(Lb)
173        c_nbf = int_nbf_x(Lc)
174        d_nbf = int_nbf_x(Ld)
175        leri=a_nbf*b_nbf*c_nbf*d_nbf*a_gen*b_gen*c_gen*d_gen
176        call hnd_giahnd(
177     &      coords(1,a_cent,ab_geom), dbl_mb(mb_exndcf(a_iexp,bra)),
178     &      dbl_mb(mb_exndcf(a_icfp,bra)), a_prim, a_gen, La,
179     &      coords(1,b_cent,ab_geom), dbl_mb(mb_exndcf(b_iexp,bra)),
180     &      dbl_mb(mb_exndcf(b_icfp,bra)), b_prim, b_gen, Lb,
181     &      coords(1,c_cent,cd_geom), dbl_mb(mb_exndcf(c_iexp,ket)),
182     &      dbl_mb(mb_exndcf(c_icfp,ket)), c_prim, c_gen, Lc,
183     &      coords(1,d_cent,cd_geom), dbl_mb(mb_exndcf(d_iexp,ket)),
184     &      dbl_mb(mb_exndcf(d_icfp,ket)), d_prim,d_gen,Ld,
185     &      ieqj, keql, eri, leri, scr, lscr)
186c
187c eri has cartesian block of integrals (llo:lhi,klo:khi,jlo:jhi,ilo:ihi,6)
188c
189        if (any_spherical) then
190          a_nbf_s = int_nbf_s(La)
191          b_nbf_s = int_nbf_s(Lb)
192          c_nbf_s = int_nbf_s(Lc)
193          d_nbf_s = int_nbf_s(Ld)
194          cd_gen = c_gen*d_gen
195          ab_gen = a_gen*b_gen
196          ab_cmp = a_nbf*b_nbf
197          cd_cmp = c_nbf*d_nbf
198          do_bra=bas_spherical(bra)
199          do_ket=bas_spherical(ket)
200          call giao_to_sph(eri,leri,scr,La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf,
201     &                     d_nbf,a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s,
202     &                     a_gen,b_gen,c_gen,d_gen,ab_gen,ab_cmp,
203     &                     cd_gen,cd_cmp,do_bra,do_ket)
204        endif
205      else
206        write(luout,*)'int_giao_2e: cannot do sp integrals'
207        write(luout,*)' brain :',brain
208        write(luout,*)' ketin :',ketin
209        write(luout,*)' ish   :',ish
210        write(luout,*)' jsh   :',jsh
211        write(luout,*)' ksh   :',ksh
212        write(luout,*)' lsh   :',lsh
213        call errquit('int_giao_2e: fatal error ',0, INT_ERR)
214      endif
215      end
216c
217      subroutine giao_to_sph(eri,leri,scr,La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf,
218     &                     d_nbf,a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s,
219     &                     a_gen,b_gen,c_gen,d_gen,ab_gen,ab_cmp,
220     &                     cd_gen,cd_cmp,do_bra,do_ket)
221      implicit none
222c
223      integer leri
224      double precision eri(leri,6), scr(*)
225      integer La,Lb,Lc,Ld,a_nbf,b_nbf,c_nbf,d_nbf
226      integer a_nbf_s,b_nbf_s,c_nbf_s,d_nbf_s,a_gen,b_gen,c_gen,d_gen
227      integer ab_gen,ab_cmp,cd_gen,cd_cmp
228      logical do_bra,do_ket
229c
230      integer ityp
231c
232      ab_cmp = a_nbf_s*b_nbf_s
233      do ityp = 1, 6
234         call spcart_bra2etran(eri(1,ityp),scr,
235     &        b_nbf,a_nbf,b_nbf_s,a_nbf_s,
236     &        Lb, La, b_gen, a_gen,
237     &        cd_gen*cd_cmp,.false.)
238         call spcart_ket2etran(eri(1,ityp),scr,
239     &        d_nbf,c_nbf,d_nbf_s,c_nbf_s,
240     &        Ld, Lc, d_gen, c_gen,
241     &        ab_gen*ab_cmp,.false.)
242      enddo
243      a_nbf = a_nbf_s
244      b_nbf = b_nbf_s
245      c_nbf = c_nbf_s
246      d_nbf = d_nbf_s
247      end
248c
249C> @}
250