TABLE OF CONTENTS


ABINIT/m_ksdiago [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ksdiago

FUNCTION

  Direct diagonalization of the KS Hamiltonian H_k(G,G')

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group (MG)
  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_ksdiago
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35  use m_hamiltonian
36 
37  use m_fstrings,          only : toupper
38  use m_geometry,          only : metric
39  use m_hide_lapack,       only : xheev, xhegv, xheevx, xhegvx
40  use m_kg,                only : mkkin, mkkpg
41  use m_fftcore,           only : kpgsph
42  use m_fft,               only : fftpac
43  use m_cgtools,           only : set_istwfk
44  use m_electronpositron,  only : electronpositron_type
45  use m_mpinfo,            only : destroy_mpi_enreg, initmpi_seq
46  use m_pawtab,            only : pawtab_type
47  use m_paw_ij,            only : paw_ij_type
48  use m_pawcprj,           only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_reorder
49  use m_pawfgr,            only : pawfgr_type
50  use m_initylmg,          only : initylmg
51  use m_mkffnl,            only : mkffnl
52  use m_getghc,            only : getghc
53  use m_fourier_interpol,  only : transgrid
54  use m_cgprj,             only : getcprj
55 
56  implicit none
57 
58  private

m_ksdiago/ddiago_ctl_type [ Types ]

[ Top ] [ m_ksdiago ] [ Types ]

NAME

  ddiago_ctl_type

FUNCTION

  Structure storing the variables controlling the direct diagonalization of the Kohn-Sham Hamiltonian.
  Mainly used for debugging (and in the KSS code!)

SOURCE

 71  type, public :: ddiago_ctl_type
 72 
 73   integer :: isppol
 74    ! The spin component of the Hamiltonian (1 if nspinor==1 or nsppol==1).
 75 
 76   integer :: istwf_k
 77    ! Option defining whether time-reversal symmetry is used at particular k-points
 78    ! If 0, the code will automatically use TR symmetry if possible (depending on the k-point)
 79 
 80   integer :: nband_k
 81    ! Number of bands to be calculated.
 82 
 83   integer :: npw_k
 84   ! The number of planes waves for the wavefunctions taking into account time-reversal symmetry.
 85 
 86   integer :: npwtot
 87   ! The number of planes waves in the Hamiltonian without taking into account istwf_k
 88 
 89   integer :: nspinor
 90   ! Number of spinorial components.
 91 
 92   integer :: prtvol
 93    ! Flag controlling the verbosity.
 94 
 95   integer :: use_scalapack
 96   ! 0 if diagonalization is done in sequential on each node.
 97   ! 1 to use scalapack
 98 
 99   real(dp) :: abstol
100    ! used fro RANGE="V","I", and "A" when do_full_diago=.FALSE.
101    ! The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted
102    ! as converged when it is determined to lie in an interval [a,b] of width less than or equal to
103    !
104    !         ABSTOL + EPS *   max( |a|,|b| ) ,
105    !
106    ! where EPS is the machine precision.  If ABSTOL is less than or equal to zero, then  EPS*|T|  will be used in its place,
107    ! where |T| is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
108    !
109    ! Eigenvalues will be computed most accurately when ABSTOL is
110    ! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
111    ! If this routine returns with INFO>0, indicating that some
112    ! eigenvectors did not converge, try setting ABSTOL to 2*DLAMCH('S').
113 
114   real(dp) :: ecut
115    ! The cutoff energy for the plane wave basis set.
116 
117   real(dp) :: ecutsm
118    ! Smearing energy for plane wave kinetic energy (Ha)
119 
120   real(dp) :: effmass_free
121    ! Effective mass for electrons (usually one).
122 
123   logical :: do_full_diago
124   ! Specifies whether direct or partial diagonalization will be performed.
125   ! Meaningful only if RANGE='A'.
126 
127   integer :: ilu(2)
128    ! If RANGE='I', the indices (in ascending order) of the smallest and largest eigenvalues to be returned.
129    ! il=ilu(1), iu=ilu(2) where
130    ! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. NOT used if RANGE = 'A' or 'V'.
131 
132   integer :: nloalg(3)
133 
134   real(dp) :: kpoint(3)
135    ! The k-point in reduced coordinates at which the Hamiltonian is diagonalized.
136 
137   real(dp) :: vlu(2)
138    ! If RANGE='V', the lower and upper bounds of the interval to
139    ! be searched for eigenvalues. vl=vlu(1) and vu=vlu(2) with VL < VU.
140    ! Not referenced if RANGE = 'A' or 'I'.
141 
142   character(len=1) :: jobz
143    ! character defining whether wavefunctions are required (lapack option).
144    ! "N":  Compute eigenvalues only;
145    ! "V":  Compute eigenvalues and eigenvectors.
146 
147   character(len=1) :: range
148    ! character defining the subset of eigenstates that will be calculated (lapack option).
149    ! "A": all eigenvalues will be found.
150    ! "V": all eigenvalues in the half-open interval (VL,VU] will be found.
151    ! "I": the IL-th through IU-th eigenvalues will be found.
152 
153   !$character(len=fnlen) :: fname
154   ! The name of the file storing the eigenvectors and eigenvalues (only if jobz="V")
155 
156  end type ddiago_ctl_type
157 
158  public ::  ksdiago
159  public ::  init_ddiago_ctl

m_ksdiago/init_ddiago_ctl [ Functions ]

[ Top ] [ m_ksdiago ] [ Functions ]

NAME

  init_ddiago_ctl

FUNCTION

INPUTS

OUTPUT

PARENTS

      m_shirley,outkss

CHILDREN

      destroy_mpi_enreg,initmpi_seq,kpgsph,wrtout

SOURCE

711 subroutine init_ddiago_ctl(Dctl,jobz,isppol,nspinor,ecut,kpoint,nloalg,gmet,&
712 & nband_k,istwf_k,ecutsm,effmass_free,abstol,range,ilu,vlu,use_scalapack,prtvol)
713 
714 
715 !This section has been created automatically by the script Abilint (TD).
716 !Do not modify the following lines by hand.
717 #undef ABI_FUNC
718 #define ABI_FUNC 'init_ddiago_ctl'
719 !End of the abilint section
720 
721  implicit none
722 
723 !Arguments ------------------------------------
724 !scalars
725  integer,intent(in) :: isppol,nspinor
726  integer,optional,intent(in) :: istwf_k,prtvol,use_scalapack,nband_k
727  real(dp),intent(in) :: ecut
728  real(dp),optional,intent(in) :: ecutsm,effmass_free
729  real(dp),optional,intent(in) :: abstol
730  character(len=*),intent(in) :: jobz
731  character(len=*),optional,intent(in) :: range
732  type(ddiago_ctl_type),intent(out) :: Dctl
733 !arrays
734  integer,intent(in) :: nloalg(3)
735  integer,optional,intent(in) :: ilu(2)
736  real(dp),intent(in) :: kpoint(3)
737  real(dp),optional,intent(in) :: vlu(2)
738  real(dp),intent(in) :: gmet(3,3)
739 
740 !Local variables-------------------------------
741 !scalars
742  integer :: npw_k
743  logical :: ltest
744  character(len=500) :: msg
745  type(MPI_type) :: MPI_enreg_seq
746 !arrays
747  integer,allocatable :: kg_k(:,:)
748 
749 ! *************************************************************************
750 
751  call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type.
752 
753  Dctl%isppol  = isppol
754  Dctl%nspinor = nspinor
755  Dctl%kpoint  = kpoint
756 
757  if (PRESENT(istwf_k)) then
758   Dctl%istwf_k = istwf_k
759  else
760   Dctl%istwf_k = set_istwfk(kpoint)
761  end if
762 
763  ABI_CHECK(Dctl%istwf_k==1,"istwf_k/=1 not coded")
764 
765  Dctl%jobz   = toupper(jobz(1:1))
766  Dctl%range  = "A"
767  if (PRESENT(range)) then
768   Dctl%range = toupper(range)
769  end if
770 
771  Dctl%ecut = ecut
772  Dctl%ecutsm = zero; if (PRESENT(ecutsm)) Dctl%ecutsm = ecutsm
773  Dctl%effmass_free = one; if (PRESENT(effmass_free)) Dctl%effmass_free = effmass_free
774  Dctl%nloalg  = nloalg
775  Dctl%prtvol = 0; if (PRESENT(prtvol)) Dctl%prtvol = prtvol
776  Dctl%abstol = -tol8; if (PRESENT(abstol)) Dctl%abstol = abstol
777 
778  ABI_ALLOCATE(kg_k,(3,0))
779 
780 ! * Total number of G-vectors for this k-point with istwf_k=1.
781  call kpgsph(ecut,0,gmet,0,0,1,kg_k,kpoint,0,MPI_enreg_seq,0,Dctl%npwtot)
782 
783 ! * G-vectors taking into account time-reversal symmetry.
784  call kpgsph(ecut,0,gmet,0,0,istwf_k,kg_k,kpoint,0,MPI_enreg_seq,0,npw_k)
785 
786  Dctl%npw_k = npw_k
787  ABI_DEALLOCATE(kg_k)
788 
789  Dctl%do_full_diago = .FALSE.
790 
791  SELECT CASE (Dctl%range)
792   CASE ("A")
793 
794   ! Check on the number of stored bands.
795   Dctl%nband_k=-1
796   if (PRESENT(nband_k)) Dctl%nband_k=nband_k
797 
798   if (Dctl%nband_k==-1.or.Dctl%nband_k>=npw_k*nspinor) then
799     Dctl%nband_k=npw_k*nspinor
800     write(msg,'(4a)')ch10,&
801 &    'Since the number of bands to be computed was (-1) or',ch10,&
802 &    'too large, it has been set to the max. value npw_k*nspinor. '
803     if (Dctl%prtvol>0) then
804       call wrtout(std_out,msg,'COLL')
805     end if
806   else
807     Dctl%nband_k=nband_k
808   end if
809 
810   Dctl%do_full_diago = (Dctl%nband_k==npw_k*nspinor)
811 
812   if (Dctl%do_full_diago) then
813     write(msg,'(6a)')ch10,&
814 &    'Since the number of bands to be computed',ch10,&
815 &    'is equal to the number of G-vectors found for this k-point,',ch10,&
816 &    'the program will perform complete diagonalization.'
817   else
818     write(msg,'(6a)')ch10,&
819 &     'Since the number of bands to be computed',ch10,&
820 &     'is less than the number of G-vectors found,',ch10,&
821 &     'the program will perform partial diagonalization.'
822   end if
823   if (Dctl%prtvol>0) then
824     call wrtout(std_out,msg,'COLL')
825   end if
826 
827  CASE ("I")
828   if (.not.PRESENT(ilu)) then
829     MSG_ERROR(" ilu must be specified when range=I ")
830   end if
831   Dctl%ilu = ilu
832 
833   ltest = ( ( ilu(2)>=ilu(1) ) .and. ilu(1)>=1 .and. ilu(2)<=Dctl%npwtot )
834   write(msg,'(a,2i0)')" Illegal value for ilu: ",ilu
835   ABI_CHECK(ltest,msg)
836   Dctl%nband_k= ilu(2)-ilu(1)+1
837 
838  CASE ("V")
839   if (.not.PRESENT(vlu)) then
840     MSG_ERROR(" vlu must be specified when range=V ")
841   end if
842   Dctl%vlu = vlu
843 
844   Dctl%nband_k=-1 !??
845 
846   ltest = (vlu(2)>vlu(1))
847   write(msg,'(a,2f0.3)')" Illegal value for vlu: ",vlu
848   ABI_CHECK(ltest,msg)
849 
850  CASE DEFAULT
851    msg = " Unknown value for range: "//TRIM(Dctl%range)
852    MSG_ERROR(msg)
853  END SELECT
854 
855  ! Consider the case in which we asked for the entire set of eigenvectors
856  ! but the number of bands is less that npw_k. Therefore have to prepare the call to ZHEEVX.
857  ! TODO this has to be done in a cleaner way.
858  if (Dctl%range=="A".and. (.not.Dctl%do_full_diago)) then
859    Dctl%range="I"
860    Dctl%ilu(1) = 1
861    Dctl%ilu(2) = npw_k*nspinor
862    Dctl%nband_k= npw_k*nspinor
863  end if
864 
865  Dctl%use_scalapack=0
866  if (PRESENT(use_scalapack)) then
867    Dctl%use_scalapack=use_scalapack
868  end if
869  ABI_CHECK(Dctl%use_scalapack==0," scalapack mode not coded")
870 
871  call destroy_mpi_enreg(MPI_enreg_seq)
872 
873 end subroutine init_ddiago_ctl

m_ksdiago/ksdiago [ Functions ]

[ Top ] [ m_ksdiago ] [ Functions ]

NAME

 ksdiago

FUNCTION

  This routine performs the direct diagonalization of the Kohn-Sham Hamiltonian
  for a given k-point and spin. The routine drives the following operations:

    1) Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
       starting from the knowledge of the local potential on the real-space FFT mesh.

    2) Diagonalizing H in the plane-wave basis.

  It is called in outkss.F90 during the generation of the KSS file
  needed for a GW post-treatment. Since many-body calculations usually
  require a large number of eigenstates eigen-functions, a direct
  diagonalization of the Hamiltonian might reveal more stable than iterative
  techniques that might be problematic when several high energy states are required.
  The main drawback of the direct diagonalization is the bad scaling with the size
  of the basis set (npw**3) and the large memory requirements.
  At present, only norm-conserving pseudopotentials are implemented.

