TABLE OF CONTENTS


ABINIT/m_newrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_newrho

FUNCTION

COPYRIGHT

  Copyright (C) 2005-2022 ABINIT group (MT).
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_newrho
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_errors
27  use m_abicore
28  use m_ab7_mixing
29  use m_abi2big
30  use m_dtset
31 
32  use defs_datatypes, only : pseudopotential_type
33  use defs_abitypes,     only : MPI_type
34  use m_time,     only : timab
35  use m_geometry, only : metric
36  use m_pawtab,   only : pawtab_type
37  use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter
38  use m_prcref,   only : prcref
39  use m_wvl_rho, only : wvl_prcref
40  use m_fft,     only : fourdp
41 
42  implicit none
43 
44  private

ABINIT/newrho [ Functions ]

[ Top ] [ Functions ]

NAME

 newrho

FUNCTION

 Compute new trial density by mixing new and old values.
 Call prcref to compute preconditioned residual density and forces,
 Then, call one of the self-consistency drivers, then update density.

INPUTS

  atindx(natom)=index table for atoms (see gstate.f)
  dielar(7)=input parameters for dielectric matrix:
                diecut,dielng,diemac,diemix,diegap,dielam,diemixmag.
  dielinv(2,npwdiel,nspden,npwdiel,nspden)=
                              inverse of the dielectric matrix in rec. space
  dielstrt=number of the step at which the dielectric preconditioning begins.
  dtset <type(dataset_type)>=all input variables in this dataset
   | densfor_pred= governs the preconditioning of the atomic charges
   | iprcel= governs the preconditioning of the density residual
   | iprcfc= governs the preconditioning of the forces
   | iscf=( <= 0 =>non-SCF), >0 => SCF)
   |  iscf =11 => determination of the largest eigenvalue of the SCF cycle
   |  iscf =12 => SCF cycle, simple mixing
   |  iscf =13 => SCF cycle, Anderson mixing
   |  iscf =14 => SCF cycle, Anderson mixing (order 2)
   |  iscf =15 => SCF cycle, CG based on the minimization of the energy
   |  iscf =17 => SCF cycle, Pulay mixing
   | isecur=level of security of the computation
   | mffmem=governs the number of FFT arrays which are fit in core memory
   |          it is either 1, in which case the array f_fftgr is used,
   |          or 0, in which case the array f_fftgr_disk is used
   | natom=number of atoms
   | nspden=number of spin-density components
   | pawoptmix=-PAW- 1 if the computed residuals include the PAW (rhoij) part
   | prtvol=control print volume and debugging
  etotal=the total energy obtained from the input density
  fnametmp_fft=name of _FFT file
  fcart(3,natom)=cartesian forces (hartree/bohr)
  ffttomix(nfft*(1-nfftmix/nfft))=Index of the points of the FFT (fine) grid on the grid used for mixing (coarse)
  gmet(3,3)=metrix tensor in G space in Bohr**-2.
  grhf(3,natom)=Hellman-Feynman derivatives of the total energy
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  initialized= if 0, the initialization of the gstate run is not yet finished
  ispmix=1 if mixing is done in real space, 2 if mixing is done in reciprocal space
  istep= number of the step in the SCF cycle
  kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
  kxc(nfft,nkxc)=exchange-correlation kernel, needed only for electronic
     dielectric matrix
  mgfft=maximum size of 1D FFTs
  mixtofft(nfftmix*(1-nfftmix/nfft))=Index of the points of the FFT grid used for mixing (coarse) on the FFT (fine) grid
  moved_atm_inside= if 1, then the preconditioned forces
    as well as the preconditioned density residual must be computed;
    otherwise, compute only the preconditioned density residual.
  mpi_enreg=information about MPI parallelization
  my_natom=number of atoms treated by current processor
  nattyp(ntypat)=number of atoms of each type in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  nfftmix=dimension of FFT grid used to mix the densities (used in PAW only)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  ngfftmix(18)=contain all needed information about 3D FFT, for the grid corresponding to nfftmix
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  npawmix=-PAW only- number of spherical part elements to be mixed
  npwdiel=number of planewaves for dielectric matrix
  nresid(nfft,nspden)=array for the residual of the density
  ntypat=number of types of atoms in cell.
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
                                         Use here rhoij residuals (and gradients)
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  susmat(2,npwdiel,nspden,npwdiel,nspden)=
   the susceptibility (or density-density response) matrix in reciprocal space
  usepaw= 0 for non paw calculation; =1 for paw calculation
  vtrial(nfft,nspden)=the trial potential that gave vresid.
  xred(3,natom)=reduced dimensionless atomic coordinates
  tauresid(nfft,nspden*dtset%usekden)=array for kinetic energy density residue (out - in)

