TABLE OF CONTENTS


ABINIT/mkrho [ Functions ]

[ Top ] [ Functions ]

NAME

 mkrho

FUNCTION

 Depending on option argument value:
 --Compute charge density rho(r) and rho(G) in electrons/bohr**3
   from input wavefunctions, band occupations, and k point wts.
 --Compute kinetic energy density tau(r) and tau(G) in bohr**-5
   from input wavefunctions, band occupations, and k point wts.
 --Compute a given element of the kinetic energy density tensor
   tau_{alpha,beta}(r) and tau_{alpha,beta}(G) in bohr**-5
   from input wavefunctions, band occupations, and k point wts.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LSI, AR, MB)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mcg)=wf in G space
  dtset <type(dataset_type)>=all input variables for this dataset
   | istwfk(nkpt)=input option parameter that describes the storage of wfs
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs
   | mkmem=Number of k points treated by this node
   | mpw=maximum allowed value for npw
   | nband(nkpt*nsppol)=number of bands to be included in summation
   |  at each k point for each spin channel
   | nfft=(effective) number of FFT grid points (for this processor)
   | ngfft(18)=contain all needed information about 3D FFT,
   |  see ~abinit/doc/variables/vargs.htm#ngfft
   | nkpt=number of k points
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in group (at least 1 for identity)
   | symafm(nsym)=(anti)ferromagnetic part of symmetry operations
   | symrel(3,3,nsym)=symmetry matrices in real space (integers)
   | wtk(nkpt)=k point weights (they sum to 1.0)
  gprimd(3,3)=dimensional reciprocal space primitive translations
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  kg(3,mpw*mkmem)=reduced planewave coordinates
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mpi_enreg=informations about MPI parallelization
  npwarr(nkpt)=number of planewaves and boundary planewaves at each k
  occ(mband*nkpt*nsppol)=
          occupation numbers for each band (usually 2.0) at each k point
  option if 0: compute rhor (electron density)
         if 1: compute taur (kinetic energy density)
               (i.e. Trace over the kinetic energy density tensor)
         if 2: compute taur_{alpha,beta} !!NOT YET IMPLEMENTED
               (a given element of the kinetic energy density tensor)
  paw_dmft  <type(paw_dmft_type)>= paw+dmft related data
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  rprimd(3,3)=dimensional real space primitive translations
  tim_mkrho=timing code of the calling routine(can be set to 0 if not attributed)
  ucvol=unit cell volume (Bohr**3)
  wvl_den <type(wvl_denspot_type)>=density informations for wavelets
  wvl_wfs <type(wvl_projector_type)>=wavefunctions informations for wavelets

OUTPUT

 rhog(2,nfft)=total electron density in G space
 rhor(nfft,nspden)=electron density in r space
   (if spin polarized, array contains total density in first half and spin-up density in second half)
   (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2)

PARENTS

      afterscfloop,energy,gstate,respfn,scfcv,vtorho

CHILDREN

      bandfft_kpt_set_ikpt,fftpac,fourwf,prep_fourwf,prtrhomxmn
      sphereboundary,symrhg,timab,wrtout,wvl_mkrho,xmpi_sum

SOURCE

 79 #if defined HAVE_CONFIG_H
 80 #include "config.h"
 81 #endif
 82 
 83 #include "abi_common.h"
 84 
 85 subroutine mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
 86 &                rhog,rhor,rprimd,tim_mkrho,ucvol,wvl_den,wvl_wfs,&
 87 &                option) !optional
 88 
 89  use defs_basis
 90  use defs_abitypes
 91  use defs_wvltypes
 92  use m_profiling_abi
 93  use m_xmpi
 94  use m_errors
 95 
 96  use m_hamiltonian,  only : gs_hamiltonian_type
 97  use m_bandfft_kpt,  only : bandfft_kpt_set_ikpt
 98  use m_paw_dmft,     only : paw_dmft_type
 99 
