1#define STAG 1
2* $Id$
3c     Modified for construction of multiple sets of perturbed densities
4c     and explicit derivative densities wrt nuclei
5c
6c     BGJ - 8/98
7c
8C> \ingroup nwdft_xc
9C> @{
10C>
11C> \file xc_rhogen.F
12C> Evaluate density and derivatives
13C>
14C> @}
15C>
16C> \ingroup nwdft_xc
17C> @{
18C>
19C> \brief Evaluate the density and their derivatives on a grid
20C>
21C>
22      Subroutine xc_rhogen(what_in,
23     T     tol_rho, basis, g_dens, max_at_bf,
24     N     natoms, curatoms, ncuratoms, npert,
25     I     ipol, nq, nbf, mbf, GRAD, ipol2,
26     &     F, Pmat, ff, ffd,
27     C     chi, delchi, heschi,
28     I     ibf, iniz, ifin,
29     &     rho, delrho, lap, rchi_atom,
30     &     rdelchi_atom, rdens_atom, cetobfr, wmax,
31     &     ttau, kske, dolap)
32      implicit none
33#include "global.fh"
34#include "errquit.fh"
35#include "mafdecls.fh"
36#include "dftpara.fh"
37#include "dist.fh"
38#include "dft_fdist.fh"
39
40      Logical GRAD !< [Input] .true. when using gradient corrected
41                   !< functional
42      Logical kske !< [Input] .true. when using kinetic energy density
43                   !< functional
44      integer what_in !< [Input] What to calculate
45                   !< - what=0: calculate density
46                   !< - what=1: calculate perturbed density (derivative
47                   !<   with respect to electron position)
48                   !< - what=2: calculate derivative with respect to
49                   !<   nuclear coordinates
50       logical dolap       !< [Input] .true. if the Laplacian of the
51                           !< density is required
52c what=-1 dplot aka density contours
53c what=0 dens
54c what=1 pert
55c what=2 nucder
56      integer basis        !< [Input] basis set handle
57      integer ipol         !< [Input] no. of spin channels
58      integer ipol2        !< [Input] no. of spin channels in density
59                           !< - 1 if closed shell
60                           !< - 3 if open shell
61      integer nbf          !< [Input] no. of basis functions
62      integer mbf          !< [Input] "restricted" no. of basis functions
63      integer max_at_bf    !< [Input] max no. bf per atom
64      integer nq           !< [Input] no. of quadrature points
65      integer natoms       !< [Input] no. of atoms
66      double precision wmax !< [Input] max weight
67      integer npert        !< [Input] number of perturbed densities
68      integer curatoms(*) !< [Input] indexing array for "active" atoms
69      integer ncuratoms        !< [Input] number of currently active atoms
70      integer g_dens(*) !< [Input] GA handle for DM
71      integer ibf(mbf)     !< [Input] mapping of nbf_ao -> mbf
72      integer iniz(natoms) !< [Input] mapping of nbf_ao -> mbf
73      integer ifin(natoms) !< [Input] mapping of nbf_ao -> mbf
74      double precision tol_rho !< [Input] accuracy for rho evaluation
75      double precision chi(nq,mbf)     !< [Input] function values
76      double precision delchi(nq,3,mbf)!< [Input] function gradients
77      double precision heschi(nq,6,mbf)!< [Input] function hessians
78c
79      double precision ttau(nq,ipol,*) !< [Output] Total Kohn-Sham K.E.density
80c
81       double precision lap(nq,ipol2,*)
82c
83      double precision delrho(nq,3,ipol,*) !< [Output] Derivative of density
84      double precision Pmat(*) !< [Scratch] scratch vector
85      double precision F(max_at_bf*max_at_bf) !< [Scratch] scratch vector
86      double precision ff(nq,*)   !< [Scratch] scratch array
87      double precision ffd(nq,3,*) !< [Scratch] scratch array
88      double precision rho(nq,ipol2,*) !< [Output] The density
89      double precision rchi_atom(natoms) !< [Input] Screening parameters
90      double precision rdelchi_atom(natoms) !< [Input] Screening parameters
91      double precision rdens_atom(natoms,natoms,ipol) !< [Input] Screening parameters
92      integer cetobfr(2,natoms) !< [Input] Centers to basis functions
93c
94c     local declarations
95c
96      integer i0, ii, mu, n, npol
97      integer ipert        ! perturbation loop index
98      integer iat, inizia, ifirst, ilast, nbfia, nnia, iat_in
99      integer ifinia, ifinja
100      integer jat, inizja, jfirst, jlast, nbfja, nnja
101      integer iatcur, jatcur
102      double precision FUNC_MAX, DELFUNC_MAX, FUNC_MAXI, FUNC_MAXJ
103      double precision P_MAXJ, P_MAXJ_A, P_MAXJ_B, P_MAXIJ
104      double precision dabsmax
105      integer g_keepd(2)
106      integer nbhandl
107      integer jj
108      logical doffd,doitt
109      external dabsmax
110      external xc_rhoscreen
111      integer xc_rhoscreen
112      integer nonzero,nonz0
113      integer i_nz,l_nz
114      integer sizeblk, gindx
115      integer what
116      logical zapnegatives
117#ifdef DEBUG
118      integer ga_nodeid
119      external ga_nodeid
120#endif
121      g_keepd(1)=0
122      g_keepd(2)=0
123      call starttimer(monitor_xcrho)
124      if(what_in.ge.0) then
125         what=what_in
126         zapnegatives=.true.
127      else
128         what=-what_in
129         zapnegatives=.false.
130      endif
131c
132c     Evaluate the charge density and its gradient at each of the
133c     sampling points
134c
135      doffd=(what.eq.2.and.grad).or.(what.eq.0.and.kske).or.
136     &      (what.eq.1.and.kske).or.(what.eq.0.and.dolap)
137
138      npol = (ipol*(ipol+1))/2
139c     to keep compilers quiet
140      iatcur=1
141      jatcur=1
142c
143      if(what.eq.0) then
144         call dcopy(nq*npol,0d0,0,rho,1)
145         if (grad) call dcopy(3*nq*ipol,0d0,0,delrho,1)
146         if (kske) call dcopy(nq*ipol,0d0,0,ttau,1)  ! total
147         if (dolap) call dcopy(nq*ipol2,0d0,0,lap,1)  ! total
148c
149c     repl stuff
150c
151         g_keepd(1)=g_dens(1)
152         if(ipol.eq.2) g_keepd(2)=g_dens(2)
153         if(xcreplicated.and.dorepdm) then
154            g_dens(1)=g_repdm(1)
155            if(ipol.eq.2) g_dens(2)=g_repdm(2)
156         endif
157      elseif(what.eq.1) then
158         call dcopy(nq*ipol*npert,0d0,0,rho,1)
159         if (grad) call dcopy(3*nq*ipol*npert,0d0,0,delrho,1)
160         if (kske) call dcopy(nq*ipol*npert,0d0,0,ttau,1)   ! total
161      elseif(what.eq.2) then
162         call dcopy(nq*ipol*3*ncuratoms,0d0,0,rho,1)
163         if (grad)call dcopy(3*nq*ipol*3*ncuratoms,0d0,0,delrho,1)
164      else
165         call errquit('wrong what value',0,0)
166      endif
167c
168c     Screening is accomplished by:  p(r) <= |Xi(r)|*|Xj(r)|*|Dij|
169c     Xi(r) is screened on desired accuracy/max(|Xj(r)|)*max(|Dij|)
170c     Dij is screened on desired accuracy/max(|Xi(r)|)*max(|Xj(r)|)
171c
172      i0=ipol-1
173c
174      FUNC_MAX = dabsmax(natoms,rchi_atom)
175      DELFUNC_MAX=0d0
176      if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom)
177c
178      nonzero=0
179      if(dftnbget) then
180         if (.not.ma_push_get
181     &        (mt_int,(ipol*npert*natoms*(natoms+1))/2,
182     N        'nzmap map',l_nz,i_nz))
183     &        call errquit('xcrho:push_get failed', 13, MA_ERR)
184         nonzero=xc_rhoscreen(grad,ipol,natoms,npert,
185     I        iniz,
186     W        tol_rho,wmax,
187     O        int_mb(i_nz),
188     R        rchi_atom,rdelchi_atom,rdens_atom)
189
190         if(nonzero.eq.0) goto 1688
191c
192c     prefetch first DM block
193c
194         nonz0=1
195         call xc_getdmblock(int_mb(i_nz),nonz0,natoms,cetobfr,
196     G        g_dens(1),
197     A        Pmat,nbhandl)
198      endif
199
200#ifdef STAG
201      do 230 iat_in = 1+ga_nodeid(), natoms+ga_nodeid()
202         iat=mod(iat_in,natoms)
203         if(iat.eq.0) iat=natoms
204#else
205      do 230 iat = 1, natoms
206#endif
207         inizia = iniz(iat)
208         if (inizia.eq.0)goto 230
209         if(what.eq.2) then
210            iatcur = curatoms(iat)
211         endif
212         ifinia = ifin(iat)
213         ifirst = cetobfr(1,iat)
214         ilast = cetobfr(2,iat)
215         nbfia = ilast-ifirst+1
216         nnia = ifinia-inizia+1
217c
218c        screening parameters
219c
220         FUNC_MAXI = rchi_atom(iat)
221         if(grad)
222     .     FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat))
223         FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX)
224         if(what.lt.2) then
225            if (ipol.gt.1)then
226              P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1))
227              P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2))
228              P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B)
229            else
230              P_MAXJ=0d0
231              do jat=1,iat
232                 if(iniz(jat).ne.0)
233     .                P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1))
234              enddo
235            endif
236            if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.lt.tol_rho) goto 225
237         endif
238         do 220 jat = 1, iat
239            inizja = iniz(jat)
240            if (inizja.eq.0)goto 220
241            if(what.eq.2) then
242              jatcur = curatoms(jat)
243              if (jatcur.eq.0.and.iatcur.eq.0) goto 220
244            endif
245            call starttimer(monitor_rscreen0)
246            ifinja = ifin(jat)
247            jfirst = cetobfr(1,jat)
248            jlast = cetobfr(2,jat)
249            nbfja = jlast-jfirst+1
250            nnja = ifinja-inizja+1
251c
252c           screening parameters
253c
254            FUNC_MAXJ=rchi_atom(jat)
255            if(grad)
256     .        FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat))
257            doitt=.true.
258            if(what.lt.2) then
259              P_MAXIJ = rdens_atom(iat,jat,1)
260              if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ,
261     &             rdens_atom(iat,jat,2))
262              doitt=(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho)
263            endif
264            call endtimer(monitor_rscreen0)
265            if (doitt) then
266              do ii = 1, ipol
267c
268c               screening parameters
269c
270                if((what.gt.1).or.(rdens_atom(iat,jat,ii)*
271     R               wmax*FUNC_MAXI*FUNC_MAXJ.ge.tol_rho)) then
272c
273c                  Loop over perturbations
274c
275                   do 215 ipert = 1,npert
276                     sizeblk=nbfia*nbfja
277                     call updist(monitor_size_ga_get, sizeblk)
278c
279                     if(dftnbget) then
280                       call starttimer(monitor_wait3)
281                       call ga_nbwait(nbhandl)
282                       call endtimer(monitor_wait3)
283                       call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja,
284     I                      ifirst, jfirst, ibf(inizia), ibf(inizja))
285                       nonz0=nonz0+1
286                       if((npert*ipol).gt.1) then
287                         gindx=ipert+(ii-1)*npert+1
288
289                         if(gindx.gt.npert*ipol)
290     +                      gindx=mod(gindx,npert*ipol)
291                       else
292                         gindx=1
293                       endif
294                       if(nonz0.le.nonzero)
295     C                   call xc_getdmblock(int_mb(i_nz),nonz0,natoms,
296     C                   cetobfr,g_dens(gindx),
297     A                      Pmat,nbhandl)
298                       if (wmax*FUNC_MAXI*FUNC_MAXJ*
299     .                     dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215
300                     else
301                       call starttimer(monitor_gaget)
302                       if(truerepdm) then
303                         call xc_dmget(
304     +                      dbl_mb(k_repdm(ipert+(ii-1)*npert)),nbf_ld,
305     %                      ifirst, ilast, jfirst, jlast, Pmat,nbfia)
306                       else
307                         call ga_get(g_dens(ipert+(ii-1)*npert),
308     %                        ifirst, ilast, jfirst, jlast, Pmat,nbfia)
309                       endif
310                       call endtimer(monitor_gaget)
311                       if (wmax*FUNC_MAXI*FUNC_MAXJ*
312     .                     dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215
313                       call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja,
314     +                      ifirst, jfirst, ibf(inizia), ibf(inizja))
315                     endif
316c
317                     call starttimer(monitor_rhocomp)
318                     if(iat.ne.jat) call dscal(nnia*nnja,2d0,F,1)
319c
320c                    Compute Xiat(r)*Xjat(r)*Diat,jat
321c
322                     call dgemm('n','n',nq,nnja,nnia,1d0,
323     A                    chi(1,inizia),nq,F,nnia,0d0,ff,nq)
324                     if(what.lt.2) then
325                        jj=i0+ii
326                        if(what.eq.1) jj=ii
327                        do mu=inizja,ifinja
328                           do n=1,nq
329                              rho(n,jj,ipert) = rho(n,jj,ipert) +
330     P                             chi(n,mu)*ff(n,mu-inizja+1)
331                           enddo
332                        enddo
333                     endif
334                     if(doffd) then
335                       call dgemm('n','n',nq*3,nnja,nnia,1d0,
336     A                      delchi(1,1,inizia),nq*3,F,nnia,0d0,ffd,nq*3)
337c
338c                      build tau for meta GGA
339c
340                       if(kske) then
341                         do mu=inizja,ifinja
342                           do n=1,nq
343                             ttau(n,ii,ipert) = ttau(n,ii,ipert)+0.5d0*(
344     &                         delchi(n,1,mu)*ffd(n,1,mu-inizja+1)+
345     &                         delchi(n,2,mu)*ffd(n,2,mu-inizja+1)+
346     &                         delchi(n,3,mu)*ffd(n,3,mu-inizja+1))
347                           enddo
348                         enddo
349                       endif
350
351                       if (dolap) then
352                         do mu=inizja,ifinja
353                           do n=1,nq
354                              ! total
355                              lap(n,ii,ipert) = lap(n,ii,ipert) + 2d0*
356     &                            (delchi(n,1,mu)*ffd(n,1,mu-inizja+1) +
357     &                             delchi(n,2,mu)*ffd(n,2,mu-inizja+1) +
358     &                             delchi(n,3,mu)*ffd(n,3,mu-inizja+1))
359cold     &                             delchi(n,3,mu)*ffd(n,3,mu-inizja+1) +
360cold     &                             heschi(n,1,mu)*ff(n,mu-inizja+1) +
361cold     &                             heschi(n,4,mu)*ff(n,mu-inizja+1) +
362cold     &                             heschi(n,6,mu)*ff(n,mu-inizja+1))
363                           enddo
364                         enddo
365                       endif
366                     endif
367c
368                     if((what.eq.2.and.jatcur.ne.0).or.
369c
370c                      We need the "grad" code at zero order in the
371c                      nuclear deriv case, but we can skip this part if
372c                      iat is not active
373c
374     O                  (what.ne.2.and.grad)) then
375                       if(what.eq.0.or.what.eq.1) then
376                         call xc_dchiff(nq,inizja,ifinja,
377     P                        delrho(1,1,ii,ipert),delchi,
378     F                        ff)
379cc DMR/Begin
380                         if (dolap) then
381                            do mu=inizja,ifinja
382                               do n=1,nq
383                                  lap(n,ii,ipert) = lap(n,ii,ipert) +
384     &                            ff(n,mu-inizja+1)*(heschi(n,1,mu) +
385     &                            heschi(n,4,mu) + heschi(n,6,mu))
386                               enddo
387                            enddo
388                         endif
389cc DMR/End
390                       elseif(what.eq.2) then
391                         call xc_dchiffp(nq,ipol2,inizja,ifinja,
392     P                        rho,delchi,ff,
393     M                        ii,jat)
394                         if (grad) then
395c
396c                          Compute nuclear gradient of delrho
397c
398                           call xc_drhonuc(nq,ipol,inizja,ifinja,
399     D                          delrho,heschi,delchi,ff,ffd,
400     I                          ii,jat)
401                         endif
402                       endif
403                     endif
404c
405                     if((what.eq.2.and.iatcur.ne.0).or.
406     O                  (what.ne.2.and.grad)) then
407c
408c                      Compute delXiat(r)*Xjat(r)*Diat,jat
409c
410                       call dgemm('n','t',nq,nnia,nnja,1d0,
411     A                      chi(1,inizja),nq,F,nnia,0d0,ff,nq)
412                       if(what.lt.2) then
413                         call xc_dchiff(nq,inizia,ifinia,
414     P                        delrho(1,1,ii,ipert),delchi,
415     F                        ff)
416cc DMR/Begin
417                         if (dolap) then
418                            do mu=inizia,ifinia
419                               do n=1,nq
420                                  lap(n,ii,ipert) = lap(n,ii,ipert) +
421     &                              ff(n,mu-inizia+1)*(heschi(n,1,mu) +
422     &                              heschi(n,4,mu) + heschi(n,6,mu))
423                               enddo
424                            enddo
425                         endif
426cc DMR/End
427                       elseif(what.eq.2) then
428                         call xc_dchiffp(nq,ipol2,inizia,ifinia,
429     P                        rho,delchi,ff,
430     M                        ii,iat)
431                         if(grad) then
432                           call dgemm('n','t',nq*3,nnia,nnja,1d0,
433     A                          delchi(1,1,inizja),nq*3,F,nnia,0d0,ffd,
434     +                          nq*3)
435c
436c                          Compute nuclear gradient of delrho
437c
438                           call xc_drhonuc(nq,ipol,inizia,ifinia,
439     D                          delrho,heschi,delchi,ff,ffd,
440     I                          ii,iat)
441
442                         endif
443                       endif
444                     endif
445                     call endtimer(monitor_rhocomp)
446  215              continue
447                endif
448              enddo
449            endif
450  220    continue
451  225    continue
452  230 continue
453      if(zapnegatives) then
454c
455c     Enforce results that are compatible with the laws of physics.
456c     I.e. Rho >= 0, Tau >= 0, Rho=0 ==> Grad(Rho)=0, and
457c     Tau=0 ==> Grad(Rho)=0. The reason for doing this is that density
458c     functionals have been designed to be valid for physically
459c     sensible data points. They may fail spectacularly for data points
460c     outside the physically sensible range. Alternatively to screening
461c     the data here it might be done in the functionals themselves.
462c     However, this requires a great deal of care to make sure that the
463c     right limits are taken.
464c
465c     Rho should be non-negative. Numerical errors in the construction
466c     may leave small negative values that are formally incorrect.
467c     Hence filter negative values out.
468c
469c     Also note that rho.eq.0 implies delrho=0, so this is enforced as
470c     well. Proof: Rho is a non-negative quantity. Hence if rho.eq.0
471c     then rho is at a minimum. A minimum being an extremum it follows
472c     that its gradient must be zero. QED.
473c
474      if (what.eq.0) then
475        do ii = 1, ipol
476          jj = i0 + ii
477          do n = 1, nq
478            if (rho(n,jj,1).le.0.0d0) then
479              rho(n,jj,1) = 0.0d0
480              if (grad) then
481                delrho(n,1,ii,1) = 0.0d0
482                delrho(n,2,ii,1) = 0.0d0
483                delrho(n,3,ii,1) = 0.0d0
484              endif
485            endif
486          enddo
487        enddo
488      else if (what.eq.1) then
489        do ii = 1, ipol
490          jj = ii
491          do n = 1, nq
492            if (rho(n,jj,npert).le.0.0d0) then
493              rho(n,jj,npert) = 0.0d0
494              if (grad) then
495                delrho(n,1,ii,npert) = 0.0d0
496                delrho(n,2,ii,npert) = 0.0d0
497                delrho(n,3,ii,npert) = 0.0d0
498              endif
499            endif
500          enddo
501        enddo
502      endif
503c
504c     Tau should be non-negative. Numerical errors in the construction
505c     may leave small negative values that can cause trouble in the
506c     density functionals. Hence filter negative values out.
507c
508c     Also note that tau.eq.0 implies delrho=0, so this is enforced as
509c     well. Proof: tau =1/2 \sum_i (\nabla\phi_i)\cdot(\nabla\phi_i)
510c     hence tau.eq.0 iff (\nabla\phi_i)=0 for all i. Delrho is
511c     defined as delrho = \sum_i\nabla\phi_i^2 hence if the gradient
512c     of the orbital is zero for all orbitals i then delrho must be
513c     zero as well. QED.
514c
515      if (kske) then
516        if (what.eq.0) then
517          do ii = 1, ipol
518            do n = 1, nq
519              if (ttau(n,ii,1).le.0.0d0) then
520                ttau(n,ii,1) = 0.0d0
521                if (grad) then
522                  delrho(n,1,ii,1) = 0.0d0
523                  delrho(n,2,ii,1) = 0.0d0
524                  delrho(n,3,ii,1) = 0.0d0
525                endif
526              endif
527            enddo
528          enddo
529        else if (what.eq.1) then
530          do ii = 1, ipol
531            do n = 1, nq
532              if (ttau(n,ii,npert+1).le.0.0d0) then
533                ttau(n,ii,npert+1) = 0.0d0
534                if (grad) then
535                  delrho(n,1,ii,npert+1) = 0.0d0
536                  delrho(n,2,ii,npert+1) = 0.0d0
537                  delrho(n,3,ii,npert+1) = 0.0d0
538                endif
539              endif
540            enddo
541          enddo
542        endif
543      endif
544      endif
545c
546      call starttimer(monitor_rhocomp2)
547      if(what.eq.0) then
548c
549c       Only construct total densities for regular case
550c
551         if (ipol.eq.2)then
552            call dcopy(nq, rho(1,2,1), 1, rho(1,1,1), 1)
553            call daxpy(nq, 1.d0, rho(1,3,1), 1, rho(1,1,1), 1)
554         endif
555      endif
556      if(what.eq.0) then
557         if(xcreplicated.and.dorepdm) then
558            g_dens(1)=g_keepd(1)
559            if(ipol.eq.2)  g_dens(2)=g_keepd(2)
560         endif
561      endif
562      call endtimer(monitor_rhocomp2)
563c
564 1688 continue
565      if(dftnbget) then
566         if (.not.ma_pop_stack(l_nz))
567     &        call errquit('xcrho:pop_stack failed', 13, MA_ERR)
568      endif
569
570      call endtimer(monitor_xcrho)
571      return
572      end
573c
574      subroutine xc_drhonuc(nq,ipol,n0,n1,
575     D     delrho,heschi,delchi,ff,ffd,
576     I     ii,iat)
577      implicit none
578      integer nq,ipol,ii,n0,n1,iat
579      double precision delrho(nq,3,ipol,3,*)
580      double precision heschi(nq,6,*)
581      double precision delchi(nq,3,*)
582      double precision ff(nq,*)
583      double precision ffd(nq,3,*)
584c
585      integer n,mu,mu1
586c
587      do mu=n0,n1
588         mu1=mu-n0+1
589         do n = 1, nq
590            delrho(n,1,ii,1,iat) = delrho(n,1,ii,1,iat) -
591     &           heschi(n,1,mu)*ff(n,mu1) -
592     -           delchi(n,1,mu)*ffd(n,1,mu1)
593            delrho(n,2,ii,1,iat) = delrho(n,2,ii,1,iat) -
594     &           heschi(n,2,mu)*ff(n,mu1) -
595     -           delchi(n,1,mu)*ffd(n,2,mu1)
596            delrho(n,3,ii,1,iat) = delrho(n,3,ii,1,iat) -
597     &           heschi(n,3,mu)*ff(n,mu1) -
598     -           delchi(n,1,mu)*ffd(n,3,mu1)
599            delrho(n,1,ii,2,iat) = delrho(n,1,ii,2,iat) -
600     &           heschi(n,2,mu)*ff(n,mu1) -
601     -           delchi(n,2,mu)*ffd(n,1,mu1)
602            delrho(n,2,ii,2,iat) = delrho(n,2,ii,2,iat) -
603     &           heschi(n,4,mu)*ff(n,mu1) -
604     -           delchi(n,2,mu)*ffd(n,2,mu1)
605            delrho(n,3,ii,2,iat) = delrho(n,3,ii,2,iat) -
606     &           heschi(n,5,mu)*ff(n,mu1) -
607     -           delchi(n,2,mu)*ffd(n,3,mu1)
608            delrho(n,1,ii,3,iat) = delrho(n,1,ii,3,iat) -
609     &           heschi(n,3,mu)*ff(n,mu1) -
610     -           delchi(n,3,mu)*ffd(n,1,mu1)
611            delrho(n,2,ii,3,iat) = delrho(n,2,ii,3,iat) -
612     &           heschi(n,5,mu)*ff(n,mu1) -
613     -           delchi(n,3,mu)*ffd(n,2,mu1)
614            delrho(n,3,ii,3,iat) = delrho(n,3,ii,3,iat) -
615     &           heschi(n,6,mu)*ff(n,mu1) -
616     -           delchi(n,3,mu)*ffd(n,3,mu1)
617         enddo
618      enddo
619      return
620      end
621      subroutine xc_dchiff(nq,n0,n1,
622     P     delrho,delchi,ff)
623      implicit none
624      integer nq,n0,n1
625      double precision delrho(nq,3)
626      double precision delchi(nq,3,*)
627      double precision ff(nq,*)
628c
629      integer n,mu,mu1
630c
631      do mu=n0,n1
632         mu1=mu-n0+1
633         do n = 1, nq
634            delrho(n,1) = delrho(n,1) + delchi(n,1,mu)*ff(n,mu1)
635            delrho(n,2) = delrho(n,2) + delchi(n,2,mu)*ff(n,mu1)
636            delrho(n,3) = delrho(n,3) + delchi(n,3,mu)*ff(n,mu1)
637         enddo
638      enddo
639      return
640      end
641      subroutine xc_dchiffp(nq,ipol2,n0,n1,
642     P     rho,delchi,ff,
643     M     ii,iat)
644      implicit none
645      integer nq,ipol2,n0,n1
646      double precision rho(nq,ipol2,3,*)
647      double precision delchi(nq,3,*)
648      double precision ff(nq,*)
649      integer ii,iat
650c
651      integer n,mu,mu1
652c
653      do mu=n0,n1
654         mu1=mu-n0+1
655         do n = 1, nq
656            rho(n,ii,1,iat) = rho(n,ii,1,iat)-delchi(n,1,mu)*ff(n,mu1)
657            rho(n,ii,2,iat) = rho(n,ii,2,iat)-delchi(n,2,mu)*ff(n,mu1)
658            rho(n,ii,3,iat) = rho(n,ii,3,iat)-delchi(n,3,mu)*ff(n,mu1)
659         enddo
660      enddo
661      return
662      end
663
664c@@@@@@@@@@@@@@@@@
665      integer function xc_rhoscreen(grad,ipol,natoms,npert,
666     I     iniz,
667     W     tol_rho,wmax,
668     O     nz,
669     R     rchi_atom,rdelchi_atom,rdens_atom)
670      implicit none
671#include "global.fh"
672      logical grad
673      integer ipol,natoms,npert
674      integer iniz(natoms) ! mapping of nbf_ao -> mbf
675      double precision tol_rho,wmax
676      double precision rchi_atom(*)
677      double precision rdelchi_atom(*)
678      double precision rdens_atom(natoms,natoms,ipol)
679c
680c      output
681c
682      integer nz(*)
683c
684      integer iat,jat,inizia,inizja,iat_in
685      integer ipert
686      double precision FUNC_MAXI,FUNC_MAXJ,DELFUNC_MAX,FUNC_MAX,
687     ,     P_MAXJ,P_MAXJ_A,P_MAXJ_B,P_MAXIJ
688      integer nonzero,ii
689      double precision dabsmax
690      external dabsmax
691c
692      FUNC_MAX = dabsmax(natoms,rchi_atom)
693      DELFUNC_MAX=0d0
694      if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom)
695
696      nonzero=0
697cscat
698#ifdef STAG
699      do iat_in = 1+ga_nodeid(), natoms+ga_nodeid()
700         iat=mod(iat_in,natoms)
701         if(iat.eq.0) iat=natoms
702#else
703      do iat = 1, natoms
704#endif
705         inizia = iniz(iat)
706         if (inizia.ne.0) then
707c
708c     screening parameters
709c
710            FUNC_MAXI = rchi_atom(iat)
711            if(grad)
712     .           FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat))
713            FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX)
714            if (ipol.gt.1)then
715               P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1))
716               P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2))
717               P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B)
718            else
719               P_MAXJ=0d0
720               do jat=1,iat
721                  if(iniz(jat).ne.0)
722     .                 P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1))
723               enddo
724            endif
725            if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.ge.tol_rho) then
726               do jat = 1, iat
727                  inizja = iniz(jat)
728                  if (inizja.ne.0) then
729c
730c     screening parameters
731c
732                     FUNC_MAXJ=rchi_atom(jat)
733                     if(grad)
734     .                    FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat))
735                     P_MAXIJ = rdens_atom(iat,jat,1)
736                     if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ,
737     &                    rdens_atom(iat,jat,2))
738                     if(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho)then
739c
740                        do ipert = 1,npert
741                           do ii = 1, ipol
742c
743c     screening parameters
744c
745                              P_MAXIJ = rdens_atom(iat,jat,ii)
746                              if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.
747     P                             ge.tol_rho) then
748                                 nonzero=nonzero+1
749                                 nz(nonzero)=(iat-1)*natoms+jat
750                              endif
751                           enddo
752                        enddo
753                     endif
754                  endif
755               enddo            ! loop jat
756            endif
757         endif
758      enddo                     ! loop iat
759      xc_rhoscreen=nonzero
760      return
761      end
762      subroutine xc_getdmblock(nz,nonz0,natoms,cetobfr,g_dens,
763     A     Pmat,nbhandl)
764      implicit none
765#include "dist.fh"
766#include "dft_fdist.fh"
767      integer nz(*)
768      integer natoms
769      integer nbhandl
770      integer g_dens
771      integer cetobfr(2,*)
772      integer nonz0
773      double precision pmat(*)
774c
775      integer ij0
776      integer iat0,jat0,ifirst,ilast,jfirst,jlast,nbfia
777c
778      ij0=nz(nonz0)
779      iat0=(ij0-1)/natoms+1
780      jat0=ij0-(iat0-1)*natoms
781      ifirst = cetobfr(1,iat0)
782      ilast = cetobfr(2,iat0)
783      jfirst = cetobfr(1,jat0)
784      jlast = cetobfr(2,jat0)
785      nbfia= ilast-ifirst+1
786      call starttimer(monitor_ganbget)
787      call ga_nbget(g_dens,
788     %     ifirst, ilast, jfirst, jlast, Pmat,nbfia,nbhandl)
789      call endtimer(monitor_ganbget)
790      return
791      end
792      subroutine xc_dmget(repdm, nbf_ld,
793     %                 ilo, ihi, jlo, jhi, Pmat,nbfia)
794      implicit none
795      integer  ilo, ihi, jlo, jhi, nbfia,nbf_ld
796      double precision pmat(nbfia,*)
797      double precision repdm(*)
798      integer ij,nnn
799      integer i,j
800c
801      nnn=nbfia
802      if(ilo.ne.jlo) then
803         do j=jlo,jhi
804            ij=((j-1)*(2*(nbf_ld+1)-j)+1)/2+ilo-j+1
805            call dcopy(nnn,repdm(ij),1,pmat(1,j-jlo+1),1)
806         enddo
807      else
808c diag block: copy only lower tr
809         do j=jlo,jhi
810            ij=((j-1)*(2*(nbf_ld+1)-j)+1)/2+1
811            call dcopy(nnn,repdm(ij),1,pmat(j-jlo+1,j-jlo+1),1)
812            nnn=nnn-1
813         enddo
814c     copy offdiag terms (aka transp)
815         do j=1,nbfia
816            do i=j+1,nbfia
817               pmat(j,i)=pmat(i,j)
818            enddo
819         enddo
820      endif
821      return
822      end
823c
824c    Essentially a clone of the old xc_rhogen.F without the new negative density screening
825c    We need this for transition densities which can be -ve
826c
827      subroutine td_rhogen(what,
828     T     tol_rho, basis, g_dens, max_at_bf,
829     N     natoms, curatoms, ncuratoms, npert,
830     I     ipol, nq, nbf, mbf, GRAD, ipol2,
831     &     F, Pmat, ff, ffd,
832     C     chi, delchi, heschi,
833     I     ibf, iniz, ifin,
834     &     rho, delrho, lap, rchi_atom,
835     &     rdelchi_atom, rdens_atom, cetobfr, wmax,
836     &     ttau, kske, dolap)
837c
838      implicit none
839c
840#include "errquit.fh"
841#include "mafdecls.fh"
842#include "dftpara.fh"
843#include "dist.fh"
844#include "dft_fdist.fh"
845
846      Logical GRAD !< [Input] .true. when using gradient corrected
847                   !< functional
848      Logical kske !< [Input] .true. when using kinetic energy density
849                   !< functional
850      integer what !< [Input] What to calculate
851                   !< - what=0: calculate density
852                   !< - what=1: calculate perturbed density (derivative
853                   !<   with respect to electron position)
854                   !< - what=2: calculate derivative with respect to
855                   !<   nuclear coordinates
856       logical dolap        ! true if lap is required [in]
857
858c what=0 dens
859c what=1 pert
860c what=2 nucder
861      integer basis        !< [Input] basis set handle
862      integer ipol         !< [Input] no. of spin channels
863      integer ipol2        !< [Input] no. of spin channels in density
864                           !< - 1 if closed shell
865                           !< - 3 if open shell
866      integer nbf          !< [Input] no. of basis functions
867      integer mbf          !< [Input] "restricted" no. of basis functions
868      integer max_at_bf    !< [Input] max no. bf per atom
869      integer nq           !< [Input] no. of quadrature points
870      integer natoms       !< [Input] no. of atoms
871      double precision wmax !< [Input] max weight
872      integer npert        !< [Input] number of perturbed densities
873      integer curatoms(*) !< [Input] indexing array for "active" atoms
874      integer ncuratoms        !< [Input] number of currently active atoms
875      integer g_dens(*) !< [Input] GA handle for DM
876      integer ibf(mbf)     !< [Input] mapping of nbf_ao -> mbf
877      integer iniz(natoms) !< [Input] mapping of nbf_ao -> mbf
878      integer ifin(natoms) !< [Input] mapping of nbf_ao -> mbf
879      double precision tol_rho !< [Input] accuracy for rho evaluation
880      double precision chi(nq,mbf)     !< [Input] function values
881      double precision delchi(nq,3,mbf)!< [Input] function gradients
882      double precision heschi(nq,6,mbf)!< [Input] function hessians
883c
884      double precision ttau(nq,ipol,*) !< [Output] Total Kohn-Sham K.E.density
885c
886       double precision lap(nq,ipol2,*)
887c
888      double precision delrho(nq,3,ipol,*) !< [Output] Derivative of density
889      double precision Pmat(*) !< [Scratch] scratch vector
890      double precision F(max_at_bf*max_at_bf) !< [Scratch] scratch vector
891      double precision ff(nq,*)   !< [Scratch] scratch array
892      double precision ffd(nq,3,*) !< [Scratch] scratch array
893      double precision rho(nq,ipol2,*) !< [Output] The density
894      double precision rchi_atom(natoms) !< [Input] Screening parameters
895      double precision rdelchi_atom(natoms) !< [Input] Screening parameters
896      double precision rdens_atom(natoms,natoms,ipol) !< [Input] Screening parameters
897      integer cetobfr(2,natoms) !< [Input] Centers to basis functions
898c
899c     local declarations
900c
901      integer i0, ii, mu, n, npol
902      integer ipert        ! perturbation loop index
903      integer iat, inizia, ifirst, ilast, nbfia, nnia
904      integer ifinia, ifinja
905      integer jat, inizja, jfirst, jlast, nbfja, nnja
906      integer iatcur, jatcur
907      double precision FUNC_MAX, DELFUNC_MAX, FUNC_MAXI, FUNC_MAXJ
908      double precision P_MAXJ, P_MAXJ_A, P_MAXJ_B, P_MAXIJ
909      double precision dabsmax
910      integer g_keepd(2)
911      integer nbhandl
912      integer jj
913      logical doffd,doitt
914      external dabsmax
915      external xc_rhoscreen
916      integer xc_rhoscreen
917      integer nonzero,nonz0
918      integer i_nz,l_nz
919      integer sizeblk, gindx
920#ifdef DEBUG
921      integer ga_nodeid
922      external ga_nodeid
923#endif
924      call starttimer(monitor_xcrho)
925      g_keepd(1)=0
926      g_keepd(2)=0
927c
928c     Evaluate the charge density and its gradient at each of the
929c     sampling points
930c
931      doffd=(what.eq.2.and.grad).or.(what.eq.0.and.kske).or.
932     & (what.eq.1.and.kske).or.(what.eq.0.and.dolap)
933
934      npol = (ipol*(ipol+1))/2
935c     to keep compilers quiet
936      iatcur=1
937      jatcur=1
938c
939      if(what.eq.0) then
940         call dcopy(nq*npol,0.D0,0,rho,1)
941         if (grad) call dcopy(3*nq*ipol,0.D0,0,delrho,1)
942         if (kske) call dcopy(nq*ipol,0.D0,0,ttau,1)  ! total
943         if (dolap) call dfill(nq*ipol2,0.D0,lap,1)  ! total
944c
945c     repl stuff
946c
947         g_keepd(1)=g_dens(1)
948         if(ipol.eq.2) g_keepd(2)=g_dens(2)
949         if(xcreplicated.and.dorepdm) then
950            g_dens(1)=g_repdm(1)
951            if(ipol.eq.2) g_dens(2)=g_repdm(2)
952         endif
953      elseif(what.eq.1) then
954         call dfill(nq*ipol*npert,0.D0,rho,1)
955         if (grad) call dfill(3*nq*ipol*npert,0.D0,delrho,1)
956         if (kske) call dfill(nq*ipol*npert,0.D0,ttau,1)   ! total
957      elseif(what.eq.2) then
958         call dfill(nq*ipol*3*ncuratoms,0.D0,rho,1)
959         if (grad)call dfill(3*nq*ipol*3*ncuratoms,0.D0,delrho,1)
960      else
961         call errquit('wrong what value',0,0)
962      endif
963c
964c     Screening is accomplished by:  p(r) <= |Xi(r)|*|Xj(r)|*|Dij|
965c     Xi(r) is screened on desired accuracy/max(|Xj(r)|)*max(|Dij|)
966c     Dij is screened on desired accuracy/max(|Xi(r)|)*max(|Xj(r)|)
967c
968      i0=ipol-1
969c
970      FUNC_MAX = dabsmax(natoms,rchi_atom)
971      DELFUNC_MAX=0d0
972      if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom)
973c
974      nonzero=0
975      if(dftnbget) then
976         if (.not.ma_push_get
977     &        (mt_int,(ipol*npert*natoms*(natoms+1))/2,
978     N        'nzmap map',l_nz,i_nz))
979     &        call errquit('xcrho:push_get failed', 13, MA_ERR)
980         nonzero=xc_rhoscreen(grad,ipol,natoms,npert,
981     I        iniz,
982     W        tol_rho,wmax,
983     O        int_mb(i_nz),
984     R        rchi_atom,rdelchi_atom,rdens_atom)
985
986         if(nonzero.eq.0) goto 1688
987c
988c     prefetch first DM block
989c
990         nonz0=1
991         call xc_getdmblock(int_mb(i_nz),nonz0,natoms,cetobfr,
992     G        g_dens(1),
993     A        Pmat,nbhandl)
994         endif
995
996      do 230 iat = 1, natoms
997         inizia = iniz(iat)
998         if (inizia.eq.0)goto 230
999         if(what.eq.2) then
1000            iatcur = curatoms(iat)
1001         endif
1002         ifinia = ifin(iat)
1003         ifirst = cetobfr(1,iat)
1004         ilast = cetobfr(2,iat)
1005         nbfia = ilast-ifirst+1
1006         nnia = ifinia-inizia+1
1007c
1008c        screening parameters
1009c
1010         FUNC_MAXI = rchi_atom(iat)
1011         if(grad)
1012     .   FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat))
1013         FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX)
1014         if(what.lt.2) then
1015            if (ipol.gt.1)then
1016            P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1))
1017            P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2))
1018            P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B)
1019         else
1020            P_MAXJ=0d0
1021            do jat=1,iat
1022               if(iniz(jat).ne.0)
1023     .              P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1))
1024            enddo
1025         endif
1026         if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.lt.tol_rho) goto 225
1027      endif
1028         do 220 jat = 1, iat
1029            inizja = iniz(jat)
1030            if (inizja.eq.0)goto 220
1031            if(what.eq.2) then
1032               jatcur = curatoms(jat)
1033               if (jatcur.eq.0.and.iatcur.eq.0) goto 220
1034            endif
1035            call starttimer(monitor_rscreen0)
1036            ifinja = ifin(jat)
1037            jfirst = cetobfr(1,jat)
1038            jlast = cetobfr(2,jat)
1039            nbfja = jlast-jfirst+1
1040            nnja = ifinja-inizja+1
1041c
1042c           screening parameters
1043c
1044            FUNC_MAXJ=rchi_atom(jat)
1045            if(grad)
1046     .      FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat))
1047            doitt=.true.
1048            if(what.lt.2) then
1049               P_MAXIJ = rdens_atom(iat,jat,1)
1050               if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ,
1051     &              rdens_atom(iat,jat,2))
1052               doitt=(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho)
1053            endif
1054            call endtimer(monitor_rscreen0)
1055            if (doitt) then
1056            do ii = 1, ipol
1057c
1058c              screening parameters
1059c
1060            if((what.gt.1).or.(rdens_atom(iat,jat,ii)*
1061     R              wmax*FUNC_MAXI*FUNC_MAXJ.ge.tol_rho)) then
1062c
1063c           Loop over perturbations
1064c
1065            do 215 ipert = 1,npert
1066               sizeblk=nbfia*nbfja
1067               call updist(monitor_size_ga_get, sizeblk)
1068c
1069               if(dftnbget) then
1070                  call starttimer(monitor_wait3)
1071                  call ga_nbwait(nbhandl)
1072                  call endtimer(monitor_wait3)
1073                  call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja,
1074     I                 ifirst, jfirst, ibf(inizia), ibf(inizja))
1075                  nonz0=nonz0+1
1076                  if((npert*ipol).gt.1) then
1077                     gindx=ipert+(ii-1)*npert+1
1078
1079                     if(gindx.gt.npert*ipol)gindx=mod(gindx,npert*ipol)
1080                  else
1081                     gindx=1
1082                  endif
1083                  if(nonz0.le.nonzero)
1084     C                 call xc_getdmblock(int_mb(i_nz),nonz0,natoms,
1085     C                 cetobfr,g_dens(gindx),
1086     A                    Pmat,nbhandl)
1087               if (wmax*FUNC_MAXI*FUNC_MAXJ*
1088     .              dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215
1089               else
1090                  call starttimer(monitor_gaget)
1091                  if(truerepdm) then
1092                  call xc_dmget(dbl_mb(k_repdm(ipert+(ii-1)*npert)),
1093     &                    nbf_ld,
1094     %                 ifirst, ilast, jfirst, jlast, Pmat,nbfia)
1095                  else
1096                  call ga_get(g_dens(ipert+(ii-1)*npert),
1097     %                 ifirst, ilast, jfirst, jlast, Pmat,nbfia)
1098                  endif
1099                  call endtimer(monitor_gaget)
1100               if (wmax*FUNC_MAXI*FUNC_MAXJ*
1101     .              dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215
1102               call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja, ifirst,
1103     &                      jfirst, ibf(inizia), ibf(inizja))
1104               endif
1105
1106c
1107               call starttimer(monitor_rhocomp)
1108               if(iat.ne.jat) call dscal(nnia*nnja,2d0,F,1)
1109c
1110c              Compute Xiat(r)*Xjat(r)*Diat,jat
1111c
1112               call dgemm('n','n',nq,nnja,nnia,1d0,
1113     A              chi(1,inizia),nq,F,nnia,0d0,ff,nq)
1114               if(what.lt.2) then
1115                  jj=i0+ii
1116                  if(what.eq.1) jj=ii
1117                  do mu=inizja,ifinja
1118                     do n=1,nq
1119                        rho(n,jj,ipert) = rho(n,jj,ipert) +
1120     P                       chi(n,mu)*ff(n,mu-inizja+1)
1121                     enddo
1122                  enddo
1123               endif
1124               if(doffd) then
1125                  call dgemm('n','n',nq*3,nnja,nnia,1d0,
1126     A                 delchi(1,1,inizia),nq*3,F,nnia,0d0,ffd,nq*3)
1127c
1128c              build tau for meta GGA
1129c
1130                  if(kske) then
1131                     do mu=inizja,ifinja
1132                        do n=1,nq
1133                           ttau(n,ii,ipert) = ttau(n,ii,ipert)+0.5d0*(
1134     &                       delchi(n,1,mu)*ffd(n,1,mu-inizja+1)+
1135     &                       delchi(n,2,mu)*ffd(n,2,mu-inizja+1)+
1136     &                       delchi(n,3,mu)*ffd(n,3,mu-inizja+1))
1137                        enddo
1138                     enddo
1139                  endif
1140
1141                  if (dolap) then
1142                     do mu=inizja,ifinja
1143                        do n=1,nq
1144                           ! total
1145                           lap(n,ii,ipert) = lap(n,ii,ipert) + 2d0*
1146     &                          (delchi(n,1,mu)*ffd(n,1,mu-inizja+1) +
1147     &                          delchi(n,2,mu)*ffd(n,2,mu-inizja+1) +
1148     &                          delchi(n,3,mu)*ffd(n,3,mu-inizja+1) +
1149     &                          heschi(n,1,mu)*ff(n,mu-inizja+1) +
1150     &                          heschi(n,4,mu)*ff(n,mu-inizja+1) +
1151     &                          heschi(n,6,mu)*ff(n,mu-inizja+1))
1152                        enddo
1153                     enddo
1154                  endif
1155               endif
1156c
1157               if((what.eq.2.and.jatcur.ne.0).or.
1158c
1159c     We need the "grad" code at zero order in the nuclear deriv case,
1160c     but we can skip this part if iat is not active
1161c
1162     O              (what.ne.2.and.grad)) then
1163                  if(what.eq.0.or.what.eq.1) then
1164                     call xc_dchiff(nq,inizja,ifinja,
1165     P                    delrho(1,1,ii,ipert),delchi,
1166     F                    ff)
1167                  elseif(what.eq.2) then
1168                     call xc_dchiffp(nq,ipol2,inizja,ifinja,
1169     P                    rho,delchi,ff,
1170     M                    ii,jat)
1171                     if (grad) then
1172c
1173c     Compute nuclear gradient of delrho
1174c
1175                        call xc_drhonuc(nq,ipol,inizja,ifinja,
1176     D                       delrho,heschi,delchi,ff,ffd,
1177     I                       ii,jat)
1178                     endif
1179                  endif
1180               endif
1181c
1182               if((what.eq.2.and.iatcur.ne.0).or.
1183     O              (what.ne.2.and.grad)) then
1184c
1185c                 Compute delXiat(r)*Xjat(r)*Diat,jat
1186c
1187               call dgemm('n','t',nq,nnia,nnja,1d0,
1188     A              chi(1,inizja),nq,F,nnia,0d0,ff,nq)
1189                  if(what.lt.2) then
1190                     call xc_dchiff(nq,inizia,ifinia,
1191     P                    delrho(1,1,ii,ipert),delchi,
1192     F                    ff)
1193                  elseif(what.eq.2) then
1194                     call xc_dchiffp(nq,ipol2,inizia,ifinia,
1195     P                    rho,delchi,ff,
1196     M                    ii,iat)
1197                     if(grad) then
1198                        call dgemm('n','t',nq*3,nnia,nnja,1d0,
1199     A                    delchi(1,1,inizja),nq*3,F,nnia,0d0,ffd,nq*3)
1200c
1201c     Compute nuclear gradient of delrho
1202c
1203                        call xc_drhonuc(nq,ipol,inizia,ifinia,
1204     D                       delrho,heschi,delchi,ff,ffd,
1205     I                       ii,iat)
1206
1207                     endif
1208                  endif
1209            endif
1210               call endtimer(monitor_rhocomp)
1211  215       continue
1212         endif
1213         enddo
1214         endif
1215  220    continue
1216  225    continue
1217  230 continue
1218
1219      call starttimer(monitor_rhocomp2)
1220      if(what.eq.0) then
1221c
1222c     Only construct total densities for regular case
1223c
1224         if (ipol.eq.2)then
1225            call dcopy(nq, rho(1,2,1), 1, rho(1,1,1), 1)
1226            call daxpy(nq, 1.d0, rho(1,3,1), 1, rho(1,1,1), 1)
1227         endif
1228      endif
1229      if(what.eq.0) then
1230         if(xcreplicated.and.dorepdm) then
1231            g_dens(1)=g_keepd(1)
1232            if(ipol.eq.2)  g_dens(2)=g_keepd(2)
1233         endif
1234      endif
1235      call endtimer(monitor_rhocomp2)
1236c
1237 1688 continue
1238      if(dftnbget) then
1239         if (.not.ma_pop_stack(l_nz))
1240     &        call errquit('xcrho:pop_stack failed', 13, MA_ERR)
1241      endif
1242
1243      call endtimer(monitor_xcrho)
1244
1245      return
1246      end
1247
1248      subroutine xc_delrhodotr(nq,ipol,qxyz,rho,delrho,delrhoq)
1249      implicit none
1250      integer nq,ipol                 ![in]
1251      double precision rho(nq,*)      ! [in]
1252      double precision delrho(nq,3,*) ! [in]
1253      double precision qxyz(3,*)      ! [in]
1254      double precision delrhoq(nq,*)  ! [out]
1255c
1256      integer ii,n
1257
1258      if(ipol.eq.1) then
1259         do n = 1, nq
1260            delrhoq(n,1) = 3d0*rho(n,1)+
1261     D           delrho(n,1,1)*qxyz(1,n) +
1262     D           delrho(n,2,1)*qxyz(2,n) +
1263     D           delrho(n,3,1)*qxyz(3,n)
1264         enddo
1265      else
1266      do ii=2,3
1267         do n = 1, nq
1268            delrhoq(n,ii-1) = 3d0*rho(n,ii)+
1269     D           delrho(n,1,ii-1)*qxyz(1,n) +
1270     D           delrho(n,2,ii-1)*qxyz(2,n) +
1271     D           delrho(n,3,ii-1)*qxyz(3,n)
1272         enddo
1273      enddo
1274      endif
1275      return
1276      end
1277C>
1278C> @}
1279
1280
1281
1282
1283