TABLE OF CONTENTS


ABINIT/cgq_builder [ Functions ]

[ Top ] [ Functions ]

NAME

 cgq_builder

FUNCTION

 This routine locates cgq for efield calculations, especially for parallel case

INPUTS

  berryflag = logical flag determining use of electric field variables
  cg(2,mcg)=planewave coefficients of wavefunctions.
  dtset <type(dataset_type)>=all input variables for this dataset
  ikpt=index of current k kpt
  ikpt_loc=index of k point on current processor (see vtorho.F90)
  isspol=value of spin polarization currently treated
  me_distrb=current value from spaceComm_distrb (see vtorho.F90)
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcgq=size of cgq array (see vtorho.F90)
  mkgq=size of pwnsfacq array (see vtorho.F90)
  my_nspinor=nspinor value determined by current // set up
  nband_k=number of bands at each k point
  nproc_distrb=nproc from spaceComm_distrb (see vtorho.F90)
  npwarr(nkpt)=number of planewaves in basis at this k point
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
                           (see initberry.f)
  pwind_alloc = first dimension of pwind
  spaceComm_distrb=comm_cell from mpi_enreg

OUTPUT

  cgq(2,mcgq)=planewave coefficients of wavenfunctions adjacent to cg at ikpt
  pwnsfacq(2,mkgq)=phase factors for non-symmorphic translations for cg's adjacent to cg(ikpt)

SIDE EFFECTS

 Input/Output
   dtefield <type(efield_type)> = efield variables
   mpi_enreg=information about MPI parallelization

TODO

NOTES

PARENTS

      vtorho

CHILDREN

      timab,xmpi_recv,xmpi_send

SOURCE

2345 subroutine cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,&
2346 &                      me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,&
2347 &                      npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb)
2348 
2349 
2350 !This section has been created automatically by the script Abilint (TD).
2351 !Do not modify the following lines by hand.
2352 #undef ABI_FUNC
2353 #define ABI_FUNC 'cgq_builder'
2354 !End of the abilint section
2355 
2356  implicit none
2357 
2358 !Arguments ------------------------------------
2359  integer,intent(in) :: ikpt,ikpt_loc,isppol,me_distrb,mcg,mcgq,mkgq,my_nspinor,nband_k
2360  integer,intent(in) :: nproc_distrb,pwind_alloc,spaceComm_distrb
2361  logical,intent(in) :: berryflag
2362  type(dataset_type), intent(in) :: dtset
2363  type(efield_type), intent(inout) :: dtefield
2364  type(MPI_type), intent(in) :: mpi_enreg
2365 !arrays
2366  integer,intent(in) :: npwarr(dtset%nkpt)
2367  real(dp),intent(in) :: cg(2,mcg),pwnsfac(2,pwind_alloc)
2368  real(dp),intent(out) :: cgq(2,mcgq),pwnsfacq(2,mkgq)
2369 
2370 !Local variables -------------------------
2371 !scalars
2372  integer :: count,count1,icg1,icg2,dest,his_source
2373  integer :: idir,ierr,ifor,ikg1,ikg2,ikptf,ikpt1f,ikpt1i
2374  integer :: jkpt,jkpt1i,jkptf,jkpt1f,jsppol,my_source,npw_k1,tag
2375 !arrays
2376  integer,allocatable :: flag_send(:,:), flag_receive(:)
2377  real(dp) :: tsec(2)
2378  real(dp),allocatable :: buffer(:,:)
2379 
2380 ! *************************************************************************
2381 
2382 !DEBUG
2383 !write(std_out,'(a)')'cgq_builder enter'
2384 !DEBUG
2385 
2386  if (mcgq==0.or.mkgq==0) return
2387 
2388  call timab(983,1,tsec)
2389 
2390 !Test compatbility of berryflag
2391  if (berryflag) then
2392    ABI_ALLOCATE(flag_send,(0:nproc_distrb-1,dtefield%fnkpt))
2393  end if
2394  ABI_ALLOCATE(flag_receive,(dtset%nkpt))
2395  flag_send(:,:) = 0
2396  flag_receive(:) = 0
2397 
2398  if (berryflag) ikptf = dtefield%i2fbz(ikpt)
2399 
2400  do idir = 1, 3
2401 
2402 !  skip idir values for which efield_dot(idir) = 0
2403    if (berryflag .and. abs(dtefield%efield_dot(idir)) < tol12 ) cycle
2404 
2405    do ifor = 1, 2
2406 
2407      if(berryflag) then
2408        dtefield%sflag(:,ikpt + dtset%nkpt*(isppol - 1),ifor,idir) = 0
2409        ikpt1f = dtefield%ikpt_dk(ikptf,ifor,idir)
2410        ikpt1i = dtefield%indkk_f2ibz(ikpt1f,1)
2411      end if
2412 
2413      npw_k1 = npwarr(ikpt1i)
2414      count = npw_k1*my_nspinor*nband_k
2415      my_source = mpi_enreg%proc_distrb(ikpt1i,1,isppol)
2416 
2417      do dest = 0, nproc_distrb-1
2418 
2419        if ((dest==me_distrb).and.(ikpt_loc <= dtset%mkmem)) then
2420 !        I am dest and have something to do
2421 
2422          if ( my_source == me_distrb ) then
2423 !          I am destination and source
2424 
2425            if(berryflag) then
2426              ikg1 = dtefield%fkgindex(ikpt1f)
2427              ikg2 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2428              icg1 = dtefield%cgindex(ikpt1i,isppol)
2429              icg2 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2430            end if
2431 
2432            pwnsfacq(:,ikg2 + 1:ikg2 + npw_k1) = pwnsfac(:,ikg1 + 1:ikg1 + npw_k1)
2433            cgq(:,icg2 + 1:icg2 + count) = cg(:,icg1 + 1:icg1 + count)
2434 
2435          else !  I am the destination but not the source -> receive
2436 !          receive pwnsfacq
2437            if(berryflag) then
2438              tag = ikpt1f + (isppol - 1)*dtefield%fnkpt
2439              ikg1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2440            end if
2441            ABI_ALLOCATE(buffer,(2,npw_k1))
2442            call xmpi_recv(buffer,my_source,tag,spaceComm_distrb,ierr)
2443            pwnsfacq(:,ikg1+1:ikg1+npw_k1) = buffer(:,1:npw_k1)
2444            ABI_DEALLOCATE(buffer)
2445 
2446 !          receive cgq if necessary
2447            if(flag_receive(ikpt1i) == 0) then
2448              ABI_ALLOCATE(buffer,(2,count))
2449              tag = ikpt1i + (isppol - 1)*dtset%nkpt
2450              call xmpi_recv(buffer,my_source,tag,spaceComm_distrb,ierr)
2451              if(berryflag) icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*dtset%nkpt)
2452              cgq(:,icg1+1:icg1+count) = buffer(:,1:count)
2453              ABI_DEALLOCATE(buffer)
2454              flag_receive(ikpt1i) = 1
2455            end if ! end if flag_receive == 0
2456          end if ! end tasks if I am the destination
2457 
2458        else if (ikpt_loc <= mpi_enreg%mkmem(dest)) then  ! dest != me and the dest has a k-point to treat
2459 
2460 !        jkpt is the kpt which is being treated by dest (in ibz)
2461 !        jsppol is his isppol
2462          jkpt = mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,1)
2463          jsppol = mpi_enreg%kpt_loc2ibz_sp(dest, ikpt_loc,2)
2464 
2465          if(jkpt > 0 .and. jsppol > 0) then
2466 
2467            if(berryflag) then
2468              jkptf = dtefield%i2fbz(jkpt)
2469              jkpt1f = dtefield%ikpt_dk(jkptf,ifor,idir)
2470              jkpt1i = dtefield%indkk_f2ibz(jkpt1f,1)
2471            end if
2472            his_source = mpi_enreg%proc_distrb(jkpt1i,1,jsppol)
2473 
2474            if (his_source == me_distrb) then
2475 
2476 !            send
2477 !            pwnsfacq
2478              if(berryflag) then
2479                ikg1 = dtefield%fkgindex(jkpt1f)
2480                tag = jkpt1f + (jsppol - 1)*dtefield%fnkpt
2481              end if
2482              count1 = npwarr(jkpt1i)
2483              ABI_ALLOCATE(buffer,(2,count1))
2484              buffer(:,1:count1)  = pwnsfac(:,ikg1+1:ikg1+count1)
2485              call xmpi_send(buffer,dest,tag,spaceComm_distrb,ierr)
2486              ABI_DEALLOCATE(buffer)
2487 
2488 !            send cgq if necessary
2489              if(flag_send(dest, jkpt1i)==0) then
2490                if(berryflag) icg1 = dtefield%cgindex(jkpt1i,jsppol)
2491                tag = jkpt1i + (jsppol - 1)*dtset%nkpt
2492                count1 = npwarr(jkpt1i)*nband_k*my_nspinor
2493                ABI_ALLOCATE(buffer,(2,count1))
2494                buffer(:,1:count1)  = cg(:,icg1+1:icg1+count1)
2495                call xmpi_send(buffer,dest,tag,spaceComm_distrb,ierr)
2496                ABI_DEALLOCATE(buffer)
2497                flag_send(dest, jkpt1i)=1
2498              end if ! if send cgq
2499 
2500            end if ! end check that his_source == me
2501 
2502          end if ! end check on jkpt > 0 and jsppol > 0
2503 
2504        end if ! end check on me = dest else if me != dest
2505 
2506      end do ! end loop over dest = 0, nproc-1
2507 
2508    end do !end loop over ifor
2509 
2510  end do !end loop over idir
2511 
2512  call timab(983,2,tsec)
2513 
2514  ABI_DEALLOCATE(flag_send)
2515  ABI_DEALLOCATE(flag_receive)
2516 
2517 !DEBUG
2518 !write(std_out,'(a)')'cgq_builder exit'
2519 !END_DEBUG
2520 
2521 end subroutine cgq_builder

ABINIT/e_eigen [ Functions ]

[ Top ] [ Functions ]

NAME

  e_eigen

FUNCTION

  Computes eigenvalues energy from eigen, occ, kpt, wtk

INPUTS

  eigen(nkpt*nsppol)=eigenvalues
  mband= maximum number of bands
  nband(nkpt*nsppol)= number of bands for each k-point and spin
  nkpt= number of k-points
  nsppol= number of spin polarization
  occ(mband*nkpt*nsppol)=occupations
  wtk(nkpt)= k-point weights

OUTPUT

  e_eigenvalues= eigenvalues energy

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2052 subroutine e_eigen(eigen,e_eigenvalues,mband,nband,nkpt,nsppol,occ,wtk)
2053 
2054 
2055 !This section has been created automatically by the script Abilint (TD).
2056 !Do not modify the following lines by hand.
2057 #undef ABI_FUNC
2058 #define ABI_FUNC 'e_eigen'
2059 !End of the abilint section
2060 
2061  implicit none
2062 
2063 !Arguments ------------------------------------
2064  integer , intent(in)  :: mband,nkpt,nsppol
2065  integer , intent(in)  :: nband(nkpt*nsppol)
2066  real(dp) , intent(in)  :: eigen(mband*nkpt*nsppol)
2067  real(dp) , intent(in)  :: occ(mband*nkpt*nsppol)
2068  real(dp) , intent(in)  :: wtk(nkpt)
2069  real(dp) , intent(out) :: e_eigenvalues
2070 
2071 !Local variables-------------------------------
2072  integer :: ib,iband,ii,ikpt,isppol,nband_k
2073  real(dp) :: wtk_k
2074 
2075 ! *************************************************************************
2076 
2077    DBG_ENTER("COLL")
2078    ii=0;ib=0
2079    do isppol=1,nsppol
2080      do ikpt=1,nkpt
2081        ii=ii+1
2082        nband_k=nband(ii) ;  wtk_k=wtk(ii)
2083        do iband=1,nband_k
2084          ib=ib+1
2085          if(abs(occ(ib)) > tol8) then
2086            e_eigenvalues = e_eigenvalues + wtk_k*occ(ib)*eigen(ib)
2087          end if
2088        end do
2089      end do
2090    end do
2091 
2092    DBG_EXIT("COLL")
2093 
2094  end subroutine e_eigen

ABINIT/m_vtorho [ Modules ]

[ Top ] [ Modules ]

NAME

  m_vtorho

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MF, AR, MM, MT, FJ, MB, MT, TR)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_vtorho
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_wvltypes
33  use m_abicore
34  use m_xmpi
35  use m_ab7_mixing
36  use m_errors
37  use m_wffile
38  use m_efield
39  use m_cgtools
40  use m_gemm_nonlop
41 
42  use m_time,               only : timab
43  use m_geometry,           only : xred2xcart
44  use m_occ,                only : newocc
45  use m_dtset,              only : testsusmat
46  use m_pawang,             only : pawang_type
47  use m_pawtab,             only : pawtab_type
48  use m_paw_ij,             only : paw_ij_type
49  use m_pawfgrtab,          only : pawfgrtab_type
50  use m_pawrhoij,           only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_io, pawrhoij_get_nspden
51  use m_pawcprj,            only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_getdim
52  use m_pawfgr,             only : pawfgr_type
53  use m_energies,           only : energies_type
54  use m_hamiltonian,        only : init_hamiltonian,destroy_hamiltonian, &
55                                   load_spin_hamiltonian, load_k_hamiltonian, gs_hamiltonian_type
56  use m_bandfft_kpt,        only : bandfft_kpt, bandfft_kpt_type, bandfft_kpt_set_ikpt, &
57                                   bandfft_kpt_savetabs, bandfft_kpt_restoretabs, prep_bandfft_tabs
58  use m_electronpositron,   only : electronpositron_type,electronpositron_calctype
59  use m_paw_dmft,           only : paw_dmft_type,init_dmft,destroy_dmft,print_dmft,saveocc_dmft
60  use m_paw_correlations,   only : setnoccmmp
61  use m_paw_occupancies,   only : pawmkrhoij
62  use m_paw_mkrho,          only : pawmkrho
63  use m_crystal,            only : crystal_init, crystal_free, crystal_t
64  use m_oper,               only : oper_type,init_oper,destroy_oper
65  use m_io_tools,           only : flush_unit
66  use m_abi2big,            only : wvl_occ_abi2big, wvl_rho_abi2big, wvl_occopt_abi2big, wvl_eigen_abi2big
67  use m_fock,               only : fock_type, fock_ACE_type, fock_updateikpt, fock_calc_ene
68  use m_invovl,             only : make_invovl
69  use m_tddft,              only : tddft
70  use m_kg,                 only : mkkin, mkkpg, mknucdipmom_k
71  use m_suscep_stat,        only : suscep_stat
72  use m_fft,                only : fftpac
73  use m_spacepar,           only : symrhg
74  use m_vtowfk,             only : vtowfk
75  use m_mkrho,              only : mkrho, prtrhomxmn
76  use m_mkffnl,             only : mkffnl
77  use m_mpinfo,             only : proc_distrb_cycle
78  use m_common,             only : prteigrs
79  use m_dmft,               only : dmft_solve
80  use m_datafordmft,        only : datafordmft
81  use m_fourier_interpol,   only : transgrid
82  use m_cgprj,              only : ctocprj
83  use m_wvl_rho,            only : wvl_mkrho
84  use m_wvl_psi,            only : wvl_hpsitopsi, wvl_psitohpsi, wvl_nl_gradient
85 
86 #if defined HAVE_BIGDFT
87  use BigDFT_API,           only : last_orthon, evaltoocc, write_energies, eigensystem_info
88 #endif
89 
90  implicit none
91 
92  private

