1C> \ingroup nwint
2C> @{
3C>
4C> \brief Compute 1-electron Gaussian periodic nuclear attraction
5C> integral derivatives
6C>
7C> See [1] for details.
8C>
9C> [1] JE Jaffe, AC Hess,
10C>     <i>"Gaussian basis density functional theory for systems
11C>     periodic in two or three dimensions: Energy and forces"</i>,
12C>    J.Chem.Phys. <b>105</b>, 10983-10998 (1996), DOI:
13C>    <a href="https://doi.org/10.1063/1.472866">
14C>    10.1063/1.472866</a>
15C>
16      subroutine intpd_1epe(i_basis,ish,Ri,j_basis,jsh,Rj,
17     &    lscr,scr,lpea,Pea)
18*
19* $Id$
20*
21      implicit none
22#include "stdio.fh"
23#include "errquit.fh"
24#include "nwc_const.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"
33c::external subroutines used
34c errquit
35c::functions
36      logical int_chk_init
37      logical int_chk_sh
38      logical cando_nw_1e
39      logical cando_nw
40      external int_chk_init
41      external int_chk_sh
42      external cando_nw_1e
43      external cando_nw
44      integer int_nint_cart
45      external int_nint_cart
46c::passed
47      integer i_basis           !< [Input] basis set handle for ish functions
48      integer j_basis           !< [Input] basis set handle for jsh functions
49      integer ish               !< [Input] lexical contraction/shell index
50      integer jsh               !< [Input] lexical contraction/shell index
51      integer lscr              !< [Input] length of the scratch array
52      integer lpea              !< [Input] length of potential energy integral array
53      double precision Pea(lpea) !< [Output] potential energy integral array
54      double precision scr(lscr) !< [Scratch] scratch array
55      double precision Ri(3)    !< [Input] translation vector for ish center (fractional coordinates)
56      double precision Rj(3)    !< [Input] translation vector for jsh center (fractional coordinates)
57c::local
58      logical shells_ok
59      integer i_geom, j_geom, ibas, jbas, ucont, mynint
60      integer Li, i_prim, i_gen, i_iexp, i_icfp, i_cent
61      integer Lj, j_prim, j_gen, j_iexp, j_icfp, j_cent
62      double precision xyz_new_i(3) ! new coordinates for ish function center
63      double precision xyz_new_j(3) ! new coordinates for jsh function center
64c
65      integer WarnP
66      save WarnP
67      data WarnP /0/
68c
69#include "bas_exndcf_sfn.fh"
70#include "bas_ibs_sfn.fh"
71c
72c check initialization and shells
73c
74      if (.not.int_chk_init('intpd_1epe'))
75     &    call errquit('intpd_1epe: int_init was not called' ,0,
76     &              INT_ERR)
77c
78      shells_ok = int_chk_sh(i_basis,ish)
79      shells_ok = shells_ok .and. int_chk_sh(j_basis,jsh)
80      if (.not.shells_ok)
81     &    call errquit('intpd_1epe: invalid contraction/shell',0,
82     &             BASIS_ERR)
83c
84      ibas = i_basis + BASIS_HANDLE_OFFSET
85      jbas = j_basis + BASIS_HANDLE_OFFSET
86c
87      ucont   = (sf_ibs_cn2ucn(ish,ibas))
88      Li      = infbs_cont(CONT_TYPE ,ucont,ibas)
89      i_prim  = infbs_cont(CONT_NPRIM,ucont,ibas)
90      i_gen   = infbs_cont(CONT_NGEN ,ucont,ibas)
91      i_iexp  = infbs_cont(CONT_IEXP ,ucont,ibas)
92      i_icfp  = infbs_cont(CONT_ICFP ,ucont,ibas)
93      i_cent  = (sf_ibs_cn2ce(ish,ibas))
94      i_geom  = ibs_geom(ibas)
95c
96      ucont   = (sf_ibs_cn2ucn(jsh,jbas))
97      Lj      = infbs_cont(CONT_TYPE ,ucont,jbas)
98      j_prim  = infbs_cont(CONT_NPRIM,ucont,jbas)
99      j_gen   = infbs_cont(CONT_NGEN ,ucont,jbas)
100      j_iexp  = infbs_cont(CONT_IEXP ,ucont,jbas)
101      j_icfp  = infbs_cont(CONT_ICFP ,ucont,jbas)
102      j_cent  = (sf_ibs_cn2ce(jsh,jbas))
103      j_geom  = ibs_geom(jbas)
104c
105      if (i_geom.ne.j_geom.and.WarnP.eq.0) then
106        write(luout,*)
107     &      'intpd_1epe: WARNING: possible geometry inconsistency'
108        write(luout,*)'i_basis geometry handle:',i_geom
109        write(luout,*)'j_basis geometry handle:',j_geom
110        WarnP = 1
111      endif
112c.. translate ish center coordinates based on Ri
113      call intp_txyz(i_cent,i_geom,Ri,xyz_new_i)
114c.. translate jsh center coordinates based on Rj
115      call intp_txyz(j_cent,j_geom,Rj,xyz_new_j)
116c
117      mynint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
118c
119      if (cando_nw(i_basis,ish,0).and.cando_nw(j_basis,jsh,0)) then
120        call hf1d(
121     &      xyz_new_i,dbl_mb(mb_exndcf(i_iexp,ibas)),
122     &      dbl_mb(mb_exndcf(i_icfp,ibas)), i_prim, i_gen, Li,
123     &      i_cent,
124     &      xyz_new_j,dbl_mb(mb_exndcf(j_iexp,jbas)),
125     &      dbl_mb(mb_exndcf(j_icfp,jbas)), j_prim, j_gen, Lj,
126     &      j_cent,
127     &      coords(1,1,i_geom),charge(1,i_geom),
128     &      geom_invnucexp(1,i_geom),ncenter(i_geom),
129c............................. doS     doT     doV    canonical
130     &      scr,scr,Pea,mynint,.false.,.false.,.true.,.false.,
131c.............. dryrun
132     &      .false.,scr,lscr)
133        if (bas_spherical(ibas).or.bas_spherical(jbas)) then
134          if (Li.eq.-1) i_gen = 1
135          if (Lj.eq.-1) j_gen = 1
136          call spcart_2cBtran(Pea,scr,lscr,
137     &        int_nbf_x(Li),int_nbf_s(Li),Li,i_gen,bas_spherical(ibas),
138     &        int_nbf_x(Lj),int_nbf_s(Lj),Lj,j_gen,bas_spherical(jbas),
139     &        (3*ncenter(i_geom)),.false.)
140        endif
141      else
142        call errquit('intpd_1epe: could not do sp or nw integrals',0,
143     &               INT_ERR)
144      endif
145c
146      end
147C> @}
148