100 !This section has been created automatically by the script Abilint (TD).
101 !Do not modify the following lines by hand.
102 #undef ABI_FUNC
103 #define ABI_FUNC 'mkrho'
104  use interfaces_14_hidewrite
105  use interfaces_18_timing
106  use interfaces_32_util
107  use interfaces_52_fft_mpi_noabirule
108  use interfaces_53_ffts
109  use interfaces_66_wfs
110  use interfaces_67_common, except_this_one => mkrho
111 !End of the abilint section
112 
113  implicit none
114 
115 !Arguments ------------------------------------
116 !scalars
117  integer,intent(in) :: mcg,tim_mkrho
118  integer,intent(in),optional :: option
119  real(dp),intent(in) :: ucvol
120  type(MPI_type),intent(inout) :: mpi_enreg
121  type(dataset_type),intent(in) :: dtset
122  type(paw_dmft_type), intent(in)  :: paw_dmft
123  type(wvl_wf_type),intent(inout) :: wvl_wfs
124  type(wvl_denspot_type), intent(inout) :: wvl_den
125 !no_abirules
126 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
127  integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,  &
128 &               (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
129  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
130  real(dp), intent(in) :: gprimd(3,3)
131  real(dp), intent(in) :: cg(2,mcg)
132  real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
133 !nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
134  real(dp), intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym),  &
135 &                                 (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
136  real(dp), intent(in) :: rprimd(3,3)
137  real(dp), intent(out) :: rhor(dtset%nfft,dtset%nspden),rhog(2,dtset%nfft)
138 
139 !Local variables-------------------------------
140 !scalars
141  integer,save :: nskip=0
142  integer :: alpha,use_nondiag_occup_dmft,bdtot_index,beta,blocksize,iband,iband1,ibandc1,ib,iblock,icg,ierr
143  integer :: ifft,ikg,ikpt,ioption,ipw,ipwsp,ishf,ispden,ispinor,ispinor_index
144  integer :: isppol,istwf_k,jspinor_index
145  integer :: me,my_nspinor,n1,n2,n3,n4,n5,n6,nalpha,nband_k,nbandc1,nbdblock,nbeta
146  integer :: ndat,nfftot,npw_k,spaceComm,tim_fourwf
147  real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp
148  real(dp) :: weight,weight_i
149  character(len=500) :: message
150 !arrays
151  integer,allocatable :: gbound(:,:),kg_k(:,:)
152  logical :: locc_test,nspinor1TreatedByThisProc,nspinor2TreatedByThisProc
153  real(dp) :: dummy(2,1),tsec(2)
154  real(dp),allocatable :: cwavef(:,:,:),cwavefb(:,:,:),cwavef_x(:,:)
155  real(dp),allocatable :: cwavef_y(:,:),cwavefb_x(:,:),cwavefb_y(:,:),kg_k_cart_block(:)
156  real(dp),allocatable :: occ_k(:),rhoaug(:,:,:),rhoaug_down(:,:,:)
157  real(dp),allocatable :: rhoaug_mx(:,:,:),rhoaug_my(:,:,:),rhoaug_up(:,:,:)
158  real(dp),allocatable :: taur_alphabeta(:,:,:,:),wfraug(:,:,:,:)
159 
160 ! *************************************************************************
161 
162  DBG_ENTER("COLL")
163 
164  call timab(790+tim_mkrho,1,tsec)
165  call timab(799,1,tsec)
166 
167  if(.not.(present(option))) then
168    ioption=0
169  else
170    ioption=option
171  end if
172 
173  if(ioption/=0.and.paw_dmft%use_sc_dmft==1) then
174    message = ' option argument value of this routines should be 0 if usedmft=1. '
175    MSG_ERROR(message)
176  end if
177  if(paw_dmft%use_sc_dmft/=0) then
178    nbandc1=(paw_dmft%mbandc-1)*paw_dmft%use_sc_dmft+1
179  else
180    nbandc1=1
181  end if
182  use_nondiag_occup_dmft=0
183 
184 !if(dtset%nspinor==2.and.paw_dmft%use_sc_dmft==1) then
185 !write(message, '(a,a,a,a)' )ch10,&
186 !&   ' mkrho : ERROR -',ch10,&
187 !&   '  nspinor argument value of this routines should be 1 if usedmft=1. '
188 !call wrtout(std_out,message,'COLL')
189 !end if
190 
191  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
192  if (mpi_enreg%paral_spinor==0) then
193    ispinor_index=1;jspinor_index=1
194    nspinor1TreatedByThisProc=.true.
195    nspinor2TreatedByThisProc=(dtset%nspinor==2)
196  else
197    ispinor_index=mpi_enreg%me_spinor+1;jspinor_index=3-ispinor_index
198    nspinor1TreatedByThisProc=(mpi_enreg%me_spinor==0)
199    nspinor2TreatedByThisProc=(mpi_enreg%me_spinor==1)
200  end if
201 
202 !Set local variable which depend on option argument
203 
204 !nalpha*nbeta is the number of element of the kinetic energy density tensor
205 !to be computed in the irreducible Brillouin Zone (BZ) to get the result in the full BZ.
206 !In case of electron density calculation, nalpha=nbeta=1
207  select case (ioption)
208  case (0)
209    nalpha = 1
210    nbeta = 1
211  case (1)
212    nalpha = 3
213    nbeta = 1
214    ABI_ALLOCATE(taur_alphabeta,(dtset%nfft,dtset%nspden,3,1))
215  case (2)
216    nalpha = 3
217    nbeta = 3
218    ABI_ALLOCATE(taur_alphabeta,(dtset%nfft,dtset%nspden,3,3))
219  case default
220    MSG_BUG(' ioption argument value should be 0,1 or 2.')
221  end select
222 
223 !Init me
224  me=mpi_enreg%me_kpt
225 !zero the charge density array in real space
226 !$OMP PARALLEL DO COLLAPSE(2)
227  do ispden=1,dtset%nspden
228    do ifft=1,dtset%nfft
229      rhor(ifft,ispden)=zero
230    end do
231  end do
232 
233 !WVL - Branching with a separate mkrho procedure
234 !in wavelet.
235  if (dtset%usewvl == 1) then
236    select case(ioption)
237    case (0)
238      call wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den)
239      return
240    case (1)
241 !      call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
242      message = ' Sorry, kinetic energy density (taur) is not yet implemented in wavelet formalism.'
243      MSG_ERROR(message)
244    case (2)
245 !      call wvl_mkrho(dtset, mpi_enreg, occ, rhor, wvl_wfs, wvl_den)
246      message = '  Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented in wavelet formalism.'
247      MSG_BUG(message)
248    end select
249  end if
250 !WVL - Following is done in plane waves.
251 
252 !start loop over alpha and beta
253 
254  do alpha=1,nalpha
255    do beta=1,nbeta
256 
257 !    start loop over spin and k points
258      bdtot_index=0
259      icg=0
260 
261 !    n4,n5,n6 are FFT dimensions, modified to avoir cache trashing
262      n1 = dtset%ngfft(1) ; n2 = dtset%ngfft(2) ; n3 = dtset%ngfft(3)
263      n4 = dtset%ngfft(4) ; n5 = dtset%ngfft(5) ; n6 = dtset%ngfft(6)
264      ndat = 1 ; if (mpi_enreg%paral_kgb==1) ndat = mpi_enreg%bandpp
265      ABI_ALLOCATE(cwavef,(2,dtset%mpw,my_nspinor))
266      ABI_ALLOCATE(rhoaug,(n4,n5,n6))
267      ABI_ALLOCATE(wfraug,(2,n4,n5,n6*ndat))
268      ABI_ALLOCATE(cwavefb,(2,dtset%mpw*paw_dmft%use_sc_dmft,my_nspinor))
269      if(dtset%nspden==4) then
270        ABI_ALLOCATE(rhoaug_up,(n4,n5,n6))
271        ABI_ALLOCATE(rhoaug_down,(n4,n5,n6))
272        ABI_ALLOCATE(rhoaug_mx,(n4,n5,n6))
273        ABI_ALLOCATE(rhoaug_my,(n4,n5,n6))
274        rhoaug_up(:,:,:)=zero
275        rhoaug_down(:,:,:)=zero
276        rhoaug_mx(:,:,:)=zero
277        rhoaug_my(:,:,:)=zero
278      end if
279 
280      do isppol=1,dtset%nsppol
281        ikg=0
282 
283        rhoaug(:,:,:)=zero
284        do ikpt=1,dtset%nkpt
285 
286          nband_k = dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
287          npw_k=npwarr(ikpt)
288          istwf_k = dtset%istwfk(ikpt)
289 
290          if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
291            bdtot_index=bdtot_index+nband_k
292            cycle
293          end if
294 
295          ABI_ALLOCATE(gbound,(2*dtset%mgfft+8,2))
296          ABI_ALLOCATE(kg_k,(3,npw_k))
297 
298          kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
299          call sphereboundary(gbound,istwf_k,kg_k,dtset%mgfft,npw_k)
300 
301 !        Loop over bands to fft and square for rho(r)
302 !        Shoulb be changed to treat bands by batch always
303 
304          if(mpi_enreg%paral_kgb /= 1) then  ! Not yet parallelized on spinors
305            do iband=1,nband_k
306 !            if(paw_dmft%use_sc_dmft==1) then
307 !            write(std_out,*) 'iband  ',iband,occ(iband+bdtot_index),paw_dmft%occnd(iband,iband,ikpt,isppol)
308 !            else
309 !            write(std_out,*) 'iband  ',iband,occ(iband+bdtot_index)
310 !            endif
311              if(mpi_enreg%paralbd==1)then
312                if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me)) cycle
313              end if
314              do ibandc1=1,nbandc1 ! in case of DMFT
315 !              Check if DMFT and only treat occupied states (check on occ.)
316                if(paw_dmft%use_sc_dmft == 1) then
317                  iband1 = paw_dmft%include_bands(ibandc1)
318                  if(paw_dmft%band_in(iband)) then
319                    if(.not. paw_dmft%band_in(iband1))  stop
320                    use_nondiag_occup_dmft = 1
321                    locc_test = abs(paw_dmft%occnd(1,iband,iband1,ikpt,isppol)) +&
322 &                   abs(paw_dmft%occnd(2,iband,iband1,ikpt,isppol))>tol8
323 !                  write(std_out,*) "mkrho,ikpt,iband,use_occnd",ikpt,iband
324                  else
325                    use_nondiag_occup_dmft = 0
326                    locc_test = abs(occ(iband+bdtot_index))>tol8
327                    if(ibandc1 /=1 .and. .not. paw_dmft%band_in(iband)) cycle
328                  end if
329                else
330                  use_nondiag_occup_dmft = 0
331                  locc_test = abs(occ(iband+bdtot_index))>tol8
332                end if
333 
334                if (locc_test) then
335 !                Obtain Fourier transform in fft box and accumulate the density or the kinetic energy density
336 !                Not yet parallise on nspinor if paral_kgb non equal to 1
337                  ipwsp=(iband-1)*npw_k*my_nspinor +icg
338                  cwavef(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k)
339                  if (my_nspinor==2) cwavef(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k)
340 
341                  if(ioption==1)then
342 !                  Multiplication by 2pi i (k+G)_alpha
343                    gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
344                    kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
345                    do ispinor=1,my_nspinor
346                      do ipw=1,npw_k
347                        kg_k_cart=gp2pi1*kg_k(1,ipw)+gp2pi2*kg_k(2,ipw)+gp2pi3*kg_k(3,ipw)+kpt_cart
348                        cwftmp=-cwavef(2,ipw,ispinor)*kg_k_cart
349                        cwavef(2,ipw,ispinor)=cwavef(1,ipw,ispinor)*kg_k_cart
350                        cwavef(1,ipw,ispinor)=cwftmp
351                      end do
352                    end do
353                  else if(ioption==2)then
354                    message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.'
355                    MSG_ERROR(message)
356                  end if
357 
358 !                Non diag occupation in DMFT.
359                  if(use_nondiag_occup_dmft==1) then
360                    ipwsp=(iband1-1)*npw_k*my_nspinor +icg
361                    cwavefb(:,1:npw_k,1)=cg(:,1+ipwsp:ipwsp+npw_k)
362                    if (my_nspinor==2) cwavefb(:,1:npw_k,2)=cg(:,ipwsp+npw_k+1:ipwsp+2*npw_k)
363                    weight  =paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
364                    if(dtset%nspinor==1) weight_i=zero
365                    if(dtset%nspinor==2) weight_i=paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*dtset%wtk(ikpt)/ucvol
366 
367                  else
368                    weight=occ(iband+bdtot_index)*dtset%wtk(ikpt)/ucvol
369                    weight_i=weight
370                  end if
371 
372                  if(mpi_enreg%paralbd==0) tim_fourwf=3
373                  if(mpi_enreg%paralbd==1) tim_fourwf=6
374 
375 !                The same section of code is also found in vtowfk.F90 : should be rationalized !
376 
377                  call fourwf(1,rhoaug,cwavef(:,:,1),dummy,wfraug,gbound,gbound,&
378 &                 istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
379 &                 npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,&
380 &                 use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,1),&
381 &                 use_gpu_cuda=dtset%use_gpu_cuda)
382 
383 
384                  if(dtset%nspinor==2)then
385 !                  DEBUG GZ !To obtain a x-directed magnetization(test)
386 !                  cwavef1(1,1:npw_k)=-cwavef(2,1:npw_k)
387 !                  cwavef1(2,1:npw_k)= cwavef(1,1:npw_k)
388 !                  ENDDEBUG
389 
390                    if(dtset%nspden==1) then
391 !                    We need only the total density : accumulation continues on top of rhoaug
392 
393                      call fourwf(1,rhoaug,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
394 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
395 &                     npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,&
396 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda)
397 
398 
399                    else if(dtset%nspden==4) then
400 !                    Build the four components of rho. We use only norm quantities and, so fourwf.
401 !$\sum_{n} f_n \Psi^{* \alpha}_n \Psi^{\alpha}_n =\rho^{\alpha \alpha}$
402 !$\sum_{n} f_n (\Psi^{1}+\Psi^{2})^*_n (\Psi^{1}+\Psi^{2})_n=rho+m_x$
403 !$\sum_{n} f_n (\Psi^{1}-i \Psi^{2})^*_n (\Psi^{1}-i \Psi^{2})_n=rho+m_y$
404                      ABI_ALLOCATE(cwavef_x,(2,npw_k))
405                      ABI_ALLOCATE(cwavef_y,(2,npw_k))
406                      ABI_ALLOCATE(cwavefb_x,(2,npw_k*paw_dmft%use_sc_dmft))
407                      ABI_ALLOCATE(cwavefb_y,(2,npw_k*paw_dmft%use_sc_dmft))
408 !$(\Psi^{1}+\Psi^{2})$
409                      cwavef_x(:,:)=cwavef(:,1:npw_k,1)+cwavef(:,1:npw_k,2)
410 !$(\Psi^{1}-i \Psi^{2})$
411                      cwavef_y(1,:)=cwavef(1,1:npw_k,1)+cwavef(2,1:npw_k,2)
412                      cwavef_y(2,:)=cwavef(2,1:npw_k,1)-cwavef(1,1:npw_k,2)
413                      if(use_nondiag_occup_dmft==1) then
414                        cwavefb_x(:,:)=cwavefb(:,1:npw_k,1)+cwavefb(:,1:npw_k,2)
415                        cwavefb_y(1,:)=cwavefb(1,1:npw_k,1)+cwavefb(2,1:npw_k,2)
416                        cwavefb_y(2,:)=cwavefb(2,1:npw_k,1)-cwavefb(1,1:npw_k,2)
417                      end if
418                      rhoaug_up(:,:,:)=rhoaug(:,:,:) !Already computed
419 
420                      call fourwf(1,rhoaug_down,cwavef(:,:,2),dummy,wfraug,gbound,gbound,&
421 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
422 &                     npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,&
423 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb(:,:,2),use_gpu_cuda=dtset%use_gpu_cuda)
424 
425                      call fourwf(1,rhoaug_mx,cwavef_x,dummy,wfraug,gbound,gbound,&
426 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
427 &                     npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,&
428 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_x,use_gpu_cuda=dtset%use_gpu_cuda)
429 
430                      call fourwf(1,rhoaug_my,cwavef_y,dummy,wfraug,gbound,gbound,&
431 &                     istwf_k,kg_k,kg_k,dtset%mgfft,mpi_enreg,1,dtset%ngfft,&
432 &                     npw_k,1,n4,n5,n6,1,dtset%paral_kgb,tim_fourwf,weight,weight_i,&
433 &                     use_ndo=use_nondiag_occup_dmft,fofginb=cwavefb_y,use_gpu_cuda=dtset%use_gpu_cuda)
434 
435                      ABI_DEALLOCATE(cwavef_x)
436                      ABI_DEALLOCATE(cwavef_y)
437                      ABI_DEALLOCATE(cwavefb_x)
438                      ABI_DEALLOCATE(cwavefb_y)
439 
440                    end if ! dtset%nspden/=4
441 
442 
443                  end if
444 
445                else
446 !                Accumulate the number of one-way 3D ffts skipped
447                  nskip=nskip+1
448                end if ! abs(occ(iband+bdtot_index))>tol8
449 !              End loop on iband
450              end do ! iband1=1,(nband_k-1)*paw_dmft%use_sc_dmft+1
451            end do ! iband=1,nband_k
452 
453          else !paral_kgb==1
454            if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
455 
456            call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
457 
458            nbdblock=nband_k/(mpi_enreg%nproc_band * mpi_enreg%bandpp)
459            blocksize=nband_k/nbdblock
460            if(allocated(cwavef))  then
461              ABI_DEALLOCATE(cwavef)
462            end if
463            ABI_ALLOCATE(cwavef,(2,npw_k*blocksize,dtset%nspinor))
464            if(ioption==1)  then
465              ABI_ALLOCATE(kg_k_cart_block,(npw_k))
466            end if
467            ABI_ALLOCATE(occ_k,(nband_k))
468            occ_k(:)=occ(bdtot_index+1:bdtot_index+nband_k)
469 
470            do iblock=1,nbdblock
471              if (dtset%nspinor==1) then
472                cwavef(:,1:npw_k*blocksize,1)=cg(:,1+(iblock-1)*npw_k*blocksize+icg:iblock*npw_k*blocksize+icg)
473              else
474                if (mpi_enreg%paral_spinor==0) then
475                  ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg
476                  do ib=1,blocksize
477                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,1)=cg(:,1+(2*ib-2)*npw_k+ishf:(2*ib-1)*npw_k+ishf)
478                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,2)=cg(:,1+(2*ib-1)*npw_k+ishf:ib*2*npw_k+ishf)
479                  end do
480                else
481                  ishf=(iblock-1)*npw_k*my_nspinor*blocksize+icg
482                  do ib=1,blocksize
483                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,ispinor_index)=&
484 &                   cg(:,1+(ib-1)*npw_k+ishf:ib*npw_k+ishf)
485                    cwavef(:,(ib-1)*npw_k+1:ib*npw_k,jspinor_index)=zero
486                  end do
487                  call xmpi_sum(cwavef,mpi_enreg%comm_spinor,ierr)
488                end if
489              end if
490              if(ioption==1)then
491 !              Multiplication by 2pi i (k+G)_alpha
492                gp2pi1=gprimd(alpha,1)*two_pi ; gp2pi2=gprimd(alpha,2)*two_pi ; gp2pi3=gprimd(alpha,3)*two_pi
493                kpt_cart=gp2pi1*dtset%kptns(1,ikpt)+gp2pi2*dtset%kptns(2,ikpt)+gp2pi3*dtset%kptns(3,ikpt)
494                kg_k_cart_block(1:npw_k)=gp2pi1*kg_k(1,1:npw_k)+gp2pi2*kg_k(2,1:npw_k)+gp2pi3*kg_k(3,1:npw_k)+kpt_cart
495                do ib=1,blocksize
496                  do ipw=1,npw_k
497                    cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
498                    cwavef(2,ipw,1)=cwavef(1,ipw+(ib-1)*npw_k,1)*kg_k_cart_block(ipw)
499                    cwavef(1,ipw,1)=cwftmp
500                    if (my_nspinor==2) then
501                      cwftmp=-cwavef(2,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
502                      cwavef(2,ipw,2)=cwavef(1,ipw+(ib-1)*npw_k,2)*kg_k_cart_block(ipw)
503                      cwavef(1,ipw,2)=cwftmp
504                    end if
505                  end do
506                end do
507              else if(ioption==2)then
508                message = '  Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.'
509                MSG_ERROR(message)
510              end if
511 
512              call timab(538,1,tsec)
513              if (nspinor1TreatedByThisProc) then
514                call prep_fourwf(rhoaug,blocksize,cwavef(:,:,1),wfraug,iblock,istwf_k,dtset%mgfft,mpi_enreg,&
515 &               nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
516 &               dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
517              end if
518              call timab(538,2,tsec)
519              if(dtset%nspinor==2)then
520                if (dtset%nspden==1) then
521                  if (nspinor2TreatedByThisProc) then
522                    call prep_fourwf(rhoaug,blocksize,cwavef(:,:,2),wfraug,&
523 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
524 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
525 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
526                  end if
527                else if(dtset%nspden==4 ) then
528                  ABI_ALLOCATE(cwavef_x,(2,npw_k*blocksize))
529                  ABI_ALLOCATE(cwavef_y,(2,npw_k*blocksize))
530                  cwavef_x(:,:)=cwavef(:,:,1)+cwavef(:,:,2)
531                  cwavef_y(1,:)=cwavef(1,:,1)+cwavef(2,:,2)
532                  cwavef_y(2,:)=cwavef(2,:,1)-cwavef(1,:,2)
533                  call timab(538,1,tsec)
534                  if (nspinor1TreatedByThisProc) then
535                    call prep_fourwf(rhoaug_down,blocksize,cwavef(:,:,2),wfraug,&
536 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
537 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
538 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
539                  end if
540                  if (nspinor2TreatedByThisProc) then
541                    call prep_fourwf(rhoaug_mx,blocksize,cwavef_x,wfraug,&
542 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
543 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
544 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
545                    call prep_fourwf(rhoaug_my,blocksize,cwavef_y,wfraug,&
546 &                   iblock,istwf_k,dtset%mgfft,mpi_enreg,&
547 &                   nband_k,ndat,dtset%ngfft,npw_k,n4,n5,n6,occ_k,1,ucvol,&
548 &                   dtset%wtk(ikpt),use_gpu_cuda=dtset%use_gpu_cuda)
549                  end if
550                  call timab(538,2,tsec)
551                  ABI_DEALLOCATE(cwavef_x)
552                  ABI_DEALLOCATE(cwavef_y)
553                end if
554              end if
555            end do !iblock
556            if(ioption==1)  then
557              ABI_DEALLOCATE(kg_k_cart_block)
558            end if
559            if (allocated(cwavef))  then
560              ABI_DEALLOCATE(cwavef)
561            end if
562            ABI_DEALLOCATE(occ_k)
563          end if
564 
565          ABI_DEALLOCATE(gbound)
566          ABI_DEALLOCATE(kg_k)
567 
568          bdtot_index=bdtot_index+nband_k
569 
570          if (dtset%mkmem/=0) then
571            icg=icg+npw_k*my_nspinor*nband_k
572            ikg=ikg+npw_k
573          end if
574 
575 !        End loop on ikpt:
576        end do
577 
578 
579        if(mpi_enreg%paral_kgb == 1) then
580          call bandfft_kpt_set_ikpt(-1,mpi_enreg)
581          if (dtset%nspden==4) then
582 !          Sum the contribution of the band and of the FFT
583            call xmpi_sum(rhoaug     ,mpi_enreg%comm_bandspinorfft, ierr)
584            call xmpi_sum(rhoaug_down,mpi_enreg%comm_bandspinorfft, ierr)
585            call xmpi_sum(rhoaug_mx ,mpi_enreg%comm_bandspinorfft, ierr)
586            call xmpi_sum(rhoaug_my ,mpi_enreg%comm_bandspinorfft, ierr)
587            rhoaug_up(:,:,:) = rhoaug(:,:,:)
588          else
589            call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr)
590          end if
591        end if
592 
593 !      Transfer density on augmented fft grid to normal fft grid in real space
594 !      Take also into account the spin, to place it correctly in rhor.
595        if(dtset%nspden==1 .or. dtset%nspden==2) then
596          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug,1)
597        else if(dtset%nspden==4) then
598          ispden=1
599          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_up,1)
600          ispden=2
601          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_mx,1)
602          ispden=3
603          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_my,1)
604          ispden=4
605          call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug_down,1)
606          ABI_DEALLOCATE(rhoaug_up)
607          ABI_DEALLOCATE(rhoaug_down)
608          ABI_DEALLOCATE(rhoaug_mx)
609          ABI_DEALLOCATE(rhoaug_my)
610        end if
611 
612      end do !  isppol=1,dtset%nsppol
613 
614      if(allocated(cwavef))  then
615        ABI_DEALLOCATE(cwavef)
616      end if
617      ABI_DEALLOCATE(cwavefb)
618      ABI_DEALLOCATE(rhoaug)
619      ABI_DEALLOCATE(wfraug)
620 
621 !    Recreate full rhor on all proc.
622      call timab(48,1,tsec)
623      call timab(71,1,tsec)
624      spaceComm=mpi_enreg%comm_cell
625      if (mpi_enreg%paral_hf==1)spaceComm=mpi_enreg%comm_kpt
626      if(mpi_enreg%paral_kgb==1)spaceComm=mpi_enreg%comm_kpt
627      call xmpi_sum(rhor,spaceComm,ierr)
628      call timab(71,2,tsec)
629      call timab(48,2,tsec)
630 
631      call timab(799,2,tsec)
632      call timab(549,1,tsec)
633 
634      if(ioption==1 .or. ioption==2) then
635 !$OMP PARALLEL DO COLLAPSE(2)
636        do ispden=1,dtset%nspden
637          do ifft=1,dtset%nfft
638            taur_alphabeta(ifft,ispden,alpha,beta) = rhor(ifft,ispden)
639          end do
640        end do
641      end if
642 
643    end do !  beta=1,nbeta
644  end do !  alpha=1,nalpha
645 
646 !Compute the trace over the kinetic energy density tensor. i.e. Sum of the 3 diagonal elements.
647  if(ioption==1)then
648 !  zero rhor array in real space
649    do ispden=1,dtset%nspden
650 !$OMP PARALLEL DO
651      do ifft=1,dtset%nfft
652        rhor(ifft,ispden)=zero
653      end do
654    end do
655    do alpha = 1, nalpha
656 !$OMP PARALLEL DO COLLAPSE(2)
657      do ispden=1,dtset%nspden
658        do ifft=1,dtset%nfft
659          rhor(ifft,ispden) = rhor(ifft,ispden) + taur_alphabeta(ifft,ispden,alpha,1)
660        end do
661      end do
662    end do
663  end if
664 
665  nfftot=dtset%ngfft(1) * dtset%ngfft(2) * dtset%ngfft(3)
666 
667  select case (ioption)
668  case(0,1)
669    call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
670    dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
671    if(ioption==1)then
672 !$OMP PARALLEL DO
673      do ifft=1,dtset%nfft
674        do ispden=1,dtset%nspden
675          rhor(ifft,ispden)=1.0d0/2.0d0*rhor(ifft,ispden)
676        end do
677        rhog(:,ifft)=1.0d0/2.0d0*rhog(:,ifft)
678      end do
679    end if
680  case(2)
681    message = ' Sorry, kinetic energy density tensor (taur_(alpha,beta)) is not yet implemented.'
682    MSG_BUG(message)
683 
684 !    call symtaug(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,&
685 !    dtset%paral_kgb,phnons,rhog,rhor,rprimd,dtset%symafm,dtset%symrel)
686  end select
687 
688  call timab(549,2,tsec)
689 
690 !We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
691 !we also have the spin-up density, symmetrized, in rhor(:,2).
692 !In case of non collinear magnetism, we have rho,mx,my,mz. No symmetry is applied
693 
694  call timab(799,1,tsec)
695 
696  if(ioption==1 .or. ioption==2)  then
697    ABI_DEALLOCATE(taur_alphabeta)
698  end if
699 
700 !Find and print minimum and maximum total electron density
701 !(or total kinetic energy density, or total element of kinetic energy density tensor) and locations
702  call wrtout(std_out,'mkrho: echo density (plane-wave part only)','COLL')
703  call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,optrhor=ioption,ucvol=ucvol)
704 
705  call timab(799,2,tsec)
706  call timab(790+tim_mkrho,2,tsec)
707 
708  DBG_EXIT("COLL")
709 
710 end subroutine mkrho