1c
2c     == EFGZ4 scalar and spin-orbit components (x,y,z)
3      subroutine zora_getv_EFGZ4_SO(rtdb,
4     &                              g_dens,        ! in : atomic density
5     &                              zora_calc_type,! in : =4 EFG-NUM, =3 EFGZ4-SR
6     &                              zora_Qpq,      ! in : type of EFG potential
7     &                              xyz_EFGcoords, ! in : EFG-nuclear coordinates
8     &                              g_efgz4_sf,    ! out: munu matrix
9     &                              g_efgz4_so,    ! out: munu matrix
10     &                              nexc)
11c
12C$Id$
13c Adapted from zora_getv_so
14
15      implicit none
16#include "rtdb.fh"
17#include "bas.fh"
18#include "cdft.fh"
19#include "errquit.fh"
20#include "mafdecls.fh"
21#include "global.fh"
22#include "geom.fh"
23#include "msgtypesf.h"
24#include "msgids.fh"
25#include "stdio.fh"
26#include "cgridfile.fh"
27#include "grid_cube.fh"
28#include "modelpotential.fh"
29c
30c     == arguments ==
31      integer rtdb
32      integer g_dens(2)
33      integer g_efgz4_sf(2)
34      integer g_efgz4_so(3)
35      integer zora_Qpq,zora_calc_type
36      double precision xyz_EFGcoords(3)
37      integer nexc
38c
39c     == local variables ==
40      integer i,j,k,ind,nij
41      double precision rho_n
42      double precision tmat
43      double precision dummy(2)
44      integer iqsh, istep, nxyz, ncontrset
45      integer ixyz, lxyz, icharge, lcharge, itags, ltags
46      integer lrqbuf,irqbuf,lqxyz,iqxyz,lqwght,iqwght,nqpts,ncube,
47     &     ictr_buf,iqpts
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     &        lscal0, iscal0,
57     &        lscalx, iscalx,
58     &        lscaly, iscaly,
59     &        lscalz, iscalz
60      integer lscal(4),iscal(4)
61      integer inntsize,ddblsize,ok
62      logical grid_file_rewind
63      external grid_file_rewind,calc_zora_EFGZ4_SO,calc_EFG,
64     &         ga_antisymmetrize
65c
66c     model potential parameters
67      character*2 gelem(ncenters)
68      double precision gexpo(ncenters,50)
69      double precision gcoef(ncenters,50)
70c
71c     == allocate memory ==
72      do i=1,4
73       if (.not.MA_Push_Get(mt_dbl,nbf_ao*nbf_ao,
74     &                     'scal0',lscal(i),iscal(i)))
75     &    call errquit('zora_getv_so: scali',0, MA_ERR)
76      enddo
77c     == preliminaries ==
78      do i= 1, nbf_ao*nbf_ao
79       do j=1,4
80         dbl_mb(iscal(j)+i-1)=0.d0
81       enddo
82      enddo
83c
84c     get zora model potential parameters of geometry
85      if (use_modelpotential)
86     &  call get_modelpotential_params(rtdb,ncenters,gelem,gexpo,gcoef)
87c
88c     == generate the grid ==
89      dummy(1) = 0.d0
90      dummy(2) = 0.d0
91      call grid_quadv0(rtdb,g_dens,g_efgz4_sf,nexc,rho_n,dummy,tmat)
92c     == ao basis set info used by xc_eval_basis ==
93      if (.not.bas_numcont(AO_bas_han, ncontrset))
94     &     call errquit('zora_getv_sf:bas_numcont',0, BASIS_ERR)
95      if (.not.MA_Push_Get(mt_int, 3*ncenters, 'bas_cent_info',
96     &     lbas_cent_info, ibas_cent_info))
97     &     call errquit('zora_getv_sf: cannot allocate bas_cent_info',0,
98     &       MA_ERR)
99      if (.not.MA_Push_Get(mt_int, 6*ncontrset, 'bas_cset_info',
100     &     lbas_cset_info, ibas_cset_info))
101     &     call errquit('zora_getv_sf: cannot allocate bas_cset_info',0,
102     &       MA_ERR)
103      call xc_make_basis_info(AO_bas_han, int_mb(ibas_cent_info),
104     &     int_mb(ibas_cset_info), ncenters)
105      if (.not.MA_Push_Get(mt_log, ncontrset, 'docset',
106     &     ldocset, idocset))
107     &     call errquit('zora_getv_sf: cannot allocate ccdocset',
108     .     ncontrset, MA_ERR)
109      do i=1,ncontrset
110         log_mb(idocset+i-1)=.true.
111      enddo
112c
113      if(.not.MA_push_get(MT_int, ncenters, 'iniz',
114     &     lniz, iniz))
115     &     call errquit("zora_getv_sf:iniz",0, MA_ERR)
116      do i= 1, ncenters
117         int_mb(iniz+i-1)=1
118      enddo
119c
120      nxyz = 3*ncenters
121      if (.not.MA_push_Get(MT_Dbl,nxyz,'xyz',lxyz,ixyz))
122     &   call errquit('zora_getv_sf: cannot allocate xyz',0, MA_ERR)
123      if (.not.MA_Push_Get(MT_Dbl,ncenters,'charge',lcharge,icharge))
124     &   call errquit('zora_getv_sf: cannot allocate charge',0, MA_ERR)
125      if (.not.MA_Push_Get(MT_Byte,ncenters*16,'tags',ltags,itags))
126     &   call errquit('zora_getv_sf: cannot allocate tags',0, MA_ERR)
127      if (.not. geom_cart_get(geom, ncenters, Byte_MB(itags),
128     &                        Dbl_MB(ixyz), Dbl_MB(icharge)))
129     &   call errquit('zora_getv_sf: geom_cart_get failed',74, GEOM_ERR)
130
131      if (.not.MA_Push_get(mt_dbl,3*n_per_rec,'qxyz',lqxyz,iqxyz))
132     &   call errquit('zora_getv_sf: cannot allocate qxyz',0, MA_ERR)
133      if (.not.MA_Push_get(mt_dbl,n_per_rec,'qwght',lqwght,iqwght))
134     &   call errquit('zora_getv_sf: cannot allocate qwght',0, MA_ERR)
135      if (.not.MA_Push_get(MT_dbl, 4*buffer_size+4,
136     &     'quad pts buffer', lrqbuf, irqbuf))
137     &     call errquit('zora_getv_sf: quad buffer', 3, MA_ERR)
138
139      if (.not. grid_file_rewind())
140     $     call errquit('zora_getv_sf: rewinding gridpts?', 0,
141     &       UNKNOWN_ERR)
142c
143c     == loop over records in the grid file ==
144      do iqsh = 1, n_rec_in_file
145c
146c       == define the current range of radial shells and integration center ==
147        call grid_file_read(n_per_rec, nqpts, ictr_buf,
148     &        rad,dbl_mb(irqbuf),nsubb)
149
150        if(nqpts.gt.buffer_size)
151     &    call errquit(' buffersize exceed by qpts ',nqpts, UNKNOWN_ERR)
152c
153c        == loop over a subset of the grid ==
154         istep=0
155         do  ncube=1,nsubb
156c
157c           put buf into currently used arrays qxyz and qwght
158            call grid_repack(dbl_mb(irqbuf), dbl_mb(iqxyz),
159     &           dbl_mb(iqwght), nqpts, rad,istep)
160
161            if(nqpts.ne.0) then
162c
163c              == compute the basis functions over the grid ==
164               if(.not.MA_Push_get(MT_dbl, ncenters, 'rchi_atom',
165     &             l_rchi_atom,i_rchi_atom))
166     &             call errquit("zora_getv:rchi_atom",0, MA_ERR)
167c
168               if(.not.MA_Push_get(MT_dbl, nqpts*ncenters, 'rq',
169     &             l_rq,i_rq))
170     &             call errquit("zora_getv_sf:rq",0, MA_ERR)
171c
172c              == delchi ==
173               if (.not.MA_Push_Get(mt_dbl, 3*nqpts*nbf_ao,
174     &             'delchi_ao', ldelchi_ao, idelchi_ao))
175     &             call errquit('zora_getv: delchi_ao',0, MA_ERR)
176c
177c              == chi ==
178               if (.not.MA_Push_Get(mt_dbl, nqpts*nbf_ao,
179     &             'chi_ao', lchi_ao, ichi_ao))
180     &             call errquit('zora_getv: chi_ao',0, MA_ERR)
181               call qdist(dbl_mb(i_rchi_atom), dbl_mb(i_rq),
182     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters)
183               call xc_eval_basis(ao_bas_han, 1, dbl_mb(ichi_ao),
184     &              dbl_mb(idelchi_ao), 0d0, 0d0, dbl_mb(i_rq),
185     &              dbl_mb(iqxyz), dbl_mb(ixyz), nqpts, ncenters,
186     &              int_mb(iniz), log_mb(idocset),
187     &              int_mb(ibas_cent_info), int_mb(ibas_cset_info))
188               if (zora_calc_type.eq.3) then
189c                write(*,*) 'enter calc_zora_EFGZ4_so() ...'
190                call calc_zora_EFGZ4_so(ao_bas_han,geom,ipol,g_dens,
191     &                                  dbl_mb(idelchi_ao),
192     &                                  dbl_mb(iqxyz),dbl_mb(iqwght),
193     &                                  nbf_ao,nqpts,ncenters,
194     &                                  zora_Qpq,xyz_EFGcoords,
195     &                                  use_modelpotential,
196     &                                  gexpo, gcoef,
197     &                                  dbl_mb(iscal(1)), ! out
198     &                                  dbl_mb(iscal(2)), ! out
199     &                                  dbl_mb(iscal(3)), ! out
200     &                                  dbl_mb(iscal(4))) ! out
201               else if (zora_calc_type.eq.4) then
202c                write(*,*) 'enter calc_EFG() ...'
203                call calc_EFG(ao_bas_han,geom,ipol,g_dens,
204     &                        dbl_mb(ichi_ao),
205     &                        dbl_mb(iqxyz),dbl_mb(iqwght),
206     &                        nbf_ao,nqpts,ncenters,
207     &                        zora_Qpq,xyz_EFGcoords,
208     &                        dbl_mb(iscal(1))) ! out
209               endif
210c              == delete memory ==
211               if(.not.MA_chop_stack(l_rchi_atom))
212     &            call errquit("zora_getv: pop rchi_atom",100,MA_ERR)
213            endif ! nqpts
214         enddo ! ncube
215      end do ! iqsh
216c
217c     == delete memory ==
218      if(.not.MA_chop_stack(lbas_cent_info))
219     &     call errquit("zora_getv_sf: pop bas_cent_info", 100, MA_ERR)
220c
221c     == tally up over all the nodes ==
222      do k=1,4
223         call ga_dgop(msg_excrho, dbl_mb(iscal(k)), nbf_ao*nbf_ao, '+')
224      enddo
225c
226c     == pack into a ga ==
227c
228c     == scalar contribution ==
229      do i=1,ipol
230       call ga_zero(g_efgz4_sf(i))
231       call ga_put(g_efgz4_sf(i),1,nbf_ao,1,nbf_ao,dbl_mb(iscal(1))
232     &            ,nbf_ao)
233       call ga_symmetrize(g_efgz4_sf(i))
234      enddo
235c     == spin-orbit contributions ==
236      ind=4 ! = 4,3,2 = z,y,x in g_zora_so(i)
237      do i=1,3 ! = z,y,x
238       call ga_zero(g_efgz4_so(i))
239       call ga_put(g_efgz4_so(i),
240     &             1,nbf_ao,1,nbf_ao,dbl_mb(iscal(ind)),nbf_ao)
241       call ga_antisymmetrize(g_efgz4_so(i))
242       ind=ind-1
243      enddo
244      if (.not.MA_chop_stack(lscal(1)))
245     &    call errquit('zora_getv_so: scali chop_stack',0, MA_ERR)
246      call ga_sync()
247      return
248      end
249