1      subroutine int_giao_1ega(ibas,jbas,g,integ_type,xyzpt,nat,
2     &                         oskel)
3C$Id$
4      implicit none
5#include "errquit.fh"
6#include "mafdecls.fh"
7#include "global.fh"
8#include "rtdb.fh"
9#include "inp.fh"
10#include "apiP.fh"
11#include "bas.fh"
12#include "sym.fh"
13#include "geom.fh"
14c
15c     Compute the desired type of 1e GIAO integrals
16c     and ADD them into the given global array.
17c     This version computes the full square of integrals and should work
18c     OK even if ibas != jbas.
19c
20c     Oskel indicates that the skeleton (petite-list symmetry) matrix should be
21c     built ... requires that ibas = jbas.
22c
23c     arguments
24c
25      integer ibas, jbas            ! [input] bra and ket basis sets
26      integer g                     ! [output] GA handle to array, one for each field direction (if needed)
27      character*(*) integ_type      ! [input] Name of integrals to compute
28      logical oskel                 ! [input] If true generate symmetry unique list
29      integer nat                   ! [input] number of atoms for which we get the integrals (if needed)
30      double precision xyzpt(3,*)   ! [input] coordinates of requested atoms (if needed)
31c
32c     local variables
33c
34      integer mone, two
35      parameter (mone=-1, two=2)
36      integer nshell_i, nshell_j
37      integer ishell, jshell, iproc, nproc, mem1, max1e
38      integer ijshell, ilo, ihi, jlo, jhi, ilen, jlen
39      integer l_buf, l_scr
40      integer k_buf, k_scr
41      integer alo(3), ahi(3), ld(2)
42      integer type
43      logical odoit
44      double precision q2
45      integer nblocks
46c
47      logical odbug
48      logical osome
49c
50      integer irtdb
51      integer int_get_rtdb
52      external int_get_rtdb
53      integer iextbq
54      integer nbq,nextbq,ncosbq,ntmp
55      integer l_qbq,l_cbq
56      integer k_qbq,k_cbq
57      integer k_qextbq,k_cextbq
58c
59      odbug=.false.
60      osome=.false.
61      osome=osome.or.odbug
62      odbug=odbug.and.(ga_nodeid().eq.0)
63      osome=osome.and.(ga_nodeid().eq.0)
64      if(osome) then
65         write(6,*) 'in -int_giao_1ega- ... integ_type = ',
66     $              integ_type,ga_nodeid()
67         call util_flush(6)
68      endif
69c
70      call ga_sync()
71c
72      if (oskel) then
73         if (ibas.ne.jbas) call errquit
74     $      ('int_giao_1ega: use of symmetry requires ibas=jbas', ibas,
75     &       BASIS_ERR)
76      end if
77c
78      if (inp_compare(.false., integ_type, 's10')) then
79         type = 1
80         nblocks = 3
81      elseif (inp_compare(.false., integ_type, 'srxRb')) then
82         type = 2
83         nblocks = 3
84      else if (inp_compare(.false., integ_type, 'l10')) then
85         type = 3
86         nblocks = 3
87      else if (inp_compare(.false., integ_type, 'tv10')) then
88         type = 4
89         nblocks = 3
90      else if (inp_compare(.false., integ_type, 'h01')) then
91         type = 5
92         nblocks = 3 * nat
93      else if (inp_compare(.false., integ_type, 'h11 para'))then
94         type = 6
95         nblocks = 9 * nat
96      else if (inp_compare(.false., integ_type, 'h11 dia'))then
97         type = 7
98         nblocks = 9 * nat
99      else if (inp_compare(.false., integ_type, 'h11 all'))then
100         type = 8
101         nblocks = 9 * nat
102      else if (inp_compare(.false., integ_type, 'dso'))then
103         type = 9
104         nblocks = 9 * nat ! nat is number of pairs in this case : nat * (nat-1)
105      else if (inp_compare(.false., integ_type, 'pso'))then
106         type = 10
107         nblocks = 3 * nat
108      else if (inp_compare(.false., integ_type, 'fc'))then
109         type = 11
110         nblocks = nat
111      else if (inp_compare(.false., integ_type, 'sd+fc'))then
112         type = 12
113         nblocks = 6 * nat
114      else if (inp_compare(.false., integ_type, 'velocity'))then
115         type = 13
116         nblocks = 3
117      else if (inp_compare(.false., integ_type, 'angmom'))then
118         type = 14
119         nblocks = 3
120      else if (inp_compare(.false., integ_type, 'bq10')) then
121         type = 15
122         nblocks = 3
123      else
124         write(6,*) ' integ_type = ', integ_type,ga_nodeid()
125         call errquit('int_giao_1ega: unknown integ_type', 0, INT_ERR)
126      end if
127c
128c     Get info about the basis sets
129c
130      if (.not. bas_numcont(ibas, nshell_i)) call errquit
131     $     ('rhf_fock_1e: bas_numcont failed for ibas', ibas,
132     &       BASIS_ERR)
133      if (.not. bas_numcont(jbas, nshell_j)) call errquit
134     $     ('rhf_fock_1e: bas_numcont failed for jbas', jbas,
135     &       BASIS_ERR)
136c
137c     allocate necessary local temporary arrays on the stack
138c
139c     l_buf ... buffer to hold shell block of matrix
140c     l_s   ... buffer to hold shell block of matrix
141c     l_scr ... workspace for integral routines
142c
143c     k_* are the offsets corrsponding to the l_* handles
144c
145c
146      if (type.lt.9) then
147         call int_mem_1e(max1e, mem1)
148         max1e = max1e*nblocks
149      elseif (type.eq.9) then
150         call int_init_dso(max1e,mem1,ibas,nat)
151      elseif (type.eq.10) then
152         call int_init_pso(max1e,mem1,ibas,nat)
153      elseif (type.eq.11.or.type.eq.12) then
154         call int_init_1eelec(max1e,mem1,ibas,2,nat)
155      elseif (type.eq.15) then
156         call int_mem_1e(max1e, mem1)
157         max1e = max1e*nblocks
158      else
159         call int_init_dip(max1e,mem1,ibas)
160         max1e = max1e*nblocks
161      endif
162      mem1 = max(mem1,max1e)
163c
164      if(.not.MA_push_get(MT_DBL,max1e,'int_giao_1ega:buf',l_buf,k_buf))
165     $     call errquit('int_giao_1ega: ma failed', max1e, MA_ERR)
166      if(.not.MA_push_get(MT_DBL, mem1,'int_giao_1ega:scr',l_scr,k_scr))
167     $     call errquit('int_giao_1ega: ma failed', mem1, MA_ERR)
168c
169c     Get bq charges (external bq and cosmo)
170c
171      if (type.eq.15) then
172        nbq=0
173        nextbq = 0
174        ncosbq = 0
175        if(geom_extbq_on()) then
176           nextbq = geom_extbq_ncenter()
177           k_cextbq = geom_extbq_coord()
178           k_qextbq = geom_extbq_charge()
179           nbq = nextbq  ! external bq centers
180        end if
181c
182c       Get rtdb handle
183c
184        irtdb = int_get_rtdb()   ! get rtdb handle
185        if (rtdb_get(irtdb,'cosmo:nefc',mt_int,1,ncosbq))
186     &       nbq = ncosbq  ! cosmo bq centers
187c
188        if (nextbq.gt.0.and.ncosbq.gt.0)
189     &       nbq = nextbq + ncosbq  ! all bq centers
190c
191c       Allocate memory for all the bq charges
192c
193        if(.not.MA_push_get(MT_DBL,nbq,'qbq',l_qbq,k_qbq))
194     &    call errquit('int_giao_1ega: ma failed', l_qbq, MA_ERR)
195        call dfill(nbq,0.d0,dbl_mb(k_qbq),1)
196        if(.not.MA_push_get(MT_DBL,3*nbq,'cbq',l_cbq,k_cbq))
197     $    call errquit('int_giao_1ega: ma failed', l_cbq, MA_ERR)
198        call dfill(3*nbq,0.d0,dbl_mb(k_cbq),1)
199c
200c       Assign external bq charges and coordinates
201c
202        if (nextbq.gt.0) then
203         call dcopy(nextbq,dbl_mb(k_qextbq),1,dbl_mb(k_qbq),1)
204         call dcopy(3*nextbq,dbl_mb(k_cextbq),1,dbl_mb(k_cbq),1)
205        end if ! nextbq
206c
207c       Get cosmo charges and coordinates
208c
209        if (ncosbq.gt.0) then
210         call bq_fromrtdb(irtdb,'cosmo:efcz','cosmo:efcc',
211     &       ntmp,dbl_mb(k_qbq+nextbq),dbl_mb(k_cbq+3*nextbq))
212        end if ! ncosbq gt 0
213      end if ! type .eq. 15
214c
215c     Loop thru shells with static parallel work decomposition
216c
217      iproc = ga_nodeid()
218      nproc = ga_nnodes()
219      ijshell = 0
220      q2 = 1.0d0
221      do jshell = 1, nshell_j
222         do ishell = 1, nshell_i
223c
224            if (mod(ijshell, nproc) .eq. iproc) then
225               odoit = .true.
226               if (oskel)
227     $              odoit = sym_shell_pair(ibas, ishell, jshell, q2)
228c
229               if (odoit) then
230                  if (.not. bas_cn2bfr(ibas, ishell, ilo, ihi))
231     $                 call errquit('int_1e_ga: bas_cn2bfr ?', ibas,
232     &       BASIS_ERR)
233                  if (.not. bas_cn2bfr(jbas, jshell, jlo, jhi))
234     $                 call errquit('int_1e_ga: bas_cn2bfr ?', jbas,
235     &       BASIS_ERR)
236c
237                  ilen = ihi-ilo+1
238                  jlen = jhi-jlo+1
239c
240c     Generate the integrals
241c
242                  if (type .eq. 1) then      ! 3
243                     call int_giaos10 (jbas, jshell, ibas, ishell,
244     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
245                  else if (type .eq. 2) then ! 3
246                     call int_giaos100(jbas, jshell, ibas, ishell,
247     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
248                  else if (type .eq. 3) then ! 3
249                     call int_giaol10 (jbas, jshell, ibas, ishell,
250     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
251                  else if (type .eq. 4) then ! 3
252                     call int_giaotv10 (jbas, jshell, ibas, ishell,
253     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
254                  else if (type .eq. 5) then ! 3*nat
255                     call int_giaoh01 (jbas, jshell, ibas, ishell,
256     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
257     $                    xyzpt, nat)
258                  else if (type .eq. 6) then ! 9*nat
259                     call int_giaoh11 (jbas, jshell, ibas, ishell,
260     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
261     $                    xyzpt, nat, .true.,.false.)
262                  else if (type .eq. 7) then ! 9*nat
263                     call int_giaoh11 (jbas, jshell, ibas, ishell,
264     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
265     $                    xyzpt, nat, .false.,.true.)
266                  else if (type .eq. 8) then ! 9*nat
267                     call int_giaoh11 (jbas, jshell, ibas, ishell,
268     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
269     $                    xyzpt, nat, .true.,.true.)
270                  else if (type .eq. 9) then ! 9*nat*(nat-1)
271                     call int_dso (jbas, jshell, ibas, ishell,
272     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
273     $                    xyzpt, nat)
274                  else if (type .eq. 10) then ! 3*nat
275                     call int_pso (jbas, jshell, ibas, ishell,
276     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
277     $                    xyzpt, nat)
278                  else if (type .eq. 11) then ! nat
279                     call int_1eelec(jbas, jshell, ibas,ishell,
280     &                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
281     &                    mone,xyzpt,nat)
282                  else if (type .eq. 12) then ! 6*nat
283                     call int_1eelec(jbas, jshell, ibas,ishell,
284     &                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
285     &                    two,xyzpt,nat)
286                  else if (type .eq. 13) then ! 3
287                     call int_veloc(jbas, jshell, ibas, ishell,
288     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
289                  else if (type .eq. 14) then ! 3
290                     call int_angmom(jbas, jshell, ibas, ishell,
291     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf))
292                  else if (type .eq. 15) then ! 3
293                     call int_giaobq10(jbas, jshell, ibas, ishell,
294     $                    mem1, dbl_mb(k_scr), max1e, dbl_mb(k_buf),
295     &                    dbl_mb(k_qbq),dbl_mb(k_cbq),nbq)
296                  else
297                     call errquit('int_giao_1ega: invalid type?', type,
298     &       GA_ERR)
299                  end if
300c
301c     Add the integrals into the global array
302c
303                  alo(1) = jlo
304                  ahi(1) = jhi
305                  alo(2) = ilo
306                  ahi(2) = ihi
307                  alo(3) = 1
308                  ahi(3) = nblocks
309                  ld(1) = jlen
310                  ld(2) = ilen
311                  call nga_acc(g,alo,ahi,dbl_mb(k_buf),ld,1.0d0)
312               end if
313            endif
314            ijshell = ijshell + 1
315         end do
316      end do
317c
318c     chop stack at first item allocated
319c
320      if (type.eq.15) then  ! for bq's
321        if (.not. MA_pop_stack(l_cbq)) call errquit
322     $     ('int_giao_1ega: pop failed', l_cbq, MA_ERR)
323        if (.not. MA_pop_stack(l_qbq)) call errquit
324     $     ('int_giao_1ega: pop failed', l_qbq, MA_ERR)
325      end if ! type 15
326      if (.not. MA_pop_stack(l_scr)) call errquit
327     $     ('int_giao_1ega: pop failed', l_scr, MA_ERR)
328      if (.not. MA_pop_stack(l_buf)) call errquit
329     $     ('int_giao_1ega: pop failed', l_buf, MA_ERR)
330
331      call ga_sync()            ! So that no nasty races can result
332c
333      end
334