1C> \ingroup nwint
2C> @{
3C>
4C> \brief Compute the 1-electron derivative integrals
5C>
6C> This routine compute the 1-electron derivative integrals. The routine
7C> computes the nuclear attraction integral derivatives
8C> \f{eqnarray*}{
9C>   \sum_A\frac{\partial (\mu|Z_AR_A^{-1}|\nu)}{\partial R_x} &=&
10C>      \sum_A\int \frac{\partial [g_\mu(X_\mu,r_1)Z_AR^{-1}_{A1}g_\nu(X_\nu,r_1)]}
11C>                      {\partial X_x}dr_1
12C> \f}
13C> It also adds the kinetic energy integral derivatives
14C> \f{eqnarray*}{
15C>   -\frac{1}{2}\frac{\partial (\mu|\nabla_1^2|\nu)}{\partial R_x} &=&
16C>   -\frac{1}{2}\int \frac{\partial [g_\mu(X_\mu,r_1)\nabla_1^2g_\nu(X_\nu,r_1)]}
17C>                         {\partial X_x}dr_1
18C> \f}
19C> The results are stored in `H1a` as if that variable was declared as
20C> `H1a(nint,ncoord,natom)` where `nint` is the number of integrals for a given
21C> pair of shells, `ncoord` is the number of Cartesian coordinates, and `natom`
22C> is the number of atoms.
23C>
24      subroutine intd_1eh1(i_basis,ish,j_basis,jsh,lscr,scr,
25     &       lH1a,H1a)
26C $Id$
27      implicit none
28#include "stdio.fh"
29#include "errquit.fh"
30#include "nwc_const.fh"
31#include "basP.fh"
32#include "basdeclsP.fh"
33#include "geomP.fh"
34#include "geobasmapP.fh"
35c
36c layer routine to compute the derivative 1 electron hamiltonian integrals
37c for shells/contractions ish,jsh
38c
39c Order is...   nint*3*nat (3=> xyz, nat=number of atoms)
40c
41c  /                   |
42c | nint,d <ij>        |
43c |      --------------|
44c  \     d[idatom(1),x]|
45c                          |
46c       nint,d <ij>        |
47c            --------------|
48c            d[idatom(1),y]|
49c                              |
50c           nint,d <ij>        |
51c                --------------|
52c                d[idatom(1),z]|
53c                                  |
54c               nint,d <ij>        |
55c                    --------------|
56c                    d[idatom(2),x]|
57c                                      |
58c                   nint,d <ij>        |
59c                        --------------|
60c                        d[idatom(2),y]|
61c                                           |
62c                       nint,d <ij>         |
63c                            -------------- |
64c                            d[idatom(2),z] |
65c
66c                                  . . .
67c                                                            |
68c                                         nint,d <ij>        |
69c                                              --------------|
70c                                            d[idatom(nat),x]|
71c                                                                |
72c                                             nint,d <ij>        |
73c                                                  --------------|
74c                                                d[idatom(nat),y]|
75c                                                                    \
76c                                                 nint,d <ij>         |
77c                                                      -------------- |
78c                                                    d[idatom(nat),z]/
79c
80c::functions
81      integer int_nint_cart
82      external int_nint_cart
83c::passed
84      integer i_basis   !< [Input] ish basis set handle
85      integer ish       !< [Input] ``i'' contraction index
86      integer j_basis   !< [Input] jsh basis set handle
87      integer jsh       !< [Input] ``j'' contraction index
88      integer lscr      !< [Input] length of scratch space
89      integer lH1a      !< [Input] number of h1 integral derivatives in shells ish and jsh
90c                       ! NOTE: nint*3 integral derivatives returned per unique center
91      double precision scr(lscr) !< [Input] scratch array
92      double precision H1a(*)    !< [Output] derivative integrals
93c
94c::local
95      integer nint, offset, scrsize, nat
96c
97      nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset)))
98c
99      nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
100      if (nint*3*nat.gt.lH1a) then
101        write(luout,*) 'nint*3*nat = ',nint*3*nat
102        write(luout,*) 'lH1a       = ',lH1a
103        call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR)
104      endif
105c
106      offset = nint*3*2       ! scratch for Ta array in intd_1eh1P routine
107      scrsize = lscr - offset ! new scratch array size
108      offset = offset + 1     ! increment for passing to intd_1eh1P
109c
110      call intd_1eh1P(i_basis,ish,j_basis,jsh,
111cc AJL/Begin/SPIN ECPs
112     &       scrsize,scr(offset),nint,H1a,scr,1)
113c                                             ^ ECP Integral (ALPHA)
114c
115      end
116      subroutine intd_1eh1P(i_basis,ish,j_basis,jsh,lscr,scr,
117     &       nint,H1a,Ta,ecp_channel)
118cc AJL/End
119      implicit none
120#include "stdio.fh"
121#include "errquit.fh"
122#include "apiP.fh"
123#include "nwc_const.fh"
124#include "int_nbf.fh"
125#include "basP.fh"
126#include "basdeclsP.fh"
127#include "geomP.fh"
128#include "geom.fh"
129#include "geobasmapP.fh"
130#include "mafdecls.fh"
131#include "bas_exndcf_dec.fh"
132#include "bas_ibs_dec.fh"
133#include "rel_nwc.fh"
134c::external subroutines used
135c... errquit
136c::functions
137      logical cando_hnd_1e
138      logical cando_nw
139      logical cando_sim
140      external cando_hnd_1e
141      external cando_sim
142c::passed
143      integer i_basis   ! [input] ish basis set handle
144      integer ish       ! [input] ``i'' contraction index
145      integer j_basis   ! [input] jsh basis set handle
146      integer jsh       ! [input] ``j'' contraction index
147      integer lscr      ! [input] length of scratch space
148      integer nint      ! [input] number of integrals in shells ish and jsh
149c                       ! NOTE: nint*3 integral derivatives returned per unique center
150      double precision scr(lscr) ! [input] scratch array
151      double precision H1a(nint,3,*)    ! [output] derivative integrals (nint,3,n_atoms)
152      double precision Ta(nint,3,2)     ! [scratch] space for kinetic integrals
153cc AJL/Begin/SPIN ECPs
154      integer ecp_channel !< [Input] ECP Channel (1 = Alpha, 2 = Beta)
155cc AJL/End
156c::local
157      logical doT
158      integer ucont,uconts
159      integer ibas,iatom,inp,igen,iexp,icf,itype,igeom,isbas,icfS
160      integer jbas,jatom,jnp,jgen,jexp,jcf,jtype,jgeom,jsbas,jcfS
161      integer nat
162      integer nintV
163      integer offset
164c
165      logical any_spherical
166      logical orel, oirel, ojrel, oNR
167      logical ohnd_ok, onw_ok,osim_ok
168      integer i_nbf_x, j_nbf_x
169      integer i_nbf_s, j_nbf_s
170      integer nint_x, nint_s
171      integer zatom, zyx
172      integer lbas, sbas, abas
173c
174#include "bas_exndcf_sfn.fh"
175#include "bas_ibs_sfn.fh"
176c
177c  check if gencon/sp shells
178c
179      call int_nogencont_check(i_basis,'intd_1eh1P:i_basis')
180      call int_nogencont_check(j_basis,'intd_1eh1P:j_basis')
181      call int_nospshell_check(i_basis,'intd_1eh1P:i_basis')
182      call int_nospshell_check(j_basis,'intd_1eh1P:j_basis')
183c
184      ibas = i_basis + BASIS_HANDLE_OFFSET
185      jbas = j_basis + BASIS_HANDLE_OFFSET
186      oNR = .true.
187      oirel = .false.
188      ojrel = .false.
189      orel = .false.
190c
191      if (dyall_mod_dir) then
192c
193c     get basis set handles; relativistic integral option only valid
194c     if both ibas and jbas are the ao basis.
195c
196        lbas = lc_bsh + BASIS_HANDLE_OFFSET
197        sbas = sc_bsh + BASIS_HANDLE_OFFSET
198        abas = ao_bsh + BASIS_HANDLE_OFFSET
199        orel = ibas .eq. abas .and. jbas .eq. abas
200      end if
201c
202c   i shell
203c
204      ucont = (sf_ibs_cn2ucn(ish,ibas))
205c
206c     check for relativistic shell
207c
208      if (orel .and. (infbs_cont(CONT_RELLS ,ucont,ibas) .ne. 0)) then
209        oirel = .true.
210        isbas = sbas
211        uconts = ao_to_ls(ucont)
212        if (uconts .eq. 0) call errquit (
213     &      'intd_1eh1: no relativistic pointer',911, INT_ERR)
214        if (nesc_1e_approx) then
215          ibas = lbas
216          ucont = uconts
217        end if
218      else
219        uconts = ucont
220        isbas = ibas
221      end if
222c
223      inp   = infbs_cont(CONT_NPRIM,ucont,ibas)
224      igen  = infbs_cont(CONT_NGEN,ucont,ibas)
225      iexp  = infbs_cont(CONT_IEXP,ucont,ibas)
226      icf   = infbs_cont(CONT_ICFP,ucont,ibas)
227      itype = infbs_cont(CONT_TYPE,ucont,ibas)
228      igeom = ibs_geom(ibas)
229      iatom = (sf_ibs_cn2ce(ish,ibas))
230      icfS  = infbs_cont(CONT_ICFP ,uconts,isbas)
231c
232c   j shell
233c
234      ucont = (sf_ibs_cn2ucn(jsh,jbas))
235c
236c     check for relativistic shell
237c
238      if (orel .and. (infbs_cont(CONT_RELLS ,ucont,jbas) .ne. 0)) then
239        ojrel = .true.
240        jsbas = sbas
241        uconts = ao_to_ls(ucont)
242        if (uconts .eq. 0) call errquit (
243     &      'intd_1eh1: no relativistic pointer',911, INT_ERR)
244        if (nesc_1e_approx) then
245          jbas = lbas
246          ucont = uconts
247        end if
248      else
249        uconts = ucont
250        jsbas = jbas
251      end if
252c
253      jnp   = infbs_cont(CONT_NPRIM,ucont,jbas)
254      jgen  = infbs_cont(CONT_NGEN,ucont,jbas)
255      jexp  = infbs_cont(CONT_IEXP,ucont,jbas)
256      jcf   = infbs_cont(CONT_ICFP,ucont,jbas)
257      jtype = infbs_cont(CONT_TYPE,ucont,jbas)
258      jgeom = ibs_geom(jbas)
259      jatom = (sf_ibs_cn2ce(jsh,jbas))
260      jcfS  = infbs_cont(CONT_ICFP ,uconts,jsbas)
261c
262      oNR = .not.(oirel.and.ojrel)
263      orel = oirel.or.ojrel
264c
265      if (igeom.ne.jgeom) then
266        write(luout,*)'intd_1eh1P.F: two different geometries for',
267     &         ' derivatives?'
268        call errquit('intd_1eh1P: geom error ',911, INT_ERR)
269      endif
270c
271      if (iatom.eq.jatom) then
272        doT = .false.
273      else
274        doT = .true.
275      endif
276c
277      ohnd_ok = cando_hnd_1e(i_basis,ish,0)
278     &    .and. cando_hnd_1e(j_basis,jsh,0)
279     &    .and. (.not.geom_any_finuc (igeom))
280     &    .and. (.not.geom_any_finuc (jgeom))
281      onw_ok = cando_nw(i_basis,ish,0) .and. cando_nw(j_basis,jsh,0)
282      osim_ok = cando_sim(i_basis,ish,0) .and. cando_sim(j_basis,jsh,0)
283c
284      if (orel) then
285        call rel_oneld (
286     &      coords(1,iatom,igeom),
287     &      dbl_mb(mb_exndcf(iexp,ibas)),
288     &      dbl_mb(mb_exndcf(icf,ibas)),
289     &      dbl_mb(mb_exndcf(icfS,isbas)),inp,igen,itype,iatom,
290     &      coords(1,jatom,jgeom),
291     &      dbl_mb(mb_exndcf(jexp,jbas)),
292     &      dbl_mb(mb_exndcf(jcf,jbas)),
293     &      dbl_mb(mb_exndcf(jcfS,jsbas)),jnp,jgen,jtype,jatom,
294     &      coords(1,1,igeom),charge(1,igeom),
295     &      geom_invnucexp(1,igeom),ncenter(igeom),
296c........................     doS   doT  doV    canAB
297     &      scr,Ta,H1a,nint,.false.,doT,.true.,.false.,onw_ok,
298c...........       nonrel dryrun
299     &      ohnd_ok,oNR,.false.,scr,lscr,rel_dbg,rel_typ)
300c     send it here for SImint bit
301      else if (onw_ok.or.osim_ok) then
302        call hf1d(
303     &      coords(1,iatom,igeom),
304     &      dbl_mb(mb_exndcf(iexp,ibas)),
305     &      dbl_mb(mb_exndcf(icf,ibas)),
306     &      inp,igen,itype,iatom,
307c
308     &      coords(1,jatom,jgeom),
309     &      dbl_mb(mb_exndcf(jexp,jbas)),
310     &      dbl_mb(mb_exndcf(jcf,jbas)),
311     &      jnp,jgen,jtype,jatom,
312c
313     &      coords(1,1,igeom),charge(1,igeom),
314     &      geom_invnucexp(1,igeom),ncenter(igeom),
315     &      scr,Ta,H1a,nint,
316c..............overlap, k-e,  pot-e,  canab,   dryrun
317     &      .false., doT, .true., .false., .false.,
318     &      scr,lscr)
319      elseif (ohnd_ok) then
320        call hnd_stvintd(
321     &      coords(1,iatom,igeom),
322     &      dbl_mb(mb_exndcf(iexp,ibas)),
323     &      dbl_mb(mb_exndcf(icf,ibas)),
324     &      inp,igen,itype,iatom,
325c
326     &      coords(1,jatom,jgeom),
327     &      dbl_mb(mb_exndcf(jexp,jbas)),
328     &      dbl_mb(mb_exndcf(jcf,jbas)),
329     &      jnp,jgen,jtype,jatom,
330c
331     &      coords(1,1,igeom),charge(1,igeom),ncenter(igeom),
332     &      scr,Ta,H1a,nint,
333c............overlap, k-e,     pot-e,
334     &      .false.,  doT, .true.,
335     &      scr,lscr)
336      else
337        call errquit('intd_1eh1: could not do hnd or nw integrals',
338     &                0, INT_ERR)
339      endif
340c
341c if needed add in Ta derivative integrals
342c
343      if (doT) then
344        call daxpy(nint*3,1.0d00,Ta(1,1,1),1,H1a(1,1,iatom),1)
345        call daxpy(nint*3,1.0d00,Ta(1,1,2),1,H1a(1,1,jatom),1)
346      endif
347c
348c check for ecp
349c
350*
351* this should move to hf1dsp when sp is enabled.
352*
353      nat = ncenter(igeom)  ! needed for both ecp and spherical
354      if (any_ecp) then
355        nintV = int_nbf_x(itype)*int_nbf_x(jtype)
356        offset = nintV*3*nat + 1
357*       write(luout,*)' lscr to ecp_hf1:',(lscr-nintV)
358*       call dcopy(nintV*3*nat,0.0d00,0,scr,1) ! buffer zeroed in ecp_gradient
359c
360cc AJL/Begin/SPIN ECPs
361        if (ecp_channel.eq.1) then
362          call intd_ecp_hf1(
363     &        coords(1,iatom,igeom),
364     &        dbl_mb(mb_exndcf(iexp,ibas)),
365     &        dbl_mb(mb_exndcf(icf,ibas)),
366     &        inp,igen,itype,iatom,
367c
368     &        coords(1,jatom,jgeom),
369     &        dbl_mb(mb_exndcf(jexp,jbas)),
370     &        dbl_mb(mb_exndcf(jcf,jbas)),
371     &        jnp,jgen,jtype,jatom,
372c
373     &        scr,nintV,nat,
374     &        scr(offset),(lscr-offset-1),
375     &        .false.)
376c
377        else
378          call intd_ecp_hf1_beta(
379     &        coords(1,iatom,igeom),
380     &        dbl_mb(mb_exndcf(iexp,ibas)),
381     &        dbl_mb(mb_exndcf(icf,ibas)),
382     &        inp,igen,itype,iatom,
383c
384     &        coords(1,jatom,jgeom),
385     &        dbl_mb(mb_exndcf(jexp,jbas)),
386     &        dbl_mb(mb_exndcf(jcf,jbas)),
387     &        jnp,jgen,jtype,jatom,
388c
389     &        scr,nintV,nat,
390     &        scr(offset),(lscr-offset-1),
391     &        .false.)
392c
393        endif
394c*... sum ecp into derivative H1 block
395        call daxpy(nintV*3*nat,1.0d00,scr,1,H1a,1)
396cc AJL/End
397c
398      endif
399c
400*     H1a now has the cartesian integral block  (jlo:jhi,ilo:ihi)
401*
402      any_spherical = bas_spherical(ibas).or.bas_spherical(jbas)
403      if (.not.any_spherical) return
404c
405c ... reset general contractions for sp shells to 1 since they are handled
406c     as a block of 4.
407c
408      if (itype.eq.-1) igen = 1
409      if (jtype.eq.-1) jgen = 1
410c
411      if (bas_spherical(ibas).and.bas_spherical(jbas)) then
412*... transform both i and j integrals
413        i_nbf_x = int_nbf_x(Itype)
414        i_nbf_s = int_nbf_s(Itype)
415        j_nbf_x = int_nbf_x(Jtype)
416        j_nbf_s = int_nbf_s(Jtype)
417c
418        do zatom = 1,nat
419          do zyx = 1,3
420            call spcart_tran1e(H1a(1,zyx,zatom),scr,
421     &          j_nbf_x,i_nbf_x,Jtype,jgen,
422     &          j_nbf_s,i_nbf_s,Itype,igen,
423     &          .false.)
424          enddo
425        enddo
426      else if (bas_spherical(ibas)) then
427*.. transform on i component
428        i_nbf_x = int_nbf_x(Itype)
429        i_nbf_s = int_nbf_s(Itype)
430        j_nbf_x = int_nbf_x(Jtype)
431        j_nbf_s = j_nbf_x
432        do zatom = 1,nat
433          do zyx = 1,3
434            call spcart_tran1e(H1a(1,zyx,zatom),scr,
435     &          j_nbf_x,i_nbf_x,0,jgen,
436     &          j_nbf_s,i_nbf_s,Itype,igen,
437     &          .false.)
438          enddo
439        enddo
440      else if (bas_spherical(jbas)) then
441*.. transform on j component
442        i_nbf_x = int_nbf_x(Itype)
443        i_nbf_s = i_nbf_x
444        j_nbf_x = int_nbf_x(Jtype)
445        j_nbf_s = int_nbf_s(Jtype)
446        do zatom = 1,nat
447          do zyx = 1,3
448            call spcart_tran1e(H1a(1,zyx,zatom),scr,
449     &        j_nbf_x,i_nbf_x,Jtype,jgen,
450     &        j_nbf_s,i_nbf_s,0,igen,
451     &        .false.)
452          enddo
453        enddo
454      else
455        call errquit(
456     &        'int_1eh1: should never reach transform blocked else',
457     &        911, INT_ERR)
458      endif
459c
460c now shuffle transformed buffers to contiguous space
461c
462      nint_x = i_nbf_x*j_nbf_x
463      nint_s = i_nbf_s*j_nbf_s
464      if (nint_s.gt.nint_x) then
465        call errquit
466     &      ('intd_1eh1: nint_s >.nint_x diff=',(nint_s-nint_x),
467     &      INT_ERR)
468      elseif (nint_s.eq.nint_x) then
469        return
470      else
471        call int_c2s_mv
472     &      (H1a,nint_x,nint_s,(3*nat),scr,lscr,'intd_1eh1')
473      endif
474c
475
476      end
477cc AJL/Begin/SPIN ECPs
478cc This is a wrapper, identical to intd_1epot except we have
479cc separated out for alpha and beta channels
480      subroutine intd_1eh1_beta(i_basis,ish,j_basis,jsh,lscr,scr,
481     &       lH1a,H1a)
482      implicit none
483#include "stdio.fh"
484#include "errquit.fh"
485#include "nwc_const.fh"
486#include "basP.fh"
487#include "basdeclsP.fh"
488#include "geomP.fh"
489#include "geobasmapP.fh"
490c::functions
491      integer int_nint_cart
492      external int_nint_cart
493c::passed
494      integer i_basis   !< [Input] ish basis set handle
495      integer ish       !< [Input] ``i'' contraction index
496      integer j_basis   !< [Input] jsh basis set handle
497      integer jsh       !< [Input] ``j'' contraction index
498      integer lscr      !< [Input] length of scratch space
499      integer lH1a      !< [Input] number of h1 integral derivatives in shells ish and jsh
500c                       ! NOTE: nint*3 integral derivatives returned per unique center
501      double precision scr(lscr) !< [Input] scratch array
502      double precision H1a(*)    !< [Output] derivative integrals
503c
504c::local
505      integer nint, offset, scrsize, nat
506c
507      nat = ncenter(ibs_geom((i_basis + Basis_Handle_Offset)))
508c
509      nint = int_nint_cart(i_basis,ish,j_basis,jsh,0,0,0,0)
510      if (nint*3*nat.gt.lH1a) then
511        write(luout,*) 'nint*3*nat = ',nint*3*nat
512        write(luout,*) 'lH1a       = ',lH1a
513        call errquit('intd_1eh1: nint>lH1a error',911, INT_ERR)
514      endif
515c
516      offset = nint*3*2       ! scratch for Ta array in intd_1eh1P routine
517      scrsize = lscr - offset ! new scratch array size
518      offset = offset + 1     ! increment for passing to intd_1eh1P
519c
520      call intd_1eh1P(i_basis,ish,j_basis,jsh,
521     &       scrsize,scr(offset),nint,H1a,scr,2)
522c                                             ^ ECP Integral (BETA)
523      end
524cc AJL/End
525C> @}
526