TABLE OF CONTENTS


ABINIT/m_newrho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_newrho

FUNCTION

COPYRIGHT

  Copyright (C) 2005-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_newrho
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32 
33  use m_wvl_rho, only : wvl_prcref
34  use m_fft,     only : fourdp
35 
36  implicit none
37 
38  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
  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 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*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.

PARENTS

      scfcv

CHILDREN

      ab7_mixing_copy_current_step,ab7_mixing_eval,ab7_mixing_eval_allocate
      ab7_mixing_eval_deallocate,ab7_mixing_use_moving_atoms,fourdp,metric
      prcref,timab,wvl_prcref,wvl_rho_abi2big

SOURCE

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