Line data Source code
1 : !!****m* ABINIT/m_ksdiago
2 : !! NAME
3 : !! m_ksdiago
4 : !!
5 : !! FUNCTION
6 : !! Direct diagonalization of the KS Hamiltonian H_k(G,G')
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_ksdiago
23 :
24 : use, intrinsic :: iso_c_binding
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_xmpi
29 : use m_xomp
30 : use m_hamiltonian
31 : use m_distribfft
32 : use libxc_functionals
33 : use m_ebands
34 : use m_nctk
35 : use m_dtfil
36 : use m_hdr
37 : use m_wfk
38 :
39 : use defs_datatypes, only : pseudopotential_type
40 : use defs_abitypes, only : MPI_type
41 : use m_gwdefs, only : GW_TOLQ0, GW_Q0_DEFAULT !, cone_gw, czero_gw, j_gw
42 : use m_dtset, only : dataset_type
43 : use m_fstrings, only : toupper, ktoa, itoa, sjoin, ftoa, ltoa
44 : use m_io_tools, only : iomode_from_fname, get_unit
45 : use m_yaml, only : yamldoc_t, yamldoc_open
46 : use m_numeric_tools, only : blocked_loop
47 : use m_time, only : cwtime, cwtime_report, timab
48 : use m_geometry, only : metric, normv
49 : use m_hide_lapack, only : xhegv_cplex, xheev_cplex, xheevx_cplex, xhegvx_cplex
50 : use m_slk, only : slkmat_dp_t, slk_processor_t, block_dist_1d, &
51 : compute_eigen_problem, compute_generalized_eigen_problem
52 : use m_bz_mesh, only : findnq, findq, findqg0, identk
53 : use m_kg, only : mkkin, mkkpg
54 : use m_crystal, only : crystal_t
55 : use m_fftcore, only : kpgsph, get_kg
56 : use m_fft_mesh, only : calc_ceigr, get_gfft
57 : use m_fft, only : fftpac, uplan_t, fftbox_plan3_t, zerosym
58 : use m_cgtools, only : set_istwfk
59 : use m_electronpositron, only : electronpositron_type
60 : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
61 : use m_pawtab, only : pawtab_type
62 : use m_paw_ij, only : paw_ij_type
63 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_reorder, &
64 : pawcprj_set_zero, pawcprj_mpi_sum, pawcprj_copy
65 : use m_cgprj, only : getcprj
66 : use m_pawfgr, only : pawfgr_type
67 : use m_initylmg, only : initylmg
68 : use m_mkffnl, only : mkffnl
69 : use m_getghc, only : getghc, multithreaded_getghc
70 : use m_wfd, only : wfd_t
71 : use m_vcoul, only : vcgen_t
72 : use m_occ, only : get_fact_spin_tol_empty
73 : use m_pstat, only : pstat_proc
74 :
75 : implicit none
76 :
77 : private
78 : !!***
79 :
80 : !!****t* m_ksdiago/ugb_t
81 : !! NAME
82 : !! ugb_t
83 : !!
84 : !! FUNCTION
85 : !! This object stores the wavefunctions for given (k-point, spin) in a
86 : !! PBLAS matrix distributed over bands (columns)
87 : !!
88 : !! SOURCE
89 :
90 : type, public :: ugb_t
91 :
92 : integer :: istwf_k = -1
93 : ! Storage mode of cg_k.
94 :
95 : integer :: nspinor = -1
96 : ! Number of spinors.
97 :
98 : integer :: npw_k = -1
99 : ! Number of planewaves.
100 :
101 : integer :: npwsp = -1
102 : ! nnpw_k * nspinor.
103 :
104 : integer :: nband_k = - 1
105 : ! Total number of bands (global)
106 :
107 : integer :: my_bstart = -1, my_bstop = - 1, my_nband = - 1
108 : ! 1) Initial band
109 : ! 2) Last band
110 : ! 3) Number of bands treated by this proc. 0 if idle proc.
111 :
112 : logical :: has_idle_procs
113 : ! True if there are procs in comm who don't own any column.
114 :
115 : integer,pointer :: comm
116 : ! pointer to MPI communicator in mat
117 :
118 : type(slk_processor_t) :: processor
119 :
120 : type(slkmat_dp_t) :: mat
121 : ! PBLAS matrix with MPI-distributed Fourier components (double precision)
122 : ! Local buffer: (2, npwsp * my_nband)
123 : ! Global matrix: (npwsp, nband_k)
124 :
125 : integer, allocatable :: kg_k(:,:)
126 : ! (3, npw_k)
127 : ! G-vectors in reduced coordinates.
128 :
129 : real(dp), contiguous, pointer :: cg_k(:,:,:)
130 : ! (2, npwsp * my_nband)
131 : ! NB: This is a pointer to mat%buffer_cplx
132 :
133 : type(pawcprj_type),allocatable :: cprj_k(:,:)
134 : ! (natom, nspinor * my_nband))
135 : ! PAW projections ordered according to natom and NOT according to typat.
136 : ! NOTE my_nband
137 :
138 : contains
139 :
140 : procedure :: from_diago => ugb_from_diago
141 : ! Build object by direct diagonalization of the KS Hamiltonian.
142 :
143 : procedure :: from_wfk_file => ugb_from_wfk_file
144 : ! Build object from WFK file.
145 :
146 : procedure :: free => ugb_free
147 : ! Free memory.
148 :
149 : procedure :: print => ugb_print
150 : ! Print info on object.
151 :
152 : procedure :: collect_cprj => ugb_collect_cprj
153 : ! Collect a subset of PAW cprj on all processors.
154 :
155 : end type ugb_t
156 : !!***
157 :
158 : !!****t* m_ksdiago/hyb_t
159 : !! NAME
160 : !! hyb_t
161 : !!
162 : !! FUNCTION
163 : !!
164 : !! SOURCE
165 :
166 : type, public :: hyb_t
167 :
168 : integer :: nkibz = -1, nkbz = -1
169 : integer :: nqibz = -1, nqbz = -1
170 :
171 : integer :: mg0(3) = [2, 2, 2]
172 : ! Max shifts to account for umklapps.
173 :
174 : type(wfd_t) :: wfd
175 : type(vcgen_t) :: vcgen
176 : type(ebands_t) :: ebands
177 :
178 : real(dp),allocatable :: kibz(:,:), kbz(:,:)
179 : real(dp),allocatable :: qibz(:,:), qbz(:,:), wtq(:)
180 :
181 : integer,allocatable :: kbz2ibz(:,:)
182 : ! kbz2ibz(6, nkbz))
183 :
184 : integer,allocatable :: kbz2ibz_symrel(:,:)
185 : ! kbz2ibz_symrel(6, nkbz))
186 :
187 : integer,allocatable :: qbz2ibz(:,:)
188 : ! qbz2ibz(6, nqbz))
189 :
190 : contains
191 :
192 : procedure :: from_wfk_file => hyb_from_wfk_file
193 : ! Build object from WFK file
194 :
195 : !procedure :: print => hyb_print
196 : ! Print info on object.
197 :
198 : procedure :: free => hyb_free
199 : ! Free memory.
200 :
201 : end type hyb_t
202 : !!***
203 :
204 : !!****t* m_ksdiago/ddiago_ctl_type
205 : !! NAME
206 : !! ddiago_ctl_type
207 : !!
208 : !! FUNCTION
209 : !! Structure storing the variables controlling the direct diagonalization of the Kohn-Sham Hamiltonian.
210 : !! Mainly used for debugging (and in the KSS code!)
211 : !!
212 : !! SOURCE
213 :
214 : type, public :: ddiago_ctl_type
215 :
216 : integer :: spin
217 : ! The spin component of the Hamiltonian (1 if nspinor==1 or nsppol==1).
218 :
219 : integer :: istwf_k
220 : ! Option defining whether time-reversal symmetry is used at particular k-points
221 : ! If 0, the code will automatically use TR symmetry if possible (depending on the k-point)
222 :
223 : integer :: nband_k
224 : ! Number of bands to be calculated.
225 :
226 : integer :: npw_k
227 : ! The number of planes waves for the wavefunctions taking into account time-reversal symmetry.
228 :
229 : integer :: npwtot
230 : ! The number of planes waves in the Hamiltonian without taking into account istwf_k
231 :
232 : integer :: nspinor
233 : ! Number of spinorial components.
234 :
235 : integer :: prtvol
236 : ! Flag controlling the verbosity.
237 :
238 : integer :: use_scalapack
239 : ! 0 if diagonalization is done in sequential on each node.
240 : ! 1 to use scalapack
241 : ! TODO Not implemented
242 :
243 : real(dp) :: abstol
244 : ! used for RANGE= "V", "I", and "A" when do_full_diago=.FALSE.
245 : ! The absolute error tolerance for the eigenvalues. An approximate eigenvalue is accepted
246 : ! as converged when it is determined to lie in an interval [a,b] of width less than or equal to
247 : !
248 : ! ABSTOL + EPS * max( |a|,|b| ) ,
249 : !
250 : ! where EPS is the machine precision. If ABSTOL is less than or equal to zero, then EPS*|T| will be used in its place,
251 : ! where |T| is the 1-norm of the tridiagonal matrix obtained by reducing A to tridiagonal form.
252 : !
253 : ! Eigenvalues will be computed most accurately when ABSTOL is
254 : ! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
255 : ! If this routine returns with INFO>0, indicating that some
256 : ! eigenvectors did not converge, try setting ABSTOL to 2*DLAMCH('S').
257 :
258 : real(dp) :: ecut
259 : ! The cutoff energy for the plane wave basis set.
260 :
261 : real(dp) :: ecutsm
262 : ! Smearing energy for plane wave kinetic energy (Ha)
263 :
264 : real(dp) :: effmass_free
265 : ! Effective mass for electrons (usually one).
266 :
267 : logical :: do_full_diago
268 : ! Specifies whether direct or partial diagonalization will be performed.
269 : ! Meaningful only if RANGE='A'.
270 :
271 : integer :: ilu(2)
272 : ! If RANGE='I', the indices (in ascending order) of the smallest and largest eigenvalues to be returned.
273 : ! il=ilu(1), iu=ilu(2) where
274 : ! 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. NOT used if RANGE = 'A' or 'V'.
275 :
276 : integer :: nloalg(3)
277 :
278 : real(dp) :: kpoint(3)
279 : ! The k-point in reduced coordinates at which the Hamiltonian is diagonalized.
280 :
281 : real(dp) :: vlu(2)
282 : ! If RANGE='V', the lower and upper bounds of the interval to
283 : ! be searched for eigenvalues. vl=vlu(1) and vu=vlu(2) with VL < VU.
284 : ! Not referenced if RANGE = 'A' or 'I'.
285 :
286 : character(len=1) :: jobz
287 : ! character defining whether wavefunctions are required (lapack option).
288 : ! "N": Compute eigenvalues only;
289 : ! "V": Compute eigenvalues and eigenvectors.
290 :
291 : character(len=1) :: range
292 : ! character defining the subset of eigenstates that will be calculated (lapack option).
293 : ! "A": all eigenvalues will be found.
294 : ! "V": all eigenvalues in the half-open interval (VL,VU] will be found.
295 : ! "I": the IL-th through IU-th eigenvalues will be found.
296 :
297 : !$character(len=fnlen) :: fname
298 : ! The name of the file storing the eigenvectors and eigenvalues (only if jobz="V")
299 :
300 : end type ddiago_ctl_type
301 : !!***
302 :
303 :
304 : !!****t* m_ksdiago/psbands_t
305 : !! NAME
306 : !! psbands_t
307 : !!
308 : !! FUNCTION
309 : !!
310 : !! SOURCE
311 :
312 : type, public :: psbands_t
313 :
314 : integer :: nb_tot = -1
315 : ! Total number of states (protected + pseudo bands)
316 :
317 : integer :: nb_protected = -1
318 : ! Number of protected bands.
319 :
320 : integer :: nslices = -1
321 : ! Number of slices.
322 :
323 : integer :: maxsto_per_slice = -1
324 : ! Max number of pseudo bands per slice.
325 :
326 : real(dp) :: efrac
327 :
328 : integer,allocatable :: subspace(:,:)
329 : ! (3, nslices)
330 : ! For each slice, the first and last band index and the number of pseudo bands in the slice.
331 :
332 : real(dp),allocatable :: ps_eig(:)
333 : ! (nb_tot)
334 : ! eigenvalues (KS + pseudo energies)
335 :
336 : contains
337 :
338 : procedure :: init => psbands_init
339 : ! Initialize the object
340 :
341 : procedure :: band2slice => psbands_band2slice
342 : ! Return the slice index from the band index.
343 :
344 : procedure :: free => psbands_free
345 : ! Free memory.
346 :
347 : end type psbands_t
348 : !!***
349 :
350 : public :: ksdiago
351 : public :: init_ddiago_ctl
352 :
353 : !!***
354 :
355 : contains
356 : !!***
357 :
358 : !!****f* m_ksdiago/ksdiago
359 : !! NAME
360 : !! ksdiago
361 : !!
362 : !! FUNCTION
363 : !! This routine performs the direct diagonalization of the Kohn-Sham Hamiltonian
364 : !! for a given k-point and spin. The routine drives the following operations:
365 : !!
366 : !! 1) Re-computing <G|H|G_prim> matrix elements for all (G, G_prim).
367 : !! starting from the knowledge of the local potential on the real-space FFT mesh.
368 : !!
369 : !! 2) Diagonalizing H in the plane-wave basis.
370 : !!
371 : !! It is called in outkss.F90 during the generation of the KSS file
372 : !! needed for a GW post-treatment. Since many-body calculations usually
373 : !! require a large number of eigenstates eigen-functions, a direct
374 : !! diagonalization of the Hamiltonian might reveal more stable than iterative
375 : !! techniques that might be problematic when several high energy states are required.
376 : !! The main drawback of the direct diagonalization is the bad scaling with the size
377 : !! of the basis set (npw**3) and the large memory requirements.
378 : !!
379 : !! INPUTS
380 : !! kpoint(3)
381 : !! prtvol=Integer Flags defining verbosity level
382 : !! ecut=cut-off energy for plane wave basis sphere (Ha)
383 : !! mgfftc=maximum size of 1D FFTs (coarse mesh).
384 : !! natom=number of atoms in cell.
385 : !! nfftf=(effective) number of FFT grid points in the dense FFT mesh (for this processor)
386 : !! (nfftf=nfft for norm-conserving potential runs)
387 : !! nspinor=number of spinorial components of the wavefunctions
388 : !! nsppol=1 for unpolarized, 2 for spin-polarized
389 : !! nspden=number of density components
390 : !! pawtab(psps%ntypat*psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
391 : !! pawfgr<pawfgr_type>=fine grid parameters and related data
392 : !! paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
393 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
394 : !! rprimd(3,3)=dimensional primitive translations for real space (bohr)
395 : !! vtrial(nfftf,nspden)=the trial potential
396 : !! xred(3,natom)=reduced dimensionless atomic coordinates
397 : !! comm=MPI communicator.
398 : !! [electronpositron] <electronpositron_type>=quantities for the electron-positron annihilation.
399 : !! nfftc=Number of points in the coarse FFT mesh.
400 : !! ngfftc(18)=Info about 3D FFT for the coarse mesh, see ~abinit/doc/variables/vargs.htm#ngfft
401 : !! Diago_ctl<ddiago_ctl_type>=Datatype storing variables and options controlling the direct diagonalization.
402 : !!
403 : !! OUTPUT
404 : !! ierr=Status error.
405 : !! onband_diago
406 : !!
407 : !! SIDE EFFECTS
408 : !! eig_ene(:)=Pointer used for allocating and storing the eigenvalues (hartree)
409 : !! input: pointer to NULL
410 : !! output: eig_ene(onband_diago)=The calculatated eigenvalues in ascending order.
411 : !!
412 : !! eig_vec(:,:,:)=Pointer used for allocating and holding the wave functions at this k-point and spin.
413 : !! input: pointer to NULL
414 : !! output: eig_vec(2,npw_k*nspinor,onband_diago)=The calculated eigenvectors.
415 : !!
416 : !! cprj_k(natom,nspinor*onband_diago) PAW only===
417 : !! input: pointer to NULL
418 : !! output: Projected eigenstates <Proj_i|Cnk> from output eigenstates.
419 : !!
420 : !! NOTES
421 : !! * The routine can be time consuming (in particular when computing <G1|H|G2> elements for all (G1,G2)).
422 : !! So, it is recommended to call it once per run.
423 : !!
424 : !! * The routine RE-compute all Hamiltonian terms. So it is equivalent to an additional electronic SCF cycle.
425 : !! (This has no effect is convergence was reached.
426 : !! If not, eigenvalues/vectors may differs from the conjugate gradient ones)
427 : !!
428 : !! * Please, do NOT pass Dtset% to this routine. Either use a local variable properly initialized
429 : !! or add the additional variable to ddiago_ctl_type and change the creation method accordingly.
430 : !! ksdiago is designed such that it is possible to diagonalize the Hamiltonian at an arbitrary k-point
431 : !! or spin (not efficient but easy to code). Therefore ksdiago is useful non only for
432 : !! the KSS generation but also for testing more advanced iterative algorithms as well as interpolation techniques.
433 : !!
434 : !! SOURCE
435 :
436 10 : subroutine ksdiago(Diago_ctl, nband_k, nfftc, mgfftc, ngfftc, natom, &
437 10 : typat, nfftf, nspinor, nspden, nsppol, pawtab, pawfgr, paw_ij,&
438 10 : psps, rprimd, vtrial, xred, onband_diago, eig_ene, eig_vec, cprj_k, comm, ierr,&
439 : electronpositron) ! Optional arguments
440 :
441 : !Arguments ------------------------------------
442 : !scalars
443 : integer,intent(in) :: mgfftc,natom,comm,nband_k,nfftf,nsppol,nspden,nspinor,nfftc
444 : integer,intent(out) :: ierr, onband_diago
445 : type(pseudopotential_type),intent(in) :: psps
446 : type(pawfgr_type),intent(in) :: pawfgr
447 : type(ddiago_ctl_type),intent(in) :: Diago_ctl
448 : !arrays
449 : integer,intent(in) :: typat(natom), ngfftc(18)
450 : real(dp),intent(in) :: rprimd(3,3)
451 : real(dp),intent(inout) :: vtrial(nfftf,nspden)
452 : real(dp),intent(in) :: xred(3,natom)
453 : real(dp),pointer :: eig_ene(:),eig_vec(:,:,:)
454 : type(pawcprj_type),pointer :: cprj_k(:,:)
455 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
456 : type(paw_ij_type),intent(in) :: paw_ij(natom*psps%usepaw)
457 : type(electronpositron_type),optional,pointer :: Electronpositron
458 :
459 : !Local variables-------------------------------
460 : !scalars
461 : integer,parameter :: mkmem1 = 1, tim_getghc = 4, paral_kgb0 = 0, master = 0, ndat1 = 1, ncomp1 = 1
462 : integer :: cprj_choice,cpopt,dimffnl,ib,ider,idir,spin,npw_k
463 : integer :: ikg,istwf_k,exchn2n3d,prtvol
464 : integer :: jj,n1,n2,n3,n4,n5,n6,negv,nkpg,nproc,npw_k_test,my_rank,optder
465 : integer :: type_calc,sij_opt,igsp2,cplex_ghg,iband,ibs1,ibs2
466 : real(dp),parameter :: lambda0 = zero
467 : real(dp) :: ucvol,ecutsm,effmass_free,size_mat,ecut
468 : logical :: do_full_diago
469 : character(len=50) :: jobz,range
470 : character(len=80) :: frmt1,frmt2
471 : character(len=10) :: stag(2)
472 : character(len=500) :: msg
473 10 : type(MPI_type) :: mpi_enreg_seq
474 10 : type(gs_hamiltonian_type) :: gs_hamk
475 : !arrays
476 : integer :: nloalg(3)
477 10 : integer,allocatable :: kg_k(:,:)
478 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),kptns_(3,1),kpoint(3),ylmgr_dum(1,1,1)
479 10 : real(dp),allocatable :: ph3d(:,:,:),bras(:,:),ffnl(:,:,:,:),kinpw(:),kpg_k(:,:)
480 10 : real(dp),allocatable :: vlocal(:,:,:,:),ylm_k(:,:),dum_ylm_gr_k(:,:,:), vxctaulocal(:,:,:,:,:)
481 10 : real(dp),allocatable :: ghc(:,:),gvnlxc(:,:),gsc(:,:),ghg_mat(:,:,:),gsg_mat(:,:,:)
482 10 : real(dp),pointer :: cwavef(:,:)
483 10 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
484 : ! *********************************************************************
485 :
486 10 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
487 :
488 10 : if (nproc > 1) then
489 0 : ABI_WARNING("ksdiago not supported in parallel. Running in sequential.")
490 : end if
491 :
492 10 : call initmpi_seq(mpi_enreg_seq) ! Fake MPI_type for sequential part.
493 10 : call mpi_enreg_seq%distribfft%init_seq('c', ngfftc(2), ngfftc(3), 'all')
494 10 : if (pawfgr%usefinegrid /= 0) then
495 0 : call mpi_enreg_seq%distribfft%init_seq('f', pawfgr%ngfft(2), pawfgr%ngfft(3), 'all')
496 : end if
497 :
498 10 : spin = Diago_ctl%spin
499 40 : kpoint = Diago_ctl%kpoint
500 10 : istwf_k = Diago_ctl%istwf_k
501 : !% nband_k = Diago_ctl%nband_k
502 10 : npw_k = Diago_ctl%npw_k
503 40 : nloalg = Diago_ctl%nloalg
504 10 : ecut = Diago_ctl%ecut
505 10 : ecutsm = Diago_ctl%ecutsm
506 10 : effmass_free = Diago_ctl%effmass_free
507 10 : prtvol = Diago_ctl%prtvol
508 :
509 10 : call metric(gmet, gprimd, -1, rmet, rprimd, ucvol)
510 :
511 30 : if (nsppol == 1) stag = [' ',' ']
512 10 : if (nsppol == 2) stag = ['SPIN UP: ','SPIN DOWN:']
513 :
514 : ! The coarse FFT mesh.
515 10 : n1 = ngfftc(1); n2 = ngfftc(2); n3 = ngfftc(3)
516 10 : n4 = ngfftc(4); n5 = ngfftc(5); n6 = ngfftc(6)
517 :
518 : !====================
519 : !=== Check input ====
520 : !====================
521 10 : ierr = 0
522 :
523 : ! istwfk must be 1 for each k-point
524 10 : if (istwf_k/=1) then
525 : write(msg,'(7a)')&
526 0 : ' istwfk /= 1 not allowed:',ch10,&
527 0 : ' States output not programmed for time-reversal symmetry.',ch10,&
528 0 : ' Action: change istwfk in input file (put it to 1 for all kpt).',ch10,&
529 0 : ' Program does not stop but _KSS file will not be created...'
530 0 : ABI_WARNING(msg)
531 0 : ierr = ierr + 1
532 : end if
533 :
534 10 : if (ierr /= 0) RETURN ! Houston we have a problem!
535 :
536 : ! Initialize the Hamiltonian datatype on the coarse FFT mesh.
537 10 : if (present(electronpositron)) then
538 : call gs_hamk%init(psps, pawtab, nspinor, nsppol, nspden, natom, typat, xred, nfftc, &
539 0 : mgfftc, ngfftc, rprimd, nloalg, paw_ij=paw_ij, usecprj=0, electronpositron=electronpositron)
540 : else
541 : call gs_hamk%init(psps, pawtab, nspinor, nsppol, nspden, natom, typat, xred, nfftc, &
542 10 : mgfftc, ngfftc, rprimd, nloalg, paw_ij=paw_ij, usecprj=0)
543 : end if
544 :
545 : ! Check on the number of stored bands.
546 10 : onband_diago = nband_k
547 10 : if (nband_k==-1 .or. nband_k >= npw_k*nspinor) then
548 0 : onband_diago = npw_k*nspinor
549 0 : write(msg,'(4a,i0)')ch10,&
550 0 : ' Since the number of bands to be computed was -1 or',ch10,&
551 0 : ' too large, it has been set to the maximum value npw_k*nspinor: ',npw_k*nspinor
552 0 : call wrtout(std_out, msg)
553 : end if
554 :
555 : !do_full_diago = (onband_diago==npw_k*nspinor)
556 10 : do_full_diago = Diago_ctl%do_full_diago
557 :
558 10 : if (do_full_diago) then
559 0 : write(msg,'(6a)')ch10,&
560 0 : ' Since the number of bands to be computed',ch10,&
561 0 : ' is equal to the number of G-vectors found for this kpt,',ch10,&
562 0 : ' the program will perform complete diagonalization.'
563 : else
564 10 : write(msg,'(6a)')ch10,&
565 10 : ' Since the number of bands to be computed',ch10,&
566 10 : ' is less than the number of G-vectors found,',ch10,&
567 20 : ' the program will perform partial diagonalization.'
568 : end if
569 10 : if (prtvol > 0) call wrtout(std_out, msg)
570 :
571 : ! Set up local potential vlocal with proper dimensioning, from vtrial.
572 : ! Select spin component of interest if nspden<=2 as nvloc==1, for nspden==4, nvloc==4
573 : ! option=2: vtrial(n1*n2*n3,ispden) --> vlocal(nd1,nd2,nd3) real case
574 :
575 60 : ABI_MALLOC(vlocal, (n4, n5, n6, gs_hamk%nvloc))
576 : !if (with_vxctau) then
577 : ! ABI_MALLOC(vxctaulocal,(n4,n5,n6,gs_hamk%nvloc,4))
578 : !end if
579 :
580 : ! Set up local potential vlocal on the coarse FFT mesh from vtrial taking into account the spin.
581 :
582 : call gspot_transgrid_and_pack(spin, psps%usepaw, paral_kgb0, nfftc, ngfftc, nfftf, &
583 10 : nspden, gs_hamk%nvloc, ncomp1, pawfgr, mpi_enreg_seq, vtrial, vlocal)
584 10 : call gs_hamk%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
585 :
586 : ! This for meta-gga.
587 : !if (with_vxctau) then
588 : ! call gspot_transgrid_and_pack(spin, psps%usepaw, paral_kgb0, nfftc, ngfftc, nfftf, &
589 : ! nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
590 : ! call gs_hamk%load_spin(spin, vxctaulocal=vxctaulocal)
591 : !end if
592 :
593 : ! Calculate G-vectors, for this k-point. Count also the number of planewaves as a check.
594 10 : exchn2n3d = 0; ikg = 0
595 30 : ABI_MALLOC(kg_k, (3, npw_k))
596 :
597 10 : call kpgsph(ecut, exchn2n3d, gmet, ikg, 0, istwf_k, kg_k, kpoint, 0, mpi_enreg_seq, 0, npw_k_test)
598 10 : ABI_CHECK(npw_k_test == npw_k, "npw_k_test/=npw_k")
599 10 : call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,kg_k,kpoint,mkmem1,mpi_enreg_seq,npw_k,npw_k_test)
600 :
601 : !========================
602 : !==== Kinetic energy ====
603 : !========================
604 30 : ABI_MALLOC(kinpw, (npw_k))
605 10 : call mkkin(ecut, ecutsm, effmass_free, gmet, kg_k, kinpw, kpoint, npw_k, 0, 0)
606 :
607 : !================================
608 : !==== Non-local form factors ====
609 : !================================
610 40 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*psps%useylm))
611 :
612 10 : if (psps%useylm == 1) then
613 0 : optder = 0
614 0 : ABI_MALLOC(dum_ylm_gr_k, (npw_k, 3+6*(optder/2),psps%mpsang**2))
615 0 : kptns_(:,1) = kpoint
616 :
617 : ! Here mband is not used if paral_compil_kpt=0
618 : call initylmg(gprimd, kg_k, kptns_, mkmem1, mpi_enreg_seq, psps%mpsang, npw_k, [nband_k], 1, &
619 0 : [npw_k], 1, optder, rprimd, ylm_k, dum_ylm_gr_k)
620 :
621 0 : ABI_FREE(dum_ylm_gr_k)
622 : end if
623 :
624 : ! Compute (k+G) vectors (only if useylm=1)
625 10 : nkpg = 3 * nloalg(3)
626 40 : ABI_MALLOC(kpg_k, (npw_k, nkpg))
627 10 : if (nkpg > 0) call mkkpg(kg_k, kpg_k, kpoint, nkpg, npw_k)
628 :
629 : ! Compute nonlocal form factors ffnl at all (k+G):
630 10 : idir=0; ider=0; dimffnl=1+ider ! Now the derivative is not needed anymore.
631 50 : ABI_MALLOC(ffnl, (npw_k, dimffnl, psps%lmnmax, psps%ntypat))
632 :
633 : call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, gmet, gprimd, ider, idir, psps%indlmn, &
634 : kg_k, kpg_k, kpoint, psps%lmnmax, psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw_k, &
635 10 : psps%ntypat, psps%pspso, psps%qgrid_ff, rmet, psps%usepaw, psps%useylm, ylm_k, ylmgr_dum)
636 :
637 10 : ABI_FREE(ylm_k)
638 :
639 : ! Load k-dependent part in the Hamiltonian datastructure
640 40 : ABI_MALLOC(ph3d, (2, npw_k, gs_hamk%matblk))
641 : call gs_hamk%load_k(kpt_k=kpoint, istwf_k=istwf_k, npw_k=npw_k, kinpw_k=kinpw, &
642 10 : kg_k=kg_k, kpg_k=kpg_k, ffnl_k=ffnl, ph3d_k=ph3d, compute_ph3d=.true., compute_gbound=.true.)
643 :
644 : ! Prepare call to getghc.
645 10 : type_calc = 0 ! For applying the whole Hamiltonian
646 10 : sij_opt = 0; if (psps%usepaw==1) sij_opt = 1 ! For PAW, <k+G|S|k+G"> is also needed.
647 :
648 10 : cpopt = -1 ! If cpopt=-1, <p_lmn|in> (and derivatives) are computed here (and not saved)
649 : if (psps%usepaw==1.and..FALSE.) then ! TODO Calculate <p_lmn|k+G>.
650 : cpopt = 0 ! <p_lmn|in> are computed here and saved
651 : end if
652 :
653 30 : ABI_MALLOC(ghc, (2, npw_k*nspinor*ndat1))
654 20 : ABI_MALLOC(gvnlxc, (2, npw_k*nspinor*ndat1))
655 30 : ABI_MALLOC(gsc, (2, npw_k*nspinor*ndat1*(sij_opt+1)/2))
656 :
657 10 : cplex_ghg = 2
658 10 : size_mat = cplex_ghg*(npw_k*nspinor)**2*dp*b2Mb
659 10 : write(msg,'(a,f0.3,a)')" Out-of-memory in ghg_mat. Memory required by the Hamiltonian matrix: ",size_mat," [Mb]."
660 40 : ABI_STAT_MALLOC(ghg_mat, (cplex_ghg, npw_k*nspinor, npw_k*nspinor), ierr)
661 10 : ABI_CHECK(ierr == 0, msg)
662 10 : write(msg,'(a,f0.3,a)')" Out-of-memory in gsg_mat. Memory required by the PAW overlap operator: ",size_mat," [Mb]."
663 40 : ABI_STAT_MALLOC(gsg_mat, (cplex_ghg, npw_k*nspinor, npw_k*nspinor*psps%usepaw), ierr)
664 10 : ABI_CHECK(ierr == 0, msg)
665 :
666 : ! cwaveprj is ordered by atom type, see nonlop_ylm.
667 20 : ABI_MALLOC(cwaveprj, (natom, nspinor*(1+cpopt)*gs_hamk%usepaw))
668 : if (cpopt == 0) call pawcprj_alloc(cwaveprj, 0, gs_hamk%dimcprj)
669 :
670 : ! Initialize plane-wave array with zeros
671 8708 : ABI_CALLOC(bras, (2, npw_k*nspinor))
672 10 : if (prtvol > 0) call wrtout(std_out, ' Calculating <G|H|G''> elements')
673 :
674 : ! Loop over the |beta,G''> component.
675 2906 : do igsp2=1,npw_k*nspinor
676 2896 : bras(1, igsp2) = one
677 :
678 : ! Get <:|H|beta,G''> and <:|S_{PAW}|beta,G''>
679 : call getghc(cpopt, bras, cwaveprj, ghc, gsc, gs_hamk, gvnlxc, lambda0, mpi_enreg_seq, ndat1, &
680 2896 : prtvol, sij_opt, tim_getghc, type_calc)
681 :
682 : ! Fill the upper triangle.
683 1265485 : ghg_mat(:,1:igsp2,igsp2) = ghc(:,1:igsp2)
684 2896 : if (psps%usepaw == 1) gsg_mat(:,1:igsp2,igsp2) = gsc(:,1:igsp2)
685 :
686 : ! Reset the |G,beta> component that has been treated.
687 2906 : bras(1, igsp2) = zero
688 : end do
689 :
690 : ! Free workspace memory allocated so far.
691 10 : ABI_FREE(bras)
692 10 : ABI_FREE(kinpw)
693 10 : ABI_FREE(vlocal)
694 10 : ABI_FREE(ghc)
695 10 : ABI_FREE(gvnlxc)
696 10 : ABI_FREE(gsc)
697 : ABI_SFREE(vxctaulocal)
698 :
699 : if (psps%usepaw == 1 .and. cpopt == 0) call pawcprj_free(Cwaveprj)
700 10 : ABI_FREE(cwaveprj)
701 :
702 : !===========================================
703 : !=== Diagonalization of <G|H|G''> matrix ===
704 : !===========================================
705 30 : ABI_MALLOC(eig_ene, (onband_diago))
706 40 : ABI_MALLOC(eig_vec, (cplex_ghg, npw_k*nspinor, onband_diago))
707 :
708 10 : jobz = Diago_ctl%jobz !jobz="Vectors"
709 :
710 10 : if (do_full_diago) then
711 : ! Full diagonalization
712 0 : write(msg,'(6a,i0)')ch10,&
713 0 : ' Begin full diagonalization for kpt: ',trim(ktoa(kpoint)), stag(spin), ch10,&
714 0 : ' Matrix size: ', npw_k*nspinor
715 0 : call wrtout(std_out, msg)
716 :
717 0 : if (psps%usepaw == 0) then
718 0 : call xheev_cplex(jobz, "Upper", cplex_ghg, npw_k*nspinor, ghg_mat, eig_ene, msg, ierr)
719 : else
720 0 : call xhegv_cplex(1, jobz, "Upper", cplex_ghg, npw_k*nspinor, ghg_mat, gsg_mat, eig_ene, msg, ierr)
721 : end if
722 0 : ABI_CHECK(ierr == 0, msg)
723 0 : eig_vec(:,:,:)= ghg_mat
724 :
725 : else
726 : ! Partial diagonalization
727 10 : range = Diago_ctl%range !range="Irange"
728 :
729 10 : write(msg,'(2a,3es16.8,3a,i0,a,i0)')ch10,&
730 10 : ' Begin partial diagonalization for kpt= ',kpoint, stag(spin),ch10,&
731 20 : ' - Size of mat.=',npw_k*nspinor,' - # out_nband: ',onband_diago
732 10 : call wrtout(std_out, msg)
733 :
734 10 : if (psps%usepaw == 0) then
735 : call xheevx_cplex(jobz, range, "Upper", cplex_ghg, npw_k*nspinor, ghg_mat, zero, zero,&
736 10 : 1, onband_diago, -tol8, negv, eig_ene, eig_vec, npw_k*nspinor, msg, ierr)
737 : else
738 : call xhegvx_cplex(1, jobz, range, "Upper", cplex_ghg, npw_k*nspinor, ghg_mat, gsg_mat, zero, zero,&
739 0 : 1, onband_diago, -tol8, negv, eig_ene, eig_vec, npw_k*nspinor, msg, ierr)
740 : end if
741 10 : ABI_CHECK(ierr == 0, msg)
742 : end if
743 :
744 10 : ABI_FREE(ghg_mat)
745 10 : ABI_FREE(gsg_mat)
746 :
747 10 : if (prtvol > 0 .and. my_rank == master) then
748 : ! Write eigenvalues.
749 0 : frmt1 = '(8x,9(1x,f7.2))'; frmt2 = '(8x,9(1x,f7.2))'
750 0 : write(msg,'(2a,3x,a)')' Eigenvalues in eV for kpt: ', trim(ktoa(kpoint)), stag(spin)
751 0 : call wrtout(std_out, msg)
752 :
753 0 : write(msg,frmt1)(eig_ene(ib)*Ha_eV,ib=1,MIN(9,onband_diago))
754 0 : call wrtout(std_out, msg)
755 0 : if (onband_diago >9 ) then
756 0 : do jj=10,onband_diago,9
757 0 : write(msg, frmt2) (eig_ene(ib)*Ha_eV,ib=jj,MIN(jj+8,onband_diago)); call wrtout(std_out, msg)
758 : end do
759 : end if
760 : end if
761 :
762 : !========================================================
763 : !==== Calculate <Proj_i|Cnk> from output eigenstates ====
764 : !========================================================
765 10 : if (psps%usepaw == 1) then
766 :
767 0 : ABI_MALLOC(cprj_k,(natom, nspinor*onband_diago))
768 0 : call pawcprj_alloc(cprj_k, 0, gs_hamk%dimcprj)
769 :
770 0 : idir = 0; cprj_choice = 1 ! Only projected wave functions.
771 :
772 0 : do iband=1,onband_diago
773 0 : ibs1 = nspinor * (iband - 1) + 1
774 0 : ibs2 = ibs1; if (nspinor == 2) ibs2=ibs2+1
775 0 : cwavef => eig_vec(1:2,1:npw_k,iband)
776 :
777 : call getcprj(cprj_choice, 0, cwavef, cprj_k(:,ibs1:ibs2), &
778 : gs_hamk%ffnl_k, idir, gs_hamk%indlmn, gs_hamk%istwf_k, gs_hamk%kg_k, &
779 : gs_hamk%kpg_k, gs_hamk%kpt_k, gs_hamk%lmnmax, gs_hamk%mgfft, mpi_enreg_seq, 1, &
780 : gs_hamk%natom, gs_hamk%nattyp, gs_hamk%ngfft, gs_hamk%nloalg, gs_hamk%npw_k, gs_hamk%nspinor, &
781 0 : gs_hamk%ntypat, gs_hamk%phkxred, gs_hamk%ph1d, gs_hamk%ph3d_k, gs_hamk%ucvol, gs_hamk%useylm)
782 : end do
783 :
784 : ! Reorder the cprj (order is now the same as in input file)
785 0 : call pawcprj_reorder(cprj_k, gs_hamk%atindx1)
786 : end if ! usepaw
787 :
788 : ! Free memory.
789 10 : ABI_FREE(kpg_k)
790 10 : ABI_FREE(kg_k)
791 10 : ABI_FREE(ph3d)
792 10 : ABI_FREE(ffnl)
793 :
794 10 : call destroy_mpi_enreg(mpi_enreg_seq)
795 10 : call gs_hamk%free()
796 10 : call xmpi_barrier(comm)
797 :
798 20 : end subroutine ksdiago
799 : !!***
800 : !----------------------------------------------------------------------
801 :
802 : !!****f* m_ksdiago/init_ddiago_ctl
803 : !! NAME
804 : !! init_ddiago_ctl
805 : !!
806 : !! FUNCTION
807 : !!
808 : !! INPUTS
809 : !!
810 : !! OUTPUT
811 : !!
812 : !! SOURCE
813 :
814 10 : subroutine init_ddiago_ctl(Dctl, jobz, spin, nspinor, ecut, kpoint, nloalg, gmet, &
815 : nband_k, istwf_k, ecutsm, effmass_free, abstol, range, ilu, vlu, use_scalapack, prtvol)
816 :
817 : !Arguments ------------------------------------
818 : !scalars
819 : integer,intent(in) :: spin,nspinor
820 : integer,optional,intent(in) :: istwf_k,prtvol,use_scalapack,nband_k
821 : real(dp),intent(in) :: ecut
822 : real(dp),optional,intent(in) :: ecutsm,effmass_free
823 : real(dp),optional,intent(in) :: abstol
824 : character(len=*),intent(in) :: jobz
825 : character(len=*),optional,intent(in) :: range
826 : type(ddiago_ctl_type),intent(out) :: Dctl
827 : !arrays
828 : integer,intent(in) :: nloalg(3)
829 : integer,optional,intent(in) :: ilu(2)
830 : real(dp),intent(in) :: kpoint(3)
831 : real(dp),optional,intent(in) :: vlu(2)
832 : real(dp),intent(in) :: gmet(3,3)
833 :
834 : !Local variables-------------------------------
835 : !scalars
836 : integer :: npw_k
837 : logical :: ltest
838 : character(len=500) :: msg
839 10 : type(MPI_type) :: mpi_enreg_seq
840 : !arrays
841 10 : integer,allocatable :: kg_k(:,:)
842 : ! *************************************************************************
843 :
844 10 : call initmpi_seq(mpi_enreg_seq) ! Fake MPI_type.
845 :
846 10 : Dctl%spin = spin
847 10 : Dctl%nspinor = nspinor
848 40 : Dctl%kpoint = kpoint
849 :
850 10 : if (PRESENT(istwf_k)) then
851 10 : Dctl%istwf_k = istwf_k
852 : else
853 0 : Dctl%istwf_k = set_istwfk(kpoint)
854 : end if
855 :
856 10 : ABI_CHECK(Dctl%istwf_k == 1, "istwf_k/=1 not coded")
857 :
858 10 : Dctl%jobz = toupper(jobz(1:1))
859 10 : Dctl%range = "A"
860 10 : if (PRESENT(range)) Dctl%range = toupper(range)
861 :
862 10 : Dctl%ecut = ecut
863 10 : Dctl%ecutsm = zero; if (PRESENT(ecutsm)) Dctl%ecutsm = ecutsm
864 10 : Dctl%effmass_free = one; if (PRESENT(effmass_free)) Dctl%effmass_free = effmass_free
865 40 : Dctl%nloalg = nloalg
866 10 : Dctl%prtvol = 0; if (PRESENT(prtvol)) Dctl%prtvol = prtvol
867 10 : Dctl%abstol = -tol8; if (PRESENT(abstol)) Dctl%abstol = abstol
868 :
869 10 : ABI_MALLOC(kg_k,(3,0))
870 :
871 : ! Total number of G-vectors for this k-point with istwf_k=1.
872 10 : call kpgsph(ecut,0,gmet,0,0,1,kg_k,kpoint,0,mpi_enreg_seq,0,Dctl%npwtot)
873 :
874 : ! G-vectors taking into account time-reversal symmetry.
875 10 : call kpgsph(ecut,0,gmet,0,0,istwf_k,kg_k,kpoint,0,mpi_enreg_seq,0,npw_k)
876 :
877 10 : Dctl%npw_k = npw_k
878 10 : ABI_FREE(kg_k)
879 :
880 10 : Dctl%do_full_diago = .FALSE.
881 :
882 10 : SELECT CASE (Dctl%range)
883 : CASE ("A")
884 :
885 : ! Check on the number of stored bands.
886 10 : Dctl%nband_k=-1
887 10 : if (PRESENT(nband_k)) Dctl%nband_k=nband_k
888 :
889 10 : if (Dctl%nband_k==-1.or.Dctl%nband_k>=npw_k*nspinor) then
890 0 : Dctl%nband_k=npw_k*nspinor
891 0 : write(msg,'(4a)')ch10,&
892 0 : 'Since the number of bands to be computed was (-1) or',ch10,&
893 0 : 'too large, it has been set to the max. value npw_k*nspinor. '
894 0 : if (Dctl%prtvol>0) call wrtout(std_out, msg)
895 : else
896 10 : Dctl%nband_k=nband_k
897 : end if
898 :
899 10 : Dctl%do_full_diago = (Dctl%nband_k==npw_k*nspinor)
900 :
901 10 : if (Dctl%do_full_diago) then
902 0 : write(msg,'(6a)')ch10,&
903 0 : 'Since the number of bands to be computed',ch10,&
904 0 : 'is equal to the number of G-vectors found for this k-point,',ch10,&
905 0 : 'the program will perform complete diagonalization.'
906 : else
907 10 : write(msg,'(6a)')ch10,&
908 10 : 'Since the number of bands to be computed',ch10,&
909 10 : 'is less than the number of G-vectors found,',ch10,&
910 20 : 'the program will perform partial diagonalization.'
911 : end if
912 10 : if (Dctl%prtvol>0) call wrtout(std_out, msg)
913 :
914 : CASE ("I")
915 0 : if (.not.PRESENT(ilu)) then
916 0 : ABI_ERROR(" ilu must be specified when range=I ")
917 : end if
918 0 : Dctl%ilu = ilu
919 :
920 0 : ltest = ( ( ilu(2)>=ilu(1) ) .and. ilu(1)>=1 .and. ilu(2)<=Dctl%npwtot )
921 0 : write(msg,'(a,2i0)')" Illegal value for ilu: ",ilu
922 0 : ABI_CHECK(ltest,msg)
923 0 : Dctl%nband_k= ilu(2)-ilu(1)+1
924 :
925 : CASE ("V")
926 0 : if (.not.PRESENT(vlu)) then
927 0 : ABI_ERROR(" vlu must be specified when range=V ")
928 : end if
929 0 : Dctl%vlu = vlu
930 :
931 0 : Dctl%nband_k=-1 !??
932 :
933 0 : ltest = (vlu(2)>vlu(1))
934 0 : write(msg,'(a,2f0.3)')" Illegal value for vlu: ",vlu
935 0 : ABI_CHECK(ltest,msg)
936 :
937 : CASE DEFAULT
938 10 : ABI_ERROR(" Unknown value for range: "//TRIM(Dctl%range))
939 : END SELECT
940 :
941 : ! Consider the case in which we asked for the entire set of eigenvectors
942 : ! but the number of bands is less that npw_k. Therefore have to prepare the call to ZHEEVX.
943 : ! TODO this has to be done in a cleaner way.
944 10 : if (Dctl%range == "A" .and. .not. dctl%do_full_diago) then
945 10 : Dctl%range="I"
946 10 : Dctl%ilu(1) = 1
947 10 : Dctl%ilu(2) = npw_k*nspinor
948 10 : Dctl%nband_k= npw_k*nspinor
949 : end if
950 :
951 10 : Dctl%use_scalapack=0
952 10 : if (PRESENT(use_scalapack)) Dctl%use_scalapack=use_scalapack
953 10 : ABI_CHECK(Dctl%use_scalapack==0," scalapack mode not coded yet")
954 :
955 10 : call destroy_mpi_enreg(mpi_enreg_seq)
956 :
957 10 : end subroutine init_ddiago_ctl
958 : !!***
959 :
960 : !!****f* m_ksdiago/ugb_from_diago
961 : !! NAME
962 : !! ugb_from_diago
963 : !!
964 : !! FUNCTION
965 : !! This routine performs the direct diagonalization of the Kohn-Sham Hamiltonian
966 : !! for a given k-point and spin using Scalapack/ELPA.
967 : !!
968 : !! INPUTS
969 : !! spin= spin index.
970 : !! istwf_k= Storage mode for wavefunctions.
971 : !! kpoint(3)= k-point in reduced coordinates
972 : !! ecut= Cutoff energy
973 : !! gs_fermie=Fermi level as computed from the previous GS run.
974 : !! nband_k=Number of bands
975 : !! prtvol=Verbosity level
976 : !! nfftf=(effective) number of FFT grid points in the dense FFT mesh (for this processor)
977 : !! (nfftf=nfft for norm-conserving potential runs)
978 : !! pawtab(psps%ntypat*psps%usepaw) <type(pawtab_type)>=paw tabulated starting data
979 : !! pawfgr<pawfgr_type>=fine grid parameters and related data
980 : !! paw_ij(natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
981 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
982 : !! vtrial(nfftf,nspden)=the trial potential
983 : !! comm=MPI communicator.
984 : !! nfftc=Number of points in the coarse FFT mesh.
985 : !! ngfftc(18)=Info about 3D FFT for the coarse mesh, see ~abinit/doc/variables/vargs.htm#ngfft
986 : !! [electronpositron] <electronpositron_type>=quantities for the electron-positron annihilation.
987 : !!
988 : !! OUTPUT
989 : !! eig_k(1:nband_k)=The calculatated eigenvalues in ascending order.
990 : !!
991 : !! SOURCE
992 :
993 0 : subroutine ugb_from_diago(ugb, spin, istwf_k, kpoint, ecut, gs_fermie, nband_k, ngfftc, nfftf, &
994 0 : dtset, pawtab, pawfgr, paw_ij, cryst, psps, vtrial, eig_k, hyb, comm, &
995 : electronpositron) ! Optional arguments
996 :
997 : !Arguments ------------------------------------
998 : !scalars
999 : class(ugb_t),target,intent(out) :: ugb
1000 : integer,intent(in) :: spin, istwf_k
1001 : real(dp),intent(in) :: kpoint(3), ecut, gs_fermie
1002 : type(dataset_type),intent(in) :: dtset
1003 : integer,intent(in) :: comm,nfftf
1004 : integer,intent(inout) :: nband_k
1005 : type(crystal_t),intent(in) :: cryst
1006 : type(pseudopotential_type),intent(in) :: psps
1007 : type(pawfgr_type),intent(in) :: pawfgr
1008 : type(hyb_t),intent(inout) :: hyb
1009 : !arrays
1010 : integer,intent(in) :: ngfftc(18)
1011 : real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden)
1012 : !real(dp),intent(inout) :: vxctau(nfftf, dtset%nspden, 4*usevxctau)
1013 : real(dp),allocatable,intent(out) :: eig_k(:)
1014 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
1015 : type(paw_ij_type),intent(in) :: paw_ij(cryst%natom*psps%usepaw)
1016 : type(electronpositron_type),optional,pointer :: electronpositron
1017 :
1018 : !Local variables-------------------------------
1019 : !scalars
1020 : integer,parameter :: mkmem1 = 1, tim_getghc = 4, paral_kgb0 = 0, master = 0, ncomp1 = 1
1021 : integer :: cprj_choice,cpopt,dimffnl,ib,ider,idir,npw_k,nfftc,mgfftc, igs, ige, omp_nt
1022 : integer :: jj,n1,n2,n3,n4,n5,n6,nkpg,nproc,my_rank,optder, ib_glob, nb_glob, ib_loc
1023 : integer :: type_calc,sij_opt,igsp2_start,ig, my_ib, ibs1, ipwsp, islice, igsp_loc, my_npwsp
1024 : integer :: npwsp, col_bsize, nsppol, nspinor, nspden, loc2_size, il_g1, il_g2, ig1, ig2, ierr, min_my_nband, band_sum
1025 : integer :: idat, ndat, batch_size, h_size !, mene_found
1026 : integer :: ik_ibz, ik_bz, isym_k, trev_k, g0_k(3), g0(3)
1027 : integer :: iq_bz !, isym_q, trev_q !, g0_q(3) iq_ibz
1028 : real(dp),parameter :: lambda0 = zero
1029 : real(dp) :: cpu, wall, gflops, mem_mb, f_bsum, fact_spin, tol_empty_in, tol_empty, gsq_max, inv_sqrt_ucvol, rcut
1030 : logical :: do_full_diago, haveit, isirr_k, q_is_gamma ! isirr_q,
1031 : character(len=80) :: frmt1
1032 : character(len=10) :: stag(2)
1033 : character(len=500) :: msg
1034 0 : type(MPI_type) :: mpi_enreg_seq
1035 0 : type(gs_hamiltonian_type) :: gs_hamk
1036 0 : type(slkmat_dp_t) :: ghg_mat, gsg_mat, ghg_4diag, gsg_4diag, eigvec
1037 : type(slk_processor_t) :: proc_1d, proc_4diag
1038 0 : type(uplan_t) :: uplan_k
1039 : type(fftbox_plan3_t) :: box_plan
1040 0 : type(psbands_t) :: psb
1041 : !arrays
1042 0 : integer,allocatable :: gfft(:,:)
1043 : real(dp) :: kptns_(3,1), ylmgr_dum(1,1,1), tsec(2), ksum(3), kk_ibz(3), kgw_m_ksum(3), qq_bz(3), my_gw_qlwl(3) ! q0(3),
1044 0 : real(dp),allocatable :: ph3d(:,:,:), ffnl(:,:,:,:), kinpw(:), kpg_k(:,:), thetas(:,:)
1045 0 : real(dp),allocatable :: vlocal(:,:,:,:), ylm_k(:,:), dum_ylm_gr_k(:,:,:), eig_ene(:), ghc(:,:), gvnlxc(:,:), gsc(:,:), vcg_qbz(:,:)
1046 0 : real(dp),target,allocatable :: bras(:,:)
1047 0 : complex(dp),allocatable :: ps_ug(:,:,:)
1048 0 : complex(gwp),allocatable :: cbras_box(:,:), cbras_g(:,:), vc_sqrt(:), ur(:), rfg_box(:,:)
1049 0 : type(pawcprj_type),allocatable :: cwaveprj(:,:)
1050 : ! *********************************************************************
1051 :
1052 0 : call timab(1919, 1, tsec)
1053 0 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1054 :
1055 : ! See sequence of calls in vtorho.
1056 : ! Check that usekden is not 0 if want to use vxctau
1057 : !with_vxctau = dtset%usekden/=0
1058 :
1059 0 : if (dtset%usekden/=0) then
1060 0 : ABI_ERROR("nscf_init with mgga not yet coded")
1061 : end if
1062 : ! Check if want to use vxctau
1063 : !with_vxctau = (present(vxctau).and.usevxctau/=0)
1064 :
1065 : !====================
1066 : !=== Check input ====
1067 : !====================
1068 0 : if (all(istwf_k /= [1, 2])) then
1069 0 : ABI_ERROR(sjoin("istwfk:", itoa(istwf_k), "not allowed:"))
1070 : end if
1071 :
1072 0 : if (istwf_k == 2) then
1073 0 : ABI_ERROR("istwfk == 2 with direct diago is still under development")
1074 : !ABI_WARNING("istwfk == 2 with direct diago is still under development")
1075 : end if
1076 :
1077 0 : if (dtset%ixc < 0) then
1078 0 : if (libxc_functionals_ismgga() .and. .not. libxc_functionals_is_potential_only()) then
1079 0 : ABI_ERROR("meta-gga functionals are not compatible with direct diagonalization!")
1080 : end if
1081 : end if
1082 :
1083 : ! MPI_type for sequential part.
1084 0 : call initmpi_seq(mpi_enreg_seq)
1085 0 : call mpi_enreg_seq%distribfft%init_seq('c', ngfftc(2), ngfftc(3), 'all')
1086 0 : if (pawfgr%usefinegrid /= 0) then
1087 0 : call mpi_enreg_seq%distribfft%init_seq('f', pawfgr%ngfft(2), pawfgr%ngfft(3), 'all')
1088 : end if
1089 :
1090 0 : nspinor = dtset%nspinor; nsppol = dtset%nsppol; nspden = dtset%nspden
1091 0 : if (nsppol == 1) stag = [' ',' ']
1092 0 : if (nsppol == 2) stag = ['SPIN UP: ','SPIN DOWN:']
1093 :
1094 : ! Get g-vectors from kpt and ecut.
1095 0 : call get_kg(kpoint, istwf_k, ecut, cryst%gmet, npw_k, ugb%kg_k)
1096 0 : npwsp = npw_k * nspinor
1097 :
1098 : ! The coarse FFT mesh for the application of the Hamiltonian.
1099 0 : n1 = ngfftc(1); n2 = ngfftc(2); n3 = ngfftc(3)
1100 0 : n4 = ngfftc(4); n5 = ngfftc(5); n6 = ngfftc(6)
1101 0 : nfftc = product(ngfftc(1:3)); mgfftc = maxval(ngfftc(1:3))
1102 :
1103 : ! Initialize the Hamiltonian on the coarse FFT mesh.
1104 0 : if (present(electronpositron)) then
1105 : call gs_hamk%init(psps, pawtab, nspinor, nsppol, nspden, cryst%natom, cryst%typat, cryst%xred, nfftc, &
1106 0 : mgfftc, ngfftc, cryst%rprimd, dtset%nloalg, paw_ij=paw_ij, usecprj=0, gpu_option=dtset%gpu_option, electronpositron=electronpositron)
1107 : else
1108 : call gs_hamk%init(psps, pawtab, nspinor, nsppol, nspden, cryst%natom, cryst%typat, cryst%xred, nfftc, &
1109 0 : mgfftc, ngfftc, cryst%rprimd, dtset%nloalg, paw_ij=paw_ij, usecprj=0, gpu_option=dtset%gpu_option)
1110 : end if
1111 :
1112 : ! Check on the number of stored bands.
1113 0 : if (nband_k == -1 .or. nband_k >= npwsp) then
1114 0 : nband_k = npwsp
1115 0 : write(msg,'(4a, i0)')ch10,&
1116 0 : ' Since the number of bands to be computed was -1 or',ch10,&
1117 0 : ' too large, it has been set to the maximum value. npw_k*nspinor: ',npwsp
1118 0 : call wrtout(std_out, msg)
1119 : end if
1120 :
1121 0 : do_full_diago = nband_k == npwsp
1122 :
1123 : ! Set up local potential vlocal with proper dimensioning, from vtrial.
1124 : ! Select spin component of interest if nspden<=2 as nvloc==1, for nspden==4, nvloc==4
1125 : ! option=2: vtrial(n1*n2*n3,ispden) --> vlocal(nd1,nd2,nd3) real case
1126 :
1127 0 : ABI_MALLOC(vlocal, (n4, n5, n6, gs_hamk%nvloc))
1128 : call gspot_transgrid_and_pack(spin, psps%usepaw, paral_kgb0, nfftc, ngfftc, nfftf, &
1129 0 : nspden, gs_hamk%nvloc, ncomp1, pawfgr, mpi_enreg_seq, vtrial, vlocal)
1130 0 : call gs_hamk%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
1131 :
1132 : ! TODO: This for meta-gga.
1133 : !if (with_vxctau) then
1134 : ! call gspot_transgrid_and_pack(spin, psps%usepaw, paral_kgb0, nfftc, ngfftc, nfftf, &
1135 : ! nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal)
1136 : ! call gs_hamk%load_spin(spin, vxctaulocal=vxctaulocal)
1137 : !end if
1138 :
1139 : !========================
1140 : !==== Kinetic energy ====
1141 : !========================
1142 0 : ABI_MALLOC(kinpw, (npw_k))
1143 0 : call mkkin(ecut, dtset%ecutsm, dtset%effmass_free, cryst%gmet, ugb%kg_k, kinpw, kpoint, npw_k, 0, 0)
1144 :
1145 : !================================
1146 : !==== Non-local form factors ====
1147 : !================================
1148 0 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2*psps%useylm))
1149 :
1150 0 : if (psps%useylm == 1) then
1151 0 : optder = 0
1152 0 : ABI_MALLOC(dum_ylm_gr_k, (npw_k, 3+6*(optder/2),psps%mpsang**2))
1153 0 : kptns_(:,1) = kpoint
1154 :
1155 : ! NB: Here mband is not used if paral_compil_kpt = 0
1156 : call initylmg(cryst%gprimd, ugb%kg_k, kptns_, mkmem1, mpi_enreg_seq, psps%mpsang, npw_k, [nband_k], 1, &
1157 0 : [npw_k], 1, optder, cryst%rprimd, ylm_k, dum_ylm_gr_k)
1158 :
1159 0 : ABI_FREE(dum_ylm_gr_k)
1160 : end if
1161 :
1162 : ! Compute (k+G) vectors (only if useylm=1)
1163 0 : nkpg = 3 * dtset%nloalg(3)
1164 0 : ABI_MALLOC(kpg_k, (npw_k, nkpg))
1165 0 : if (nkpg > 0) call mkkpg(ugb%kg_k, kpg_k, kpoint, nkpg, npw_k)
1166 :
1167 : ! Compute nonlocal form factors ffnl at all (k+G):
1168 0 : idir=0; ider=0; dimffnl=1+ider ! Now the derivative is not needed anymore.
1169 0 : ABI_MALLOC(ffnl, (npw_k, dimffnl, psps%lmnmax, psps%ntypat))
1170 :
1171 : call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl, psps%ffspl, cryst%gmet, cryst%gprimd, ider, idir, psps%indlmn, &
1172 : ugb%kg_k, kpg_k, kpoint, psps%lmnmax, psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw_k, &
1173 0 : psps%ntypat, psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm_k, ylmgr_dum)
1174 :
1175 0 : ABI_FREE(ylm_k)
1176 :
1177 : ! Load k-dependent part of the Hamiltonian.
1178 0 : ABI_MALLOC(ph3d, (2, npw_k, gs_hamk%matblk))
1179 : call gs_hamk%load_k(kpt_k=kpoint, istwf_k=istwf_k, npw_k=npw_k, kinpw_k=kinpw, &
1180 0 : kg_k=ugb%kg_k, kpg_k=kpg_k, ffnl_k=ffnl, ph3d_k=ph3d, compute_ph3d=.true., compute_gbound=.true.)
1181 :
1182 : ! Prepare call to getghc.
1183 0 : type_calc = 0 ! For applying the whole Hamiltonian
1184 0 : sij_opt = 0; if (psps%usepaw==1) sij_opt = 1 ! For PAW, <k+G|S|k+G"> is also needed.
1185 :
1186 0 : cpopt = -1
1187 0 : if (psps%usepaw == 1) cpopt = 0 ! <p_lmn|in> are computed here and saved
1188 :
1189 : ! Init 1D PBLAS grid to block-distribute H along columns.
1190 0 : call proc_1d%init(comm, grid_dims=[1, nproc])
1191 0 : h_size = npwsp; if (istwf_k == 2) h_size = 2*npwsp - 1
1192 :
1193 0 : ABI_CHECK(block_dist_1d(h_size, nproc, col_bsize, msg), msg)
1194 0 : call ghg_mat%init(h_size, h_size, proc_1d, istwf_k, size_blocs=[h_size, col_bsize])
1195 0 : if (psps%usepaw == 1) call gsg_mat%init(h_size, h_size, proc_1d, istwf_k, size_blocs=[h_size, col_bsize])
1196 :
1197 : ! Estimate memory
1198 0 : mem_mb = ghg_mat%locmem_mb()
1199 0 : mem_mb = two * (psps%usepaw + 1) * mem_mb + mem_mb ! last term for eigvec matrix
1200 0 : call wrtout(std_out, sjoin(" Local memory for scalapack matrices:", ftoa(mem_mb, fmt="(f8.1)"), ' [Mb] <<< MEM'))
1201 :
1202 : ! Define batch size for the application of the Hamiltonian
1203 : ! This is useful if OpenMP is activated thus we use multiples of omp_nt.
1204 0 : omp_nt = xomp_get_num_threads(open_parallel=.True.)
1205 0 : batch_size = 8 * omp_nt
1206 0 : if (istwf_k == 2) batch_size = 1 ! FIXME
1207 : !batch_size = 1
1208 0 : if (gs_hamk%gpu_option == ABI_GPU_OPENMP) then
1209 0 : batch_size = 32
1210 : end if
1211 :
1212 0 : call wrtout(std_out, sjoin(" Building H^KS with batch_size:", itoa(batch_size)))
1213 :
1214 0 : ABI_MALLOC(bras, (2, npwsp * batch_size))
1215 : ! cwaveprj is ordered by atom type, see nonlop_ylm.
1216 0 : ABI_MALLOC(cwaveprj, (cryst%natom, nspinor*(1+cpopt)*gs_hamk%usepaw*batch_size))
1217 0 : if (cpopt == 0) call pawcprj_alloc(cwaveprj, 0, gs_hamk%dimcprj)
1218 0 : ABI_MALLOC(ghc, (2, npwsp * batch_size))
1219 0 : ABI_MALLOC(gvnlxc, (2, npwsp * batch_size))
1220 0 : ABI_MALLOC(gsc, (2, npwsp * batch_size*(sij_opt+1)/2))
1221 :
1222 : ! Loop over the |beta,G''> component.
1223 0 : call cwtime(cpu, wall, gflops, "start")
1224 0 : loc2_size = ghg_mat%size_local(2)
1225 :
1226 0 : if (my_rank == master) call pstat_proc%print(_PSTAT_ARGS_)
1227 :
1228 0 : do il_g2=1, loc2_size, batch_size
1229 : ! Operate on ndat g-vectors starting at the igsp2_start global index.
1230 0 : igsp2_start = ghg_mat%loc2gcol(il_g2)
1231 0 : ndat = blocked_loop(il_g2, loc2_size, batch_size)
1232 :
1233 0 : bras = zero
1234 0 : if (istwf_k == 1) then
1235 0 : do idat=0,ndat-1
1236 0 : bras(1, igsp2_start + idat * npwsp + idat) = one
1237 : end do
1238 : else
1239 : ! only istwf_k == 2 is coded here. NB: there's a check at the beginning of this routine.
1240 0 : do idat=0,ndat-1
1241 0 : if (igsp2_start + idat <= npwsp) then
1242 : ! Cosine term
1243 0 : bras(1, igsp2_start + idat*npwsp + idat) = half
1244 0 : if (igsp2_start == 1) bras(1, igsp2_start + idat*npwsp + idat) = one
1245 : else
1246 : ! Sine term
1247 : !ig = igsp2_start - npwsp + 1
1248 0 : ig = igsp2_start - npwsp + 1 + 1 ! This should be OK
1249 0 : bras(2, ig + idat*npwsp + idat) = half
1250 : end if
1251 : end do
1252 : end if
1253 :
1254 : ! Get <:|H|beta,G''> and <:|S_{PAW}|beta,G''>
1255 : call multithreaded_getghc(cpopt, bras, cwaveprj, ghc, gsc, gs_hamk, gvnlxc, lambda0, mpi_enreg_seq, ndat, &
1256 0 : dtset%prtvol, sij_opt, tim_getghc, type_calc)
1257 :
1258 : ! Now fill my local buffer of ghg/gsg.
1259 0 : if (istwf_k == 1) then
1260 : ! Complex wavefunctions.
1261 0 : do idat=0,ndat-1
1262 0 : igs = 1 + idat * npwsp; ige = igs + npwsp - 1
1263 0 : ghg_mat%buffer_cplx(:, il_g2+idat) = cmplx(ghc(1, igs:ige), ghc(2, igs:ige), kind=dp)
1264 : end do
1265 0 : if (psps%usepaw == 1) then
1266 0 : do idat=0,ndat-1
1267 0 : igs = 1 + idat * npwsp; ige = igs + npwsp - 1
1268 0 : gsg_mat%buffer_cplx(:, il_g2+idat) = cmplx(gsc(1,igs:ige), gsc(2,igs:ige), kind=dp)
1269 : end do
1270 : end if
1271 :
1272 : else
1273 : ! Real wavefunctions.
1274 0 : do idat=0,ndat-1
1275 0 : igs = 1 + idat*npwsp; ige = igs + npwsp - 1
1276 : !if (igsp2_start == 1 .or. igsp2_start == npwsp + 1 .and. idat == 0) then
1277 : ! ghc(:, igs:ige) = tol3 !; print *, ghc(:, igs:ige)
1278 : !end if
1279 0 : ghg_mat%buffer_real(1:npwsp, il_g2+idat) = ghc(1, igs:ige) ! CC or CS
1280 0 : ghg_mat%buffer_real(npwsp+1:, il_g2+idat) = -ghc(2, igs+1:ige) ! SC or SS. Note igs+1
1281 : end do
1282 0 : if (psps%usepaw == 1) then
1283 0 : NOT_IMPLEMENTED_ERROR()
1284 : !gsg_mat%buffer_real(...)
1285 : end if
1286 : end if ! istwf_k
1287 : end do ! il_g2
1288 :
1289 : ! MG: DEBUG
1290 : !call wrtout(std_out, " WARNING: Setting H_KS to zero for debugging purposes!"); ghg_mat%buffer_cplx = czero
1291 0 : call cwtime_report(" build H^KS_g1g2", cpu, wall, gflops)
1292 :
1293 : ! Free workspace memory allocated so far.
1294 0 : ABI_FREE(bras)
1295 0 : ABI_FREE(kinpw)
1296 0 : ABI_FREE(vlocal)
1297 0 : ABI_FREE(ghc)
1298 0 : ABI_FREE(gvnlxc)
1299 0 : ABI_FREE(gsc)
1300 0 : if (psps%usepaw == 1 .and. cpopt == 0) call pawcprj_free(cwaveprj)
1301 0 : ABI_FREE(cwaveprj)
1302 :
1303 : ! ==================================
1304 : ! Compute Fock operator F^k_{g1,g2}
1305 : ! ==================================
1306 0 : if (dtset%usefock == 1) then
1307 : !if (.False.) then
1308 0 : call cwtime(cpu, wall, gflops, "start")
1309 0 : call wrtout(std_out, sjoin(" Building Fock operator F^k_{g1,g2} with batch_size:", itoa(batch_size)))
1310 0 : ABI_CHECK(dtset%usepaw == 0, "DIRECT DIAGO OF FOCK OPERATOR WITH PAW IS NOT CODED!")
1311 0 : inv_sqrt_ucvol = one / sqrt(cryst%ucvol)
1312 :
1313 0 : call hyb%wfd%change_ngfft(cryst, psps, ngfftc)
1314 :
1315 0 : ABI_MALLOC(ur, (nfftc*nspinor))
1316 0 : ABI_MALLOC(cbras_g, (npw_k*nspinor, batch_size))
1317 0 : ABI_MALLOC(cbras_box, (nfftc*nspinor, batch_size))
1318 0 : ABI_MALLOC(rfg_box, (nfftc*nspinor, batch_size))
1319 0 : ABI_MALLOC(vc_sqrt, (nfftc))
1320 0 : ABI_MALLOC(vcg_qbz, (nfftc, hyb%nqbz))
1321 :
1322 : ! Set tolerance used to decide if a band is empty.
1323 0 : tol_empty_in = 0.01_dp
1324 0 : call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
1325 :
1326 : ! Precompute the Coulomb term here to avoid tons of calls inside the loop over ig2.
1327 : ! Get g-vectors in the FFT box for vcoul.
1328 0 : ABI_MALLOC(gfft, (3, nfftc))
1329 0 : call get_gfft(ngfftc, kpoint, cryst%gmet, gsq_max, gfft)
1330 :
1331 0 : my_gw_qlwl(:) = GW_Q0_DEFAULT; if (dtset%gw_nqlwl > 0) my_gw_qlwl = dtset%gw_qlwl(:,1)
1332 : !my_gw_qlwl = zero
1333 0 : do ik_bz=1,hyb%nkbz
1334 0 : ksum = hyb%kbz(:, ik_bz)
1335 0 : kgw_m_ksum = kpoint - ksum
1336 : !print *, "kpoint", kpoint, "ksum:", ksum
1337 0 : call findqg0(iq_bz, g0, kgw_m_ksum, hyb%nqbz, hyb%qbz, hyb%mG0)
1338 0 : ABI_CHECK(all(g0 == 0), sjoin("g0 = ", ltoa(g0)))
1339 0 : qq_bz = hyb%qbz(:,iq_bz)
1340 0 : q_is_gamma = normv(qq_bz, cryst%gmet, "G") < GW_TOLQ0
1341 0 : call hyb%vcgen%get_vc_sqrt(qq_bz, nfftc, gfft, my_gw_qlwl, cryst, vc_sqrt, comm, vc=vcg_qbz(:,iq_bz))
1342 : ! A non-positive value of rcut activates the recipe of Spencer & Alavi, PRB 77, 193110 (2008) [[cite:Spencer2008]].
1343 0 : rcut = (cryst%ucvol * hyb%nkbz * 3.d0 / four_pi) ** third
1344 : !vcgen%i_sz = two_pi * rcut**2
1345 0 : if (q_is_gamma) then
1346 : !vcg_qbz(1,iq_bz) = two_pi * rcut**2 ! FIXME: This is used in GW
1347 : !vcg_qbz(1,iq_bz) = hyb%vcgen%i_sz
1348 0 : vcg_qbz(1,iq_bz) = two_pi/three * rcut**2 ! FIXME: This is used in m_fock
1349 : !vcg_qbz(1,iq_bz) = zero
1350 : end if
1351 : !vcg_qbz(2:,iq_bz) = vcg_qbz(2:,iq_bz) * (inv_sqrt_ucvol**2)
1352 0 : vcg_qbz(:,iq_bz) = vcg_qbz(:,iq_bz) * (inv_sqrt_ucvol**2)
1353 0 : vcg_qbz(:,iq_bz) = one
1354 0 : call zerosym(vcg_qbz(:,iq_bz), 1, n1, n2, n3)
1355 : end do ! ik_bz
1356 0 : ABI_FREE(gfft)
1357 :
1358 : ! Build plans for (dense, g-sphere) FFTs.
1359 0 : call box_plan%from_ngfft(ngfftc, nspinor*batch_size, dtset%gpu_option)
1360 0 : call uplan_k%init(npw_k, nspinor, batch_size, ngfftc, istwf_k, ugb%kg_k, gwp, dtset%gpu_option)
1361 :
1362 : ! Blocked loop over the columns of F^k_{g1,g2}.
1363 0 : do ig2=1, npwsp, batch_size
1364 0 : ndat = blocked_loop(ig2, npwsp, batch_size)
1365 : ! Fill cbras_box(r) with e^{ig2.r}.
1366 0 : do idat=1,ndat
1367 0 : call calc_ceigr(ugb%kg_k(:,ig2+idat-1), nfftc, nspinor, ngfftc, cbras_box(:,idat))
1368 0 : cbras_box(:,idat) = cbras_box(:,idat) * inv_sqrt_ucvol
1369 : end do
1370 :
1371 : ! ==============================
1372 : ! ==== Sum over k in the BZ ====
1373 : ! ==============================
1374 0 : rfg_box = zero
1375 0 : do ik_bz=1,hyb%nkbz
1376 0 : ksum = hyb%kbz(:, ik_bz)
1377 : ! Parallelism over k-points.
1378 : !if (.not. hyb%wfd%ihave_ug(0, ik_ibz, spin)) cycle
1379 :
1380 : ! Find the symmetrical image of ksum in the IBZ
1381 : ! FIXME: Be careful with the symmetry conventions here and the interplay between umklapp in q and FFT
1382 0 : ik_ibz = hyb%kbz2ibz_symrel(1, ik_bz); isym_k = hyb%kbz2ibz_symrel(2, ik_bz)
1383 0 : trev_k = hyb%kbz2ibz_symrel(6, ik_bz); g0_k = hyb%kbz2ibz_symrel(3:5, ik_bz)
1384 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
1385 0 : kk_ibz = hyb%kibz(:, ik_ibz)
1386 :
1387 : ! Identify q and G0 where q + G0 = k_GW - ksum
1388 : !kgw_m_ksum = kpoint - ksum
1389 : !call findqg0(iq_bz, g0, kgw_m_ksum, hyb%nqbz, hyb%qbz, hyb%mG0)
1390 : !ABI_CHECK(all(g0 == 0), sjoin("g0 = ", ltoa(g0)))
1391 :
1392 : !qq_bz = hyb%qbz(:, iq_bz)
1393 : !iq_ibz = hyb%qbz2ibz(1, iq_bz); isym_q = hyb%qbz2ibz(2, iq_bz)
1394 : !trev_q = hyb%qbz2ibz(6, iq_bz); g0_q = hyb%qbz2ibz(3:5, iq_bz)
1395 : !isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0))
1396 :
1397 : !! Find the corresponding irreducible q-point.
1398 : !! NB: non-zero umklapp G_o is not allowed. There's a check in setup_sigma
1399 : !!call qmesh%get_BZ_item(iq_bz, qbz, iq_ibz, isym_q, itim_q)
1400 : !q_is_gamma = normv(qq_bz, cryst%gmet, "G") < GW_TOLQ0
1401 :
1402 : ! ==========================
1403 : ! Sum over (occupied) bands
1404 : ! ==========================
1405 0 : do band_sum=1, hyb%wfd%nband(ik_ibz, spin)
1406 : ! MPI parallelism over bands.
1407 0 : if (.not. hyb%wfd%ihave_ug(band_sum, ik_ibz, spin)) cycle
1408 0 : f_bsum = hyb%ebands%occ(band_sum, ik_ibz, spin) * fact_spin; if (abs(f_bsum) <= tol_empty) cycle
1409 :
1410 : !print *, "band_sum, ik_ibz, spin", band_sum, ik_ibz, spin
1411 0 : call hyb%wfd%get_ur(band_sum, ik_ibz, spin, ur)
1412 0 : ur = ur * inv_sqrt_ucvol
1413 0 : do idat=1,ndat
1414 0 : cbras_box(:,idat) = f_bsum * conjg(ur) * cbras_box(:,idat)
1415 : end do
1416 :
1417 : ! FFT r --> g and multiply by v(g,q) on the FFT box.
1418 0 : call box_plan%execute(cbras_box(:,1), -1, ndat=ndat)
1419 0 : do idat=1,ndat
1420 0 : cbras_box(:,idat) = cbras_box(:,idat) * vcg_qbz(:, iq_bz)
1421 : end do
1422 : ! FFT g --> r, multiply by u(r) and accumulate in rfg_box
1423 0 : call box_plan%execute(cbras_box(:,1), +1, ndat=ndat)
1424 0 : do idat=1,ndat
1425 0 : rfg_box(:,idat) = rfg_box(:,idat) + cbras_box(:,idat) * ur
1426 : end do
1427 : end do ! band_sum
1428 : end do ! ik_bz
1429 :
1430 : ! FFT r --> g_sphere and MPI sum partial contributions.
1431 0 : call uplan_k%execute_rg(ndat, rfg_box(:,1), cbras_g(:,1))
1432 0 : call xmpi_sum(cbras_g, hyb%wfd%comm_spin(spin), ierr)
1433 0 : cbras_g = - cbras_g * sqrt(cryst%ucvol)
1434 : !cbras_g = - half * cbras_g * sqrt(cryst%ucvol)
1435 : !cbras_g = - 10000 * cbras_g / (hyb%nkbz*cryst%ucvol) ! * alpha_hyb
1436 : !cbras_g = - sqrt(cryst%ucvol) * cbras_g ! / (hyb%nkbz*cryst%ucvol) ! * alpha_hyb
1437 : !cbras_g = -half * inv_sqrt_ucvol * cbras_g / (hyb%nkbz*cryst%ucvol) ! * alpha_hyb
1438 : !cbras_g = - cbras_g / (hyb%nkbz * cryst%ucvol) ! * alpha_hyb
1439 : !cbras_g = -half * sqrt(cryst%ucvol) * cbras_g !/ (hyb%nkbz*cryst%ucvol) ! * alpha_hyb
1440 :
1441 : ! Update my local buffer of ghg_mat.
1442 0 : do idat=1,ndat
1443 0 : do ig1=1,npwsp
1444 : ! From global to local indices.
1445 0 : call ghg_mat%glob2loc(ig1, ig2+idat-1, il_g1, il_g2, haveit); if (.not. haveit) cycle
1446 : !print *, "ig1, idat, cbras_g", ig1, idat, cbras_g(ig1, idat)
1447 0 : if (istwf_k == 1) then
1448 : ! Complex wavefunctions.
1449 0 : ghg_mat%buffer_cplx(il_g1, il_g2) = ghg_mat%buffer_cplx(il_g1, il_g2) + cbras_g(ig1,idat)
1450 : else
1451 : ! Real wavefunctions.
1452 0 : NOT_IMPLEMENTED_ERROR()
1453 : end if ! istwf_k
1454 : end do ! ig1
1455 : end do ! idat
1456 :
1457 : end do ! ig2
1458 :
1459 0 : ABI_FREE(vc_sqrt)
1460 0 : ABI_FREE(vcg_qbz)
1461 0 : ABI_FREE(cbras_box)
1462 0 : ABI_FREE(rfg_box)
1463 0 : ABI_FREE(cbras_g)
1464 0 : ABI_FREE(ur)
1465 0 : call uplan_k%free(); call box_plan%free()
1466 0 : call cwtime_report(" build Fock_g1g2", cpu, wall, gflops)
1467 : end if ! usefock
1468 :
1469 : !===========================================
1470 : !=== Diagonalization of <G|H|G''> matrix ===
1471 : !===========================================
1472 0 : ABI_MALLOC(eig_ene, (h_size))
1473 :
1474 : ! Change size block. Use 2D rectangular grid of processors for diagonalization, if possible.
1475 0 : call proc_4diag%init(comm)
1476 0 : call ghg_mat%change_size_blocs(ghg_4diag, processor=proc_4diag, free=.True.)
1477 0 : if (psps%usepaw == 1) call gsg_mat%change_size_blocs(gsg_4diag, processor=proc_4diag, free=.True.)
1478 : !call ghg_mat%copy(ghg_4diag); call ghg_mat%free()
1479 :
1480 : ! NB: global H shape is (h_size, h_size) even for partial diago.
1481 : ! then one extracts the (hsize, nband_k) sub-matrix before returning.
1482 0 : call ghg_4diag%copy(eigvec)
1483 0 : if (my_rank == master) call pstat_proc%print(_PSTAT_ARGS_)
1484 :
1485 : #ifndef HAVE_LINALG_ELPA
1486 : call wrtout([std_out, ab_out], &
1487 0 : "- WARNING: Using ScaLAPACK for diagonalization, but ELPA library is highly recommended for both efficiency and memory reasons.")
1488 : #endif
1489 :
1490 0 : if (do_full_diago) then
1491 0 : write(msg,'(5a, (a,i0), 2a)')ch10,&
1492 0 : ' Begin full diagonalization for kpt: ',trim(ktoa(kpoint)), stag(spin), ch10,&
1493 0 : " H_gg' Matrix size: ",npwsp, ", Scalapack grid: ", trim(ltoa(ghg_4diag%processor%grid%dims))
1494 0 : call wrtout(std_out, msg)
1495 0 : call cwtime(cpu, wall, gflops, "start")
1496 0 : if (psps%usepaw == 0) then
1497 : !call ghg_4diag%pzheev("V", "U", eigvec, eig_ene)
1498 0 : call compute_eigen_problem(ghg_4diag%processor, ghg_4diag, eigvec, eig_ene, comm, istwf_k)
1499 : else
1500 0 : call compute_generalized_eigen_problem(ghg_4diag%processor, ghg_4diag, gsg_4diag, eigvec, eig_ene, comm, istwf_k)
1501 : end if
1502 0 : call cwtime_report(" full_diago", cpu, wall, gflops)
1503 :
1504 : else
1505 0 : write(msg,'(6a,i0,(a,i0), 2a)') ch10,&
1506 0 : ' Begin partial diagonalization for kpt: ',trim(ktoa(kpoint)), stag(spin), ch10,&
1507 0 : " H_gg' Matrix size: ",npwsp,', nband_k: ', nband_k,", Scalapack grid: ", trim(ltoa(ghg_4diag%processor%grid%dims))
1508 0 : call wrtout(std_out, msg)
1509 :
1510 0 : call cwtime(cpu, wall, gflops, "start")
1511 0 : if (psps%usepaw == 0) then
1512 : !call ghg_4diag%pzheevx("V", "I", "U", zero, zero, 1, nband_k, -tol8, eigvec, mene_found, eig_ene)
1513 0 : call compute_eigen_problem(ghg_4diag%processor, ghg_4diag, eigvec, eig_ene, comm, istwf_k, nev=nband_k)
1514 : else
1515 : !call ghg_4diag%pzhegvx(1, "V", "I", "U", gsg_4diag, zero, zero, 1, nband_k, -tol8, eigvec, mene_found, eig_ene)
1516 : call compute_generalized_eigen_problem(ghg_4diag%processor, ghg_4diag, gsg_4diag, eigvec, eig_ene, comm, istwf_k, &
1517 0 : nev=nband_k)
1518 : end if
1519 0 : call cwtime_report(" partial_diago", cpu, wall, gflops)
1520 : end if
1521 :
1522 0 : if (my_rank == master) then
1523 0 : call pstat_proc%print(_PSTAT_ARGS_)
1524 : ! Write eigenvalues.
1525 0 : frmt1 = '(8x,*(1x,f7.3))'
1526 0 : write(msg, '(2a,3x,a)')' Eigenvalues in eV for kpt: ', trim(ktoa(kpoint)), stag(spin); call wrtout(std_out, msg)
1527 0 : write(msg, frmt1)(eig_ene(ib)*Ha_eV,ib=1,min(9,nband_k)); call wrtout(std_out, msg)
1528 : ! HYB DEBUG
1529 : !call wrtout(std_out, "hyb%ebands")
1530 : !write(msg, frmt1)(hyb%ebands%eig(ib,1,spin)*Ha_eV, ib=1,min(9,hyb%ebands%mband)); call wrtout(std_out, msg)
1531 0 : if (nband_k > 9 .and. dtset%prtvol > 0) then
1532 0 : do jj=10,nband_k,9
1533 0 : write(msg, frmt1) (eig_ene(ib)*Ha_eV,ib=jj,min(jj+8,nband_k)); call wrtout(std_out, msg)
1534 : end do
1535 : end if
1536 : end if
1537 :
1538 : ! Free memory
1539 0 : call ghg_4diag%free(); call gsg_4diag%free(); call proc_1d%free()
1540 :
1541 : ! ================
1542 : ! Stochastic bands
1543 : ! ================
1544 : !if (dtset%nb_protected /= 0) then
1545 : if (.False.) then
1546 : call wrtout(std_out, " Generating stochastic bands...")
1547 : ! Initial setup.
1548 : call psb%init(dtset, h_size, eig_ene, gs_fermie) !, nband_k)
1549 : my_npwsp = eigvec%size_local(1)
1550 : nb_glob = eigvec%size_global(2)
1551 : ABI_CALLOC(ps_ug, (my_npwsp, psb%maxsto_per_slice, psb%nslices))
1552 : ABI_MALLOC(thetas, (nb_glob, psb%maxsto_per_slice))
1553 :
1554 : ! Loop over global bands.
1555 : do ib_glob=1, nb_glob
1556 : ! Need the same random phases on all MPI procs.
1557 : if (eigvec%processor%my_rank == master) call random_number(thetas)
1558 : call xmpi_bcast(thetas, master, eigvec%processor%comm, ierr)
1559 :
1560 : ! Get slice index from ib_glob.
1561 : islice = psb%band2slice(ib_glob); if (islice == -1) cycle
1562 : !band_block = psb%subspace(1:2, islice)
1563 : !nb_in_slice = psb%subspace(3,islice)
1564 :
1565 : ! Loop over global PW index.
1566 : do ipwsp=1,npwsp
1567 : call eigvec%glob2loc(ipwsp, ib_glob, igsp_loc, ib_loc, haveit); if (.not. haveit) cycle
1568 : do ib=1,psb%subspace(3,islice)
1569 : ps_ug(igsp_loc, ib, islice) = ps_ug(igsp_loc, ib, islice) + &
1570 : eigvec%buffer_cplx(igsp_loc, ib_loc) * exp(j_dpc*two_pi*thetas(ib_glob,ib))
1571 : end do
1572 : end do
1573 : end do ! ib_glob
1574 : ABI_FREE(thetas)
1575 :
1576 : ! Normalize
1577 : ! TODO: Need MPI communicator over columns here.
1578 : !call xmpi_sum(ps_ug, eigvec%column_comm, ierr)
1579 : do islice=1,psb%nslices
1580 : do ib=1,psb%subspace(3,islice)
1581 : if (psb%subspace(3,islice) == 1) cycle
1582 : ps_ug(:,ib,islice) = ps_ug(:,ib,islice) / sqrt(one * psb%subspace(3,islice))
1583 : end do
1584 : end do
1585 :
1586 : ! Now insert ps_ug in the right position in eigevec
1587 : do ib_glob=1, nb_glob
1588 : islice = psb%band2slice(ib_glob); if (islice == -1) cycle
1589 : !band_start = 1 + (islice - 1) * psb%nb_per_slice
1590 : ! Loop over global PW index.
1591 : do ipwsp=1,npwsp
1592 : call eigvec%glob2loc(ipwsp, ib_glob, igsp_loc, ib_loc, haveit); if (.not. haveit) cycle
1593 : !do ib=1,psb%nb_per_slice
1594 : ! eigvec%buffer_cplx(igsp_loc, ib_loc) = ps_ug(igsp_loc, ib, islice)
1595 : !end do
1596 : end do
1597 : end do
1598 :
1599 : ABI_FREE(ps_ug)
1600 :
1601 : ! here we change the value of nband_k and eig_k.
1602 : nband_k = psb%nb_tot
1603 : ABI_MALLOC(eig_k, (nband_k))
1604 : eig_k = psb%ps_eig
1605 :
1606 : else
1607 : ! No pseudo bands.
1608 0 : ABI_MALLOC(eig_k, (nband_k))
1609 0 : eig_k(:) = eig_ene(1:nband_k)
1610 : end if
1611 :
1612 : ! Now transfer eigvec to the ugb datastructure using 1d grid (block column distribution).
1613 0 : call wrtout(std_out, " Moving to PBLAS block column distribution...")
1614 0 : call cwtime(cpu, wall, gflops, "start")
1615 :
1616 0 : call ugb%processor%init(comm, grid_dims=[1, nproc])
1617 0 : ABI_CHECK(block_dist_1d(nband_k, nproc, col_bsize, msg), msg)
1618 0 : call eigvec%cut(h_size, nband_k, ugb%mat, size_blocs=[h_size, col_bsize], processor=ugb%processor, free=.True.)
1619 0 : call proc_4diag%free()
1620 :
1621 : ! =================
1622 : ! Build ugb object
1623 : ! =================
1624 0 : ugb%istwf_k = istwf_k
1625 0 : ugb%nspinor = nspinor
1626 0 : ugb%npw_k = npw_k
1627 0 : ugb%npwsp = npwsp
1628 0 : ugb%nband_k = nband_k
1629 0 : ugb%comm => ugb%mat%processor%comm
1630 :
1631 0 : ugb%my_bstart = ugb%mat%loc2gcol(1)
1632 0 : ugb%my_bstop = ugb%mat%loc2gcol(ugb%mat%size_local(2))
1633 0 : ugb%my_nband = ugb%my_bstop - ugb%my_bstart + 1
1634 :
1635 0 : if (ugb%my_nband > 0) then
1636 0 : call c_f_pointer(c_loc(ugb%mat%buffer_cplx), ugb%cg_k, shape=[2, npwsp, ugb%my_nband])
1637 : else
1638 0 : ugb%my_nband = 0; ugb%cg_k => null()
1639 : end if
1640 :
1641 0 : call xmpi_min(ugb%my_nband, min_my_nband, comm, ierr)
1642 0 : ugb%has_idle_procs = min_my_nband == 0
1643 :
1644 0 : if (psps%usepaw == 1 .and. ugb%my_nband > 0) then
1645 : ! Calculate <Proj_i|Cnk> from output eigenstates. Note array allocated with ugb%my_nband
1646 0 : ABI_MALLOC(ugb%cprj_k, (cryst%natom, nspinor * ugb%my_nband))
1647 0 : call pawcprj_alloc(ugb%cprj_k, 0, gs_hamk%dimcprj)
1648 0 : idir = 0; cprj_choice = 1 ! Only projected wave functions.
1649 :
1650 0 : do my_ib=1,ugb%my_nband
1651 0 : ibs1 = nspinor * (my_ib - 1) + 1
1652 : call getcprj(cprj_choice, 0, ugb%cg_k(:,:,my_ib), ugb%cprj_k(:,ibs1), &
1653 : gs_hamk%ffnl_k, idir, gs_hamk%indlmn, gs_hamk%istwf_k, gs_hamk%kg_k, &
1654 : gs_hamk%kpg_k, gs_hamk%kpt_k, gs_hamk%lmnmax, gs_hamk%mgfft, mpi_enreg_seq, 1, &
1655 : gs_hamk%natom, gs_hamk%nattyp, gs_hamk%ngfft, gs_hamk%nloalg, gs_hamk%npw_k, gs_hamk%nspinor, &
1656 0 : gs_hamk%ntypat, gs_hamk%phkxred, gs_hamk%ph1d, gs_hamk%ph3d_k, gs_hamk%ucvol, gs_hamk%useylm)
1657 : end do
1658 :
1659 : ! Reorder the cprj (order is now the same as in the input file)
1660 0 : call pawcprj_reorder(ugb%cprj_k, gs_hamk%atindx1)
1661 : end if ! usepaw
1662 :
1663 0 : call cwtime_report(" block column distribution completed", cpu, wall, gflops)
1664 :
1665 : ! Free memory.
1666 :
1667 0 : ABI_FREE(eig_ene)
1668 0 : ABI_FREE(kpg_k)
1669 0 : ABI_FREE(ph3d)
1670 0 : ABI_FREE(ffnl)
1671 0 : call destroy_mpi_enreg(mpi_enreg_seq); call gs_hamk%free(); call psb%free()
1672 :
1673 0 : if (my_rank == master) call pstat_proc%print(_PSTAT_ARGS_)
1674 :
1675 0 : call timab(1919, 2, tsec)
1676 :
1677 0 : end subroutine ugb_from_diago
1678 : !!***
1679 :
1680 : !!****f* m_ksdiago/ugb_from_wfk_file
1681 : !! NAME
1682 : !! ugb_from_wfk_file
1683 : !!
1684 : !! FUNCTION
1685 : !! Initialize an ugb_t instance from a WFK file.
1686 : !!
1687 : !! INPUTS
1688 : !! spin: spin index.
1689 : !! kpoint(3)
1690 : !! comm=MPI communicator.
1691 : !!
1692 : !! OUTPUT
1693 : !! eig_k(1:nband_k)=The calculatated eigenvalues in ascending order.
1694 : !!
1695 : !! SOURCE
1696 :
1697 0 : subroutine ugb_from_wfk_file(ugb, ik_ibz, spin, istwf_k, kpoint, nband_k, &
1698 : dtset, dtfil, cryst, eig_k, comm)
1699 :
1700 : use m_wfk
1701 :
1702 : !Arguments ------------------------------------
1703 : !scalars
1704 : class(ugb_t),target,intent(out) :: ugb
1705 : integer,intent(in) :: ik_ibz, spin, istwf_k
1706 : real(dp),intent(in) :: kpoint(3)
1707 : type(dataset_type),intent(in) :: dtset
1708 : type(datafiles_type),intent(in) :: dtfil
1709 : integer,intent(in) :: nband_k, comm
1710 : type(crystal_t),intent(in) :: cryst
1711 : !arrays
1712 : real(dp),allocatable,intent(out) :: eig_k(:)
1713 :
1714 : !Local variables-------------------------------
1715 : !scalars
1716 : integer,parameter :: master = 0, formeig0 = 0
1717 : integer :: ierr, bcast_comm, color, min_my_nband
1718 : integer :: nprocs, my_rank, nbsum, npwsp, bstart, bstop, band_step, nb, npw_k, col_bsize, band, ib, il_b, iloc
1719 : logical :: have_band
1720 0 : type(ebands_t) :: wfk_ebands
1721 0 : type(wfk_t) :: wfk
1722 0 : type(hdr_type) :: wfk_hdr
1723 : character(len=fnlen) :: wfk_path
1724 : character(len=500) :: msg
1725 : !type(gs_hamiltonian_type) :: gs_hamk
1726 : !arrays
1727 : integer :: units(2)
1728 0 : real(dp),target,allocatable :: cg_work(:,:,:)
1729 0 : real(dp),contiguous, pointer :: cg_k(:,:)
1730 : ! *********************************************************************
1731 :
1732 0 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1733 0 : units(:) = [std_out, ab_out]
1734 :
1735 0 : wfk_path = dtfil%fnamewffk
1736 0 : if (my_rank == master) then
1737 0 : if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then
1738 0 : ABI_ERROR(sjoin("Cannot find HYBRYD WFK file:", wfk_path, ". Error:", msg))
1739 : end if
1740 0 : call wrtout(units, sjoin("- Reading HYBRID orbitals from WFK file:", wfk_path), pre_newlines=2)
1741 : end if
1742 :
1743 : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
1744 0 : call xmpi_bcast(wfk_path, master, comm, ierr)
1745 :
1746 : ! Read energies and performs some basic consistency checks.
1747 0 : wfk_ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
1748 0 : call wfk_hdr%vs_dtset(dtset)
1749 0 : ABI_CHECK_IEQ(dtset%ixc, wfk_hdr%ixc, "dtset%ixc /= wfk_hdr%ixc")
1750 0 : ABI_CHECK(all(abs(wfk_hdr%kptns(:,ik_ibz) - kpoint) < tol6), "Different kpoint")
1751 0 : ABI_CHECK_IRANGE(nband_k, 1, wfk_ebands%mband, "nband_k > mband.")
1752 :
1753 0 : ABI_MALLOC(eig_k, (nband_k))
1754 0 : eig_k = wfk_ebands%eig(1:nband_k, ik_ibz, spin)
1755 :
1756 0 : npw_k = wfk_hdr%npwarr(ik_ibz)
1757 0 : npwsp = npw_k * wfk_hdr%nspinor
1758 0 : ABI_CHECK_IEQ(istwf_k, wfk_hdr%istwfk(ik_ibz), "different istwfk_k")
1759 :
1760 : ! Init scalapack matrix
1761 0 : call ugb%processor%init(comm, grid_dims=[1, nprocs])
1762 0 : ABI_CHECK(block_dist_1d(nband_k, nprocs, col_bsize, msg), msg)
1763 0 : call ugb%mat%init(npwsp, nband_k, ugb%processor, 1, size_blocs=[-1, col_bsize])
1764 :
1765 0 : ABI_MALLOC(ugb%kg_k, (3, npw_k))
1766 :
1767 : ! Master reads and broadcasts. Much faster on lumi
1768 0 : if (my_rank == master) then
1769 0 : call wfk%open_read(wfk_path, formeig0, iomode_from_fname(wfk_path), get_unit(), xmpi_comm_self)
1770 : end if
1771 :
1772 : ! TODO: Optimize this part
1773 : ! Find band_step that gives good compromise between memory and efficiency.
1774 : !band_step = memb_limited_step(1, nbsum, 2*npwsp, xmpi_bsize_dp, 1024.0_dp)
1775 0 : band_step = 200
1776 : !band_step = 100
1777 0 : nbsum = nband_k
1778 0 : do bstart=1, nbsum, band_step
1779 0 : bstop = min(bstart + band_step - 1, nbsum); nb = bstop - bstart + 1
1780 :
1781 0 : ABI_MALLOC(cg_work, (2, npwsp, nb)) ! This array is always dp
1782 0 : if (my_rank == master) then
1783 0 : call c_f_pointer(c_loc(cg_work), cg_k, shape=[2, npwsp * nb])
1784 0 : call wfk%read_band_block([bstart, bstop], ik_ibz, spin, xmpio_single, kg_k=ugb%kg_k, cg_k=cg_k)
1785 : end if
1786 :
1787 0 : call xmpi_bcast(ugb%kg_k, master, comm, ierr)
1788 :
1789 : ! Create communicator with master and all procs requiring this set of bands block (color == 1)
1790 0 : color = 0
1791 0 : do band=bstart, bstop
1792 0 : call ugb%mat%glob2loc(1, band, iloc, il_b, have_band)
1793 0 : if (have_band) then
1794 0 : color = 1; exit
1795 : end if
1796 : end do
1797 0 : if (my_rank == master) color = 1
1798 0 : call xmpi_comm_split(comm, color, my_rank, bcast_comm, ierr)
1799 :
1800 0 : if (color == 1) then
1801 0 : call xmpi_bcast(cg_work, master, bcast_comm, ierr)
1802 : endif
1803 0 : call xmpi_comm_free(bcast_comm)
1804 :
1805 : ! Copy my portion of cg_work to buffer_cplx (here we have dp --> sp conversion).
1806 0 : if (color == 1) then
1807 0 : do band=bstart, bstop
1808 0 : ib = band - bstart + 1
1809 0 : call ugb%mat%glob2loc(1, band, iloc, il_b, have_band); if (.not. have_band) cycle
1810 0 : ugb%mat%buffer_cplx(:, il_b) = cmplx(cg_work(1,:,ib), cg_work(2,:,ib), kind=gwp)
1811 : end do
1812 : end if
1813 0 : ABI_FREE(cg_work)
1814 : end do ! bstart
1815 :
1816 0 : if (my_rank == master) call wfk%close()
1817 :
1818 : ! =================
1819 : ! Build ugb object
1820 : ! =================
1821 0 : ugb%istwf_k = istwf_k
1822 0 : ugb%nspinor = wfk_hdr%nspinor
1823 0 : ugb%npw_k = npw_k
1824 0 : ugb%npwsp = npwsp
1825 0 : ugb%nband_k = nband_k
1826 0 : ugb%comm => ugb%mat%processor%comm
1827 :
1828 0 : ugb%my_bstart = ugb%mat%loc2gcol(1)
1829 0 : ugb%my_bstop = ugb%mat%loc2gcol(ugb%mat%size_local(2))
1830 0 : ugb%my_nband = ugb%my_bstop - ugb%my_bstart + 1
1831 :
1832 0 : if (ugb%my_nband > 0) then
1833 0 : call c_f_pointer(c_loc(ugb%mat%buffer_cplx), ugb%cg_k, shape=[2, ugb%npwsp, ugb%my_nband])
1834 : else
1835 0 : ugb%my_nband = 0
1836 0 : ugb%cg_k => null()
1837 : end if
1838 :
1839 0 : call xmpi_min(ugb%my_nband, min_my_nband, comm, ierr)
1840 0 : ugb%has_idle_procs = min_my_nband == 0
1841 :
1842 : ! TODO
1843 0 : if (dtset%usepaw == 1 .and. ugb%my_nband > 0) then
1844 0 : ABI_ERROR("ugb_from_wfk does not support PAW")
1845 : ! Calculate <Proj_i|Cnk> from output eigenstates. Note array allocated with ugb%my_nband
1846 0 : ABI_MALLOC(ugb%cprj_k, (cryst%natom, ugb%nspinor * ugb%my_nband))
1847 : !call pawcprj_alloc(ugb%cprj_k, 0, gs_hamk%dimcprj)
1848 : !idir = 0; cprj_choice = 1 ! Only projected wave functions.
1849 :
1850 : !do my_ib=1,ugb%my_nband
1851 : ! ibs1 = nspinor * (my_ib - 1) + 1
1852 : ! call getcprj(cprj_choice, 0, ugb%cg_k(:,:,my_ib), ugb%cprj_k(:,ibs1), &
1853 : ! gs_hamk%ffnl_k, idir, gs_hamk%indlmn, gs_hamk%istwf_k, gs_hamk%kg_k, &
1854 : ! gs_hamk%kpg_k, gs_hamk%kpt_k, gs_hamk%lmnmax, gs_hamk%mgfft, mpi_enreg_seq, &
1855 : ! gs_hamk%natom, gs_hamk%nattyp, gs_hamk%ngfft, gs_hamk%nloalg, gs_hamk%npw_k, gs_hamk%nspinor, &
1856 : ! gs_hamk%ntypat, gs_hamk%phkxred, gs_hamk%ph1d, gs_hamk%ph3d_k, gs_hamk%ucvol, gs_hamk%useylm)
1857 : !end do
1858 :
1859 : !! Reorder the cprj (order is now the same as in the input file)
1860 : !call pawcprj_reorder(ugb%cprj_k, gs_hamk%atindx1)
1861 : end if ! usepaw
1862 :
1863 0 : call ugb%print(units, dtset%prtvol)
1864 :
1865 0 : call wfk_hdr%free(); call wfk_ebands%free()
1866 :
1867 0 : end subroutine ugb_from_wfk_file
1868 : !!***
1869 :
1870 : !----------------------------------------------------------------------
1871 :
1872 : !!****f* m_gwr/ugb_free
1873 : !! NAME
1874 : !! ugb_free
1875 : !!
1876 : !! FUNCTION
1877 : !! Free dynamic memory.
1878 : !!
1879 : !! SOURCE
1880 :
1881 0 : subroutine ugb_free(ugb)
1882 :
1883 : !Arguments ------------------------------------
1884 : class(ugb_t),intent(inout) :: ugb
1885 : ! *************************************************************************
1886 :
1887 0 : call ugb%mat%free()
1888 0 : call ugb%processor%free()
1889 0 : ABI_SFREE(ugb%kg_k)
1890 0 : ugb%cg_k => null()
1891 0 : ugb%comm => null()
1892 :
1893 0 : if (allocated(ugb%cprj_k)) then
1894 0 : call pawcprj_free(ugb%cprj_k)
1895 0 : ABI_FREE(ugb%cprj_k)
1896 : end if
1897 :
1898 0 : end subroutine ugb_free
1899 : !!***
1900 :
1901 : !----------------------------------------------------------------------
1902 :
1903 : !!****f* m_gwr/ugb_print
1904 : !! NAME
1905 : !! ugb_print
1906 : !!
1907 : !! FUNCTION
1908 : !! Print info on the object.
1909 : !!
1910 : !! SOURCE
1911 :
1912 0 : subroutine ugb_print(ugb, units, prtvol, header)
1913 :
1914 : !Arguments ------------------------------------
1915 : class(ugb_t),intent(in) :: ugb
1916 : integer,intent(in) :: units(:), prtvol
1917 : character(len=*),optional,intent(in) :: header
1918 :
1919 : !Local variables-------------------------------
1920 : character(len=500) :: msg
1921 : type(yamldoc_t) :: ydoc
1922 : ! *************************************************************************
1923 :
1924 : ABI_UNUSED(prtvol)
1925 :
1926 0 : msg = ' ==== Info on the ugb_t object ==== '; if (present(header)) msg = ' ==== '//trim(adjustl(header))//' ==== '
1927 0 : call wrtout(units, msg)
1928 :
1929 0 : ydoc = yamldoc_open('ugb_t') !, width=11, real_fmt='(3f8.3)')
1930 0 : call ydoc%add_int("istwf_k", ugb%istwf_k)
1931 0 : call ydoc%add_int("nspinor", ugb%nspinor)
1932 0 : call ydoc%add_int("npw_k", ugb%npw_k)
1933 0 : call ydoc%add_int("nband_k", ugb%nband_k)
1934 0 : call ydoc%add_int("my_bstart", ugb%my_bstart)
1935 0 : call ydoc%add_int("my_bstop", ugb%my_bstop)
1936 0 : call ydoc%add_int("my_nband", ugb%my_nband)
1937 0 : call ydoc%write_units_and_free(units)
1938 :
1939 0 : end subroutine ugb_print
1940 : !!***
1941 : !----------------------------------------------------------------------
1942 :
1943 : !!****f* m_ksdiago/ugb_collect_cprj
1944 : !! NAME
1945 : !! ugb_collect_cprj
1946 : !!
1947 : !! FUNCTION
1948 : !! This is a collective routine that returns in `out_cprj` the PAW projections
1949 : !! for `nb` bands starting at `band_start` NB: `out_cprj` is supposed to be allocated in the parent
1950 : !!
1951 : !! SOURCE
1952 :
1953 0 : subroutine ugb_collect_cprj(ugb, nspinor, nb, band_start, out_cprj)
1954 :
1955 : !Arguments ------------------------------------
1956 : class(ugb_t),intent(in) :: ugb
1957 : integer,intent(in) :: nspinor, nb, band_start
1958 : type(pawcprj_type),intent(inout) :: out_cprj(:,:)
1959 :
1960 : !Local variables-------------------------------
1961 : integer :: ierr, my_ibs, out_ibs, band, cnt
1962 : ! *************************************************************************
1963 :
1964 0 : ABI_CHECK_IEQ(size(ugb%cprj_k, dim=1), size(out_cprj, dim=1), "size1 should be the same")
1965 0 : ABI_CHECK_IGEQ(size(out_cprj, dim=2), nb*nspinor, "size2 too small!")
1966 :
1967 : ! TODO: Numb algorithm based on xmpi_sum. Might be optimized.
1968 0 : call pawcprj_set_zero(out_cprj)
1969 :
1970 0 : cnt = nspinor - 1
1971 0 : do band=band_start, band_start+nb-1
1972 0 : if (band >= ugb%my_bstart .and. band <= ugb%my_bstop) then
1973 0 : my_ibs = 1 + (band - ugb%my_bstart) * nspinor
1974 0 : out_ibs = 1 + (band - band_start) * nspinor
1975 0 : call pawcprj_copy(ugb%cprj_k(:,my_ibs:my_ibs+cnt), out_cprj(:,out_ibs:out_ibs+cnt))
1976 : end if
1977 : end do
1978 :
1979 0 : call pawcprj_mpi_sum(out_cprj, ugb%comm, ierr)
1980 :
1981 0 : end subroutine ugb_collect_cprj
1982 : !!***
1983 :
1984 : !----------------------------------------------------------------------
1985 :
1986 : !!****f* m_gwr/hyb_from_wfk_file
1987 : !! NAME
1988 : !! hyb_from_wfk_file
1989 : !!
1990 : !! FUNCTION
1991 : !! Read the WFK file compute with HYBRID functional
1992 : !!
1993 : !! SOURCE
1994 :
1995 0 : subroutine hyb_from_wfk_file(hyb, cryst, dtfil, dtset, psps, pawtab, ngfftc, diago_pool, comm)
1996 :
1997 : use m_krank
1998 : use m_kpts
1999 :
2000 : !Arguments ------------------------------------
2001 : class(hyb_t),intent(out) :: hyb
2002 : type(crystal_t),intent(in) :: cryst
2003 : type(datafiles_type),intent(in) :: dtfil
2004 : type(dataset_type),intent(in) :: dtset
2005 : type(pseudopotential_type),intent(in) :: psps
2006 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
2007 : type(xmpi_pool2d_t),intent(in) :: diago_pool
2008 : integer,intent(in) :: ngfftc(18), comm
2009 :
2010 : !Local variables ------------------------------
2011 : integer,parameter :: master = 0
2012 : integer :: nprocs, my_rank, ierr, mband, nkibz, nsppol, spin, ik_ibz, ebands_kptopt ! b1, b2,
2013 : real(dp) :: vc_ecut
2014 : character(len=5000) :: msg
2015 0 : type(hdr_type) :: wfk_hdr
2016 0 : type(crystal_t) :: wfk_cryst
2017 0 : type(krank_t) :: krank_ibz ! qrank,
2018 : character(len=fnlen) :: wfk_path
2019 : integer :: nqbzX
2020 : integer :: units(2)
2021 0 : integer,allocatable :: nband(:,:), wfd_istwfk(:), qtab(:), qtabi(:), qtabo(:)
2022 0 : real(dp),allocatable :: qbz(:,:), wtk(:), wtq(:)
2023 0 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
2024 : !************************************************************************
2025 :
2026 0 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2027 0 : units(:) = [std_out, ab_out]
2028 :
2029 0 : wfk_path = dtfil%fnamewffk
2030 0 : if (my_rank == master) then
2031 0 : if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then
2032 0 : ABI_ERROR(sjoin("Cannot find HYBRYD WFK file:", wfk_path, ". Error:", msg))
2033 : end if
2034 0 : call wrtout(units, sjoin("- Reading HYBRID orbitals from WFK file:", wfk_path), pre_newlines=2)
2035 : end if
2036 :
2037 : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
2038 0 : call xmpi_bcast(wfk_path, master, comm, ierr)
2039 :
2040 : ! Construct crystal and hyb%ebands from the GS WFK file.
2041 0 : hyb%ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
2042 0 : call wfk_hdr%vs_dtset(dtset)
2043 0 : ABI_CHECK_IEQ(dtset%ixc, wfk_hdr%ixc, "dtset%ixc /= wfk_hdr%ixc")
2044 :
2045 0 : wfk_cryst = wfk_hdr%get_crystal()
2046 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
2047 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
2048 : end if
2049 : !call wfk_cryst%print(header="crystal structure from WFK file")
2050 0 : call wfk_cryst%free()
2051 : ! TODO: Add more consistency checks e.g. nkibz,...
2052 : !cryst = wfk_hdr%get_crystal()
2053 : !call cryst%print(header="crystal structure from WFK file")
2054 :
2055 0 : nkibz = hyb%ebands%nkpt; nsppol = hyb%ebands%nsppol
2056 0 : mband = hyb%ebands%mband
2057 :
2058 : ! Initialize the wave function descriptor.
2059 : ! Only wavefunctions for the symmetrical imagine of the k wavevectors
2060 : ! treated by this MPI rank are stored.
2061 0 : ABI_MALLOC(nband, (nkibz, nsppol))
2062 0 : ABI_MALLOC(bks_mask, (mband, nkibz, nsppol))
2063 0 : ABI_MALLOC(keep_ur, (mband, nkibz, nsppol))
2064 0 : nband = mband; bks_mask = .False.; keep_ur = .False.
2065 :
2066 : ! Set tolerance used to decide if a band is empty
2067 : !tol_empty_in = 0.01_dp
2068 : !call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
2069 :
2070 : !do hyb_ik_ibz=1,nkibz
2071 0 : do spin=1,nsppol
2072 0 : if (all(.not. diago_pool%treats(:, spin))) cycle ! MPI distribution of collinear spins.
2073 0 : do ik_ibz=1,nkibz
2074 0 : bks_mask(:, ik_ibz, spin) = .True.
2075 : !bks_mask(b1:b2, ik_ibz, spin) = .True.
2076 : end do
2077 : end do
2078 : !end do
2079 :
2080 : ! Impose istwfk = 1 for all k-points.
2081 : ! wfd_read_wfk will handle a possible conversion if the WFK contains istwfk /= 1.
2082 0 : ABI_MALLOC(wfd_istwfk, (nkibz))
2083 0 : wfd_istwfk = 1 !; wfd_istwfk = wfk_hdr%istwf_k
2084 :
2085 : call hyb%wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkibz, dtset%nsppol, bks_mask, &
2086 : dtset%nspden, dtset%nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, hyb%ebands%kptns, ngfftc, &
2087 0 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
2088 :
2089 0 : call hyb%wfd%print([std_out], header="Wavefunctions for Hybrid WKF file")
2090 :
2091 0 : ABI_FREE(nband)
2092 0 : ABI_FREE(keep_ur)
2093 0 : ABI_FREE(wfd_istwfk)
2094 0 : ABI_FREE(bks_mask)
2095 :
2096 0 : call wfk_hdr%free()
2097 :
2098 : ! Read wavefunctions.
2099 0 : call hyb%wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path))
2100 :
2101 : ! This piece of code is taken from m_gwr.
2102 :
2103 : ! =======================
2104 : ! Setup k-mesh and q-mesh
2105 : ! =======================
2106 : ! Get full kBZ associated to hyb%ebands
2107 : call kpts_ibz_from_kptrlatt(cryst, hyb%ebands%kptrlatt, hyb%ebands%kptopt, hyb%ebands%nshiftk, hyb%ebands%shiftk, &
2108 0 : hyb%nkibz, hyb%kibz, wtk, hyb%nkbz, hyb%kbz) !, bz2ibz=bz2ibz)
2109 : !new_kptrlatt=gwr%kptrlatt, new_shiftk=gwr%kshift,
2110 : !bz2ibz=new%ind_qbz2ibz) # FIXME
2111 0 : ABI_FREE(wtk)
2112 :
2113 : ! In principle kibz should be equal to hyb%ebands%kptns.
2114 0 : ABI_CHECK_IEQ(hyb%nkibz, hyb%ebands%nkpt, "nkibz != hyb%ebands%nkpt")
2115 0 : ABI_CHECK(all(abs(hyb%ebands%kptns - hyb%kibz) < tol12), "hyb%ebands%kibz != hyb%kibz")
2116 :
2117 : ! Note symrec convention.
2118 0 : ebands_kptopt = hyb%ebands%kptopt
2119 0 : call krank_ibz%from_kptrlatt(hyb%nkibz, hyb%kibz, hyb%ebands%kptrlatt, compute_invrank=.False.)
2120 :
2121 0 : ABI_MALLOC(hyb%kbz2ibz, (6, hyb%nkbz))
2122 0 : if (kpts_map("symrec", ebands_kptopt, cryst, krank_ibz, hyb%nkbz, hyb%kbz, hyb%kbz2ibz) /= 0) then
2123 0 : ABI_ERROR("Cannot map kBZ to IBZ!")
2124 : end if
2125 :
2126 : ! Order kbz by stars and rearrange entries in kbz2ibz table.
2127 0 : call kpts_pack_in_stars(hyb%nkbz, hyb%kbz, hyb%kbz2ibz)
2128 :
2129 0 : if (my_rank == master) then
2130 0 : call kpts_map_print(units, " Mapping kBZ --> kIBZ", "symrec", hyb%kbz, hyb%kibz, hyb%kbz2ibz, dtset%prtvol)
2131 : end if
2132 :
2133 : ! Table with symrel conventions for the symmetrization of the wfs.
2134 0 : ABI_MALLOC(hyb%kbz2ibz_symrel, (6, hyb%nkbz))
2135 0 : if (kpts_map("symrel", ebands_kptopt, cryst, krank_ibz, hyb%nkbz, hyb%kbz, hyb%kbz2ibz_symrel) /= 0) then
2136 0 : ABI_ERROR("Cannot map kBZ to IBZ!")
2137 : end if
2138 0 : call krank_ibz%free()
2139 :
2140 : ! Setup qIBZ, weights and BZ.
2141 : ! Always use q --> -q symmetry even in systems without inversion
2142 : ! TODO: Might add input variable to rescale the q-mesh.
2143 :
2144 : ! Find the number of q-points such that q = k1-k2.
2145 0 : call findnq(hyb%nkbz, hyb%kbz, cryst%nsym, cryst%symrec, cryst%symafm, hyb%nqibz, cryst%timrev)
2146 :
2147 : ! Find the coordinates of the q-points in the IBZ.
2148 0 : ABI_MALLOC(hyb%qibz, (3, hyb%nqibz))
2149 0 : call findq(hyb%nkbz, hyb%kbz, cryst%nsym, cryst%symrec, cryst%symafm, cryst%gprimd, hyb%nqibz, hyb%qibz, cryst%timrev)
2150 0 : ABI_CHECK(all(abs(hyb%qibz(:,1)) < tol16), "First qpoint in qibz should be Gamma!")
2151 :
2152 : ! HM: the bz2ibz produced above is incomplete, I do it here using listkk
2153 : !ABI_MALLOC(hyb%qbz2ibz, (6, hyb%nqbz))
2154 : !qrank = krank_from_kptrlatt(hyb%nqibz, hyb%qibz, qptrlatt, compute_invrank=.False.)
2155 :
2156 : !if (kpts_map("symrec", qtimrev1, cryst, qrank, hyb%nqbz, hyb%qbz, hyb%qbz2ibz) /= 0) then
2157 : ! ABI_ERROR("Cannot map qBZ to IBZ!")
2158 : !end if
2159 : !call qrank%free()
2160 :
2161 : ! Order qbz by stars and rearrange entries in qbz2ibz table.
2162 : !call kpts_pack_in_stars(hyb%nqbz, hyb%qbz, hyb%qbz2ibz)
2163 : !if (my_rank == master) then
2164 : ! call kpts_map_print(units, " Mapping qBZ --> qIBZ", "symrec", hyb%qbz, hyb%qibz, hyb%qbz2ibz, dtset%prtvol)
2165 : !end if
2166 :
2167 0 : nqbzX = hyb%nqibz*cryst%nsym*cryst%timrev ! Maximum possible number
2168 0 : ABI_MALLOC(qbz, (3, nqbzX))
2169 0 : ABI_MALLOC(wtq, (hyb%nqibz))
2170 0 : ABI_MALLOC(qtab, (nqbzX))
2171 0 : ABI_MALLOC(qtabi, (nqbzX))
2172 0 : ABI_MALLOC(qtabo, (nqbzX))
2173 :
2174 0 : call identk(hyb%qibz, hyb%nqibz, nqbzX, cryst%nsym, cryst%timrev, cryst%symrec, cryst%symafm, qbz, qtab, qtabi, qtabo, hyb%nqbz, wtq)
2175 :
2176 0 : ABI_MALLOC(hyb%qbz, (3, hyb%nqibz))
2177 0 : hyb%qbz = qbz(:,1:hyb%nqibz)
2178 :
2179 0 : ABI_FREE(qbz)
2180 0 : ABI_FREE(wtq)
2181 0 : ABI_FREE(qtab)
2182 0 : ABI_FREE(qtabi)
2183 0 : ABI_FREE(qtabo)
2184 :
2185 : ! TODO: MC technique does not seem to work as expected, even in the legacy code.
2186 0 : vc_ecut = dtset%ecut ! * four
2187 : call hyb%vcgen%init(cryst, hyb%ebands%kptrlatt, hyb%nkbz, hyb%nqibz, hyb%nqbz, hyb%qbz, &
2188 0 : dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, vc_ecut, comm)
2189 :
2190 0 : end subroutine hyb_from_wfk_file
2191 : !!***
2192 :
2193 : !----------------------------------------------------------------------
2194 :
2195 : !!****f* m_gwr/hyb_free
2196 : !! NAME
2197 : !! hyb_free
2198 : !!
2199 : !! FUNCTION
2200 : !! Free dynamic memory.
2201 : !!
2202 : !! SOURCE
2203 :
2204 0 : subroutine hyb_free(hyb)
2205 :
2206 : !Arguments ------------------------------------
2207 : class(hyb_t),intent(inout) :: hyb
2208 : ! *************************************************************************
2209 :
2210 0 : ABI_SFREE(hyb%kibz)
2211 0 : ABI_SFREE(hyb%kbz)
2212 0 : ABI_SFREE(hyb%qibz)
2213 0 : ABI_SFREE(hyb%qbz)
2214 0 : ABI_SFREE(hyb%wtq)
2215 0 : ABI_SFREE(hyb%kbz2ibz)
2216 0 : ABI_SFREE(hyb%kbz2ibz_symrel)
2217 0 : ABI_SFREE(hyb%qbz2ibz)
2218 :
2219 : ! Free datatypes
2220 0 : call hyb%wfd%free(); call hyb%vcgen%free(); call hyb%ebands%free()
2221 :
2222 0 : end subroutine hyb_free
2223 : !!***
2224 :
2225 : !!****f* m_ksdiago/psbands_init
2226 : !! NAME
2227 : !! psbands_init
2228 : !!
2229 : !! FUNCTION
2230 : !!
2231 : !! INPUTS
2232 : !!
2233 : !! OUTPUT
2234 : !!
2235 : !! SOURCE
2236 :
2237 0 : subroutine psbands_init(psb, dtset, eig_size, eig_k, gs_fermie)
2238 :
2239 : !Arguments ------------------------------------
2240 : class(psbands_t),intent(out) :: psb
2241 : class(dataset_type),target,intent(in) :: dtset
2242 : integer,intent(in) :: eig_size
2243 : real(dp),intent(in) :: gs_fermie
2244 : !arrays
2245 : real(dp),intent(in) :: eig_k(eig_size)
2246 :
2247 : !Local variables-------------------------------
2248 : !scalars
2249 : integer :: islice, ib, cnt, units(2), first_band, last_band, nb
2250 : real(dp) :: first_eig, last_eig
2251 0 : real(dp),allocatable :: tmp_eig_k(:)
2252 : ! *********************************************************************
2253 :
2254 : ABI_UNUSED(dtset%natom)
2255 :
2256 : ! Shift energies wrt the input Fermi level.
2257 0 : ABI_MALLOC(tmp_eig_k, (eig_size))
2258 0 : tmp_eig_k = eig_k - gs_fermie
2259 :
2260 0 : psb%nb_protected = huge(1)
2261 : !psb%nb_protected = dtset%nb_protected
2262 : !psb%maxsto_per_slice = dtset%nb_per_slice
2263 : ! TODO
2264 0 : psb%efrac = 0.02_dp ! dtset%efrac
2265 :
2266 : ! Compute nslices and subspace
2267 : ! TODO: Add possibility of treating occupied states as well?
2268 0 : ABI_MALLOC(psb%subspace, (3, eig_size))
2269 0 : first_band = psb%nb_protected + 1
2270 0 : psb%nslices = 0
2271 :
2272 : do while (first_band > 0)
2273 : first_eig = tmp_eig_k(first_band)
2274 : last_eig = first_eig + (first_eig * psb%efrac)
2275 : last_band = get_band_with_energy_small_than(first_band+1, eig_size, last_eig)
2276 : psb%nslices = psb%nslices + 1
2277 : psb%subspace(1, psb%nslices) = first_band
2278 : if (last_band == -1) then
2279 : psb%subspace(2, psb%nslices) = eig_size
2280 : else
2281 : psb%subspace(2, psb%nslices) = last_band
2282 : end if
2283 : nb = psb%subspace(2, psb%nslices) - psb%subspace(1, psb%nslices) + 1
2284 : if (last_band == first_band) then
2285 : ! Won't use pseudo bands in this case.
2286 : psb%subspace(3, psb%nslices) = 1
2287 : else
2288 : !psb%subspace(3, psb%nslices) = min(dtset%nb_per_slice, nb)
2289 : end if
2290 : first_band = last_band + 1
2291 : !write(std_out,'(a,i0,a,*(1x,i0))')" islice: ", psb%nslices, " subspace:", psb%subspace(:, psb%nslices)
2292 : end do
2293 :
2294 : ! Copy eigenvalues of the protected states.
2295 0 : psb%nb_tot = psb%nb_protected + sum(psb%subspace(3,1:psb%nslices))
2296 0 : ABI_MALLOC(psb%ps_eig, (psb%nb_tot))
2297 0 : psb%ps_eig(1:psb%nb_protected) = tmp_eig_k(1:psb%nb_protected)
2298 :
2299 : cnt = 0
2300 : do islice=1,psb%nslices
2301 : first_band = psb%subspace(1, islice)
2302 : last_band = psb%subspace(2, islice)
2303 : ! Take average of eigenvalues inside the slice.
2304 : do ib=1,psb%subspace(3, islice)
2305 : cnt = cnt + 1
2306 : psb%ps_eig(psb%nb_protected + cnt) = sum(tmp_eig_k(first_band:last_band)) / dble(last_band - first_band + 1)
2307 : end do
2308 : end do
2309 0 : ABI_FREE(tmp_eig_k)
2310 :
2311 0 : psb%ps_eig = psb%ps_eig + gs_fermie
2312 :
2313 0 : units = [std_out, ab_out]
2314 0 : call wrtout(units, ' Stochastic pseudobands setup:', pre_newlines=1)
2315 0 : call wrtout(units, sjoin(' Number of stochastic subspaces: ', itoa(psb%nslices)))
2316 : !call wrtout(units, sjoin(' Number of stochastic pseudobands per subspace: ', itoa(dtset%nb_per_slice)))
2317 0 : call wrtout(units, sjoin(' Original number of bands: ', itoa(eig_size)))
2318 0 : call wrtout(units, sjoin(' Number of bands in the protection window: ', itoa(psb%nb_protected)))
2319 0 : call wrtout(units, sjoin(' Final number of bands: ', itoa(psb%nb_tot)), newlines=1)
2320 :
2321 : !if (dtset%prtvol > 5) then
2322 : ! do islice=1,psb%nslices
2323 : ! write(msg,'(a,i0,a,*(1x,i0))')" islice: ", psb%nslices, " subspace:", psb%subspace(:, psb%nslices)
2324 : ! call wrtout(units, msg)
2325 : ! end do
2326 : !end if
2327 :
2328 : contains
2329 :
2330 : integer function get_band_with_energy_small_than(idx_start, idx_end, energy) result(band)
2331 : integer, intent(in) :: idx_start, idx_end
2332 : integer :: ib
2333 : real(dp), intent(in) :: energy
2334 :
2335 : band = -1
2336 : do ib=idx_start,idx_end
2337 : if (tmp_eig_k(ib) > energy) then
2338 : band = ib - 1; return
2339 : end if
2340 : end do
2341 : end function get_band_with_energy_small_than
2342 :
2343 : end subroutine psbands_init
2344 : !!***
2345 :
2346 : !!****f* m_ksdiago/psbands_band2slice
2347 : !! NAME
2348 : !! psbands_band2slice
2349 : !!
2350 : !! FUNCTION
2351 : !! Return the slice index from the band index. -1 if band is protected.
2352 : !!
2353 : !! SOURCE
2354 :
2355 0 : integer function psbands_band2slice(psb, band) result(islice)
2356 :
2357 : !Arguments ------------------------------------
2358 : class(psbands_t),intent(in) :: psb
2359 : integer,intent(in) :: band
2360 : ! *********************************************************************
2361 :
2362 0 : do islice=1,psb%nslices
2363 0 : if (band >= psb%subspace(1,islice) .and. &
2364 0 : band <= psb%subspace(2,islice)) return
2365 : end do
2366 0 : islice = -1
2367 :
2368 : end function psbands_band2slice
2369 : !!***
2370 :
2371 : !!****f* m_ksdiago/psbands_free
2372 : !! NAME
2373 : !! psbands_free
2374 : !!
2375 : !! FUNCTION
2376 : !! Free memory
2377 : !!
2378 : !! SOURCE
2379 :
2380 0 : subroutine psbands_free(psb)
2381 :
2382 : !Arguments ------------------------------------
2383 : class(psbands_t),intent(inout) :: psb
2384 : ! *********************************************************************
2385 :
2386 0 : ABI_SFREE(psb%ps_eig)
2387 0 : ABI_SFREE(psb%subspace)
2388 :
2389 0 : end subroutine psbands_free
2390 : !!***
2391 :
2392 0 : end module m_ksdiago
2393 : !!***
|
