TABLE OF CONTENTS


ABINIT/m_newvtr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_newvtr

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, 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_newvtr
28 
29  use m_fft,     only : fourdp
30 
31  implicit none
32 
33  private

ABINIT/newvtr [ Functions ]

[ Top ] [ Functions ]

NAME

 newvtr

FUNCTION

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

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
   | spinmagntarget=input variable that governs fixed moment calculation
   | intxc=control xc quadrature
   | densfor_pred= governs the preconditioning of the atomic charges
   | iprcel= governs the preconditioning of the potential residual
   | iprcfc=governs the preconditioning of the forces
   | iscf=( <= 0 =>non-SCF), >0 => SCF)
   |  iscf =1 => determination of the largest eigenvalue of the SCF cycle
   |  iscf =2 => SCF cycle, simple mixing
   |  iscf =3 => SCF cycle, Anderson mixing
   |  iscf =4 => SCF cycle, Anderson mixing (order 2)
   |  iscf =5 => SCF cycle, CG based on the minimization of the energy
   |  iscf =7 => SCF cycle, Pulay mixing
   | isecur=level of security of the computation
   | ixc=exchange-correlation choice parameter.
   | 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
   | occopt=option for occupancies
   | paral_kgb=option for (kpt,g vectors,bands) parallelism
   | pawoptmix= - PAW only - 1 if the computed residuals include the PAW (rhoij) part
   | prtvol=control print volume and debugging
   | typat(natom)=integer type for each atom in cell
  etotal=the total energy obtained from the input vtrial
  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!
  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 potential residual must be computed;
    otherwise, compute only the preconditioned potential 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
  nstep=number of steps expected in iterations.
  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)
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
  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
  vhartr(nfft)=array for holding Hartree potential
  vnew_mean(nspden)=constrained mean value of the future trial potential (might be
    spin-polarized
  vpsp(nfft)=array for holding local psp
  vresid(nfft,nspden)=array for the residual of the potential
  vxc(nfft,nspden)=exchange-correlation potential (hartree)
  xred(3,natom)=reduced dimensionless atomic coordinates

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
  vtrial(nfft,nspden)= at input, it is the "in" trial potential that gave vresid=(v_out-v_in)
       at output, it is an updated "mixed" trial potential
  ===== 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,rhoijselect,rhoijp= several arrays
                containing new values of rhoij (augmentation occupancies)

WARNINGS

 depending on the value of densfor_pred and moved_atm_inside,
 the xc potential or the Hxc potential may have been subtracted from vtrial !

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.
  In case of norm-conserving calculations the FFT grid is the usual FFT grid.

  Subtility in PAW and non-collinear magnetism:
    Potentials are stored in (up-up,dn-dn,Re[up-dn],Im[up-dn]) format
    On-site occupancies (rhoij) are stored in (n,mx,my,mz)
    This is compatible provided that the mixing factors for n and m are identical
    and that the residual is not a combination of V_res and rhoij_res (pawoptmix=0).

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,mean_fftr
      metric,prcref_pma,timab,wvl_prcref,wvl_vtrial_abi2big

SOURCE

173 subroutine newvtr(atindx,dbl_nnsclo,dielar,dielinv,dielstrt,&
174      &  dtn_pc,dtset,etotal,fcart,ffttomix,&
175      &  gmet,grhf,gsqcut,&
176      &  initialized,ispmix,&
177      &  istep,&
178      &  kg_diel,kxc,mgfft,mix,mixtofft,&
179      &  moved_atm_inside,mpi_enreg,my_natom,nattyp,nfft,nfftmix,&
180      &  ngfft,ngfftmix,nkxc,npawmix,npwdiel,&
181      &  nstep,ntypat,n1xccc,&
182      &  pawrhoij,&
183      &  ph1d,&
184      &  psps,rhor,rprimd,susmat,usepaw,&
185      &  vhartr,vnew_mean,vpsp,vresid,&
186      &  vtrial,vxc,xred,&
187      &  nfftf,&
188      &  pawtab,&
189      &  rhog,&
190      &  wvl)
191 
192  use defs_basis
193  use defs_datatypes
194  use defs_abitypes
195  use defs_wvltypes
196  use m_abicore
197  use m_errors
198  use m_abi2big
199  use m_ab7_mixing
200  use m_cgtools
201 
202  use m_time,     only : timab
203  use m_geometry, only : metric
204  use m_pawtab,   only : pawtab_type
205  use m_pawrhoij, only : pawrhoij_type
206  use m_prcref,   only : prcref_PMA
207  use m_wvl_rho,  only : wvl_prcref
208 
209 !This section has been created automatically by the script Abilint (TD).
210 !Do not modify the following lines by hand.
211 #undef ABI_FUNC
212 #define ABI_FUNC 'newvtr'
213 !End of the abilint section
214 
215  implicit none
216 
217 !Arguments-------------------------------
218   ! WARNING
219   ! BEWARE THERE IS TWO DIFFERENT SIZE DECLARED FOR ARRAY NHAT IN RHOTOV AND RHOHXC
220   ! THIS MIGHT RESULT IN A BUG
221 !scalars
222  integer,intent(in) :: dielstrt,initialized,ispmix,istep,mgfft
223  integer,intent(in) :: moved_atm_inside,my_natom,n1xccc,nfft
224  integer,intent(in) :: nfftf,nfftmix,nkxc,npawmix,npwdiel,nstep
225  integer,intent(in) :: ntypat,usepaw
226  integer,intent(inout) :: dbl_nnsclo
227  real(dp),intent(in) :: etotal,gsqcut
228  type(MPI_type),intent(in) :: mpi_enreg
229  type(dataset_type),intent(in) :: dtset
230  type(ab7_mixing_object),intent(inout) :: mix
231  type(pseudopotential_type),intent(in) :: psps
232  type(wvl_data), intent(inout) :: wvl
233 !arrays
234  integer,intent(in) :: atindx(dtset%natom)
235  integer,intent(in) :: ffttomix(nfft*(1-nfftmix/nfft))
236  integer,intent(in) :: kg_diel(3,npwdiel)
237  integer,intent(in) :: mixtofft(nfftmix*(1-nfftmix/nfft)),nattyp(ntypat)
238  integer,intent(in) :: ngfft(18),ngfftmix(18)
239  real(dp),intent(in) :: dielar(7)
240  real(dp),intent(in) :: fcart(3,dtset%natom),grhf(3,dtset%natom)
241  real(dp),intent(in) :: rprimd(3,3)
242  real(dp),intent(in) :: susmat(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
243  real(dp),intent(in) :: vhartr(nfft),vnew_mean(dtset%nspden)
244  real(dp),intent(in) :: vxc(nfft,dtset%nspden)
245  real(dp),intent(inout) :: dielinv(2,npwdiel,dtset%nspden,npwdiel,dtset%nspden)
246  real(dp),intent(inout), target :: dtn_pc(3,dtset%natom)
247  real(dp),intent(inout) :: gmet(3,3)
248  real(dp),intent(inout) :: kxc(nfft,nkxc),ph1d(2,3*(2*mgfft+1)*dtset%natom)
249  real(dp),intent(inout) :: rhog(2,nfftf),vpsp(nfft)
250  real(dp),intent(inout), target :: rhor(nfft,dtset%nspden)
251  real(dp),intent(inout) :: vresid(nfft,dtset%nspden),vtrial(nfft,dtset%nspden)
252  real(dp),intent(inout), target :: xred(3,dtset%natom)
253  type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw)
254  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
255 
256 !Local variables-------------------------------
257 !scalars
258  integer :: cplex,dplex,i_vresid1,i_vrespc1
259 ! integer :: i1,i2,i3,ifft2,ifft3,ifft4,ifft5,ii1,ii2,ii3,ii4,ii5
260  integer :: iatom,ifft,indx,irhoij,mpicomm,mpi_comm_sphgrid
261  integer :: ispden,jfft,jrhoij,klmn,kmix,n1,n2,n3,nfftot,nselect
262  integer :: errid,tim_fourdp
263  logical :: mpi_summarize,reset
264  real(dp) :: dielng,diemix,fact,ucvol,ucvol_local,vme
265  character(len=500) :: message
266 !arrays
267  real(dp),parameter :: identity(4)=(/one,one,zero,zero/)
268  real(dp) :: gprimd(3,3),rmet(3,3),tsec(2),vmean(dtset%nspden)
269  real(dp),allocatable :: rhoijrespc(:)
270  real(dp),allocatable :: rhoijtmp(:,:)
271  real(dp),allocatable :: vresid0(:,:),vrespc(:,:),vreswk(:,:),vtrialg(:,:,:)
272  real(dp),pointer :: vtrial0(:,:),vpaw(:)
273 
274 ! *************************************************************************
275 
276 !DEBUG
277 !write(std_out,*)' newvtr : enter '
278 !write(std_out,*)' newvtr : ispmix,nfft,nfftmix=',ispmix,nfft,nfftmix
279 !ENDDEBUG
280 
281  call timab(93,1,tsec)
282  call timab(901,1,tsec)
283  tim_fourdp=8
284 
285 !mpicomm over spherical grid:
286  mpi_comm_sphgrid=mpi_enreg%comm_fft
287  if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
288 
289 !Compatibility tests
290  if(nfftmix>nfft) then
291    MSG_BUG('  nfftmix>nfft not allowed !')
292  end if
293 
294  if(ispmix/=2.and.nfftmix/=nfft) then
295    message = '  nfftmix/=nfft allowed only when ispmix=2 !'
296    MSG_BUG(message)
297  end if
298 
299  if(dtset%usewvl==1) then
300    if(dtset%wvl_bigdft_comp==1) then
301      message = 'newvtr: usewvl == 1 and wvl_bigdft_comp==1 not allowed (use wvl_newtr() instead)!'
302      MSG_BUG(message)
303    end if
304    if(ispmix/=1 .or. nfftmix/=nfft) then
305      MSG_BUG('newvtr: nfftmix/=nfft, ispmix/=1 not allowed for wavelets')
306    end if
307  end if
308 
309  if(usepaw==1.and.dtset%nspden==4.and.dtset%pawoptmix==1) then
310    message = ' pawoptmix=1 is not compatible with nspden=4 !'
311    MSG_ERROR(message)
312  end if
313 
314  dielng=dielar(2)
315  diemix=dielar(4)
316  n1=ngfft(1)
317  n2=ngfft(2)
318  n3=ngfft(3)
319  if (usepaw==1.and.my_natom>0) then
320    cplex=pawrhoij(1)%cplex;dplex=cplex-1
321  else
322    cplex=0;dplex=0
323  end if
324 
325 !Get size of FFT grid
326  nfftot=PRODUCT(ngfft(1:3))
327 
328 !Compute different geometric tensor, as well as ucvol, from rprimd
329  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
330 
331  if(dtset%usewvl==0) then
332    ucvol_local=ucvol
333 #if defined HAVE_BIGDFT
334  else
335    ucvol_local = product(wvl%den%denspot%dpbox%hgrids) * real(nfftot, dp)
336 #endif
337  end if
338 
339 !------Treat the mean of potentiel residual
340 
341 !Special care must be taken with components of the
342 !potential that are associated with NO density change.
343 !In general, only the global mean of the potential has
344 !such an anomalous feature. However, in the spin
345 !polarized cas with fixed occupancies, also the
346 !mean of each spin-potential (independently of the other)
347 !has such a behaviour. The trick is to remove these
348 !variables before going in the predictive routines,
349 !then to put them back
350 
351 !Compute the mean of the old vtrial
352  call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid)
353 
354 !When (collinear) spin-polarized and fixed occupation numbers,
355 !treat separately spin up and spin down.
356 !Otherwise, use only global mean
357  do ispden=1,dtset%nspden
358    if (dtset%nspden==2.and.dtset%occopt>=3.and. &
359 &   abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
360      vme=(vmean(1)+vmean(2))*half
361    else
362      vme=vmean(ispden)
363    end if
364    vtrial(:,ispden)=vtrial(:,ispden)-vme
365  end do
366 
367  call timab(901,2,tsec)
368 
369 !Select components of potential to be mixed
370  ABI_ALLOCATE(vtrial0,(ispmix*nfftmix,dtset%nspden))
371  ABI_ALLOCATE(vresid0,(ispmix*nfftmix,dtset%nspden))
372  if (ispmix==1.and.nfft==nfftmix) then
373    vtrial0=vtrial;vresid0=vresid
374  else if (nfft==nfftmix) then
375    do ispden=1,dtset%nspden
376      call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
377      call fourdp(1,vresid0(:,ispden),vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
378    end do
379  else
380    ABI_ALLOCATE(vtrialg,(2,nfft,dtset%nspden))
381    ABI_ALLOCATE(vreswk,(2,nfft))
382    do ispden=1,dtset%nspden
383      fact=dielar(4);if (ispden>1) fact=dielar(7)
384      call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
385      call fourdp(1,vreswk,vresid(:,ispden),-1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
386      do ifft=1,nfft
387        if (ffttomix(ifft)>0) then
388          jfft=2*ffttomix(ifft)
389          vtrial0(jfft-1,ispden)=vtrialg(1,ifft,ispden)
390          vtrial0(jfft  ,ispden)=vtrialg(2,ifft,ispden)
391          vresid0(jfft-1,ispden)=vreswk(1,ifft)
392          vresid0(jfft  ,ispden)=vreswk(2,ifft)
393        else
394          vtrialg(:,ifft,ispden)=vtrialg(:,ifft,ispden)+fact*vreswk(:,ifft)
395        end if
396      end do
397    end do
398    ABI_DEALLOCATE(vreswk)
399  end if
400 
401  call timab(902,1,tsec)
402 
403 !Choice of preconditioner governed by iprcel, densfor_pred and iprcfc
404  ABI_ALLOCATE(vrespc,(ispmix*nfftmix,dtset%nspden))
405  ABI_ALLOCATE(vpaw,(npawmix*usepaw))
406  if (usepaw==1)  then
407    ABI_ALLOCATE(rhoijrespc,(npawmix))
408  else
409    ABI_ALLOCATE(rhoijrespc,(0))
410  end if
411 
412  call timab(902,2,tsec)
413  call timab(903,1,tsec)
414 
415  if(dtset%usewvl==0) then
416    call prcref_PMA(atindx,dielar,dielinv,dielstrt,dtn_pc,dtset,fcart,ffttomix,gmet,gsqcut,&
417 &   istep,kg_diel,kxc,mgfft,moved_atm_inside,mpi_enreg,my_natom,&
418 &   nattyp,nfft,nfftmix,ngfft,ngfftmix,nkxc,npawmix,npwdiel,ntypat,n1xccc,&
419 &   ispmix,0,pawrhoij,ph1d,psps,rhog,rhoijrespc,rhor,rprimd,susmat,&
420 &   vhartr,vpsp,vresid0,vrespc,vxc,xred,&
421 &   etotal,pawtab,wvl)
422  else
423    call wvl_prcref(dielar,dtset%iprcel,my_natom,nfftmix,npawmix,dtset%nspden,pawrhoij,&
424 &   rhoijrespc,psps%usepaw,vresid0,vrespc)
425  end if
426 
427  call timab(903,2,tsec)
428  call timab(904,1,tsec)
429 
430 !------Compute new vtrial and eventual new atomic positions
431 
432  i_vresid1=mix%i_vresid(1)
433  i_vrespc1=mix%i_vrespc(1)
434 
435 !Initialise working arrays for the mixing object.
436  if (moved_atm_inside == 1) then
437    call ab7_mixing_use_moving_atoms(mix, dtset%natom, xred, dtn_pc)
438  end if
439  call ab7_mixing_eval_allocate(mix, istep)
440 !Copy current step arrays.
441  if (moved_atm_inside == 1) then
442    call ab7_mixing_copy_current_step(mix, vresid0, errid, message, &
443 &   arr_respc = vrespc, arr_atm = grhf)
444  else
445    call ab7_mixing_copy_current_step(mix, vresid0, errid, message, &
446 &   arr_respc = vrespc)
447  end if
448  if (errid /= AB7_NO_ERROR) then
449    MSG_ERROR(message)
450  end if
451  ABI_DEALLOCATE(vresid0)
452 
453 !PAW: either use the array f_paw or the array f_paw_disk
454  if (usepaw==1) then
455    indx=-dplex
456    do iatom=1,my_natom
457      do ispden=1,pawrhoij(iatom)%nspden
458        ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,1))
459        rhoijtmp=zero
460        jrhoij=1
461        do irhoij=1,pawrhoij(iatom)%nrhoijsel
462          klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
463          rhoijtmp(klmn:klmn+dplex,1)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
464          jrhoij=jrhoij+cplex
465        end do
466        do kmix=1,pawrhoij(iatom)%lmnmix_sz
467          indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
468          vpaw(indx:indx+dplex)=rhoijtmp(klmn:klmn+dplex,1)-pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
469          mix%f_paw(indx:indx+dplex,i_vresid1)=pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
470          mix%f_paw(indx:indx+dplex,i_vrespc1)=rhoijrespc(indx:indx+dplex)
471        end do
472        ABI_DEALLOCATE(rhoijtmp)
473      end do
474    end do
475  end if
476 
477 !------Prediction of the components of the potential associated with a density change
478 
479 !Init mpicomm
480  if(mpi_enreg%paral_kgb==1)then
481    mpicomm=mpi_enreg%comm_fft
482    mpi_summarize=.true.
483  else
484    mpicomm=0
485    mpi_summarize=.false.
486  end if
487 
488  reset = .false.
489  if (initialized == 0) reset = .true.
490  call ab7_mixing_eval(mix, vtrial0, istep, nfftot, ucvol_local, &
491 & mpicomm, mpi_summarize, errid, message, &
492 & reset = reset, isecur = dtset%isecur, &
493 & pawopt = dtset%pawoptmix, pawarr = vpaw, etotal = etotal, potden = rhor, &
494 & comm_atom=mpi_enreg%comm_atom)
495 
496  if (errid == AB7_ERROR_MIXING_INC_NNSLOOP) then
497    dbl_nnsclo = 1
498  else if (errid /= AB7_NO_ERROR) then
499    MSG_ERROR(message)
500  end if
501 
502 !PAW: apply a simple mixing to rhoij (this is temporary)
503  if(dtset%iscf==5 .or. dtset%iscf==6)then
504    if (usepaw==1) then
505      indx=-dplex
506      do iatom=1,my_natom
507        ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
508        rhoijtmp=zero
509        if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
510          do ispden=1,pawrhoij(iatom)%nspden
511            do kmix=1,pawrhoij(iatom)%lmnmix_sz
512              indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
513              rhoijtmp(klmn:klmn+dplex,ispden)=rhoijrespc(indx:indx+dplex) &
514 &             -pawrhoij(iatom)%rhoijres(klmn:klmn+dplex,ispden)
515            end do
516          end do
517        end if
518        do ispden=1,pawrhoij(iatom)%nspden
519          jrhoij=1
520          do irhoij=1,pawrhoij(iatom)%nrhoijsel
521            klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
522            rhoijtmp(klmn:klmn+dplex,ispden)=rhoijtmp(klmn:klmn+dplex,ispden) &
523 &           +pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
524            jrhoij=jrhoij+cplex
525          end do
526        end do
527        nselect=0
528        do klmn=1,pawrhoij(iatom)%lmn2_size
529          if (any(abs(rhoijtmp(cplex*klmn-dplex:cplex*klmn,:))>tol10)) then
530            nselect=nselect+1
531            pawrhoij(iatom)%rhoijselect(nselect)=klmn
532            do ispden=1,pawrhoij(iatom)%nspden
533              pawrhoij(iatom)%rhoijp(cplex*nselect-dplex:cplex*nselect,ispden)=&
534 &             rhoijtmp(cplex*klmn-dplex:cplex*klmn,ispden)
535            end do
536          end if
537        end do
538        pawrhoij(iatom)%nrhoijsel=nselect
539        ABI_DEALLOCATE(rhoijtmp)
540      end do
541    end if
542  end if
543 
544  ABI_DEALLOCATE(rhoijrespc)
545 
546 !PAW: restore rhoij from compact storage
547  if (usepaw==1.and.dtset%iscf/=5.and.dtset%iscf/=6) then
548    indx=-dplex
549    do iatom=1,my_natom
550      ABI_ALLOCATE(rhoijtmp,(cplex*pawrhoij(iatom)%lmn2_size,pawrhoij(iatom)%nspden))
551      rhoijtmp=zero
552      if (pawrhoij(iatom)%lmnmix_sz<pawrhoij(iatom)%lmn2_size) then
553        do ispden=1,pawrhoij(iatom)%nspden
554          jrhoij=1
555          do irhoij=1,pawrhoij(iatom)%nrhoijsel
556            klmn=cplex*pawrhoij(iatom)%rhoijselect(irhoij)-dplex
557            rhoijtmp(klmn:klmn+dplex,ispden)=pawrhoij(iatom)%rhoijp(jrhoij:jrhoij+dplex,ispden)
558            jrhoij=jrhoij+cplex
559          end do
560        end do
561      end if
562      do ispden=1,pawrhoij(iatom)%nspden
563        do kmix=1,pawrhoij(iatom)%lmnmix_sz
564          indx=indx+cplex;klmn=cplex*pawrhoij(iatom)%kpawmix(kmix)-dplex
565          rhoijtmp(klmn:klmn+dplex,ispden)=vpaw(indx:indx+dplex)
566        end do
567      end do
568      nselect=0
569      if (cplex==1) then
570        do klmn=1,pawrhoij(iatom)%lmn2_size
571          if (any(abs(rhoijtmp(klmn,:))>tol10)) then
572            nselect=nselect+1
573            pawrhoij(iatom)%rhoijselect(nselect)=klmn
574            do ispden=1,pawrhoij(iatom)%nspden
575              pawrhoij(iatom)%rhoijp(nselect,ispden)=rhoijtmp(klmn,ispden)
576            end do
577          end if
578        end do
579      else
580        do klmn=1,pawrhoij(iatom)%lmn2_size
581          if (any(abs(rhoijtmp(2*klmn-1:2*klmn,:))>tol10)) then
582            nselect=nselect+1
583            pawrhoij(iatom)%rhoijselect(nselect)=klmn
584            do ispden=1,pawrhoij(iatom)%nspden
585              pawrhoij(iatom)%rhoijp(2*nselect-1:2*nselect,ispden)=rhoijtmp(2*klmn-1:2*klmn,ispden)
586            end do
587          end if
588        end do
589      end if
590      pawrhoij(iatom)%nrhoijsel=nselect
591      ABI_DEALLOCATE(rhoijtmp)
592    end do
593  end if
594  ABI_DEALLOCATE(vpaw)
595  ABI_DEALLOCATE(vrespc)
596 
597 !Eventually write the data on disk and deallocate f_fftgr_disk
598  call ab7_mixing_eval_deallocate(mix)
599 
600  call timab(904,2,tsec)
601 
602 !Restore potential
603  if (ispmix==1.and.nfft==nfftmix) then
604    vtrial=vtrial0
605  else if (nfft==nfftmix) then
606    do ispden=1,dtset%nspden
607      call fourdp(1,vtrial0(:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
608    end do
609  else
610    do ispden=1,dtset%nspden
611      do ifft=1,nfftmix
612        jfft=mixtofft(ifft)
613        vtrialg(1,jfft,ispden)=vtrial0(2*ifft-1,ispden)
614        vtrialg(2,jfft,ispden)=vtrial0(2*ifft  ,ispden)
615      end do
616      call fourdp(1,vtrialg(:,:,ispden),vtrial(:,ispden),+1,mpi_enreg,nfft,ngfft,dtset%paral_kgb,tim_fourdp)
617    end do
618    ABI_DEALLOCATE(vtrialg)
619  end if
620  ABI_DEALLOCATE(vtrial0)
621 
622  call timab(905,1,tsec)
623 
624 !------Treat the mean of the potential
625 
626 !Compute the mean of the new vtrial
627  call mean_fftr(vtrial,vmean,nfft,nfftot,dtset%nspden,mpi_comm_sphgrid)
628 
629 !Reset the mean of the new vtrial, to the value vnew_mean
630 !When spin-polarized and fixed occupation numbers,
631 !treat separately spin up and spin down.
632 !Otherwise, use only global mean
633  do ispden=1,dtset%nspden
634    if (dtset%nspden==2.and.dtset%occopt>=3.and. &
635 &   abs(dtset%spinmagntarget+99.99_dp)<1.0d-10)then
636      vme=(vnew_mean(1)+vnew_mean(2)-vmean(1)-vmean(2))*half
637    else
638      vme=vnew_mean(ispden)-vmean(ispden)
639    end if
640    vtrial(:,ispden)=vtrial(:,ispden)+vme
641  end do
642 
643  if(moved_atm_inside==1 .and. istep/=nstep )then
644    if(abs(dtset%densfor_pred)==1.or.abs(dtset%densfor_pred)==4)then
645 !    Subtract current local psp, but also vxc (for core charges)
646      do ispden=1,dtset%nspden
647        vtrial(:,ispden)=vtrial(:,ispden)-vpsp(:)*identity(ispden)-vxc(:,ispden)
648      end do
649    else if(abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)then
650 !    Subtract current vpsp+Hxc from vtrial. This should be rationalized later
651      do ispden=1,dtset%nspden
652        vtrial(:,ispden)=vtrial(:,ispden)-(vpsp(:)+vhartr(:))*identity(ispden)-vxc(:,ispden)
653      end do
654    end if
655  end if
656 
657 !In WVL: copy vtrial to BigDFT object:
658  if(dtset%usewvl==1) then
659    call wvl_vtrial_abi2big(1,vtrial,wvl%den)
660  end if
661 
662  call timab(905,2,tsec)
663  call timab(93,2,tsec)
664 
665 end subroutine newvtr