OUTPUT

  dbl_nnsclo=1 if nnsclo has to be doubled to secure the convergence.

SIDE EFFECTS

  dtn_pc(3,natom)=preconditioned change of atomic position,
                                          in reduced coordinates
  mix<type(ab7_mixing_object)>=all data defining the mixing algorithm for the density
  rhor(nfft,nspden)= at input, it is the "out" trial density that gave nresid=(rho_out-rho_in)
                     at output, it is an updated "mixed" trial density
  rhog(2,nfft)= Fourier transform of the new trial density
  ===== if usekden==1 =====
  [mix_mgga<type(ab7_mixing_object)>]=all data defining the mixing algorithm
     for the kinetic energy density
  ===== if densfor_pred==3 .and. moved_atm_inside==1 =====
    ph1d(2,3*(2*mgfft+1)*natom)=1-dim structure factor phases
  ==== if usepaw==1
    pawrhoij(natom)%nrhoijsel=number of non-zero values of rhoij
    pawrhoij(natom)%rhoijp(cplex_rhoij*lmn2_size,nspden)= new (mixed) value of rhoij quantities in PACKED STORAGE
    pawrhoij(natom)%rhoijselect(lmn2_size)=select the non-zero values of rhoij
  taug(2,nfft*dtset%usekden)=array for Fourier transform of kinetic
     energy density
  taur(nfft,nspden*dtset%usekden)=array for kinetic energy density

NOTES

  In case of PAW calculations:
    Computations are done either on the fine FFT grid or the coarse grid (depending on dtset%pawmixdg)
    All variables (nfft,ngfft,mgfft) refer to the fine FFT grid.
    All arrays (densities/potentials...) are computed on this fine FFT grid.
  ! Developpers have to be careful when introducing others arrays:
      they have to be stored on the fine FFT grid (except f_fftgr).
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

SOURCE

