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