1C> \ingroup nwint
2C> @{
3C>
4C> \brief Compute the nuclear attaction integrals for pseudo nucleii
5C>
6C> This subroutine is modified from int_1epe. It evaluates the nuclear
7C> attraction integrals for each pseudo nucleus at positions specified by the
8C> array qxyz, which are just the quardrature points. The effective charges
9C> of these pseudo nuclei are determined by the transition densities and/or
10C> the quardrature weights.
11C>
12c $Id$
13      subroutine int_vstat1e(i_basis,ish,j_basis,jsh,lscr,scr,nqpts,
14     &     vint,qxyz,zq,dens)
15c
16c     This subroutine is modified from int_1epe. It evaluates the nuclear
17c     attraction integrals for each pseudo nucleus at positions specified by the
18c     array qxyz, which are just the quardrature points. The effective charges
19c     of these pseudo nuclei are determined by the transition densities and/or
20c     the quardrature weights.
21c
22      implicit none
23#include "nwc_const.fh"
24#include "errquit.fh"
25#include "basP.fh"
26#include "basdeclsP.fh"
27#include "geomP.fh"
28#include "geobasmapP.fh"
29#include "mafdecls.fh"
30#include "bas_exndcf_dec.fh"
31#include "bas_ibs_dec.fh"
32#include "int_nbf.fh"
33#include "stdio.fh"
34#include "apiP.fh"
35c
36c     external subroutines used
37c
38      logical cando_hnd_1e
39      logical cando_nw_1e
40      logical cando_nw
41      logical int_chk_init
42      logical int_chk_sh
43      external int_chk_init
44      external int_chk_sh
45      external cando_hnd_1e
46      external cando_nw_1e
47      external cando_nw
48c
49      integer i_basis !< [Input] basis set handle for ish
50      integer ish     !< [Input] i shell/contraction
51      integer j_basis !< [Input] basis set handle for jsh
52      integer jsh     !< [Input] j shell/contraction
53      integer lscr    !< [Input] length of scratch array
54      double precision scr(lscr)    !< [Scratch] scratch array
55      integer nqpts                 !< [Input] length of potential buffer
56      double precision vint(nqpts)  !< [Output] potential integrals
57      double precision qxyz(3,nqpts) !< [Input] coordinates of q points
58      double precision zq(nqpts)    !< [Input] effective charges on q points
59      double precision dens(*)      !< [Input] elements of density matrix
60c::local
61      logical shells_ok, orel, oirel, ojrel, oNR
62      integer i_geom, j_geom, ibas, jbas, ucont
63      integer lbas, sbas, isbas, jsbas
64      integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent, i_icfpS
65      integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent, j_icfpS
66c
67      logical any_spherical
68      integer i_nbf_x, j_nbf_x
69      integer i_nbf_s, j_nbf_s
70      integer rel_dbg, rel_typ
71      integer ig,jg
72c
73      integer WarnP
74      save WarnP
75      data WarnP /0/
76c
77#include "bas_exndcf_sfn.fh"
78#include "bas_ibs_sfn.fh"
79c
80c check initialization and shells
81c
82      if (.not.int_chk_init('int_vstat1e'))
83     &       call errquit('int_vstat1e: int_init was not called' ,0,
84     &            INT_ERR)
85c
86      shells_ok = int_chk_sh(i_basis,ish)
87      shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh)
88      if (.not.shells_ok)
89     &       call errquit('int_vstat1e: invalid contraction/shell',0,
90     &               BASIS_ERR)
91c
92      ibas = i_basis + BASIS_HANDLE_OFFSET
93      jbas = j_basis + BASIS_HANDLE_OFFSET
94      isbas = ibas
95      jsbas = jbas
96      oNR = .true.
97c
98      ucont   = (sf_ibs_cn2ucn(ish,ibas))
99      Li      = infbs_cont(CONT_TYPE ,ucont,ibas)
100      i_prim  = infbs_cont(CONT_NPRIM,ucont,ibas)
101      i_gen   = infbs_cont(CONT_NGEN ,ucont,ibas)
102      i_iexp  = infbs_cont(CONT_IEXP ,ucont,ibas)
103      i_icfp  = infbs_cont(CONT_ICFP ,ucont,ibas)
104      i_cent  = (sf_ibs_cn2ce(ish,ibas))
105      i_geom  = ibs_geom(ibas)
106      i_icfpS  = infbs_cont(CONT_ICFP ,ucont,isbas)
107c
108      ucont   = (sf_ibs_cn2ucn(jsh,jbas))
109      Lj      = infbs_cont(CONT_TYPE ,ucont,jbas)
110      j_prim  = infbs_cont(CONT_NPRIM,ucont,jbas)
111      j_gen   = infbs_cont(CONT_NGEN ,ucont,jbas)
112      j_iexp  = infbs_cont(CONT_IEXP ,ucont,jbas)
113      j_icfp  = infbs_cont(CONT_ICFP ,ucont,jbas)
114      j_cent  = (sf_ibs_cn2ce(jsh,jbas))
115      j_geom  = ibs_geom(jbas)
116      j_icfpS  = infbs_cont(CONT_ICFP ,ucont,jsbas)
117c
118      if (i_geom.ne.j_geom.and.WarnP.eq.0) then
119        write(luout,*)
120     &      'int_vstat1e: WARNING: possible geometry inconsistency'
121        write(luout,*)'i_basis geometry handle:',i_geom
122        write(luout,*)'j_basis geometry handle:',j_geom
123        WarnP = 1
124      endif
125c
126      if (cando_hnd_1e(i_basis,ish,0).and.cando_hnd_1e(j_basis,jsh,0))
127     &     then
128         call hnd_vstat(
129     &        coords(1,i_cent,i_geom),dbl_mb(mb_exndcf(i_iexp,ibas)),
130     &        dbl_mb(mb_exndcf(i_icfp,ibas)),
131     &        i_prim, i_gen, Li,
132     &        coords(1,j_cent,j_geom),dbl_mb(mb_exndcf(j_iexp,jbas)),
133     &        dbl_mb(mb_exndcf(j_icfp,jbas)),
134     &        j_prim, j_gen, Lj,
135     &        qxyz, zq, nqpts,dens,
136c..............................doS     doT     doV
137     &        scr,scr,vint,nqpts,.false.,.false.,.true.,
138     &        scr,lscr)
139c
140      else
141         call errquit('int_vstat1e: could not do hnd, try sp or nw
142     &        integrals?', 0, INT_ERR)
143      endif
144c      write(*,*)"vints in int_vstat1e"
145c      write(*,"(4f12.8)")((qxyz(ig,jg),ig=1,3),vint(jg),jg=1,nqpts)
146c
147*     vint now has the cartesian integral block  (jlo:jhi,ilo:ihi)
148*
149      any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
150      if (.not.any_spherical) return
151      if (any_spherical)
152     &     call errquit("int_vstat1e: spherical basis not supported
153     &     yet", 100, BASIS_ERR)
154c
155c ... reset general contractions for sp shells to 1 since they are handled
156c     as a block of 4. Since int_nbf_* arrays are set to the appropriate size.
157c
158      if (li.eq.-1) i_gen = 1
159      if (lj.eq.-1) j_gen = 1
160c
161      if (bas_spherical(ibas).and.bas_spherical(jbas)) then
162*...  transform both i and j integrals
163         i_nbf_x = int_nbf_x(Li)
164         i_nbf_s = int_nbf_s(Li)
165         j_nbf_x = int_nbf_x(Lj)
166         j_nbf_s = int_nbf_s(Lj)
167c
168         call spcart_tran1e(vint,scr,
169     &        j_nbf_x,i_nbf_x,Lj,j_gen,
170     &        j_nbf_s,i_nbf_s,Li,i_gen,
171     &        .false.)
172      else if (bas_spherical(ibas)) then
173*..   transform on i component
174         i_nbf_x = int_nbf_x(Li)
175         i_nbf_s = int_nbf_s(Li)
176         j_nbf_x = int_nbf_x(Lj)
177         j_nbf_s = j_nbf_x
178         call spcart_tran1e(vint,scr,
179     &        j_nbf_x,i_nbf_x,0,j_gen,
180     &        j_nbf_s,i_nbf_s,Li,i_gen,
181     &        .false.)
182      else if (bas_spherical(jbas)) then
183*..   transform on j component
184         i_nbf_x = int_nbf_x(Li)
185         i_nbf_s = i_nbf_x
186         j_nbf_x = int_nbf_x(Lj)
187         j_nbf_s = int_nbf_s(Lj)
188         call spcart_tran1e(vint,scr,
189     &        j_nbf_x,i_nbf_x,Lj,j_gen,
190     &        j_nbf_s,i_nbf_s,0,i_gen,
191     &        .false.)
192      else
193         call errquit(
194     &        'int_vstat1e: should never reach transform blocked else',
195     &        911, INT_ERR)
196      endif
197c
198      end
199C> @}
200