TABLE OF CONTENTS


ABINIT/dfpt_vtorho [ Functions ]

[ Top ] [ Functions ]

NAME

 dfpt_vtorho

FUNCTION

 This routine compute the new 1-density from a fixed 1-potential (vtrial1)
 but might also simply compute eigenvectors and eigenvalues.
 The main part of it is a wf update over all k points

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, DRH, MB, XW, 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 .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions
  cgq(2,mpw1*nspinor*mband*mkqmem*nsppol)=pw coefficients of GS wavefunctions at k+q.
  cg1(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF wavefunctions at k,q.
  cplex: if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX
  cprj(natom,nspinor*mband*mkmem*nsppol*usecprj)= wave functions at k
              projected with non-local projectors: cprj=<p_i|Cnk>
  cprjq(natom,nspinor*mband*mkqmem*nsppol*usecprj)= wave functions at k+q
              projected with non-local projectors: cprjq=<p_i|Cnk+q>
  dbl_nnsclo=if 1, will double the value of dtset%nnsclo
  doccde_rbz(mband*nkpt_rbz*nsppol)=derivative of occ_rbz wrt the energy
  docckqde(mband*nkpt_rbz*nsppol)=derivative of occkq wrt the energy
  dtefield = variables related to response Berry-phase calculation
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
  eigenq(mband*nkpt_rbz*nsppol)=GS eigenvalues at k+q (hartree)
  eigen0(mband*nkpt_rbz*nsppol)=GS eigenvalues at k (hartree)
  fermie1=derivative of fermi energy wrt (strain) perturbation
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
  idir=direction of the perturbation
  indsy1(4,nsym1,natom)=indirect indexing array for atom labels
  ipert=type of the perturbation
  irrzon1(nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  istwfk_rbz(nkpt_rbz)=input option parameter that describes the storage of wfs
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg1(3,mpw1*mk1mem)=reduced planewave coordinates at k+q, with RF k points
  kpt_rbz(3,nkpt_rbz)=reduced coordinates of k points.
  mband=maximum number of bands
  mkmem =number of k points treated by this node (GS data).
  mkqmem =number of k+q points treated by this node (GS data)
  mk1mem =number of k points treated by this node (RF data)
  mpw=maximum dimensioned size of npw or wfs at k
  mpw1=maximum dimensioned size of npw for wfs at k+q (also for 1-order wfs).
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nband_rbz(nkpt_rbz*nsppol)=number of bands at each RF k point for each spin
  ncpgr=number of gradients stored in cprj array (cprj=<p_i|Cnk>)
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs - see comment in respfn.F90)
  nkpt_rbz=number of k points in the IBZ for this perturbation
  mpi_enreg=information about MPI parallelization
  npwarr(nkpt_rbz)=number of planewaves in basis at this GS k point
  npwar1(nkpt_rbz)=number of planewaves in basis at this RF k+q point
  nspden=number of spin-density components
  nsppol=1 for unpolarized, 2 for spin-polarized
  nsym1=number of symmetry elements in space group consistent with
   perturbation
  ntypat=number of types of atoms in unit cell.
  occkq(mband*nkpt_rbz*nsppol)=occupation number for each band (often 2)
   at each k+q point of the reduced Brillouin zone.
  occ_rbz(mband*nkpt_rbz*nsppol)=occupation number for each band and k
   (usually 2)
  optres=0: the new value of the density is computed in place of the input value
         1: only the density residual is computed ; the input density is kept
  paw_ij(natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS
  paw_ij1(natom) <type(paw_ij_type)>=1st-order paw arrays given on (i,j) channels
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawang1 <type(pawang_type)>=pawang datastructure containing only the symmetries preserving the perturbation
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid for the GS
  pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
  pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data for the GS
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
  ph1d(2,3*(2*dtset%mgfft+1)*natom)=one-dimensional structure factor information
  prtvol=control print volume and debugging output
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices
  qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
  inverse of the overlap matrix
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional real space primitive translations
  symaf1(nsym1)=(anti)ferromagnetic part of symmetry operations
  symrc1(3,3,nsym1)=symmetry operations in reciprocal space
  symrl1(3,3,nsym1)=symmetry operations in real space
  ucvol=unit cell volume in bohr**3.
  usecprj= 1 if cprj, cprjq, cprj1 arrays are stored in memory
  useylmgr1= 1 if ylmgr1 array is allocated
  ddk<wfk_t)=struct info DDK file
  vtrial(nfftf,nspden)=GS Vtrial(r).
  vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r).
  wtk_rbz(nkpt_rbz)=weight assigned to each k point.
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylm1(mpw1*mk1mem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k+g point
  ylmgr1(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics for each G and k+g point

OUTPUT

  cg1(2,mpw*nspinor*mband*mk1mem*nsppol)=updated wavefunctions, orthogonalized to the occupied states
  cg1_active(2,mpw1*nspinor*mband*mk1mem*nsppol)=pw coefficients of RF
    wavefunctions at k,q. They are orthogonalized to the active.
  eigen1(2*mband*mband*nkpt_rbz*nsppol)=array for holding eigenvalues
    (hartree)
  edocc=correction to 2nd-order total energy coming from changes of occupation
  eeig0=0th-order eigenenergies part of 2nd-order total energy
  ek0=0th-order kinetic energy part of 2nd-order total energy.
  ek1=1st-order kinetic energy part of 2nd-order total energy
    (not for phonons)
  eloc0=0th-order local (psp+vxc+Hart) part of 2nd-order total energy
  enl0=0th-order nonlocal pseudopot. part of 2nd-order total energy.
  enl1=1st-order nonlocal pseudopot. part of 2nd-order total energy.
  gh1c_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(1)}|nK>
  gh0c1_set(2,mpw1*nspinor*mband*mk1mem*nsppol*dim_eig2rf)= set of <G|H^{(0)}|\Psi^{(1)}>
      The wavefunction is orthogonal to the active space (for metals). It is not
      coherent with cg1.
  resid(mband*nkpt_rbz*nsppol)=residuals for each band over all k points.
  residm=maximum value from resid array (except for nbdbuf highest bands)
  rhog1(2,nfftf)=RF electron density in reciprocal space
  ==== if optres==1
    nres2=square of the norm of the residual
    nvresid1(cplex*nfftf,nspden)=1st-order density residual
  ==== if psps%usepaw==1
    cprj1(natom,nspinor*mband*mk1mem*nsppol*usecprj)=
              1st-order wave functions at k,q projected with non-local projectors:
                       cprj1=<p_i|C1nk,q> where p_i is a non-local projector
    nhat1(cplex*nfftf,nspden*psps%usepaw)=1st-order compensation charge density

SIDE EFFECTS

  pawrhoij1(natom) <type(pawrhoij_type)>= 1st-order paw rhoij occupancies and related data
  rhor1(cplex*nfftf,nspden)=RF electron density in electrons/bohr**3.

PARENTS

      dfpt_scfcv

CHILDREN

      destroy_hamiltonian,destroy_rf_hamiltonian,dfpt_vtowfk,dfptff_gbefd
      dfptff_gradberry,fftpac,getgh1c_setup,init_hamiltonian
      init_rf_hamiltonian,load_spin_hamiltonian,load_spin_rf_hamiltonian
      occeig,pawmkrho,pawrhoij_alloc,pawrhoij_free,pawrhoij_init_unpacked
      pawrhoij_mpisum_unpacked,rf_transgrid_and_pack,sqnorm_v,symrhg,timab
      xmpi_sum

SOURCE

154 #if defined HAVE_CONFIG_H
155 #include "config.h"
156 #endif
157 
158 #include "abi_common.h"
159 
160 subroutine dfpt_vtorho(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dbl_nnsclo,&
161 & dim_eig2rf,doccde_rbz,docckqde,dtefield,dtfil,dtset,qphon,&
162 & edocc,eeig0,eigenq,eigen0,eigen1,ek0,ek1,eloc0,enl0,enl1,&
163 & fermie1,gh0c1_set,gh1c_set,gmet,gprimd,idir,indsy1,&
164 & ipert,irrzon1,istwfk_rbz,kg,kg1,kpt_rbz,mband,&
165 & mkmem,mkqmem,mk1mem,mpi_enreg,mpw,mpw1,my_natom,&
166 & natom,nband_rbz,ncpgr,nfftf,nhat1,nkpt_rbz,npwarr,npwar1,nres2,nspden,&
167 & nsppol,nsym1,ntypat,nvresid1,occkq,occ_rbz,optres,&
168 & paw_ij,paw_ij1,pawang,pawang1,pawfgr,pawfgrtab,pawrhoij,pawrhoij1,pawtab,&
169 & phnons1,ph1d,prtvol,psps,pwindall,qmat,resid,residm,rhog1,rhor1,rmet,rprimd,symaf1,symrc1,symrl1,ucvol,&
170 & usecprj,useylmgr1,ddk_f,vtrial,vtrial1,wtk_rbz,xred,ylm,ylm1,ylmgr1,cg1_out)
171 
172  use defs_basis
173  use defs_datatypes
174  use defs_abitypes
175  use m_profiling_abi
176  use m_xmpi
177  use m_errors
178  use m_efield
179  use m_hamiltonian
180  use m_wfk
181  use m_cgtools
182 
183  use m_occ,      only : occeig
184  use m_hdr,      only : hdr_skip, hdr_io
185  use m_pawang,   only : pawang_type
186  use m_pawtab,   only : pawtab_type
187  use m_paw_ij,   only : paw_ij_type
188  use m_pawfgrtab,only : pawfgrtab_type
189  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, &
190 &                       pawrhoij_init_unpacked, pawrhoij_free_unpacked, &
191 &                       pawrhoij_mpisum_unpacked
192  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get
193  use m_pawfgr,   only : pawfgr_type
194 
195 !This section has been created automatically by the script Abilint (TD).
196 !Do not modify the following lines by hand.
197 #undef ABI_FUNC
198 #define ABI_FUNC 'dfpt_vtorho'
199  use interfaces_18_timing
200  use interfaces_32_util
201  use interfaces_53_ffts
202  use interfaces_65_paw
203  use interfaces_66_wfs
204  use interfaces_67_common
205  use interfaces_72_response, except_this_one => dfpt_vtorho
206 !End of the abilint section
207 
208  implicit none
209 
210 !Arguments -------------------------------
211 !scalars
212  integer,intent(in) :: cplex,dbl_nnsclo,dim_eig2rf,idir,ipert,mband,mk1mem,mkmem
213  integer,intent(in) :: mkqmem,mpw,mpw1,my_natom,natom,ncpgr,nfftf,nkpt_rbz,nspden
214  integer,intent(in) :: nsppol,nsym1,ntypat,optres,prtvol,usecprj,useylmgr1
215  integer,optional,intent(in) :: cg1_out
216  real(dp),intent(in) :: fermie1,ucvol
217  real(dp),intent(out) :: edocc,eeig0,ek0,ek1,eloc0,enl0,enl1,nres2,residm
218  type(MPI_type),intent(in) :: mpi_enreg
219  type(datafiles_type),intent(in) :: dtfil
220  type(dataset_type),intent(in) :: dtset
221  type(efield_type),intent(in) :: dtefield
222  type(pawang_type),intent(in) :: pawang,pawang1
223  type(pawfgr_type),intent(in) :: pawfgr
224  type(pseudopotential_type),intent(in) :: psps
225 
226 !arrays
227  integer,intent(in) :: indsy1(4,nsym1,natom)
228 !                      nfft**(1-1/nsym1) is 1 if nsym1==1, and nfft otherwise
229  integer,intent(in) :: irrzon1(dtset%nfft**(1-1/nsym1),2,(nspden/nsppol)-3*(nspden/4))
230  integer,intent(in) :: istwfk_rbz(nkpt_rbz)
231  integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mk1mem)
232  integer,intent(in) :: nband_rbz(nkpt_rbz*nsppol),npwar1(nkpt_rbz,2)
233  integer,intent(in) :: npwarr(nkpt_rbz,2),symaf1(nsym1),symrc1(3,3,nsym1),symrl1(3,3,nsym1)
234  real(dp),intent(in) :: qphon(3)
235  real(dp),intent(in) :: cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)
236  real(dp),intent(inout) :: cg1(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol)
237  real(dp),intent(inout):: cg1_active(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf)
238  real(dp),intent(out) :: gh1c_set(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf)
239  real(dp),intent(out) :: gh0c1_set(2,mpw1*dtset%nspinor*mband*mk1mem*nsppol*dim_eig2rf)
240  real(dp),intent(in) :: cgq(2,mpw1*dtset%nspinor*mband*mkqmem*nsppol)
241  real(dp),intent(in) :: doccde_rbz(mband*nkpt_rbz*nsppol)
242  real(dp),intent(in) :: docckqde(mband*nkpt_rbz*nsppol)
243  real(dp),intent(in) :: eigen0(mband*nkpt_rbz*nsppol)
244  real(dp),intent(out) :: eigen1(2*mband*mband*nkpt_rbz*nsppol)
245  real(dp),intent(in) :: eigenq(mband*nkpt_rbz*nsppol),gmet(3,3),gprimd(3,3)
246  real(dp),intent(in) :: kpt_rbz(3,nkpt_rbz),occ_rbz(mband*nkpt_rbz*nsppol)
247  real(dp),intent(in) :: occkq(mband*nkpt_rbz*nsppol),ph1d(2,3*(2*dtset%mgfft+1)*natom)
248  real(dp),intent(in) :: phnons1(2,dtset%nfft**(1-1/nsym1),(nspden/nsppol)-3*(nspden/4))
249  real(dp), intent(out) :: nhat1(cplex*nfftf,dtset%nspden*psps%usepaw)
250  real(dp),intent(out) :: resid(mband*nkpt_rbz*nsppol),rhog1(2,nfftf)
251  real(dp),intent(inout) :: nvresid1(cplex*nfftf,nspden),rhor1(cplex*nfftf,nspden)
252  real(dp),intent(in) :: rmet(3,3),rprimd(3,3)
253  real(dp),intent(in),target :: vtrial(nfftf,nspden)
254  real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden)
255  real(dp),intent(in) :: wtk_rbz(nkpt_rbz),xred(3,natom)
256  real(dp),intent(in) :: ylm(mpw*mkmem,psps%mpsang*psps%mpsang*psps%useylm)
257  real(dp),intent(in) :: ylm1(mpw1*mk1mem,psps%mpsang*psps%mpsang*psps%useylm)
258  real(dp),intent(in) :: ylmgr1(mpw1*mk1mem,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
259  integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem,8,3)
260  real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt_rbz,2,3)
261  type(pawcprj_type),intent(in) :: cprj (natom,dtset%nspinor*mband*mkmem *nsppol*usecprj)
262  type(pawcprj_type),intent(in) :: cprjq(natom,dtset%nspinor*mband*mkqmem*nsppol*usecprj)
263  type(pawcprj_type),intent(inout) :: cprj1(natom,dtset%nspinor*mband*mk1mem*nsppol*usecprj)
264  type(paw_ij_type),intent(in) :: paw_ij(my_natom*psps%usepaw),paw_ij1(my_natom*psps%usepaw)
265  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
266  type(pawrhoij_type),intent(in) :: pawrhoij(my_natom*psps%usepaw)
267  type(pawrhoij_type),target,intent(inout) :: pawrhoij1(my_natom*psps%usepaw)
268  type(pawtab_type), intent(in) :: pawtab(ntypat*psps%usepaw)
269  type(wfk_t),intent(inout) :: ddk_f(4)
270 
271 !Local variables-------------------------------
272 !scalars
273  integer,parameter :: level=13
274  integer :: bd2tot_index,bdtot_index,buffer_size,counter,cplex_rhoij
275  integer :: iband,nlines_done,ibdkpt,ibg,ibg1,ibgq,icg,icg1,icgq,ierr
276  integer :: ii,ikg,ikg1,ikpt,ilm,index1,ispden,iscf_mod,isppol,istwf_k
277  integer :: mbd2kpsp,mbdkpsp,mcgq,mcgq_disk,mcprjq
278  integer :: mcprjq_disk,me,n1,n2,n3,n4,n5,n6,nband_k,nband_kq,nkpg,nkpg1
279  integer :: nnsclo_now,npw1_k,npw_k,nspden_rhoij,spaceworld,test_dot
280  logical :: paral_atom,qne0
281  real(dp) :: arg,wtk_k
282  type(gs_hamiltonian_type) :: gs_hamkq
283  type(rf_hamiltonian_type) :: rf_hamkq,rf_hamk_dir2
284 !arrays
285  integer,allocatable :: kg1_k(:,:),kg_k(:,:)
286  integer, pointer :: my_atmtab(:)
287  real(dp) :: kpoint(3),kpq(3)
288  real(dp) :: tsec(2)
289  real(dp),allocatable :: buffer1(:)
290  real(dp),allocatable :: ddkinpw(:),dkinpw(:),dkinpw2(:)
291  real(dp),allocatable :: doccde_k(:),doccde_kq(:)
292  real(dp),allocatable :: edocc_k(:),eeig0_k(:),eig0_k(:),eig0_kq(:),eig1_k(:)
293  real(dp),allocatable :: ek0_k(:),ek1_k(:),eloc0_k(:),enl0_k(:),enl1_k(:)
294  real(dp),allocatable :: ffnl1(:,:,:,:),ffnlk(:,:,:,:)
295  real(dp),allocatable :: grad_berry(:,:,:),kinpw1(:),kpg1_k(:,:)
296  real(dp),allocatable :: kpg_k(:,:),occ_k(:),occ_kq(:)
297  real(dp),allocatable :: ph3d(:,:,:),ph3d1(:,:,:),resid_k(:)
298  real(dp),allocatable :: rho1wfg(:,:),rho1wfr(:,:),rhoaug1(:,:,:,:),rocceig(:,:)
299  real(dp),allocatable :: vlocal(:,:,:,:),vlocal1(:,:,:,:)
300  real(dp),allocatable :: ylm1_k(:,:),ylm_k(:,:),ylmgr1_k(:,:,:)
301  type(pawrhoij_type),pointer :: pawrhoij1_unsym(:)
302 
303 ! *********************************************************************
304 
305  DBG_ENTER('COLL')
306 
307 !Keep track of total time spent in this routine
308  call timab(121,1,tsec)
309  call timab(124,1,tsec)
310 
311 !Retrieve parallelism data
312  spaceworld=mpi_enreg%comm_cell
313  me=mpi_enreg%me_kpt
314  paral_atom=(my_natom/=natom)
315  my_atmtab=>mpi_enreg%my_atmtab
316 
317  if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
318 & dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
319    ABI_ALLOCATE(grad_berry,(2,mpw1,dtefield%mband_occ))
320  else
321    ABI_ALLOCATE(grad_berry,(0,0,0))
322  end if
323 
324 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
325  if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
326    MSG_BUG('wrong values for nfft, nfftf!')
327  end if
328 
329 !The value of iscf must be modified if ddk perturbation, see dfpt_looppert.f
330  iscf_mod=dtset%iscf;if(ipert==natom+1.or.ipert==natom+10.or.ipert==natom+11) iscf_mod=-3
331 
332  edocc=zero ; eeig0=zero ; ek0=zero  ; ek1=zero
333  eloc0=zero ; enl0=zero ; enl1=zero
334  bdtot_index=0
335  bd2tot_index=0
336  ibg=0;icg=0
337  ibgq=0;icgq=0
338  ibg1=0;icg1=0
339  mbdkpsp=mband*nkpt_rbz*nsppol
340  mbd2kpsp=2*mband**2*nkpt_rbz*nsppol
341 
342  n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3)
343  n4=dtset%ngfft(4); n5=dtset%ngfft(5); n6=dtset%ngfft(6)
344  qne0=(qphon(1)**2+qphon(2)**2+qphon(3)**2>=tol14)
345 
346 !Initialize PW 1st-order density if needed
347 !Also store old rho1 in case of density mixing
348  if (iscf_mod>0) then
349    if (optres==1) nvresid1=rhor1
350    if (psps%usepaw==0) then
351      rhor1(:,:)=zero
352    else
353      ABI_ALLOCATE(rho1wfr,(cplex*dtset%nfft,dtset%nspden))
354      ABI_ALLOCATE(rho1wfg,(2,dtset%nfft))
355      rho1wfr(:,:)=zero
356    end if
357  end if
358 
359 !Set max number of non-self-consistent loops nnsclo_now for use in dfpt_vtowfk
360  if(iscf_mod<=0 .and. iscf_mod/=-3)then
361    nnsclo_now=dtset%nstep
362  else
363    if(dtset%nnsclo>0)then
364      nnsclo_now=dtset%nnsclo
365    else
366      nnsclo_now=1
367    end if
368    if(dbl_nnsclo==1) nnsclo_now=nnsclo_now*2
369  end if
370 
371 !Prepare GS k+q wf
372  mcgq=mpw1*dtset%nspinor*mband*mkqmem*nsppol;mcgq_disk=0
373 
374 !Prepare RF PAW files
375  if (psps%usepaw==1) then
376    mcprjq=dtset%nspinor*mband*mkqmem*nsppol*usecprj;mcprjq_disk=0
377  else
378    mcprjq=0;mcprjq_disk=0
379  end if
380 
381 !Initialisation of the wfdot file in case of electric field (or 2nd order Sternheimer equation)
382  test_dot=0
383  if (ipert==natom+2.and.sum((dtset%qptn(1:3))**2 )<=tol7.and.&
384 & (dtset%berryopt/= 4.and.dtset%berryopt/= 6.and.dtset%berryopt/= 7.and.&
385 & dtset%berryopt/=14.and.dtset%berryopt/=16.and.dtset%berryopt/=17).or.&
386 & (ipert==natom+10.or.ipert==natom+11)) then
387    test_dot=1
388  end if
389 
390 !==== Initialize most of the Hamiltonian (and derivative) ====
391 !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
392 !2) Perform the setup needed for the non-local factors:
393 !* Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
394 !* PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
395 
396  call init_hamiltonian(gs_hamkq,psps,pawtab,dtset%nspinor,nsppol,nspden,natom,&
397 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
398 & paw_ij=paw_ij,comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
399 & usecprj=usecprj,ph1d=ph1d,nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
400 
401  call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamkq,has_e1kbsc=.true.,paw_ij1=paw_ij1,&
402 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
403  if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then
404    call init_rf_hamiltonian(cplex,gs_hamkq,ipert,rf_hamk_dir2,has_e1kbsc=.true.,paw_ij1=paw_ij1,&
405 &   comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab)
406  end if
407 
408 !PAW:allocate memory for non-symetrized 1st-order occupancies matrix (pawrhoij1)
409  pawrhoij1_unsym => pawrhoij1
410  if (psps%usepaw==1.and.iscf_mod>0) then
411    if (paral_atom) then
412      ABI_DATATYPE_ALLOCATE(pawrhoij1_unsym,(natom))
413      cplex_rhoij=max(cplex,dtset%pawcpxocc);nspden_rhoij=dtset%nspden
414      call pawrhoij_alloc(pawrhoij1_unsym,cplex_rhoij,nspden_rhoij,dtset%nspinor,&
415 &     dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0,use_rhoij_=1)
416    else
417      pawrhoij1_unsym => pawrhoij1
418      call pawrhoij_init_unpacked(pawrhoij1_unsym)
419    end if
420  end if
421 
422  ABI_ALLOCATE(rhoaug1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
423  ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamkq%nvloc))
424  ABI_ALLOCATE(vlocal1,(cplex*n4,n5,n6,gs_hamkq%nvloc))
425 
426  nlines_done = 0
427 
428 !LOOP OVER SPINS
429  do isppol=1,nsppol
430 
431 !  Rewind kpgsph data file if needed:
432    ikg=0;ikg1=0
433 
434 !  Set up local potential vlocal1 with proper dimensioning, from vtrial1
435 !  Same thing for vlocal from vtrial Also take into account the spin.
436 
437    call rf_transgrid_and_pack(isppol,nspden,psps%usepaw,cplex,nfftf,dtset%nfft,dtset%ngfft,&
438 &   gs_hamkq%nvloc,pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1)
439 
440 !  Continue to initialize the Hamiltonian
441    call load_spin_hamiltonian(gs_hamkq,isppol,vlocal=vlocal,with_nonlocal=.true.)
442    call load_spin_rf_hamiltonian(rf_hamkq,isppol,vlocal1=vlocal1,with_nonlocal=.true.)
443    if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then
444      call load_spin_rf_hamiltonian(rf_hamk_dir2,isppol,with_nonlocal=.true.)
445      if (ipert==natom+11) then ! load vlocal1
446        call load_spin_rf_hamiltonian(rf_hamk_dir2,isppol,vlocal1=vlocal1)
447      end if
448    end if
449 
450    if (ipert==natom+5) then !SPr deb, in case of magnetic field perturbation, no non-local
451      call load_spin_rf_hamiltonian(rf_hamkq,isppol,vlocal1=vlocal1)
452    end if
453 
454 !  Nullify contribution to 1st-order density from this k-point
455    rhoaug1(:,:,:,:)=zero
456 
457    call timab(125,1,tsec)
458 
459 !======================================================================
460 !==============  BIG FAT K POINT LOOP  ================================
461 !======================================================================
462 
463    do ikpt=1,nkpt_rbz
464      counter=100*ikpt+isppol
465 
466      nband_k = nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
467      istwf_k = istwfk_rbz(ikpt)
468      npw_k   = npwarr(ikpt,1)
469      npw1_k  = npwar1(ikpt,1)
470      wtk_k   = wtk_rbz(ikpt)
471 
472      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then
473        eigen1(1+bd2tot_index : 2*nband_k**2+bd2tot_index) = zero
474        resid(1+bdtot_index : nband_k+bdtot_index) = zero
475        bdtot_index=bdtot_index+nband_k
476        bd2tot_index=bd2tot_index+2*nband_k**2
477 
478        cycle ! Skip the rest of the k-point loop
479      end if
480 
481      kpoint(:)=kpt_rbz(:,ikpt)
482      kpq(:)=kpoint(:);if (ipert<natom+3.or.ipert==natom+5) kpq(:)=kpq(:)+qphon(1:3)
483      ABI_ALLOCATE(kg_k,(3,npw_k))
484      ABI_ALLOCATE(kg1_k,(3,npw1_k))
485      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
486      ABI_ALLOCATE(ylm1_k,(npw1_k,psps%mpsang*psps%mpsang*psps%useylm))
487      ABI_ALLOCATE(ylmgr1_k,(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1))
488      ABI_ALLOCATE(doccde_k,(nband_k))
489      ABI_ALLOCATE(doccde_kq,(nband_k))
490      ABI_ALLOCATE(eig0_k,(nband_k))
491      ABI_ALLOCATE(eig0_kq,(nband_k))
492      ABI_ALLOCATE(eig1_k,(2*nband_k**2))
493      ABI_ALLOCATE(edocc_k,(nband_k))
494      ABI_ALLOCATE(eeig0_k,(nband_k))
495      ABI_ALLOCATE(ek0_k,(nband_k))
496      ABI_ALLOCATE(ek1_k,(nband_k))
497      ABI_ALLOCATE(eloc0_k,(nband_k))
498      ABI_ALLOCATE(occ_k,(nband_k))
499      ABI_ALLOCATE(occ_kq,(nband_k))
500      ABI_ALLOCATE(resid_k,(nband_k))
501      ABI_ALLOCATE(rocceig,(nband_k,nband_k))
502      ABI_ALLOCATE(enl0_k,(nband_k))
503      ABI_ALLOCATE(enl1_k,(nband_k))
504 
505      eig1_k(:)=zero
506      eig0_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index)
507      eig0_kq(:)=eigenq(1+bdtot_index:nband_k+bdtot_index)
508      edocc_k(:)=zero
509      eeig0_k(:)=zero ; ek0_k(:)=zero  ; ek1_k(:)=zero
510      eloc0_k(:)=zero ; enl0_k(:)=zero ; enl1_k(:)=zero
511      occ_k(:)=occ_rbz(1+bdtot_index:nband_k+bdtot_index)
512      occ_kq(:)=occkq(1+bdtot_index:nband_k+bdtot_index)
513      doccde_k(:)=doccde_rbz(1+bdtot_index:nband_k+bdtot_index)
514      doccde_kq(:)=docckqde(1+bdtot_index:nband_k+bdtot_index)
515      resid_k(:)=zero
516 
517 !    For each pair of active bands (m,n), generates the ratios
518 !    rocceig(m,n)=(occ_kq(m)-occ_k(n))/(eig0_kq(m)-eig0_k(n))
519 !    and decide to which band to attribute it.
520      call occeig(doccde_k,doccde_kq,eig0_k,eig0_kq,nband_k,dtset%occopt,occ_k,occ_kq,rocceig)
521 
522      ! These arrays are not needed anymore.
523      ABI_DEALLOCATE(doccde_k)
524      ABI_DEALLOCATE(doccde_kq)
525      ABI_DEALLOCATE(occ_kq)
526 
527      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
528      if (psps%useylm==1) then
529        do ilm=1,psps%mpsang*psps%mpsang
530          ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
531        end do
532      end if
533 
534 !    Get (k+q+G) wave vectors and associated spherical harmonics
535      kg1_k(:,1:npw1_k)=kg1(:,1+ikg1:npw1_k+ikg1)
536      if (psps%useylm==1) then
537        do ilm=1,psps%mpsang*psps%mpsang
538          ylm1_k(1:npw1_k,ilm)=ylm1(1+ikg1:npw1_k+ikg1,ilm)
539        end do
540        if (useylmgr1==1) then
541          do ilm=1,psps%mpsang*psps%mpsang
542            do ii=1,3+6*((ipert-natom)/10)
543              ylmgr1_k(1:npw1_k,ii,ilm)=ylmgr1(1+ikg1:npw1_k+ikg1,ii,ilm)
544            end do
545          end do
546        end if
547      end if
548 
549 !    Set up the ground-state Hamiltonian, and some parts of the 1st-order Hamiltonian
550      call getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,&                              ! In
551      kpoint,kpq,idir,ipert,natom,rmet,gprimd,gmet,istwf_k,&                         ! In
552      npw_k,npw1_k,useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,&                      ! In
553      dkinpw,nkpg,nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,&                 ! Out
554      ddkinpw=ddkinpw,dkinpw2=dkinpw2,rf_hamk_dir2=rf_hamk_dir2)                     ! Out
555 
556 !    Compute the gradient of the Berry-phase term
557      if (dtset%berryopt== 4.or.dtset%berryopt== 6.or.dtset%berryopt== 7.or.&
558 &     dtset%berryopt==14.or.dtset%berryopt==16.or.dtset%berryopt==17) then
559        if (ipert<=natom) then
560 !        phonon perturbation
561          call dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt_rbz,&
562 &         npwarr,npwar1,dtset%nspinor,nsppol,qmat,pwindall)
563        else
564 !        electric field perturbation
565          call dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir,ikpt,isppol,mband,mpw,mpw1,mkmem,mk1mem,nkpt_rbz,&
566 &         npwarr,npwar1,dtset%nspinor,&
567 &         nsppol,qmat,pwindall,rprimd)
568        end if
569      end if
570 
571      ! Free some memory before calling dfpt_vtowfk
572      ABI_DEALLOCATE(ylm_k)
573      ABI_DEALLOCATE(ylm1_k)
574      ABI_DEALLOCATE(ylmgr1_k)
575 
576 !    Compute the eigenvalues, wavefunction, residuals,
577 !    contributions to kinetic energy, nonlocal energy, forces,
578 !    and update of 1st-order density to this k-point and this spin polarization.
579      nband_kq = nband_k  !Note that the calculation only works for same number of bandes on all K points.
580 !    Note that dfpt_vtowfk is called with kpoint, while kpt is used inside vtowfk3
581      call dfpt_vtowfk(cg,cgq,cg1,cg1_active,cplex,cprj,cprjq,cprj1,dim_eig2rf,dtfil,&
582 &     dtset,edocc_k,eeig0_k,eig0_k,eig0_kq,eig1_k,ek0_k,ek1_k,eloc0_k,enl0_k,enl1_k,fermie1,&
583 &     gh0c1_set,gh1c_set,grad_berry,gs_hamkq,ibg,ibgq,ibg1,icg,icgq,icg1,idir,ikpt,ipert,isppol,&
584 &     mband,mcgq,mcprjq,mkmem,mk1mem,mpi_enreg,mpw,mpw1,natom,nband_k,ncpgr,nnsclo_now,&
585 &     npw_k,npw1_k,dtset%nspinor,nsppol,n4,n5,n6,occ_k,pawrhoij1_unsym,prtvol,psps,resid_k,&
586 &     rf_hamkq,rf_hamk_dir2,rhoaug1,rocceig,ddk_f,wtk_k,nlines_done,cg1_out)
587 
588 !    Free temporary storage
589      ABI_DEALLOCATE(kinpw1)
590      ABI_DEALLOCATE(kg_k)
591      ABI_DEALLOCATE(kg1_k)
592      ABI_DEALLOCATE(kpg_k)
593      ABI_DEALLOCATE(kpg1_k)
594      ABI_DEALLOCATE(dkinpw)
595      if (ipert==natom+10) then
596        ABI_DEALLOCATE(ddkinpw)
597        if (idir>3) then
598          ABI_DEALLOCATE(dkinpw2)
599        end if
600      end if
601      ABI_DEALLOCATE(ffnlk)
602      ABI_DEALLOCATE(ffnl1)
603      ABI_DEALLOCATE(eig0_k)
604      ABI_DEALLOCATE(eig0_kq)
605      ABI_DEALLOCATE(rocceig)
606      ABI_DEALLOCATE(ph3d)
607      if (allocated(ph3d1)) then
608        ABI_DEALLOCATE(ph3d1)
609      end if
610 
611 !    Save eigenvalues (hartree), residuals (hartree**2)
612      eigen1 (1+bd2tot_index : 2*nband_k**2+bd2tot_index) = eig1_k(:)
613      resid  (1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
614 
615 !    Accumulate sum over k points for nonlocal and kinetic energies,
616 !    also accumulate gradients of Enonlocal:
617      if (iscf_mod>0 .or. iscf_mod==-3 .or. iscf_mod==-2)then
618        do iband=1,nband_k
619          edocc=edocc+wtk_k*occ_k(iband)*edocc_k(iband)
620          eeig0=eeig0+wtk_k*occ_k(iband)*eeig0_k(iband)
621          ek0=ek0+wtk_k*occ_k(iband)*ek0_k(iband)
622          ek1=ek1+wtk_k*occ_k(iband)*ek1_k(iband)
623          eloc0=eloc0+wtk_k*occ_k(iband)*eloc0_k(iband)
624          enl0=enl0+wtk_k*occ_k(iband)*enl0_k(iband)
625          enl1=enl1+wtk_k*occ_k(iband)*enl1_k(iband)
626        end do
627      end if
628 
629      ABI_DEALLOCATE(eig1_k)
630      ABI_DEALLOCATE(occ_k)
631      ABI_DEALLOCATE(resid_k)
632      ABI_DEALLOCATE(edocc_k)
633      ABI_DEALLOCATE(eeig0_k)
634      ABI_DEALLOCATE(ek0_k)
635      ABI_DEALLOCATE(ek1_k)
636      ABI_DEALLOCATE(eloc0_k)
637      ABI_DEALLOCATE(enl0_k)
638      ABI_DEALLOCATE(enl1_k)
639 
640 !    Keep track of total number of bands (all k points so far, even for k points not treated by me)
641      bdtot_index=bdtot_index+nband_k
642      bd2tot_index=bd2tot_index+2*nband_k**2
643 
644 !    Shift array memory
645      if (mkmem/=0) then
646        ibg=ibg+nband_k
647        icg=icg+npw_k*dtset%nspinor*nband_k
648        ikg=ikg+npw_k
649      end if
650      if (mkqmem/=0) then
651        ibgq=ibgq+dtset%nspinor*nband_k
652        icgq=icgq+npw1_k*dtset%nspinor*nband_k
653      end if
654      if (mk1mem/=0) then
655        ibg1=ibg1+nband_k
656        icg1=icg1+npw1_k*dtset%nspinor*nband_k
657        ikg1=ikg1+npw1_k
658      end if
659 
660    end do
661 
662 !======================================================================
663 !==================  END BIG K POINT LOOP  ============================
664 !======================================================================
665 
666    call timab(125,2,tsec)
667 
668 !  Transfer density on augmented fft grid to normal fft grid in real space. Also take into account the spin.
669 ! FR EB for the non-collinear part see vtorho.F90
670    if(iscf_mod>0) then
671      if (psps%usepaw==0) then
672        call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,1),1)
673        if(nspden==4)then
674          do ispden=2,4
675            call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rhor1,rhoaug1(:,:,:,ispden),1)
676          end do
677        end if
678      else
679        call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,1),1)
680        if(nspden==4)then
681          do ispden=2,4
682            call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,dtset%ngfft,rho1wfr,rhoaug1(:,:,:,ispden),1)
683          end do
684        end if
685      end if
686    end if
687 
688  end do !  End loop over spins
689 
690 !More memory cleaning
691  call destroy_hamiltonian(gs_hamkq)
692  call destroy_rf_hamiltonian(rf_hamkq)
693  if ((ipert==natom+10.and.idir>3).or.ipert==natom+11) then
694    call destroy_rf_hamiltonian(rf_hamk_dir2)
695  end if
696  ABI_DEALLOCATE(rhoaug1)
697  ABI_DEALLOCATE(vlocal)
698  ABI_DEALLOCATE(vlocal1)
699 
700  call timab(124,2,tsec)
701 
702 !=== MPI communications ==================
703  if(xmpi_paral==1)then
704    call timab(129,1,tsec)
705 
706 !  Compute buffer size
707    buffer_size=7+mbd2kpsp+mbdkpsp
708    if (iscf_mod>0) then
709      buffer_size=buffer_size+cplex*dtset%nfft*nspden
710    end if
711    ABI_ALLOCATE(buffer1,(buffer_size))
712 
713 !  Pack rhor1,edocc,eeig0,ek0,ek1,eloc0,enl0,enl1,eigen1,resid
714    if (iscf_mod>0) then
715      index1=cplex*dtset%nfft*nspden
716      if (psps%usepaw==0) then
717        buffer1(1:index1)=reshape(rhor1  ,(/index1/))
718      else
719        buffer1(1:index1)=reshape(rho1wfr,(/index1/))
720      end if
721    else
722      index1=0
723    end if
724    buffer1(index1+1)=edocc;buffer1(index1+2)=eeig0
725    buffer1(index1+3)=ek0  ;buffer1(index1+4)=ek1
726    buffer1(index1+5)=eloc0;buffer1(index1+6)=enl0
727    buffer1(index1+7)=enl1
728    index1=index1+7
729    bdtot_index=0;bd2tot_index=0
730    do isppol=1,nsppol
731      do ikpt=1,nkpt_rbz
732        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
733        buffer1(index1+1:index1+2*nband_k**2) = eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2)
734        buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k)= resid(bdtot_index+1:bdtot_index+nband_k)
735        bdtot_index=bdtot_index+nband_k
736        bd2tot_index=bd2tot_index+2*nband_k**2
737        index1=index1+2*nband_k**2+nband_k
738      end do
739    end do
740    if(index1<buffer_size)buffer1(index1+1:buffer_size)=zero
741 
742 !  Build sum of everything
743    call timab(48,1,tsec)
744    call xmpi_sum(buffer1,buffer_size,spaceworld,ierr)
745    call timab(48,2,tsec)
746 
747 !  Unpack the final result
748    if(iscf_mod>0) then
749      index1=cplex*dtset%nfft*nspden
750      if (psps%usepaw==0) then
751        rhor1(:,:)  =reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/))
752      else
753        rho1wfr(:,:)=reshape(buffer1(1:index1),(/cplex*dtset%nfft,nspden/))
754      end if
755    else
756      index1=0
757    end if
758 
759    edocc=buffer1(index1+1);eeig0=buffer1(index1+2)
760    ek0=buffer1(index1+3)  ;ek1=buffer1(index1+4)
761    eloc0=buffer1(index1+5);enl0=buffer1(index1+6)
762    enl1=buffer1(index1+7)
763    index1=index1+7
764    bdtot_index=0;bd2tot_index=0
765    do isppol=1,nsppol
766      do ikpt=1,nkpt_rbz
767        nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
768        eigen1(bd2tot_index+1:bd2tot_index+2*nband_k**2) = buffer1(index1+1:index1+2*nband_k**2)
769        resid(bdtot_index+1:bdtot_index+nband_k)= buffer1(index1+2*nband_k**2+1:index1+2*nband_k**2+nband_k)
770        bdtot_index=bdtot_index+nband_k
771        bd2tot_index=bd2tot_index+2*nband_k**2
772        index1=index1+2*nband_k**2+nband_k
773      end do
774    end do
775    ABI_DEALLOCATE(buffer1)
776 
777 !  Accumulate PAW occupancies
778    if (psps%usepaw==1.and.iscf_mod>0) then
779      call pawrhoij_mpisum_unpacked(pawrhoij1_unsym,spaceworld)
780    end if
781 
782    call timab(129,2,tsec)
783  end if ! if kpt parallel
784 
785  call timab(127,1,tsec)
786 
787 !If needed, compute rhog1, and symmetrizes the density
788  if (iscf_mod > 0) then
789 
790 !  In order to have the symrhg working in parallel on FFT coefficients, the size
791 !  of irzzon1 and phnons1 should be set to nfftot. Therefore, nsym\=1 does not work.
792 
793    if(nspden==4) then
794 ! FR symrhg will manage correctly this rearrangement
795      rhor1(:,2)=rhor1(:,2)+(rhor1(:,1)+rhor1(:,4))    !(n+mx)
796      rhor1(:,3)=rhor1(:,3)+(rhor1(:,1)+rhor1(:,4))    !(n+my)
797      call timab(17,2,tsec)
798    end if
799 !
800    if (psps%usepaw==0) then
801      call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
802 &     nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,rhog1  ,rhor1  ,rprimd,symaf1,symrl1)
803    else
804      call symrhg(cplex,gprimd,irrzon1,mpi_enreg,dtset%nfft,dtset%nfft,dtset%ngfft,&
805 &     nspden,nsppol,nsym1,dtset%paral_kgb,phnons1,rho1wfg,rho1wfr,rprimd,symaf1,symrl1)
806    end if
807 !  We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
808 !  we also have the spin-up density, symmetrized, in rhor1(:,2).
809  end if
810 
811  ABI_DEALLOCATE(grad_berry)
812 
813 !Find largest residual over bands, k points, and spins except for nbdbuf highest bands
814  ibdkpt=1
815  residm=zero
816  do isppol=1,nsppol
817    do ikpt=1,nkpt_rbz
818      nband_k=nband_rbz(ikpt+(isppol-1)*nkpt_rbz)
819      nband_k=max(1,nband_k-dtset%nbdbuf)
820      residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_k-1)))
821      ibdkpt=ibdkpt+nband_k
822    end do
823  end do
824 
825  call timab(127,2,tsec)
826 
827  if (iscf_mod>0) then
828 
829 !  PAW: Build new 1st-order rhoij quantities then symetrize them
830 !  Compute and add the 1st-order compensation density to rho1wfr
831 !  to get the total 1st-order density
832    if (psps%usepaw==1) then
833      call pawmkrho(arg,cplex,gprimd,idir,indsy1,ipert,mpi_enreg,&
834 &     my_natom,natom,nspden,nsym1,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
835 &     dtset%pawprtvol,pawrhoij1,pawrhoij1_unsym,pawtab,dtset%qptn,rho1wfg,rho1wfr,&
836 &     rhor1,rprimd,symaf1,symrc1,dtset%typat,ucvol,dtset%usewvl,xred,&
837 &     pawang_sym=pawang1,pawnhat=nhat1,pawrhoij0=pawrhoij,rhog=rhog1)
838      ABI_DEALLOCATE(rho1wfr)
839      ABI_DEALLOCATE(rho1wfg)
840      if (paral_atom) then
841        call pawrhoij_free(pawrhoij1_unsym)
842        ABI_DATATYPE_DEALLOCATE(pawrhoij1_unsym)
843      end if
844    end if
845 
846 !  Compute density residual (if required) and its squared norm
847    if (optres==1) then
848      nvresid1=rhor1-nvresid1
849      call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid1)
850    end if
851  end if ! iscf>0
852 
853  call timab(121,2,tsec)
854 
855  DBG_EXIT('COLL')
856 
857 end subroutine dfpt_vtorho