1c
2c     == compute the matrix elements for the zora corrections ==
3      subroutine zora_getv_sf(rtdb, g_dens,
4     &                        g_zora_sf, g_zora_scale_sf,
5     &                        ofinite,zetanuc_arr,
6     &                        Knucl, nexc)
7c
8C$Id$
9c Modified from zora_getv_sf() by FA 10-31-10
10c
11      implicit none
12c
13#include "rtdb.fh"
14#include "bas.fh"
15#include "cdft.fh"
16#include "errquit.fh"
17#include "mafdecls.fh"
18#include "global.fh"
19#include "geom.fh"
20#include "msgtypesf.h"
21#include "msgids.fh"
22#include "stdio.fh"
23#include "cgridfile.fh"
24#include "grid_cube.fh"
25#include "modelpotential.fh"
26c
27c     == arguments ==
28      integer rtdb
29      integer g_dens(2)
30      integer g_zora_sf(2)
31      integer g_zora_scale_sf(2)
32      logical ofinite
33      double precision zetanuc_arr(*)!for Gaussian Nuclear Model
34      logical Knucl
35      integer nexc
36c
37c     == local variables ==
38      integer i,j,nij
39      double precision rho_n
40      double precision tmat
41      double precision dummy(2)
42      integer ilo,ihi,jlo,jhi
43      integer iqsh, istep, nxyz, ncontrset
44      integer ixyz, lxyz, icharge, lcharge, itags, ltags
45      integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
46     &     ictr_buf,iqpts
47c
48      double precision rad,ke
49      integer lbas_cset_info, ibas_cset_info,
50     &     lbas_cent_info, ibas_cent_info,
51     &     ldocset, idocset,
52     &     l_rchi_atom,i_rchi_atom,
53     &     l_rq,i_rq,lniz, iniz,
54     &     lchi_ao, ichi_ao,
55     &     ldelchi_ao, idelchi_ao,
56     &     lzora0, izora0,
57     &     lzora0scal, izora0scal
58c
59      integer inntsize,ddblsize,ok
60c
61      logical dorepl
62      integer g_densrep(2),g_dens0(2),dorep_glob
63      logical docopy,dozero,dorepd,dorepon
64      logical util_mirrmat
65      external util_mirrmat
66c
67      logical grid_file_rewind
68      external grid_file_rewind
69c
70      logical int_normalize
71      external int_normalize
72c
73c     model potential parameters
74      character*2 gelem(ncenters)
75      double precision gexpo(ncenters,50)
76      double precision gcoef(ncenters,50)
77c
78c     == allocate memory ==
79      if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao,'zora0',lzora0,izora0))
80     &   call errquit('zora_getv: zorasf',0, MA_ERR)
81      if (.not.MA_Push_Get(mt_dbl, nbf_ao*nbf_ao,'zora0scal',
82     &  lzora0scal, izora0scal))
83     &       call errquit('zora_getv: zorasf',0, MA_ERR)
84c
85c     == preliminaries ==
86      do i= 1, nbf_ao*nbf_ao
87         dbl_mb(izora0+i-1)=0.d0
88         dbl_mb(izora0scal+i-1)=0.d0
89      enddo
90c
91c mirroring
92      dorepon=.true.
93      if (.not.rtdb_get(rtdb,'fock:mirrmat',mt_log,1,dorepon)) then
94         dorepon=nbf_ao.lt.3000
95      endif
96c
97      dorepl=.false.
98      if(ga_cluster_nnodes().gt.1.and.dorepon) then
99         docopy=.true.
100         dozero=.false.
101         dorepd=util_mirrmat(ipol,g_dens,g_densrep,docopy,dozero)
102         dorepl=dorepd
103         dorep_glob=0
104         if(dorepl) dorep_glob=1
105         call ga_igop(375,dorep_glob,1, '+')
106         dorepl=dorep_glob.eq.ga_nnodes()
107         if(dorepl) then
108            g_dens0(1)=g_dens(1)
109            g_dens(1)=g_densrep(1)
110            if(ipol.eq.2) then
111               g_dens0(2)=g_dens(2)
112               g_dens(2)=g_densrep(2)
113            endif
114         else
115            if(dorepd) then
116               call util_mirrstop(g_densrep(1))
117               if(ipol.eq.2) call util_mirrstop(g_densrep(2))
118            endif
119            if(ga_nodeid().eq.0) then
120               write(6,*) ' no mirroring in zora_getv_so'
121               call util_flush(6)
122            endif
123         endif
124      endif
125c
126c     get zora model potential parameters of geometry
127      if (use_modelpotential)
128     &  call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef)
129c
130c     == generate the grid ==
131      dummy(1) = 0.d0
132      dummy(2) = 0.d0
133      call grid_quadv0(rtdb, g_dens, g_zora_sf, nexc, rho_n, dummy,
134     &  tmat)
135c
136c     == initialization ==
137      call ga_zero(g_zora_sf(1))
138      if (ipol.gt.1) call ga_zero(g_zora_sf(2))
139      call ga_zero(g_zora_scale_sf(1))
140      if (ipol.gt.1) call ga_zero(g_zora_scale_sf(2))
141c
142c     == ao basis set info used by xc_eval_basis ==
143      if(.not.int_normalize(rtdb,AO_bas_han))
144     &  call errquit('zora_getv_sf: int_normalize failed',0, BASIS_ERR)
145      if (.not.bas_numcont(AO_bas_han, ncontrset))
146     &     call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR)
147      if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
148     &     lbas_cent_info, ibas_cent_info))
149     &     call errquit('zora_getv_sf: cannot allocate bas_cent_info',0,
150     &       MA_ERR)
151      if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
152     &     lbas_cset_info, ibas_cset_info))
153     &     call errquit('zora_getv_sf: cannot allocate bas_cset_info',0,
154     &       MA_ERR)
155      call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
156     &     int_mb(ibas_cset_info), ncenters)
157c
158      if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
159     &     ldocset, idocset))
160     &     call errquit('zora_getv_sf: cannot allocate ccdocset',
161     .     ncontrset, MA_ERR)
162      do i=1,ncontrset
163         log_mb(idocset+i-1)=.true.
164      enddo
165c
166      if (.not.MA_push_get(MT_int, ncenters, 'iniz',
167     &  lniz, iniz)) call errquit("zora_getv_sf:iniz",0, MA_ERR)
168      do i= 1, ncenters
169         int_mb(iniz+i-1)=1
170      enddo
171c
172      nxyz = 3*ncenters
173      if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
174     &   call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR)
175      if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
176     &   call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR)
177      if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
178     &   call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR)
179      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
180     &                        Dbl_MB(ixyz), Dbl_MB(icharge)))
181     &   call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR)
182
183      if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
184     &   call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR)
185      if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
186     &   call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR)
187      if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
188     &     'quad pts buffer', lrqbuf, irqbuf))
189     &     call errquit('zora_getv_sf: quad buffer', 3, MA_ERR)
190c
191c     == rewind grid file ==
192      if (.not. grid_file_rewind())
193     $     call errquit('zora_getv_sf: rewinding gridpts?', 0,
194     &       UNKNOWN_ERR)
195c     == loop over records in the grid file ==
196      do iqsh = 1, n_rec_in_file
197c
198c       == define the current range of radial shells and integration center ==
199        call grid_file_read(n_per_rec, nqpts, ictr_buf,
200     &        rad,dbl_mb(irqbuf),nsubb)
201
202        if(nqpts.gt.buffer_size)
203     &    call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
204c
205c        == loop over a subset of the grid ==
206         istep=0
207         do  ncube=1,nsubb
208c
209c           put buf into currently used arrays qxyz and qwght
210            call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
211     &           dbl_mb(iqwght), nqpts, rad,istep)
212
213            if(nqpts.ne.0) then
214c
215c              == compute the basis functions over the grid ==
216               if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
217     &             l_rchi_atom,i_rchi_atom))
218     &             call errquit("zora_getv:rchi_atom",0, MA_ERR)
219c
220               if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
221     &             l_rq,i_rq))
222     &             call errquit("zora_getv_sf:rq",0, MA_ERR)
223c
224c              == delchi ==
225               if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
226     &             'delchi_ao', ldelchi_ao, idelchi_ao))
227     &             call errquit('zora_getv: delchi_ao',0, MA_ERR)
228c
229c              == chi ==
230               if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
231     &             'chi_ao', lchi_ao, ichi_ao))
232     &             call errquit('zora_getv: chi_ao',0, MA_ERR)
233c
234               call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
235     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
236               call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
237     &              dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
238     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
239     &              int_mb(iniz), log_mb(idocset),
240     &              int_mb(ibas_cent_info), int_mb(ibas_cset_info))
241c
242c              calculate numerical overlap
243c               call num_ovlp(dbl_mb(ichi_ao), dbl_mb(iqwght),
244c     &           nbf_ao, nqpts, dbl_mb(izora0))
245c
246c              == calculate spin-free zora contribution ==
247               call calc_zora_sf(ao_bas_han, geom, ipol, g_dens,
248     &              dbl_mb(ichi_ao), dbl_mb(idelchi_ao), dbl_mb(iqxyz),
249     &              dbl_mb(iqwght), nbf_ao, nqpts, ncenters,
250     &              dbl_mb(izora0), dbl_mb(izora0scal),
251     &              ofinite,zetanuc_arr,
252     &              Knucl,  ! = .true. if including ONLY nuclear part in K ZORA
253     &              use_modelpotential,gexpo,gcoef)
254c
255c              == delete memory ==
256               if(.not.MA_pop_stack(lchi_ao))
257     &            call errquit("zora_getv: pop chi_ao", 100, MA_ERR)
258               if(.not.MA_pop_stack(ldelchi_ao))
259     &            call errquit("zora_getv: pop delchi_ao", 100, MA_ERR)
260               if(.not.MA_pop_stack(l_rq))
261     &            call errquit("zora_getv: pop rq", 100, MA_ERR)
262               if(.not.MA_pop_stack(l_rchi_atom))
263     &            call errquit("zora_getv: pop rchi_atom",100,MA_ERR)
264            endif ! nqpts
265         enddo ! ncube
266      end do ! iqsh
267c
268c     == delete memory ==
269      if(.not.MA_pop_stack(lrqbuf))
270     &     call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR)
271      if(.not.MA_pop_stack(lqwght))
272     &     call errquit("zora_getv_sf: pop qwght", 100, MA_ERR)
273      if(.not.MA_pop_stack(lqxyz))
274     &     call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR)
275      if(.not.MA_pop_stack(ltags))
276     &     call errquit("zora_getv_sf: pop tags", 100, MA_ERR)
277      if(.not.MA_pop_stack(lcharge))
278     &     call errquit("zora_getv_sf: pop charge", 100, MA_ERR)
279      if(.not.MA_pop_stack(lxyz))
280     &     call errquit("zora_getv_sf: pop xyz", 100, MA_ERR)
281      if(.not.MA_pop_stack(lniz))
282     &     call errquit("zora_getv_sf: pop niz", 100, MA_ERR)
283      if(.not.MA_pop_stack(ldocset))
284     &     call errquit("zora_getv_sf: pop docset", 100, MA_ERR)
285      if(.not.MA_pop_stack(lbas_cset_info))
286     &     call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR)
287      if(.not.MA_pop_stack(lbas_cent_info))
288     &     call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR)
289c
290      if(dorepl) then
291         call util_mirrstop(g_densrep(1))
292         if(ipol.eq.2) call util_mirrstop(g_densrep(2))
293         g_dens(1)=g_dens0(1)
294         if(ipol.eq.2) g_dens(2)=g_dens0(2)
295      endif
296c
297c     == tally up over all the nodes ==
298      call ga_dgop(msg_excrho, dbl_mb(izora0), nbf_ao*nbf_ao, '+')
299      call ga_dgop(msg_excrho, dbl_mb(izora0scal), nbf_ao*nbf_ao, '+')
300c
301c     == scalar contribution ==
302      call ga_zero(g_zora_sf(1))
303      call ga_distribution(g_zora_sf(1),
304     &     ga_nodeid(), ilo, ihi, jlo, jhi)
305      if (ilo.gt.0 .and. ilo.le.ihi) then
306         do j=jlo,jhi
307            call ga_put(g_zora_sf(1),ilo, ihi, jlo, jhi,
308     &           dbl_mb(izora0+(jlo-1)*nbf_ao+ilo-1),nbf_ao)
309         enddo
310      endif
311      call ga_symmetrize(g_zora_sf(1))
312c
313      call ga_zero(g_zora_scale_sf(1))
314      call ga_distribution(g_zora_scale_sf(1),
315     &     ga_nodeid(), ilo, ihi, jlo, jhi)
316      if (ilo.gt.0 .and. ilo.le.ihi) then
317         do j=jlo,jhi
318            call ga_put(g_zora_scale_sf(1),ilo, ihi, jlo, jhi,
319     &           dbl_mb(izora0scal+(jlo-1)*nbf_ao+ilo-1),nbf_ao)
320         enddo
321      endif
322      call ga_symmetrize(g_zora_scale_sf(1))
323c
324      if(ipol.eq.2) then
325         call ga_copy(g_zora_sf(1),g_zora_sf(2))
326         call ga_copy(g_zora_scale_sf(1),g_zora_scale_sf(2))
327      endif
328c
329c     call ga_print(g_zora_sf(1))
330c
331      call ga_sync()
332c
333      return
334      end
335c
336      subroutine num_ovlp(chi_ao,wght,nbf,npts,ovlp)
337
338      implicit none
339
340#include "stdio.fh"
341
342      integer nbf, npts
343      double precision chi_ao(npts,nbf),wght(npts)
344      double precision ovlp(nbf,nbf)
345      integer i, j, k
346
347      do i = 1, nbf
348        do j = 1, nbf
349          do k = 1, npts
350            ovlp(i,j)=ovlp(i,j)+chi_ao(k,i)*wght(k)*chi_ao(k,j)
351          enddo
352        enddo
353      enddo
354c
355      return
356      end
357