165 subroutine newrho(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,&
166 &  gmet,grhf,gsqcut,initialized,ispmix,istep,kg_diel,kxc,mgfft,mix,mixtofft,&
167 &  moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,&
168 &  nfftmix,nfftmix_per_nfft,ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
169 &  nresid,ntypat,n1xccc,pawrhoij,pawtab,&
170 &  ph1d,psps,rhog,rhor,rprimd,susmat,usepaw,vtrial,wvl,wvl_den,xred,&
171 &  mix_mgga,taug,taur,tauresid)
172 
173 !Arguments-------------------------------
174 !scalars
175  integer,intent(in) :: dielstrt,initialized,ispmix,istep,my_natom,mgfft
176  integer,intent(in) :: moved_atm_inside,n1xccc,nfft
177  integer,intent(in) :: nfftmix,nfftmix_per_nfft
178  integer,intent(in) :: nkxc,npawmix,npwdiel,ntypat,usepaw
179  integer,intent(inout) :: dbl_nnsclo
180  real(dp),intent(in) :: etotal,gsqcut
181  type(MPI_type),intent(in) :: mpi_enreg
182  type(ab7_mixing_object), intent(inout) :: mix
183  type(ab7_mixing_object), intent(inout),optional :: mix_mgga
184  type(dataset_type),intent(in) :: dtset
185  type(pseudopotential_type),intent(in) :: psps
186  type(wvl_internal_type), intent(in) :: wvl
187  type(wvl_denspot_type), intent(inout) :: wvl_den
188 !arrays
189  integer,intent(in) :: atindx(dtset%natom)
190  integer,intent(in) :: ffttomix(nfft*(nfftmix_per_nfft))
191  integer,intent(in) :: kg_diel(3,npwdiel)
192  integer,intent(in) :: mixtofft(nfftmix*nfftmix_per_nfft)
193  integer,intent(in) :: nattyp(ntypat),ngfft(18),ngfftmix(18)
194  real(dp),intent(in) :: dielar(7),fcart(3,dtset%natom),grhf(3,dtset%natom)
195  real(dp),intent(inout) :: rprimd(3,3)
196  real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
197  real(dp),intent(in), target :: vtrial(nfft,dtset%nspden)
198  real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
199  real(dp),intent(inout), target :: dtn_pc(3,dtset%natom)
200  real(dp),intent(inout) :: gmet(3,3)
201 !TODO: nresid appears to be only intent in here.
202  real(dp),intent(inout) :: kxc(nfft,nkxc),nresid(nfft,dtset%nspden)
203  real(dp),intent(inout) :: ph1d(2,3*(2*mgfft+1)*dtset%natom)
204  real(dp),intent(inout) :: rhor(nfft,dtset%nspden)
205  real(dp),intent(inout), target :: xred(3,dtset%natom)
206  real(dp),intent(inout) :: rhog(2,nfft)
207  real(dp),intent(inout), optional :: taug(2,nfft*dtset%usekden)
208  real(dp),intent(inout), optional :: taur(nfft,dtset%nspden*dtset%usekden)
209  real(dp),intent(inout), optional :: tauresid(nfft,dtset%nspden*dtset%usekden)
210 
211  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*psps%usepaw)
212  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
213 
214 !Local variables-------------------------------
215 !scalars
216  integer,parameter :: tim_fourdp9=9
217  integer :: cplex,dplex,errid,i_vresid1,i_vrespc1,iatom,ifft,indx,iq,iq0,irhoij,ispden,jfft
218  integer :: jrhoij,klmn,kklmn,kmix,mpicomm,nfftot,qphase
219  logical :: mpi_summarize,reset
220  real(dp) :: fact,ucvol,ucvol_local
221  character(len=500) :: message
222 !arrays
223  real(dp) :: gprimd(3,3),rmet(3,3),ro(2),tsec(2),vhartr_dum(1),vpsp_dum(1)
224  real(dp) :: vxc_dum(1,1)
225  real(dp),allocatable :: magng(:,:,:),magntaug(:,:,:)
226  real(dp),allocatable :: nresid0(:,:),nrespc(:,:),nreswk(:,:,:)
227  real(dp),allocatable :: rhoijrespc(:),rhoijtmp(:,:)
228 ! TODO : these should be allocatables not pointers: is there some reason to
229 !  keep them this way, eg an interface somewhere?
230  real(dp), pointer :: rhomag(:,:), npaw(:)
231  real(dp),allocatable :: tauresid0(:,:),taurespc(:,:)
232  real(dp),allocatable :: taumag(:,:)
233 
234 ! *************************************************************************
235 
236  DBG_ENTER("COLL")
237  call timab(94,1,tsec)
238 
239  nfftot=PRODUCT(ngfft(1:3))
240 
241 !Compatibility tests
242  if(nfftmix>nfft) then
243    message='nfftmix>nfft not allowed!'
244    ABI_BUG(message)
245  end if
246 
247  if(dtset%usewvl==1) then
248    if( (ispmix/=1 .or. nfftmix/=nfft)) then
249      message='nfftmix/=nfft, ispmix/=1 not allowed for wavelets!'
250      ABI_BUG(message)
251    end if
252    if(dtset%wvl_bigdft_comp==1) then
253      message='usewvl == 1 and wvl_bigdft_comp==1 not allowed!'
254      ABI_BUG(message)
255    end if
256  end if
257 
258  if(ispmix/=2.and.nfftmix/=nfft) then
259    message='nfftmix/=nfft allowed only when ispmix=2!'
260    ABI_BUG(message)
261  end if
262 
263  if (dtset%usekden==1) then
264    if ((.not.present(tauresid)).or.(.not.present(taug)).or. &
265 &      (.not.present(taur)).or.(.not.present(mix_mgga))) then
266      message='Several arrays are missing!'
267      ABI_BUG(message)
268    end if
269    if (mix_mgga%iscf==AB7_MIXING_CG_ENERGY.or.mix_mgga%iscf==AB7_MIXING_CG_ENERGY_2.or.&
270 &      mix_mgga%iscf==AB7_MIXING_EIG) then
271      message='kinetic energy density cannot be mixed with the selected mixing algorithm!'
272      ABI_ERROR(message)
273    end if
274  end if
275 
276  if (usepaw==1.and.my_natom>0) then
277    cplex=pawrhoij(1)%cplex_rhoij;dplex=cplex-1
278    qphase=pawrhoij(1)%qphase
279  else
280    cplex = 0;dplex = 0 ; qphase=0
281  end if
282 
283 !Compute different geometric tensor, as well as ucvol, from rprimd
284  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
285 
286  if(dtset%usewvl==0) then
287    ucvol_local=ucvol
288 #if defined HAVE_BIGDFT
289  else
290    ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(nfftot, dp)
291 #endif
292  end if
293 
294 !Select components of density to be mixed
295  ABI_MALLOC(rhomag,(ispmix*nfftmix,dtset%nspden))
296  ABI_MALLOC(nresid0,(ispmix*nfftmix,dtset%nspden))
297  ABI_MALLOC(taumag,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
298  ABI_MALLOC(tauresid0,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
299  ! real space and all fft points are here
300  if (ispmix==1.and.nfft==nfftmix) then
301    rhomag(:,1:dtset%nspden)=rhor(:,1:dtset%nspden)
302    nresid0(:,1:dtset%nspden)=nresid(:,1:dtset%nspden)
303    if (dtset%usekden==1) then
304      taumag(:,1:dtset%nspden)=taur(:,1:dtset%nspden)
305      tauresid0(:,1:dtset%nspden)=tauresid(:,1:dtset%nspden)
306    end if
307  ! recip space and all fft points are here
308  else if (nfft==nfftmix) then
309    do ispden=1,dtset%nspden
310      call fourdp(1,nresid0(:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
311    end do
312    rhomag(:,1)=reshape(rhog,(/2*nfft/))
313    if (dtset%nspden>1) then
314      do ispden=2,dtset%nspden
315        call fourdp(1,rhomag(:,ispden),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
316      end do
317    end if
318    if (dtset%usekden==1) then
319      do ispden=1,dtset%nspden
320        call fourdp(1,tauresid0(:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
321      end do
322      taumag(:,1)=reshape(taug,(/2*nfft/))
323      if (dtset%nspden>1) then
324        do ispden=2,dtset%nspden
325          call fourdp(1,taumag(:,ispden),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
326        end do
327      end if
328    end if
329  ! not all fft points are here - presumes recip space
330  else
331    fact=dielar(4)-1._dp
332    ABI_MALLOC(nreswk,(2,nfft,dtset%nspden))
333    do ispden=1,dtset%nspden
334      call fourdp(1,nreswk(:,:,ispden),nresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
335    end do
336    do ifft=1,nfft
337      if (ffttomix(ifft)>0) then
338        jfft=2*ffttomix(ifft)
339        rhomag (jfft-1:jfft,1)=rhog(1:2,ifft)
340        nresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1)
341      else
342        rhog(:,ifft)=rhog(:,ifft)+fact*nreswk(:,ifft,1)
343      end if
344    end do
345    if (dtset%nspden>1) then
346      ABI_MALLOC(magng,(2,nfft,dtset%nspden-1))
347      do ispden=2,dtset%nspden
348        call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
349        do ifft=1,nfft
350          if (ffttomix(ifft)>0) then
351            jfft=2*ffttomix(ifft)
352            rhomag (jfft-1:jfft,ispden)=magng (1:2,ifft,ispden-1)
353            nresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden)
354          else
355            magng(:,ifft,ispden-1)=magng(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden)
356            if (dtset%nspden==2) magng(:,ifft,1)=two*magng(:,ifft,1)-rhog(:,ifft)
357          end if
358        end do
359      end do
360    end if
361    if (dtset%usekden==1) then
362      do ispden=1,dtset%nspden
363        call fourdp(1,nreswk(:,:,ispden),tauresid(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
364      end do
365      do ifft=1,nfft
366        if (ffttomix(ifft)>0) then
367          jfft=2*ffttomix(ifft)
368          taumag (jfft-1:jfft,1)=taug(1:2,ifft)
369          tauresid0(jfft-1:jfft,1)=nreswk(1:2,ifft,1)
370        else
371          taug(:,ifft)=taug(:,ifft)+fact*nreswk(:,ifft,1)
372        end if
373      end do
374      if (dtset%nspden>1) then
375        ABI_MALLOC(magntaug,(2,nfft,dtset%nspden-1))
376        do ispden=2,dtset%nspden
377          call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
378          do ifft=1,nfft
379            if (ffttomix(ifft)>0) then
380              jfft=2*ffttomix(ifft)
381              taumag (jfft-1:jfft,ispden)=magntaug(1:2,ifft,ispden-1)
382              tauresid0(jfft-1:jfft,ispden)=nreswk(1:2,ifft,ispden)
383            else
384              magntaug(:,ifft,ispden-1)=magntaug(:,ifft,ispden-1)+fact*nreswk(:,ifft,ispden)
385              if (dtset%nspden==2) magntaug(:,ifft,1)=two*magntaug(:,ifft,1)-taug(:,ifft)
386            end if
387          end do
388        end do
389      end if
390    end if
391    ABI_FREE(nreswk)
392  end if
393 
394 !Retrieve "input" density from "output" density and density residual
395  rhomag(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)-nresid0(:,1:dtset%nspden)
396  if (dtset%usekden==1) then
397    taumag(:,1:dtset%nspden)=taumag(:,1:dtset%nspden)-tauresid0(:,1:dtset%nspden)
398  end if
399 
400 !If nspden==2, separate density and magnetization
401  if (dtset%nspden==2) then
402    rhomag (:,2)=two*rhomag (:,2)-rhomag (:,1)
403    nresid0(:,2)=two*nresid0(:,2)-nresid0(:,1)
404    if (dtset%usekden==1) then
405      taumag (:,2)=two*taumag (:,2)-taumag (:,1)
406      tauresid0(:,2)=two*tauresid0(:,2)-tauresid0(:,1)
407    end if
408  end if
409 
410 !If PAW, handle occupancy matrix
411  if (usepaw==1.and.my_natom>0) then
412    if (pawrhoij(1)%nspden==2) then
413      do iatom=1,my_natom
414        do iq=1,qphase
415          iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
416          jrhoij=1+iq0
417          do irhoij=1,pawrhoij(iatom)%nrhoijsel
418            ro(1:1+dplex)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
419            pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
420            pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
421            jrhoij=jrhoij+cplex
422          end do
423          do kmix=1,pawrhoij(iatom)%lmnmix_sz
424            klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
425            ro(1:1+dplex)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)
426            pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,1)=ro(1:1+dplex)+pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)
427            pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)=ro(1:1+dplex)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,2)
428          end do
429        end do
430      end do
431    end if
432  end if
433 
434 !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc
435  ABI_MALLOC(nrespc,(ispmix*nfftmix,dtset%nspden))
436  ABI_MALLOC(taurespc,(ispmix*nfftmix,dtset%nspden*dtset%usekden))
437  ABI_MALLOC(npaw,(npawmix*usepaw))
438  if (usepaw==1)  then
439    ABI_MALLOC(rhoijrespc,(npawmix))
440  else
441    ABI_MALLOC(rhoijrespc,(0))
442  end if
443  if(dtset%usewvl==0) then
444    call prcref(atindx,dielar,dielinv,&
445 &   dielstrt,dtn_pc,dtset,etotal,fcart,ffttomix,gmet,gsqcut,&
446 &   istep,kg_diel,kxc,&
447 &   mgfft,moved_atm_inside,mpi_enreg,my_natom,&
448 &   nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
449 &   ispmix,1,pawrhoij,pawtab,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,&
450 &   susmat,vhartr_dum,vpsp_dum,nresid0,nrespc,vxc_dum,wvl,wvl_den,xred)
451  else
452    call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,&
453 &   rhoijrespc,psps%usepaw,nresid0,nrespc)
454  end if
455 !At present, only a simple precoditionning for the kinetic energy density
456 ! (is Kerker mixing valid for tau?)
457  if (dtset%usekden==1) then
458    do ispden=1,dtset%nspden
459      fact=dielar(4);if (ispden>1) fact=abs(dielar(7))
460      taurespc(1:ispmix*nfftmix,ispden)=fact*tauresid0(1:ispmix*nfftmix,ispden)
461    end do
462  end if
463 
464 !------Compute new trial density and eventual new atomic positions
465 
466  if (mix%n_fftgr>0) then
467    i_vresid1=mix%i_vresid(1)
468    i_vrespc1=mix%i_vrespc(1)
469  end if
470 
471 !Initialise working arrays for the mixing object.
472  if (moved_atm_inside == 1) then
473    call ab7_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
474  end if
475  call ab7_mixing_eval_allocate(mix, istep)
476 
477 !Copy current step arrays.
478  if (moved_atm_inside == 1) then
479    call ab7_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc, arr_atm = grhf)
480  else
481    call ab7_mixing_copy_current_step(mix, nresid0, errid, message, arr_respc = nrespc)
482  end if
483  if (errid /= AB7_NO_ERROR) then
484    ABI_ERROR(message)
485  end if
486 
487 !Same treatment for the kinetic energy density
488  if (dtset%usekden==1) then
489    call ab7_mixing_eval_allocate(mix_mgga, istep)
490    call ab7_mixing_copy_current_step(mix_mgga, tauresid0, errid, message, arr_respc = taurespc)
491    if (errid /= AB7_NO_ERROR) then
492      ABI_ERROR(message)
493    end if
494  end if
495 
496  ABI_FREE(nresid0)
497  ABI_FREE(nrespc)
498  ABI_FREE(tauresid0)
499  ABI_FREE(taurespc)
500 
501 !PAW: either use the array f_paw or the array f_paw_disk
502  if (usepaw==1) then
503    indx=-dplex
504    do iatom=1,my_natom
505      ABI_MALLOC(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1))
506      do iq=1,qphase
507        iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
508        do ispden=1,pawrhoij(iatom)%nspden
509          rhoijtmp=zero ; jrhoij=1+iq0
510          do irhoij=1,pawrhoij(iatom)%nrhoijsel
511            klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
512            rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
513            jrhoij=jrhoij+cplex
514          end do
515          do kmix=1,pawrhoij(iatom)%lmnmix_sz
516            indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex ; kklmn=klmn+iq0
517            npaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
518            mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(kklmn:kklmn+dplex,ispden)
519            mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
520          end do
521        end do
522      end do
523      ABI_FREE(rhoijtmp)
524    end do
525  end if
526 
527 !------Prediction of the components of the density
528 
529 !Init mpicomm
530  if(mpi_enreg%paral_kgb==1)then
531    mpicomm=mpi_enreg%comm_fft
532    mpi_summarize=.true.
533  else
534    mpicomm=0
535    mpi_summarize=.false.
536  end if
537  if(dtset%usewvl==1) then
538    mpicomm=mpi_enreg%comm_wvl
539    mpi_summarize=(mpi_enreg%nproc_wvl > 1)
540  end if
541 
542  reset = .false.
543  if (initialized == 0) reset = .true.
544 
545 !Electronic density mixing
546  call ab7_mixing_eval(mix, rhomag, istep, nfftot, ucvol_local, &
547 & mpicomm, mpi_summarize, errid, message, &
548 & reset = reset, isecur = dtset%isecur,&
549 & pawopt = dtset%pawoptmix, pawarr = npaw, &
550 & etotal = etotal, potden = vtrial, &
551 & comm_atom=mpi_enreg%comm_atom)
552  if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
553    dbl_nnsclo = 1
554  else if (errid /= AB7_NO_ERROR) then
555    ABI_ERROR(message)
556  end if
557 !Kinetic energy density mixing (if any)
558  if (dtset%usekden==1) then
559    call ab7_mixing_eval(mix_mgga, taumag, istep, nfftot, ucvol_local, &
560 &   mpicomm, mpi_summarize, errid, message, reset = reset)
561    if (errid /= AB7_NO_ERROR) then
562      ABI_ERROR(message)
563    end if
564  end if
565 
566 !PAW: apply a simple mixing to rhoij (this is temporary)
567  if(dtset%iscf==15 .or. dtset%iscf==16)then
568    if (usepaw==1) then
569      indx=-dplex
570      do iatom=1,my_natom
571        ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
572        rhoijtmp=zero
573        do iq=1,qphase
574          iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
575          if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
576            do ispden=1,pawrhoij(iatom)%nspden
577              do kmix=1,pawrhoij(iatom)%lmnmix_sz
578                indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
579                rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) &
580 &               -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
581              end do
582            end do
583          end if
584          if (pawrhoij(iatom)%nspden/=2) then
585            do ispden=1,pawrhoij(iatom)%nspden
586              jrhoij=iq0+1
587              do irhoij=1,pawrhoij(iatom)%nrhoijsel
588                klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
589                rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
590 &               +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
591                jrhoij=jrhoij+cplex
592              end do
593            end do
594          else
595            jrhoij=iq0+1
596            do irhoij=1,pawrhoij(iatom)%nrhoijsel
597              klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
598              ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1)
599              rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2)) &
600 &             +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,1)
601              rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2)) &
602 &             +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,2)
603              jrhoij=jrhoij+cplex
604            end do
605          end if
606        end do
607        call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
608 &           pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
609 &           pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
610        ABI_FREE(rhoijtmp)
611      end do
612    end if
613  end if
614 
615  !if (usepaw==1)  then
616  ABI_FREE(rhoijrespc)
617  !end if
618 
619 !PAW: restore rhoij from compact storage
620  if (usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16) then
621    indx=-dplex
622    do iatom=1,my_natom
623      ABI_MALLOC(rhoijtmp,(cplex*qphase*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
624      rhoijtmp=zero
625      do iq=1,qphase
626        iq0=merge(0,cplex*pawrhoij(iatom)%lmn2_size,iq==1)
627        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
628          do ispden=1,pawrhoij(iatom)%nspden
629            jrhoij=iq0+1
630            do irhoij=1,pawrhoij(iatom)%nrhoijsel
631              klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
632              rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
633              jrhoij=jrhoij+cplex
634            end do
635          end do
636        end if
637        do ispden=1,pawrhoij(iatom)%nspden
638          do kmix=1,pawrhoij(iatom)%lmnmix_sz
639            indx=indx+cplex;klmn=iq0+cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
640            rhoijtmp(klmn:klmn+dplex,ispden)=npaw(indx:indx+dplex)
641          end do
642        end do
643        if (pawrhoij(iatom)%nspden==2) then
644          jrhoij=iq0+1
645          do irhoij=1,pawrhoij(iatom)%nrhoijsel
646            klmn=iq0+cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
647            ro(1:1+dplex)=rhoijtmp(klmn:klmn+dplex,1)
648            rhoijtmp(klmn:klmn+dplex,1)=half*(ro(1:1+dplex)+rhoijtmp(klmn:klmn+dplex,2))
649            rhoijtmp(klmn:klmn+dplex,2)=half*(ro(1:1+dplex)-rhoijtmp(klmn:klmn+dplex,2))
650            jrhoij=jrhoij+cplex
651          end do
652        end if
653      end do
654      call pawrhoij_filter(pawrhoij(iatom)%rhoijp,pawrhoij(iatom)%rhoijselect,&
655 &         pawrhoij(iatom)%nrhoijsel,pawrhoij(iatom)%cplex_rhoij,pawrhoij(iatom)%qphase,&
656 &         pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden,rhoij_input=rhoijtmp)
657      ABI_FREE(rhoijtmp)
658    end do
659  end if   ! usepaw==1.and.dtset%iscf/=15.and.dtset%iscf/=16
660  ABI_FREE(npaw)
661 
662 !Eventually write the data on disk and deallocate f_fftgr_disk
663  call ab7_mixing_eval_deallocate(mix)
664  if (dtset%usekden==1) call ab7_mixing_eval_deallocate(mix_mgga)
665 
666 !Fourier transform the density
667  if (ispmix==1.and.nfft==nfftmix) then
668    !Real space mixing, no need to transform rhomag
669    rhor(:,1:dtset%nspden)=rhomag(:,1:dtset%nspden)
670    if(dtset%usewvl==0) then
671      !Get rhog from rhor(:,1)
672      call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
673    end if
674    if (dtset%usekden==1) then
675      taur(:,1:dtset%nspden)=taumag(:,1:dtset%nspden)
676      if(dtset%usewvl==0) then
677        call fourdp(1,taug,taur(:,1),-1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
678      end if
679    end if
680  else if (nfft==nfftmix) then
681    !Reciprocal mixing space mixing, need to generate rhor in real space from rhomag in reciprocal space
682    do ispden=1,dtset%nspden
683      call fourdp(1,rhomag(:,ispden),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
684    end do
685    rhog(:,:)=reshape(rhomag(:,1),(/2,nfft/))
686    if (dtset%usekden==1) then
687      do ispden=1,dtset%nspden
688        call fourdp(1,taumag(:,ispden),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
689      end do
690      taug(:,:)=reshape(taumag(:,1),(/2,nfft/))
691    end if
692  else
693    do ifft=1,nfftmix
694      jfft=mixtofft(ifft)
695      rhog(1:2,jfft)=rhomag(2*ifft-1:2*ifft,1)
696    end do
697    call fourdp(1,rhog,rhor(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
698    if (dtset%nspden>1) then
699      do ispden=2,dtset%nspden
700        do ifft=1,nfftmix
701          jfft=mixtofft(ifft)
702          magng(1:2,jfft,ispden-1)=rhomag(2*ifft-1:2*ifft,ispden)
703        end do
704        call fourdp(1,magng(:,:,ispden-1),rhor(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
705      end do
706      ABI_FREE(magng)
707    end if
708    if (dtset%usekden==1) then
709      do ifft=1,nfftmix
710        jfft=mixtofft(ifft)
711        taug(1:2,jfft)=taumag(2*ifft-1:2*ifft,1)
712      end do
713      call fourdp(1,taug,taur(:,1),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
714      if (dtset%nspden>1) then
715        do ispden=2,dtset%nspden
716          do ifft=1,nfftmix
717            jfft=mixtofft(ifft)
718            magntaug(1:2,jfft,ispden-1)=taumag(2*ifft-1:2*ifft,ispden)
719          end do
720          call fourdp(1,magntaug(:,:,ispden-1),taur(:,ispden),+1,mpi_enreg,nfft,1,ngfft,tim_fourdp9)
721        end do
722        ABI_FREE(magntaug)
723      end if
724    end if
725  end if
726  ABI_FREE(rhomag)
727  ABI_FREE(taumag)
728 
729 !Set back rho in (up+dn,up) form if nspden=2
730  if (dtset%nspden==2) then
731    rhor(:,2)=half*(rhor(:,1)+rhor(:,2))
732    if (dtset%usekden==1) taur(:,2)=half*(taur(:,1)+taur(:,2))
733  end if
734 
735 !In WVL: copy density to BigDFT object:
736  if(dtset%usewvl==1) then
737    call wvl_rho_abi2big(1,rhor,wvl_den)
738  end if
739 
740  call timab(94,2,tsec)
741 
742  DBG_EXIT("COLL")
743 
744 end subroutine newrho