1czora...Build the zora contribution to the fock matrix
2czora...This is routine is analogous to the xc_tabcd routine that
3czora...contructs the exchange-correlation contribution
4czora...TODO: Clean up and more comments
5c
6      subroutine zora_tabcd(what,l3d,
7     ,     tol_rho, Tmat, TTmat, Amat, Bmat, Cmat, Dmat,
8     N     Emat,Fmat,qxyz,xyz,
9     &     chi, delchi, heschi,
10     N     curatoms,ncuratoms,nmat,
11     I     ipol, nq, nbf, max_at_bf, max_at_bf2,
12     G     GRAD, basis, natoms, iniz, ifin,
13     &     g_zora, ibf, rchi_atom, rdelchi_atom,
14     &     rdens_atom, cetobfr,kske,Mmat,scr)
15c
16      implicit none
17#include "mafdecls.fh"
18#include "dftpara.fh"
19#include "dft2drv.fh"
20#include "dist.fh"
21#include "dft_fdist.fh"
22c
23czora...start
24#include "zora.fh"
25czora...end
26c
27c
28      Logical GRAD
29      integer what ! [in]
30cwhat=0
31cwhat=1 CPKS_LHS
32cwhat=2 CPKS_RHS
33cwhat=2 NMR_RHS
34      integer basis
35      integer max_at_bf ! [input]
36      integer max_at_bf2 ! [input]
37      integer nmat ! Number of XC matrices (alpha + beta sets) to make [input]
38c                  ! e.g. number of perturbations for CPSCF
39      integer imat ! XC matrix loop index
40      integer ipol  ! [input]
41      integer nq    ! [input]
42      integer nbf    ! [input]
43      integer natoms ! [input]
44      integer ncuratoms ! Number of current "active" atoms [input]
45      integer curatoms(*) ! Mapping array for current atoms [input]
46      integer jatcur,  nu, nu1,  indT
47      double precision tol_rho
48      double precision Tmat(*), TTmat(*)
49      double precision rchi_atom(natoms)
50      double precision rdelchi_atom(natoms)
51      double precision rdens_atom(natoms,natoms,ipol)
52      integer cetobfr(2,natoms)
53c
54czora...start
55      integer g_zora(*)
56czora...end
57c
58c     Sampling Matrices for the XC Potential & Energy
59c
60      double precision Amat(nq,ipol,*), Cmat(nq,3,ipol,*)
61      double precision Mmat(nq,ipol)
62c     nmr
63      double precision Emat(nq,max_at_bf)
64      double precision Fmat(nq,3,max_at_bf)
65      double precision qxyz(3,*)                     ! [input] grid point coordinates
66      double precision xyz(3,*)                  ! [input] nuclear coordinates
67      logical kske
68c#elif defined(TABCD_CPKS_LHS)
69c
70c     Note: Meaning of dimensioning of Amat and Cmat changes for
71c           second derivatives, simulating "overloading" of
72c           Amat and Cmat
73c
74c     Sampling Matrices for the XC part of integrand when making
75c     multiple matrices, e.g. XC part of perturbations for CPSCF
76c
77c      double precision Amat(nq,ipol,nmat), Cmat(nq,3,ipol,nmat)
78c#elif defined(TABCD_CPKS_RHS)
79c
80c     For explicit nuclear derivatives of XC matrix, the same functional
81c     derivative values are combined with different basis fn derivatives
82c
83c      double precision Amat(nq,ipol), Cmat(nq,3,ipol)
84
85c
86c     Sampling Matrices for [Products of] Basis Functions & Gradients
87c
88      double precision Bmat(nq,max_at_bf)
89      double precision Dmat(nq,3,max_at_bf)
90      integer iniz(natoms), ifin(natoms)
91c
92c     Basis Functions & Gradients
93c
94      double precision chi(nq,nbf), delchi(nq,3,nbf)
95      double precision heschi(nq,6,*)
96      integer ibf(nbf)
97      double precision A_MAX, C_MAX, AC_MAX, FUNC_MAXI,
98     &                 B_MAX, D_MAX, BD_MAX, FUNC_MAXJ
99      integer iat, inizia, ifinia, nbfia, nnia, ifirst, ilast
100      integer jat, inizja, ifinja, nbfja, nnja, jfirst, jlast
101      integer ii, mu, mu1
102      integer n,lastjat
103      double precision chi1
104      double precision scr(nq)
105      double precision dabsmax
106      double precision tolrho15
107      external dabsmax
108      logical l3d
109      integer jrsh,jrsh2,n3d,idir,jdir
110      integer g_update(2)
111cnmr
112      double precision Rij(3) ! vector R(jat) - R(iat)
113      integer inia, iq, ix, ix1, ix2
114c
115c     Indexing array for basis function hessian columns as if
116c     it were a 3x3 matrix
117c
118      integer indh(3,3)
119      logical w01,w02,w013,dofull
120      double precision ddot
121      double precision raa,rbb
122      integer nbhandl1,nbhandl2
123      logical nbfirst1,nbfirst2,doitt
124      integer sizeblk
125#include "nwc_const.fh"
126      integer nonzero(nw_max_atom),natleft,
127     A     iat0,jat0
128      external ddot
129      data indh / 1, 2, 3,
130     &            2, 4, 5,
131     &            3, 5, 6 /
132
133c         0: l3d=.f.    & n3d=1
134ccc	rhs: l3d=.true. & n3d=3
135ccc	lhs: l3d=.true. & n3d=1
136c
137      call starttimer(monitor_tabcd)
138      call starttimer(monitor_screen0)
139      natleft=0
140      do  iat = 1, natoms
141         if(iniz(iat).ne.0) then
142            natleft=natleft+1
143            nonzero(natleft)=iat
144         endif
145      enddo
146      tolrho15=tol_rho**1.25d0
147      if(what.eq.0) then
148        n3d=1
149      elseif(what.eq.1) then
150         n3d=1
151      elseif(what.eq.2) then
152         n3d=3
153      elseif(what.eq.3) then
154         n3d=3
155      else
156         call errquit(' wrong what value for xctabcd ',0,0)
157      endif
158      w01=what.eq.0.or.what.eq.1
159      w013=w01.or.what.eq.3
160      w02=what.eq.0.or.what.eq.2
161      dofull=what.ne.0.or.l3d
162      nbfirst1=.true.
163      nbfirst2=.true.
164      call endtimer(monitor_screen0)
165c
166c
167c     Beginning of loop over multiple XC matrices
168c
169      do 500 imat = 1,nmat
170c
171c     Compute the matrix product for the XC potential and energy:
172c
173c              T = transpose(A*B) + transpose(C*D)
174c
175
176         call starttimer(monitor_screen1)
177
178         A_MAX = dabsmax(nq*ipol,Amat(1,1,imat))
179         if (GRAD) then
180            C_MAX = dabsmax(nq*3*ipol,Cmat(1,1,1,imat))
181         else
182            C_MAX = 0d0
183         endif
184         AC_MAX = max(A_MAX,C_MAX)
185
186         call endtimer(monitor_screen1)
187
188c
189c     repl stuff
190c
191      if(xcreplicated.and.dorepxc) then
192         g_update(1)=k_repxc(1)
193         g_update(2)=k_repxc(2)
194      else
195         g_update(1)=g_zora(1)
196         g_update(2)=g_zora(2)
197      endif
198      do 430 iat0=1,natleft
199         call starttimer(monitor_screen2)
200         iat=nonzero(iat0)
201         inizia = iniz(iat)
202         FUNC_MAXI = rchi_atom(iat)
203         if(GRAD) FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat))
204         doitt=(AC_MAX*FUNC_MAXI).ge.tol_rho
205         call endtimer(monitor_screen2)
206         if(what.eq.2.or.doitt) then
207         ifinia = ifin(iat)
208         ifirst = cetobfr(1,iat)
209         ilast = cetobfr(2,iat)
210         nnia = ifinia - inizia + 1
211         nbfia = ilast - ifirst + 1
212         do ii = 1, ipol
213            do mu = 1, nnia
214
215               call starttimer(monitor_mult1)
216
217               mu1 = mu+inizia-1
218               if (GRAD) then
219                 do n = 1, nq
220                  chi1 = chi(n,mu1)
221
222                  Bmat(n,mu) = 0.d0
223
224                  if (so_term .eq. 1) then        ! z-component
225                    Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat)
226                    Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat)
227                    Dmat(n,3,mu) = 0.d0
228                  else if (so_term .eq. 2) then   ! y-component
229                    Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat)
230                    Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat)
231                    Dmat(n,2,mu) = 0.d0
232                  else if (so_term .eq. 3) then   ! x-component
233                    Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat)
234                    Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat)
235                    Dmat(n,1,mu) = 0.d0
236                  else if (so_term .eq. 0) then   ! spin-free
237                    Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat)
238                    Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat)
239                    Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat)
240                  end if
241
242                 enddo   ! nq: over grid points
243               else
244                  do n = 1, nq
245                     Bmat(n,mu) = 0.d0
246                  enddo   ! nq: over grid points
247               endif
248
249            enddo  ! mu: over the basis set
250
251c     Monitoring
252
253            call endtimer(monitor_mult1)
254c
255            call starttimer(monitor_screen3)
256            B_MAX = dabsmax(nnia*nq,Bmat)
257            if (GRAD) then
258               D_MAX = dabsmax(nnia*nq*3,Dmat)
259            else
260               D_MAX = 0d0
261            endif
262            BD_MAX = max(B_MAX,D_MAX)
263c
264            lastjat=iat0
265            if(what.eq.2) lastjat=natleft
266            if(what.eq.3) lastjat=iat0-1
267            call endtimer(monitor_screen3)
268            do 168 jat0=1,lastjat
269               jat=nonzero(jat0)
270               if(what.eq.2) then
271c
272c     To fit better into existing structure, loop over full square
273c     of atom pairs and only compute nuclear derivative contribution
274c     from jat.  Also, this way we only need check jatcur once and
275c     for all, and don't have to check iatcur at all.
276c
277                  jatcur = curatoms(jat)
278                  if (jatcur.eq.0) goto 168
279               endif
280               call starttimer(monitor_screen4)
281               inizja = iniz(jat)
282               FUNC_MAXJ = rchi_atom(jat)
283               if(grad) FUNC_MAXJ = max(rchi_atom(jat),FUNC_MAXJ)
284               doitt=(BD_MAX*FUNC_MAXJ).ge.tol_rho
285               call endtimer(monitor_screen4)
286               if(what.eq.2.or.doitt) then
287
288c     Monitoring
289
290            call starttimer(monitor_mult2)
291
292               if(what.eq.3) then
293                     Rij(1) = 0.5d0*(xyz(1,jat)-xyz(1,iat))
294                     Rij(2) = 0.5d0*(xyz(2,jat)-xyz(2,iat))
295                     Rij(3) = 0.5d0*(xyz(3,jat)-xyz(3,iat))
296               endif
297               ifinja = ifin(jat)
298               jfirst = cetobfr(1,jat)
299               jlast = cetobfr(2,jat)
300               nbfja = jlast - jfirst + 1
301               nnja = ifinja - inizja + 1
302               sizeblk=n3d*nbfia*nbfja
303               if(what.eq.2.or.what.eq.3)
304     &            call dcopy(sizeblk, 0d0,0, TTmat,1)
305c
306c              Loop over x, y, z directions for derivative XC mats
307c
308               do jdir = 1,n3d
309c
310                  if(what.eq.3) then
311                     ix1 = mod(jdir,3)+1
312                     ix2 = mod(jdir+1,3)+1
313                     raa=rij(ix1)
314                     rbb=rij(ix2)
315                     do iq = 1, nq
316                        scr(iq) = raa*qxyz(ix2,iq) - rbb*qxyz(ix1,iq)
317                     enddo
318                     do inia = 1, nnia
319                        do iq = 1, nq
320                           Emat(iq,inia) = scr(iq)*Bmat(iq,inia)
321                        enddo
322                        if (GRAD) then
323                           do iq = 1, nq
324                               Emat(iq,inia) = 0.d0
325                           enddo
326                        endif
327                     enddo
328                     if (GRAD) then
329                        do inia = 1, nnia
330                              do iq = 1, nq
331                               Fmat(iq,1,inia) = 0.d0
332                               Fmat(iq,2,inia) = 0.d0
333                               Fmat(iq,3,inia) = 0.d0
334                           enddo
335                        enddo
336                     endif
337                  endif
338c
339                  if (GRAD) then
340                   indT = 0
341                   do nu = 1, nnja
342                     nu1 = nu+inizja-1
343                     do mu = 1, nnia
344                      indT = indT + 1
345
346czora...               build Tmat
347
348                      if (so_term .eq. 1) then       ! z-component
349                       Tmat(indT) =
350     &                   ddot(nq,Dmat(1,1,mu),1,delchi(1,2,nu1),1)
351     &                 - ddot(nq,Dmat(1,2,mu),1,delchi(1,1,nu1),1)
352                      else if (so_term .eq. 2) then  ! y-component
353                       Tmat(indT) =
354     &                   ddot(nq,Dmat(1,3,mu),1,delchi(1,1,nu1),1)
355     &                 - ddot(nq,Dmat(1,1,mu),1,delchi(1,3,nu1),1)
356                      else if (so_term .eq. 3) then  ! x-component
357                       Tmat(indT) =
358     &                   ddot(nq,Dmat(1,2,mu),1,delchi(1,3,nu1),1)
359     &                 - ddot(nq,Dmat(1,3,mu),1,delchi(1,2,nu1),1)
360                      else if (so_term .eq. 0) then  ! spin-free
361                       Tmat(indT) =
362     &                   ddot(nq,Dmat(1,1,mu),1,delchi(1,1,nu1),1)
363     &                 + ddot(nq,Dmat(1,2,mu),1,delchi(1,2,nu1),1)
364     &                 + ddot(nq,Dmat(1,3,mu),1,delchi(1,3,nu1),1)
365                      end if
366
367                     enddo ! mu: over the basis set
368                   enddo ! nu: over the basis set
369
370                  endif
371
372                  if(n3d.eq.1) then
373                     call dfill(max_at_bf2, 0.d0, TTmat, 1)
374                     if(dofull) then
375                        call scat_mat(TTmat, Tmat, nbfia, nbfja, nnia,
376     &                     nnja,ifirst,jfirst,ibf(inizia),ibf(inizja))
377                     else
378                      if (so_term.eq.0) then
379                        call scat_matup(TTmat, Tmat, nbfia, nbfja, nnia,
380     &                     nnja,ifirst,jfirst,ibf(inizia),ibf(inizja))
381                      else
382                        call scat_mat(TTmat, Tmat, nbfia, nbfja, nnia,
383     &                     nnja,ifirst,jfirst,ibf(inizia),ibf(inizja))
384                      end if
385                     endif
386                  else
387                     call scat_mat3(n3d,jdir,
388     &                    TTmat, Tmat, nbfia, nbfja, nnia,
389     &                    nnja,ifirst,jfirst,ibf(inizia),ibf(inizja))
390                  endif
391
392               enddo            ! jdir (loop over x, y, z directions for nmr)
393
394c     Monitoring
395
396            call endtimer(monitor_mult2)
397
398                  doitt=.true.
399            call starttimer(monitor_screen5)
400               if(what.eq.0) then
401                  doitt=dabsmax(sizeblk,ttmat).ge.tolrho15
402                  jrsh=ii
403               elseif(what.eq.1) then
404                  doitt=.true.
405                  jrsh=imat+(ii-1)*nmat
406               elseif(what.eq.3) then
407                  jrsh=(ii-1)*n3d+1
408               else
409                  doitt=dabsmax(sizeblk,ttmat).ge.tol_rho
410                  jrsh=1+(jat-1)*3+(ii-1)*3*natoms
411               endif
412               call endtimer(monitor_screen5)
413               if(doitt) then
414
415                  jrsh2=jrsh+n3d-1
416
417c     Monitoring
418
419               call updist(monitor_size_ga_acc1, sizeblk)
420               call starttimer( monitor_comm_ga_acc1)
421
422                  if(l3d) then
423                     call dft_3dacc(g_zora, ttmat,
424     &                    jrsh,jrsh2,
425     &                 ifirst, ilast, jfirst, jlast, nbfia)
426                  else
427                     if(dftnbacc) then
428
429                        if(.not.nbfirst1) then
430                           call starttimer( monitor_wait1)
431                           call ga_nbwait(nbhandl1)
432                           call endtimer( monitor_wait1)
433                        endif
434                        nbfirst1=.false.
435                        call upd_atombl_nb(g_update(ii),
436     &                       basis,iat,jat,ttmat,nbhandl1)
437                     else
438                        call upd_atom_block(g_update(ii),
439     &                       basis,iat,jat,ttmat)
440
441
442                     endif
443                  endif
444
445
446
447c     Monitoring
448
449               call endtimer( monitor_comm_ga_acc1)
450
451
452               if(what.ne.0.or.l3d) then
453c
454c                 check to see if can skip and use ga_symmetrize
455c
456                  if ((w013.and.iat.ne.jat).or.what.eq.2) then
457c     For CPKS RHS, we update with transpose even for iat = jat,
458c     since that is necessary to get both contributions
459c     mu * del(nu) and del(mu) * nu
460
461
462                     call starttimer(monitor_comp_transp)
463
464                     if(n3d.eq.1) then
465                        call transp_mat(TTmat, Tmat,
466     ,                       nbfia, nbfja)
467                     else
468                        if(what.eq.3) then
469                           call dscal(n3d*nbfia*nbfja,-1.0d0,TTmat,1)
470                        endif
471                        call transp_mat3(n3d,TTmat, Tmat,
472     ,                       nbfia, nbfja)
473                     endif
474
475                     call endtimer(monitor_comp_transp)
476
477
478c     Monitoring
479
480                     call starttimer(monitor_comm_ga_acc2)
481
482
483                     if(l3d) then
484
485                        call dft_3dacc(g_zora, tmat,
486     &                       jrsh,jrsh2,
487     %                       jfirst, jlast, ifirst, ilast, nbfja)
488                     else
489                        if(dftnbacc) then
490
491                           if(.not.nbfirst2) then
492                              call ga_nbwait(nbhandl2)
493                           endif
494                           nbfirst2=.false.
495                           call upd_atombl_nb(g_update(ii),
496     .                          basis,jat,iat,tmat,nbhandl2)
497                        else
498
499                           call upd_atom_block(g_update(ii),basis,
500     J                          jat,iat,tmat)
501                        endif
502                     endif
503c     Monitoring
504
505                     call endtimer(monitor_comm_ga_acc2)
506
507
508                  endif
509               endif
510
511            endif
512           endif
513
514 168       continue
515        enddo ! ipol loop
516      endif
517 430  continue
518 500  continue
519
520      return
521      end
522c $Id$
523