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

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

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

ABINIT/m_dfpt_vtorho [ Modules ]

[ Top ] [ Modules ]

NAME

 m_dfpt_vtorho

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_dfpt_vtorho
27 
28  use defs_basis
29  use defs_datatypes
30  use defs_abitypes
31  use m_abicore
32  use m_xmpi
33  use m_errors
34  use m_efield
35  use m_hamiltonian
36  use m_wfk
37  use m_cgtools
38 
39  use m_time,     only : timab
40  use m_occ,      only : occeig
41  use m_hdr,      only : hdr_skip, hdr_io
42  use m_pawang,   only : pawang_type
43  use m_pawtab,   only : pawtab_type
44  use m_paw_ij,   only : paw_ij_type
45  use m_pawfgrtab,only : pawfgrtab_type
46  use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, &
47 &                       pawrhoij_init_unpacked, pawrhoij_free_unpacked, &
48 &                       pawrhoij_mpisum_unpacked
49  use m_pawcprj,  only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_get
50  use m_pawfgr,   only : pawfgr_type
51  use m_paw_mkrho,only : pawmkrho
52  use m_fft,      only : fftpac
53  use m_spacepar, only : symrhg
54  use m_getgh1c,  only : rf_transgrid_and_pack, getgh1c_setup
55  use m_dfpt_vtowfk, only : dfpt_vtowfk
56  use m_dfpt_fef,    only : dfptff_gradberry, dfptff_gbefd
57  use m_mpinfo,      only : proc_distrb_cycle
58  use m_fourier_interpol, only : transgrid
59 
60  implicit none
61 
62  private