TABLE OF CONTENTS


ABINIT/m_newvtr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_newvtr

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_newvtr
23 
24  use defs_basis
25  use defs_wvltypes
26  use m_abicore
27  use m_errors
28  use m_abi2big
29  use m_ab7_mixing
30  use m_cgtools
31  use m_dtset
32 
33  use defs_datatypes, only : pseudopotential_type
34  use defs_abitypes,     only : MPI_type
35  use m_time,     only : timab
36  use m_geometry, only : metric
37  use m_pawtab,   only : pawtab_type
38  use m_pawrhoij, only : pawrhoij_type,pawrhoij_filter
39  use m_prcref,   only : prcref_PMA
40  use m_wvl_rho,  only : wvl_prcref
41  use m_fft,      only : fourdp
42  use m_xctk,     only : xcpot
43 
44  implicit none
45 
46  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
  [vtauresid(nfft,nspden*dtset%usekden)]=array for vtau residue (see vtau below))
  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)
  vres_mean(nspden)=mean value of the potential residual
  vpsp(nfft)=array for holding local psp
  vresid(nfft,nspden)=array for the residual of the potential
  vxc(nfft,nspden)=exchange-correlation potential (hartree)
  [vtau(nfftf,dtset%nspden,4*dtset%usekden)]=derivative of XC energy density
      with respect to kinetic energy density (metaGGA cases) (optional)
  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 usekden==1 =====
  [mix_mgga<type(ab7_mixing_object)>]=all data defining the mixing algorithm for
    the kinetic energy 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).

SOURCE

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