Branch data 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 : : character(len=5000) :: msg
2014 : 0 : type(hdr_type) :: wfk_hdr
2015 [ # # # # : 0 : type(crystal_t) :: wfk_cryst
# # # # #
# # # # #
# # # # ]
2016 : 0 : type(krank_t) :: krank_ibz ! qrank,
2017 : : character(len=fnlen) :: wfk_path
2018 : : integer :: nqbzX
2019 : : integer :: units(2)
2020 : 0 : integer,allocatable :: nband(:,:), wfd_istwfk(:), qtab(:), qtabi(:), qtabo(:)
2021 [ # # ]: 0 : real(dp),allocatable :: qbz(:,:), wtk(:), wtq(:)
2022 : 0 : logical,allocatable :: bks_mask(:,:,:), keep_ur(:,:,:)
2023 : : !************************************************************************
2024 : :
2025 : 0 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2026 [ # # ]: 0 : units(:) = [std_out, ab_out]
2027 : :
2028 : 0 : wfk_path = dtfil%fnamewffk
2029 [ # # ]: 0 : if (my_rank == master) then
2030 [ # # ]: 0 : if (nctk_try_fort_or_ncfile(wfk_path, msg) /= 0) then
2031 : 0 : ABI_ERROR(sjoin("Cannot find HYBRYD WFK file:", wfk_path, ". Error:", msg))
2032 : : end if
2033 : 0 : call wrtout(units, sjoin("- Reading HYBRID orbitals from WFK file:", wfk_path), pre_newlines=2)
2034 : : end if
2035 : :
2036 : : ! Broadcast filenames (needed because they might have been changed if we are using netcdf files)
2037 : 0 : call xmpi_bcast(wfk_path, master, comm, ierr)
2038 : :
2039 : : ! Construct crystal and hyb%ebands from the GS WFK file.
2040 [ # # # # : 0 : hyb%ebands = wfk_read_ebands(wfk_path, comm, out_hdr=wfk_hdr)
# # # # #
# # # # #
# # # # #
# # # ]
2041 : 0 : call wfk_hdr%vs_dtset(dtset)
2042 [ # # ]: 0 : ABI_CHECK_IEQ(dtset%ixc, wfk_hdr%ixc, "dtset%ixc /= wfk_hdr%ixc")
2043 : :
2044 : 0 : wfk_cryst = wfk_hdr%get_crystal()
2045 [ # # ]: 0 : if (cryst%compare(wfk_cryst, header=" Comparing input crystal with WFK crystal") /= 0) then
2046 : 0 : ABI_ERROR("Crystal structure from input and from WFK file do not agree! Check messages above!")
2047 : : end if
2048 : : !call wfk_cryst%print(header="crystal structure from WFK file")
2049 : 0 : call wfk_cryst%free()
2050 : : ! TODO: Add more consistency checks e.g. nkibz,...
2051 : : !cryst = wfk_hdr%get_crystal()
2052 : : !call cryst%print(header="crystal structure from WFK file")
2053 : :
2054 : 0 : nkibz = hyb%ebands%nkpt; nsppol = hyb%ebands%nsppol
2055 : 0 : mband = hyb%ebands%mband
2056 : :
2057 : : ! Initialize the wave function descriptor.
2058 : : ! Only wavefunctions for the symmetrical imagine of the k wavevectors
2059 : : ! treated by this MPI rank are stored.
2060 [ # # # # : 0 : ABI_MALLOC(nband, (nkibz, nsppol))
# # # # #
# # # # #
# # ]
2061 [ # # # # : 0 : ABI_MALLOC(bks_mask, (mband, nkibz, nsppol))
# # # # #
# # # # #
# # # # #
# ]
2062 [ # # # # : 0 : ABI_MALLOC(keep_ur, (mband, nkibz, nsppol))
# # # # #
# # # # #
# # ]
2063 [ # # # # : 0 : nband = mband; bks_mask = .False.; keep_ur = .False.
# # # # #
# # # # #
# # ]
2064 : :
2065 : : ! Set tolerance used to decide if a band is empty
2066 : : !tol_empty_in = 0.01_dp
2067 : : !call get_fact_spin_tol_empty(nsppol, nspinor, tol_empty_in, fact_spin, tol_empty)
2068 : :
2069 : : !do hyb_ik_ibz=1,nkibz
2070 [ # # ]: 0 : do spin=1,nsppol
2071 [ # # # # : 0 : if (all(.not. diago_pool%treats(:, spin))) cycle ! MPI distribution of collinear spins.
# # ]
2072 [ # # ]: 0 : do ik_ibz=1,nkibz
2073 [ # # ]: 0 : bks_mask(:, ik_ibz, spin) = .True.
2074 : : !bks_mask(b1:b2, ik_ibz, spin) = .True.
2075 : : end do
2076 : : end do
2077 : : !end do
2078 : :
2079 : : ! Impose istwfk = 1 for all k-points.
2080 : : ! wfd_read_wfk will handle a possible conversion if the WFK contains istwfk /= 1.
2081 [ # # # # : 0 : ABI_MALLOC(wfd_istwfk, (nkibz))
# # # # #
# # # ]
2082 [ # # ]: 0 : wfd_istwfk = 1 !; wfd_istwfk = wfk_hdr%istwf_k
2083 : :
2084 : : call hyb%wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkibz, dtset%nsppol, bks_mask, &
2085 : : dtset%nspden, dtset%nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, hyb%ebands%kptns, ngfftc, &
2086 : 0 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
2087 : :
2088 [ # # ]: 0 : call hyb%wfd%print([std_out], header="Wavefunctions for Hybrid WKF file")
2089 : :
2090 [ # # ]: 0 : ABI_FREE(nband)
2091 [ # # ]: 0 : ABI_FREE(keep_ur)
2092 [ # # ]: 0 : ABI_FREE(wfd_istwfk)
2093 [ # # ]: 0 : ABI_FREE(bks_mask)
2094 : :
2095 : 0 : call wfk_hdr%free()
2096 : :
2097 : : ! Read wavefunctions.
2098 : 0 : call hyb%wfd%read_wfk(wfk_path, iomode_from_fname(wfk_path))
2099 : :
2100 : : ! This piece of code is taken from m_gwr.
2101 : :
2102 : : ! =======================
2103 : : ! Setup k-mesh and q-mesh
2104 : : ! =======================
2105 : : ! Get full kBZ associated to hyb%ebands
2106 : : call kpts_ibz_from_kptrlatt(cryst, hyb%ebands%kptrlatt, hyb%ebands%kptopt, hyb%ebands%nshiftk, hyb%ebands%shiftk, &
2107 [ # # # # ]: 0 : hyb%nkibz, hyb%kibz, wtk, hyb%nkbz, hyb%kbz) !, bz2ibz=bz2ibz)
2108 : : !new_kptrlatt=gwr%kptrlatt, new_shiftk=gwr%kshift,
2109 : : !bz2ibz=new%ind_qbz2ibz) # FIXME
2110 [ # # ]: 0 : ABI_FREE(wtk)
2111 : :
2112 : : ! In principle kibz should be equal to hyb%ebands%kptns.
2113 [ # # ]: 0 : ABI_CHECK_IEQ(hyb%nkibz, hyb%ebands%nkpt, "nkibz != hyb%ebands%nkpt")
2114 [ # # # # : 0 : ABI_CHECK(all(abs(hyb%ebands%kptns - hyb%kibz) < tol12), "hyb%ebands%kibz != hyb%kibz")
# # # # ]
2115 : :
2116 : : ! Note symrec convention.
2117 : 0 : ebands_kptopt = hyb%ebands%kptopt
2118 : 0 : call krank_ibz%from_kptrlatt(hyb%nkibz, hyb%kibz, hyb%ebands%kptrlatt, compute_invrank=.False.)
2119 : :
2120 [ # # # # : 0 : ABI_MALLOC(hyb%kbz2ibz, (6, hyb%nkbz))
# # # # #
# # # #
# ]
2121 [ # # ]: 0 : if (kpts_map("symrec", ebands_kptopt, cryst, krank_ibz, hyb%nkbz, hyb%kbz, hyb%kbz2ibz) /= 0) then
2122 : 0 : ABI_ERROR("Cannot map kBZ to IBZ!")
2123 : : end if
2124 : :
2125 : : ! Order kbz by stars and rearrange entries in kbz2ibz table.
2126 : 0 : call kpts_pack_in_stars(hyb%nkbz, hyb%kbz, hyb%kbz2ibz)
2127 : :
2128 [ # # ]: 0 : if (my_rank == master) then
2129 : 0 : call kpts_map_print(units, " Mapping kBZ --> kIBZ", "symrec", hyb%kbz, hyb%kibz, hyb%kbz2ibz, dtset%prtvol)
2130 : : end if
2131 : :
2132 : : ! Table with symrel conventions for the symmetrization of the wfs.
2133 [ # # # # : 0 : ABI_MALLOC(hyb%kbz2ibz_symrel, (6, hyb%nkbz))
# # # # #
# # # #
# ]
2134 [ # # ]: 0 : if (kpts_map("symrel", ebands_kptopt, cryst, krank_ibz, hyb%nkbz, hyb%kbz, hyb%kbz2ibz_symrel) /= 0) then
2135 : 0 : ABI_ERROR("Cannot map kBZ to IBZ!")
2136 : : end if
2137 : 0 : call krank_ibz%free()
2138 : :
2139 : : ! Setup qIBZ, weights and BZ.
2140 : : ! Always use q --> -q symmetry even in systems without inversion
2141 : : ! TODO: Might add input variable to rescale the q-mesh.
2142 : :
2143 : : ! Find the number of q-points such that q = k1-k2.
2144 : 0 : call findnq(hyb%nkbz, hyb%kbz, cryst%nsym, cryst%symrec, cryst%symafm, hyb%nqibz, cryst%timrev)
2145 : :
2146 : : ! Find the coordinates of the q-points in the IBZ.
2147 [ # # # # : 0 : ABI_MALLOC(hyb%qibz, (3, hyb%nqibz))
# # # # #
# # # #
# ]
2148 : 0 : call findq(hyb%nkbz, hyb%kbz, cryst%nsym, cryst%symrec, cryst%symafm, cryst%gprimd, hyb%nqibz, hyb%qibz, cryst%timrev)
2149 [ # # # # : 0 : ABI_CHECK(all(abs(hyb%qibz(:,1)) < tol16), "First qpoint in qibz should be Gamma!")
# # ]
2150 : :
2151 : : ! HM: the bz2ibz produced above is incomplete, I do it here using listkk
2152 : : !ABI_MALLOC(hyb%qbz2ibz, (6, hyb%nqbz))
2153 : : !qrank = krank_from_kptrlatt(hyb%nqibz, hyb%qibz, qptrlatt, compute_invrank=.False.)
2154 : :
2155 : : !if (kpts_map("symrec", qtimrev1, cryst, qrank, hyb%nqbz, hyb%qbz, hyb%qbz2ibz) /= 0) then
2156 : : ! ABI_ERROR("Cannot map qBZ to IBZ!")
2157 : : !end if
2158 : : !call qrank%free()
2159 : :
2160 : : ! Order qbz by stars and rearrange entries in qbz2ibz table.
2161 : : !call kpts_pack_in_stars(hyb%nqbz, hyb%qbz, hyb%qbz2ibz)
2162 : : !if (my_rank == master) then
2163 : : ! call kpts_map_print(units, " Mapping qBZ --> qIBZ", "symrec", hyb%qbz, hyb%qibz, hyb%qbz2ibz, dtset%prtvol)
2164 : : !end if
2165 : :
2166 : 0 : nqbzX = hyb%nqibz*cryst%nsym*cryst%timrev ! Maximum possible number
2167 [ # # # # : 0 : ABI_MALLOC(qbz, (3, nqbzX))
# # # # #
# # # ]
2168 [ # # # # : 0 : ABI_MALLOC(wtq, (hyb%nqibz))
# # # # #
# # # ]
2169 [ # # # # : 0 : ABI_MALLOC(qtab, (nqbzX))
# # # # #
# # # ]
2170 [ # # # # : 0 : ABI_MALLOC(qtabi, (nqbzX))
# # # # ]
2171 [ # # # # : 0 : ABI_MALLOC(qtabo, (nqbzX))
# # # # ]
2172 : :
2173 : 0 : call identk(hyb%qibz, hyb%nqibz, nqbzX, cryst%nsym, cryst%timrev, cryst%symrec, cryst%symafm, qbz, qtab, qtabi, qtabo, hyb%nqbz, wtq)
2174 : :
2175 [ # # # # : 0 : ABI_MALLOC(hyb%qbz, (3, hyb%nqibz))
# # # # #
# # # #
# ]
2176 [ # # # # : 0 : hyb%qbz = qbz(:,1:hyb%nqibz)
# # ]
2177 : :
2178 [ # # ]: 0 : ABI_FREE(qbz)
2179 [ # # ]: 0 : ABI_FREE(wtq)
2180 [ # # ]: 0 : ABI_FREE(qtab)
2181 [ # # ]: 0 : ABI_FREE(qtabi)
2182 [ # # ]: 0 : ABI_FREE(qtabo)
2183 : :
2184 : : ! TODO: MC technique does not seem to work as expected, even in the legacy code.
2185 : : call hyb%vcgen%init(cryst, hyb%ebands%kptrlatt, hyb%nkbz, hyb%nqibz, hyb%nqbz, hyb%qbz, &
2186 : 0 : dtset%rcut, dtset%gw_icutcoul, dtset%vcutgeo, comm)
2187 : :
2188 [ # # # # : 0 : end subroutine hyb_from_wfk_file
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
2189 : : !!***
2190 : :
2191 : : !----------------------------------------------------------------------
2192 : :
2193 : : !!****f* m_gwr/hyb_free
2194 : : !! NAME
2195 : : !! hyb_free
2196 : : !!
2197 : : !! FUNCTION
2198 : : !! Free dynamic memory.
2199 : : !!
2200 : : !! SOURCE
2201 : :
2202 : 0 : subroutine hyb_free(hyb)
2203 : :
2204 : : !Arguments ------------------------------------
2205 : : class(hyb_t),intent(inout) :: hyb
2206 : : ! *************************************************************************
2207 : :
2208 [ # # ]: 0 : ABI_SFREE(hyb%kibz)
2209 [ # # ]: 0 : ABI_SFREE(hyb%kbz)
2210 [ # # ]: 0 : ABI_SFREE(hyb%qibz)
2211 [ # # ]: 0 : ABI_SFREE(hyb%qbz)
2212 [ # # ]: 0 : ABI_SFREE(hyb%wtq)
2213 [ # # ]: 0 : ABI_SFREE(hyb%kbz2ibz)
2214 [ # # ]: 0 : ABI_SFREE(hyb%kbz2ibz_symrel)
2215 [ # # ]: 0 : ABI_SFREE(hyb%qbz2ibz)
2216 : :
2217 : : ! Free datatypes
2218 : 0 : call hyb%wfd%free(); call hyb%vcgen%free(); call hyb%ebands%free()
2219 : :
2220 : 0 : end subroutine hyb_free
2221 : : !!***
2222 : :
2223 : : !!****f* m_ksdiago/psbands_init
2224 : : !! NAME
2225 : : !! psbands_init
2226 : : !!
2227 : : !! FUNCTION
2228 : : !!
2229 : : !! INPUTS
2230 : : !!
2231 : : !! OUTPUT
2232 : : !!
2233 : : !! SOURCE
2234 : :
2235 : 0 : subroutine psbands_init(psb, dtset, eig_size, eig_k, gs_fermie)
2236 : :
2237 : : !Arguments ------------------------------------
2238 : : class(psbands_t),intent(out) :: psb
2239 : : class(dataset_type),target,intent(in) :: dtset
2240 : : integer,intent(in) :: eig_size
2241 : : real(dp),intent(in) :: gs_fermie
2242 : : !arrays
2243 : : real(dp),intent(in) :: eig_k(eig_size)
2244 : :
2245 : : !Local variables-------------------------------
2246 : : !scalars
2247 : : integer :: islice, ib, cnt, units(2), first_band, last_band, nb
2248 : : real(dp) :: first_eig, last_eig
2249 [ # # ]: 0 : real(dp),allocatable :: tmp_eig_k(:)
2250 : : ! *********************************************************************
2251 : :
2252 : : ABI_UNUSED(dtset%natom)
2253 : :
2254 : : ! Shift energies wrt the input Fermi level.
2255 [ # # # # : 0 : ABI_MALLOC(tmp_eig_k, (eig_size))
# # # # #
# # # ]
2256 [ # # # # ]: 0 : tmp_eig_k = eig_k - gs_fermie
2257 : :
2258 : 0 : psb%nb_protected = huge(1)
2259 : : !psb%nb_protected = dtset%nb_protected
2260 : : !psb%maxsto_per_slice = dtset%nb_per_slice
2261 : : ! TODO
2262 : 0 : psb%efrac = 0.02_dp ! dtset%efrac
2263 : :
2264 : : ! Compute nslices and subspace
2265 : : ! TODO: Add possibility of treating occupied states as well?
2266 [ # # # # : 0 : ABI_MALLOC(psb%subspace, (3, eig_size))
# # # # #
# # # #
# ]
2267 : 0 : first_band = psb%nb_protected + 1
2268 : 0 : psb%nslices = 0
2269 : :
2270 : : do while (first_band > 0)
2271 : : first_eig = tmp_eig_k(first_band)
2272 : : last_eig = first_eig + (first_eig * psb%efrac)
2273 : : last_band = get_band_with_energy_small_than(first_band+1, eig_size, last_eig)
2274 : : psb%nslices = psb%nslices + 1
2275 : : psb%subspace(1, psb%nslices) = first_band
2276 : : if (last_band == -1) then
2277 : : psb%subspace(2, psb%nslices) = eig_size
2278 : : else
2279 : : psb%subspace(2, psb%nslices) = last_band
2280 : : end if
2281 : : nb = psb%subspace(2, psb%nslices) - psb%subspace(1, psb%nslices) + 1
2282 : : if (last_band == first_band) then
2283 : : ! Won't use pseudo bands in this case.
2284 : : psb%subspace(3, psb%nslices) = 1
2285 : : else
2286 : : !psb%subspace(3, psb%nslices) = min(dtset%nb_per_slice, nb)
2287 : : end if
2288 : : first_band = last_band + 1
2289 : : !write(std_out,'(a,i0,a,*(1x,i0))')" islice: ", psb%nslices, " subspace:", psb%subspace(:, psb%nslices)
2290 : : end do
2291 : :
2292 : : ! Copy eigenvalues of the protected states.
2293 : 0 : psb%nb_tot = psb%nb_protected + sum(psb%subspace(3,1:psb%nslices))
2294 [ # # # # ]: 0 : ABI_MALLOC(psb%ps_eig, (psb%nb_tot))
2295 [ # # ]: 0 : psb%ps_eig(1:psb%nb_protected) = tmp_eig_k(1:psb%nb_protected)
2296 : :
2297 : : cnt = 0
2298 : : do islice=1,psb%nslices
2299 : : first_band = psb%subspace(1, islice)
2300 : : last_band = psb%subspace(2, islice)
2301 : : ! Take average of eigenvalues inside the slice.
2302 : : do ib=1,psb%subspace(3, islice)
2303 : : cnt = cnt + 1
2304 : : psb%ps_eig(psb%nb_protected + cnt) = sum(tmp_eig_k(first_band:last_band)) / dble(last_band - first_band + 1)
2305 : : end do
2306 : : end do
2307 [ # # ]: 0 : ABI_FREE(tmp_eig_k)
2308 : :
2309 [ # # ]: 0 : psb%ps_eig = psb%ps_eig + gs_fermie
2310 : :
2311 [ # # ]: 0 : units = [std_out, ab_out]
2312 : 0 : call wrtout(units, ' Stochastic pseudobands setup:', pre_newlines=1)
2313 : 0 : call wrtout(units, sjoin(' Number of stochastic subspaces: ', itoa(psb%nslices)))
2314 : : !call wrtout(units, sjoin(' Number of stochastic pseudobands per subspace: ', itoa(dtset%nb_per_slice)))
2315 : 0 : call wrtout(units, sjoin(' Original number of bands: ', itoa(eig_size)))
2316 : 0 : call wrtout(units, sjoin(' Number of bands in the protection window: ', itoa(psb%nb_protected)))
2317 : 0 : call wrtout(units, sjoin(' Final number of bands: ', itoa(psb%nb_tot)), newlines=1)
2318 : :
2319 : : !if (dtset%prtvol > 5) then
2320 : : ! do islice=1,psb%nslices
2321 : : ! write(msg,'(a,i0,a,*(1x,i0))')" islice: ", psb%nslices, " subspace:", psb%subspace(:, psb%nslices)
2322 : : ! call wrtout(units, msg)
2323 : : ! end do
2324 : : !end if
2325 : :
2326 : : contains
2327 : :
2328 : : integer function get_band_with_energy_small_than(idx_start, idx_end, energy) result(band)
2329 : : integer, intent(in) :: idx_start, idx_end
2330 : : integer :: ib
2331 : : real(dp), intent(in) :: energy
2332 : :
2333 : : band = -1
2334 : : do ib=idx_start,idx_end
2335 : : if (tmp_eig_k(ib) > energy) then
2336 : : band = ib - 1; return
2337 : : end if
2338 : : end do
2339 : : end function get_band_with_energy_small_than
2340 : :
2341 : : end subroutine psbands_init
2342 : : !!***
2343 : :
2344 : : !!****f* m_ksdiago/psbands_band2slice
2345 : : !! NAME
2346 : : !! psbands_band2slice
2347 : : !!
2348 : : !! FUNCTION
2349 : : !! Return the slice index from the band index. -1 if band is protected.
2350 : : !!
2351 : : !! SOURCE
2352 : :
2353 : 0 : integer function psbands_band2slice(psb, band) result(islice)
2354 : :
2355 : : !Arguments ------------------------------------
2356 : : class(psbands_t),intent(in) :: psb
2357 : : integer,intent(in) :: band
2358 : : ! *********************************************************************
2359 : :
2360 [ # # ]: 0 : do islice=1,psb%nslices
2361 [ # # # # ]: 0 : if (band >= psb%subspace(1,islice) .and. &
2362 : 0 : band <= psb%subspace(2,islice)) return
2363 : : end do
2364 : 0 : islice = -1
2365 : :
2366 : : end function psbands_band2slice
2367 : : !!***
2368 : :
2369 : : !!****f* m_ksdiago/psbands_free
2370 : : !! NAME
2371 : : !! psbands_free
2372 : : !!
2373 : : !! FUNCTION
2374 : : !! Free memory
2375 : : !!
2376 : : !! SOURCE
2377 : :
2378 : 0 : subroutine psbands_free(psb)
2379 : :
2380 : : !Arguments ------------------------------------
2381 : : class(psbands_t),intent(inout) :: psb
2382 : : ! *********************************************************************
2383 : :
2384 [ # # ]: 0 : ABI_SFREE(psb%ps_eig)
2385 [ # # ]: 0 : ABI_SFREE(psb%subspace)
2386 : :
2387 : 0 : end subroutine psbands_free
2388 : : !!***
2389 : :
2390 [ # # # # : 0 : end module m_ksdiago
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
2391 : : !!***
|