ABINIT/vtorho [ Functions ]

[ Top ] [ Functions ]

NAME

 vtorho

FUNCTION

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

INPUTS

  afford=used to dimension susmat
  atindx(natom)=index table for atoms (see gstate.f)
  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cpus= cpu time limit in seconds
  dbl_nnsclo=if 1, will double the value of dtset%nnsclo
  dielop= if positive, the dielectric matrix must be computed.
  dielstrt=number of the step at which the dielectric preconditioning begins.
  dmatpawu= fixed occupation matrix of correlated orbitals (DFT+U or DMFT only)
  dtefield <type(efield_type)> = variables related to Berry phase
      calculations (see initberry.f)
  dtfil <type(datafiles_type)>=variables related to files
  dtset <type(dataset_type)>=all input variables for this dataset
   | mband=maximum number of bands
   | mgfft=maximum size of 1D FFTs
   | mkmem =number of k points treated by this node.
   | mpw=maximum dimensioned size of npw
   | nfft=(effective) number of FFT grid points (for this processor)
   | nkpt=number of k points.
   | nspden=number of spin-density components
   | nsppol=1 for unpolarized, 2 for spin-polarized
   | nsym=number of symmetry elements in space group
   | typat= array of types of the natoms
  electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
  etotal=total energy (Ha) - only needed for tddft
  fock <type(fock_type)>= quantities to calculate Fock exact exchange
  gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  gprimd(3,3)=dimensional reciprocal space primitive translations
   (3x3 tensor) and grads wrt atomic coordinates (3*natom)
  gsqcut=cutoff on (k+G)^2 (bohr^-2)
  hdr <type(hdr_type)>=the header of wf, den and pot files
  indsym(4,nsym,natom)=indirect indexing array for atom labels
  irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
  irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data for diel matrix
                                     nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  istep=index of the number of steps in the routine scfcv
  istep_mix=index of the number of steps for the SCF mixing (can be <istep)
  kg(3,mpw*mkmem)=reduced planewave coordinates.
  kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
  kxc(nfftf,nkxc)=exchange-correlation kernel, needed only if nkxc/=0 .
  lmax_diel=1+max. value of l angular momentum used for dielectric matrix
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
  mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
  mpi_enreg=informations about MPI parallelization
  my_natom=number of atoms treated by current processor
  natom=number of atoms in cell.
  nattyp(ntypat)= # atoms of each type.
  nfftf= -PAW ONLY- number of FFT grid points for the fine grid
         (nfftf=nfft for norm-conserving potential runs)
  nfftdiel=number of fft grid points for the computation of the diel matrix
  ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix,
                see ~abinit/doc/variables/vargs.htm#ngfft
  nkxc=second dimension of the array kxc, see rhotoxc.f for a description
  npwarr(nkpt)=number of planewaves in basis at this k point
  npwdiel=size of the susmat array.
  ntypat=number of types of atoms in unit cell.
  optforces=option for the computation of forces (0: no force;1: forces)
  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_dmft  <type(paw_dmft_type)>= paw+dmft related data
  paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  pawang <type(pawang)>=paw angular mesh and related data
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  pawfgrtab(my_natom*usepaw) <type(pawfgrtab_type)>=atomic data given on fine rectangular grid
  pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
  phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
                                    nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  phnonsdiel(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases,
   for diel matr
                                     nfft**(1-1/nsym) is 1 if nsym==1, and nfft otherwise
  ph1d(2,3*(2*mgfft+1)*natom)=one-dimensional structure factor information
  ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information
                                             for the dielectric matrix
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
  pwind_alloc = first dimension of pwind
  pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
                           (see initberry.f)
  rmet(3,3)=real space metric (bohr**2)
  rprimd(3,3)=dimensional primitive vectors
  symrec(3,3,nsym)=symmetry operations in reciprocal space
  ucvol=unit cell volume in bohr**3.
  usecprj=1 if cprj datastructure is stored in memory
  wffnew,unit numbers for wf disk files.
  vtrial(nfftf,nspden)=INPUT potential Vtrial(r).
  vxctau=(only for meta-GGA): derivative of XC energy density with respect to
    kinetic energy density (depsxcdtau). The arrays vxctau(nfft,nspden,4) contains also
    the gradient of vxctau (gvxctau) in vxctau(:,:,2:4)
  xred(3,natom)=reduced dimensionless atomic coordinates
  ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
  ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
  ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point
                                 for the dielectric matrix

OUTPUT

  compch_fft=-PAW only- compensation charge inside spheres computed over fine fft grid
  dphase(3) : dphase(idir) = accumulated change in the string-averaged
     Zak phase along the idir-th direction caused by the update of all
     the occupied Bloch states at all the k-points (only if finite electric field)
  eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
  resid(mband*nkpt*nsppol)=residuals for each band over all k points.
  residm=maximum value from resid array (except for nbdbuf highest bands)
  susmat(2,npwdiel*afford,nspden,npwdiel,nspden)=
   the susceptibility (or density-density response) matrix in reciprocal space
  === if optforces>0 ===
    grnl(3*natom)=stores grads of nonlocal energy wrt length scales
  ==== if optres==1
    nres2=square of the norm of the residual
    nvresid(nfftf,nspden)=density residual
  ==== if psps%usepaw==1
    cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
                               cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.
    nhat(nfftf,nspden*psps%usepaw)=compensation charge density on rectangular grid in real space

SIDE EFFECTS

  cg(2,mpw*dtset%nspinor*mband*mkmem*nsppol)=planewave coefficients of wavefunctions.
   At output contains updated wavefunctions coefficients;
    if nkpt>1, these are kept in a disk file.
  energies <type(energies_type)>=storage for energies computed here :
   | e_eigenvalues=Sum of the eigenvalues - Band energy (Hartree)
   | e_kinetic=kinetic energy part of total energy
   | e_nonlocalpsp=nonlocal pseudopotential part of total energy
   | e_fermie=fermi energy (Hartree)
  occ(mband*nkpt*nsppol)=occupation number for each band for each k.
      (input if insulator - occopt<3 - ; output if metallic)
  pawrhoij(my_natom*usepaw) <type(pawrhoij_type)>= paw rhoij occupancies and related data
  rhog(2,nfftf)=Fourier transform of total electron density
  rhor(nfftf,nspden)=total electron density (el/bohr**3)
  taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
  taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
  tauresid(nfftf,nspden*dtset%usekden)=array for kinetic energy density residual
  wvl <type(wvl_data)>=wavelets structures in case of wavelets basis.
  ==== if (usepaw==1) ====
    cprj(natom,mcprj*usecprj)= wave functions projected with non-local projectors:
                               cprj(n,k,i)=<p_i|Cnk> where p_i is a non-local projector.

PARENTS

      scfcv

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

NOTES

  Be careful to the meaning of nfft (size of FFT grids):
   - In case of norm-conserving calculations the FFT grid is the usual FFT grid.
   - In case of PAW calculations:
     Two FFT grids are used; one with nfft points (coarse grid) for
     the computation of wave functions ; one with nfftf points
     (fine grid) for the computation of total density.

  The total electronic density (rhor,rhog) is divided into two terms:
   - The density related to WFs =Sum[Psi**2]
   - The compensation density (nhat) - only in PAW

  The parallelisation needed for the electric field should be
  made an independent subroutine, so that this routine could be put
  back in the 95_drive directory.

SOURCE

 273 subroutine vtorho(afford,atindx,atindx1,cg,compch_fft,cprj,cpus,dbl_nnsclo,&
 274 &           dielop,dielstrt,dmatpawu,dphase,dtefield,dtfil,dtset,&
 275 &           eigen,electronpositron,energies,etotal,gbound_diel,&
 276 &           gmet,gprimd,grnl,gsqcut,hdr,indsym,irrzon,irrzondiel,&
 277 &           istep,istep_mix,kg,kg_diel,kxc,lmax_diel,mcg,mcprj,mgfftdiel,mpi_enreg,&
 278 &           my_natom,natom,nattyp,nfftf,nfftdiel,ngfftdiel,nhat,nkxc,&
 279 &           npwarr,npwdiel,nres2,ntypat,nvresid,occ,optforces,&
 280 &           optres,paw_dmft,paw_ij,pawang,pawfgr,pawfgrtab,pawrhoij,pawtab,&
 281 &           phnons,phnonsdiel,ph1d,ph1ddiel,psps,fock,&
 282 &           pwind,pwind_alloc,pwnsfac,resid,residm,rhog,rhor,&
 283 &           rmet,rprimd,susmat,symrec,taug,taur,tauresid,&
 284 &           ucvol,usecprj,wffnew,vtrial,vxctau,wvl,xred,ylm,ylmgr,ylmdiel)
 285 
 286 
 287 !This section has been created automatically by the script Abilint (TD).
 288 !Do not modify the following lines by hand.
 289 #undef ABI_FUNC
 290 #define ABI_FUNC 'vtorho'
 291 !End of the abilint section
 292 
 293  implicit none
 294 
 295 !Arguments -------------------------------
 296  integer, intent(in) :: afford,dbl_nnsclo,dielop,dielstrt,istep,istep_mix,lmax_diel,mcg,mcprj,mgfftdiel
 297  integer, intent(in) :: my_natom,natom,nfftf,nfftdiel,nkxc,npwdiel
 298  integer, intent(in) :: ntypat,optforces,optres,pwind_alloc,usecprj
 299  real(dp), intent(in) :: cpus,etotal,gsqcut,ucvol
 300  real(dp), intent(out) :: compch_fft,nres2,residm
 301  type(MPI_type), intent(inout) :: mpi_enreg
 302  type(datafiles_type), intent(in) :: dtfil
 303  type(dataset_type), intent(inout) :: dtset
 304  type(efield_type), intent(inout) :: dtefield
 305  type(electronpositron_type),pointer :: electronpositron
 306  type(energies_type), intent(inout) :: energies
 307  type(hdr_type), intent(inout) :: hdr
 308  type(paw_dmft_type), intent(inout)  :: paw_dmft
 309  type(pawang_type), intent(in) :: pawang
 310  type(pawfgr_type), intent(in) :: pawfgr
 311  type(pseudopotential_type), intent(in) :: psps
 312  type(fock_type),pointer, intent(inout) :: fock
 313  type(wffile_type), intent(inout) :: wffnew
 314  type(wvl_data), intent(inout) :: wvl
 315  integer, intent(in) :: atindx(natom),atindx1(natom),gbound_diel(2*mgfftdiel+8,2)
 316  integer, intent(in) :: indsym(4,dtset%nsym,natom)
 317  integer, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 318  integer, intent(in) :: irrzondiel(nfftdiel**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 319  integer, intent(in) :: kg(3,dtset%mpw*dtset%mkmem),kg_diel(3,npwdiel),nattyp(ntypat),ngfftdiel(18),npwarr(dtset%nkpt)
 320  integer, intent(in) :: pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
 321  real(dp), intent(in) :: dmatpawu(:,:,:,:),gmet(3,3),gprimd(3,3),ph1d(2,3*(2*dtset%mgfft+1)*natom)
 322  real(dp), intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*psps%usepaw)
 323  real(dp), intent(in) :: phnons(2,dtset%nfft**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 324  real(dp), intent(in) :: phnonsdiel(2,nfftdiel**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))
 325  real(dp), intent(in) :: pwnsfac(2,pwind_alloc),rmet(3,3),rprimd(3,3)
 326  real(dp), intent(inout) :: vtrial(nfftf,dtset%nspden)
 327  real(dp), intent(inout) :: xred(3,natom)
 328  real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
 329  real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
 330  real(dp), intent(in) :: ylmdiel(npwdiel,lmax_diel**2)
 331  real(dp), intent(out) :: dphase(3),grnl(3*natom)
 332  real(dp), intent(inout) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
 333  real(dp), intent(out) :: nhat(nfftf,dtset%nspden*psps%usepaw)
 334  real(dp), intent(inout) :: resid(dtset%mband*dtset%nkpt*dtset%nsppol)
 335  real(dp), intent(out) :: nvresid(nfftf,dtset%nspden)
 336  real(dp), intent(out) :: susmat(2,npwdiel*afford,dtset%nspden,npwdiel,dtset%nspden)
 337  real(dp), intent(inout) :: cg(2,mcg)
 338  real(dp), intent(inout) :: kxc(nfftf,nkxc),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
 339  real(dp), intent(inout) :: rhog(2,nfftf),rhor(nfftf,dtset%nspden)
 340  real(dp), intent(inout) :: taug(2,nfftf*dtset%usekden),taur(nfftf,dtset%nspden*dtset%usekden)
 341  real(dp), intent(inout) :: tauresid(nfftf,dtset%nspden*dtset%usekden)
 342  real(dp), intent(inout),optional :: vxctau(nfftf,dtset%nspden,4*dtset%usekden)
 343  type(pawcprj_type),pointer,intent(inout) :: cprj(:,:)
 344  type(paw_ij_type),intent(inout) :: paw_ij(my_natom*psps%usepaw)
 345  type(pawfgrtab_type),intent(inout) :: pawfgrtab(my_natom*psps%usepaw)
 346  type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom*psps%usepaw)
 347  type(pawtab_type),intent(in) :: pawtab(ntypat*psps%usepaw)
 348 
 349 !Local variables-------------------------------
 350 !scalars
 351  integer,parameter :: level=111,tim_mkrho=2
 352  integer,save :: nwarning=0
 353  integer :: bdtot_index,counter,cplex,dimffnl,enunit,iband,iband1,ibdkpt
 354  integer :: ibg,icg,ider,idir,ierr,ifft,ifor,ifor1,ii,ikg,ikpt
 355  integer :: ikpt_loc,ikpt1,my_ikpt,ikxc,ilm,imagn,index1,iorder_cprj,ipert,iplex
 356  integer :: iscf,ispden,isppol,istwf_k,mband_cprj,mbdkpsp,mb2dkpsp
 357  integer :: mcgq,mcprj_local,mcprj_tmp,me_distrb,mkgq,mpi_comm_sphgrid
 358  integer :: mwarning,my_nspinor,n1,n2,n3,n4,n5,n6,nband_eff
 359  integer :: nband_k,nband_cprj_k,nbuf,neglect_pawhat,nfftot,nkpg,nkpt1,nnn,nnsclo_now
 360  integer :: nproc_distrb,npw_k,nspden_rhoij,option,prtvol
 361  integer :: spaceComm_distrb,usecprj_local,usefock_ACE,usetimerev
 362  logical :: berryflag,computesusmat,fixed_occ
 363  logical :: locc_test,paral_atom,remove_inv,usefock,with_vxctau
 364  logical :: do_last_ortho,wvlbigdft=.false.
 365  real(dp) :: dmft_ldaocc
 366  real(dp) :: edmft,ebandlda,ebanddmft,ebandldatot,ekindmft,ekindmft2,ekinlda
 367  real(dp) :: min_occ,vxcavg_dum,strsxc(6)
 368  character(len=500) :: message
 369  type(bandfft_kpt_type),pointer :: my_bandfft_kpt => null()
 370  type(gs_hamiltonian_type) :: gs_hamk
 371 !arrays
 372  integer,allocatable :: kg_k(:,:)
 373  real(dp) :: dielar(7),dphase_k(3),kpoint(3),qpt(3),rhodum(1),tsec(2),ylmgr_dum(0,0,0)
 374  real(dp),allocatable :: EigMin(:,:),buffer1(:),buffer2(:),cgq(:,:)
 375  real(dp),allocatable :: cgrkxc(:,:),cgrvtrial(:,:),doccde(:)
 376  real(dp),allocatable :: dphasek(:,:),eig_k(:),ek_k(:),ek_k_nd(:,:,:),eknk(:),eknk_nd(:,:,:,:,:)
 377  real(dp),allocatable :: enl_k(:),enlnk(:),focknk(:),fockfornk(:,:,:),ffnl(:,:,:,:),grnl_k(:,:), xcart(:,:)
 378  real(dp),allocatable :: grnlnk(:,:),kinpw(:),kpg_k(:,:),occ_k(:),ph3d(:,:,:)
 379  real(dp),allocatable :: pwnsfacq(:,:),resid_k(:),rhoaug(:,:,:,:),rhowfg(:,:),rhowfr(:,:)
 380  real(dp),allocatable :: vlocal(:,:,:,:),vlocal_tmp(:,:,:),vxctaulocal(:,:,:,:,:),ylm_k(:,:),zshift(:)
 381  complex(dpc),target,allocatable :: nucdipmom_k(:)
 382  type(pawcprj_type),allocatable :: cprj_tmp(:,:)
 383  type(pawcprj_type),allocatable,target:: cprj_local(:,:)
 384  type(oper_type) :: lda_occup
 385  type(pawrhoij_type),pointer :: pawrhoij_unsym(:)
 386 
 387  type(crystal_t) :: cryst_struc
 388  integer :: idum1(0),idum3(0,0,0)
 389  real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
 390 
 391 !Variables for BigDFT
 392 #if defined HAVE_BIGDFT
 393  integer :: occopt_bigdft
 394 #endif
 395 
 396 ! *********************************************************************
 397 
 398  DBG_ENTER("COLL")
 399 
 400 !Keep track of total time spent in vtorho
 401  call timab(980,1,tsec)
 402  call timab(981,1,tsec)
 403 
 404 !Structured debugging if prtvol==-level
 405  prtvol=dtset%prtvol
 406 
 407 !If usewvl: wvlbigdft indicates that the BigDFT workflow will be followed
 408  wvlbigdft=(dtset%usewvl==1.and.dtset%wvl_bigdft_comp==1)
 409 
 410 !Several inits
 411  n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
 412  n4=dtset%ngfft(4) ; n5=dtset%ngfft(5) ; n6=dtset%ngfft(6)
 413  usecprj_local=0;if (psps%usepaw==1) usecprj_local=1
 414  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 415  paral_atom=(my_natom/=natom)
 416  compch_fft=-1.d5
 417 
 418 !Check that usekden is not 0 if want to use vxctau
 419  with_vxctau = (dtset%usekden/=0)
 420 
 421 !Check that fock is present if want to use fock option
 422  usefock = (dtset%usefock==1 .and. associated(fock))
 423  usefock_ACE=0
 424  if (usefock) usefock_ACE=fock%fock_common%use_ACE
 425 
 426 
 427 !Init MPI
 428  spaceComm_distrb=mpi_enreg%comm_cell
 429  if (mpi_enreg%paral_kgb==1) spaceComm_distrb=mpi_enreg%comm_kpt
 430  if (mpi_enreg%paral_hf ==1) spaceComm_distrb=mpi_enreg%comm_kpt
 431  nproc_distrb=xmpi_comm_size(spaceComm_distrb)
 432  me_distrb=xmpi_comm_rank(spaceComm_distrb)
 433  mpi_comm_sphgrid=mpi_enreg%comm_fft
 434  if(dtset%usewvl==1) mpi_comm_sphgrid=mpi_enreg%comm_wvl
 435  if (mpi_enreg%me_img/=0) nwarning=nwarning+1
 436 
 437 !Test size of FFT grids (1 grid in norm-conserving, 2 grids in PAW)
 438  if ((psps%usepaw==1.and.pawfgr%nfft/=nfftf).or.(psps%usepaw==0.and.dtset%nfft/=nfftf)) then
 439    MSG_BUG('wrong values for nfft, nfftf!')
 440  end if
 441 
 442 !Test optforces (to prevent memory overflow)
 443  if (optforces/=0.and.optforces/=1) then
 444    write(message,'(a,i0)')' wrong value for optforces = ',optforces
 445    MSG_BUG(message)
 446  end if
 447 
 448  iscf=dtset%iscf
 449  fixed_occ=(dtset%occopt<3.or.electronpositron_calctype(electronpositron)==1)
 450  if(.not. wvlbigdft) then
 451    energies%e_eigenvalues = zero
 452    energies%e_kinetic     = zero
 453    energies%e_nonlocalpsp = zero
 454    if (usefock) then
 455      energies%e_fock=zero
 456      energies%e_fockdc=zero
 457    end if
 458    grnl(:)=zero
 459    resid(:) = zero ! JWZ 13 May 2010. resid and eigen need to be fully zeroed each time before use
 460    eigen(:) = zero
 461    bdtot_index=0
 462    ibg=0;icg=0
 463    mbdkpsp=dtset%mband*dtset%nkpt*dtset%nsppol
 464    if(paw_dmft%use_dmft==1) mb2dkpsp=2*dtset%mband*dtset%mband*dtset%nkpt*dtset%nsppol
 465  end if
 466 
 467  if(dtset%usewvl==0) then
 468    ABI_ALLOCATE(eknk,(mbdkpsp))
 469    ABI_ALLOCATE(enlnk,(mbdkpsp))
 470    ABI_ALLOCATE(eknk_nd,(dtset%nsppol,dtset%nkpt,2,dtset%mband,dtset%mband*paw_dmft%use_dmft))
 471    ABI_ALLOCATE(EigMin,(2,dtset%mband))
 472    ABI_ALLOCATE(grnlnk,(3*natom,mbdkpsp*optforces))
 473    if (usefock) then
 474      ABI_ALLOCATE(focknk,(mbdkpsp))
 475      focknk=zero
 476      if (optforces>0)then
 477        ABI_ALLOCATE(fockfornk,(3,natom,mbdkpsp))
 478        fockfornk=zero
 479      end if
 480    end if
 481    eknk(:)=zero;enlnk(:)=zero
 482    if (optforces>0) grnlnk(:,:)=zero
 483    if(paw_dmft%use_dmft==1) eknk_nd=zero
 484  end if !usewvl==0
 485 
 486 !Initialize rhor if needed; store old rhor
 487  if(iscf>=0 .or. iscf==-3) then
 488    if (optres==1) then
 489      nvresid=rhor
 490      tauresid=taur
 491    end if
 492 !  NC and plane waves
 493    if (psps%usepaw==0 .and. dtset%usewvl==0) then
 494      rhor=zero
 495 !    PAW
 496    elseif(psps%usepaw==1) then
 497      ABI_ALLOCATE(rhowfr,(dtset%nfft,dtset%nspden))
 498      ABI_ALLOCATE(rhowfg,(2,dtset%nfft))
 499      rhowfr(:,:)=zero
 500    end if
 501  end if
 502 
 503 !Set max number of non-self-consistent loops nnsclo_now for use in vtowfk
 504  if(iscf<0)then
 505    ! ===== Non self-consistent =====
 506    nnsclo_now=dtset%nstep
 507  else
 508    ! ===== Self-consistent =====
 509    if(dtset%nnsclo>0) then
 510    ! ===== Self-consistent + imposed =====
 511      nnsclo_now=dtset%nnsclo
 512    else
 513    ! ===== Self-consistent + default =====
 514      nnsclo_now=1
 515      if (dtset%usewvl==0) then
 516      ! ----- Plane waves -----
 517        if (istep<=2.and.iscf/=0) nnsclo_now=2
 518      else
 519      ! ----- Wavelets -----
 520        if (iscf==0) then
 521          nnsclo_now=0
 522        else if (istep<=2) then
 523          nnsclo_now=3
 524        else if (istep<=4) then
 525          nnsclo_now=2
 526        end if
 527      end if
 528    end if
 529    ! ===== Double is required =====
 530    if(dbl_nnsclo==1)then
 531      nnsclo_now=nnsclo_now*2
 532    end if
 533  end if
 534  if(dtset%wfoptalg==2)nnsclo_now=40  ! UNDER DEVELOPMENT
 535 
 536  write(message, '(a,i0,a,3(i0,1x))' ) ' vtorho : nnsclo_now=',nnsclo_now,&
 537 & ', note that nnsclo,dbl_nnsclo,istep=',dtset%nnsclo,dbl_nnsclo,istep
 538  call wrtout(std_out,message,'COLL')
 539 
 540 !==== Initialize most of the Hamiltonian ====
 541 !Allocate all arrays and initialize quantities that do not depend on k and spin.
 542  call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,natom,&
 543 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,&
 544 & paw_ij=paw_ij,ph1d=ph1d,usecprj=usecprj_local,electronpositron=electronpositron,fock=fock,&
 545 & comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab,mpi_spintab=mpi_enreg%my_isppoltab,&
 546 & nucdipmom=dtset%nucdipmom,use_gpu_cuda=dtset%use_gpu_cuda)
 547 
 548 !Initializations for PAW (projected wave functions)
 549  mcprj_local=0 ; mband_cprj=0
 550  if (psps%usepaw==1) then
 551    mband_cprj=dtset%mband
 552    if (dtset%paral_kgb/=0) mband_cprj=mband_cprj/mpi_enreg%nproc_band
 553    iorder_cprj=0 ; mcprj_local=mcprj
 554    if (usecprj==0) then
 555      mcprj_local=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
 556      !This is a check but should always be true since scfcv allocated cprj anyway
 557      if (allocated(cprj_local)) then
 558        !Was allocated in scfcv so we just destroy and reconstruct it as desired
 559        call pawcprj_free(cprj_local)
 560        ABI_DATATYPE_DEALLOCATE(cprj_local)
 561      end if
 562      ABI_DATATYPE_ALLOCATE(cprj_local,(dtset%natom,mcprj_local))
 563      call pawcprj_alloc(cprj_local,0,gs_hamk%dimcprj)
 564      cprj=> cprj_local
 565    end if
 566  end if
 567 
 568  call timab(981,2,tsec)
 569 
 570 !===================================================================
 571 ! WAVELETS - Branching with a separate VTORHO procedure
 572 !===================================================================
 573 
 574  if (dtset%usewvl == 1) then
 575 #ifndef HAVE_BIGDFT
 576    BIGDFT_NOTENABLED_ERROR()
 577 #else
 578 
 579 !  do_last_ortho in case of diagonalization scheme
 580    if (     wvlbigdft) do_last_ortho=(dtset%iscf/=0)
 581    if (.not.wvlbigdft) do_last_ortho=(.true.)
 582 
 583    ABI_ALLOCATE(xcart,(3, dtset%natom))
 584    call xred2xcart(dtset%natom, rprimd, xcart, xred)
 585 
 586    if(wvlbigdft) then
 587 !    NSCF loop for wvlbigdt:
 588      call wvl_nscf_loop_bigdft()
 589    else
 590 !    NSCF loop for WVL: (not wvlbigdft)
 591      call wvl_nscf_loop()
 592    end if
 593 
 594 !  Eventually orthogonalize WFs now
 595    if (do_last_ortho) then
 596      call write_energies(ii,0,wvl%e%energs,0.d0,0.d0,"final")
 597      call last_orthon(me_distrb, nproc_distrb, ii, wvl%wfs%ks, wvl%e%energs%evsum, .true.)
 598      if(wvlbigdft) energies%e_xcdc = wvl%e%energs%evxc
 599 !    If occupation numbers are not changed...
 600      if (fixed_occ .or. (iscf<0 .and. iscf/=-3)) then
 601        call wvl_comm_eigen()
 602      end if
 603 !    ... or update occupations:
 604      if( ( .not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
 605        if(wvlbigdft) then
 606          call wvl_occ_bigdft()
 607        else
 608 !        Communicate eigenvalues:
 609          call wvl_comm_eigen()
 610 !        Update occ and Fermi level
 611          call wvl_occ()
 612        end if
 613      end if
 614 !    This might accelerate convergence:
 615      wvl%wfs%ks%diis%energy_min=one
 616      wvl%wfs%ks%diis%alpha=two
 617    end if !do_last_ortho
 618 
 619 !  Compute eigenvalues energy
 620    if(.not. wvlbigdft .and. nnsclo_now>0) then
 621      call e_eigen(eigen,energies%e_eigenvalues,dtset%mband,dtset%nband,dtset%nkpt,&
 622 &     dtset%nsppol,occ,dtset%wtk)
 623    else
 624      energies%e_eigenvalues = energies%e_kinetic + energies%e_localpsp &
 625 &     + energies%e_xcdc  + two*energies%e_hartree +energies%e_nonlocalpsp
 626    end if
 627 
 628    if (optforces == 1) then ! not compatible with iscf=0 and wvlbigdftcomp=1 + PAW
 629      call wvl_nl_gradient(grnl, mpi_enreg, dtset%natom, rprimd, wvl, xcart)
 630    end if
 631 
 632 !  For iscf<0 we do not update the density
 633    if (dtset%iscf>=0) then !(dtset%iscf>=0 .and. .not. wvlbigdft ) then
 634      call wvl_mkrho(dtset,irrzon,mpi_enreg,phnons,rhor,wvl%wfs,wvl%den)
 635    end if
 636    ABI_DEALLOCATE(xcart)
 637 
 638 !  Note in WVL+NC: the rest will be skipped.
 639 !  For PAW: we will compute Rho_ij at the end.
 640    !if(wvlbigdft) return
 641 #endif
 642  else
 643 
 644 !===================================================================
 645 ! PLANE WAVES - Standard VTORHO procedure
 646 !===================================================================
 647 
 648 !  Electric fields: set flag to turn on various behaviors
 649    berryflag = (dtset%berryopt == 4 .or. dtset%berryopt == 14 .or. &
 650 &   dtset%berryopt == 6 .or. dtset%berryopt == 16 .or. &
 651 &   dtset%berryopt == 7 .or. dtset%berryopt == 17)
 652 
 653 !  Electric field: allocate dphasek
 654    nkpt1 = dtset%nkpt
 655    if ( berryflag ) then
 656      ABI_ALLOCATE(dphasek,(3,dtset%nkpt*dtset%nsppol))
 657      dphasek(:,:) = zero
 658      nkpt1 = dtefield%mkmem_max
 659    end if
 660 
 661    ABI_ALLOCATE(rhoaug,(n4,n5,n6,gs_hamk%nvloc))
 662    ABI_ALLOCATE(vlocal,(n4,n5,n6,gs_hamk%nvloc))
 663    if(with_vxctau) then
 664      ABI_ALLOCATE(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
 665    end if
 666 
 667 !  LOOP OVER SPINS
 668    do isppol=1,dtset%nsppol
 669      call timab(982,1,tsec)
 670 
 671      ikpt_loc = 0
 672      ikg=0
 673 
 674 !    Set up local potential vlocal with proper dimensioning, from vtrial
 675 !    Also take into account the spin.
 676      if(dtset%nspden/=4)then
 677        if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then
 678          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal,2)
 679          if(with_vxctau) then
 680            do ii=1,4
 681              call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,&
 682 &             vxctau(:,:,ii),vxctaulocal(:,:,:,:,ii),2)
 683            end do
 684          end if
 685        else
 686          ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
 687          call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
 688          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal,2)
 689          ABI_DEALLOCATE(cgrvtrial)
 690        end if
 691      else
 692        ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6))
 693        if (psps%usepaw==0.or.pawfgr%usefinegrid==0) then
 694          do ispden=1,dtset%nspden
 695            call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,vtrial,vlocal_tmp,2)
 696            vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
 697          end do
 698        else
 699          ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
 700          call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
 701          do ispden=1,dtset%nspden
 702            call fftpac(ispden,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,cgrvtrial,vlocal_tmp,2)
 703            vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
 704          end do
 705          ABI_DEALLOCATE(cgrvtrial)
 706        end if
 707        ABI_DEALLOCATE(vlocal_tmp)
 708      end if ! nspden
 709 
 710      rhoaug(:,:,:,:)=zero
 711 
 712 !    Continue to initialize the Hamiltonian
 713      call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
 714      if (with_vxctau) then
 715        call load_spin_hamiltonian(gs_hamk,isppol,vxctaulocal=vxctaulocal)
 716      end if
 717 
 718      call timab(982,2,tsec)
 719 
 720 !    BIG FAT k POINT LOOP
 721 !    MVeithen: I had to modify the structure of this loop in order to implement MPI // of the electric field
 722 !    Note that the loop here differs from the similar one in berryphase_new.F90.
 723 !    here, ikpt_loc numbers the kpts treated by the current processor.
 724 !    in berryphase_new.F90, ikpt_loc ALSO includes info about value of isppol.
 725 
 726      ikpt = 0
 727      do while (ikpt_loc < nkpt1)
 728 
 729        call timab(997,1,tsec)
 730 
 731        if ( .not.berryflag ) then
 732          ikpt_loc = ikpt_loc + 1
 733          ikpt = ikpt_loc
 734          my_ikpt = mpi_enreg%my_kpttab(ikpt)
 735        else
 736          if (ikpt_loc < dtset%mkmem) ikpt = ikpt + 1
 737          if ((ikpt > dtset%nkpt).and.(ikpt_loc < dtset%mkmem)) exit
 738          my_ikpt=ikpt
 739        end if
 740 
 741        dphase_k(:) = zero
 742        counter=100*ikpt+isppol
 743        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
 744        nband_cprj_k=nband_k/mpi_enreg%nproc_band
 745        istwf_k=dtset%istwfk(ikpt)
 746        npw_k=npwarr(ikpt)
 747 
 748        mcgq=1 ; mkgq=1
 749        if (.not.berryflag) then
 750          if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) then
 751            eigen(1+bdtot_index : nband_k+bdtot_index) = zero
 752            resid(1+bdtot_index : nband_k+bdtot_index) = zero
 753            bdtot_index=bdtot_index+nband_k
 754            cycle
 755          end if
 756        else
 757          if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_distrb)) &
 758 &         .and.(ikpt_loc <= dtset%mkmem)) then
 759            eigen(1+bdtot_index : nband_k+bdtot_index) = zero
 760            resid(1+bdtot_index : nband_k+bdtot_index) = zero
 761            bdtot_index = bdtot_index + nband_k
 762            cycle
 763          end if
 764          ikpt_loc = ikpt_loc + 1
 765          mcgq = dtset%mpw*my_nspinor*nband_k*dtefield%nneigh(ikpt)
 766          ikg = dtefield%kgindex(ikpt)
 767          mkgq = 6*dtset%mpw
 768        end if ! berryflag
 769 
 770        call timab(997,2,tsec)
 771 
 772 !      In case of MPI // of a finite field calculation
 773 !      build the cgq array that stores the wavefunctions for the
 774 !      neighbours of ikpt, and the pwnsfacq array that stores the
 775 !      corresponding phase factors (in case of tnons)
 776        ABI_ALLOCATE(cgq,(2,mcgq))
 777        ABI_ALLOCATE(pwnsfacq,(2,mkgq))
 778        if ( berryflag ) then
 779          call cgq_builder(berryflag,cg,cgq,dtefield,dtset,ikpt,ikpt_loc,isppol,mcg,mcgq,&
 780 &         me_distrb,mkgq,mpi_enreg,my_nspinor,nband_k,nproc_distrb,&
 781 &         npwarr,pwnsfac,pwnsfacq,pwind_alloc,spaceComm_distrb)
 782          if (ikpt_loc > dtset%mkmem) then
 783            ABI_DEALLOCATE(cgq)
 784            ABI_DEALLOCATE(pwnsfacq)
 785            cycle
 786          end if
 787        end if !berryopt
 788 
 789        call timab(984,1,tsec)
 790 
 791        if (mpi_enreg%paral_kgb==1) my_bandfft_kpt => bandfft_kpt(my_ikpt)
 792        call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
 793 !      my_ikpt = ikpt
 794 !       if (mpi_enreg%paral_kgb==1) then
 795 !        my_ikpt = mpi_enreg%my_kpttab(ikpt)
 796 !         my_bandfft_kpt => bandfft_kpt(my_ikpt)
 797 !         call bandfft_kpt_set_ikpt(ikpt,mpi_enreg)
 798 !       end if
 799 
 800        ABI_ALLOCATE(eig_k,(nband_k))
 801        ABI_ALLOCATE(ek_k,(nband_k))
 802        ABI_ALLOCATE(enl_k,(nband_k))
 803        ABI_ALLOCATE(ek_k_nd,(2,nband_k,nband_k*paw_dmft%use_dmft))
 804        ABI_ALLOCATE(occ_k,(nband_k))
 805        ABI_ALLOCATE(resid_k,(nband_k))
 806        ABI_ALLOCATE(zshift,(nband_k))
 807        ABI_ALLOCATE(grnl_k,(3*natom,nband_k*optforces))
 808 
 809        eig_k(:)=zero
 810        ek_k(:)=zero
 811        enl_k(:)=zero
 812        if(paw_dmft%use_dmft==1) ek_k_nd(:,:,:)=zero
 813        if (optforces>0) grnl_k(:,:)=zero
 814        kpoint(:)=dtset%kptns(:,ikpt)
 815        occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index)
 816        resid_k(:)=zero
 817        zshift(:)=dtset%eshift
 818 
 819        ABI_ALLOCATE(kg_k,(3,npw_k))
 820        ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
 821        kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
 822        if (psps%useylm==1) then
 823          do ilm=1,psps%mpsang*psps%mpsang
 824            ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
 825          end do
 826        end if
 827 
 828 !      Set up remaining of the Hamiltonian
 829 
 830 !      Compute (1/2) (2 Pi)**2 (k+G)**2:
 831        ABI_ALLOCATE(kinpw,(npw_k))
 832 !       call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k)
 833        call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
 834 
 835 !      Compute (k+G) vectors (only if useylm=1)
 836        nkpg=3*optforces*dtset%nloalg(3)
 837        ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
 838        if ((mpi_enreg%paral_kgb/=1.or.istep<=1).and.nkpg>0) then
 839          call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
 840        end if
 841 
 842 !      Compute nonlocal form factors ffnl at all (k+G):
 843        ider=0;idir=0;dimffnl=1
 844        ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,ntypat))
 845        if (mpi_enreg%paral_kgb/=1.or.istep<=1) then
 846          call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
 847 &         gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
 848 &         psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
 849 &         npw_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,&
 850 &         psps%usepaw,psps%useylm,ylm_k,ylmgr)
 851        end if
 852 
 853 !     compute and load nuclear dipole Hamiltonian at current k point
 854        if(any(abs(gs_hamk%nucdipmom)>0.0)) then
 855          if(allocated(nucdipmom_k)) then
 856            ABI_DEALLOCATE(nucdipmom_k)
 857          end if
 858          ABI_ALLOCATE(nucdipmom_k,(npw_k*(npw_k+1)/2))
 859          call mknucdipmom_k(gmet,kg_k,kpoint,natom,gs_hamk%nucdipmom,nucdipmom_k,npw_k,rprimd,ucvol,xred)
 860          if(allocated(gs_hamk%nucdipmom_k)) then
 861             ABI_DEALLOCATE(gs_hamk%nucdipmom_k)
 862          end if
 863          ABI_ALLOCATE(gs_hamk%nucdipmom_k,(npw_k*(npw_k+1)/2))
 864          call load_k_hamiltonian(gs_hamk,nucdipmom_k=nucdipmom_k)
 865        end if
 866 
 867 
 868 !      Load k-dependent part in the Hamiltonian datastructure
 869 !       - Compute 3D phase factors
 870 !       - Prepare various tabs in case of band-FFT parallelism
 871 !       - Load k-dependent quantities in the Hamiltonian
 872        ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
 873 
 874        if(usefock_ACE/=0) then
 875          call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
 876 &         kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,fockACE_k=fock%fockACE(ikpt,isppol),ph3d_k=ph3d,&
 877 &         compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),&
 878 &         compute_gbound=(mpi_enreg%paral_kgb/=1))
 879        else
 880          call load_k_hamiltonian(gs_hamk,kpt_k=dtset%kptns(:,ikpt),istwf_k=istwf_k,npw_k=npw_k,&
 881 &         kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
 882 &         compute_ph3d=(mpi_enreg%paral_kgb/=1.or.istep<=1),&
 883 &         compute_gbound=(mpi_enreg%paral_kgb/=1))
 884        end if
 885 
 886 !      Load band-FFT tabs (transposed k-dependent arrays)
 887        if (mpi_enreg%paral_kgb==1) then
 888          if (istep<=1) then
 889            call prep_bandfft_tabs(gs_hamk,ikpt,dtset%mkmem,mpi_enreg)
 890          end if
 891          call load_k_hamiltonian(gs_hamk,npw_fft_k=my_bandfft_kpt%ndatarecv, &
 892 &         gbound_k =my_bandfft_kpt%gbound, &
 893 &         kinpw_k  =my_bandfft_kpt%kinpw_gather, &
 894 &         kg_k     =my_bandfft_kpt%kg_k_gather, &
 895 &         kpg_k    =my_bandfft_kpt%kpg_k_gather, &
 896          ffnl_k   =my_bandfft_kpt%ffnl_gather, &
 897          ph3d_k   =my_bandfft_kpt%ph3d_gather)
 898        end if
 899 
 900 !      Build inverse of overlap matrix for chebfi
 901        if(psps%usepaw == 1 .and. dtset%wfoptalg == 1 .and. istep <= 1) then
 902          call make_invovl(gs_hamk, dimffnl, ffnl, ph3d, mpi_enreg)
 903        end if
 904 
 905        ! Setup gemm_nonlop
 906        if (gemm_nonlop_use_gemm) then
 907          !set the global variable indicating to gemm_nonlop where to get its data from
 908          gemm_nonlop_ikpt_this_proc_being_treated = my_ikpt
 909          if (istep <= 1) then
 910            !Init the arrays
 911            call make_gemm_nonlop(my_ikpt,gs_hamk%npw_fft_k,gs_hamk%lmnmax, &
 912 &           gs_hamk%ntypat, gs_hamk%indlmn, gs_hamk%nattyp, gs_hamk%istwf_k, gs_hamk%ucvol, gs_hamk%ffnl_k,&
 913 &           gs_hamk%ph3d_k)
 914          end if
 915        end if
 916 
 917 #if defined HAVE_GPU_CUDA
 918        if (dtset%use_gpu_cuda==1) then
 919          call gpu_update_ffnl_ph3d(ph3d,size(ph3d),ffnl,size(ffnl))
 920        end if
 921 #endif
 922 
 923        call timab(984,2,tsec)
 924 
 925 !      Update the value of ikpt,isppol in fock_exchange and allocate the memory space to perform HF calculation.
 926        if (usefock) then
 927          call fock_updateikpt(fock%fock_common,ikpt,isppol)
 928        end if
 929        if ((psps%usepaw==1).and.(usefock)) then
 930          if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then
 931            fock%fock_common%forces_ikpt=zero
 932          end if
 933        end if
 934 
 935 !      Compute the eigenvalues, wavefunction, residuals,
 936 !      contributions to kinetic energy, nonlocal energy, forces,
 937 !      and update of rhor to this k-point and this spin polarization.
 938        call vtowfk(cg,cgq,cprj,cpus,dphase_k,dtefield,dtfil,&
 939 &       dtset,eig_k,ek_k,ek_k_nd,enl_k,fixed_occ,grnl_k,gs_hamk,&
 940 &       ibg,icg,ikpt,iscf,isppol,kg_k,kinpw,mband_cprj,mcg,mcgq,mcprj_local,mkgq,&
 941 &       mpi_enreg,dtset%mpw,natom,nband_k,dtset%nkpt,nnsclo_now,npw_k,npwarr,&
 942 &       occ_k,optforces,prtvol,pwind,pwind_alloc,pwnsfac,pwnsfacq,resid_k,&
 943 &       rhoaug,paw_dmft,dtset%wtk(ikpt),zshift)
 944 
 945        call timab(985,1,tsec)
 946 
 947 #if defined HAVE_GPU_CUDA
 948        if(dtset%use_gpu_cuda==1) then
 949          call gpu_finalize_ffnl_ph3d()
 950        end if
 951 #endif
 952        ABI_DEALLOCATE(ffnl)
 953        ABI_DEALLOCATE(kinpw)
 954        ABI_DEALLOCATE(kg_k)
 955        ABI_DEALLOCATE(kpg_k)
 956        if(allocated(nucdipmom_k)) then
 957          ABI_DEALLOCATE(nucdipmom_k)
 958        end if
 959        ABI_DEALLOCATE(ylm_k)
 960        ABI_DEALLOCATE(ph3d)
 961        ABI_DEALLOCATE(cgq)
 962        ABI_DEALLOCATE(pwnsfacq)
 963 
 964 !      electric field
 965        if (berryflag) then
 966          dphasek(:,ikpt + (isppol - 1)*dtset%nkpt) = dphase_k(:)
 967 
 968 !        The overlap matrices for all first neighbours of ikpt are no more up to date
 969          do idir = 1, 3
 970            do ifor = 1, 2
 971              ikpt1 = dtefield%ikpt_dk(dtefield%i2fbz(ikpt),ifor,idir)
 972              ikpt1 = dtefield%indkk_f2ibz(ikpt1,1)
 973              ifor1 = -1*ifor + 3   ! ifor = 1 -> ifor1 = 2 & ifor = 2 -> ifor1 = 1
 974              dtefield%sflag(:,ikpt1+(isppol-1)*dtset%nkpt,ifor1,idir) = 0
 975            end do
 976          end do
 977        end if  ! berryflag
 978 
 979 !      Save eigenvalues (hartree), residuals (hartree**2)
 980        eigen(1+bdtot_index : nband_k+bdtot_index) = eig_k(:)
 981        eknk (1+bdtot_index : nband_k+bdtot_index) = ek_k (:)
 982        if(usefock) then
 983          focknk (1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%eigen_ikpt (:)
 984          if (optforces>0) fockfornk(:,:,1+bdtot_index : nband_k+bdtot_index) = fock%fock_common%forces_ikpt(:,:,:)
 985        end if
 986        if(paw_dmft%use_dmft==1) eknk_nd(isppol,ikpt,:,:,:) = ek_k_nd(:,:,:)
 987        resid(1+bdtot_index : nband_k+bdtot_index) = resid_k(:)
 988        if (optforces>0) grnlnk(:,1+bdtot_index : nband_k+bdtot_index) = grnl_k(:,:)
 989        enlnk(1+bdtot_index : nband_k+bdtot_index) = enl_k(:)
 990 
 991        if(iscf>0 .or. iscf==-3)then
 992 !        Accumulate sum over k points for band, nonlocal and kinetic energies,
 993 !        also accumulate gradients of Enonlocal:
 994          do iband=1,nband_k
 995            if (abs(occ_k(iband))>tol8) then
 996              energies%e_kinetic     = energies%e_kinetic     + dtset%wtk(ikpt)*occ_k(iband)*ek_k(iband)
 997              energies%e_eigenvalues = energies%e_eigenvalues + dtset%wtk(ikpt)*occ_k(iband)*eig_k(iband)
 998              energies%e_nonlocalpsp = energies%e_nonlocalpsp + dtset%wtk(ikpt)*occ_k(iband)*enl_k(iband)
 999              if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ_k(iband)*grnl_k(:,iband)
1000              if (usefock) energies%e_fock=energies%e_fock + half*fock%fock_common%eigen_ikpt(iband)*occ_k(iband)*dtset%wtk(ikpt)
1001            end if
1002          end do
1003 
1004 !        Calculate Fock contribution to the total energy if required
1005          if ((psps%usepaw==1).and.(usefock)) then
1006            if ((fock%fock_common%optfor).and.(usefock_ACE==0)) then
1007              call fock_calc_ene(dtset,fock%fock_common,energies%e_exactX,ikpt,nband_k,occ_k)
1008            end if
1009          end if
1010        end if
1011 
1012        call timab(985,2,tsec)
1013 
1014        ABI_DEALLOCATE(eig_k)
1015        ABI_DEALLOCATE(ek_k)
1016        ABI_DEALLOCATE(ek_k_nd)
1017        ABI_DEALLOCATE(grnl_k)
1018        ABI_DEALLOCATE(occ_k)
1019        ABI_DEALLOCATE(resid_k)
1020        ABI_DEALLOCATE(zshift)
1021        ABI_DEALLOCATE(enl_k)
1022 
1023 !      Keep track of total number of bands (all k points so far, even for k points not treated by me)
1024        bdtot_index=bdtot_index+nband_k
1025 
1026 !      Also shift array memory if dtset%mkmem/=0
1027        if (dtset%mkmem/=0) then
1028          ibg=ibg+my_nspinor*nband_cprj_k
1029          icg=icg+npw_k*my_nspinor*nband_k
1030          ikg=ikg+npw_k
1031        end if
1032 
1033      end do ! End big k point loop
1034 
1035      call timab(986,1,tsec)
1036 
1037      if (fixed_occ .and. mpi_enreg%paral_kgb==1) then
1038        call xmpi_sum(rhoaug,mpi_enreg%comm_bandspinorfft,ierr) !Sum the contributions over bands/FFT/spinors
1039      end if
1040 
1041 !    Transfer density on augmented fft grid to normal fft grid in real space
1042 !    Also take into account the spin.
1043      if(iscf>0.or.iscf==-3)then
1044        if (psps%usepaw==0) then
1045          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,1),1)
1046          if(dtset%nspden==4)then
1047            do imagn=2,4
1048              call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhor,rhoaug(:,:,:,imagn),1)
1049            end do
1050          end if
1051        else
1052          call fftpac(isppol,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,1),1)
1053          if(dtset%nspden==4)then
1054            do imagn=2,4
1055              call fftpac(imagn,mpi_enreg,dtset%nspden,n1,n2,n3,n4,n5,n6,dtset%ngfft,rhowfr,rhoaug(:,:,:,imagn),1)
1056            end do
1057          end if
1058        end if
1059      end if
1060 
1061      call timab(986,2,tsec)
1062 
1063    end do ! End loop over spins
1064 
1065    call timab(988,1,tsec)
1066    if (usefock) then
1067      if(fock%fock_common%optfor) then
1068        call xmpi_sum(fock%fock_common%forces,mpi_enreg%comm_kpt,ierr)
1069      end if
1070    end if
1071 !  Electric field: compute string-averaged change in Zak phase
1072 !  along each direction, store it in dphase(idir)
1073 !  ji: it is not convenient to do this anymore. Remove. Set dphase(idir)=0.0_dp.
1074 !  eventually, dphase(idir) will have to go...
1075    if (berryflag) then
1076      dphase(:) = zero
1077 !    In case of MPI // of a finite field calculation, send dphasek to all cpus
1078      call xmpi_sum(dphasek,spaceComm_distrb,ierr)
1079      ABI_DEALLOCATE(dphasek)
1080    end if ! berryflag
1081    ABI_DEALLOCATE(rhoaug)
1082    ABI_DEALLOCATE(vlocal)
1083    if(with_vxctau) then
1084      ABI_DEALLOCATE(vxctaulocal)
1085    end if
1086 
1087    call timab(988,2,tsec)
1088 
1089    ABI_ALLOCATE(doccde,(dtset%mband*dtset%nkpt*dtset%nsppol))
1090    doccde(:)=zero
1091 
1092 !  Treat now varying occupation numbers, in the self-consistent case
1093    if((.not.fixed_occ) .and. (iscf>0.or.iscf==-3)) then
1094 
1095 !    Parallel case
1096      if (mpi_enreg%nproc_kpt>1) then
1097 
1098        call timab(989,1,tsec)
1099 
1100 !      If needed, exchange the values of eigen,resid,eknk,enlnk,grnlnk
1101        ABI_ALLOCATE(buffer1,((4+3*natom*optforces+dtset%usefock+3*natom*dtset%usefock*optforces)*mbdkpsp))
1102        if(paw_dmft%use_dmft==1) then
1103          ABI_ALLOCATE(buffer2,(mb2dkpsp*paw_dmft%use_dmft))
1104        end if
1105 !      Pack eigen,resid,eknk,enlnk,grnlnk in buffer1
1106        buffer1(1          :  mbdkpsp)=eigen(:)
1107        buffer1(1+  mbdkpsp:2*mbdkpsp)=resid(:)
1108        buffer1(1+2*mbdkpsp:3*mbdkpsp)=eknk(:)
1109        buffer1(1+3*mbdkpsp:4*mbdkpsp)=enlnk(:)
1110        index1=4*mbdkpsp
1111        if (optforces>0) then
1112          buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(grnlnk,(/(3*natom)*mbdkpsp/) )
1113          index1=index1+3*natom*mbdkpsp
1114        end if
1115        if (usefock) then
1116          buffer1(1+index1:index1+mbdkpsp)=focknk(:)
1117          if (optforces>0) then
1118            index1=index1+mbdkpsp
1119            buffer1(index1+1:index1+3*natom*mbdkpsp)=reshape(fockfornk,(/(3*natom)*mbdkpsp/) )
1120          end if
1121        end if
1122        if(paw_dmft%use_dmft==1) then
1123          nnn=0
1124          do ikpt=1,dtset%nkpt
1125            do isppol=1,dtset%nsppol
1126              do iband=1,dtset%mband
1127                do iband1=1,dtset%mband
1128                  do iplex=1,2
1129                    nnn=nnn+1
1130                    buffer2(nnn)=eknk_nd(isppol,ikpt,iplex,iband,iband1)
1131                  end do
1132                end do
1133              end do
1134            end do
1135          end do
1136          if(nnn.ne.mb2dkpsp)  then
1137            write(message,*)' BUG in vtorho2, buffer2',nnn,mb2dkpsp
1138            MSG_BUG(message)
1139          end if
1140        end if
1141 !      Build sum of everything
1142        call timab(48,1,tsec)
1143        call xmpi_sum(buffer1,mpi_enreg%comm_kpt,ierr)
1144        if(mpi_enreg%paral_kgb/=1.and.paw_dmft%use_dmft==1) then
1145          call xmpi_sum(buffer2,mpi_enreg%comm_kpt,ierr)
1146        end if
1147        call timab(48,2,tsec)
1148 
1149 !      Unpack eigen,resid,eknk,enlnk,grnlnk in buffer1
1150        eigen(:) =buffer1(1          :  mbdkpsp)
1151        resid(:) =buffer1(1+  mbdkpsp:2*mbdkpsp)
1152        eknk(:)  =buffer1(1+2*mbdkpsp:3*mbdkpsp)
1153        enlnk(:) =buffer1(1+3*mbdkpsp:4*mbdkpsp)
1154        index1=4*mbdkpsp
1155        if (optforces>0) then
1156          grnlnk(:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3*natom,mbdkpsp/) )
1157        end if
1158        if (usefock) then
1159          focknk(:)=buffer1(1+index1:index1+mbdkpsp)
1160          if (optforces>0) then
1161            index1=index1+mbdkpsp
1162            fockfornk(:,:,:)=reshape(buffer1(index1+1:index1+3*natom*mbdkpsp),(/3,natom,mbdkpsp/) )
1163          end if
1164        end if
1165        if(paw_dmft%use_dmft==1) then
1166          nnn=0
1167          do ikpt=1,dtset%nkpt
1168            do isppol=1,dtset%nsppol
1169              do iband=1,dtset%mband
1170                do iband1=1,dtset%mband
1171                  do iplex=1,2
1172                    nnn=nnn+1
1173                    eknk_nd(isppol,ikpt,iplex,iband,iband1)=buffer2(nnn)
1174                  end do
1175                end do
1176              end do
1177            end do
1178          end do
1179        end if
1180        if(allocated(buffer2))  then
1181          ABI_DEALLOCATE(buffer2)
1182        end if
1183        ABI_DEALLOCATE(buffer1)
1184        call timab(989,2,tsec)
1185 
1186      end if ! nproc_kpt>1
1187 
1188 !    Compute the new occupation numbers from eigen
1189      call timab(990,1,tsec)
1190      call newocc(doccde,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,&
1191 &     dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,&
1192 &     dtset%nsppol,occ,dtset%occopt,prtvol,dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk)
1193      call timab(990,2,tsec)
1194 
1195 !    !=========  DMFT call begin ============================================
1196      dmft_ldaocc=0
1197      if(paw_dmft%use_dmft==1.and.psps%usepaw==1.and.dtset%nbandkss==0) then
1198        call timab(991,1,tsec)
1199 
1200        if (dtset%dmftcheck>=0.and.dtset%usedmft>=1.and.(sum(pawtab(:)%upawu)>=tol8.or.  &
1201 &       sum(pawtab(:)%jpawu)>tol8).and.dtset%dmft_entropy==0) energies%entropy=zero
1202 
1203 !      ==  0 to a dmft calculation and do not use lda occupations
1204 !      ==  1 to a lda calculation with the dmft loop
1205        if(dtset%dmftcheck==-1) dmft_ldaocc=1
1206 
1207 !      ==  initialise occnd
1208        paw_dmft%occnd=zero
1209 
1210        if(dmft_ldaocc==0) then
1211          if(dtset%occopt/=3) then
1212            message = ' occopt should be equal to 3 in dmft'
1213            MSG_ERROR(message)
1214          end if
1215 !        ==  initialise edmft
1216          if(paw_dmft%use_dmft>=1) edmft = zero
1217 
1218          !  Compute residm to check the value
1219          ibdkpt=1
1220          residm=zero
1221          do isppol=1,dtset%nsppol
1222            do ikpt=1,dtset%nkpt
1223              nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1224              nband_eff=max(1,nband_k-dtset%nbdbuf)
1225              residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
1226              ibdkpt=ibdkpt+nband_k
1227            end do
1228          end do
1229 
1230          !  Test residm
1231          if (paw_dmft%use_dmft>0 .and. residm>tol4 .and. dtset%dmftcheck>=0) then
1232            if(dtset%dmft_entropy>0)  then
1233              write(message,'(a,e12.3)')&
1234 &             ' WARNING: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm
1235              call wrtout(std_out,message,'COLL')
1236            else
1237              write(message,'(a,e12.3)')&
1238 &             ' ERROR: Wavefunctions not converged : DFT+DMFT calculation cannot be carried out safely ',residm
1239              call wrtout(std_out,message,'COLL')
1240              write(message,'(a,i0)')'  Action: increase nline and nnsclo',dtset%nstep
1241              MSG_ERROR(message)
1242            end if
1243 
1244          else if (paw_dmft%use_dmft>0 .and. residm>tol10.and. dtset%dmftcheck>=0) then
1245            write(message,'(3a)')ch10,&
1246 &           '  Wavefunctions not converged : DFT+DMFT calculation might not be carried out safely ',ch10
1247            MSG_WARNING(message)
1248          end if
1249 
1250 !        ==  allocate paw_dmft%psichi and paw_dmft%eigen_lda
1251          call init_dmft(dmatpawu,dtset,energies%e_fermie,dtfil%fnameabo_app,dtset%nspinor,paw_dmft,pawtab,psps,dtset%typat)
1252          call print_dmft(paw_dmft,dtset%pawprtvol)
1253 
1254 !        ==  gather crystal structure date into data "cryst_struc"
1255          remove_inv=.false.
1256          if(dtset%nspden==4) remove_inv=.true.
1257          call crystal_init(dtset%amu_orig(:,1),cryst_struc,dtset%spgroup,natom,dtset%npsp,ntypat, &
1258 &         dtset%nsym,rprimd,dtset%typat,xred,dtset%ziontypat,dtset%znucl,1,&
1259 &         dtset%nspden==2.and.dtset%nsppol==1,remove_inv,hdr%title,&
1260 &         dtset%symrel,dtset%tnons,dtset%symafm)
1261 
1262 !        ==  compute psichi
1263          call xmpi_barrier(spaceComm_distrb)
1264          call init_oper(paw_dmft,lda_occup)
1265          call flush_unit(std_out)
1266          call timab(620,1,tsec)
1267          call datafordmft(cryst_struc,cprj,gs_hamk%dimcprj,dtset,eigen,energies%e_fermie,&
1268 &         lda_occup,dtset%mband,mband_cprj,dtset%mkmem,mpi_enreg,&
1269 &         dtset%nkpt,my_nspinor,dtset%nsppol,occ,&
1270 &         paw_dmft,paw_ij,pawang,pawtab,psps,usecprj_local,dtfil%unpaw)
1271          call timab(620,2,tsec)
1272          call flush_unit(std_out)
1273 
1274 !        ==  solve dmft loop
1275          call xmpi_barrier(spaceComm_distrb)
1276 
1277          call dmft_solve(cryst_struc,istep,lda_occup,paw_dmft,pawang,pawtab,dtset%pawprtvol)
1278          edmft=paw_dmft%edmft
1279          energies%e_paw=energies%e_paw+edmft
1280          energies%e_pawdc=energies%e_pawdc+edmft
1281          call flush_unit(std_out)
1282 !        paw_dmft%occnd(:,:,:,:,:)=0.5_dp
1283 
1284 !        call print_dmft(paw_dmft,dtset%pawprtvol)
1285          if(dtset%paral_kgb==1) then
1286            write(message,'(5a)')ch10,&
1287 &           ' Parallelization over bands is not yet compatible with self-consistency in DMFT ',ch10,&
1288 &           ' Calculation of density does not taken into account non diagonal occupations',ch10
1289            call wrtout(std_out,message,'COLL')
1290            call wrtout(ab_out,message,'COLL')
1291 !          MSG_ERROR(message)
1292            if(dtset%nstep>1) then
1293              write(message,'(a,i0)')'  Action: use nstep=1 instead of nstep=',dtset%nstep
1294              MSG_ERROR(message)
1295            end if
1296            residm=zero
1297          end if
1298          if(dtset%nspinor==2) then
1299 !          call flush_unit(ab_out)
1300 !          write(message,'(3a)')&
1301 !          &         ' Self consistent DFT+DMFT with nspinor==2 is not possible yet ',ch10,&
1302 !          &         ' Calculation are restricted to nstep =1'
1303 !          !         MSG_ERROR(message)
1304 !          if(dtset%nstep>1) then
1305 !          write(message,'(a,i0)')' Action: use nstep=1 instead of nstep=',dtset%nstep
1306 !          !           MSG_ERROR(message)
1307 !          endif
1308          end if
1309 
1310          if(me_distrb==0) then
1311            call saveocc_dmft(paw_dmft)
1312          end if
1313          call destroy_dmft(paw_dmft)
1314 
1315 !        ==  destroy crystal_t cryst_struc
1316          call crystal_free(cryst_struc)
1317          call destroy_oper(lda_occup)
1318        end if ! dmft_ldaocc
1319        call timab(991,2,tsec)
1320 
1321      end if ! usedmft
1322 
1323      if(dtset%nbandkss/=0) then
1324        write(message,'(a,i3,2a,i3,4a)') &
1325 &       " dtset%nbandkss ==",dtset%nbandkss,ch10,&
1326 &       " and dtset%usedmft ==",dtset%usedmft,ch10,&
1327 &       " a DFT loop is carried out without DMFT.",ch10,&
1328 &       " Only psichi's will be written at convergence of the DFT loop."
1329        call wrtout(std_out,message,'COLL')
1330      end if
1331 !    !=========  DMFT call end   ============================================
1332 
1333      call timab(992,1,tsec)
1334 
1335 !    Compute eeig, ek,enl and grnl from the new occ, and the shared eknk,enlnk,grnlnk
1336      energies%e_eigenvalues = zero
1337      energies%e_kinetic     = zero
1338      energies%e_nonlocalpsp = zero
1339      if (usefock) then
1340        energies%e_fock     = zero
1341        if (optforces>0) fock%fock_common%forces=zero
1342      end if
1343      if (optforces>0) grnl(:)=zero
1344      if(paw_dmft%use_dmft>=1) then
1345        ebandlda               = zero
1346        ebanddmft              = zero
1347        ebandldatot            = zero
1348        ekindmft               = zero
1349        ekindmft2              = zero
1350        ekinlda                = zero
1351      end if
1352 
1353 !    Compute new energy terms due to non diagonal occupations and DMFT.
1354      bdtot_index=1
1355      do isppol=1,dtset%nsppol
1356        do ikpt=1,dtset%nkpt
1357          nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1358          do iband=1,nband_k
1359 
1360            locc_test = abs(occ(bdtot_index))>tol8
1361 !          dmft
1362            if(paw_dmft%use_dmft>=1.and.dtset%nbandkss==0) then
1363              if(paw_dmft%band_in(iband)) then
1364                if( paw_dmft%use_dmft == 1 .and. dmft_ldaocc == 1 ) then ! test of the code
1365                  paw_dmft%occnd(1,iband,iband,ikpt,isppol)=occ(bdtot_index)
1366                end if
1367                locc_test = abs(paw_dmft%occnd(1,iband,iband,ikpt,isppol))+&
1368 &               abs(paw_dmft%occnd(2,iband,iband,ikpt,isppol))>tol8
1369              end if
1370            end if
1371 
1372            if (locc_test) then
1373 !            dmft
1374              if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then
1375                ebandldatot=ebandldatot+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1376                if(paw_dmft%band_in(iband)) then
1377                  ebandlda=ebandlda+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1378                  ekinlda=ekinlda+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1379                  occ(bdtot_index)=paw_dmft%occnd(1,iband,iband,ikpt,isppol)
1380                  ebanddmft=ebanddmft+dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1381                  ekindmft=ekindmft+dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1382                end if
1383              end if
1384 
1385              energies%e_eigenvalues = energies%e_eigenvalues + &
1386 &             dtset%wtk(ikpt)*occ(bdtot_index)*eigen(bdtot_index)
1387              energies%e_kinetic = energies%e_kinetic + &
1388 &             dtset%wtk(ikpt)*occ(bdtot_index)*eknk(bdtot_index)
1389              energies%e_nonlocalpsp = energies%e_nonlocalpsp + &
1390 &             dtset%wtk(ikpt)*occ(bdtot_index)*enlnk(bdtot_index)
1391              if (usefock) then
1392                energies%e_fock=energies%e_fock + half*focknk(bdtot_index)*occ(bdtot_index)*dtset%wtk(ikpt)
1393                if (optforces>0) fock%fock_common%forces(:,:)=fock%fock_common%forces(:,:)+&
1394 &               dtset%wtk(ikpt)*occ(bdtot_index)*fockfornk(:,:,bdtot_index)
1395              end if
1396              if (optforces>0) grnl(:)=grnl(:)+dtset%wtk(ikpt)*occ(bdtot_index)*grnlnk(:,bdtot_index)
1397            end if
1398            bdtot_index=bdtot_index+1
1399            if(paw_dmft%use_dmft==1.and.dtset%nbandkss==0) then
1400              do iband1=1,nband_k
1401                if(paw_dmft%band_in(iband).and.paw_dmft%band_in(iband1)) then
1402 !                write(std_out,*) "II+", isppol,ikpt,iband,iband1
1403                  ekindmft2=ekindmft2  +  dtset%wtk(ikpt)*paw_dmft%occnd(1,iband,iband1,ikpt,isppol)*&
1404 &                 eknk_nd(isppol,ikpt,1,iband,iband1)
1405                  if(dtset%nspinor==2) then
1406                    ekindmft2=ekindmft2  -  dtset%wtk(ikpt)*paw_dmft%occnd(2,iband,iband1,ikpt,isppol)*&
1407 &                   eknk_nd(isppol,ikpt,2,iband,iband1)
1408                  end if
1409 !                write(std_out,*) "II", occnd(1,iband,iband1,ikpt,isppol),eknk_nd(isppol,ikpt,iband,iband1)
1410                end if
1411              end do
1412            end if
1413          end do
1414        end do
1415      end do
1416 
1417      if(paw_dmft%use_dmft==1) then
1418        energies%e_kinetic = energies%e_kinetic -ekindmft+ekindmft2
1419        if(abs(dtset%pawprtvol)>=2) then
1420          write(message,'(4a,7(2x,a,2x,e14.7,a),a)') &
1421 &         "-----------------------------------------------",ch10,&
1422 &         "--- Energy for DMFT and tests (in Ha)  ",ch10,&
1423 &         "--- Ebandldatot    (Ha.) = ",ebandldatot,ch10,&
1424 &         "--- Ebandlda       (Ha.) = ",ebandlda,ch10,&
1425 &         "--- Ebanddmft      (Ha.) = ",ebanddmft,ch10,&
1426 &         "--- Ekinlda        (Ha.) = ",ekinlda,ch10, &
1427 &         "--- Ekindmftdiag   (Ha.) = ",ekindmft,ch10,&
1428 &         "--- Ekindmftnondiag(Ha.) = ",ekindmft2,ch10,&
1429 &         "--- Edmft=         (Ha.) = ",edmft,ch10,&
1430 &         "-----------------------------------------------"
1431          call wrtout(std_out,message,'COLL')
1432        end if
1433        if(paw_dmft%use_dmft==1.and.mpi_enreg%paral_kgb==1) paw_dmft%use_dmft=0
1434      end if
1435 
1436      if (psps%usepaw==0) then
1437        call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1438 &       rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1439      else
1440        call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1441 &       rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs)
1442      end if
1443      call timab(992,2,tsec)
1444 
1445 !    Treat fixed occupation numbers or non-self-consistent case
1446    else
1447 
1448      if (mpi_enreg%nproc_kpt>1) then
1449 
1450        call timab(989,1,tsec)
1451 
1452        nbuf=2*mbdkpsp+dtset%nfft*dtset%nspden+3+3*natom*optforces
1453 !      * If Hartree-Fock calculation, the exact exchange energy is k-dependent.
1454        if(dtset%usefock==1) then
1455          nbuf=nbuf+1
1456          if (optforces>0) nbuf=nbuf+3*natom
1457        end if
1458        if(iscf==-1 .or. iscf==-2)nbuf=2*mbdkpsp
1459        ABI_ALLOCATE(buffer1,(nbuf))
1460 !      Pack eigen,resid,rho[wf]r,grnl,enl,ek
1461        buffer1(1:mbdkpsp)=eigen(:)
1462        buffer1(1+mbdkpsp:2*mbdkpsp)=resid(:)
1463        index1=2*mbdkpsp
1464        if(iscf/=-1 .and. iscf/=-2)then
1465          if (psps%usepaw==0) then
1466            buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhor,(/dtset%nfft*dtset%nspden/))
1467          else
1468            buffer1(index1+1:index1+dtset%nfft*dtset%nspden)=reshape(rhowfr,(/dtset%nfft*dtset%nspden/))
1469          end if
1470          index1=index1+dtset%nfft*dtset%nspden
1471          buffer1(index1+1) = energies%e_kinetic
1472          buffer1(index1+2) = energies%e_eigenvalues
1473          buffer1(index1+3) = energies%e_nonlocalpsp
1474          index1=index1+3
1475 !        * If Hartree-Fock calculation, save e_fock in buffer1
1476          if (dtset%usefock==1) then
1477            buffer1(index1+1) = energies%e_fock
1478            index1=index1+1
1479            if (optforces>0)then
1480              buffer1(index1+1:index1+3*natom)=reshape(fock%fock_common%forces,(/3*natom/))
1481              index1=index1+3*natom
1482            end if
1483          end if
1484          if (optforces>0) buffer1(index1+1:index1+3*natom)=grnl(1:3*natom)
1485        end if
1486 
1487 !      Build sum of everything
1488        call timab(48,1,tsec)
1489        call xmpi_sum(buffer1,nbuf,mpi_enreg%comm_kpt ,ierr)
1490        call timab(48,2,tsec)
1491 
1492 !      Unpack the final result
1493        eigen(:)=buffer1(1:mbdkpsp)
1494        resid(:)=buffer1(1+mbdkpsp:2*mbdkpsp)
1495        index1=2*mbdkpsp
1496        if(iscf/=-1 .and. iscf/=-2)then
1497          if (psps%usepaw==0) then
1498            ii=1
1499            do ispden=1,dtset%nspden
1500              do ifft=1,dtset%nfft
1501                rhor(ifft,ispden)=buffer1(index1+ii)
1502                ii=ii+1
1503              end do
1504            end do
1505          else
1506            ii=1
1507            do ispden=1,dtset%nspden
1508              do ifft=1,dtset%nfft
1509                rhowfr(ifft,ispden)=buffer1(index1+ii)
1510                ii=ii+1
1511              end do
1512            end do
1513          end if
1514          index1=index1+dtset%nfft*dtset%nspden
1515          energies%e_kinetic = buffer1(index1+1)
1516          energies%e_eigenvalues = buffer1(index1+2)
1517          energies%e_nonlocalpsp = buffer1(index1+3)
1518          index1=index1+3
1519 !        * If Hartree-Fock calculation, save e_fock in buffer1
1520          if (dtset%usefock==1) then
1521            energies%e_fock = buffer1(index1+1)
1522            index1=index1+1
1523            if (optforces>0) then
1524              fock%fock_common%forces(:,:)=reshape(buffer1(index1+1:index1+3*natom),(/3,natom/))
1525              index1=index1+3*natom
1526            end if
1527          end if
1528          if (optforces>0) grnl(1:3*natom)=buffer1(index1+1:index1+3*natom)
1529        end if
1530        ABI_DEALLOCATE(buffer1)
1531        call timab(989,2,tsec)
1532 
1533      end if ! nproc_kpt>1
1534 
1535 
1536 !    Compute the highest occupied eigenenergy
1537      if(iscf/=-1 .and. iscf/=-2)then
1538        call timab(993,1,tsec)
1539        energies%e_fermie = -huge(one)
1540        bdtot_index=1
1541        do isppol=1,dtset%nsppol
1542          do ikpt=1,dtset%nkpt
1543            nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1544            do iband=1,nband_k
1545              if(abs(occ(bdtot_index))>tol8 .and. eigen(bdtot_index)>energies%e_fermie+tol10) then
1546                energies%e_fermie=eigen(bdtot_index)
1547              end if
1548              bdtot_index=bdtot_index+1
1549            end do
1550          end do
1551        end do
1552        call xmpi_max(energies%e_fermie,spaceComm_distrb,ierr)
1553        call timab(993,2,tsec)
1554      end if
1555 
1556 !    If needed, compute rhog, and symmetrizes the density
1557      if (iscf > 0 .or. iscf==-3 ) then
1558 !      energies%e_fermie=zero  ! Actually, should determine the maximum of the valence band XG20020802
1559        nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
1560 
1561        call timab(994,1,tsec)
1562        if (psps%usepaw==0) then
1563          call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
1564 &         dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhog  ,rhor  ,rprimd,dtset%symafm,dtset%symrel)
1565        else
1566          call symrhg(1,gprimd,irrzon,mpi_enreg,dtset%nfft,nfftot,dtset%ngfft,dtset%nspden,&
1567 &         dtset%nsppol,dtset%nsym,dtset%paral_kgb,phnons,rhowfg,rhowfr,rprimd,dtset%symafm,dtset%symrel)
1568        end if
1569        call timab(994,2,tsec)
1570 !      We now have both rho(r) and rho(G), symmetrized, and if dtset%nsppol=2
1571 !      we also have the spin-up density, symmetrized, in rhor(:,2).
1572      end if
1573 
1574    end if !  End of test on varying or fixed occupation numbers
1575 
1576    call timab(994,1,tsec)
1577 
1578 
1579 !  Compute the kinetic energy density
1580    if(dtset%usekden==1 .and. (iscf > 0 .or. iscf==-3 ) )then
1581      call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,mpi_enreg,npwarr,occ,paw_dmft,phnons,&
1582 &     taug,taur,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs,option=1)
1583    end if
1584 
1585    ABI_DEALLOCATE(eknk)
1586    if (usefock) then
1587      ABI_DEALLOCATE(focknk)
1588      if (optforces>0)then
1589        ABI_DEALLOCATE(fockfornk)
1590      end if
1591    end if
1592    ABI_DEALLOCATE(eknk_nd)
1593    ABI_DEALLOCATE(grnlnk)
1594    ABI_DEALLOCATE(enlnk)
1595 
1596 !  In the non-self-consistent case, print eigenvalues and residuals
1597    if(iscf<=0 .and. me_distrb==0)then
1598      option=2 ; enunit=1 ; vxcavg_dum=zero
1599      call prteigrs(eigen,enunit,energies%e_fermie,dtfil%fnameabo_app_eig,&
1600 &     ab_out,iscf,dtset%kptns,dtset%kptopt,dtset%mband,dtset%nband,&
1601 &     dtset%nkpt,nnsclo_now,dtset%nsppol,occ,dtset%occopt,option,&
1602 &     dtset%prteig,prtvol,resid,dtset%tolwfr,vxcavg_dum,dtset%wtk)
1603    end if
1604 
1605 !  Find largest residual over bands, k points, and spins, except for nbdbuf highest bands
1606    ibdkpt=1
1607    residm=zero
1608    do isppol=1,dtset%nsppol
1609      do ikpt=1,dtset%nkpt
1610        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1611        nband_eff=max(1,nband_k-dtset%nbdbuf)
1612        residm=max(residm,maxval(resid(ibdkpt:ibdkpt+nband_eff-1)))
1613        ibdkpt=ibdkpt+nband_k
1614      end do
1615    end do
1616 
1617  end if !usewvl==0
1618 
1619 !===================================================================
1620 ! End of PLANE WAVES section
1621 !===================================================================
1622 
1623 !In the self-consistent case, diagnose lack of unoccupied state (for each spin and k-point).
1624 !Print a warning if the number of such messages already written does not exceed mwarning.
1625  mwarning=5
1626  if(nwarning<mwarning .and. iscf>=0)then
1627    nwarning=nwarning+1
1628    bdtot_index=1
1629    do isppol=1,dtset%nsppol
1630      do ikpt=1,dtset%nkpt
1631        min_occ=two
1632        nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
1633        do iband=1,nband_k
1634          if(occ(bdtot_index)<min_occ)min_occ=occ(bdtot_index)
1635          bdtot_index=bdtot_index+1
1636        end do
1637        if(min_occ>0.01_dp)then
1638          if(dtset%nsppol==1)then
1639            write(message, '(a,i0,3a,f7.3,5a)' )&
1640 &           'For k-point number ',ikpt,',',ch10,&
1641 &           'The minimal occupation factor is',min_occ,'.',ch10,&
1642 &           'An adequate monitoring of convergence requires it to be  at most 0.01_dp.',ch10,&
1643 &           'Action: increase slightly the number of bands.'
1644          else
1645            write(message, '(a,i0,3a,i0,a,f7.3,5a)' )&
1646 &           'For k-point number ',ikpt,', and',ch10,&
1647 &           'for spin polarization',isppol,'the minimal occupation factor is',min_occ,'.',ch10,&
1648 &           'An adequate monitoring of convergence requires it to be at most 0.01_dp.',ch10,&
1649 &           'Action: increase slightly the number of bands.'
1650          end if
1651          MSG_WARNING(message)
1652          exit ! It is enough if one lack of adequate occupation is identified, so exit.
1653        end if
1654      end do
1655    end do
1656  end if
1657 
1658  if (iscf>0.or.iscf==-3 .or. (dtset%usewvl==1 .and. iscf==0)) then
1659 
1660 !  PAW: Build new rhoij quantities from new occ then symetrize them
1661 !  Compute and add the compensation density to rhowfr to get the total density
1662    if (psps%usepaw==1) then
1663      call timab(555,1,tsec)
1664      if (paral_atom) then
1665        ABI_DATATYPE_ALLOCATE(pawrhoij_unsym,(natom))
1666        nspden_rhoij=pawrhoij_get_nspden(dtset%nspden,dtset%nspinor,dtset%pawspnorb)
1667        call pawrhoij_alloc(pawrhoij_unsym,dtset%pawcpxocc,nspden_rhoij,dtset%nspinor,&
1668 &       dtset%nsppol,dtset%typat,pawtab=pawtab,use_rhoijp=0)
1669      else
1670        pawrhoij_unsym => pawrhoij
1671      end if
1672      if (usecprj_local==1) then
1673        call pawmkrhoij(atindx,atindx1,cprj,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,dtset%mband,mband_cprj,&
1674 &       mcprj_local,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,dtset%nspinor,dtset%nsppol,&
1675 &       occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,dtfil%unpaw,&
1676 &       dtset%usewvl,dtset%wtk)
1677      else
1678        mcprj_tmp=my_nspinor*mband_cprj*dtset%mkmem*dtset%nsppol
1679        ABI_DATATYPE_ALLOCATE(cprj_tmp,(natom,mcprj_tmp))
1680        call pawcprj_alloc(cprj_tmp,0,gs_hamk%dimcprj)
1681        call ctocprj(atindx,cg,1,cprj_tmp,gmet,gprimd,0,0,0,dtset%istwfk,kg,dtset%kptns,&
1682 &       mcg,mcprj_tmp,dtset%mgfft,dtset%mkmem,mpi_enreg,psps%mpsang,dtset%mpw,&
1683 &       dtset%natom,nattyp,dtset%nband,dtset%natom,dtset%ngfft,dtset%nkpt,dtset%nloalg,&
1684 &       npwarr,dtset%nspinor,dtset%nsppol,ntypat,dtset%paral_kgb,ph1d,psps,rmet,dtset%typat,&
1685 &       ucvol,dtfil%unpaw,xred,ylm,ylmgr_dum)
1686        call pawmkrhoij(atindx,atindx1,cprj_tmp,gs_hamk%dimcprj,dtset%istwfk,dtset%kptopt,&
1687 &       dtset%mband,mband_cprj,mcprj_tmp,dtset%mkmem,mpi_enreg,natom,dtset%nband,dtset%nkpt,&
1688 &       dtset%nspinor,dtset%nsppol,occ,dtset%paral_kgb,paw_dmft,dtset%pawprtvol,pawrhoij_unsym,&
1689 &       dtfil%unpaw,dtset%usewvl,dtset%wtk)
1690        call pawcprj_free(cprj_tmp)
1691        ABI_DATATYPE_DEALLOCATE(cprj_tmp)
1692      end if
1693      call timab(555,2,tsec)
1694 !    Build symetrized packed rhoij and compensated pseudo density
1695      cplex=1;ipert=0;idir=0;qpt(:)=zero
1696      if(dtset%usewvl==0) then
1697        call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1698 &       my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1699 &       dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
1700 &       symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat,rhog=rhog)
1701      else
1702 !      here do not pass rhog, we do not use it
1703        call pawmkrho(1,compch_fft,cplex,gprimd,idir,indsym,ipert,mpi_enreg,&
1704 &       my_natom,natom,dtset%nspden,dtset%nsym,ntypat,dtset%paral_kgb,pawang,pawfgr,pawfgrtab,&
1705 &       dtset%pawprtvol,pawrhoij,pawrhoij_unsym,pawtab,qpt,rhowfg,rhowfr,rhor,rprimd,dtset%symafm,&
1706 &       symrec,dtset%typat,ucvol,dtset%usewvl,xred,pawnhat=nhat)
1707 !      In WVL: copy density to BigDFT object:
1708        call wvl_rho_abi2big(1,rhor,wvl%den)
1709      end if
1710      if (paral_atom) then
1711        call pawrhoij_free(pawrhoij_unsym)
1712        ABI_DATATYPE_DEALLOCATE(pawrhoij_unsym)
1713      end if
1714    end if
1715 
1716    if(paw_dmft%use_dmft==1) then
1717 !    == check noccmmp
1718      call setnoccmmp(1,0,rdum4,0,0,idum3,my_natom,natom,0,1,dtset%nsppol,0,ntypat,&
1719 &     paw_ij,pawang,dtset%pawprtvol,pawrhoij,pawtab,rdum2,idum1,dtset%typat,0,dtset%usepawu,&
1720 &     comm_atom=mpi_enreg%comm_atom,mpi_atmtab=mpi_enreg%my_atmtab)
1721    end if
1722 
1723 !  Find and print minimum and maximum total electron density and locations
1724 !  Compute density residual (if required) and its squared norm
1725    if (iscf>=0) then
1726      if (psps%usepaw==0) then
1727        call prtrhomxmn(std_out,mpi_enreg,dtset%nfft,dtset%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
1728      else
1729        call prtrhomxmn(std_out,mpi_enreg,nfftf,pawfgr%ngfft,dtset%nspden,1,rhor,ucvol=ucvol)
1730      end if
1731      if (optres==1) then
1732        nvresid=rhor-nvresid
1733        call sqnorm_v(1,nfftf,nres2,dtset%nspden,optres,nvresid,mpi_comm_sphgrid=mpi_comm_sphgrid)
1734        tauresid=taur-tauresid
1735      end if
1736    end if
1737 
1738  end if ! iscf>0 or iscf=-3
1739 
1740  if(psps%usepaw==1.and.(iscf>=0.or.iscf==-3))  then
1741    ABI_DEALLOCATE(rhowfr)
1742    ABI_DEALLOCATE(rhowfg)
1743  end if
1744 
1745  call timab(994,2,tsec)
1746 
1747  if (iscf==-1) then
1748 !  Eventually compute the excited states within tddft
1749    call timab(995,1,tsec)
1750    if (psps%usepaw==1) then
1751 !    In case of PAW calculation, have to transfer kxc from the fine to the coarse grid:
1752      ABI_ALLOCATE(cgrkxc,(dtset%nfft,nkxc))
1753      do ikxc=1,nkxc
1754        call transgrid(1,mpi_enreg,1,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrkxc(:,ikxc),kxc(:,ikxc))
1755      end do
1756      call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
1757 &     kg,cgrkxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
1758 &     ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
1759      ABI_DEALLOCATE(cgrkxc)
1760    else
1761      call tddft(cg,dtfil,dtset,eigen,etotal,gmet,gprimd,gsqcut,&
1762 &     kg,kxc,dtset%mband,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nfft,&
1763 &     ngfftdiel,dtset%nkpt,nkxc,npwarr,dtset%nspinor,dtset%nsppol,occ,ucvol,wffnew)
1764    end if
1765    call timab(995,2,tsec)
1766 
1767  else
1768 !  Eventually compute the susceptibility matrix and the
1769 !  dielectric matrix when istep_mix is equal to 1 or dielstrt
1770    call timab(996,1,tsec)
1771    call testsusmat(computesusmat,dielop,dielstrt,dtset,istep_mix) !test if the matrix is to be computed
1772    if(computesusmat) then
1773      dielar(1)=dtset%diecut;dielar(2)=dtset%dielng
1774      dielar(3)=dtset%diemac;dielar(4)=dtset%diemix
1775      dielar(5)=dtset%diegap;dielar(6)=dtset%dielam
1776      dielar(7)=dtset%diemix;if (iscf>=10) dielar(7)=dtset%diemixmag
1777      usetimerev=1
1778      if (psps%usepaw==1.and.dtset%pawspnorb>0.and.dtset%kptopt/=1.and.dtset%kptopt/=2) usetimerev=0
1779      neglect_pawhat=1-dtset%pawsushat
1780      call suscep_stat(atindx,atindx1,cg,cprj,dielar,&
1781 &     gs_hamk%dimcprj,doccde,eigen,gbound_diel,gprimd,&
1782 &     irrzondiel,dtset%istwfk,kg,kg_diel,lmax_diel,&
1783 &     dtset%mband,mcg,mcprj_local,mgfftdiel,dtset%mkmem,mpi_enreg,dtset%mpw,natom,dtset%nband,&
1784 &     neglect_pawhat,nfftdiel,ngfftdiel,&
1785 &     dtset%nkpt,npwarr,npwdiel,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%nsym,ntypat,&
1786 &     occ,dtset%occopt,pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,&
1787 &     susmat,dtset%symafm,dtset%symrel,dtset%tnons,dtset%typat,ucvol,&
1788 &     dtfil%unpaw,usecprj_local,psps%usepaw,usetimerev,dtset%wtk,ylmdiel)
1789    end if
1790    call timab(996,2,tsec)
1791 
1792  end if ! end condition on iscf
1793 
1794  call destroy_hamiltonian(gs_hamk)
1795 
1796  if (psps%usepaw==1) then
1797    if (usecprj==0) then
1798      call pawcprj_free(cprj_local)
1799      ABI_DATATYPE_DEALLOCATE(cprj_local)
1800    end if
1801  end if
1802 
1803  if(dtset%usewvl==0) then
1804    ABI_DEALLOCATE(EigMin)
1805    ABI_DEALLOCATE(doccde)
1806 #if defined HAVE_GPU_CUDA
1807    if(dtset%use_gpu_cuda==1) then
1808      call gpu_finalize_ham_data()
1809    end if
1810 #endif
1811  end if
1812 
1813  call timab(980,2,tsec)
1814 
1815  DBG_EXIT("COLL")
1816 
1817  contains

ABINIT/wvl_comm_eigen [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_comm_eigen

FUNCTION

  Computes occupations for the wavelet case
  Using BigDFT routines

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2241 subroutine wvl_comm_eigen()
2242 
2243 
2244 !This section has been created automatically by the script Abilint (TD).
2245 !Do not modify the following lines by hand.
2246 #undef ABI_FUNC
2247 #define ABI_FUNC 'wvl_comm_eigen'
2248 !End of the abilint section
2249 
2250  implicit none
2251 
2252 !Arguments ------------------------------------
2253 
2254 !Local variables-------------------------------
2255 #if defined HAVE_BIGDFT
2256  integer:: ikpt,norb,shift
2257 #endif
2258 
2259 ! *************************************************************************
2260 
2261    DBG_ENTER("COLL")
2262 
2263 #if defined HAVE_BIGDFT
2264    if(wvlbigdft) then
2265 !  Communicates eigenvalues to all procs.
2266 !  This will print out the eigenvalues and Fermi level.
2267      call eigensystem_info(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl,0.d0,&
2268 &     wvl%wfs%ks%Lzd%Glr%wfd%nvctr_c+7*wvl%wfs%ks%Lzd%Glr%wfd%nvctr_f,&
2269 &     wvl%wfs%ks%orbs,wvl%wfs%ks%psi)
2270    else
2271 !  Send all eigenvalues to all procs.
2272 !  I simply communicate eigenvalues: I do not print them into screen, nor calculate Fermi-level.
2273      norb=wvl%wfs%ks%orbs%norb
2274      if (mpi_enreg%nproc_wvl > 1) then
2275        shift=1
2276        do ikpt = 1, wvl%wfs%ks%orbs%nkpts
2277          call xmpi_bcast(wvl%wfs%ks%orbs%eval(shift:shift+norb-1),wvl%wfs%ks%orbs%ikptproc(ikpt),mpi_enreg%comm_wvl,ierr)
2278          shift=shift+norb
2279        end do
2280      end if
2281    end if
2282 
2283 !Copy eigenvalues from BigDFT object to "eigen"
2284    call wvl_eigen_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,eigen,2,wvl%wfs)
2285 
2286 #else
2287    BIGDFT_NOTENABLED_ERROR()
2288 #endif
2289 
2290    DBG_EXIT("COLL")
2291 
2292  end subroutine wvl_comm_eigen
2293 
2294 end subroutine vtorho

ABINIT/wvl_nscf_loop [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_nscf_loop

FUNCTION

  Non-self-consistent field cycle in Wavelets
  See also "wvl_nscf_loop_bigdft"

INPUTS

  nnsclo= number of non-self consistent field iterations

OUTPUT

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

1841 subroutine wvl_nscf_loop()
1842 
1843 
1844 !This section has been created automatically by the script Abilint (TD).
1845 !Do not modify the following lines by hand.
1846 #undef ABI_FUNC
1847 #define ABI_FUNC 'wvl_nscf_loop'
1848 !End of the abilint section
1849 
1850  implicit none
1851 
1852 !Arguments ------------------------------------
1853 ! integer, intent(in)                    :: istep,mcprj,nfft,nnsclo
1854 ! real(dp), intent(inout)                :: residm
1855 ! type(dataset_type), intent(in)         :: dtset
1856 ! type(MPI_type), intent(in)             :: mpi_enreg
1857 ! type(energies_type), intent(inout)     :: energies
1858 ! type(wvl_data), intent(inout)          :: wvl
1859 ! !arrays
1860 ! real(dp), intent(inout)                :: xcart(3, dtset%natom)
1861 ! real(dp), dimension(6), intent(out)    :: strsxc
1862 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
1863 
1864 !Local variables-------------------------------
1865  integer :: inonsc,ii
1866  integer,parameter :: iscf_=-1       !do not do a SCF cycle
1867  logical,parameter :: do_scf=.false. !do not do a SCF cycle
1868  logical,parameter :: wvlbigdft=.false.
1869  real(dp) :: dum,eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
1870 
1871 ! *************************************************************************
1872 
1873    DBG_ENTER("COLL")
1874 
1875    if(nnsclo_now>0) then
1876      do inonsc=1,nnsclo_now
1877        call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1878 &       istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1879 &       nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1880 &       dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1881        call wvl_hpsitopsi(cprj,dtset,energies,inonsc,mcprj_local,mpi_enreg, &
1882 &       residm,wvl,xcart)
1883        if(residm<dtset%tolwfr) exit !Exit loop if converged
1884      end do
1885 
1886    else
1887      do ii=1, dtset%nline
1888 !      Direct minimization technique: no diagonalization
1889        call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1890 &       istep,ii,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1891 &       nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1892 &       dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1893        call wvl_hpsitopsi(cprj,dtset,energies,ii,mcprj_local,mpi_enreg, &
1894 &       residm,wvl,xcart)
1895        if(residm<dtset%tolwfr) exit !Exit loop if converged
1896      end do
1897    end if
1898 
1899 !  Update energies depending on new WF
1900    energies%e_kinetic=ekin
1901    energies%e_nonlocalpsp=enl
1902    energies%e_exactX=eexctx
1903    energies%e_sicdc=esicdc
1904 
1905 !  Eventually update energies depending on density
1906    if (dtset%iscf<10) then
1907      energies%e_localpsp=eloc
1908      energies%e_hartree=eh
1909      energies%e_xc=exc ; energies%e_xcdc=evxc
1910    else if (nnsclo_now==0) then
1911      energies%e_localpsp=eloc
1912    end if
1913 
1914 !  End of nscf iterations
1915    if (do_last_ortho) then
1916 !    !Don't update energies (nscf cycle has been done); just recompute potential
1917      inonsc=nnsclo_now;if (nnsclo_now==0) inonsc=dtset%nline
1918      call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc,&
1919 &     istep,inonsc,iscf_,mpi_enreg%me_wvl,dtset%natom,&
1920 &     nfftf,mpi_enreg%nproc_wvl,dtset%nspden,&
1921 &     dum,do_scf,evxc,wvl,wvlbigdft,xcart,strsxc)
1922    end if
1923 
1924    DBG_EXIT("COLL")
1925 
1926  end subroutine wvl_nscf_loop

ABINIT/wvl_nscf_loop_bigdft [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_nscf_loop_bigdft

FUNCTION

  Non-self-consistent field cycle in Wavelets
  It follows the BigDFT scheme.
  See also "wvl_nscf_loop"

INPUTS

  nnsclo= number of non-self consistent field iterations

OUTPUT

  argout(sizeout)=description

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

1952 subroutine wvl_nscf_loop_bigdft()
1953 
1954 
1955 !This section has been created automatically by the script Abilint (TD).
1956 !Do not modify the following lines by hand.
1957 #undef ABI_FUNC
1958 #define ABI_FUNC 'wvl_nscf_loop_bigdft'
1959 !End of the abilint section
1960 
1961  implicit none
1962 
1963 !Arguments ------------------------------------
1964 ! integer, intent(in)                    :: istep,mcprj,nfft,nnsclo
1965 ! real(dp), intent(inout)                :: residm
1966 ! real(dp), intent(out)                  :: nres2
1967 ! type(dataset_type), intent(in)         :: dtset
1968 ! type(MPI_type), intent(in)             :: mpi_enreg
1969 ! type(energies_type), intent(inout)     :: energies
1970 ! type(wvl_data), intent(inout)          :: wvl
1971  !arrays
1972 ! real(dp), intent(inout)                :: xcart(3, dtset%natom)
1973 ! real(dp), dimension(6), intent(out)    :: strsxc
1974 ! type(pawcprj_type),dimension(dtset%natom,mcprj),intent(out)::cprj
1975 
1976 !Local variables-------------------------------
1977  integer :: inonsc
1978  integer,parameter :: iscf_=-1       !do not do a SCF cycle
1979  logical,parameter :: do_scf=.false. !do not do a SCF cycle
1980  logical,parameter :: wvlbigdft=.true.
1981  real(dp) :: eexctx,eh,ekin,eloc,enl,esicdc,evxc,exc
1982 
1983 ! *************************************************************************
1984 
1985    DBG_ENTER("COLL")
1986 
1987    call wvl_hpsitopsi(cprj,dtset, energies, istep, mcprj_local,mpi_enreg, &
1988 &   residm, wvl,xcart)
1989 
1990    if (nnsclo_now>2) then
1991      do inonsc = 2, nnsclo_now-1
1992        call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
1993 &       energies%e_hartree, energies%e_kinetic, energies%e_localpsp, &
1994 &       energies%e_nonlocalpsp, energies%e_sicdc, istep, inonsc, iscf_, &
1995 &       mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
1996 &       dtset%nspden, nres2, do_scf,energies%e_xcdc, &
1997 &       wvl, wvlbigdft, xcart, strsxc)
1998        call wvl_hpsitopsi(cprj,dtset, energies, inonsc, mcprj_local,mpi_enreg, &
1999 &       residm, wvl,xcart)
2000      end do
2001    end if
2002 
2003 !  End of nscf iterations
2004    if (do_last_ortho.and.nnsclo_now<=1) then
2005 !    !Don't update energies (nscf cycle has been done); just recompute potential
2006      call wvl_psitohpsi(dtset%diemix,eexctx,exc,eh,ekin,eloc,enl,esicdc, &
2007 &     istep, 1, iscf_, mpi_enreg%me_wvl, dtset%natom, nfftf, &
2008 &     mpi_enreg%nproc_wvl,dtset%nspden, nres2, do_scf,evxc, &
2009 &     wvl, wvlbigdft, xcart, strsxc)
2010    else if (do_last_ortho.and.nnsclo_now>1) then
2011 !    !Update energies and potential (nscf cycles are not finished)
2012      call wvl_psitohpsi(dtset%diemix, energies%e_exactX, energies%e_xc, &
2013 &     energies%e_hartree,energies%e_kinetic, energies%e_localpsp, &
2014 &     energies%e_nonlocalpsp, energies%e_sicdc, istep, nnsclo_now, iscf_, &
2015 &     mpi_enreg%me_wvl, dtset%natom, nfftf, mpi_enreg%nproc_wvl,&
2016 &     dtset%nspden, nres2, do_scf,energies%e_xcdc, &
2017 &     wvl, wvlbigdft, xcart, strsxc)
2018    end if
2019 
2020    DBG_EXIT("COLL")
2021 
2022  end subroutine wvl_nscf_loop_bigdft

ABINIT/wvl_occ [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_occ

FUNCTION

  Computes occupations for the wavelet case

INPUTS

OUTPUT

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2119 subroutine wvl_occ()
2120 
2121 
2122 !This section has been created automatically by the script Abilint (TD).
2123 !Do not modify the following lines by hand.
2124 #undef ABI_FUNC
2125 #define ABI_FUNC 'wvl_occ'
2126 !End of the abilint section
2127 
2128  implicit none
2129 
2130 !Local variables-------------------------------
2131  real(dp):: doccde_(dtset%mband*dtset%nkpt*dtset%nsppol)
2132 ! *************************************************************************
2133 
2134    DBG_ENTER("COLL")
2135 
2136 !  Compute the new occupation numbers from eigen
2137    call newocc(doccde_,eigen,energies%entropy,energies%e_fermie,dtset%spinmagntarget,&
2138 &   dtset%mband,dtset%nband,dtset%nelect,dtset%nkpt,dtset%nspinor,&
2139 &   dtset%nsppol,occ,dtset%occopt,prtvol,&
2140 &   dtset%stmbias,dtset%tphysel,dtset%tsmear,dtset%wtk)
2141 
2142 ! Copy occupations and efermi to BigDFT variables
2143    call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,1,wvl%wfs)
2144 
2145 #if defined HAVE_BIGDFT
2146 !  Copy Fermi level to BigDFT variable:
2147    wvl%wfs%ks%orbs%efermi=energies%e_fermie
2148 #endif
2149 
2150    DBG_EXIT("COLL")
2151 
2152  end subroutine wvl_occ

ABINIT/wvl_occ_bigdft [ Functions ]

[ Top ] [ Functions ]

NAME

  wvl_occ_bigdft

FUNCTION

  Computes occupations for the wavelet case
  Using BigDFT routines

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 for the wvlbigdft case, see the routine 'wvl_occ_bigdft'

PARENTS

      vtorho

CHILDREN

      eigensystem_info,wvl_eigen_abi2big,xmpi_bcast

SOURCE

2180 subroutine wvl_occ_bigdft()
2181 
2182 
2183 !This section has been created automatically by the script Abilint (TD).
2184 !Do not modify the following lines by hand.
2185 #undef ABI_FUNC
2186 #define ABI_FUNC 'wvl_occ_bigdft'
2187 !End of the abilint section
2188 
2189  implicit none
2190 
2191 ! *************************************************************************
2192 
2193    DBG_ENTER("COLL")
2194 
2195 ! Transfer occopt from ABINIT to BigDFT
2196 #if defined HAVE_BIGDFT
2197    occopt_bigdft=dtset%occopt
2198    call wvl_occopt_abi2big(occopt_bigdft,occopt_bigdft,1)
2199 
2200 !Calculate occupations using BigDFT routine
2201    call evaltoocc(mpi_enreg%me_wvl, mpi_enreg%nproc_wvl, .false., &
2202 &   dtset%tsmear, wvl%wfs%ks%orbs,  occopt_bigdft)
2203 
2204 !Pass e_fermi from BigDFT object to ABINIT variable:
2205    energies%e_fermie = wvl%wfs%ks%orbs%efermi
2206 
2207 !Copy occupations from BigDFT to ABINIT variables
2208    call wvl_occ_abi2big(dtset%mband,dtset%nkpt,dtset%nsppol,occ,2,wvl%wfs)
2209 #endif
2210 
2211    DBG_EXIT("COLL")
2212 
2213  end subroutine wvl_occ_bigdft