1      subroutine zora_getv_HFine_slow(
2     &                           rtdb,
3     &                           g_dens,        ! in  : atomic density
4     &                           ofinite,       ! in  : = .true. if Gaussian Nucl. Model requested
5     &                           atmass,        ! in  : atomic mass
6     &                           xyz_NMRcoords, ! in  : nuclear coordinates
7     &                           g_zpso,        ! out : ZPSO term
8     &                           g_fcsd,        ! out : FC+SD (v,u) term
9     &                           nexc)
10c
11C$Id$
12c Adapted from zora_getv_so
13
14      implicit none
15#include "rtdb.fh"
16#include "bas.fh"
17#include "cdft.fh"
18#include "errquit.fh"
19#include "mafdecls.fh"
20#include "global.fh"
21#include "geom.fh"
22#include "msgtypesf.h"
23#include "msgids.fh"
24#include "stdio.fh"
25#include "cgridfile.fh"
26#include "grid_cube.fh"
27#include "modelpotential.fh"
28c
29c     == arguments ==
30      integer rtdb
31      integer g_dens(2)
32      integer g_zpso(3),g_fcsd(3,3)
33      integer nexc
34c
35c     == local variables ==
36      integer i,j,k,t,n,ind,nij
37      double precision rho_n
38      double precision tmat
39      double precision dummy(2)
40      integer iqsh, istep, nxyz, ncontrset
41      integer ixyz, lxyz, icharge, lcharge, itags, ltags
42      integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
43     &     ictr_buf,iqpts
44      double precision rad,ke
45      integer lbas_cset_info, ibas_cset_info,
46     &        lbas_cent_info, ibas_cent_info,
47     &        ldocset, idocset,
48     &        l_rchi_atom,i_rchi_atom,
49     &        l_rq,i_rq,l_iniz, k_iniz,
50     &        lchi_ao, ichi_ao,
51     &        ldelchi_ao, idelchi_ao
52      integer lzpso(3),izpso(3),
53     &        lfcsd(3,3),ifcsd(3,3)
54      integer inntsize,ddblsize,ok
55      double precision xyz_NMRcoords(3),atmass
56      double precision chi_cntr(3,nbf_ao)
57
58      logical grid_file_rewind,ofinite
59      external grid_file_rewind,calc_zora_HFine_slow,
60     &         ga_antisymmetrize
61c
62c     model potential parameters
63      character*2 gelem(ncenters)
64      double precision gexpo(ncenters,50)
65      double precision gcoef(ncenters,50)
66c
67c     == allocate memory ==
68      do i=1,3
69       if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao,
70     &                      'lnmr',lzpso(i),izpso(i)))
71     &    call errquit('zora_getv_HFine: zpso',911,MA_ERR)
72       enddo
73       do i=1,3
74          do j=i,3
75             if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao,
76     &            'lfcsd',lfcsd(i,j),ifcsd(i,j)))
77     &            call errquit('zora_getv_HFine: fcsd',911,MA_ERR)
78          enddo
79      enddo ! end-loop-i
80c     == preliminaries ==
81      call dfill(3*nbf_ao*nbf_ao,0d0,dbl_mb(izpso(j)),1)
82      do k=1,3
83         do j=k,3
84            call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(ifcsd(k,j)),1)
85         enddo
86      enddo
87c
88c     get zora model potential parameters of geometry
89      if (use_modelpotential)
90     &  call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef)
91c
92c     == generate the grid ==
93      dummy(1) = 0.d0
94      dummy(2) = 0.d0
95      call grid_quadv0(rtdb,g_dens,g_zpso,nexc,rho_n,dummy,tmat)
96c     == ao basis set info used by xc_eval_basis ==
97      if (.not.bas_numcont(AO_bas_han, ncontrset))
98     &     call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR)
99      if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
100     &     lbas_cent_info, ibas_cent_info))
101     &     call errquit('zora_getv_sf: cannot allocate bas_cent_info',0,
102     &       MA_ERR)
103      if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
104     &     lbas_cset_info, ibas_cset_info))
105     &     call errquit('zora_getv_sf: cannot allocate bas_cset_info',0,
106     &       MA_ERR)
107      call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
108     &     int_mb(ibas_cset_info), ncenters)
109      if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
110     &     ldocset, idocset))
111     &     call errquit('zora_getv_sf: cannot allocate ccdocset',
112     .     ncontrset, MA_ERR)
113      do i=1,ncontrset
114         log_mb(idocset+i-1)=.true.
115      enddo
116      if(.not.MA_push_get(MT_int, ncenters, 'iniz',
117     &     l_iniz, k_iniz))
118     &     call errquit("zora_getv_sf:iniz",0, MA_ERR)
119      do i= 1, ncenters
120         int_mb(k_iniz+i-1)=1
121      enddo
122      nxyz = 3*ncenters
123      if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
124     &   call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR)
125      if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
126     &   call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR)
127      if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
128     &   call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR)
129      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
130     &                        Dbl_MB(ixyz), Dbl_MB(icharge)))
131     &   call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR)
132
133      if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
134     &   call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR)
135      if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
136     &   call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR)
137      if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
138     &     'quad pts buffer', lrqbuf, irqbuf))
139     &     call errquit('zora_getv_sf: quad buffer', 3, MA_ERR)
140
141      if (.not. grid_file_rewind())
142     $     call errquit('zora_getv_sf: rewinding gridpts?', 0,
143     &       UNKNOWN_ERR)
144c
145c     == loop over records in the grid file ==
146      do iqsh = 1, n_rec_in_file
147c
148c       == define the current range of radial shells and integration center ==
149        call grid_file_read(n_per_rec, nqpts, ictr_buf,
150     &        rad,dbl_mb(irqbuf),nsubb)
151
152        if(nqpts.gt.buffer_size)
153     &    call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
154c
155c        == loop over a subset of the grid ==
156         istep=0
157         do  ncube=1,nsubb
158c
159c           put buf into currently used arrays qxyz and qwght
160            call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
161     &           dbl_mb(iqwght), nqpts, rad,istep)
162
163            if(nqpts.ne.0) then
164c
165c              == compute the basis functions over the grid ==
166               if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
167     &             l_rchi_atom,i_rchi_atom))
168     &             call errquit("zora_getv:rchi_atom",0, MA_ERR)
169c
170               if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
171     &             l_rq,i_rq))
172     &             call errquit("zora_getv_sf:rq",0, MA_ERR)
173c
174c              == delchi ==
175               if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
176     &             'delchi_ao', ldelchi_ao, idelchi_ao))
177     &             call errquit('zora_getv: delchi_ao',0, MA_ERR)
178c
179c              == chi ==
180               if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
181     &             'chi_ao', lchi_ao, ichi_ao))
182     &             call errquit('zora_getv: chi_ao',0, MA_ERR)
183               call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
184     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
185               call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
186     &              dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
187     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
188     &              int_mb(k_iniz), log_mb(idocset),
189     &              int_mb(ibas_cent_info), int_mb(ibas_cset_info))
190                call calc_zora_HFine_slow(
191     &                                  ao_bas_han,geom,ipol,g_dens,
192     &                                  dbl_mb(ichi_ao),
193     &                                  dbl_mb(idelchi_ao),
194     &                                  dbl_mb(iqxyz),dbl_mb(iqwght),
195     &                                  nbf_ao,nqpts,ncenters,
196     &                                  ofinite,
197     &                                  atmass,
198     &                                  xyz_NMRcoords,
199     &                                  use_modelpotential,
200     &                                  gexpo,
201     &                                  gcoef,
202     &                                  dbl_mb(izpso(1)),  ! out
203     &                                  dbl_mb(izpso(2)),  ! out
204     &                                  dbl_mb(izpso(3)),  ! out
205     &                                  dbl_mb(ifcsd(1,1)),! out
206     &                                  dbl_mb(ifcsd(1,2)),! out
207     &                                  dbl_mb(ifcsd(1,3)),! out
208c     &                                  dbl_mb(ifcsd(2,1)),! out
209     &                                  dbl_mb(ifcsd(2,2)),! out
210     &                                  dbl_mb(ifcsd(2,3)),! out
211c     &                                  dbl_mb(ifcsd(3,1)),! out
212c     &                                  dbl_mb(ifcsd(3,2)),! out
213     &                                  dbl_mb(ifcsd(3,3)))! out
214
215c              == delete memory ==
216               if(.not.MA_chop_stack(l_rchi_atom))
217     &            call errquit("zora_getv: pop rchi_atom",100,MA_ERR)
218            endif ! nqpts
219         enddo ! ncube
220      end do ! iqsh
221c
222c     == delete memory ==
223      if(.not.MA_chop_stack(lbas_cent_info))
224     &     call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR)
225c
226c     == tally up over all the nodes ==
227      do k=1,3
228         call ga_dgop(msg_excrho+k,
229     D        dbl_mb(izpso(k)), nbf_ao*nbf_ao, '+')
230         do t=k,3
231c     do t=1,3
232          call ga_dgop(msg_excrho+t*k,
233     D           dbl_mb(ifcsd(k,t)),nbf_ao*nbf_ao, '+')
234         enddo                  ! end-loop-t
235      enddo                     ! end-loop-k
236c
237c     == pack into a ga all HFine AOs ==
238      if(ga_nodeid().eq.0) then
239      do i=1,3
240       call ga_zero(g_zpso(i))
241       call ga_put(g_zpso(i),
242     &             1,nbf_ao,1,nbf_ao,dbl_mb(izpso(i)),nbf_ao)
243       do j=i,3
244        call ga_zero(g_fcsd(i,j))
245        call ga_put(g_fcsd(i,j),
246     &              1,nbf_ao,1,nbf_ao,dbl_mb(ifcsd(i,j)),nbf_ao)
247        if(i.ne.j) then
248           call ga_copy(g_fcsd(i,j),g_fcsd(j,i))
249        endif
250       enddo ! end-loop-j
251      enddo ! end-loop-i
252      endif
253      call ga_sync()
254c ---- Free memory
255      if (.not.ma_chop_stack(lzpso(1))) call errquit
256     &     ('zora_getv_HFine: ma_chop_stack of lzpso failed',
257     &     911,MA_ERR)
258      return
259      end
260c ++++++++++++++++++++++++++++++++
261c ++++++++++++++++++++++++++++++++
262      subroutine zora_getv_HFine_fast(
263     &                           rtdb,
264     &                           g_dens,        ! in  : atomic density
265     &                           ofinite,       ! in  : = .true. if Gaussian Nucl. Model requested
266     &                           zetanuc_arr,   ! in  : sqrt(zetanuc(i)) i=1,natoms
267     &                           zetanuc_slc,   ! in  : zetanuc(i)
268     &                           Knucl,
269     &                           xyz_NMRcoords, ! in  : nuclear coordinates
270     &                           g_zpso,        ! out : ZPSO term
271     &                           g_fcsd,        ! out : FC+SD (v,u) term
272     &                           nexc)
273c
274C$Id$
275c Adapted from zora_getv_so
276
277      implicit none
278#include "rtdb.fh"
279#include "bas.fh"
280#include "cdft.fh"
281#include "errquit.fh"
282#include "mafdecls.fh"
283#include "global.fh"
284#include "geom.fh"
285#include "msgtypesf.h"
286#include "msgids.fh"
287#include "stdio.fh"
288#include "cgridfile.fh"
289#include "grid_cube.fh"
290#include "modelpotential.fh"
291c
292c     == arguments ==
293      integer rtdb
294      integer g_dens(2)
295      integer g_zpso(6),g_fcsd(3,3)
296      integer nexc
297c
298c     == local variables ==
299      integer i,j,k,t,m,n,ind,nij,ac
300      double precision rho_n
301      double precision tmat
302      double precision dummy(2)
303      integer iqsh, istep, nxyz, ncontrset
304      integer ixyz, lxyz, icharge, lcharge, itags, ltags
305      integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
306     &     ictr_buf,iqpts
307      double precision rad,ke
308      double precision zetanuc_arr(*) ! For Gaussian Nuclear Model
309      double precision zetanuc_slc    ! For Gaussian Nuclear Model
310      integer lbas_cset_info, ibas_cset_info,
311     &        lbas_cent_info, ibas_cent_info,
312     &        ldocset, idocset,
313     &        l_rchi_atom,i_rchi_atom,
314     &        l_rq,i_rq,
315     &        lchi_ao, ichi_ao,
316     &        ldelchi_ao, idelchi_ao
317      integer lzpso(3,3),izpso(3,3),
318     &        lfcsd(3,3),ifcsd(3,3)
319      integer inntsize,ddblsize,ok,atn
320      double precision xyz_NMRcoords(3),atmass
321      double precision chi_cntr(3,nbf_ao)
322      logical grid_file_rewind,ofinite,Knucl
323      integer ind_tmn(2,3)
324      data ind_tmn / 2, 3,  ! tmn=123
325     &               3, 1,  ! tmn=231
326     &               1, 2 / ! tmn=312
327      logical dft_mirrdens_start,dorepl
328      integer g_dens0(2),g_densrep(2),ii
329      external dft_mirrdens_start
330      external grid_file_rewind
331c
332c     model potential parameters
333      character*2 gelem(ncenters)
334      double precision gexpo(ncenters,50)
335      double precision gcoef(ncenters,50)
336c     mbf
337      integer grid_nbfm
338      external grid_nbfm
339      character*32 pname
340      double precision acc_ao_gauss
341      integer mbf_ao,k_expo,l_expo,l_ifin,k_ifin,l_iniz,k_iniz,
342     k     k_ibf_ao,l_ibf_ao
343      pname = 'zora_getv_hfine: '
344c
345c     == allocate memory ==
346       do t=1,3
347        m=ind_tmn(1,t)
348        n=ind_tmn(2,t)
349        if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao,
350     &                'lnmr',lzpso(m,n),izpso(m,n)))
351     &    call errquit('zora_getv_NMR: zpso',911,MA_ERR)
352        if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao,
353     &                'lnmr',lzpso(n,m),izpso(n,m)))
354     &    call errquit('zora_getv_NMR: zpso',911,MA_ERR)
355          call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(izpso(m,n)),1)
356          call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(izpso(n,m)),1)
357       enddo ! end-loop-t
358       do i=1,3
359        do j=i,3
360          if (.not. ma_push_get(mt_dbl,nbf_ao*nbf_ao,
361     &                      'lfcsd',lfcsd(i,j),ifcsd(i,j)))
362     &    call errquit('zora_getv_HFine: fcsd',911,MA_ERR)
363          call dfill(nbf_ao*nbf_ao,0d0,dbl_mb(ifcsd(i,j)),1)
364        enddo
365       enddo ! end-loop-i
366c
367c     get zora model potential parameters of geometry
368      if (use_modelpotential)
369     &  call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef)
370c
371c     == generate the grid ==
372      dummy(1) = 0.d0
373      dummy(2) = 0.d0
374      dorepl = dft_mirrdens_start(g_dens,g_densrep,g_dens0,
375     i     ipol)
376      call grid_quadv0(rtdb,g_dens,g_zpso,nexc,rho_n,dummy,tmat)
377c     == ao basis set info used by xc_eval_basis ==
378      if (.not.bas_numcont(AO_bas_han, ncontrset))
379     &     call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR)
380      if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
381     &     lbas_cent_info, ibas_cent_info))
382     &     call errquit('zora_getv_sf: cannot allocate bas_cent_info',0,
383     &       MA_ERR)
384      if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
385     &     lbas_cset_info, ibas_cset_info))
386     &     call errquit('zora_getv_sf: cannot allocate bas_cset_info',0,
387     &       MA_ERR)
388      call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
389     &     int_mb(ibas_cset_info), ncenters)
390      nxyz = 3*ncenters
391      if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
392     &   call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR)
393      if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
394     &   call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR)
395      if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
396     &   call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR)
397      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
398     &                        Dbl_MB(ixyz), Dbl_MB(icharge)))
399     &   call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR)
400
401      if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
402     &   call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR)
403      if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
404     &   call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR)
405      if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
406     &     'quad pts buffer', lrqbuf, irqbuf))
407     &     call errquit('zora_getv_sf: quad buffer', 3, MA_ERR)
408
409      if (.not. grid_file_rewind())
410     $     call errquit('zora_getv_sf: rewinding gridpts?', 0,
411     &       UNKNOWN_ERR)
412      if (.not.MA_Push_Get(MT_Dbl,nbf_ao_mxprim,'expo',l_expo,k_expo))
413     &   call errquit(pname//'cannot allocate expo',0, MA_ERR)
414      if (.not.MA_Push_Get(mt_int, nbf_ao, 'ibf_ao', l_ibf_ao,
415     &     k_ibf_ao))
416     &     call errquit(pname//'cannot allocate ibf_ao',2,
417     &       MA_ERR)
418      if (.not.MA_Push_get(MT_int,ncenters,'atom list',l_iniz,k_iniz))
419     &   call errquit(pname//'cannot allocate iniz',0, MA_ERR)
420      if (.not.MA_Push_get(MT_int,ncenters,'atom list',l_ifin,k_ifin))
421     &   call errquit(pname//'cannot allocate fin',0, MA_ERR)
422      do i= 1, ncenters
423         int_mb(k_iniz+i-1)=1
424      enddo
425      if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
426     &     ldocset, idocset))
427     &     call errquit('zora_getv_sf: cannot allocate ccdocset',
428     .     ncontrset, MA_ERR)
429      do i=1,ncontrset
430         log_mb(idocset+i-1)=.true.
431      enddo
432c
433c     == loop over records in the grid file ==
434      do iqsh = 1, n_rec_in_file
435c
436c       == define the current range of radial shells and integration center ==
437        call grid_file_read(n_per_rec, nqpts, ictr_buf,
438     &        rad,dbl_mb(irqbuf),nsubb)
439
440        if(nqpts.gt.buffer_size)
441     &    call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
442c
443c        == loop over a subset of the grid ==
444         istep=0
445         do  ncube=1,nsubb
446c
447c           put buf into currently used arrays qxyz and qwght
448            call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
449     &           dbl_mb(iqwght), nqpts, rad,istep)
450
451            if(nqpts.ne.0) then
452               mbf_ao=nbf_ao
453
454crestrict to mbf_ao subset
455c               acc_ao_gauss = dble(max(iaoacc,25))
456               acc_ao_gauss = dble(iaoacc)
457                 call icopy(mbf_ao, 0,0, int_mb(k_ibf_ao), 1)
458                  mbf_ao=grid_nbfm(ao_bas_han,  ncenters,
459     &                 ictr_buf,rad,
460     Q                 dbl_mb(ixyz),dbl_mb(iqxyz),nqpts,
461     &                 int_mb(k_ibf_ao), log_mb(idocset),
462     I                 int_mb(k_iniz), int_mb(k_ifin), dbl_mb(k_expo),
463     &                 minexp,ldiff,acc_ao_gauss,iatype_pt_chg)
464
465c
466c              == compute the basis functions over the grid ==
467               if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
468     &             l_rchi_atom,i_rchi_atom))
469     &             call errquit("zora_getv:rchi_atom",0, MA_ERR)
470c
471               if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
472     &             l_rq,i_rq))
473     &             call errquit("zora_getv_sf:rq",0, MA_ERR)
474c
475c              == delchi ==
476               if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
477     &             'delchi_ao', ldelchi_ao, idelchi_ao))
478     &             call errquit('zora_getv: delchi_ao',0, MA_ERR)
479c
480c              == chi ==
481               if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
482     &             'chi_ao', lchi_ao, ichi_ao))
483     &             call errquit('zora_getv: chi_ao',0, MA_ERR)
484               call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
485     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
486               call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
487     &              dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
488     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
489     &              int_mb(k_iniz), log_mb(idocset),
490     &              int_mb(ibas_cent_info), int_mb(ibas_cset_info))
491                call calc_zora_HFine_fast(
492     &                                ao_bas_han,geom,ipol,g_dens,
493     &                                dbl_mb(ichi_ao),
494     &                                dbl_mb(idelchi_ao),
495     &                                dbl_mb(iqxyz),dbl_mb(iqwght),
496     &                                nbf_ao,mbf_ao,int_mb(k_ibf_ao),
497     N              nqpts,ncenters,
498     &                                ofinite,
499     &                                zetanuc_arr,
500     &                                zetanuc_slc,
501     &                                Knucl,
502     &                                xyz_NMRcoords,
503     &                                use_modelpotential,
504     &                                gexpo,
505     &                                gcoef,
506     &                                dbl_mb(izpso(1,2)),! out
507     &                                dbl_mb(izpso(2,3)),! out
508     &                                dbl_mb(izpso(3,1)),! out
509     &                                dbl_mb(ifcsd(1,1)),! out
510     &                                dbl_mb(ifcsd(1,2)),! out
511     &                                dbl_mb(ifcsd(1,3)),! out
512c     &                                dbl_mb(ifcsd(2,1)),! out
513     &                                dbl_mb(ifcsd(2,2)),! out
514     &                                dbl_mb(ifcsd(2,3)),! out
515c     &                                dbl_mb(ifcsd(3,1)),! out
516c     &                                dbl_mb(ifcsd(3,2)),! out
517     &                                dbl_mb(ifcsd(3,3)))! out
518c              == delete memory ==
519               if(.not.MA_chop_stack(l_rchi_atom))
520     &            call errquit("zora_getv: pop rchi_atom",100,MA_ERR)
521            endif ! nqpts
522         enddo ! ncube
523      end do ! iqsh
524c
525c     == delete memory ==
526      if(.not.MA_chop_stack(lbas_cent_info))
527     &     call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR)
528      call ga_sync()
529      if(dorepl) then
530         do ii=1,ipol
531            call util_mirrstop(g_densrep(ii))
532            g_dens(ii)=g_dens0(ii)
533         enddo
534      endif
535c
536c     == tally up over all the nodes ==
537      do k=1,3
538         m=ind_tmn(1,k)
539         n=ind_tmn(2,k)
540         if(k.eq.1) then
541            call ga_mask_sync(.true.,.false.)
542         else
543            call ga_mask_sync(.false.,.false.)
544         endif
545         call ga_dgop(msg_excrho,dbl_mb(izpso(m,n)),nbf_ao*nbf_ao,'+')
546         call dcopy(
547     D          nbf_ao*nbf_ao,dbl_mb(izpso(n,m)),1,dbl_mb(izpso(m,n)),1)
548c         call ga_dgop(msg_excrho,dbl_mb(izpso(n,m)),nbf_ao*nbf_ao,'+')
549         do t=k,3
550         if(t.eq.3) then
551            call ga_mask_sync(.false.,.true.)
552         else
553            call ga_mask_sync(.false.,.false.)
554         endif
555           call ga_dgop(msg_excrho,dbl_mb(ifcsd(k,t)),nbf_ao*nbf_ao,'+')
556c ifcsd upper triangle gone
557c           if(t.ne.k) call dcopy(
558c     D          nbf_ao*nbf_ao,dbl_mb(ifcsd(k,t)),1,dbl_mb(ifcsd(t,k)),1)
559         enddo                  ! end-loop-t
560      enddo                     ! end-loop-k
561c
562c     == pack into a ga all HFine AOs ==
563      ac=1
564      do t=1,3
565       m=ind_tmn(1,t)
566       n=ind_tmn(2,t)
567       call ga_zero(g_zpso(ac))
568       if(ga_nodeid().eq.0)
569     c      call ga_put(g_zpso(ac),
570     &             1,nbf_ao,1,nbf_ao,dbl_mb(izpso(m,n)),nbf_ao)
571       ac=ac+1
572      enddo ! end-loop-t
573      do t=1,3
574       m=ind_tmn(2,t)
575       n=ind_tmn(1,t)
576       call ga_zero(g_zpso(ac))
577       if(ga_nodeid().eq.0)
578     c      call ga_put(g_zpso(ac),
579     &             1,nbf_ao,1,nbf_ao,dbl_mb(izpso(m,n)),nbf_ao)
580       ac=ac+1
581      enddo ! end-loop-t
582      call ga_sync()
583      do i=1,3
584       do j=i,3
585        call ga_zero(g_fcsd(i,j))
586       if(ga_nodeid().eq.0)
587     c       call ga_put(g_fcsd(i,j),
588     &              1,nbf_ao,1,nbf_ao,dbl_mb(ifcsd(i,j)),nbf_ao)
589       if(i.ne.j) call ga_copy(g_fcsd(i,j),g_fcsd(j,i))
590       enddo ! end-loop-j
591      enddo ! end-loop-i
592      call ga_sync()
593c ----- free memory ----------- START
594       m=ind_tmn(1,1)
595       n=ind_tmn(2,1)
596      if (.not.ma_chop_stack(lzpso(m,n))) call errquit
597     &   ('zora_getv_HFine: ma_chop_stack of lzpso failed',
598     &    911,MA_ERR)
599
600c ----- free memory ----------- END
601      return
602      end
603
604      subroutine zora_getv_NMRHFine_F1ji(
605     &          rtdb,
606     &          g_dens,       ! in : atomic density
607     &          g_hfine,      ! out:
608     &          natoms,       ! in: nr. atoms
609     &          ofinite,      ! in: = .true. if Gaussian Nucl. Model of charges requested
610     &          zetanuc_arr,  ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model
611     &          Knucl,        ! in: = .true. if K_ZORA(V=Nuclear pot. only)
612     &          nexc)
613c
614C$Id$
615c Adapted from zora_getv_sf
616
617      implicit none
618#include "rtdb.fh"
619#include "bas.fh"
620#include "cdft.fh"
621#include "errquit.fh"
622#include "mafdecls.fh"
623#include "global.fh"
624#include "geom.fh"
625#include "msgtypesf.h"
626#include "msgids.fh"
627#include "stdio.fh"
628#include "cgridfile.fh"
629#include "grid_cube.fh"
630#include "modelpotential.fh"
631c
632c     == arguments ==
633      integer rtdb
634      integer g_dens(2)
635      integer g_hfine(3)
636      integer nexc
637c     == local variables ==
638      integer i,j,k,t,a,m,n,ind,nij,ac
639      double precision rho_n
640      double precision tmat
641      double precision dummy(2)
642      integer iqsh, istep, nxyz, ncontrset
643      integer ixyz, lxyz, icharge, lcharge, itags, ltags
644      integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
645     &     ictr_buf,iqpts
646      double precision rad,ke
647      integer lbas_cset_info, ibas_cset_info,
648     &        lbas_cent_info, ibas_cent_info,
649     &        ldocset, idocset,
650     &        l_rchi_atom,i_rchi_atom,
651     &        l_rq,i_rq,lniz, iniz,
652     &        lchi_ao, ichi_ao,
653     &        ldelchi_ao, idelchi_ao
654      double precision xyz_NMRcoords(3)
655      integer lhfine(3)   ,ihfine(3)
656      integer inntsize,ddblsize,ok
657      logical grid_file_rewind
658      integer natoms
659      logical ofinite,Knucl
660      double precision zetanuc_arr(natoms)
661      external grid_file_rewind,calc_NMRHFine_F1ij,
662     &         ga_antisymmetrize,grid_quadv0
663c
664c     model potential parameters
665      character*2 gelem(ncenters)
666      double precision gexpo(ncenters,50)
667      double precision gcoef(ncenters,50)
668c
669c     == allocate memory ==
670      do i=1,3
671      if (.not. ma_alloc_get(mt_dbl,nbf_ao*nbf_ao,
672     &                      'lhfine',lhfine(i),ihfine(i)))
673     &    call errquit('zora_getv_NMR: hfine',911,MA_ERR)
674      enddo ! end-loop-i
675c     == preliminaries ==
676       do j=1,3
677         do i= 1, nbf_ao*nbf_ao
678          dbl_mb(ihfine(j)+i-1) =0.d0
679         enddo
680       enddo
681c
682c     get zora model potential parameters of geometry
683      if (use_modelpotential)
684     &  call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef)
685c
686c     == generate the grid ==
687      dummy(1) = 0.d0
688      dummy(2) = 0.d0
689      call grid_quadv0(rtdb,g_dens,g_hfine(1),nexc,rho_n,dummy,tmat)
690c     == ao basis set info used by xc_eval_basis ==
691      if (.not.bas_numcont(AO_bas_han, ncontrset))
692     &     call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR)
693      if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
694     &     lbas_cent_info, ibas_cent_info))
695     &     call errquit('zora_getv_sf: cannot allocate bas_cent_info',0,
696     &       MA_ERR)
697      if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
698     &     lbas_cset_info, ibas_cset_info))
699     &     call errquit('zora_getv_sf: cannot allocate bas_cset_info',0,
700     &       MA_ERR)
701      call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
702     &     int_mb(ibas_cset_info), ncenters)
703      if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
704     &     ldocset, idocset))
705     &     call errquit('zora_getv_sf: cannot allocate ccdocset',
706     .     ncontrset, MA_ERR)
707      do i=1,ncontrset
708         log_mb(idocset+i-1)=.true.
709      enddo
710      if(.not.MA_push_get(MT_int, ncenters, 'iniz',
711     &     lniz, iniz))
712     &     call errquit("zora_getv_sf:iniz",0, MA_ERR)
713      do i= 1, ncenters
714         int_mb(iniz+i-1)=1
715      enddo
716      nxyz = 3*ncenters
717      if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
718     &   call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR)
719      if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
720     &   call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR)
721      if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
722     &   call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR)
723      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
724     &                        Dbl_MB(ixyz), Dbl_MB(icharge)))
725     &   call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR)
726
727      if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
728     &   call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR)
729      if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
730     &   call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR)
731      if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
732     &     'quad pts buffer', lrqbuf, irqbuf))
733     &     call errquit('zora_getv_sf: quad buffer', 3, MA_ERR)
734
735      if (.not. grid_file_rewind())
736     $     call errquit('zora_getv_sf: rewinding gridpts?', 0,
737     &       UNKNOWN_ERR)
738c
739c     == loop over records in the grid file ==
740      do iqsh = 1, n_rec_in_file
741c       == define the current range of radial shells and integration center ==
742        call grid_file_read(n_per_rec, nqpts, ictr_buf,
743     &        rad,dbl_mb(irqbuf),nsubb)
744        if(nqpts.gt.buffer_size)
745     &    call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
746c        == loop over a subset of the grid ==
747         istep=0
748         do  ncube=1,nsubb
749c        put buf into currently used arrays qxyz and qwght
750            call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
751     &           dbl_mb(iqwght), nqpts, rad,istep)
752            if(nqpts.ne.0) then
753c        == compute the basis functions over the grid ==
754               if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
755     &             l_rchi_atom,i_rchi_atom))
756     &             call errquit("zora_getv:rchi_atom",0, MA_ERR)
757               if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
758     &             l_rq,i_rq))
759     &             call errquit("zora_getv_sf:rq",0, MA_ERR)
760c              == delchi ==
761               if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
762     &             'delchi_ao', ldelchi_ao, idelchi_ao))
763     &             call errquit('zora_getv: delchi_ao',0, MA_ERR)
764c              == chi ==
765               if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
766     &             'chi_ao', lchi_ao, ichi_ao))
767     &             call errquit('zora_getv: chi_ao',0, MA_ERR)
768               call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
769     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
770               call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
771     &              dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
772     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
773     &              int_mb(iniz), log_mb(idocset),
774     &              int_mb(ibas_cent_info), int_mb(ibas_cset_info))
775               call calc_NMRHFine_F1ij(ao_bas_han,geom,ipol,g_dens,
776     &                                 dbl_mb(idelchi_ao),
777     &                                 dbl_mb(iqxyz),dbl_mb(iqwght),
778     &                                 nbf_ao,nqpts,ncenters,
779     &                                 ofinite,      ! in: = .true. if Gaussian Nucl. Model of charges requested
780     &                                 zetanuc_arr,  ! in: sqrt(zetanuc(i)) i=1,natoms for Gaussian Nuclear Model
781     &                                 Knucl,
782     &                                 use_modelpotential,
783     &                                 gexpo,
784     &                                 gcoef,
785     &                                 dbl_mb(ihfine(1)),  ! out
786     &                                 dbl_mb(ihfine(2)),  ! out
787     &                                 dbl_mb(ihfine(3)))  ! out
788c              == delete memory ==
789               if(.not.MA_pop_stack(lchi_ao))
790     &            call errquit("zora_getv: pop chi_ao", 100, MA_ERR)
791               if(.not.MA_pop_stack(ldelchi_ao))
792     &            call errquit("zora_getv: pop delchi_ao", 100, MA_ERR)
793               if(.not.MA_pop_stack(l_rq))
794     &            call errquit("zora_getv: pop rq", 100, MA_ERR)
795               if(.not.MA_pop_stack(l_rchi_atom))
796     &            call errquit("zora_getv: pop rchi_atom",100,MA_ERR)
797            endif ! nqpts
798         enddo ! ncube
799      end do ! iqsh
800c     == delete memory ==
801      if(.not.MA_pop_stack(lrqbuf))
802     &     call errquit("zora_getv_sf: pop rqbuf", 100, MA_ERR)
803      if(.not.MA_pop_stack(lqwght))
804     &     call errquit("zora_getv_sf: pop qwght", 100, MA_ERR)
805      if(.not.MA_pop_stack(lqxyz))
806     &     call errquit("zora_getv_sf: pop qxyz", 100, MA_ERR)
807      if(.not.MA_pop_stack(ltags))
808     &     call errquit("zora_getv_sf: pop tags", 100, MA_ERR)
809      if(.not.MA_pop_stack(lcharge))
810     &     call errquit("zora_getv_sf: pop charge", 100, MA_ERR)
811      if(.not.MA_pop_stack(lxyz))
812     &     call errquit("zora_getv_sf: pop xyz", 100, MA_ERR)
813      if(.not.MA_pop_stack(lniz))
814     &     call errquit("zora_getv_sf: pop niz", 100, MA_ERR)
815      if(.not.MA_pop_stack(ldocset))
816     &     call errquit("zora_getv_sf: pop docset", 100, MA_ERR)
817      if(.not.MA_pop_stack(lbas_cset_info))
818     &     call errquit("zora_getv_sf: pop bas_cset_info", 100, MA_ERR)
819      if(.not.MA_pop_stack(lbas_cent_info))
820     &     call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR)
821c
822c     == tally up over all the nodes ==
823      do k=1,3
824         call ga_dgop(msg_excrho,dbl_mb(ihfine(k)), nbf_ao*nbf_ao, '+')
825      enddo                     ! end-loop-k
826c
827c     == pack into a ga all NMR AOs ==
828      do i=1,3
829       call ga_zero(g_hfine(i))
830       call ga_put(g_hfine(i),
831     &             1,nbf_ao,1,nbf_ao,dbl_mb(ihfine(i)),nbf_ao)
832       call ga_antisymmetrize(g_hfine(i))
833      enddo ! end-loop-i
834      call ga_sync()
835c ----- free memory ----------- START
836      do i=1,3
837      if (.not.ma_free_heap(lhfine(i))) call errquit
838     &   ('zora_getv_NMR: ma_free_heap of lhfine failed',911,MA_ERR)
839      enddo
840c ----- free memory ----------- END
841      return
842      end
843      logical function dft_mirrdens_start(g_dens,g_densrep,g_dens0,
844     i     ipol)
845      implicit none
846#include "global.fh"
847      integer g_dens(2) ! [in/out]
848      integer g_dens0(2),g_densrep(2) ! [in/out]
849      integer ipol ! [in]
850c
851      integer dorep_glob,ii
852      logical dorepd,dorepl,docopy,dozero
853      logical util_mirrmat
854      external util_mirrmat
855      dorepl=.false.
856      call ga_sync()
857      if(ga_cluster_nnodes().gt.1) then
858         docopy=.true.
859         dozero=.false.
860         dorepd=util_mirrmat(ipol,g_dens,g_densrep,docopy,dozero)
861         dorep_glob=0
862         if(dorepd) dorep_glob=1
863         call ga_igop(375,dorep_glob,1, '+')
864         dorepl=dorep_glob.eq.ga_nnodes()
865         if(dorepl) then
866            do ii=1,ipol
867               g_dens0(ii)=g_dens(ii)
868               g_dens(ii)=g_densrep(ii)
869            enddo
870         else
871            if(dorepd) then
872               call util_mirrstop(g_densrep(1))
873               if(ipol.eq.2) call util_mirrstop(g_densrep(2))
874            endif
875            if(ga_nodeid().eq.0) then
876               write(6,*) ' no DM mirroring in zora_getv'
877               call util_flush(6)
878            endif
879         endif
880      endif
881      dft_mirrdens_start=dorepl
882      return
883      end
884