INPUTS

  kpoint(3)
  prtvol=Integer Flags  defining verbosity level
  ecut=cut-off energy for plane wave basis sphere (Ha)
  mgfftc=maximum size of 1D FFTs (coarse mesh).
  natom=number of atoms in cell.
  nfftf=(effective) number of FFT grid points in the dense FFT mesh (for this processor)
         (nfftf=nfft for norm-conserving potential runs)
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  nspden=number of density components
  Pawtab(Psps%ntypat*Psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
  Pawfgr<pawfgr_type>=fine grid parameters and related data
  Paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  vtrial(nfftf,nspden)=the trial potential
  xred(3,natom)=reduced dimensionless atomic coordinates
  comm=MPI communicator.
  [Electronpositron] <electronpositron_type>=quantities for the electron-positron annihilation.
  nfftc=Number of points in the coarse FFT mesh.
  ngfftc(18)=Info about 3D FFT for the coarse mesh, see ~abinit/doc/variables/vargs.htm#ngfft
  Diago_ctl<ddiago_ctl_type>=Datatype storing variables and options controlling the direct diagonalization.

OUTPUT

  ierr=Status error.
  onband_diago

SIDE EFFECTS

  eig_ene(:)=Pointer used for allocating and storing the eigenvalues (hartree)
    input: pointer to NULL
     output: eig_ene(onband_diago)=The calculatated eigenvalues in ascending order.

  eig_vec(:,:,:)=Pointer used for allocating and holding the wave functions at this k-point and spin.
    input: pointer to NULL
    output: eig_vec(2,npw_k*nspinor,onband_diago)=The calculated eigenvectors.

  Cprj_k(natom,nspinor*onband_diago) PAW only===
   input: pointer to NULL
   output: Projected eigenstates <Proj_i|Cnk> from output eigenstates.

NOTES

 * The routine can be time consuming (in particular when computing <G1|H|G2> elements for all (G1,G2)).
   So, it is recommended to call it once per run.

 * The routine RE-compute all Hamiltonian terms. So it is equivalent to an additional electronic SC cycle.
   (This has no effect is convergence was reach. If not, eigenvalues/vectors may differs from the conjugate gradient ones)

 * Please, do NOT pass Dtset% to this routine. Either use a local variable properly initialized
   or add the additional variable to ddiago_ctl_type and change the creation method accordingly.
   ksdiago is designed such that it is possible to diagonalize the Hamiltonian at an arbitrary k-point
   or spin (not efficient but easy to code). Therefore ksdiago is useful non only for
   the KSS generation but also for testing more advanced iterative algorithms as well as interpolation techniques.

PARENTS

      m_shirley,outkss

CHILDREN

      destroy_hamiltonian,destroy_mpi_enreg,fftpac,getcprj,getghc
      init_distribfft_seq,init_hamiltonian,initmpi_seq,initylmg,kpgsph
      load_k_hamiltonian,load_spin_hamiltonian,metric,mkffnl,mkkin,mkkpg
      pawcprj_alloc,pawcprj_free,pawcprj_reorder,transgrid,wrtout,xheev
      xheevx,xhegv,xhegvx,xmpi_barrier

SOURCE

253 subroutine ksdiago(Diago_ctl,nband_k,nfftc,mgfftc,ngfftc,natom,&
254 & typat,nfftf,nspinor,nspden,nsppol,Pawtab,Pawfgr,Paw_ij,&
255 & Psps,rprimd,vtrial,xred,onband_diago,eig_ene,eig_vec,Cprj_k,comm,ierr,&
256 & Electronpositron) ! Optional arguments
257 
258 
259 !This section has been created automatically by the script Abilint (TD).
260 !Do not modify the following lines by hand.
261 #undef ABI_FUNC
262 #define ABI_FUNC 'ksdiago'
263 !End of the abilint section
264 
265  implicit none
266 
267 !Arguments ------------------------------------
268 !scalars
269  integer,intent(in) :: mgfftc,natom,comm,nband_k
270  integer,intent(in) :: nfftf,nsppol,nspden,nspinor,nfftc
271  integer,intent(out) :: ierr,onband_diago
272  type(pseudopotential_type),intent(in) :: Psps
273  type(pawfgr_type),intent(in) :: Pawfgr
274  type(ddiago_ctl_type),intent(in) :: Diago_ctl
275 !arrays
276  integer,intent(in) :: typat(natom)
277  integer,intent(in) :: ngfftc(18)
278  real(dp),intent(in) :: rprimd(3,3)
279  real(dp),intent(inout) :: vtrial(nfftf,nspden)
280  real(dp),intent(in) :: xred(3,natom)
281  real(dp),pointer :: eig_ene(:),eig_vec(:,:,:)
282  type(pawcprj_type),pointer :: Cprj_k(:,:)
283  type(pawtab_type),intent(in) :: Pawtab(Psps%ntypat*Psps%usepaw)
284  type(paw_ij_type),intent(in) :: Paw_ij(natom*Psps%usepaw)
285  type(electronpositron_type),optional,pointer :: Electronpositron
286 
287 !Local variables-------------------------------
288 !scalars
289  integer,parameter :: mkmem_=1,tim_getghc=4,paral_kgb=0
290  integer :: cprj_choice,cpopt,dimffnl,ib,ider,idir,isppol,npw_k
291  integer :: ikg,master,istwf_k,exchn2n3d,prtvol
292  integer :: jj,n1,n2,n3,n4,n5,n6,negv,nkpg,nprocs,npw_k_test
293  integer :: my_rank,optder
294  integer :: ispden,ndat,type_calc,sij_opt,igsp2,cplex_ghg
295  integer :: iband,iorder_cprj,ibs1,ibs2
296  real(dp) :: cfact,ucvol,ecutsm,effmass_free,lambda,size_mat,ecut
297  logical :: do_full_diago
298  character(len=50) :: jobz,range
299  character(len=80) :: frmt1,frmt2
300  character(len=10) :: stag(2)
301  character(len=500) :: msg
302  type(MPI_type) :: MPI_enreg_seq
303  type(gs_hamiltonian_type) :: gs_hamk
304 !arrays
305  integer :: nloalg(3)
306  integer,allocatable :: kg_k(:,:)
307  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),kptns_(3,1)
308  real(dp) :: kpoint(3),ylmgr_dum(1,1,1)
309  real(dp) :: rhodum(1)
310  real(dp),allocatable :: ph3d(:,:,:),pwave(:,:)
311  real(dp),allocatable :: ffnl(:,:,:,:),kinpw(:),kpg_k(:,:)
312  real(dp),allocatable :: vlocal(:,:,:,:),ylm_k(:,:),dum_ylm_gr_k(:,:,:),vlocal_tmp(:,:,:)
313  real(dp),allocatable :: ghc(:,:),gvnlc(:,:),gsc(:,:)
314  real(dp),allocatable :: ghg_mat(:,:,:),gtg_mat(:,:,:)
315  real(dp),allocatable :: cgrvtrial(:,:)
316  real(dp),pointer :: cwavef(:,:)
317  type(pawcprj_type),allocatable :: Cwaveprj(:,:)
318 
319 ! *********************************************************************
320 
321  DBG_ENTER("COLL")
322 
323  nprocs  = xmpi_comm_size(comm)
324  my_rank = xmpi_comm_rank(comm)
325  master=0
326 
327  if (nprocs>1) then
328    MSG_WARNING("ksdiago not supported in parallel. Running in sequential.")
329  end if
330 
331  call initmpi_seq(MPI_enreg_seq) ! Fake MPI_type for sequential part.
332  call init_distribfft_seq(MPI_enreg_seq%distribfft,'c',ngfftc(2),ngfftc(3),'all')
333  if (Pawfgr%usefinegrid/=0) then
334    call init_distribfft_seq(MPI_enreg_seq%distribfft,'f',Pawfgr%ngfft(2),pawfgr%ngfft(3),'all')
335  end if
336 
337  isppol  = Diago_ctl%isppol
338  kpoint  = Diago_ctl%kpoint
339  istwf_k = Diago_ctl%istwf_k
340 !% nband_k = Diago_ctl%nband_k
341  npw_k   = Diago_ctl%npw_k
342  nloalg  = Diago_ctl%nloalg
343  ecut    = Diago_ctl%ecut
344  ecutsm  = Diago_ctl%ecutsm
345  effmass_free = Diago_ctl%effmass_free
346  prtvol  = Diago_ctl%prtvol
347 
348  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
349 
350  if (nsppol==1) stag=(/'          ','          '/)
351  if (nsppol==2) stag=(/'SPIN UP:  ','SPIN DOWN:'/)
352 !
353 !The coarse FFT mesh.
354  n1=ngfftc(1); n2=ngfftc(2); n3=ngfftc(3)
355  n4=ngfftc(4); n5=ngfftc(5); n6=ngfftc(6)
356 !
357 !====================
358 !=== Check input ====
359 !====================
360  ierr=0
361 !
362 !* istwfk must be 1 for each k-point
363  if (istwf_k/=1) then
364    write(msg,'(7a)')&
365 &   ' istwfk/=1 not allowed:',ch10,&
366 &   ' States output not programmed for time-reversal symmetry.',ch10,&
367 &   ' Action: change istwfk in input file (put it to 1 for all kpt).',ch10,&
368 &   ' Program does not stop but _KSS file will not be created...'
369    MSG_WARNING(msg)
370    ierr=ierr+1
371  end if
372 
373  if (paral_kgb/=0) then
374    write(msg,'(3a)')&
375 &   ' paral_kgb/=0 not allowed:',ch10,&
376 &   ' Program does not stop but _KSS file will not be created...'
377    MSG_WARNING(msg)
378    ierr=ierr+1
379  end if
380 
381  if (ierr/=0) RETURN ! Houston we have a problem!
382 !
383 !Initialize the Hamiltonian datatype on the coarse FFT mesh.
384  if (PRESENT(Electronpositron)) then
385    call init_hamiltonian(gs_hamk,Psps,pawtab,nspinor,nsppol,nspden,natom,typat,xred,nfftc,&
386 &   mgfftc,ngfftc,rprimd,nloalg,paw_ij=Paw_ij,usecprj=0,Electronpositron=Electronpositron)
387  else
388    call init_hamiltonian(gs_hamk,Psps,pawtab,nspinor,nsppol,nspden,natom,typat,xred,nfftc,&
389 &   mgfftc,ngfftc,rprimd,nloalg,paw_ij=Paw_ij,usecprj=0)
390  end if
391 
392 !Check on the number of stored bands.
393  if (nband_k==-1.or.nband_k>=npw_k*nspinor) then
394    onband_diago=npw_k*nspinor
395    write(msg,'(4a)')ch10,&
396 &   ' Since the number of bands to be computed was (-1) or',ch10,&
397 &   ' too large, it has been set to the max. value npw_k*nspinor. '
398    call wrtout(std_out,msg,'COLL')
399  else
400    onband_diago=nband_k
401  end if
402 
403 !do_full_diago = (onband_diago==npw_k*nspinor)
404  do_full_diago = Diago_ctl%do_full_diago
405 
406  if (do_full_diago) then
407    write(msg,'(6a)')ch10,&
408 &   ' Since the number of bands to be computed',ch10,&
409 &   ' is equal to the nb of G-vectors found for this k-pt,',ch10,&
410 &   ' the program will perform complete diagonalizations.'
411  else
412    write(msg,'(6a)')ch10,&
413 &   ' Since the number of bands to be computed',ch10,&
414 &   ' is less than the number of G-vectors found,',ch10,&
415 &   ' the program will perform partial diagonalizations.'
416  end if
417  if (prtvol>0) then
418    call wrtout(std_out,msg,'COLL')
419  end if
420 !
421 
422 !* Set up local potential vlocal with proper dimensioning, from vtrial.
423 !* Select spin component of interest if nspden<=2 as nvloc==1, for nspden==4, nvloc==4
424 !* option=2: vtrial(n1*n2*n3,ispden) --> vlocal(nd1,nd2,nd3) real case
425 
426 !nvloc=1; if (nspden==4) nvloc=4
427  ABI_MALLOC(vlocal,(n4,n5,n6,gs_hamk%nvloc))
428 
429  if (nspden/=4)then
430    if (Psps%usepaw==0.or.Pawfgr%usefinegrid==0) then
431      call fftpac(isppol,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal,2)
432    else ! Move from fine to coarse FFT mesh (PAW)
433      ABI_MALLOC(cgrvtrial,(nfftc,nspden))
434      call transgrid(1,MPI_enreg_seq,nspden,-1,0,0,paral_kgb,Pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
435      call fftpac(isppol,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal,2)
436      ABI_FREE(cgrvtrial)
437    end if
438  else
439    ABI_MALLOC(vlocal_tmp,(n4,n5,n6))
440    if (Psps%usepaw==0.or.Pawfgr%usefinegrid==0) then
441      do ispden=1,nspden
442        call fftpac(ispden,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,vtrial,vlocal_tmp,2)
443        vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
444      end do
445    else ! Move from fine to coarse FFT mesh (PAW)
446      ABI_MALLOC(cgrvtrial,(nfftc,nspden))
447      call transgrid(1,MPI_enreg_seq,nspden,-1,0,0,paral_kgb,Pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
448      do ispden=1,nspden
449        call fftpac(ispden,MPI_enreg_seq,nspden,n1,n2,n3,n4,n5,n6,ngfftc,cgrvtrial,vlocal_tmp,2)
450        vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
451      end do
452      ABI_FREE(cgrvtrial)
453    end if
454    ABI_FREE(vlocal_tmp)
455  end if
456 
457 !Continue to initialize the Hamiltonian (spin-dependent part)
458  call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
459 !
460 !* Calculate G-vectors, for this k-point.
461 !* Count the number of planewaves as a check.
462  exchn2n3d=0; ikg=0
463  ABI_MALLOC(kg_k,(3,npw_k))
464 
465  call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,kg_k,kpoint,0,MPI_enreg_seq,0,npw_k_test)
466  ABI_CHECK(npw_k_test==npw_k,"npw_k_test/=npw_k")
467 
468  call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,kg_k,kpoint,mkmem_,MPI_enreg_seq,npw_k,npw_k_test)
469 
470 !========================
471 !==== Kinetic energy ====
472 !========================
473  ABI_MALLOC(kinpw,(npw_k))
474  call mkkin(ecut,ecutsm,effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
475 !
476 !================================
477 !==== Non-local form factors ====
478 !================================
479  ABI_MALLOC(ylm_k,(npw_k,Psps%mpsang**2*Psps%useylm))
480 
481  if (Psps%useylm==1) then
482    optder=0
483    ABI_MALLOC(dum_ylm_gr_k,(npw_k,3+6*(optder/2),Psps%mpsang**2))
484    kptns_(:,1) = kpoint
485 
486 !  Here mband is not used if paral_compil_kpt=0
487    call initylmg(gprimd,kg_k,kptns_,mkmem_,MPI_enreg_seq,Psps%mpsang,npw_k,(/nband_k/),1,&
488 &   (/npw_k/),1,optder,rprimd,ylm_k,dum_ylm_gr_k)
489 
490    ABI_FREE(dum_ylm_gr_k)
491  end if
492 
493 !Compute (k+G) vectors (only if useylm=1)
494  nkpg=3*nloalg(3)
495  ABI_MALLOC(kpg_k,(npw_k,nkpg))
496  if (nkpg>0) then
497    call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
498  end if
499 
500 !Compute nonlocal form factors ffnl at all (k+G):
501  idir=0; ider=0; dimffnl=1+ider ! Now the derivative is not needed anymore.
502  ABI_MALLOC(ffnl,(npw_k,dimffnl,Psps%lmnmax,Psps%ntypat))
503 
504  call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,ffnl,Psps%ffspl,gmet,gprimd,ider,idir,Psps%indlmn,&
505 & kg_k,kpg_k,kpoint,Psps%lmnmax,Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,&
506 & Psps%ntypat,Psps%pspso,Psps%qgrid_ff,rmet,Psps%usepaw,Psps%useylm,ylm_k,ylmgr_dum)
507 
508  ABI_FREE(ylm_k)
509 
510 !Load k-dependent part in the Hamiltonian datastructure
511  ABI_MALLOC(ph3d,(2,npw_k,gs_hamk%matblk))
512  call load_k_hamiltonian(gs_hamk,kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,kinpw_k=kinpw,&
513 & kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl,ph3d_k=ph3d,&
514 & compute_ph3d=.true.,compute_gbound=.true.)
515 
516 !Prepare the call to getghc.
517  ndat=1; lambda=zero; type_calc=0         ! For applying the whole Hamiltonian
518  sij_opt=0; if (Psps%usepaw==1) sij_opt=1 ! For PAW, <k+G|1+S|k+G"> is also needed.
519 
520  cpopt=-1    ! If cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
521  if (Psps%usepaw==1.and..FALSE.) then ! TODO Calculate <p_lmn|k+G>.
522    cpopt = 0  ! <p_lmn|in> are computed here and saved
523  end if
524 
525  ABI_MALLOC(ghc  ,(2,npw_k*nspinor*ndat))
526  ABI_MALLOC(gvnlc,(2,npw_k*nspinor*ndat))
527  ABI_MALLOC(gsc  ,(2,npw_k*nspinor*ndat*(sij_opt+1)/2))
528 
529  cplex_ghg=2
530  size_mat = cplex_ghg*(npw_k*nspinor)**2*dp*b2Mb
531 !; if (Psps%usepaw==1) size_mat=two*size_mat
532 !write(msg,'(a,f0.3,a)')" Memory required by the Hamiltonian matrix: ",size_mat," [Mb]."
533 !call wrtout(std_out,msg,"COLL")
534 
535  write(msg,'(a,f0.3,a)')" Out-of-memory in ghg_mat. Memory required by the Hamiltonian matrix: ",size_mat," [Mb]."
536  ABI_STAT_MALLOC(ghg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor), ierr)
537  ABI_CHECK(ierr==0, msg)
538 
539  write(msg,'(a,f0.3,a)')" Out-of-memory in gtg_mat. Memory required by the PAW overlap operator: ",size_mat," [Mb]."
540  ABI_STAT_MALLOC(gtg_mat,(cplex_ghg,npw_k*nspinor,npw_k*nspinor*Psps%usepaw), ierr)
541  ABI_CHECK(ierr==0, msg)
542 
543  ABI_DT_MALLOC(Cwaveprj,(natom,nspinor*(1+cpopt)*gs_hamk%usepaw))
544  if (cpopt==0) then  ! Cwaveprj is ordered, see nonlop_ylm.
545    call pawcprj_alloc(Cwaveprj,0,gs_hamk%dimcprj)
546  end if
547 
548  ABI_MALLOC(pwave,(2,npw_k*nspinor))
549  pwave=zero ! Initialize plane-wave array:
550 
551  if (prtvol>0) then
552    call wrtout(std_out,' Calculating <G|H|G''> elements','PERS')
553  end if
554 
555  do igsp2=1,npw_k*nspinor ! Loop over the |beta,G''> component.
556 
557    pwave(1,igsp2)=one      ! Get <:|H|beta,G''> and <:|T_{PAW}|beta,G''>
558 
559    call getghc(cpopt,pwave,Cwaveprj,ghc,gsc,gs_hamk,gvnlc,lambda,MPI_enreg_seq,ndat,&
560 &   prtvol,sij_opt,tim_getghc,type_calc)
561 
562 !  Fill the upper triangle.
563    ghg_mat(:,1:igsp2,igsp2) = ghc(:,1:igsp2)
564    if (Psps%usepaw==1) gtg_mat(:,1:igsp2,igsp2) = gsc(:,1:igsp2)
565 
566    pwave(1,igsp2)=zero ! Reset the |G,beta> component that has been treated.
567  end do
568 
569 !Free workspace memory allocated so far.
570  ABI_FREE(pwave)
571  ABI_FREE(kinpw)
572  ABI_FREE(vlocal)
573  ABI_FREE(ghc)
574  ABI_FREE(gvnlc)
575  ABI_FREE(gsc)
576 
577  if (Psps%usepaw==1.and.cpopt==0) then
578    call pawcprj_free(Cwaveprj)
579  end if
580  ABI_DT_FREE(Cwaveprj)
581 !
582 !===========================================
583 !=== Diagonalization of <G|H|G''> matrix ===
584 !===========================================
585 !
586 !*** Allocate the pointers ***
587  ABI_MALLOC(eig_ene,(onband_diago))
588  ABI_MALLOC(eig_vec,(cplex_ghg,npw_k*nspinor,onband_diago))
589 
590  jobz =Diago_ctl%jobz  !jobz="Vectors"
591 
592  if (do_full_diago) then ! * Complete diagonalization
593 
594    write(msg,'(2a,3es16.8,3x,3a,i5)')ch10,&
595 &   ' Begin complete diagonalization for kpt= ',kpoint(:),stag(isppol),ch10,&
596 &   ' - Size of mat.=',npw_k*nspinor
597    if (prtvol>0) then
598      call wrtout(std_out,msg,'PERS')
599    end if
600 
601    if (Psps%usepaw==0) then
602      call xheev(  jobz,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,eig_ene)
603    else
604      call xhegv(1,jobz,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,gtg_mat,eig_ene)
605    end if
606 
607    eig_vec(:,:,:)=  ghg_mat
608 
609  else ! * Partial diagonalization
610 
611    range=Diago_ctl%range !range="Irange"
612 
613    write(msg,'(2a,3es16.8,3a,i5,a,i5)')ch10,&
614 &   ' Begin partial diagonalization for kpt= ',kpoint,stag(isppol),ch10,&
615 &   ' - Size of mat.=',npw_k*nspinor,' - # bnds=',onband_diago
616    if (prtvol>0) then
617      call wrtout(std_out,msg,'PERS')
618    end if
619 
620    if (Psps%usepaw==0) then
621      call xheevx(jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,zero,zero,&
622 &     1,onband_diago,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor)
623    else
624      call xhegvx(1,jobz,range,"Upper",cplex_ghg,npw_k*nspinor,ghg_mat,gtg_mat,zero,zero,&
625 &     1,onband_diago,-tol8,negv,eig_ene,eig_vec,npw_k*nspinor)
626    end if
627 
628  end if
629 
630  ABI_FREE(ghg_mat)
631  ABI_FREE(gtg_mat)
632 !
633 !========================================================
634 !==== Calculate <Proj_i|Cnk> from output eigenstates ====
635 !========================================================
636  if (Psps%usepaw==1) then
637 
638    iorder_cprj=1 !  Ordered (order does change wrt input file); will be changed later
639    ABI_DT_MALLOC(Cprj_k,(natom,nspinor*onband_diago))
640    call pawcprj_alloc(Cprj_k,0,gs_hamk%dimcprj)
641 
642    idir=0; cprj_choice=1  ! Only projected wave functions.
643 
644    do iband=1,onband_diago
645      ibs1=nspinor*(iband-1)+1
646      ibs2=ibs1; if (nspinor==2) ibs2=ibs2+1
647      cwavef => eig_vec(1:2,1:npw_k,iband)
648 
649      call getcprj(cprj_choice,0,cwavef,Cprj_k(:,ibs1:ibs2),&
650 &     gs_hamk%ffnl_k,idir,gs_hamk%indlmn,gs_hamk%istwf_k,gs_hamk%kg_k,&
651 &     gs_hamk%kpg_k,gs_hamk%kpt_k,gs_hamk%lmnmax,gs_hamk%mgfft,MPI_enreg_seq,&
652 &     gs_hamk%natom,gs_hamk%nattyp,gs_hamk%ngfft,gs_hamk%nloalg,gs_hamk%npw_k,gs_hamk%nspinor,&
653 &     gs_hamk%ntypat,gs_hamk%phkxred,gs_hamk%ph1d,gs_hamk%ph3d_k,gs_hamk%ucvol,gs_hamk%useylm)
654    end do
655 
656 !  Reorder the cprj (order is now the same as in input file)
657    call pawcprj_reorder(Cprj_k,gs_hamk%atindx1)
658 
659 !  deallocate(cwavef)
660  end if !usepaw==1
661 
662  if (prtvol>0) then ! Write eigenvalues.
663    cfact=Ha_eV ; frmt1='(8x,9(1x,f7.2))' ; frmt2='(8x,9(1x,f7.2))'
664    write(msg,'(a,3es16.8,3x,a)')' Eigenvalues in eV for kpt= ',kpoint,stag(isppol)
665    call wrtout(std_out,msg,'COLL')
666 
667    write(msg,frmt1)(eig_ene(ib)*cfact,ib=1,MIN(9,onband_diago))
668    call wrtout(std_out,msg,'COLL')
669    if (onband_diago>9) then
670      do jj=10,onband_diago,9
671        write(msg,frmt2) (eig_ene(ib)*cfact,ib=jj,MIN(jj+8,onband_diago))
672        call wrtout(std_out,msg,'COLL')
673      end do
674    end if
675  end if
676 
677  ABI_FREE(kpg_k)
678  ABI_FREE(kg_k)
679  ABI_FREE(ph3d)
680  ABI_FREE(ffnl)
681 
682  call destroy_mpi_enreg(MPI_enreg_seq)
683  call destroy_hamiltonian(gs_hamk)
684 
685  call xmpi_barrier(comm)
686 
687  DBG_EXIT("COLL")
688 
689 end subroutine ksdiago