1c $Id$
2*
3C> \ingroup nwint
4C> @{
5C>
6C> \brief Compute the 1-electron GIAO overlap integrals
7C> perturbed by the magnetic field
8C>
9C> Compute the the 1-electron GIAO overlap integrals
10C> perturbed by the magnetic field.
11C> For a definition of the GIAO basis functions see the
12C> `int_giao_2e` routine (see also [1]).
13C>
14C> [1] M. Dupuis,
15C>     "New integral transforms for molecular properties and application
16C>     to a massively parallel GIAO-SCF implementation",
17C>     Comp. Phys. Comm. <b>134</b>, 150-166 (2001), DOI:
18C>     <a href="https://doi.org/10.1016/S0010-4655(00)00195-8">
19C>     10.1016/S0010-4655(00)00195-8</a>.
20C>
21c:tex-% this is part of the API Standard Integral routines.
22c:tex-\subsection{int\_giaos10}
23c:tex-This routine computes the 1-elec GIAO integrals of s perturbed by the magnetic field
24c:tex-
25c:tex-{\it Syntax:}
26c:tex-\begin{verbatim}
27      subroutine int_giaos10(i_basis,ish,j_basis,jsh,lscr,scr,ls10,s10)
28c:tex-\end{verbatim}
29      implicit none
30#include "nwc_const.fh"
31#include "errquit.fh"
32#include "basP.fh"
33#include "basdeclsP.fh"
34#include "geomP.fh"
35#include "geobasmapP.fh"
36#include "mafdecls.fh"
37#include "bas_exndcf_dec.fh"
38#include "bas_ibs_dec.fh"
39#include "int_nbf.fh"
40#include "stdio.fh"
41#include "apiP.fh"
42#include "util.fh"
43c::external subroutines used
44c... errquit
45c::functions
46      logical cando_hnd_1e_prp
47      logical int_chk_init
48      logical int_chk_sh
49      external int_chk_init
50      external int_chk_sh
51      external cando_hnd_1e_prp
52c::passed
53c:tex-\begin{verbatim}
54      integer i_basis !< [Input] basis set handle for ish
55      integer ish     !< [Input] i shell/contraction
56      integer j_basis !< [Input] basis set handle for jsh
57      integer jsh     !< [Input] j shell/contraction
58      integer lscr    !< [Input] length of scratch array
59      double precision scr(lscr) !< [Scratch] scratch array
60      integer ls10               !< [Input] length of s10 buffer
61      double precision s10(ls10) !< [Output] s10 integrals
62c:tex-\end{verbatim}
63c::local
64      integer igeom, jgeom, ibas, jbas, ucont
65      integer itype, inp, igen, iexp, icent, icf, iatom
66      integer jtype, jnp, jgen, jexp, jcent, jcf, jatom
67c
68      logical any_spherical, trani, tranj, shells_ok
69      integer i_nbf_x, j_nbf_x
70      integer i_nbf_s, j_nbf_s
71      integer ipts, ncartint ,i,j
72c
73#include "bas_exndcf_sfn.fh"
74#include "bas_ibs_sfn.fh"
75c
76c check initialization and shells
77c
78      if (.not.int_chk_init('int_giaos10'))
79     &       call errquit('int_giaos10: int_init was not called' ,0,
80     &       INT_ERR)
81c
82      shells_ok = int_chk_sh(i_basis,ish)
83      shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh)
84      if (.not.shells_ok)
85     &       call errquit('int_giaos10: invalid contraction/shell',0,
86     &       INT_ERR)
87c
88c  check if gencont
89c
90      call int_nogencont_check(i_basis,'int_giaos10:i_basis')
91      call int_nogencont_check(j_basis,'int_giaos10:j_basis')
92c
93      ibas = i_basis + basis_handle_offset
94      jbas = j_basis + basis_handle_offset
95c
96      ucont = (sf_ibs_cn2ucn(ish,ibas))
97      itype = infbs_cont(CONT_TYPE ,ucont,ibas)
98      inp   = infbs_cont(CONT_NPRIM,ucont,ibas)
99      igen  = infbs_cont(CONT_NGEN ,ucont,ibas)
100      iexp  = infbs_cont(CONT_IEXP ,ucont,ibas)
101      icf   = infbs_cont(CONT_ICFP ,ucont,ibas)
102      iatom = (sf_ibs_cn2ce(ish,ibas))
103      igeom = ibs_geom(ibas)
104c
105      ucont = (sf_ibs_cn2ucn(jsh,jbas))
106      jtype = infbs_cont(CONT_TYPE ,ucont,jbas)
107      jnp   = infbs_cont(CONT_NPRIM,ucont,jbas)
108      jgen  = infbs_cont(CONT_NGEN ,ucont,jbas)
109      jexp  = infbs_cont(CONT_IEXP ,ucont,jbas)
110      jcf   = infbs_cont(CONT_ICFP ,ucont,jbas)
111      jatom = (sf_ibs_cn2ce(jsh,jbas))
112      jgeom = ibs_geom(jbas)
113c
114      if (igeom.ne.jgeom) then
115        write(luout,*)'int_giaos10: two different geometries for',
116     &         ' properties?'
117        call errquit('int_giaos10: geom error ',911, GEOM_ERR)
118      endif
119c
120c     Determine # of cartesian integrals in block
121c
122      ncartint = int_nbf_x(itype)*int_nbf_x(jtype)
123c
124      call hnd_giaos10(
125     &       coords(1,iatom,igeom),
126     &       dbl_mb(mb_exndcf(iexp,ibas)),
127     &       dbl_mb(mb_exndcf(icf,ibas)),
128     &       inp,igen,itype,
129c
130     &       coords(1,jatom,jgeom),
131     &       dbl_mb(mb_exndcf(jexp,jbas)),
132     &       dbl_mb(mb_exndcf(jcf,jbas)),
133     &       jnp,jgen,jtype,
134c
135     &       ncartint,s10,scr,lscr)
136c
137c     h01 now has three blocks, one for each component x, y, z
138c
139      any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
140      if (.not.any_spherical) return
141c
142      i_nbf_x = int_nbf_x(itype)
143      j_nbf_x = int_nbf_x(jtype)
144c
145c... assume we need to transform both i and j integrals
146c
147      trani = .true.
148      tranj = .true.
149*.. do not tranform i component
150      if (.not.bas_spherical(ibas)) trani = .false.
151*.. do not tranform j component
152      if (.not.bas_spherical(jbas)) tranj = .false.
153c
154c ... reset general contractions for sp shells to 1 since they are handled
155c     as a block of 4.
156c
157      if (itype.eq.-1) igen = 1
158      if (jtype.eq.-1) jgen = 1
159      call spcart_2cBtran(s10,scr,lscr,
160     &    j_nbf_x,int_nbf_s(jtype),jtype,jgen,tranj,
161     &    i_nbf_x,int_nbf_s(itype),itype,igen,trani,
162     &    3,.false.)
163c
164c     We now have the integrals in array (nsph_ints,3)
165c
166      return
167      end
168C> @}
169