TABLE OF CONTENTS
- ABINIT/m_gstore
- m_gstore/gqk_dbldelta_qpt
- m_gstore/gqk_free
- m_gstore/gqk_myqpt
- m_gstore/gqk_t
- m_gstore/gstore_calc_my_phonons
- m_gstore/gstore_compute
- m_gstore/gstore_distribute_spins
- m_gstore/gstore_fill_bks_mask
- m_gstore/gstore_filter_erange__
- m_gstore/gstore_filter_fs_tetra__
- m_gstore/gstore_free
- m_gstore/gstore_from_ncpath
- m_gstore/gstore_get_a2fw
- m_gstore/gstore_get_lambda_iso_iw
- m_gstore/gstore_get_mpw_gmax
- m_gstore/gstore_init
- m_gstore/gstore_malloc__
- m_gstore/gstore_print
- m_gstore/gstore_set_mpi_grid__
- m_gstore/gstore_spin2my_is
- m_gstore/gstore_t
- m_gstore/priority_from_eph_task
- m_gstore/recompute_select_qbz_spin
ABINIT/m_gstore [ Modules ]
NAME
m_gstore
FUNCTION
This module implements the gstore_t object that allows one to **precompute"" the e-ph matrix elements g and store them in memory with a MPI-distributed data structure. This approach is the most CPU-efficient one when one has to deal with algorithms in which the same g(q, k) is required several (many) times. Typical examples are iterative solvers for non-linear equations that are called inside a loop over T. At each iteration, indeed, we need g(q, k) and computing these quantities from scratch would be very expensive. Note that g depends on the two wave vectors (q, k), two electron band indices (m, n), phonon mode nu with momentum q and spin index if nsppol == 2 (collinear case). g(q, k) is therefore a sloppy notation for: g(q, k) = <k + q, m, spin| \Delta_{q, \nu} V^{spin}_{scf} | k, n, spin> There are lots of technical details that should be discussed but, roughly speaking, the gqk API allows one to: - select whether q or k should be in the IBZ or in the BZ. NB: It is not possible to use the IBZ both for q and k as g(Sk, q) = g(k, S^{-1}q) thus one has to select the appropriate zones beforehand. - filter bands and/or k/q wavevectors according to some criterion. In superconductors, for instance, only k/k+q states on the Fermi surface are needed. In semiconductors, one can include only k, k+q inside an energy window around the band edge. - whether the code should compute and store the complex valued g or |g|^2. Expression depending of |g|^2 are gauge-invariant provided that all degenerate states are summed over. On the contrary, the complex valued g is gauge-dependent and hic sunt leones. In Abinit, the g elements are computed within the same gauge by reconstructing Bloch states in the BZ from the IBZ by using a deterministic symmetrization algorithm Client code reading the e-ph matrix elements produced by ABINIT is expected to follow the same conventions, especially if one needs to mix g with wavefunctions in the BZ. At the level of the API, we have three different routines. 1) gstore_new builds the object, defines the BZ sampling type (e.g. k in the IBZ, q in the BZ) and implements filtering techniques. The MPI grid is automatically generated at this level 2) gstore_compute evaluates the e-ph matrix elements in parallel and dumps the results to GSTORE.nc 3) gstore%from_ncpath reconstructs the object from GSTORE.nc In a typical scenario, one uses eph_task 11 to generate GSTORE.nc i.e. steps 1) and 2). Then one introduces a new value of eph_task in which we read the object from file and call a specialized routine that implements the "post-processing" steps needed to compute the physical properties of interest. Now, let's discuss the MPI-distribution. The (q, k) matrix is distributed inside a 2D cartesian grid using block distribution. This is schematic representation for 4 procs with 2 procs for k and 2 procs for q: k-axis (kpt_comm) |-------------------- | | | | P00 | P01 | | | | q-axis |-------------------- (qpt_comm) | | | | P10 | P11 | | | | |-------------------- Each processor stores all the (band_1, band_2) transitions for a given (q, k) pair. Perturbations can be optionally distributed along a third axis (pert_comm). Note, however, that the parallelism over perturbations is not expected to be the most efficient although it allows one to reduce the memory required to store the scattering potential in the supercell as we can distribute W(r, R, 3 * natom) over the last dimension. For electronic properties, one usually uses k-points in the IBZ and q-points in the BZ. hence the parallelism over q-points is the most efficient one in terms of wall-time. Keep in mind, however, that the k-point parallelism allows one to reduce the memory allocated for the wavefunctions. Using some procs for k-point is also beneficial in terms of performance as we reduce load imbalance with the number of procs in qpt_comm does not divide nqbz. NB: If nsppol == 2, we create two gqk objects, one for each spin. The reason is that dimensions such as the number of effective bands/q-points/k-points depends on spin if filters are employed.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
97 #if defined HAVE_CONFIG_H 98 #include "config.h" 99 #endif 100 101 #include "abi_common.h" 102 103 module m_gstore 104 105 use defs_basis 106 use m_abicore 107 use m_xmpi 108 use m_errors 109 use m_htetra 110 use libtetrabz 111 use m_ebands 112 use, intrinsic :: iso_c_binding 113 use netcdf 114 use m_nctk 115 use m_wfk 116 use m_ddb 117 use m_ddk 118 use m_dvdb 119 use m_crystal 120 use m_fft 121 use m_hamiltonian 122 use m_pawcprj 123 use m_dtset 124 use m_dtfil 125 use m_wfd 126 use m_ephtk 127 use m_mkffnl 128 129 use defs_abitypes, only : mpi_type 130 use m_time, only : cwtime, cwtime_report, sec2str 131 use m_fstrings, only : tolower, itoa, ftoa, sjoin, ktoa, ltoa, strcat, replace_ch0 132 !use m_yaml, only : yamldoc_t 133 use m_numeric_tools, only : arth, get_diag, isdiagmat 134 use m_krank, only : krank_t, krank_new, krank_from_kptrlatt, get_ibz2bz, star_from_ibz_idx 135 use m_io_tools, only : iomode_from_fname 136 use m_special_funcs, only : gaussian 137 use m_copy, only : alloc_copy 138 use m_fftcore, only : ngfft_seq, get_kg 139 use m_cgtools, only : cg_zdotc 140 use m_kg, only : getph, mkkpg 141 use defs_datatypes, only : ebands_t, pseudopotential_type 142 use m_hdr, only : hdr_type, fform_from_ext, hdr_ncread 143 use m_symtk, only : matr3inv 144 use m_kpts, only : kpts_ibz_from_kptrlatt, kpts_timrev_from_kptopt, kpts_map, kpts_sort, kpts_pack_in_stars 145 use m_getgh1c, only : getgh1c, rf_transgrid_and_pack, getgh1c_setup 146 use m_ifc, only : ifc_type 147 use m_phonons, only : pheigvec_rotate 148 use m_pawang, only : pawang_type 149 use m_pawrad, only : pawrad_type 150 use m_pawtab, only : pawtab_type 151 use m_pawfgr, only : pawfgr_type 152 153 implicit none 154 155 private
m_gstore/gqk_dbldelta_qpt [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gqk_dbldelta_qpt
FUNCTION
Note that k/q weights are not included in dbldelta_q
INPUTS
OUTPUT
SOURCE
2281 subroutine gqk_dbldelta_qpt(gqk, my_iq, gstore, eph_intmeth, eph_fsmear, qpt, weight_q, dbldelta_q) 2282 2283 !Arguments ------------------------------------ 2284 class(gqk_t),intent(in) :: gqk 2285 class(gstore_t),target,intent(inout) :: gstore 2286 integer,intent(in) :: my_iq, eph_intmeth 2287 real(dp),intent(in) :: eph_fsmear 2288 real(dp),intent(out) :: qpt(3), weight_q, dbldelta_q(gqk%nb, gqk%nb, gqk%my_nk) 2289 2290 !Local variables ------------------------------ 2291 !scalars 2292 real(dp), parameter :: min_smear = tol9 2293 integer :: nb, nkbz, spin, my_ik, ib, ib1, ib2, band1, band2, nesting, ebands_timrev 2294 integer :: ik_ibz, isym_k, trev_k, tsign_k, g0_k(3) 2295 integer :: ikq_ibz, isym_kq, trev_kq, tsign_kq, g0_kq(3) 2296 integer :: ii, i1, i2, i3, cnt, ik_bz, ltetra 2297 real(dp) :: g1, g2, sigma !, weight_k !, cpu, wall, gflops 2298 logical :: isirr_k, isirr_kq, use_adaptive 2299 type(ebands_t), pointer :: ebands 2300 type(crystal_t), pointer :: cryst 2301 type(krank_t) :: my_krank 2302 !arrays 2303 integer :: nge(3), ngw(3) 2304 integer,allocatable :: my_kqmap(:,:), kmesh_map(:,:) 2305 real(dp) :: kk(3), kmesh_cartvec(3,3), rlatt(3,3), klatt(3,3), vb_k(3, gqk%nb), vb_kq(3, gqk%nb) 2306 real(dp),allocatable :: eig_k(:,:), eig_kq(:,:), kmesh(:,:), wght_bz(:,:,:) 2307 2308 !---------------------------------------------------------------------- 2309 2310 nb = gqk%nb; nkbz = gstore%nkbz; spin = gqk%spin 2311 2312 ebands => gstore%ebands; cryst => gstore%cryst 2313 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 2314 2315 call gqk%myqpt(my_iq, gstore, weight_q, qpt) 2316 2317 ! The double delta with tetra is ill-defined for q == 0. In this case we fall back to gaussian. 2318 nesting = merge(1, 0, abs(eph_intmeth) == 2 .and. all(abs(qpt) < tol12)) 2319 2320 rlatt = gstore%ebands%kptrlatt; call matr3inv(rlatt, klatt) 2321 kmesh_cartvec(:, 1) = cryst%gprimd(:,1)*klatt(1,1) + cryst%gprimd(:,2)*klatt(2,1) + cryst%gprimd(:,3)*klatt(3,1) 2322 kmesh_cartvec(:, 2) = cryst%gprimd(:,1)*klatt(1,2) + cryst%gprimd(:,2)*klatt(2,2) + cryst%gprimd(:,3)*klatt(3,2) 2323 kmesh_cartvec(:, 3) = cryst%gprimd(:,1)*klatt(1,3) + cryst%gprimd(:,2)*klatt(2,3) + cryst%gprimd(:,3)*klatt(3,3) 2324 ! TODO: It seems that two_pi is not needed here! 2325 2326 if (abs(eph_intmeth) == 1 .or. nesting /= 0) then 2327 use_adaptive = eph_fsmear < zero .or. abs(eph_intmeth) == 2 2328 if (use_adaptive) then 2329 ABI_CHECK(allocated(gqk%vk_cart_ibz), "vk_cart_ibz should be allocated when use_adaptive is .True.") 2330 end if 2331 2332 ! Find k + q in the IBZ for all my k-points. 2333 ABI_MALLOC(my_kqmap, (6, gqk%my_nk)) 2334 if (kpts_map("symrel", ebands_timrev, cryst, gstore%krank_ibz, gqk%my_nk, gqk%my_kpts, my_kqmap, qpt=qpt) /= 0) then 2335 ABI_ERROR(sjoin("Cannot map k+q to IBZ with qpt:", ktoa(qpt))) 2336 end if 2337 2338 ! Init default sigma 2339 sigma = eph_fsmear 2340 2341 do my_ik=1,gqk%my_nk 2342 kk = gqk%my_kpts(:, my_ik) 2343 2344 ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik) 2345 trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5, my_ik) 2346 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 2347 tsign_k = 1; if (trev_k == 1) tsign_k = -1 2348 2349 ikq_ibz = my_kqmap(1, my_ik); isym_kq = my_kqmap(2, my_ik) 2350 trev_kq = my_kqmap(6, my_ik); g0_kq = my_kqmap(3:5, my_ik) 2351 isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0)) 2352 tsign_kq = 1; if (trev_kq == 1) tsign_kq = -1 2353 2354 if (use_adaptive) then 2355 ! If k or k+q is not in the IBZ, we need to recostruct the value by symmetry using v(Sq) = S v(q). 2356 ! Use transpose(R) because we are using the tables for the wavefunctions 2357 ! In this case listkk has been called with symrec and use_symrec=False 2358 ! so q_bz = S^T q_ibz where S is the isym_kq symmetry 2359 vb_k = gqk%vk_cart_ibz(:,:,ik_ibz) 2360 vb_kq = gqk%vk_cart_ibz(:,:,ikq_ibz) 2361 2362 if (.not. isirr_k) then 2363 do ib=1,nb 2364 vb_k(:,ib) = tsign_k * matmul(transpose(cryst%symrel_cart(:,:,isym_k)), vb_k(:,ib)) 2365 end do 2366 end if 2367 if (.not. isirr_kq) then 2368 do ib=1,nb 2369 vb_kq(:,ib) = tsign_kq * matmul(transpose(cryst%symrel_cart(:,:,isym_kq)), vb_kq(:,ib)) 2370 end do 2371 end if 2372 end if 2373 2374 do ib2=1,nb 2375 band2 = ib2 + gqk%bstart - 1 2376 if (use_adaptive) then 2377 sigma = max(maxval([(abs(dot_product(vb_k(:, ib2), kmesh_cartvec(:,ii))), ii=1,3)]), min_smear) 2378 !write(std_out, *)"sigma:", sigma * Ha_eV 2379 end if 2380 g2 = gaussian(ebands%eig(band2, ik_ibz, spin) - ebands%fermie, sigma) 2381 2382 do ib1=1,nb 2383 band1 = ib1 + gqk%bstart - 1 2384 if (use_adaptive) then 2385 sigma = max(maxval([(abs(dot_product(vb_kq(:, ib1), kmesh_cartvec(:,ii))), ii=1,3)]), min_smear) 2386 end if 2387 g1 = gaussian(ebands%eig(band1, ikq_ibz, spin) - ebands%fermie, sigma) 2388 dbldelta_q(ib1, ib2, my_ik) = g1 * g2 ! / fs%nktot 2389 end do 2390 2391 end do 2392 end do 2393 2394 ABI_FREE(my_kqmap) 2395 2396 else if (abs(eph_intmeth) == 2) then 2397 2398 ABI_CHECK(isdiagmat(ebands%kptrlatt), "kptrlatt must be diagonal when tetra is used.") 2399 ABI_CHECK(ebands%nshiftk == 1, "nshiftk must be 1 when tetra is used") 2400 nge = get_diag(ebands%kptrlatt); ngw = nge 2401 ABI_CHECK(nkbz == product(nge(1:3)), "Wrong nge") 2402 2403 ! Compute eig_k and eig_kq in full BZ for the relevant bands around Ef. 2404 ABI_MALLOC(kmesh, (3, nkbz)) 2405 ABI_MALLOC(eig_k, (nb, nkbz)) 2406 ABI_MALLOC(eig_kq, (nb, nkbz)) 2407 2408 ! Technical problems: 2409 ! 2410 ! 1) libtetrabz works with the BZ and assume a certaing ordering of the k-points (see below) 2411 ! so we have to fill the array with eig_k and eig_kq from the IBZ by remapping the libtetra kk 2412 ! to the Abinit IBZ 2413 2414 ! 2) The dbldelta weights are given in the BZ, while the caller requires weights for k in the IBZ 2415 ! and moreover only for the IBZ k-point treated by this MPI proc. 2416 2417 ik_bz = 0 2418 do i3=0,nge(3) - 1 2419 do i2=0,nge(2) - 1 2420 do i1=0,nge(1) - 1 2421 ik_bz = ik_bz + 1 2422 kk = ([i1, i2, i3] + ebands%shiftk(:, 1)) / nge(:) 2423 kmesh(:, ik_bz) = kk 2424 end do 2425 end do 2426 end do 2427 2428 ! Map libtetra BZ mesh to IBZ and fill eig_k 2429 !call cwtime(cpu, wall, gflops, "start") 2430 ABI_MALLOC(kmesh_map, (6, nkbz)) 2431 2432 ! Find correspondence between libtetra mesh and the IBZ. 2433 if (kpts_map("symrec", ebands_timrev, cryst, gstore%krank_ibz, nkbz, kmesh, kmesh_map) /= 0) then 2434 ABI_ERROR("Cannot map libtetra mesh to IBZ") 2435 end if 2436 2437 do ik_bz=1,nkbz 2438 ik_ibz = kmesh_map(1, ik_bz) 2439 eig_k(:, ik_bz) = ebands%eig(gqk%bstart:gqk%bstop, ik_ibz, spin) - ebands%fermie 2440 end do 2441 2442 ! Map libtetra BZ mesh + q to IBZ and fill eig_kq 2443 if (kpts_map("symrec", ebands_timrev, cryst, gstore%krank_ibz, nkbz, kmesh, kmesh_map, qpt=qpt) /= 0) then 2444 ABI_ERROR(sjoin("Cannot map libtetra k+q to IBZ with qpt:", ktoa(qpt))) 2445 end if 2446 2447 do ik_bz=1,nkbz 2448 ikq_ibz = kmesh_map(1, ik_bz) 2449 eig_kq(:, ik_bz) = ebands%eig(gqk%bstart:gqk%bstop, ikq_ibz, spin) - ebands%fermie 2450 end do 2451 2452 ABI_FREE(kmesh_map) 2453 !call cwtime_report(" kmesh_map", cpu, wall, gflops) 2454 2455 ! Call libtetra routine to compute weights for double delta integration. 2456 ! Note that libtetra assumes Ef set to zero. 2457 ! TODO: Average weights over degenerate states? 2458 ! NB: This is a bootleneck, can pass comm_kp 2459 2460 ! Select option for double delta with tetra. 2461 ! 2 for the optimized tetrahedron method. 2462 ! -2 for the linear tetrahedron method. 2463 ltetra = 0 2464 if (eph_intmeth == 2) ltetra = 2 2465 if (eph_intmeth == -2) ltetra = 1 2466 2467 ABI_MALLOC(wght_bz, (nb, nb, nkbz)) 2468 call libtetrabz_dbldelta(ltetra, gstore%cryst%gprimd, nb, nge, eig_k, eig_kq, ngw, wght_bz) !, comm=comm) 2469 !call cwtime_report(" libtetrabz_dbldelta", cpu, wall, gflops) 2470 2471 my_krank = krank_new(gqk%my_nk, gqk%my_kpts) 2472 2473 ! Reindex from full BZ to my set of kpoints and rescale weights. 2474 cnt = 0 2475 do ik_bz=1,nkbz 2476 my_ik = my_krank%get_index(kmesh(:, ik_bz)) 2477 if (my_ik /= -1) then 2478 dbldelta_q(:,:,my_ik) = wght_bz(:,:,ik_bz) * gstore%nkbz 2479 cnt = cnt + 1 2480 end if 2481 end do 2482 2483 ! FIXME: bug if k-point (and q-point) parallelism. 2484 ABI_CHECK(cnt == gqk%my_nk, sjoin("cnt != my_nk, ", itoa(cnt), itoa(gqk%my_nk))) 2485 call my_krank%free() 2486 !call cwtime_report(" transfer", cpu, wall, gflops) 2487 2488 ABI_FREE(wght_bz) 2489 ABI_FREE(kmesh) 2490 ABI_FREE(eig_k) 2491 ABI_FREE(eig_kq) 2492 2493 else 2494 ABI_ERROR(sjoin("Invalid eph_intmeth:", itoa(eph_intmeth))) 2495 end if 2496 2497 end subroutine gqk_dbldelta_qpt
m_gstore/gqk_free [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gqk_free
FUNCTION
Free dynamic memory.
SOURCE
2510 subroutine gqk_free(gqk) 2511 2512 !Arguments ------------------------------------ 2513 class(gqk_t),intent(inout) :: gqk 2514 !---------------------------------------------------------------------- 2515 2516 ABI_SFREE(gqk%my_k2ibz) 2517 ABI_SFREE(gqk%my_kpts) 2518 ABI_SFREE(gqk%my_wtk) 2519 ABI_SFREE(gqk%my_q2ibz) 2520 ABI_SFREE(gqk%my_q2bz) 2521 ABI_SFREE(gqk%my_wnuq) 2522 ABI_SFREE(gqk%my_displ_cart) 2523 ABI_SFREE(gqk%my_g) 2524 ABI_SFREE(gqk%my_g2) 2525 ABI_SFREE(gqk%my_iperts) 2526 ABI_SFREE(gqk%vk_cart_ibz) 2527 ABI_SFREE(gqk%vkmat_cart_ibz) 2528 2529 ! Free communicators 2530 call gqk%kpt_comm%free() 2531 call gqk%qpt_comm%free() 2532 call gqk%pert_comm%free() 2533 call gqk%qpt_pert_comm%free() 2534 call gqk%grid_comm%free() 2535 2536 end subroutine gqk_free
m_gstore/gqk_myqpt [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gqk_myqpt
FUNCTION
Return the weight and the reduced coordinates of the q-point from the local index my_iq
INPUTS
OUTPUT
SOURCE
2233 pure subroutine gqk_myqpt(gqk, my_iq, gstore, weight, qpt) 2234 2235 !Arguments ------------------------------------ 2236 class(gqk_t),intent(in) :: gqk 2237 class(gstore_t),intent(in) :: gstore 2238 integer,intent(in) :: my_iq 2239 real(dp),intent(out) :: weight, qpt(3) 2240 2241 !Local variables ------------------------------ 2242 !scalars 2243 integer :: iq_ibz, isym_q, trev_q, tsign_q, g0_q(3) 2244 logical :: isirr_q 2245 2246 !---------------------------------------------------------------------- 2247 2248 iq_ibz = gqk%my_q2ibz(1, my_iq); isym_q = gqk%my_q2ibz(2, my_iq) 2249 trev_q = gqk%my_q2ibz(6, my_iq); g0_q = gqk%my_q2ibz(3:5, my_iq) 2250 isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0)) 2251 tsign_q = 1; if (trev_q == 1) tsign_q = -1 2252 2253 ! NB: Use symrec convention for q 2254 qpt = tsign_q * matmul(gstore%cryst%symrec(:,:,isym_q), gstore%qibz(:, iq_ibz)) + g0_q 2255 2256 select case(gstore%qzone) 2257 case ("ibz") 2258 weight = gstore%wtq(iq_ibz) 2259 case ("bz") 2260 weight = one / gstore%nqbz 2261 end select 2262 2263 end subroutine gqk_myqpt
m_gstore/gqk_t [ Types ]
[ Top ] [ m_gstore ] [ Types ]
NAME
gqk_t
FUNCTION
This object stores MPI-distributed e-ph matrix elements for a given spin index (if collinear magnetism i.e. nsppol 2). local dimensions and arrays start with `my_`, global dimensions start with `glob_
SOURCE
172 type, public :: gqk_t 173 174 integer :: cplex = -1 175 ! 1 if |g|^2 should be stored 176 ! 2 if complex valued g (mind the gauge) 177 178 integer :: spin = -1 179 ! Spin index. 180 181 integer :: natom3 = -1 182 ! 3 * natom 183 ! Mainly used to dimension arrays 184 185 integer :: nb = -1 186 ! Number of bands included in the calculation for this spin 187 ! Global as this dimension is not MPI-distributed due to (m, n) pairs. 188 ! NB: nb is not necessarily equal to nband. 189 190 integer :: bstart = -1, bstop = -1 191 ! The first band starts at bstart. 192 ! The last band is bstop (global indices) 193 194 integer :: my_npert = -1 195 ! Number of perturbations treated by this MPI rank. 196 197 integer :: glob_nk = -1, glob_nq = -1 198 ! Total number of k/q points in global matrix. 199 ! Note that k-points/q-points can be filtered. 200 ! Use kzone, qzone and kfilter to interpret these dimensions. 201 202 integer :: my_nk = -1, my_nq = -1 203 ! Number of k/q points treated by this MPI proc. 204 ! Used to loop and allocate local arrays. 205 206 integer :: my_kstart = -1, my_qstart = -1 207 ! Index of the first k/q point in the global matrix treated by this MPI proc 208 209 integer,allocatable :: my_k2ibz(:,:) 210 ! (6, my_nk) 211 ! Mapping my_kpoints --> kibz 212 213 !integer,allocatable :: my_k2bz(:,:) 214 ! (my_nk) 215 ! Mapping my_kpoints --> ik_bz 216 217 real(dp),allocatable :: my_kpts(:,:) 218 ! (3, my_nkpt) 219 ! k-points treated by this MPI proc. 220 221 real(dp),allocatable :: my_wtk(:) 222 ! (my_nkpt) 223 ! k-point weights 224 225 integer,allocatable :: my_q2ibz(:,:) 226 ! (6, my_nq) 227 ! Mapping my_qpoints --> qibz 228 229 integer,allocatable :: my_q2bz(:) 230 ! (my_nq) 231 ! Mapping my_iq index --> iq_bz index in the full BZ 232 233 !integer,allocatable :: my_q2glob(:) 234 ! (my_nq) 235 ! Mapping my_iq index --> global index in the g(q, k) matrix. 236 237 real(dp),allocatable :: vk_cart_ibz(:,:,:) 238 ! (3, nb, nkibz) 239 ! Diagonal v_{m, m,k} for k in the IBZ. 240 ! Values in the BZ can be reconstructed by symmetry. 241 ! Allocated if gstore%with_vk == 1 242 243 real(dp),allocatable :: vkmat_cart_ibz(:,:,:,:,:) 244 ! (3, nb, nb, nkibz) 245 ! v_{m, n,k} for the k in the IBZ 246 ! Allocated if gstore%with_vk in (1, 2) 247 248 integer,allocatable :: my_iperts(:) 249 ! (my_npert) 250 ! List of perturbation indices treated by this MPI proc. 251 ! Contiguous indices. 252 253 complex(dp), allocatable :: my_g(:,:,:,:,:) 254 ! (my_npert, nb, my_nq, nb, my_nk) 255 ! ( p, b1, q, b2, k) --> <k+q, b1| D_{q,p}H |k, b2> 256 ! e-ph matrix elements g (local buffer). Allocated if cplex == 2 257 258 real(dp), allocatable :: my_g2(:,:,:,:,:) 259 ! (my_npert, nb, my_nq, nb, my_nk) 260 ! |g|^2 (local buffer). Allocated if cplex == 1 261 262 integer :: coords_qkp(3) 263 ! Coordinates of this processor in the (q, k, pert) Cartesian grid. 264 265 type(xcomm_t) :: kpt_comm 266 ! MPI communicator over k-points 267 268 type(xcomm_t) :: qpt_comm 269 ! MPI communicator over q-points 270 271 type(xcomm_t) :: pert_comm 272 ! MPI communicator over atomic perturbations. 273 274 type(xcomm_t) :: qpt_pert_comm 275 ! MPI communicator over the 2d grid (qpt, atomic perturbations) 276 277 type(xcomm_t) :: grid_comm 278 ! MPI communicator for full (q,k,pert) grid. 279 280 real(dp),allocatable :: my_wnuq(:,:) 281 ! (my_npert, my_nq) 282 ! Phonon frequencies in Ha (MPI distributed) 283 284 real(dp),allocatable :: my_displ_cart(:,:,:,:,:) 285 ! (2, 3, cryst%natom, my_npert, my_nq)) 286 ! Phonon displacements (MPI distributed) 287 288 contains 289 290 procedure :: myqpt => gqk_myqpt 291 ! Return the q-point and the weight from my local index my_iq 292 293 !procedure :: my_qweight => gqk_my_qweight 294 ! Return the q-point weight from my local index my_iq 295 296 procedure :: dbldelta_qpt => gqk_dbldelta_qpt 297 ! Compute weights for the double delta. 298 299 procedure :: free => gqk_free 300 ! Free memory 301 302 end type gqk_t
m_gstore/gstore_calc_my_phonons [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_calc_my_phonons
FUNCTION
Compute and store ph frequencies and optionally the displacements for all the q-points treated by this MPI rank. TODO: Rewrite this part, Read stuff from the IBZ and symmetrize.
INPUTS
OUTPUT
SOURCE
1948 subroutine gstore_calc_my_phonons(gstore, store_phdispl) 1949 1950 !Arguments ------------------------------------ 1951 class(gstore_t),target,intent(inout) :: gstore 1952 logical,intent(in) :: store_phdispl 1953 1954 !Local variables------------------------------- 1955 integer :: my_is, my_iq, ierr 1956 real(dp) :: weight_q, cpu, wall, gflops, qpt(3) 1957 type(crystal_t),pointer :: cryst 1958 type(gqk_t),pointer :: gqk 1959 real(dp),allocatable :: phfrq(:), displ_cart(:,:,:,:) 1960 1961 !---------------------------------------------------------------------- 1962 1963 call cwtime(cpu, wall, gflops, "start") 1964 cryst => gstore%cryst 1965 1966 ABI_MALLOC(displ_cart, (2, 3, cryst%natom, 3 * cryst%natom)) 1967 ABI_MALLOC(phfrq, (3 * cryst%natom)) 1968 1969 do my_is=1,gstore%my_nspins 1970 gqk => gstore%gqk(my_is) 1971 1972 ABI_RECALLOC(gqk%my_wnuq, (gqk%my_npert, gqk%my_nq)) 1973 if (store_phdispl) then 1974 ABI_RECALLOC(gqk%my_displ_cart, (2, 3, cryst%natom, gqk%my_npert, gqk%my_nq)) 1975 end if 1976 1977 do my_iq=1,gqk%my_nq 1978 if (gqk%kpt_comm%skip(my_iq)) cycle ! MPI parallelism inside kpt_comm 1979 1980 ! Get phonon frequencies and eigenvectors for this q-point. 1981 ! FIXME: Perhaps it's better if we compute everything at q_ibz and then rotate 1982 ! so that we are guaranteed to have e(-q) == e(q)^* 1983 call gqk%myqpt(my_iq, gstore, weight_q, qpt) 1984 call gstore%ifc%fourq(cryst, qpt, phfrq, displ_cart) 1985 1986 gqk%my_wnuq(:, my_iq) = phfrq(gqk%my_iperts(:)) 1987 if (store_phdispl) gqk%my_displ_cart(:,:,:,:,my_iq) = displ_cart(:,:,:,gqk%my_iperts(:)) 1988 end do 1989 1990 ! Note that the computation is replicated across the perturbation communicator. 1991 ! This should not represent a problem since the phase in the displacement is fixed in dfpt_phfrq 1992 ! hence results should be coherent across different procs. 1993 call xmpi_sum(gqk%my_wnuq, gqk%kpt_comm%value, ierr) 1994 if (store_phdispl) call xmpi_sum(gqk%my_displ_cart, gqk%kpt_comm%value, ierr) 1995 end do ! my_is 1996 1997 ABI_FREE(displ_cart) 1998 ABI_FREE(phfrq) 1999 2000 call cwtime_report(" gstore_calc_my_phonons", cpu, wall, gflops) 2001 2002 end subroutine gstore_calc_my_phonons
m_gstore/gstore_compute [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_compute
FUNCTION
Compute MPI-distributed e-ph matrix elements
INPUTS
wk0_path=String with the path to the GS unperturbed WFK file. ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes. dtset<dataset_type>=All input variables for this dataset. ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) dvdb<dbdb_type>=Database with the DFPT SCF potentials. ifc<ifc_type>=interatomic force constants and corresponding real space grid info. pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. comm=MPI communicator.
OUTPUT
SOURCE
2566 subroutine gstore_compute(gstore, wfk0_path, ngfft, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 2567 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 2568 2569 !Arguments ------------------------------------ 2570 !scalars 2571 class(gstore_t),target,intent(inout) :: gstore 2572 character(len=*),intent(in) :: wfk0_path 2573 integer,intent(in) :: comm 2574 type(dataset_type),intent(in) :: dtset 2575 type(crystal_t),intent(in) :: cryst 2576 type(ebands_t),intent(in) :: ebands 2577 type(dvdb_t),intent(inout) :: dvdb 2578 type(pawang_type),intent(in) :: pawang 2579 type(pseudopotential_type),intent(in) :: psps 2580 type(pawfgr_type),intent(in) :: pawfgr 2581 type(ifc_type),intent(in) :: ifc 2582 type(mpi_type),intent(in) :: mpi_enreg 2583 !arrays 2584 integer,intent(in) :: ngfft(18),ngfftf(18) 2585 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 2586 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 2587 2588 !Local variables ------------------------------ 2589 !scalars 2590 integer,parameter :: tim_getgh1c = 1, berryopt0 = 0, ider0 = 0, idir0 = 0, LOG_MODQ = 5 2591 integer,parameter :: useylmgr = 0, useylmgr1 = 0, master = 0, ndat1 = 1 2592 integer :: my_rank,nproc,mband,nsppol,nkibz,idir,ipert,ebands_timrev, iq_bz 2593 integer :: cplex,natom,natom3,ipc,nspinor, nskip_tetra_kq 2594 integer :: bstart_k,bstart_kq,nband_k,nband_kq,band_k, ib_k, ib_kq !ib1,ib2, band_kq, 2595 integer :: ik_ibz,ikq_ibz,isym_k,isym_kq,trev_k,trev_kq 2596 integer :: my_ik, my_is, comm_rpt, my_npert, my_ip, my_iq, spin,istwf_k,istwf_kq,npw_k,npw_kq 2597 integer :: mpw, nb,ierr,cnt, n1,n2,n3,n4,n5,n6,nspden,ndone 2598 integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1 2599 integer :: nfft,nfftf,mgfft,mgfftf, nkpg, nkpg1, qbuf_size, iqbuf_cnt, root_ncid, spin_ncid, ncerr 2600 integer :: ii, my_nqibz, iq_start, iq_ibz, isym_q, trev_q, prev_iqbz 2601 real(dp) :: cpu, wall, gflops, cpu_q, wall_q, gflops_q, cpu_all, wall_all, gflops_all 2602 real(dp) :: ecut, eshift, eig0nk, weight_q, weight_k 2603 logical :: gen_eigenpb, isirr_k, isirr_kq, isirr_q, print_time 2604 type(wfd_t) :: wfd 2605 type(gs_hamiltonian_type) :: gs_hamkq 2606 type(rf_hamiltonian_type) :: rf_hamkq 2607 type(ddkop_t) :: ddkop 2608 type(gqk_t),pointer :: gqk 2609 character(len=500) :: msg 2610 !arrays 2611 integer :: g0_k(3), g0_kq(3), g0_q(3), work_ngfft(18),gmax(3),indkk_kq(6,1) 2612 integer(i1b),allocatable :: itreat_qibz(:) 2613 integer,allocatable :: kg_k(:,:), kg_kq(:,:), nband(:,:), wfd_istwfk(:), qselect(:) 2614 integer,allocatable :: my_pinfo(:,:), pert_table(:,:) 2615 integer,allocatable :: iq_buf(:,:), done_qbz_spin(:,:), my_iqibz_inds(:) 2616 real(dp) :: kk_bz(3),kq_bz(3),kk_ibz(3),kq_ibz(3), qq_bz(3), qq_ibz(3), vk(3) 2617 real(dp) :: phfrq(3*cryst%natom), ylmgr_dum(1,1,1) 2618 real(dp),allocatable :: displ_cart_qibz(:,:,:,:), displ_red_qibz(:,:,:,:), pheigvec_qibz(:,:,:,:) 2619 real(dp),allocatable :: displ_cart_qbz(:,:,:,:), displ_red_qbz(:,:,:,:), pheigvec_qbz(:,:,:,:) 2620 real(dp),allocatable :: grad_berry(:,:), kinpw1(:), kpg1_k(:,:), kpg_k(:,:), dkinpw(:) 2621 real(dp),allocatable :: ffnlk(:,:,:,:), ffnl1(:,:,:,:), ph3d(:,:,:), ph3d1(:,:,:) 2622 real(dp),allocatable :: v1scf(:,:,:,:), gkk_atm(:,:,:,:),gkq_nu(:,:,:,:) 2623 real(dp),allocatable :: bras_kq(:,:,:), kets_k(:,:,:), h1kets_kq(:,:,:), cgwork(:,:) 2624 real(dp),allocatable :: ph1d(:,:), vlocal(:,:,:,:), vlocal1(:,:,:,:,:) 2625 real(dp),allocatable :: ylm_kq(:,:), ylm_k(:,:), ylmgr_kq(:,:,:) 2626 real(dp),allocatable :: dummy_vtrial(:,:), gvnlx1(:,:), work(:,:,:,:) 2627 real(dp),allocatable :: gs1c(:,:), vk_cart_ibz(:,:,:) !, vkmat_cart_ibz(:,:,:,:) 2628 real(dp),allocatable :: my_gbuf(:,:,:,:,:,:), buf_wqnu(:,:), buf_eigvec_cart(:,:,:,:,:) 2629 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 2630 type(pawcprj_type),allocatable :: cwaveprj0(:,:) 2631 2632 !************************************************************************ 2633 2634 ! This parameter defines the size of the q-buffer used to store the g(k, q) e-ph matrix elements 2635 ! for all the k-point treated by this MPI rank. 2636 ! Increasing the buffer size increases the memory requirements 2637 ! but it leads to better performance as the number of IO operations is decreased. 2638 ! TODO: Should compute it on the basis of my_nkpt and my_nqpt 2639 qbuf_size = 4 2640 call wrtout(std_out, sjoin(" Begin computation of e-ph matrix elements with qbuf_size:", itoa(qbuf_size))) 2641 2642 if (psps%usepaw == 1) then 2643 ABI_ERROR("PAW not implemented") 2644 ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/)) 2645 end if 2646 2647 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 2648 call cwtime(cpu_all, wall_all, gflops_all, "start") 2649 2650 ! Copy important dimensions 2651 natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden 2652 nkibz = ebands%nkpt; mband = ebands%mband 2653 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 2654 2655 ! FFT meshes 2656 nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3)) 2657 nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3)) 2658 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3) 2659 n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6) 2660 2661 ! Open the DVDB file 2662 call dvdb%open_read(ngfftf, xmpi_comm_self) 2663 2664 do my_is=1,gstore%my_nspins 2665 spin = gstore%my_spins(my_is) 2666 gqk => gstore%gqk(my_is) 2667 if (gqk%pert_comm%nproc > 1) then 2668 ! Activate parallelism over perturbations 2669 ! Build table with list of perturbations treated by this CPU inside pert_comm 2670 ABI_WARNING("pert_comm%nproc > 1 not tested") 2671 call ephtk_set_pertables(cryst%natom, my_npert, pert_table, my_pinfo, gqk%pert_comm%value) 2672 call dvdb%set_pert_distrib(my_npert, natom3, my_pinfo, pert_table, gqk%pert_comm%value) 2673 ABI_CHECK(all(my_pinfo(3, :) == gqk%my_iperts), "my_pinfo(3, :) != gqk%my_iperts") 2674 ABI_FREE(my_pinfo) 2675 ABI_FREE(pert_table) 2676 end if 2677 end do 2678 2679 !call wrtout([std_out, ab_out], " Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.") 2680 ! Prepare Fourier interpolation of DFPT potentials. 2681 comm_rpt = xmpi_comm_self 2682 !comm_rpt = bqs_comm%value 2683 call dvdb%ftinterp_setup(dtset%ddb_ngqpt, gstore%qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt) 2684 2685 ! Build q-cache in the *dense* IBZ using the global mask qselect and itreat_qibz. 2686 ABI_MALLOC(itreat_qibz, (gstore%nqibz)) 2687 ABI_MALLOC(qselect, (gstore%nqibz)) 2688 qselect = 0; itreat_qibz = 0 2689 call dvdb%ftqcache_build(nfftf, ngfftf, gstore%nqibz, gstore%qibz, dtset%dvdb_qcache_mb, qselect, itreat_qibz, gstore%comm) 2690 ABI_FREE(itreat_qibz) 2691 ABI_FREE(qselect) 2692 2693 ! Initialize the wave function descriptor. 2694 ! Only wavefunctions for the symmetrical imagine of the k/k+q wavevectors treated by this MPI rank are stored. 2695 2696 ABI_MALLOC(nband, (nkibz, nsppol)) 2697 ABI_MALLOC(bks_mask, (mband, nkibz, nsppol)) 2698 ABI_MALLOC(keep_ur, (mband, nkibz, nsppol)) 2699 nband = mband; bks_mask = .False.; keep_ur = .False. 2700 2701 call gstore%fill_bks_mask(mband, nkibz, nsppol, bks_mask) 2702 2703 ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts) 2704 ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1. 2705 ABI_MALLOC(wfd_istwfk, (nkibz)) 2706 wfd_istwfk = 1 2707 ecut = dtset%ecut 2708 2709 call wfd_init(wfd, cryst, pawtab, psps, keep_ur, mband, nband, nkibz, nsppol, bks_mask,& 2710 nspden, nspinor, ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,& 2711 dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm) 2712 2713 call wfd%print(header="Wavefunctions for GSTORE calculation") 2714 2715 ABI_FREE(nband) 2716 ABI_FREE(keep_ur) 2717 ABI_FREE(wfd_istwfk) 2718 ABI_FREE(bks_mask) 2719 2720 ! Read wavefunctions. 2721 call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path)) 2722 2723 ! one-dimensional structure factor information on the coarse grid. 2724 ABI_MALLOC(ph1d, (2, 3*(2*mgfft+1)*natom)) 2725 call getph(cryst%atindx, natom, n1, n2, n3, ph1d, cryst%xred) 2726 2727 ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ. 2728 ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work 2729 ! that will be used to symmetrize the wavefunctions in G-space. 2730 call gstore%get_mpw_gmax(ecut, mpw, gmax) 2731 2732 ! Init work_ngfft 2733 gmax = gmax + 4 ! FIXME: this is to account for umklapp 2734 gmax = 2*gmax + 1 2735 call ngfft_seq(work_ngfft, gmax) 2736 !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3) 2737 ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6))) 2738 2739 ! Allow PW-arrays dimensioned with mpw 2740 ABI_MALLOC(kg_k, (3, mpw)) 2741 ABI_MALLOC(kg_kq, (3, mpw)) 2742 2743 ! Spherical Harmonics for useylm == 1. 2744 ABI_MALLOC(ylm_k, (mpw, psps%mpsang*psps%mpsang*psps%useylm)) 2745 ABI_MALLOC(ylm_kq, (mpw, psps%mpsang*psps%mpsang*psps%useylm)) 2746 ABI_MALLOC(ylmgr_kq, (mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 2747 2748 usecprj = 0 2749 ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj)) 2750 2751 ! Prepare call to getgh1c 2752 usevnl = 0 2753 optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C> 2754 optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C> 2755 opt_gvnlx1 = 0 ! gvnlx1 is output 2756 ABI_MALLOC(gvnlx1, (2, usevnl)) 2757 ABI_MALLOC(grad_berry, (2, nspinor*(berryopt0/4))) 2758 2759 ! This part is taken from dfpt_vtorho 2760 !==== Initialize most of the Hamiltonian (and derivative) ==== 2761 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 2762 !2) Perform the setup needed for the non-local factors: 2763 ! Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 2764 ! PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 2765 2766 call init_hamiltonian(gs_hamkq, psps, pawtab, nspinor, nsppol, nspden, natom, & 2767 dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg, & 2768 comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, & 2769 usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option) 2770 2771 ! Allocate vlocal. Note nvloc 2772 ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data) 2773 ABI_MALLOC(vlocal, (n4, n5, n6, gs_hamkq%nvloc)) 2774 vlocal = huge(one) 2775 2776 ! Allocate work space arrays. 2777 ABI_MALLOC(displ_cart_qbz, (2, 3, cryst%natom, natom3)) 2778 ABI_MALLOC(displ_cart_qibz, (2, 3, cryst%natom, natom3)) 2779 ABI_MALLOC(displ_red_qbz, (2, 3, cryst%natom, natom3)) 2780 ABI_MALLOC(displ_red_qibz, (2, 3, cryst%natom, natom3)) 2781 ABI_MALLOC(pheigvec_qbz, (2, 3, cryst%natom, 3*cryst%natom)) 2782 ABI_MALLOC(pheigvec_qibz, (2, 3, cryst%natom, 3*cryst%natom)) 2783 2784 ABI_CALLOC(dummy_vtrial, (nfftf, nspden)) 2785 2786 ! Create ddkop object to compute group velocities (if needed) 2787 ddkop = ddkop_new(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft) 2788 2789 ! Open GSTORE file, and read table used for restarting. 2790 NCF_CHECK(nctk_open_modify(root_ncid, gstore%path, gstore%comm)) 2791 2792 ABI_MALLOC(done_qbz_spin, (gstore%nqbz, nsppol)) 2793 NCF_CHECK(nf90_get_var(root_ncid, nctk_idname(root_ncid, "gstore_done_qbz_spin"), done_qbz_spin)) 2794 gstore%wfk0_path = wfk0_path 2795 if (my_rank == master) then 2796 NCF_CHECK(nf90_put_var(root_ncid, root_vid("gstore_wfk0_path"), trim(gstore%wfk0_path))) 2797 end if 2798 2799 if (my_rank == master) call gstore%print(std_out) 2800 2801 ndone = count(done_qbz_spin == 1) 2802 2803 ! NB: Write phonon data here as we are not guaranteed to have all the IBZ q-points 2804 ! inside the loop over my_iq if filtering has been used. 2805 ! TODO: Write phdata only if eph_task == -11 2806 if (ndone == 0) then 2807 call wrtout(std_out, " Computing phonon frequencies and displacements in the IBZ") 2808 call cwtime(cpu, wall, gflops, "start") 2809 2810 call xmpi_split_block(gstore%nqibz, gstore%comm, my_nqibz, my_iqibz_inds) 2811 ABI_MALLOC(buf_wqnu, (natom3, my_nqibz)) 2812 ABI_MALLOC(buf_eigvec_cart, (2, 3, natom, natom3, my_nqibz)) 2813 2814 do ii=1,my_nqibz 2815 iq_ibz = my_iqibz_inds(ii) 2816 !call ifc%fourq(cryst, gstore%qibz(:, iq_ibz), buf_wqnu(:,ii), buf_displ_cart(:,:,:,:,ii)) 2817 call ifc%fourq(cryst, gstore%qibz(:, iq_ibz), buf_wqnu(:,ii), displ_cart_qibz, & 2818 out_eigvec=buf_eigvec_cart(:,:,:,:,ii)) 2819 end do 2820 if (nproc > 1 .and. gstore%nqibz >= nproc) then 2821 NCF_CHECK(nctk_set_collective(root_ncid, root_vid("phfreqs_ibz"))) 2822 NCF_CHECK(nctk_set_collective(root_ncid, root_vid("pheigvec_cart_ibz"))) 2823 end if 2824 2825 if (my_nqibz > 0) then 2826 iq_start = my_iqibz_inds(1) 2827 ncerr = nf90_put_var(root_ncid, root_vid("phfreqs_ibz"), buf_wqnu, start=[1, iq_start], count=[natom3, my_nqibz]) 2828 NCF_CHECK(ncerr) 2829 ncerr = nf90_put_var(root_ncid, root_vid("pheigvec_cart_ibz"), buf_eigvec_cart, & 2830 start=[1,1,1,1,iq_start], count=[2, 3, natom, natom3, my_nqibz]) 2831 NCF_CHECK(ncerr) 2832 end if 2833 2834 ABI_FREE(my_iqibz_inds) 2835 ABI_FREE(buf_wqnu) 2836 ABI_FREE(buf_eigvec_cart) 2837 call cwtime_report(" phonon computation + output", cpu, wall, gflops) 2838 else 2839 call wrtout(std_out, & 2840 sjoin(" Restarting GSTORE calculation. Found: ", itoa(ndone), " (qpt, spin) entries already computed")) 2841 end if 2842 2843 if (gstore%with_vk /= 0 .and. ndone == 0) then 2844 call wrtout(std_out, " Computing and writing velocity operator matrix elements in the IBZ") 2845 call wrtout(std_out, " Note that not all the k-points in the IBZ are computed when kfilter is activated!") 2846 call cwtime(cpu, wall, gflops, "start") 2847 2848 ! On disk, we have: 2849 ! 2850 ! nctkarr_t("vk_cart_ibz", "dp", "three, nb, gstore_nkibz")) 2851 ! nctkarr_t("vkmat_cart_ibz", "dp", "two, three, nb, nb, gstore_nkibz"))) 2852 2853 ABI_MALLOC(cgwork, (2, mpw*wfd%nspinor)) 2854 2855 do my_is=1,gstore%my_nspins 2856 spin = gstore%my_spins(my_is) 2857 gqk => gstore%gqk(my_is) 2858 2859 if (gstore%with_vk == 1) then 2860 ABI_CALLOC(vk_cart_ibz, (3, gqk%nb, gstore%nkibz)) 2861 else 2862 ABI_ERROR("with_vk 2") 2863 end if 2864 2865 cnt = 0 2866 do my_ik=1,gqk%my_nk 2867 ! The k-point and the symmetries relating the BZ k-point to the IBZ. 2868 kk_bz = gqk%my_kpts(:, my_ik) 2869 weight_k = gqk%my_wtk(my_ik) 2870 2871 ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik) 2872 trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5,my_ik) 2873 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 2874 if (.not. isirr_k) cycle 2875 2876 ! parallelize inside (q, pert) so that only one proc in the 3D grid 2877 ! computes vk for this kpt in the BZ and we can use xmpi_sum_master. 2878 cnt = cnt + 1 2879 if (gqk%qpt_pert_comm%skip(cnt)) cycle 2880 2881 kk_ibz = ebands%kptns(:,ik_ibz) 2882 npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz) 2883 call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk_ibz, istwf_k, npw_k, wfd%kdata(ik_ibz)%kg_k) 2884 2885 select case (gstore%with_vk) 2886 case (1) 2887 do ib_k=1,gqk%nb 2888 band_k = ib_k + gqk%bstart - 1 2889 call wfd%copy_cg(band_k, ik_ibz, spin, cgwork) 2890 vk = ddkop%get_vdiag(ebands%eig(band_k, ik_ibz, spin), istwf_k, npw_k, wfd%nspinor, cgwork, cwaveprj0) 2891 vk_cart_ibz(:, ib_k, ik_ibz) = vk 2892 end do 2893 2894 case (2) 2895 ABI_ERROR("with_vk 2") 2896 !do ib_k=1,nband_k 2897 ! band_k = ib_k + bstart_k - 1 2898 !end do 2899 end select 2900 end do ! my_ik 2901 2902 call xmpi_sum_master(vk_cart_ibz, master, gqk%grid_comm%value, ierr) 2903 if (gqk%grid_comm%me == master) then 2904 NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid)) 2905 NCF_CHECK(nf90_put_var(spin_ncid, spin_vid("vk_cart_ibz"), vk_cart_ibz)) 2906 end if 2907 ABI_SFREE(vk_cart_ibz) 2908 end do ! my_is 2909 2910 ABI_FREE(cgwork) 2911 call cwtime_report(sjoin(" Computation of v_k group velocities with with_vk:", itoa(gstore%with_vk)), cpu, wall, gflops) 2912 end if 2913 2914 ! Loop over my spins. 2915 do my_is=1,gstore%my_nspins 2916 spin = gstore%my_spins(my_is) 2917 gqk => gstore%gqk(my_is) 2918 my_npert = gqk%my_npert 2919 NCF_CHECK(nf90_inq_ncid(root_ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid)) 2920 2921 ! Allocate workspace for wavefunctions using mpw and nb 2922 nb = gqk%nb 2923 ABI_MALLOC(bras_kq, (2, mpw*nspinor, nb)) 2924 ABI_MALLOC(kets_k, (2, mpw*nspinor, nb)) 2925 ABI_MALLOC(h1kets_kq, (2, mpw*nspinor, nb)) 2926 ABI_MALLOC(gkk_atm, (2, nb, nb, natom3)) 2927 ABI_MALLOC(gkq_nu, (2, nb, nb, natom3)) 2928 ABI_MALLOC(iq_buf, (2, qbuf_size)) 2929 2930 ! Inside the loops we compute gkq_nu(2, nb, nb, natom3) 2931 ABI_MALLOC_OR_DIE(my_gbuf, (gqk%cplex, nb, nb, natom3, gqk%my_nk, qbuf_size), ierr) 2932 2933 ! Loop over my set of q-points 2934 prev_iqbz = -1 2935 do my_iq=1,gqk%my_nq 2936 print_time = my_rank == 0 .and. (my_iq <= LOG_MODQ .or. mod(my_iq, LOG_MODQ) == 0) 2937 if (print_time) call cwtime(cpu_q, wall_q, gflops_q, "start") 2938 iq_bz = gqk%my_q2bz(my_iq) 2939 if (done_qbz_spin(iq_bz, spin) == 1) cycle 2940 2941 call gqk%myqpt(my_iq, gstore, weight_q, qq_bz) 2942 2943 iq_ibz = gqk%my_q2ibz(1, my_iq); isym_q = gqk%my_q2ibz(2, my_iq) 2944 trev_q = gqk%my_q2ibz(6, my_iq); g0_q = gqk%my_q2ibz(3:5,my_iq) 2945 ! Don't test if umklapp == 0 because we use the periodic gauge: 2946 ! 2947 ! phfreq(q+G) = phfreq(q) and eigvec(q) = eigvec(q+G) 2948 ! 2949 !isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0)) 2950 isirr_q = (isym_q == 1 .and. trev_q == 0) 2951 qq_ibz = gstore%qibz(:, iq_ibz) 2952 2953 nskip_tetra_kq = 0 2954 iqbuf_cnt = 1 + mod(my_iq - 1, qbuf_size) 2955 iq_buf(:, iqbuf_cnt) = [my_iq, iq_bz] 2956 2957 if (iq_ibz /= prev_iqbz) then 2958 ! Get phonon frequencies and eigenvectors for the corresponding q-point in the IBZ. 2959 call ifc%fourq(cryst, qq_ibz, phfrq, displ_cart_qibz, out_displ_red=displ_red_qibz, out_eigvec=pheigvec_qibz) 2960 prev_iqbz = iq_ibz 2961 end if 2962 2963 if (isirr_q) then 2964 displ_cart_qbz = displ_cart_qibz; displ_red_qbz = displ_red_qibz; pheigvec_qbz = pheigvec_qibz 2965 else 2966 ! Rotate phonon eigenvectors from q_ibz to q_bz. 2967 ! This part is needed to enforce the gauge in the ph eigenvectors, including e(-q) = e(q)^* 2968 call pheigvec_rotate(cryst, qq_ibz, isym_q, trev_q, pheigvec_qibz, pheigvec_qbz, displ_cart_qbz, & 2969 displ_red_qbz=displ_red_qbz) 2970 end if 2971 2972 !call ifc%fourq(cryst, qq_bz, phfrq, displ_cart_qbz, out_displ_red=displ_red_qbz, out_eigvec=pheigvec_qbz)) 2973 ! Use Fourier interpolation of DFPT potentials to get my_npert potentials. 2974 !cplex = 2 2975 !ABI_MALLOC(v1scf, (cplex, nfft, nspden, dvdb%my_npert)) 2976 !call dvdb%ftinterp_qpt(qq_bz, nfftf, ngfftf, v1scf, dvdb%comm_rpt) 2977 2978 ! Version with qcache. 2979 call dvdb%get_ftqbz(cryst, qq_bz, qq_ibz, gqk%my_q2ibz(:, my_iq), cplex, nfftf, ngfftf, v1scf, & 2980 gqk%pert_comm%value) 2981 2982 ! Allocate vlocal1 with correct cplex. Note nvloc and my_npert. 2983 ABI_MALLOC(vlocal1, (cplex*n4, n5, n6, gs_hamkq%nvloc, my_npert)) 2984 2985 ! Set up local potential vlocal1 with proper dimensioning from vtrial1 taking into account the spin. 2986 do my_ip=1,my_npert 2987 call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, gs_hamkq%nvloc,& 2988 pawfgr, mpi_enreg, dummy_vtrial, v1scf(:,:,:,my_ip), vlocal, vlocal1(:,:,:,:,my_ip)) 2989 end do 2990 2991 ! Continue to initialize the GS Hamiltonian 2992 call gs_hamkq%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.) 2993 2994 ! Loop over my k-points 2995 do my_ik=1,gqk%my_nk 2996 ! Set entry to zero. Important as there are cycle instructions inside these loops 2997 ! and we don't want to write random numbers to disk. 2998 my_gbuf(:,:,:,:, my_ik, iqbuf_cnt) = zero 2999 3000 ! The k-point and the symmetries relating the BZ k-point to the IBZ. 3001 kk_bz = gqk%my_kpts(:, my_ik) 3002 weight_k = gqk%my_wtk(my_ik) 3003 3004 ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik) 3005 trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5,my_ik) 3006 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 3007 kk_ibz = ebands%kptns(:,ik_ibz) 3008 3009 ! Number of bands crossing the Fermi level at k 3010 !bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz) 3011 bstart_k = gqk%bstart; nband_k = gqk%nb 3012 3013 ! ============================================ 3014 ! Find symmetrical image of k+q in in the kIBZ 3015 ! ============================================ 3016 kq_bz = kk_bz + qq_bz 3017 3018 if (kpts_map("symrel", ebands_timrev, cryst, gstore%krank_ibz, 1, kq_bz, indkk_kq) /= 0) then 3019 write(msg, '(3a)' ) & 3020 "Cannot find k+q in kmesh", ch10, 'Action: check your WFK file and the (k, q) point input variables.' 3021 ABI_ERROR(msg) 3022 end if 3023 3024 ikq_ibz = indkk_kq(1, 1); isym_kq = indkk_kq(2, 1) 3025 trev_kq = indkk_kq(6, 1); g0_kq = indkk_kq(3:5, 1) 3026 isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0)) 3027 kq_ibz = ebands%kptns(:, ikq_ibz) 3028 3029 ! If we have used the KERANGE trick, we may have k or k+q points with just one G component set to zero 3030 ! so we skip this transition immediately. This should happen only if fsewin > sigma_erange. 3031 if (wfd%npwarr(ik_ibz) == 1 .or. wfd%npwarr(ikq_ibz) == 1) cycle 3032 3033 if (gstore%kfilter == "fs_tetra") then 3034 ! Check tetra delta(e_{k+q}) and cycle if all the weights at k+q are zero. 3035 if (all(abs(gstore%delta_ef_kibz_spin(:, ikq_ibz, spin)) == zero)) then 3036 nskip_tetra_kq = nskip_tetra_kq + 1; cycle 3037 end if 3038 end if 3039 3040 ! Number of bands crossing the Fermi level at k+q 3041 !bstart_kq = fs%bstart_cnt_ibz(1, ikq_ibz); nband_kq = fs%bstart_cnt_ibz(2, ikq_ibz) 3042 bstart_kq = gqk%bstart; nband_kq = gqk%nb 3043 ABI_CHECK(nband_k <= nb .and. nband_kq <= nb, "wrong nband") 3044 3045 ! Get npw_k, kg_k and symmetrize wavefunctions from IBZ (if needed). 3046 call wfd%sym_ug_kg(ecut, kk_bz, kk_ibz, bstart_k, nband_k, spin, mpw, gqk%my_k2ibz(:, my_ik), cryst, & 3047 work_ngfft, work, istwf_k, npw_k, kg_k, kets_k) 3048 3049 ! Get npw_kq, kg_kq and symmetrize wavefunctions from IBZ (if needed). 3050 call wfd%sym_ug_kg(ecut, kq_bz, kq_ibz, bstart_kq, nband_kq, spin, mpw, indkk_kq(:,1), cryst, & 3051 work_ngfft, work, istwf_kq, npw_kq, kg_kq, bras_kq) 3052 3053 ! if PAW, one has to solve a generalized eigenproblem 3054 ! Be careful here because I will need sij_opt==-1 3055 gen_eigenpb = psps%usepaw == 1; sij_opt = 0; if (gen_eigenpb) sij_opt = 1 3056 ABI_MALLOC(gs1c, (2, npw_kq*nspinor*((sij_opt+1)/2))) 3057 3058 ! Set up the spherical harmonics (Ylm) at k and k+q. See also dfpt_looppert 3059 !if (psps%useylm == 1) then 3060 ! optder = 0; if (useylmgr == 1) optder = 1 3061 ! call initylmg(cryst%gprimd, kg_k, kk_bz, mkmem1, mpi_enreg, psps%mpsang, mpw, nband, mkmem1,& 3062 ! [npw_k], dtset%nsppol, optder, cryst%rprimd, ylm_k, ylmgr) 3063 ! call initylmg(cryst%gprimd, kg_kq, kq_bz, mkmem1, mpi_enreg, psps%mpsang, mpw, nband, mkmem1,& 3064 ! [npw_kq], dtset%nsppol, optder, cryst%rprimd, ylm_kq, ylmgr_kq) 3065 !end if 3066 3067 ! Compute k+G vectors 3068 nkpg = 3 * dtset%nloalg(3) 3069 ABI_MALLOC(kpg_k, (npw_k, nkpg)) 3070 if (nkpg > 0) call mkkpg(kg_k, kpg_k, kk_bz, nkpg, npw_k) 3071 3072 ! Compute nonlocal form factors ffnlk at (k+G) 3073 ABI_MALLOC(ffnlk, (npw_k, 1, psps%lmnmax, psps%ntypat)) 3074 call mkffnl_objs(cryst, psps, 1, ffnlk, ider0, idir0, kg_k, kpg_k, kk_bz, nkpg, npw_k, ylm_k, ylmgr_dum, & 3075 comm=gqk%pert_comm%value) !, request=ffnlk_request) 3076 3077 ! Compute k+q+G vectors 3078 nkpg1 = 3 * dtset%nloalg(3) 3079 ABI_MALLOC(kpg1_k, (npw_kq, nkpg1)) 3080 if (nkpg1 > 0) call mkkpg(kg_kq, kpg1_k, kq_bz, nkpg1, npw_kq) 3081 3082 ! Compute nonlocal form factors ffnl1 at (k+q+G) 3083 ABI_MALLOC(ffnl1, (npw_kq, 1, psps%lmnmax, psps%ntypat)) 3084 call mkffnl_objs(cryst, psps, 1, ffnl1, ider0, idir0, kg_kq, kpg1_k, kq_bz, nkpg1, npw_kq, ylm_kq, ylmgr_kq, & 3085 comm=gqk%pert_comm%value) ! request=ffnl1_request) 3086 3087 ! Loop over my atomic perturbations and compute gkk_atm. 3088 gkk_atm = zero 3089 do my_ip=1,my_npert 3090 idir = dvdb%my_pinfo(1, my_ip); ipert = dvdb%my_pinfo(2, my_ip); ipc = dvdb%my_pinfo(3, my_ip) 3091 3092 ! Prepare application of the NL part. 3093 call init_rf_hamiltonian(cplex, gs_hamkq, ipert, rf_hamkq, has_e1kbsc=.true.) 3094 3095 call rf_hamkq%load_spin(spin, vlocal1=vlocal1(:,:,:,:,my_ip), with_nonlocal=.true.) 3096 3097 ! This call is not optimal because there are quantities in out that do not depend on idir,ipert 3098 call getgh1c_setup(gs_hamkq, rf_hamkq, dtset, psps, kk_bz, kq_bz, idir, ipert, & ! In 3099 cryst%natom, cryst%rmet, cryst%gprimd, cryst%gmet, istwf_k, & ! In 3100 npw_k, npw_kq, useylmgr1, kg_k, ylm_k, kg_kq, ylm_kq, ylmgr_kq, & ! In 3101 dkinpw, nkpg, nkpg1, kpg_k, kpg1_k, kinpw1, ffnlk, ffnl1, ph3d, ph3d1, & ! Out 3102 reuse_kpg_k=1, reuse_kpg1_k=1, reuse_ffnlk=1, reuse_ffnl1=1) ! Reuse some arrays 3103 3104 ! Calculate dvscf * psi_k, results stored in h1kets_kq on the k+q sphere. 3105 ! Compute H(1) applied to GS wavefunction Psi(0) 3106 do ib_k=1,nband_k 3107 band_k = ib_k + bstart_k - 1 3108 eig0nk = ebands%eig(band_k, ik_ibz, spin) 3109 ! Use scissor shift on 0-order eigenvalue 3110 eshift = eig0nk - dtset%dfpt_sciss 3111 3112 call getgh1c(berryopt0, kets_k(:,:,ib_k), cwaveprj0, h1kets_kq(:,:,ib_k), & 3113 grad_berry, gs1c, gs_hamkq, gvnlx1, idir, ipert, eshift, mpi_enreg, optlocal, & 3114 optnl, opt_gvnlx1, rf_hamkq, sij_opt, tim_getgh1c, usevnl) 3115 end do 3116 3117 call rf_hamkq%free() 3118 ABI_FREE(kinpw1) 3119 ABI_FREE(dkinpw) 3120 ABI_FREE(ph3d) 3121 ABI_SFREE(ph3d1) 3122 3123 ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation. 3124 ! No need to handle istwf_kq because it's always 1. 3125 do ib_k=1,nband_k 3126 do ib_kq=1,nband_kq 3127 gkk_atm(:, ib_kq, ib_k, ipc) = cg_zdotc(npw_kq*nspinor, bras_kq(1,1,ib_kq), h1kets_kq(1,1,ib_k)) 3128 end do 3129 end do 3130 3131 end do ! my_ip 3132 3133 ABI_FREE(gs1c) 3134 ABI_FREE(ffnlk) 3135 ABI_FREE(ffnl1) 3136 ABI_FREE(kpg1_k) 3137 ABI_FREE(kpg_k) 3138 3139 ! Collect gkk_atm inside pert_comm so that all procs can operate on the data. 3140 if (gqk%pert_comm%nproc > 1) call xmpi_sum(gkk_atm, gqk%pert_comm%value, ierr) 3141 3142 ! Get g in the phonon representation. 3143 call ephtk_gkknu_from_atm(nb, nb, 1, natom, gkk_atm, phfrq, displ_red_qbz, gkq_nu) 3144 3145 ! Save e-ph matrix elements in the buffer. 3146 select case (gqk%cplex) 3147 case (1) 3148 my_gbuf(1,:,:,:, my_ik, iqbuf_cnt) = gkq_nu(1,:,:,:) ** 2 + gkq_nu(2,:,:,:) ** 2 3149 case (2) 3150 my_gbuf(:,:,:,:, my_ik, iqbuf_cnt) = gkq_nu 3151 end select 3152 end do ! my_ik 3153 3154 ABI_FREE(v1scf) 3155 ABI_FREE(vlocal1) 3156 3157 ! Dump buffers 3158 if (iqbuf_cnt == qbuf_size) call dump_data() 3159 3160 if (print_time) then 3161 write(msg,'(2(a,i0),a)')" My q-point [", my_iq, "/", gqk%my_nq, "]" 3162 call cwtime_report(msg, cpu_q, wall_q, gflops_q); if (my_iq == LOG_MODQ) call wrtout(std_out, "...", do_flush=.True.) 3163 end if 3164 !if (my_rank == master) then 3165 ! ! Print cache stats. 3166 ! !call dvdb%ft_qcache%report_stats() 3167 !end if 3168 end do ! my_iq 3169 3170 ! Dump the remainder. 3171 if (iqbuf_cnt /= 0) call dump_data() 3172 3173 ABI_FREE(iq_buf) 3174 ABI_FREE(my_gbuf) 3175 ABI_FREE(bras_kq) 3176 ABI_FREE(kets_k) 3177 ABI_FREE(h1kets_kq) 3178 ABI_FREE(gkk_atm) 3179 ABI_FREE(gkq_nu) 3180 end do ! my_is 3181 3182 call cwtime_report(" GSTORE computation done", cpu_all, wall_all, gflops_all, pre_str=ch10, end_str=ch10) !, comm=gstore%comm) 3183 call gstore%print(std_out, header="GSTORE at the end of gstore%compute") 3184 3185 !call xmpi_barrier(gstore%comm) 3186 NCF_CHECK(nf90_close(root_ncid)) 3187 3188 ! Free memory 3189 ABI_FREE(gvnlx1) 3190 ABI_FREE(grad_berry) 3191 ABI_FREE(dummy_vtrial) 3192 ABI_FREE(work) 3193 ABI_FREE(ph1d) 3194 ABI_FREE(vlocal) 3195 ABI_FREE(kg_k) 3196 ABI_FREE(kg_kq) 3197 ABI_FREE(ylm_k) 3198 ABI_FREE(ylm_kq) 3199 ABI_FREE(ylmgr_kq) 3200 ABI_FREE(displ_cart_qbz) 3201 ABI_FREE(displ_cart_qibz) 3202 ABI_FREE(displ_red_qbz) 3203 ABI_FREE(displ_red_qibz) 3204 ABI_FREE(pheigvec_qbz) 3205 ABI_FREE(pheigvec_qibz) 3206 ABI_FREE(done_qbz_spin) 3207 3208 call pawcprj_free(cwaveprj0) 3209 ABI_FREE(cwaveprj0) 3210 call ddkop%free() 3211 call gs_hamkq%free() 3212 call wfd%free() 3213 3214 contains 3215 3216 subroutine dump_data() 3217 3218 ! This function is called inside the double loop over (my_is, my_iq) or when we exit 3219 ! from the my_iq loop to dump the remainder that is still in the q-buffer, 3220 ! All the MPI procs in the (kpt_comm x pert_comm) grid shall call this contained routine 3221 ! as we have side-effects i.e. iqbuf_cnt set to 0. 3222 3223 ! On disk we have the global arrays: 3224 ! 3225 ! nctkarr_t("gvals", "dp", "gstore_cplex, nb, nb, natom3, glob_nk, glob_nq") 3226 ! 3227 ! while the local MPI buffers are dimensioned as follows: 3228 ! 3229 ! my_gbuf(2, nb, nb, natom3, gqk%my_nk, qbuf_size) 3230 3231 ! If parallelism over pertubation is activated, only the procs treating the first perturbation 3232 ! i.e. the procs treating different k-points for this q are involved in IO 3233 ! as all the local buffers store results for all natom3 pertubations. 3234 3235 integer :: ii, iq_bz, iq_glob, my_iq 3236 3237 if (gqk%coords_qkp(3) /= 0) goto 10 ! Yes, I'm very proud of this GOTO. 3238 3239 !iq_buf(:, iqbuf_cnt) = [my_iq, iq_bz] 3240 my_iq = iq_buf(1, 1) 3241 iq_glob = my_iq + gqk%my_qstart - 1 3242 3243 ! NB: this is an individual IO operation 3244 ncerr = nf90_put_var(spin_ncid, spin_vid("gvals"), my_gbuf, & 3245 start=[1, 1, 1, 1, gqk%my_kstart, iq_glob], & 3246 count=[gqk%cplex, gqk%nb, gqk%nb, gqk%natom3, gqk%my_nk, iqbuf_cnt] & 3247 ) 3248 NCF_CHECK(ncerr) 3249 3250 ! Only one proc sets the entry in done_qbz_spin to 1 for all the q-points in the buffer. 3251 if (all(gqk%coords_qkp(2:3) == [0, 0])) then 3252 do ii=1,iqbuf_cnt 3253 iq_bz = iq_buf(2, ii) 3254 NCF_CHECK(nf90_put_var(root_ncid, root_vid("gstore_done_qbz_spin"), 1, start=[iq_bz, spin])) 3255 end do 3256 end if 3257 3258 ! Zero the counter before returning 3259 10 iqbuf_cnt = 0 3260 3261 end subroutine dump_data 3262 3263 integer function root_vid(vname) 3264 character(len=*),intent(in) :: vname 3265 root_vid = nctk_idname(root_ncid, vname) 3266 end function root_vid 3267 3268 integer function spin_vid(vname) 3269 character(len=*),intent(in) :: vname 3270 spin_vid = nctk_idname(spin_ncid, vname) 3271 end function spin_vid 3272 3273 end subroutine gstore_compute
m_gstore/gstore_distribute_spins [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_distribute_spins
FUNCTION
INPUTS
OUTPUT
SOURCE
956 subroutine gstore_distribute_spins(gstore, mband, gstore_brange, nproc_spin, comm_spin, comm) 957 958 !Arguments ------------------------------------ 959 !scalars 960 class(gstore_t),target,intent(inout) :: gstore 961 integer,intent(in) :: mband, comm, gstore_brange(2, gstore%nsppol) 962 integer,intent(out) :: nproc_spin(gstore%nsppol), comm_spin(gstore%nsppol) 963 964 !Local variables------------------------------- 965 !scalars 966 integer :: spin, my_rank, ierr, color, nsppol, all_nproc 967 !arrays 968 integer :: buff_spin(2) 969 !---------------------------------------------------------------------- 970 971 all_nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 972 nsppol = gstore%nsppol 973 974 gstore%my_nspins = 0 975 ABI_MALLOC(gstore%brange_spin, (2, nsppol)) 976 977 do spin=1,nsppol 978 ! NB: If MPI_UNDEFINED is passed as the colour value, the subgroup in which the calling 979 ! MPI process will be placed is MPI_COMM_NULL 980 color = 1 981 if (nsppol == 2 .and. all_nproc > 1) then 982 color = xmpi_undefined 983 if (spin == 1 .and. my_rank <= (all_nproc - 1) / 2) color = 1 984 if (spin == 2 .and. my_rank > (all_nproc - 1) / 2) color = 1 985 end if 986 987 call xmpi_comm_split(comm, color, my_rank, comm_spin(spin), ierr) 988 if (comm_spin(spin) /= xmpi_undefined) then 989 gstore%my_nspins = gstore%my_nspins + 1 990 buff_spin(gstore%my_nspins) = spin 991 end if 992 993 nproc_spin(spin) = xmpi_comm_size(comm_spin(spin)) 994 995 gstore%brange_spin(:, spin) = [1, mband] 996 if (all(gstore_brange /= 0)) gstore%brange_spin(:, spin) = gstore_brange(:, spin) 997 998 ABI_CHECK_IRANGE(gstore%brange_spin(1, spin), 1, mband, "gstore_brange(1, spin)") 999 ABI_CHECK_IRANGE(gstore%brange_spin(2, spin), 1, mband, "gstore_brange(2, spin)") 1000 ABI_CHECK(gstore%brange_spin(1, spin) <= gstore%brange_spin(2, spin), "brange_spin(1, spin) <= brange_spin(2, spin)") 1001 end do 1002 1003 ABI_MALLOC(gstore%my_spins, (gstore%my_nspins)) 1004 gstore%my_spins = buff_spin(1:gstore%my_nspins) 1005 ABI_MALLOC(gstore%gqk, (gstore%my_nspins)) 1006 1007 end subroutine gstore_distribute_spins
m_gstore/gstore_fill_bks_mask [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_fill_bks_mask
FUNCTION
Fills the bks_mask array defining the set of states that should be read from file by this MPI rank when computing the e-ph matrix elements.
INPUTS
OUTPUT
SOURCE
1786 subroutine gstore_fill_bks_mask(gstore, mband, nkibz, nsppol, bks_mask) 1787 1788 !Arguments ------------------------------------ 1789 class(gstore_t),target,intent(inout) :: gstore 1790 integer,intent(in) :: mband, nkibz, nsppol 1791 logical,intent(out) :: bks_mask(mband, nkibz, nsppol) 1792 1793 !Local variables------------------------------- 1794 !scalars 1795 integer :: my_is, my_ik, my_iq, spin, ik_ibz, iqk_ibz, ebands_timrev 1796 real(dp) :: weight_q, cpu, wall, gflops 1797 type(gqk_t),pointer :: gqk 1798 type(crystal_t),pointer :: cryst 1799 !arrays 1800 integer,allocatable :: map_kq(:,:) 1801 real(dp) :: qpt(3) 1802 1803 !---------------------------------------------------------------------- 1804 1805 call cwtime(cpu, wall, gflops, "start") 1806 1807 bks_mask = .False. 1808 cryst => gstore%cryst 1809 1810 ebands_timrev = kpts_timrev_from_kptopt(gstore%ebands%kptopt) 1811 1812 do my_is=1,gstore%my_nspins 1813 gqk => gstore%gqk(my_is) 1814 spin = gstore%my_spins(my_is) 1815 1816 ! We need the image of this k in the IBZ. 1817 do my_ik=1,gqk%my_nk 1818 ik_ibz = gqk%my_k2ibz(1, my_ik) 1819 bks_mask(gqk%bstart:gqk%bstart + gqk%nb - 1, ik_ibz, spin) = .True. 1820 end do 1821 1822 ! as well as the image of k+q in the IBZ. 1823 ABI_MALLOC(map_kq, (6, gqk%my_nk)) 1824 1825 do my_iq=1,gqk%my_nq 1826 call gqk%myqpt(my_iq, gstore, weight_q, qpt) 1827 1828 if (kpts_map("symrel", ebands_timrev, cryst, gstore%krank_ibz, gqk%my_nk, gqk%my_kpts, map_kq, qpt=qpt) /= 0) then 1829 ABI_ERROR(sjoin("Cannot map k+q to IBZ with qpt:", ktoa(qpt))) 1830 end if 1831 1832 do my_ik=1,gqk%my_nk 1833 iqk_ibz = map_kq(1, my_ik) 1834 bks_mask(gqk%bstart:gqk%bstop, iqk_ibz, spin) = .True. 1835 end do 1836 end do 1837 1838 ABI_FREE(map_kq) 1839 end do 1840 1841 call cwtime_report(" gstore_fill_bks_mask", cpu, wall, gflops) 1842 1843 end subroutine gstore_fill_bks_mask
m_gstore/gstore_filter_erange__ [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_filter_erange__
FUNCTION
Filter k-points and q-points according to an energy range.
INPUTS
OUTPUT
SOURCE
1526 subroutine gstore_filter_erange__(gstore, qbz, qbz2ibz, qibz2bz, kbz, kibz, kbz2ibz, kibz2bz, & 1527 select_qbz_spin, select_kbz_spin) 1528 1529 !Arguments ------------------------------------ 1530 !scalars 1531 class(gstore_t),target,intent(inout) :: gstore 1532 real(dp),intent(in) :: qbz(3, gstore%nqbz) 1533 real(dp),target,intent(in) :: kbz(3, gstore%nqbz), kibz(3, gstore%nkibz) 1534 integer,intent(in) :: qbz2ibz(6,gstore%nqbz), qibz2bz(gstore%nqibz) 1535 integer,intent(in) :: kbz2ibz(6,gstore%nkbz), kibz2bz(gstore%nkibz) 1536 integer,intent(out) :: select_qbz_spin(gstore%nqbz, gstore%nsppol) 1537 integer,intent(out) :: select_kbz_spin(gstore%nkbz, gstore%nsppol) 1538 1539 !Local variables------------------------------- 1540 !scalars 1541 integer,parameter :: tetra_opt0 = 0 1542 integer :: nsppol, cnt, spin, ii, all_nproc, my_rank, comm, ebands_timrev, band ! ierr, 1543 integer :: ik_bz, ik_ibz, iflag, nk_in_star, gap_err 1544 logical :: assume_gap 1545 real(dp) :: ee, abs_erange1, abs_erange2, vmax, cmin 1546 type(ebands_t),pointer :: ebands 1547 type(crystal_t),pointer :: cryst 1548 type(gaps_t) :: gaps 1549 !arrays 1550 integer,allocatable :: kstar_bz_inds(:) 1551 !---------------------------------------------------------------------- 1552 1553 comm = gstore%comm 1554 all_nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1555 1556 ! filter k and k + q points according to erange and define gstore%brange_spin automatically. 1557 ! NB: here we recompute brange_spin 1558 call wrtout(std_out, sjoin(" Filtering k-points using gstore_erange:", & 1559 ltoa(reshape(gstore%erange_spin, [2 * gstore%nsppol]) * Ha_eV), "(eV)")) 1560 1561 cryst => gstore%cryst; ebands => gstore%ebands; nsppol = gstore%nsppol 1562 1563 assume_gap = .not. all(gstore%erange_spin < zero) 1564 gaps = ebands_get_gaps(ebands, gap_err) 1565 if (assume_gap) call gaps%print(unit=std_out) !, header=msg) 1566 1567 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 1568 1569 select_kbz_spin = 0; cnt = 0 1570 1571 do spin=1,nsppol 1572 gstore%brange_spin(:, spin) = [huge(1), -huge(1)] 1573 abs_erange1 = abs(gstore%erange_spin(1, spin)) 1574 abs_erange2 = abs(gstore%erange_spin(2, spin)) 1575 1576 if (assume_gap) then 1577 ! Get CBM and VBM with some tolerance 1578 vmax = gaps%vb_max(spin) + tol2 * eV_Ha 1579 cmin = gaps%cb_min(spin) - tol2 * eV_Ha 1580 else 1581 vmax = ebands%fermie 1582 cmin = ebands%fermie 1583 end if 1584 1585 do band=1, ebands%mband 1586 do ik_ibz=1,gstore%nkibz 1587 !cnt = cnt + 1; if (mod(cnt, all_nproc) /= my_rank) cycle ! MPI parallelism inside comm 1588 1589 ! Use iflag to filter k-points. 1590 ee = ebands%eig(band, ik_ibz, spin) 1591 iflag = 0 1592 1593 if (abs_erange1 > zero) then 1594 if (ee <= vmax .and. vmax - ee <= abs_erange1) then 1595 iflag = 1 !; write(std_out, *), "Adding valence band", band, " with ee [eV]: ", ee * Ha_eV 1596 end if 1597 end if 1598 if (abs_erange2 > zero) then 1599 if (ee >= cmin .and. ee - cmin <= abs_erange2) then 1600 iflag = 1 !; write(std_out, *)"Adding conduction band", band, " with ee [eV]: ", ee * Ha_eV 1601 end if 1602 end if 1603 1604 if (iflag == 1) then 1605 gstore%brange_spin(1, spin) = min(gstore%brange_spin(1, spin), band) 1606 gstore%brange_spin(2, spin) = max(gstore%brange_spin(2, spin), band) 1607 1608 select case (gstore%kzone) 1609 case ("ibz") 1610 ik_bz = kibz2bz(ik_ibz) 1611 select_kbz_spin(ik_bz, spin) = iflag 1612 1613 case ("bz") 1614 call star_from_ibz_idx(ik_ibz, gstore%nkbz, kbz2ibz, nk_in_star, kstar_bz_inds) 1615 ABI_CHECK(nk_in_star > 0, "Something wrong in star_from_ibz_idx") 1616 do ii=1,nk_in_star 1617 ik_bz = kstar_bz_inds(ii) 1618 select_kbz_spin(ik_bz, spin) = iflag 1619 end do 1620 ABI_FREE(kstar_bz_inds) 1621 end select 1622 end if 1623 1624 end do ! band 1625 end do ! ik_ibz 1626 1627 !call wrtout(std_out, sjoin("brange_spin:", ltoa(gstore%brange_spin(:, spin)))) 1628 !call wrtout(std_out, sjoin("count_select_kbz:", itoa(count(select_kbz_spin == 1)))) 1629 1630 if (any(gstore%brange_spin(:, spin) == [huge(1), -huge(1)])) then 1631 ABI_ERROR("Empty list of states inside gstore_erange") 1632 end if 1633 end do ! spin 1634 1635 !call xmpi_sum(select_kbz_spin, comm, ierr) 1636 call recompute_select_qbz_spin(gstore, qbz, qbz2ibz, qibz2bz, kbz, kibz, kbz2ibz, kibz2bz, & 1637 select_kbz_spin, select_qbz_spin) 1638 1639 call gaps%free() 1640 1641 end subroutine gstore_filter_erange__
m_gstore/gstore_filter_fs_tetra__ [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_filter_fs_tetra__
FUNCTION
Compute delta(e_k - e_F) with the tetrahedron method. Use weights to filter k-points. Include only those q-points such that there exists at least one k on the FS with k + q on the FS. Also, store and precompute gstore%delta_ef_kibz_spin(max_nb, gstore%nkibz, nsppol) to be used to filter inside gstore%compute.
INPUTS
OUTPUT
SOURCE
1395 subroutine gstore_filter_fs_tetra__(gstore, qbz, qbz2ibz, qibz2bz, kbz, kibz, kbz2ibz, kibz2bz, & 1396 select_qbz_spin, select_kbz_spin) 1397 1398 !Arguments ------------------------------------ 1399 !scalars 1400 class(gstore_t),target,intent(inout) :: gstore 1401 real(dp),intent(in) :: qbz(3, gstore%nqbz) 1402 real(dp),target,intent(in) :: kbz(3, gstore%nqbz), kibz(3, gstore%nkibz) 1403 integer,intent(in) :: qbz2ibz(6,gstore%nqbz), qibz2bz(gstore%nqibz) 1404 integer,intent(in) :: kbz2ibz(6,gstore%nkbz), kibz2bz(gstore%nkibz) 1405 integer,intent(out) :: select_qbz_spin(gstore%nqbz, gstore%nsppol) 1406 integer,intent(out) :: select_kbz_spin(gstore%nkbz, gstore%nsppol) 1407 1408 !Local variables------------------------------- 1409 !scalars 1410 integer,parameter :: tetra_opt0 = 0 1411 integer :: nsppol, ierr, cnt, spin, band, ib, ii, max_nb, all_nproc, my_rank, comm, ebands_timrev 1412 integer :: ik_bz, ik_ibz, iflag, nk_in_star 1413 real(dp) :: max_occ 1414 character(len=80) :: errorstring 1415 type(ebands_t),pointer :: ebands 1416 type(crystal_t),pointer :: cryst 1417 type(htetra_t) :: ktetra 1418 !arrays 1419 integer,allocatable :: indkk(:), kstar_bz_inds(:) 1420 real(dp):: rlatt(3,3), klatt(3,3), delta_theta_ef(2) ! qpt(3), 1421 real(dp),allocatable :: eig_ibz(:) 1422 !---------------------------------------------------------------------- 1423 1424 comm = gstore%comm 1425 all_nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1426 1427 ! Use the tetrahedron method to filter k- and k+q points on the FS in metals 1428 ! and define gstore%brange_spin automatically. 1429 call wrtout(std_out, sjoin(" Filtering k-points using:", gstore%kfilter)) 1430 1431 ! NB: here we recompute brange_spin 1432 cryst => gstore%cryst 1433 ebands => gstore%ebands 1434 nsppol = gstore%nsppol 1435 1436 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 1437 1438 call ebands_get_bands_e0(ebands, ebands%fermie, gstore%brange_spin, ierr) 1439 ABI_CHECK(ierr == 0, "Error in ebands_get_bands_e0") 1440 1441 ABI_MALLOC(indkk, (gstore%nkbz)) 1442 indkk(:) = kbz2ibz(1, :) 1443 1444 rlatt = ebands%kptrlatt; call matr3inv(rlatt, klatt) 1445 call htetra_init(ktetra, indkk, gstore%cryst%gprimd, klatt, kbz, gstore%nkbz, gstore%kibz, gstore%nkibz, & 1446 ierr, errorstring, gstore%comm) 1447 ABI_CHECK(ierr == 0, errorstring) 1448 1449 ABI_MALLOC(eig_ibz, (gstore%nkibz)) 1450 max_occ = two / (ebands%nspinor * nsppol) 1451 select_kbz_spin = 0 1452 1453 max_nb = maxval(gstore%brange_spin(2, :) - gstore%brange_spin(1, :) + 1) 1454 ABI_CALLOC(gstore%delta_ef_kibz_spin, (max_nb, gstore%nkibz, gstore%nsppol)) 1455 1456 cnt = 0 1457 do spin=1,nsppol 1458 do band=gstore%brange_spin(1, spin), gstore%brange_spin(2, spin) 1459 ib = band - gstore%brange_spin(1, spin) + 1 1460 eig_ibz = ebands%eig(band, :, spin) 1461 1462 do ik_ibz=1,gstore%nkibz 1463 cnt = cnt + 1; if (mod(cnt, all_nproc) /= my_rank) cycle ! MPI parallelism inside comm 1464 1465 call ktetra%get_onewk_wvals(ik_ibz, tetra_opt0, 1, [ebands%fermie], max_occ, & 1466 gstore%nkibz, eig_ibz, delta_theta_ef) 1467 1468 gstore%delta_ef_kibz_spin(ib, ik_ibz, spin) = delta_theta_ef(1) 1469 1470 iflag = merge(1, 0, abs(delta_theta_ef(1)) > zero) 1471 1472 ! Use iflag to filter k-points. 1473 select case (gstore%kzone) 1474 case ("ibz") 1475 ik_bz = kibz2bz(ik_ibz) 1476 select_kbz_spin(ik_bz, spin) = select_kbz_spin(ik_bz, spin) + iflag 1477 1478 case ("bz") 1479 call star_from_ibz_idx(ik_ibz, gstore%nkbz, kbz2ibz, nk_in_star, kstar_bz_inds) 1480 ABI_CHECK(nk_in_star > 0, "Something wrong in star_from_ibz_idx") 1481 do ii=1,nk_in_star 1482 ik_bz = kstar_bz_inds(ii) 1483 select_kbz_spin(ik_bz, spin) = select_kbz_spin(ik_bz, spin) + iflag 1484 end do 1485 ABI_FREE(kstar_bz_inds) 1486 end select 1487 end do 1488 1489 end do 1490 end do 1491 1492 !call ktetra%print(std_out) 1493 1494 call xmpi_sum(select_kbz_spin, comm, ierr) 1495 call xmpi_sum(gstore%delta_ef_kibz_spin, comm, ierr) 1496 1497 ! Now the tricky part as we want to remove q-points that 1498 ! do not lead to any scattering process between two states on the FS 1499 ! Remember that k+q is always a sub-mesh of the input ebands k-mesh. 1500 1501 call recompute_select_qbz_spin(gstore, qbz, qbz2ibz, qibz2bz, kbz, kibz, kbz2ibz, kibz2bz, & 1502 select_kbz_spin, select_qbz_spin) 1503 1504 ABI_FREE(eig_ibz) 1505 ABI_FREE(indkk) 1506 call ktetra%free() 1507 1508 end subroutine gstore_filter_fs_tetra__
m_gstore/gstore_free [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_free
FUNCTION
Free dynamic memory.
SOURCE
2190 subroutine gstore_free(gstore) 2191 2192 !Arguments ------------------------------------ 2193 class(gstore_t),intent(inout) :: gstore 2194 2195 !Local variables------------------------------- 2196 integer :: my_is 2197 !---------------------------------------------------------------------- 2198 2199 do my_is=1,gstore%my_nspins 2200 call gstore%gqk(my_is)%free() 2201 end do 2202 ABI_SFREE(gstore%gqk) 2203 2204 ABI_SFREE(gstore%delta_ef_kibz_spin) 2205 ABI_SFREE(gstore%qibz) 2206 ABI_SFREE(gstore%wtq) 2207 ABI_SFREE(gstore%my_spins) 2208 ABI_SFREE(gstore%brange_spin) 2209 ABI_SFREE(gstore%glob_nk_spin) 2210 ABI_SFREE(gstore%glob_nq_spin) 2211 ABI_SFREE(gstore%erange_spin) 2212 2213 call gstore%krank_ibz%free() 2214 2215 end subroutine gstore_free
m_gstore/gstore_from_ncpath [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_from_ncpath
FUNCTION
Reconstruct a gstore object from a GSTORE.nc netcdf file.
INPUTS
SOURCE
3289 subroutine gstore_from_ncpath(gstore, path, with_cplex, dtset, cryst, ebands, ifc, comm) 3290 3291 !Arguments ------------------------------------ 3292 class(gstore_t),target,intent(out) :: gstore 3293 character(len=*),intent(in) :: path 3294 integer,intent(in) :: with_cplex 3295 type(dataset_type),intent(in) :: dtset 3296 integer,intent(in) :: comm 3297 class(crystal_t),target,intent(in) :: cryst 3298 class(ebands_t),target,intent(in) :: ebands 3299 class(ifc_type),target,intent(in) :: ifc 3300 3301 !Local variables------------------------------- 3302 !scalars 3303 integer,parameter :: master = 0 3304 integer :: my_rank, ncid, spin, spin_ncid, nproc, ierr, fform, max_nb, ib, natom, natom3 3305 integer :: max_nq, max_nk, gstore_cplex, ncerr, my_is, my_iq, iq_glob, my_ik, ik_glob, my_ip, ipert 3306 integer :: iq_ibz, isym_q, trev_q, tsign_q, g0_q(3) 3307 real(dp) :: cpu, wall, gflops 3308 character(len=10) :: priority 3309 logical :: store_phdispl, isirr_q 3310 type(hdr_type) :: wfk0_hdr 3311 type(crystal_t) :: gstore_cryst 3312 type(gqk_t),pointer :: gqk 3313 !arrays 3314 integer :: ibuffer(9), nproc_spin(ebands%nsppol), comm_spin(ebands%nsppol), brange_spin(2, ebands%nsppol) 3315 integer,allocatable :: qglob2bz(:,:), kglob2bz(:,:), qbz2ibz(:,:), kbz2ibz(:,:) 3316 real(dp) :: qq_ibz(3) 3317 real(dp),allocatable :: gwork_q(:,:,:,:,:), slice_bb(:,:,:) 3318 real(dp),allocatable :: phfreqs_ibz(:,:), pheigvec_cart_ibz(:,:,:,:,:) 3319 real(dp),allocatable :: pheigvec_cart_qbz(:,:,:,:), displ_cart_qbz(:,:,:,:) 3320 3321 ! ************************************************************************* 3322 3323 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 3324 3325 ! Set basic parameters. 3326 gstore%comm = comm 3327 gstore%nsppol = dtset%nsppol 3328 gstore%path = path 3329 3330 ! Get references to other data structures. 3331 gstore%cryst => cryst 3332 gstore%ebands => ebands 3333 gstore%ifc => ifc 3334 gstore%kibz => ebands%kptns 3335 natom = cryst%natom 3336 natom3 = cryst%natom * 3 3337 3338 ABI_CALLOC(gstore%erange_spin, (2, gstore%nsppol)) 3339 ABI_MALLOC(gstore%glob_nk_spin, (gstore%nsppol)) 3340 ABI_MALLOC(gstore%glob_nq_spin, (gstore%nsppol)) 3341 3342 ! ======================================================= 3343 ! Master node reads basic objects and gstore dimensions 3344 ! ======================================================= 3345 3346 if (my_rank == master) then 3347 call wrtout([std_out, ab_out], sjoin(" Initializing gstore object from:", path, ch10)) 3348 ABI_CHECK(path /= ABI_NOFILE, "Use getgstore_filepath to specify the path to GSTORE.nc") 3349 NCF_CHECK(nctk_open_read(ncid, path, xmpi_comm_self)) 3350 !NCF_CHECK(cryst%ncwrite(ncid)) 3351 3352 ! TODO Should compare ebands_file with input ebands 3353 !NCF_CHECK(ebands_ncwrite(ebands, ncid)) 3354 3355 call hdr_ncread(wfk0_hdr, ncid, fform) 3356 ABI_CHECK(fform /= 0, sjoin("Error while reading:", path)) 3357 3358 ! Read gstore dimensions 3359 NCF_CHECK(nctk_get_dim(ncid, "gstore_cplex", gstore_cplex)) 3360 NCF_CHECK(nctk_get_dim(ncid, "gstore_nkibz", gstore%nkibz)) 3361 NCF_CHECK(nctk_get_dim(ncid, "gstore_nkbz", gstore%nkbz)) 3362 NCF_CHECK(nctk_get_dim(ncid, "gstore_nqibz", gstore%nqibz)) 3363 NCF_CHECK(nctk_get_dim(ncid, "gstore_nqbz", gstore%nqbz)) 3364 NCF_CHECK(nctk_get_dim(ncid, "gstore_max_nq", max_nq)) 3365 NCF_CHECK(nctk_get_dim(ncid, "gstore_max_nk", max_nk)) 3366 NCF_CHECK(nctk_get_dim(ncid, "gstore_max_nb", max_nb)) 3367 3368 ! Read gstore variables 3369 NCF_CHECK(nf90_get_var(ncid, vid("gstore_with_vk"), gstore%with_vk)) 3370 NCF_CHECK(nf90_get_var(ncid, vid("gstore_qptopt"), gstore%qptopt)) 3371 NCF_CHECK(nf90_get_var(ncid, vid("gstore_kzone"), gstore%kzone)) 3372 call replace_ch0(gstore%kzone) 3373 NCF_CHECK(nf90_get_var(ncid, vid("gstore_qzone"), gstore%qzone)) 3374 call replace_ch0(gstore%qzone) 3375 NCF_CHECK(nf90_get_var(ncid, vid("gstore_kfilter"), gstore%kfilter)) 3376 call replace_ch0(gstore%kfilter) 3377 NCF_CHECK(nf90_get_var(ncid, vid("gstore_wfk0_path"), gstore%wfk0_path)) 3378 call replace_ch0(gstore%wfk0_path) 3379 3380 ABI_MALLOC(gstore%qibz, (3, gstore%nqibz)) 3381 ABI_MALLOC(gstore%wtq, (gstore%nqibz)) 3382 NCF_CHECK(nf90_get_var(ncid, vid("gstore_brange_spin"), brange_spin)) 3383 NCF_CHECK(nf90_get_var(ncid, vid("gstore_erange_spin"), gstore%erange_spin)) 3384 NCF_CHECK(nf90_get_var(ncid, vid("gstore_qibz"), gstore%qibz)) 3385 NCF_CHECK(nf90_get_var(ncid, vid("gstore_wtq"), gstore%wtq)) 3386 3387 ABI_MALLOC(qbz2ibz, (6, gstore%nqbz)) 3388 ABI_MALLOC(kbz2ibz, (6, gstore%nkbz)) 3389 NCF_CHECK(nf90_get_var(ncid, vid("gstore_qbz2ibz"), qbz2ibz)) 3390 NCF_CHECK(nf90_get_var(ncid, vid("gstore_kbz2ibz"), kbz2ibz)) 3391 3392 ABI_MALLOC(qglob2bz, (max_nq, gstore%nsppol)) 3393 ABI_MALLOC(kglob2bz, (max_nk, gstore%nsppol)) 3394 NCF_CHECK(nf90_get_var(ncid, vid("gstore_qglob2bz"), qglob2bz)) 3395 NCF_CHECK(nf90_get_var(ncid, vid("gstore_kglob2bz"), kglob2bz)) 3396 NCF_CHECK(nf90_get_var(ncid, vid("gstore_glob_nq_spin"), gstore%glob_nq_spin)) 3397 NCF_CHECK(nf90_get_var(ncid, vid("gstore_glob_nk_spin"), gstore%glob_nk_spin)) 3398 3399 ! Read optional variables: 3400 if (gstore%kfilter == "fs_tetra") then 3401 ABI_MALLOC(gstore%delta_ef_kibz_spin, (max_nb, gstore%nkibz, gstore%nsppol)) 3402 NCF_CHECK(nf90_get_var(ncid, vid("gstore_delta_ef_kibz_spin"), gstore%delta_ef_kibz_spin)) 3403 end if 3404 3405 NCF_CHECK(nf90_close(ncid)) 3406 3407 ! ================= 3408 ! Consistency check 3409 ! ================= 3410 gstore_cryst = wfk0_hdr%get_crystal() 3411 !call gstore_cryst%print(header="crystal structure from WFK file") 3412 3413 if (cryst%compare(gstore_cryst, header=" Comparing input crystal with the one from GSTORE file") /= 0) then 3414 ABI_ERROR("Crystal structure from input and GSTORE do not agree! Check messages above!") 3415 end if 3416 call gstore_cryst%free() 3417 3418 if (gstore_cplex == 1 .and. with_cplex == 2) then 3419 ABI_ERROR("GSTORE file contains |g| while Abinit needs complex g. Regenerate GSTORE file with gstore_cplex 2") 3420 end if 3421 end if ! master 3422 3423 if (nproc > 1) then 3424 ! Broadcast the header. 3425 call wfk0_hdr%bcast(master, my_rank, comm) 3426 3427 ! Broadcast dimensions. 3428 if (my_rank == master) then 3429 ibuffer = [gstore_cplex, gstore%nkibz, gstore%nkbz, gstore%nqibz, gstore%nqbz, & 3430 gstore%with_vk, max_nq, max_nk, max_nb] 3431 end if 3432 call xmpi_bcast(ibuffer, master, comm, ierr) 3433 3434 if (my_rank /= master) then 3435 ! Other MPI procs need to store dims and allocate memory before the bcast. 3436 gstore_cplex = ibuffer(1) 3437 gstore%nkibz = ibuffer(2) 3438 gstore%nkbz = ibuffer(3) 3439 gstore%nqibz = ibuffer(4) 3440 gstore%nqbz = ibuffer(5) 3441 gstore%with_vk = ibuffer(6) 3442 max_nq = ibuffer(7) 3443 max_nk = ibuffer(8) 3444 max_nb = ibuffer(9) 3445 3446 ABI_MALLOC(gstore%qibz, (3, gstore%nqibz)) 3447 ABI_MALLOC(gstore%wtq, (gstore%nqibz)) 3448 ABI_MALLOC(qbz2ibz, (6, gstore%nqbz)) 3449 ABI_MALLOC(kbz2ibz, (6, gstore%nkbz)) 3450 ABI_MALLOC(qglob2bz, (max_nq, gstore%nsppol)) 3451 ABI_MALLOC(kglob2bz, (max_nk, gstore%nsppol)) 3452 end if 3453 3454 call xmpi_bcast(gstore%qptopt, master, comm, ierr) 3455 call xmpi_bcast(gstore%kzone, master, comm, ierr) 3456 call xmpi_bcast(gstore%qzone, master, comm, ierr) 3457 call xmpi_bcast(gstore%kfilter, master, comm, ierr) 3458 call xmpi_bcast(gstore%wfk0_path, master, comm, ierr) 3459 call xmpi_bcast(brange_spin, master, comm, ierr) 3460 call xmpi_bcast(gstore%erange_spin, master, comm, ierr) 3461 call xmpi_bcast(gstore%qibz, master, comm, ierr) 3462 call xmpi_bcast(gstore%wtq, master, comm, ierr) 3463 call xmpi_bcast(qbz2ibz, master, comm, ierr) 3464 call xmpi_bcast(kbz2ibz, master, comm, ierr) 3465 call xmpi_bcast(qglob2bz, master, comm, ierr) 3466 call xmpi_bcast(kglob2bz, master, comm, ierr) 3467 call xmpi_bcast(gstore%glob_nq_spin, master, comm, ierr) 3468 call xmpi_bcast(gstore%glob_nk_spin, master, comm, ierr) 3469 3470 if (gstore%kfilter == "fs_tetra") then 3471 if (my_rank /= master) then 3472 ABI_MALLOC(gstore%delta_ef_kibz_spin, (max_nb, gstore%nkibz, gstore%nsppol)) 3473 end if 3474 call xmpi_bcast(gstore%delta_ef_kibz_spin, master, comm, ierr) 3475 end if 3476 3477 end if 3478 3479 ! Construct crystal and ebands from the GS WFK file. 3480 !ebands_file = wfk_read_ebands(wfk0_path, comm, out_hdr=wfk0_hdr) 3481 call wfk0_hdr%vs_dtset(dtset) 3482 call wfk0_hdr%free() 3483 3484 ! Distribute spins, create indirect mapping to spin index and init gstore%brange_spin 3485 call gstore%distribute_spins__(ebands%mband, brange_spin, nproc_spin, comm_spin, comm) 3486 3487 ! Compute krank 3488 gstore%krank_ibz = krank_from_kptrlatt(gstore%nkibz, gstore%kibz, ebands%kptrlatt, compute_invrank=.False.) 3489 3490 ! Set MPI grid. Also define priority according to eph_task. 3491 call priority_from_eph_task(dtset%eph_task, priority) 3492 3493 call gstore%set_mpi_grid__(with_cplex, dtset%eph_np_pqbks, priority, nproc_spin, comm_spin) 3494 3495 ! At this point, we have the Cartesian grid (one per spin if any) 3496 ! and we can finally allocate and distribute other arrays. 3497 call gstore%malloc__(with_cplex, max_nq, qglob2bz, max_nk, kglob2bz, qbz2ibz, kbz2ibz) 3498 3499 ABI_FREE(qglob2bz) 3500 ABI_FREE(kglob2bz) 3501 ABI_FREE(qbz2ibz) 3502 ABI_FREE(kbz2ibz) 3503 3504 call xmpi_comm_free(comm_spin) 3505 3506 ! Now we read the big arrays with MPI-IO and hdf5 groups. 3507 ! Note the loop over spin as each gqk has its own dimensions. 3508 ! Recall the shape on the arrays: 3509 ! 3510 ! In memory, we have allocated: 3511 ! 3512 ! my_g(my_npert, nb, my_nq, nb, my_nk) if with_cplex == 2 (complex array) 3513 ! 3514 ! or 3515 ! 3516 ! my_g2(my_npert, nb, my_nq, nb, my_nk) if with_cplex == 1 3517 ! 3518 ! On disk, we have: 3519 ! 3520 ! nctkarr_t("gvals", "dp", "gstore_cplex, nb, nb, natom3, glob_nk, glob_nq") 3521 ! 3522 ! so gstore_cplex == 1 and with_cplex == 2 are not compatible. 3523 ! 3524 call cwtime(cpu, wall, gflops, "start") 3525 3526 ! ======================================== 3527 ! Load phonon frequencies and eigenvectors 3528 ! ======================================== 3529 3530 ! nctkarr_t("phfreqs_ibz", "dp", "natom3, gstore_nqibz") 3531 ! nctkarr_t("pheigvec_cart_ibz", "dp", "two, three, natom, natom3, gstore_nqibz") 3532 ABI_MALLOC(phfreqs_ibz, (natom3, gstore%nqibz)) 3533 ABI_MALLOC(pheigvec_cart_ibz, (2, 3, cryst%natom, cryst%natom * 3, gstore%nqibz)) 3534 ABI_MALLOC(pheigvec_cart_qbz, (2, 3, cryst%natom, cryst%natom * 3)) 3535 ABI_MALLOC(displ_cart_qbz, (2, 3, cryst%natom, cryst%natom * 3)) 3536 3537 NCF_CHECK(nctk_open_read(ncid, gstore%path, gstore%comm)) 3538 3539 if (nproc > 1) then 3540 NCF_CHECK(nctk_set_collective(ncid, vid("phfreqs_ibz"))) 3541 end if 3542 NCF_CHECK(nf90_get_var(ncid, vid("phfreqs_ibz"), phfreqs_ibz)) 3543 if (nproc > 1) then 3544 NCF_CHECK(nctk_set_collective(ncid, vid("pheigvec_cart_ibz"))) 3545 end if 3546 NCF_CHECK(nf90_get_var(ncid, vid("pheigvec_cart_ibz"), pheigvec_cart_ibz)) 3547 NCF_CHECK(nf90_close(ncid)) 3548 3549 store_phdispl = .True. 3550 do my_is=1,gstore%my_nspins 3551 gqk => gstore%gqk(my_is) 3552 3553 ABI_MALLOC(gqk%my_wnuq, (gqk%my_npert, gqk%my_nq)) 3554 if (store_phdispl) then 3555 ABI_MALLOC(gqk%my_displ_cart, (2, 3, cryst%natom, gqk%my_npert, gqk%my_nq)) 3556 end if 3557 3558 do my_iq=1,gqk%my_nq 3559 iq_ibz = gqk%my_q2ibz(1, my_iq); isym_q = gqk%my_q2ibz(2, my_iq) 3560 trev_q = gqk%my_q2ibz(6, my_iq); g0_q = gqk%my_q2ibz(3:5, my_iq) 3561 !isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0)) 3562 isirr_q = (isym_q == 1 .and. trev_q == 0) 3563 tsign_q = 1; if (trev_q == 1) tsign_q = -1 3564 qq_ibz = gstore%qibz(:, iq_ibz) 3565 call pheigvec_rotate(cryst, qq_ibz, isym_q, trev_q, pheigvec_cart_ibz(:,:,:,:,iq_ibz), & 3566 pheigvec_cart_qbz, displ_cart_qbz) 3567 3568 gqk%my_wnuq(:, my_iq) = phfreqs_ibz(gqk%my_iperts(:), iq_ibz) 3569 if (store_phdispl) gqk%my_displ_cart(:,:,:,:,my_iq) = displ_cart_qbz(:,:,:,gqk%my_iperts(:)) 3570 end do ! my_iq 3571 end do ! my_is 3572 3573 ABI_FREE(phfreqs_ibz) 3574 ABI_SFREE(pheigvec_cart_ibz) 3575 ABI_SFREE(displ_cart_qbz) 3576 ABI_FREE(pheigvec_cart_qbz) 3577 3578 do spin=1,gstore%nsppol 3579 my_is = gstore%spin2my_is(spin) 3580 3581 if (my_is /= 0) then 3582 gqk => gstore%gqk(my_is) 3583 3584 NCF_CHECK(nctk_open_read(ncid, gstore%path, gqk%grid_comm%value)) 3585 NCF_CHECK(nf90_inq_ncid(ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid)) 3586 3587 !NCF_CHECK(nf90_get_var(spin_ncid, spin_vid("bstart"), gqk%bstart)) 3588 !NCF_CHECK(nf90_get_var(spin_ncid, spin_vid("bstop"), gqk%bstop)) 3589 3590 ! gstore_cplex defines the data on disk while 3591 ! cplex defines what we want to store in memory 3592 ABI_MALLOC_OR_DIE(gwork_q, (gstore_cplex, gqk%nb, gqk%nb, gqk%natom3, gqk%glob_nk), ierr) 3593 ABI_MALLOC(slice_bb, (gstore_cplex, gqk%nb, gqk%nb)) 3594 3595 do my_iq=1,gqk%my_nq 3596 iq_glob = my_iq + gqk%my_qstart - 1 3597 3598 ! Read q-slice (individual IO) 3599 ncerr = nf90_get_var(spin_ncid, spin_vid("gvals"), gwork_q, start=[1, 1, 1, 1, 1, iq_glob]) ! count=[]) 3600 NCF_CHECK(ncerr) 3601 3602 do my_ik=1,gqk%my_nk 3603 ik_glob = my_ik + gqk%my_kstart - 1 3604 do my_ip=1,gqk%my_npert 3605 ipert = gqk%my_iperts(my_ip) 3606 slice_bb = gwork_q(:,:,:, ipert, ik_glob) 3607 3608 ! Put data in the right place and handle conversion g --> |g|^2 3609 if (with_cplex == gstore_cplex) then 3610 if (with_cplex == 1) gqk%my_g2(my_ip,:,my_iq,:,my_ik) = slice_bb(1,:,:) 3611 if (with_cplex == 2) gqk%my_g(my_ip,:,my_iq,:,my_ik) = slice_bb(1,:,:) + j_dpc * slice_bb(2,:,:) 3612 else 3613 if (with_cplex == 1 .and. gstore_cplex == 2) then 3614 gqk%my_g2(my_ip, :, my_iq, :, my_ik) = slice_bb(1,:,:) ** 2 + slice_bb(2,:,:) ** 2 3615 else 3616 ABI_ERROR("Conversion from g2 on file to g_cplx in memory is not possible!") 3617 end if 3618 end if 3619 3620 end do 3621 end do 3622 3623 end do ! my_iq 3624 3625 ABI_FREE(gwork_q) 3626 ABI_FREE(slice_bb) 3627 3628 ! ============================================== 3629 ! Read matrix elements of the velocity operator 3630 ! ============================================== 3631 if (gstore%with_vk == 1) then 3632 if (gqk%grid_comm%nproc > 1) then 3633 NCF_CHECK(nctk_set_collective(spin_ncid, spin_vid("vk_cart_ibz"))) 3634 end if 3635 NCF_CHECK(nf90_get_var(spin_ncid, spin_vid("vk_cart_ibz"), gqk%vk_cart_ibz)) 3636 3637 else if (gstore%with_vk == 2) then 3638 if (gqk%grid_comm%nproc > 1) then 3639 NCF_CHECK(nctk_set_collective(spin_ncid, spin_vid("vkmat_cart_ibz"))) 3640 end if 3641 NCF_CHECK(nf90_get_var(spin_ncid, spin_vid("vkmat_cart_ibz"), gqk%vkmat_cart_ibz)) 3642 3643 ! Tranfer diagonal terms to vk_cart_ibz. 3644 do ib=1,gqk%nb 3645 gqk%vk_cart_ibz(:, ib, :) = gqk%vkmat_cart_ibz(1, :, ib, ib, :) 3646 end do 3647 end if 3648 3649 NCF_CHECK(nf90_close(ncid)) 3650 end if 3651 3652 call xmpi_barrier(gstore%comm) 3653 end do ! spin 3654 3655 call cwtime_report(" gstore_from_ncpath", cpu, wall, gflops) 3656 3657 contains 3658 integer function vid(vname) 3659 character(len=*),intent(in) :: vname 3660 vid = nctk_idname(ncid, vname) 3661 end function vid 3662 3663 integer function spin_vid(vname) 3664 character(len=*),intent(in) :: vname 3665 spin_vid = nctk_idname(spin_ncid, vname) 3666 end function spin_vid 3667 3668 end subroutine gstore_from_ncpath
m_gstore/gstore_get_a2fw [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_get_a2fw
FUNCTION
Compute a^2F(omega)
INPUTS
OUTPUT
SOURCE
2103 subroutine gstore_get_a2fw(gstore, dtset, nw, wmesh, a2fw) 2104 2105 !Arguments ------------------------------------ 2106 class(gstore_t),target,intent(inout) :: gstore 2107 type(dataset_type),intent(in) :: dtset 2108 integer,intent(in) :: nw 2109 real(dp),intent(in) :: wmesh(nw) 2110 real(dp),intent(out) :: a2fw(nw) 2111 2112 !Local variables------------------------------- 2113 integer :: my_is, my_ik, my_iq, my_ip, ib_k, ib_kq, ierr 2114 real(dp) :: g2, w_nuq, weight_k, weight_q, cpu, wall, gflops 2115 type(gqk_t), pointer :: gqk 2116 !arrays 2117 real(dp) :: qpt(3) 2118 real(dp),allocatable :: dbldelta_q(:,:,:), g2_mnkp(:,:,:,:), deltaw_nuq(:) 2119 2120 !---------------------------------------------------------------------- 2121 2122 call wrtout(std_out, sjoin(" Computing a^2F(w) with ph_smear:", ftoa(dtset%ph_smear * Ha_meV), "(meV)")) 2123 ABI_CHECK(gstore%qzone == "bz", "gstore_get_lambda_iso_iw assumes qzone == `bz`") 2124 call cwtime(cpu, wall, gflops, "start") 2125 2126 ABI_MALLOC(deltaw_nuq, (nw)) 2127 a2fw = zero 2128 2129 do my_is=1,gstore%my_nspins 2130 gqk => gstore%gqk(my_is) 2131 2132 ABI_CHECK(allocated(gqk%my_g2), "my_g2 is not allocated") 2133 ABI_CHECK(allocated(gqk%my_wnuq), "my_wnuq is not allocated") 2134 2135 ! Weights for delta(e_{m k+q}) delta(e_{n k}) for my list of k-points. 2136 ABI_MALLOC(dbldelta_q, (gqk%nb, gqk%nb, gqk%my_nk)) 2137 ABI_MALLOC(g2_mnkp, (gqk%nb, gqk%nb, gqk%my_nk, gqk%my_npert)) 2138 2139 do my_iq=1,gqk%my_nq 2140 ! Compute integration weights for the double delta. 2141 call gqk%dbldelta_qpt(my_iq, gstore, dtset%eph_intmeth, dtset%eph_fsmear, qpt, weight_q, dbldelta_q) 2142 2143 ! Copy data to improve memory access in the loops below. 2144 do my_ip=1,gqk%my_npert 2145 g2_mnkp(:,:,:,my_ip) = gqk%my_g2(my_ip,:,my_iq,:,:) 2146 end do 2147 2148 do my_ip=1,gqk%my_npert 2149 w_nuq = gqk%my_wnuq(my_ip, my_iq) 2150 deltaw_nuq = gaussian(wmesh - w_nuq, dtset%ph_smear) 2151 do my_ik=1,gqk%my_nk 2152 !weight_k = gqk%my_kweight(my_ik, gstore) 2153 weight_k = gqk%my_wtk(my_ik) 2154 do ib_k=1,gqk%nb 2155 do ib_kq=1,gqk%nb 2156 g2 = g2_mnkp(ib_kq, ib_k, my_ik, my_ip) 2157 a2fw(:) = a2fw(:) + deltaw_nuq(:) * g2 * weight_k * weight_q * dbldelta_q(ib_kq, ib_k, my_ik) 2158 end do 2159 end do 2160 end do 2161 end do 2162 end do ! my_iq 2163 2164 ABI_FREE(dbldelta_q) 2165 ABI_FREE(g2_mnkp) 2166 end do ! my_is 2167 2168 ABI_FREE(deltaw_nuq) 2169 2170 ! Take into account collinear spin and N(eF) TODO 2171 a2fw = a2fw * (two / (gstore%nsppol * dtset%nspinor)) 2172 call xmpi_sum(a2fw, gstore%comm, ierr) 2173 2174 call cwtime_report(" gstore_get_a2fw", cpu, wall, gflops) 2175 2176 end subroutine gstore_get_a2fw
m_gstore/gstore_get_lambda_iso_iw [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_get_lambda_iso_iw
FUNCTION
Compute isotropic lambda along the imaginary axis
INPUTS
OUTPUT
SOURCE
2020 subroutine gstore_get_lambda_iso_iw(gstore, dtset, nw, imag_w, lambda) 2021 2022 !Arguments ------------------------------------ 2023 class(gstore_t),target,intent(inout) :: gstore 2024 type(dataset_type),intent(in) :: dtset 2025 integer,intent(in) :: nw 2026 real(dp),intent(in) :: imag_w(nw) 2027 real(dp),intent(out) :: lambda(nw) 2028 2029 !Local variables------------------------------- 2030 integer :: my_is, my_ik, my_iq, my_ip, ib_k, ib_kq, ierr 2031 real(dp) :: g2, w_nuq, weight_k, weight_q 2032 type(gqk_t), pointer :: gqk 2033 !arrays 2034 real(dp) :: qpt(3) 2035 real(dp),allocatable :: dbldelta_q(:,:,:), g2_pmnk(:,:,:,:) 2036 2037 !---------------------------------------------------------------------- 2038 2039 ABI_CHECK(gstore%qzone == "bz", "gstore_get_lambda_iso_iw assumes qzone == `bz`") 2040 lambda = zero 2041 2042 do my_is=1,gstore%my_nspins 2043 gqk => gstore%gqk(my_is) 2044 2045 ABI_CHECK(allocated(gqk%my_g2), "my_g2 is not allocated") 2046 ABI_CHECK(allocated(gqk%my_wnuq), "my_wnuq is not allocated") 2047 2048 ! Weights for delta(e_{m k+q}) delta(e_{n k}) for my list of k-points. 2049 ABI_MALLOC(dbldelta_q, (gqk%nb, gqk%nb, gqk%my_nk)) 2050 ABI_MALLOC(g2_pmnk, (gqk%my_npert, gqk%nb, gqk%nb, gqk%my_nk)) 2051 2052 do my_iq=1,gqk%my_nq 2053 ! Compute integration weights for the double delta. 2054 call gqk%dbldelta_qpt(my_iq, gstore, dtset%eph_intmeth, dtset%eph_fsmear, qpt, weight_q, dbldelta_q) 2055 2056 ! Copy data to improve memory access in the loops below. 2057 g2_pmnk = gqk%my_g2(:,:,my_iq,:,:) 2058 2059 do my_ik=1,gqk%my_nk 2060 !weight_k = gqk%my_kweight(my_ik, gstore) 2061 weight_k = gqk%my_wtk(my_ik) 2062 do ib_k=1,gqk%nb 2063 do ib_kq=1,gqk%nb 2064 do my_ip=1,gqk%my_npert 2065 g2 = g2_pmnk(my_ip, ib_kq, ib_k, my_ik) 2066 ! TODO: handle w_nuq ~ 0 2067 w_nuq = gqk%my_wnuq(my_ip, my_iq) 2068 lambda(:) = lambda(:) + & 2069 two * w_nuq / (imag_w(:) ** 2 + w_nuq ** 2) * g2 * weight_k * weight_q * dbldelta_q(ib_kq, ib_k, my_ik) 2070 end do 2071 end do 2072 end do 2073 end do 2074 end do ! my_iq 2075 2076 ABI_FREE(dbldelta_q) 2077 ABI_FREE(g2_pmnk) 2078 end do ! my_is 2079 2080 ! Take into account collinear spin 2081 lambda = lambda * (two / (gstore%nsppol * dtset%nspinor)) 2082 2083 call xmpi_sum(lambda, gstore%comm, ierr) 2084 2085 end subroutine gstore_get_lambda_iso_iw
m_gstore/gstore_get_mpw_gmax [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_get_mpw_gmax
FUNCTION
INPUTS
OUTPUT
SOURCE
1860 subroutine gstore_get_mpw_gmax(gstore, ecut, mpw, gmax) 1861 1862 !Arguments ------------------------------------ 1863 class(gstore_t),target,intent(in) :: gstore 1864 real(dp),intent(in) :: ecut 1865 integer,intent(out) :: mpw, gmax(3) 1866 1867 !Local variables------------------------------- 1868 integer,parameter :: istwfk1 = 1 1869 integer :: my_is, my_ik, my_iq, spin, onpw, ii, ipw, ierr, my_mpw 1870 real(dp) :: weight_q, cpu, wall, gflops !weight_k, 1871 type(gqk_t),pointer :: gqk 1872 !arrays 1873 integer :: my_gmax(3) 1874 integer,allocatable :: gtmp(:,:) 1875 real(dp) :: kk(3), qpt(3) 1876 1877 !---------------------------------------------------------------------- 1878 1879 mpw = 0; gmax = 0 1880 1881 ! TODO: This is an hotspot due to the double loop over k and q. 1882 ! Should use a geometrical approach to compute mpw and gmax. 1883 1884 call wrtout(std_out, " Computing mpw. This may take some time for dense k/q meshes...") 1885 call cwtime(cpu, wall, gflops, "start") 1886 1887 do my_is=1,gstore%my_nspins 1888 gqk => gstore%gqk(my_is) 1889 spin = gstore%my_spins(my_is) 1890 1891 do my_ik=1,gqk%my_nk 1892 kk = gqk%my_kpts(:, my_ik) 1893 1894 ! Compute G sphere, returning npw. Note istwfk == 1. 1895 call get_kg(kk, istwfk1, ecut, gstore%cryst%gmet, onpw, gtmp) 1896 mpw = max(mpw, onpw) 1897 do ipw=1,onpw 1898 do ii=1,3 1899 gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw))) 1900 end do 1901 end do 1902 ABI_FREE(gtmp) 1903 1904 do my_iq=1,gqk%my_nq 1905 call gqk%myqpt(my_iq, gstore, weight_q, qpt) 1906 1907 ! TODO: g0 umklapp here can enter into play! 1908 ! gmax could not be large enough! 1909 call get_kg(kk + qpt, 1, ecut, gstore%cryst%gmet, onpw, gtmp) 1910 mpw = max(mpw, onpw) 1911 do ipw=1,onpw 1912 do ii=1,3 1913 gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw))) 1914 end do 1915 end do 1916 ABI_FREE(gtmp) 1917 end do 1918 1919 end do ! my_ik 1920 end do ! my_is 1921 1922 my_mpw = mpw; call xmpi_max(my_mpw, mpw, gstore%comm, ierr) 1923 my_gmax = gmax; call xmpi_max(my_gmax, gmax, gstore%comm, ierr) 1924 1925 call wrtout(std_out, sjoin(' Optimal value of mpw: ', itoa(mpw))) 1926 call cwtime_report(" gmax and mpw", cpu, wall, gflops) 1927 1928 end subroutine gstore_get_mpw_gmax
m_gstore/gstore_init [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_init
FUNCTION
Initialize the object
INPUTS
path=Filename of the output GSTORE.nc file
OUTPUT
SOURCE
498 subroutine gstore_init(gstore, path, dtset, wfk0_hdr, cryst, ebands, ifc, comm) 499 500 !Arguments ------------------------------------ 501 !scalars 502 class(gstore_t),target,intent(out) :: gstore 503 integer,intent(in) :: comm 504 character(len=*),intent(in) :: path 505 type(dataset_type),intent(in) :: dtset 506 type(hdr_type),intent(in) :: wfk0_hdr 507 class(crystal_t),target,intent(in) :: cryst 508 class(ebands_t),target,intent(in) :: ebands 509 class(ifc_type),target,intent(in) :: ifc 510 511 !Local variables------------------------------- 512 !scalars 513 integer,parameter :: master = 0 514 integer :: all_nproc, my_rank, ierr, my_nshiftq, nsppol, spin, natom3, cnt, qtimrev 515 integer :: ik_ibz, ik_bz, iq_bz, iq_ibz, ebands_timrev, max_nq, max_nk, ncid, spin_ncid, ncerr, gstore_fform 516 real(dp) :: cpu, wall, gflops 517 character(len=10) :: priority 518 !character(len=5000) :: msg 519 type(krank_t) :: qrank 520 !arrays 521 integer :: ngqpt(3), qptrlatt(3,3), comm_spin(ebands%nsppol), nproc_spin(ebands%nsppol) 522 !integer :: glob_nk_spin(ebands%nsppol), glob_nq_spin(ebands%nsppol) 523 integer,allocatable :: qbz2ibz(:,:), kbz2ibz(:,:), kibz2bz(:), qibz2bz(:), qglob2bz(:,:), kglob2bz(:,:) 524 integer,allocatable :: select_qbz_spin(:,:), select_kbz_spin(:,:) !, done_qbz_spin(:,:) 525 real(dp):: my_shiftq(3,1) 526 real(dp),allocatable :: qbz(:,:), wtk(:), kibz(:,:), kbz(:,:) 527 !integer :: out_kptrlatt(3,3) 528 !real(dp),allocatable :: out_kibz(:,:), out_wtk(:) 529 530 !---------------------------------------------------------------------- 531 532 call cwtime(cpu, wall, gflops, "start") 533 all_nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 534 natom3 = 3 * cryst%natom; nsppol = ebands%nsppol 535 536 ! Set basic parameters. 537 gstore%comm = comm 538 gstore%nsppol = nsppol 539 gstore%path = path 540 541 ! Get references to other data structures. 542 gstore%cryst => cryst 543 gstore%ebands => ebands 544 gstore%ifc => ifc 545 gstore%kibz => ebands%kptns 546 547 ! Set metadata. 548 gstore%kzone = dtset%gstore_kzone; gstore%qzone = dtset%gstore_qzone; gstore%kfilter = dtset%gstore_kfilter 549 gstore%with_vk = dtset%gstore_with_vk 550 551 ABI_CALLOC(gstore%erange_spin, (2, nsppol)) 552 gstore%erange_spin = dtset%gstore_erange(:, 1:nsppol) 553 554 if (any(gstore%erange_spin /= zero)) then 555 ABI_CHECK(gstore%kfilter == "none", sjoin("kfilter should be none when erange is used while it is:", gstore%kfilter)) 556 gstore%kfilter = "erange" 557 end if 558 559 if (gstore%kzone == "ibz" .and. gstore%qzone == "ibz") then 560 ABI_ERROR("The combination kzone = 'ibz' and qzone = 'ibz' is not allowed!") 561 end if 562 563 ! TODO 564 !gstore%kptrlatt(3, 3) 565 !gstore%kshift(3, 1) 566 !gstore%qptrlatt(3, 3) 567 !gstore%qshift(3, 1) 568 569 ! Distribute spins, create indirect mapping to spin index and init gstore%brange_spin 570 ABI_CHECK_ILEQ(dtset%mband, ebands%mband, "dtset%mband > ebands%mband") 571 call gstore%distribute_spins__(dtset%mband, dtset%gstore_brange, nproc_spin, comm_spin, comm) 572 573 ! Define q-mesh: either from DVDB (no interpolation) or eph_ngqpt_fine (Fourier interpolation) 574 ngqpt = dtset%ddb_ngqpt; my_nshiftq = 1; my_shiftq(:,1) = dtset%ddb_shiftq 575 if (all(dtset%eph_ngqpt_fine /= 0)) then 576 ngqpt = dtset%eph_ngqpt_fine; my_shiftq = 0 577 end if 578 579 ! TODO: Should fix bz2ibz to use the same conventions as krank and listkk 580 ! NB: only sigmaph seems to be using this optional argument 581 582 ! Setup qIBZ, weights and BZ. 583 ! Assume qptopt == kptopt unless value is specified in input 584 qptrlatt = 0; qptrlatt(1, 1) = ngqpt(1); qptrlatt(2, 2) = ngqpt(2); qptrlatt(3, 3) = ngqpt(3) 585 gstore%qptopt = ebands%kptopt; if (dtset%qptopt /= 0) gstore%qptopt = dtset%qptopt 586 qtimrev = kpts_timrev_from_kptopt(gstore%qptopt) 587 588 !call wrtout(std_out, sjoin(" Generating q-IBZ for with qptopt:", itoa(gstore%qptopt))) 589 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, gstore%qptopt, my_nshiftq, my_shiftq, & 590 gstore%nqibz, gstore%qibz, gstore%wtq, gstore%nqbz, qbz) 591 !new_kptrlatt=gstore%qptrlatt, new_shiftk=gstore%qshift, 592 !bz2ibz=new%ind_qbz2ibz) # FIXME 593 594 ! HM: the bz2ibz produced above is incomplete, I do it here using listkk 595 ABI_MALLOC(qbz2ibz, (6, gstore%nqbz)) 596 597 qrank = krank_from_kptrlatt(gstore%nqibz, gstore%qibz, qptrlatt, compute_invrank=.False.) 598 599 if (kpts_map("symrec", qtimrev, cryst, qrank, gstore%nqbz, qbz, qbz2ibz) /= 0) then 600 ABI_ERROR("Cannot map qBZ to IBZ!") 601 end if 602 call qrank%free() 603 604 ! Order qbz by stars and rearrange entries in qbz2ibz table. 605 call kpts_pack_in_stars(gstore%nqbz, qbz, qbz2ibz) 606 607 !call kpts_print_kmap(std_out, qibz, qbz, qbz2ibz) 608 !do iq_bz=1,gstore%nqbz 609 ! print *, "iq_bz -> iq_ibz", qbz2ibz(1, iq_bz), qbz(:, iq_bz) 610 !end do 611 612 call get_ibz2bz(gstore%nqibz, gstore%nqbz, qbz2ibz, qibz2bz, ierr) 613 ABI_CHECK(ierr == 0, "Something wrong in symmetry tables for q-points!") 614 615 ! Get full BZ associated to ebands 616 call kpts_ibz_from_kptrlatt(cryst, ebands%kptrlatt, ebands%kptopt, ebands%nshiftk, ebands%shiftk, & 617 gstore%nkibz, kibz, wtk, gstore%nkbz, kbz) !, bz2ibz=bz2ibz) 618 !new_kptrlatt=gstore%kptrlatt, new_shiftk=gstore%kshift, 619 !bz2ibz=new%ind_qbz2ibz) # FIXME 620 621 ! In principle kibz should be equal to ebands%kptns 622 ABI_CHECK(gstore%nkibz == ebands%nkpt, "nkibz != ebands%nkpt") 623 ABI_CHECK(all(abs(gstore%kibz - kibz) < tol12), "ebands%kibz != kibz") 624 ABI_FREE(kibz) 625 626 ! Note symrel and use_symrec=.False. in get_mapping. 627 ! This means that this table can be used to symmetrize wavefunctions in cgtk_rotate. 628 ! TODO This ambiguity should be removed. Change cgtk_rotate so that we can use the symrec convention. 629 630 ABI_MALLOC(kbz2ibz, (6, gstore%nkbz)) 631 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 632 gstore%krank_ibz = krank_from_kptrlatt(gstore%nkibz, gstore%kibz, ebands%kptrlatt, compute_invrank=.False.) 633 if (kpts_map("symrel", ebands_timrev, cryst, gstore%krank_ibz, gstore%nkbz, kbz, kbz2ibz) /= 0) then 634 ABI_ERROR("Cannot map kBZ to IBZ!") 635 end if 636 637 call get_ibz2bz(gstore%nkibz, gstore%nkbz, kbz2ibz, kibz2bz, ierr) 638 ABI_CHECK(ierr == 0, "Something wrong in symmetry tables for k-points") 639 640 ! These tables are used to exclude q/k points 641 ! We use the full BZ because this mask can be also used when points are restricted to the IBZ 642 ! provided we convert from ik_ibz to ik_bz. Note that both arrays are initialized with zeros. 643 ABI_ICALLOC(select_qbz_spin, (gstore%nqbz, nsppol)) 644 ABI_ICALLOC(select_kbz_spin, (gstore%nkbz, nsppol)) 645 646 select case (gstore%kzone) 647 case ("ibz") 648 do spin=1,nsppol 649 do ik_ibz=1,gstore%nkibz 650 ik_bz = kibz2bz(ik_ibz); select_kbz_spin(ik_bz, spin) = 1 651 end do 652 end do 653 654 case ("bz") 655 select_kbz_spin = 1 656 657 case default 658 ABI_ERROR(sjoin("Invalid kzone:", gstore%kzone)) 659 end select 660 661 select case (gstore%qzone) 662 case ("ibz") 663 do spin=1,nsppol 664 do iq_ibz=1, gstore%nqibz 665 iq_bz = qibz2bz(iq_ibz); select_qbz_spin(iq_bz, spin) = 1 666 end do 667 end do 668 669 case ("bz") 670 select_qbz_spin = 1 671 672 case default 673 ABI_ERROR(sjoin("Invalid qzone:", gstore%qzone)) 674 end select 675 676 ! Here we filter the electronic wavevectors k and recompute select_qbz_spin and select_kbz_spin. 677 678 select case (gstore%kfilter) 679 case ("none") 680 continue 681 682 case ("erange") 683 call gstore%filter_erange__(qbz, qbz2ibz, qibz2bz, kbz, gstore%kibz, kbz2ibz, kibz2bz, & 684 select_qbz_spin, select_kbz_spin) 685 686 case ("fs_tetra") 687 ! Use the tetrahedron method to filter k- and k+q points on the FS in metals 688 ! and define gstore%brange_spin automatically. 689 call gstore%filter_fs_tetra__(qbz, qbz2ibz, qibz2bz, kbz, gstore%kibz, kbz2ibz, kibz2bz, & 690 select_qbz_spin, select_kbz_spin) 691 case default 692 ABI_ERROR(sjoin("Invalid kfilter:", gstore%kfilter)) 693 end select 694 695 ! Total number of k/q points for each spin after filtering (if any) 696 ABI_MALLOC(gstore%glob_nk_spin, (nsppol)) 697 ABI_MALLOC(gstore%glob_nq_spin, (nsppol)) 698 gstore%glob_nk_spin(:) = count(select_kbz_spin > 0, dim=1) 699 gstore%glob_nq_spin(:) = count(select_qbz_spin > 0, dim=1) 700 701 ! We need another table mapping the global index in the gqk matrix to the q/k index in the BZ 702 ! so that one can extract the symmetry tables computed above. 703 ! Again this is needed as the global sizes of the gqk matrix 704 ! is not necessarily equal to the size of the BZ/IBZ if we have filtered the wave vectors. 705 706 max_nq = maxval(gstore%glob_nq_spin) ! Max dim over spin 707 max_nk = maxval(gstore%glob_nk_spin) 708 ABI_ICALLOC(qglob2bz, (max_nq, nsppol)) 709 ABI_ICALLOC(kglob2bz, (max_nk, nsppol)) 710 711 do spin=1,nsppol 712 cnt = 0 713 do iq_bz=1,gstore%nqbz 714 if (select_qbz_spin(iq_bz, spin) /= 0) then 715 cnt = cnt + 1; qglob2bz(cnt, spin) = iq_bz 716 end if 717 end do 718 719 cnt = 0 720 do ik_bz=1,gstore%nkbz 721 if (select_kbz_spin(ik_bz, spin) /= 0) then 722 cnt = cnt + 1; kglob2bz(cnt, spin) = ik_bz 723 end if 724 end do 725 end do 726 727 ! ============================================= 728 ! Initialize gqk basic dimensions and MPI grid 729 ! ============================================= 730 !priority = "qk" 731 call priority_from_eph_task(dtset%eph_task, priority) 732 733 call gstore%set_mpi_grid__(dtset%gstore_cplex, dtset%eph_np_pqbks, priority, nproc_spin, comm_spin) 734 call xmpi_comm_free(comm_spin) 735 736 ! At this point, we have the Cartesian grid (one per spin if any) 737 ! and we can finally allocate and distribute other arrays. 738 call gstore%malloc__(0, max_nq, qglob2bz, max_nk, kglob2bz, qbz2ibz, kbz2ibz) 739 740 ! Initialize GSTORE.nc file i.e. define dimensions and arrays 741 ! Entries such as the e-ph matrix elements will be filled afterwards in gstore_compute. 742 ! Master node writes basic objects and gstore dimensions 743 744 if (my_rank == master) then 745 NCF_CHECK(nctk_open_create(ncid, gstore%path, xmpi_comm_self)) 746 747 ! Write the abinit header with metadata, structure and occupancies. 748 gstore_fform = fform_from_ext("GSTORE.nc") 749 NCF_CHECK(wfk0_hdr%ncwrite(ncid, gstore_fform, spinat=dtset%spinat, nc_define=.True.)) 750 !NCF_CHECK(gstore%cryst%ncwrite(ncid)) 751 752 ! Add eigenvalues and occupations. 753 NCF_CHECK(ebands_ncwrite(gstore%ebands, ncid)) 754 755 ! Write gstore dimensions 756 ncerr = nctk_def_dims(ncid, [ & 757 nctkdim_t("gstore_nkibz", gstore%nkibz), & 758 nctkdim_t("gstore_nkbz", gstore%nkbz), & 759 nctkdim_t("gstore_nqibz", gstore%nqibz), & 760 nctkdim_t("gstore_nqbz", gstore%nqbz), & 761 nctkdim_t("gstore_max_nq", max_nq), & 762 nctkdim_t("gstore_max_nk", max_nk), & 763 nctkdim_t("gstore_max_nb", maxval(gstore%brange_spin(2, :) - gstore%brange_spin(1, :) + 1) ), & 764 nctkdim_t("natom", gstore%cryst%natom), & 765 nctkdim_t("natom3", 3 * gstore%cryst%natom), & 766 nctkdim_t("gstore_cplex", dtset%gstore_cplex) & 767 ], defmode=.True.) 768 NCF_CHECK(ncerr) 769 770 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "gstore_with_vk", "gstore_qptopt"]) 771 NCF_CHECK(ncerr) 772 !ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "fermi_energy", "smearing_width"]) 773 !NCF_CHECK(ncerr) 774 775 ncerr = nctk_def_arrays(ncid, [ & 776 nctkarr_t("gstore_qibz", "dp", "three, gstore_nqibz"), & 777 nctkarr_t("gstore_qbz", "dp", "three, gstore_nqbz"), & 778 nctkarr_t("gstore_wtq", "dp", "gstore_nqibz"), & 779 nctkarr_t("gstore_kbz", "dp", "three, gstore_nkbz"), & 780 nctkarr_t("gstore_kzone", "c", "fnlen"), & 781 nctkarr_t("gstore_qzone", "c", "fnlen"), & 782 nctkarr_t("gstore_kfilter", "c", "fnlen"), & 783 nctkarr_t("gstore_wfk0_path", "c", "fnlen"), & 784 nctkarr_t("gstore_brange_spin", "i", "two, number_of_spins"), & 785 nctkarr_t("gstore_erange_spin", "dp", "two, number_of_spins"), & 786 nctkarr_t("phfreqs_ibz", "dp", "natom3, gstore_nqibz"), & 787 nctkarr_t("pheigvec_cart_ibz", "dp", "two, three, natom, natom3, gstore_nqibz"), & 788 nctkarr_t("gstore_glob_nq_spin", "i", "number_of_spins"), & 789 nctkarr_t("gstore_glob_nk_spin", "i", "number_of_spins"), & 790 nctkarr_t("gstore_done_qbz_spin", "i", "gstore_nqbz, number_of_spins"), & 791 nctkarr_t("gstore_kbz2ibz", "i", "six, gstore_nkbz"), & 792 nctkarr_t("gstore_qbz2ibz", "i", "six, gstore_nqbz"), & 793 nctkarr_t("gstore_qglob2bz", "i", "gstore_max_nq, number_of_spins"), & 794 nctkarr_t("gstore_kglob2bz", "i", "gstore_max_nk, number_of_spins") & 795 ]) 796 NCF_CHECK(ncerr) 797 798 ! internal table used to restart computation. Init with zeros. 799 NCF_CHECK(nf90_def_var_fill(ncid, vid("gstore_done_qbz_spin"), NF90_FILL, zero)) 800 801 ! Optional arrays 802 if (allocated(gstore%delta_ef_kibz_spin)) then 803 ncerr = nctk_def_arrays(ncid, & 804 nctkarr_t("gstore_delta_ef_kibz_spin", "dp", "gstore_max_nb, gstore_nkibz, number_of_spins")) 805 end if 806 NCF_CHECK(ncerr) 807 808 NCF_CHECK(nctk_set_datamode(ncid)) 809 NCF_CHECK(nf90_put_var(ncid, vid("gstore_with_vk"), gstore%with_vk)) 810 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qptopt"), gstore%qptopt)) 811 NCF_CHECK(nf90_put_var(ncid, vid("gstore_kzone"), trim(gstore%kzone))) 812 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qzone"), trim(gstore%qzone))) 813 NCF_CHECK(nf90_put_var(ncid, vid("gstore_kfilter"), trim(gstore%kfilter))) 814 815 ! Write arrays 816 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qibz"), gstore%qibz)) 817 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qbz"), qbz)) 818 NCF_CHECK(nf90_put_var(ncid, vid("gstore_wtq"), gstore%wtq)) 819 NCF_CHECK(nf90_put_var(ncid, vid("gstore_kbz"), kbz)) 820 NCF_CHECK(nf90_put_var(ncid, vid("gstore_brange_spin"), gstore%brange_spin)) 821 NCF_CHECK(nf90_put_var(ncid, vid("gstore_erange_spin"), gstore%erange_spin)) 822 NCF_CHECK(nf90_put_var(ncid, vid("gstore_glob_nq_spin"), gstore%glob_nq_spin)) 823 NCF_CHECK(nf90_put_var(ncid, vid("gstore_glob_nk_spin"), gstore%glob_nk_spin)) 824 NCF_CHECK(nf90_put_var(ncid, vid("gstore_kbz2ibz"), kbz2ibz)) 825 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qbz2ibz"), qbz2ibz)) 826 NCF_CHECK(nf90_put_var(ncid, vid("gstore_qglob2bz"), qglob2bz)) 827 NCF_CHECK(nf90_put_var(ncid, vid("gstore_kglob2bz"), kglob2bz)) 828 ! NB: kibz has been already written by ebands%ncwrite 829 830 if (allocated(gstore%delta_ef_kibz_spin)) then 831 NCF_CHECK(nf90_put_var(ncid, vid("gstore_delta_ef_kibz_spin"), gstore%delta_ef_kibz_spin)) 832 end if 833 834 do spin=1,gstore%nsppol 835 ! Create group for this spin. 836 NCF_CHECK(nf90_def_grp(ncid, strcat("gqk", "_spin", itoa(spin)), spin_ncid)) 837 838 ! Dimensions in gqk_spin group 839 ncerr = nctk_def_dims(spin_ncid, [ & 840 nctkdim_t("nb", gstore%brange_spin(2, spin) - gstore%brange_spin(1, spin) + 1), & 841 nctkdim_t("glob_nk", gstore%glob_nk_spin(spin)), & 842 nctkdim_t("glob_nq", gstore%glob_nq_spin(spin)) & 843 ], defmode=.True.) 844 NCF_CHECK(ncerr) 845 846 ! Define scalars 847 ncerr = nctk_def_iscalars(spin_ncid, [character(len=nctk_slen) :: "bstart"]) 848 NCF_CHECK(ncerr) 849 850 ! arrays in gqk_spin group with the precious stuff. Note glob dimensions 851 ncerr = nctk_def_arrays(spin_ncid, [ & 852 nctkarr_t("gvals", "dp", "gstore_cplex, nb, nb, natom3, glob_nk, glob_nq") & 853 ]) 854 NCF_CHECK(ncerr) 855 856 ! IMPORTANT: Init with zeros. 857 NCF_CHECK(nf90_def_var_fill(spin_ncid, vid_spin("gvals"), NF90_FILL, zero)) 858 859 select case(gstore%with_vk) 860 case (1) 861 ! Diagonal terms only 862 NCF_CHECK(nctk_def_arrays(spin_ncid, nctkarr_t("vk_cart_ibz", "dp", "three, nb, gstore_nkibz"))) 863 NCF_CHECK(nf90_def_var_fill(spin_ncid, vid_spin("vk_cart_ibz"), NF90_FILL, zero)) 864 case (2) 865 ! Full (nb x nb) matrix. 866 NCF_CHECK(nctk_def_arrays(spin_ncid, nctkarr_t("vkmat_cart_ibz", "dp", "two, three, nb, nb, gstore_nkibz"))) 867 NCF_CHECK(nf90_def_var_fill(spin_ncid, vid_spin("vkmat_cart_ibz"), NF90_FILL, zero)) 868 end select 869 870 ! Write (small) data 871 NCF_CHECK(nctk_set_datamode(spin_ncid)) 872 NCF_CHECK(nf90_put_var(spin_ncid, vid_spin("bstart"), gstore%brange_spin(1, spin))) 873 end do ! spin 874 875 NCF_CHECK(nf90_close(ncid)) 876 end if ! master 877 878 call xmpi_barrier(gstore%comm) 879 880 ABI_FREE(wtk) 881 ABI_FREE(kbz) 882 ABI_FREE(qbz) 883 ABI_FREE(qibz2bz) 884 ABI_FREE(kibz2bz) 885 ABI_FREE(select_qbz_spin) 886 ABI_FREE(select_kbz_spin) 887 ABI_FREE(qglob2bz) 888 ABI_FREE(kglob2bz) 889 ABI_FREE(qbz2ibz) 890 ABI_FREE(kbz2ibz) 891 892 call cwtime_report(" gstore_init:", cpu, wall, gflops) 893 894 contains 895 integer function vid(vname) 896 character(len=*),intent(in) :: vname 897 vid = nctk_idname(NCID, vname) 898 end function vid 899 integer function vid_spin(vname) 900 character(len=*),intent(in) :: vname 901 vid_spin = nctk_idname(SPIN_NCID, vname) 902 end function vid_spin 903 904 end subroutine gstore_init
m_gstore/gstore_malloc__ [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_malloc__
FUNCTION
Allocate memory.
INPUTS
OUTPUT
SOURCE
1258 subroutine gstore_malloc__(gstore, with_cplex, max_nq, qglob2bz, max_nk, kglob2bz, qbz2ibz, kbz2ibz) 1259 1260 !Arguments ------------------------------------ 1261 !scalars 1262 class(gstore_t),target,intent(inout) :: gstore 1263 integer,intent(in) :: with_cplex, max_nq, max_nk 1264 integer,intent(in) :: qglob2bz(max_nq, gstore%nsppol) 1265 integer,intent(in) :: kglob2bz(max_nk, gstore%nsppol) 1266 integer,intent(in) :: qbz2ibz(6, gstore%nqbz), kbz2ibz(6, gstore%nkbz) 1267 1268 !Local variables------------------------------- 1269 !scalars 1270 integer :: spin, my_is, ierr, my_iq, my_ik, iq_glob, iq_bz, ik_glob, ik_bz 1271 integer :: ik_ibz, isym_k, trev_k, tsign_k, g0_k(3) 1272 logical :: isirr_k 1273 real(dp) :: mem_mb 1274 type(gqk_t),pointer :: gqk 1275 !arrays 1276 integer,allocatable :: myq2glob(:), myk2glob(:) 1277 !---------------------------------------------------------------------- 1278 1279 !gstore%max_nb = maxval(gstore%brange_spin(2, :) - gstore%brange_spin(1, :) + 1) 1280 1281 do my_is=1,gstore%my_nspins 1282 gqk => gstore%gqk(my_is) 1283 spin = gstore%my_spins(my_is) 1284 1285 ! Split q-points and transfer symmetry tables. 1286 ! Note that glob_nq and glob_nk does not necessarily correspond to the size of the BZ 1287 ! First of all we have to consider kzone 1288 ! Even if kzone == "bz" we may have filtered the wavevectors e.g. Fermi surface. 1289 call xmpi_split_block(gqk%glob_nq, gqk%qpt_comm%value, gqk%my_nq, myq2glob) 1290 ABI_CHECK(gqk%my_nq > 0, "my_nq == 0") 1291 gqk%my_qstart = myq2glob(1) 1292 ABI_FREE(myq2glob) 1293 1294 ABI_MALLOC(gqk%my_q2ibz, (6, gqk%my_nq)) 1295 ABI_MALLOC(gqk%my_q2bz, (gqk%my_nq)) 1296 1297 do my_iq=1,gqk%my_nq 1298 iq_glob = my_iq + gqk%my_qstart - 1 1299 iq_bz = qglob2bz(iq_glob, spin) 1300 gqk%my_q2ibz(:, my_iq) = qbz2ibz(:, iq_bz) 1301 gqk%my_q2bz(my_iq) = iq_bz 1302 end do 1303 1304 ! Split k-points and transfer symmetry tables 1305 call xmpi_split_block(gqk%glob_nk, gqk%kpt_comm%value, gqk%my_nk, myk2glob) 1306 ABI_CHECK(gqk%my_nk > 0, "my_nk == 0") 1307 gqk%my_kstart = myk2glob(1) 1308 ABI_FREE(myk2glob) 1309 1310 ABI_MALLOC(gqk%my_k2ibz, (6, gqk%my_nk)) 1311 !ABI_MALLOC(gqk%my_k2bz, (gqk%my_nk)) 1312 ABI_MALLOC(gqk%my_kpts, (3, gqk%my_nk)) 1313 ABI_MALLOC(gqk%my_wtk, (gqk%my_nk)) 1314 1315 do my_ik=1,gqk%my_nk 1316 ik_glob = my_ik + gqk%my_kstart - 1 1317 ik_bz = kglob2bz(ik_glob, spin) 1318 gqk%my_k2ibz(:, my_ik) = kbz2ibz(:, ik_bz) 1319 !gqk%my_k2bz(my_ik) = ik_bz 1320 1321 ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik) 1322 trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5, my_ik) 1323 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 1324 tsign_k = 1; if (trev_k == 1) tsign_k = -1 1325 1326 !! symrel^T convention for k 1327 gqk%my_kpts(:, my_ik) = tsign_k * matmul(transpose(gstore%cryst%symrel(:,:,isym_k)), gstore%kibz(:, ik_ibz)) + g0_k 1328 1329 select case(gstore%kzone) 1330 case ("ibz") 1331 gqk%my_wtk(my_ik) = gstore%ebands%wtk(ik_ibz) 1332 case ("bz") 1333 gqk%my_wtk(my_ik) = one / gstore%nkbz 1334 end select 1335 end do 1336 1337 ! Allocate storage for MPI-distributed e-ph matrix elements. 1338 if (with_cplex > 0) then 1339 mem_mb = with_cplex * gqk%my_npert * gqk%nb ** 2 * gqk%my_nq * gqk%my_nk * eight * b2Mb 1340 call wrtout(std_out, sjoin(" Local memory for e-ph matrix elements:", ftoa(mem_mb, fmt="f8.1"), " [Mb] <<< MEM")) 1341 1342 ! The initialization with zero is important as not all the g are computed when we filter in k-space. 1343 ! Abinit postprocessing tools will operate of the full my_g array 1344 ! and we don't want to trigger floating point exceptions. 1345 select case (with_cplex) 1346 case (1) 1347 ABI_MALLOC_OR_DIE(gqk%my_g2, (gqk%my_npert, gqk%nb, gqk%my_nq, gqk%nb, gqk%my_nk), ierr) 1348 gqk%my_g2 = zero 1349 case (2) 1350 ABI_MALLOC_OR_DIE(gqk%my_g, (gqk%my_npert, gqk%nb, gqk%my_nq, gqk%nb, gqk%my_nk), ierr) 1351 gqk%my_g = zero 1352 case default 1353 ABI_ERROR(sjoin("Wrong with_cplex:", itoa(with_cplex))) 1354 end select 1355 1356 ! Allocate storage for MPI-distributed dH/dk matrix elements. 1357 if (any(gstore%with_vk == [1, 2])) then 1358 mem_mb = 3 * gqk%nb * gqk%my_nk * eight * b2Mb 1359 call wrtout(std_out, sjoin(" Memory for diagonal vk_cart_ibz:", ftoa(mem_mb, fmt="f8.1"), " [Mb] <<< MEM")) 1360 ABI_MALLOC_OR_DIE(gqk%vk_cart_ibz, (3, gqk%nb, gstore%nkibz), ierr) 1361 gqk%vk_cart_ibz = zero 1362 end if 1363 1364 if (gstore%with_vk == 2) then 1365 mem_mb = two * 3 * gqk%nb ** 2 * gqk%my_nk * eight * b2Mb 1366 call wrtout(std_out, sjoin(" Memory for vkmat_cart_ibz:", ftoa(mem_mb, fmt="f8.1"), " [Mb] <<< MEM")) 1367 ABI_MALLOC_OR_DIE(gqk%vkmat_cart_ibz, (2, 3, gqk%nb, gqk%nb, gstore%nkibz), ierr) 1368 gqk%vkmat_cart_ibz = zero 1369 end if 1370 end if 1371 1372 end do ! my_is 1373 1374 end subroutine gstore_malloc__
m_gstore/gstore_print [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_print
FUNCTION
Print info on the gstore object.
INPUTS
OUTPUT
SOURCE
1180 subroutine gstore_print(gstore, unit, header, prtvol) 1181 1182 !Arguments ------------------------------------ 1183 class(gstore_t),target,intent(inout) :: gstore 1184 integer,intent(in) :: unit 1185 character(len=*),optional,intent(in) :: header 1186 integer,optional,intent(in) :: prtvol 1187 1188 !Local variables ------------------------------ 1189 integer :: my_is, spin, my_prtvol 1190 type(gqk_t),pointer :: gqk 1191 !character(len=500) :: msg 1192 !type(yamldoc_t) :: ydoc 1193 !---------------------------------------------------------------------- 1194 1195 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 1196 1197 if (present(header)) then 1198 call wrtout(unit, header) 1199 else 1200 call wrtout(unit, " === Gstore parameters ===") 1201 end if 1202 !call ebands_print(gstore%ebands, header="Electron bands", unit=unit, prtvol=my_prtvol) 1203 call wrtout(unit, sjoin(" kzone:", gstore%kzone)) 1204 call wrtout(unit, sjoin(" kfilter:", gstore%kfilter)) 1205 call wrtout(unit, sjoin(" nkibz:", itoa(gstore%nkibz))) 1206 call wrtout(unit, sjoin(" nkbz:", itoa(gstore%nkbz))) 1207 call wrtout(unit, sjoin(" glob_nk_spin:", ltoa(gstore%glob_nk_spin))) 1208 call wrtout(unit, sjoin(" qzone:", gstore%qzone)) 1209 call wrtout(unit, sjoin(" nqibz:", itoa(gstore%nqibz))) 1210 call wrtout(unit, sjoin(" nqbz:", itoa(gstore%nqbz))) 1211 call wrtout(unit, sjoin(" glob_nq_spin:", ltoa(gstore%glob_nq_spin))) 1212 call wrtout(unit, sjoin(" kptopt:", itoa(gstore%ebands%kptopt))) 1213 call wrtout(unit, sjoin(" qptopt:", itoa(gstore%qptopt))) 1214 call wrtout(unit, sjoin(" with_vk:", itoa(gstore%with_vk))) 1215 1216 do my_is=1,gstore%my_nspins 1217 spin = gstore%my_spins(my_is) 1218 gqk => gstore%gqk(my_is) 1219 call wrtout(unit, sjoin(ch10, " === MPI distribution ===")) 1220 call wrtout(unit, sjoin("P Number of CPUs for parallelism over perturbations: ", itoa(gqk%pert_comm%nproc))) 1221 call wrtout(unit, sjoin("P Number of perturbations treated by this CPU: ", itoa(gqk%my_npert))) 1222 call wrtout(unit, sjoin("P Number of CPUs for parallelism over q-points: ", itoa(gqk%qpt_comm%nproc))) 1223 call wrtout(unit, sjoin("P Number of CPUs for parallelism over k-points: ", itoa(gqk%kpt_comm%nproc))) 1224 call wrtout(unit, sjoin(" gqk_cplex:", itoa(gqk%cplex)), pre_newlines=1) 1225 call wrtout(unit, sjoin(" gqk_bstart:", itoa(gqk%bstart))) 1226 call wrtout(unit, sjoin(" gqk_bstop:", itoa(gqk%bstop))) 1227 call wrtout(unit, sjoin(" gqk_nb:", itoa(gqk%nb))) 1228 call wrtout(unit, sjoin(" gqk_my_npert:", itoa(gqk%my_npert))) 1229 call wrtout(unit, sjoin(" gqk_my_nk:", itoa(gqk%my_nk))) 1230 call wrtout(unit, sjoin(" gqk_my_nq:", itoa(gqk%my_nq))) 1231 !if (gqk%cplex == 1) then 1232 ! write(msg,'(a,f8.1,a)')'- Local memory allocated for |g|^2 array: ',ABI_MEM_MB(gqk%my_g2),' [Mb] <<< MEM' 1233 ! call wrtout(unit, msg) !; print *, "my_g2 shape:", shape(gqk%my_g2) 1234 !else 1235 ! write(msg,'(a,f8.1,a)')'- Local memory allocated for g array: ',ABI_MEM_MB(gqk%my_g),' [Mb] <<< MEM' 1236 ! call wrtout(unit, msg) !; print *, "my_g shape:", shape(gqk%my_g) 1237 !end if 1238 end do 1239 1240 end subroutine gstore_print
m_gstore/gstore_set_mpi_grid__ [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_set_mpi_grid__
FUNCTION
INPUTS
OUTPUT
SOURCE
1023 subroutine gstore_set_mpi_grid__(gstore, gstore_cplex, eph_np_pqbks, priority, nproc_spin, comm_spin) 1024 1025 !Arguments ------------------------------------ 1026 !scalars 1027 class(gstore_t),target,intent(inout) :: gstore 1028 integer,intent(in) :: gstore_cplex 1029 integer,intent(in) :: eph_np_pqbks(5), nproc_spin(gstore%nsppol), comm_spin(gstore%nsppol) 1030 character(len=*),intent(in) :: priority 1031 1032 !Local variables------------------------------- 1033 !scalars 1034 integer,parameter :: master = 0, ndims = 3 1035 integer :: spin, my_is, np, my_rank, ierr !, ii 1036 type(gqk_t),pointer :: gqk 1037 character(len=5000) :: msg 1038 character(len=10) :: order 1039 integer :: comm_cart, me_cart, dims(ndims) 1040 logical :: reorder, periods(ndims), keepdim(ndims) 1041 !---------------------------------------------------------------------- 1042 1043 my_rank = xmpi_comm_rank(gstore%comm) 1044 1045 do my_is=1,gstore%my_nspins 1046 spin = gstore%my_spins(my_is) 1047 gqk => gstore%gqk(my_is) 1048 1049 gqk%spin = spin; gqk%natom3 = 3 * gstore%cryst%natom; gqk%cplex = gstore_cplex 1050 ABI_CHECK_IRANGE(gqk%cplex, 1, 2, "gstore_cplex") 1051 1052 ! Compute bstart and band size for this spin. 1053 gqk%bstart = gstore%brange_spin(1, spin) 1054 gqk%bstop = gstore%brange_spin(2, spin) 1055 gqk%nb = gstore%brange_spin(2, spin) - gstore%brange_spin(1, spin) + 1 1056 1057 ! Store global shape of the q/k matrix for this spin. 1058 gqk%glob_nq = gstore%glob_nq_spin(spin) 1059 gqk%glob_nk = gstore%glob_nk_spin(spin) 1060 end do 1061 1062 do my_is=1,gstore%my_nspins 1063 spin = gstore%my_spins(my_is) 1064 gqk => gstore%gqk(my_is) 1065 1066 ! Init for sequential execution. 1067 gqk%qpt_comm%nproc = 1; gqk%kpt_comm%nproc = 1; gqk%pert_comm%nproc = 1 1068 gqk%my_npert = gqk%natom3 1069 1070 if (any(eph_np_pqbks /= 0)) then 1071 ! Use parameters from input file. Need to perform sanity check though. 1072 gqk%pert_comm%nproc = eph_np_pqbks(1) 1073 gqk%qpt_comm%nproc = eph_np_pqbks(2) 1074 ABI_CHECK(eph_np_pqbks(3) == 1, "Band parallelism not implemented in gstore") 1075 gqk%kpt_comm%nproc = eph_np_pqbks(4) 1076 !gqk%spin_comm%nproc = eph_np_pqbks(5) 1077 gqk%my_npert = gqk%natom3 / gqk%pert_comm%nproc 1078 ABI_CHECK(gqk%my_npert > 0, "pert_comm_nproc cannot be greater than 3 * natom.") 1079 ABI_CHECK(mod(gqk%natom3, gqk%pert_comm%nproc) == 0, "pert_comm_nproc must divide 3 * natom.") 1080 1081 else 1082 ! Automatic grid generation (hopefully smart) 1083 ! Keep in mind that in gstore_build, the first loop is over q-points 1084 ! in order to reduce the number of interpolations of the DFPT potentials in q-space 1085 ! hence the q-point parallelism is expected to be more efficient. 1086 ! On the other hand, the k-point parallelism and the perturbation parallelism 1087 ! allow one to reduce the memory requirements associated to the wavefunctions (kpt) and 1088 ! the scattering potentials in the supercell (perturbations). 1089 ! Here we try to optimize performance but it's clear that for large systems the user 1090 ! should specify eph_np_pqbks in the input file. 1091 1092 select case (priority) 1093 case ("q") 1094 order = "1" 1095 case ("k") 1096 order = "2" 1097 case ("kq") 1098 order = "21" 1099 case ("qk") 1100 order = "12" 1101 case default 1102 ABI_ERROR(sjoin("Wrong priority:", priority)) 1103 end select 1104 1105 np = nproc_spin(spin) 1106 call xmpi_distrib_2d(np, order, gqk%glob_nq, gqk%glob_nk, gqk%qpt_comm%nproc, gqk%kpt_comm%nproc, ierr) 1107 ABI_CHECK(ierr == 0, sjoin("Cannot distribute nprocs:", itoa(np), " with priority: ", priority)) 1108 !if (ierr /= 0) call xmpi_distrib_2d_extra(np, gqk%qpt_comm%nproc, gqk%kpt_comm%nproc, gqk%pert_comm%nproc, ierr) 1109 end if 1110 1111 ! Consistency check. 1112 if (gqk%pert_comm%nproc * gqk%qpt_comm%nproc * gqk%kpt_comm%nproc /= nproc_spin(spin)) then 1113 write(msg, "(a,i0,3a, 4(a,1x,i0))") & 1114 "Cannot create 3d Cartesian grid with total nproc: ", nproc_spin(spin), ch10, & 1115 "Idle processes are not supported. The product of the `nproc_*` vars should be equal to nproc.", ch10, & 1116 "qpt_nproc (", gqk%qpt_comm%nproc, ") x kpt_nproc (", gqk%kpt_comm%nproc, ") x kpt_nproc", gqk%pert_comm%nproc, & 1117 ") != ", nproc_spin(spin) 1118 ABI_ERROR(msg) 1119 end if 1120 1121 end do ! my_is 1122 1123 ! For each spin treated by this rank, create 3D cartesian communicator: (q-points, k-points, perturbations) 1124 periods(:) = .False.; reorder = .False. 1125 1126 do my_is=1,gstore%my_nspins 1127 spin = gstore%my_spins(my_is) 1128 gqk => gstore%gqk(my_is) 1129 1130 dims = [gqk%qpt_comm%nproc, gqk%kpt_comm%nproc, gqk%pert_comm%nproc] 1131 1132 ! Note comm_spin(spin) 1133 gqk%grid_comm = xcomm_from_mpi_int(comm_spin(spin)) 1134 1135 #ifdef HAVE_MPI 1136 call MPI_CART_CREATE(gqk%grid_comm, ndims, dims, periods, reorder, comm_cart, ierr) 1137 1138 ! Find the index and coordinates of the current processor 1139 call MPI_COMM_RANK(comm_cart, me_cart, ierr) 1140 call MPI_CART_COORDS(comm_cart, me_cart, ndims, gqk%coords_qkp, ierr) 1141 1142 ! Create communicator for q-points 1143 keepdim = .False.; keepdim(1) = .True.; call gqk%qpt_comm%from_cart_sub(comm_cart, keepdim) 1144 ! Create communicator for k-points 1145 keepdim = .False.; keepdim(2) = .True.; call gqk%kpt_comm%from_cart_sub(comm_cart, keepdim) 1146 ! Create communicator for perturbations. 1147 keepdim = .False.; keepdim(3) = .True.; call gqk%pert_comm%from_cart_sub(comm_cart, keepdim) 1148 ! Create communicator for the (qpt, pert) 2D grid 1149 keepdim = .False.; keepdim(1) = .True.; keepdim(3) = .True.; call gqk%qpt_pert_comm%from_cart_sub(comm_cart, keepdim) 1150 call xmpi_comm_free(comm_cart) 1151 #endif 1152 1153 ! Distribute perturbations inside per_comm using block distribution. 1154 call xmpi_split_block(gqk%natom3, gqk%pert_comm%value, gqk%my_npert, gqk%my_iperts) 1155 end do ! my_is 1156 1157 if (my_rank == master) then 1158 call gstore%print(std_out) 1159 call gstore%print(ab_out) 1160 end if 1161 1162 end subroutine gstore_set_mpi_grid__
m_gstore/gstore_spin2my_is [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
gstore_spin2my_is
FUNCTION
Return the local spin index from the global spin index. 0 if this spin is not treated by this MPI proc.
INPUTS
OUTPUT
SOURCE
1754 integer pure function gstore_spin2my_is(gstore, spin) result(my_is) 1755 1756 !Arguments ------------------------------------ 1757 class(gstore_t),intent(in) :: gstore 1758 integer,intent(in) :: spin 1759 1760 !---------------------------------------------------------------------- 1761 1762 do my_is=1,gstore%my_nspins 1763 if (gstore%my_spins(my_is) == spin) return 1764 end do 1765 my_is = 0 1766 1767 end function gstore_spin2my_is
m_gstore/gstore_t [ Types ]
[ Top ] [ m_gstore ] [ Types ]
NAME
gstore_t
FUNCTION
This object stores:
- pointers to the cyrstalline structure, the KS bands, the IFCs.
- arrays that do not depend on the spin such as the IBZ and weights for k/q-points.
- metadata such as kzone, qzone and kfilter that are needed to interpret
the storage mode used for the g(k, q)
NB: the e-ph matrix element are stored in gstore%qqk(my_is)
where my_is counts the number of spins treated by this MPI processor.
NOTES
SOURCE
326 type, public :: gstore_t 327 328 integer :: nsppol 329 ! Number of independent spin polarizations. 330 331 integer :: my_nspins = 0 332 ! Number of collinear spins treated by this MPI rank 333 334 integer :: nkibz = -1, nqibz = -1 335 ! Number of k/q points in the IBZ. 336 337 integer :: nkbz = -1, nqbz = -1 338 ! Number of k/q points in the BZ. 339 340 integer :: comm 341 ! Global communicator 342 ! Inherited by the caller thus we don't free it in gstore_free. 343 344 integer :: with_vk = 0 345 ! 0 if group velocities should not be computed 346 ! 1 to compute diagonal terms only 347 ! 2 to compute diagonal and off-diagonal terms 348 349 integer :: qptopt = -1 350 ! qptopt=option for the generation of q points (defines whether spatial symmetries and/or time-reversal can be used) 351 352 character(len=fnlen) :: path = " " 353 ! Path to the nc file associated to the gstore 354 355 character(len=fnlen) :: wfk0_path = " " 356 357 character(len=fnlen) :: kzone = " ", qzone = " " 358 ! Specifies whether k- or q-points are in the BZ or in the IBZ. 359 ! Possible values are "ibz" or "bz". 360 ! Note that the combination ("ibz", "ibz") is not allowed 361 362 character(len=fnlen) :: kfilter = "none" 363 ! Specifies the tecnique used to filter k-points. 364 ! Possible values: "none", "fs_tetra", "erange". 365 366 real(dp),allocatable :: erange_spin(:, :) 367 ! (2, nsppol) 368 ! Energy window. zero if not used. Requires kfilter == "erange" 369 370 type(crystal_t), pointer :: cryst => null() 371 372 type(ebands_t), pointer :: ebands => null() 373 374 type(ifc_type), pointer :: ifc => null() 375 376 type(krank_t) :: krank_ibz !, qrank 377 ! Object used to find k-points in the IBZ and map BZ to IBZ. 378 379 integer,allocatable :: my_spins(:) 380 ! (%my_nspins) 381 ! Indirect table giving the spin indices treated by this MPI rank. 382 ! Used only in the collinear case with nsppol = 2. 383 384 integer,allocatable :: brange_spin(:, :) 385 ! (2, nsppol) 386 387 !integer :: max_nb = -1 388 ! Max number of bands over spin 389 390 integer,allocatable :: glob_nk_spin(:), glob_nq_spin(:) 391 ! (nsppol) 392 393 real(dp), contiguous, pointer :: kibz(:,:) 394 ! k-points in the IBZ. Points to ebands%kptns 395 ! (3, nkibz) 396 397 real(dp),allocatable :: delta_ef_kibz_spin(:,:,:) 398 ! (nb, gstore%nkibz, nsppol)) 399 ! Tetrahedron weights at eF in the IBZ. 400 401 real(dp), allocatable :: qibz(:,:) 402 ! (3, nqibz) 403 ! q-points in the IBZ 404 405 real(dp), allocatable :: wtq(:) 406 ! (nqibz) 407 ! q-points weights in the IBZ 408 409 !integer :: qptrlatt(3, 3) = -1 ! kptrlatt(3, 3) = -1, 410 ! k-mesh and q-mesh 411 412 !real(dp),allocatable :: kshift(:, :), qshift(:, :) 413 ! k/q-mesh shift (well, q-mesh is usually gamma-centered) 414 415 type(gqk_t), allocatable :: gqk(:) 416 ! (my_nspins) 417 ! Datastructure storing e-ph matrix elements. 418 419 !type(htetra_t) :: ktetra 420 ! Used to evaluate integrals in k-space with the tetrahedron method. 421 422 !type(htetra_t) :: qtetra 423 ! Used to evaluate integrals in q-space with the tetrahedron method. 424 425 contains 426 427 procedure :: fill_bks_mask => gstore_fill_bks_mask 428 ! Fill the table used to read (b, k, s) wavefunctions from the WFK file 429 ! keeping into account the distribution of the e-ph matrix elements. 430 431 procedure :: get_mpw_gmax => gstore_get_mpw_gmax 432 433 procedure :: spin2my_is => gstore_spin2my_is 434 ! Return the local spin index from the global spin index. 435 ! 0 if this spin is not treated by this MPI proc. 436 437 procedure :: free => gstore_free 438 ! Free memory 439 440 procedure :: print => gstore_print 441 ! Print info on the object 442 443 procedure, private :: distribute_spins__ => gstore_distribute_spins 444 ! Distribute spins (nsppol = 2) and create indirect mapping to spin index. 445 446 procedure, private :: set_mpi_grid__ => gstore_set_mpi_grid__ 447 ! Set the MPI cartesian grid 448 449 procedure, private :: malloc__ => gstore_malloc__ 450 ! Allocate local buffers once the MPI grid has been initialized. 451 452 procedure, private :: filter_fs_tetra__ => gstore_filter_fs_tetra__ 453 ! Select k-points on the FS using the tetrahedron method 454 455 procedure, private :: filter_erange__ => gstore_filter_erange__ 456 ! Select k-points inside an energy window 457 458 procedure :: compute => gstore_compute 459 ! Compute e-ph matrix elements. 460 461 procedure :: calc_my_phonons => gstore_calc_my_phonons 462 ! Helper function to compute ph quantities for all q-points treated by the MPI proc. 463 464 procedure :: get_lambda_iso_iw => gstore_get_lambda_iso_iw 465 ! Compute isotropic lamdda(iw) along the imaginary axis. 466 467 procedure :: get_a2fw => gstore_get_a2fw 468 ! Compute Eliashberg function a^2F(w). 469 470 procedure :: from_ncpath => gstore_from_ncpath 471 ! Reconstruct object from netcdf file 472 473 procedure :: init => gstore_init 474 ! Build object 475 476 end type gstore_t
m_gstore/priority_from_eph_task [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
priority_from_eph_task
FUNCTION
INPUTS
OUTPUT
SOURCE
921 subroutine priority_from_eph_task(eph_task, priority) 922 923 !Arguments ------------------------------------ 924 integer,intent(in) :: eph_task 925 character(len=*),intent(out) :: priority 926 !---------------------------------------------------------------------- 927 928 select case (eph_task) 929 case (11) 930 priority = "qk" 931 case (12, -12) 932 priority = "q" 933 case (14) 934 priority = "kq" 935 case default 936 ABI_ERROR(sjoin("Please register default priority for eph_task:", itoa(eph_task))) 937 end select 938 939 end subroutine priority_from_eph_task
m_gstore/recompute_select_qbz_spin [ Functions ]
[ Top ] [ m_gstore ] [ Functions ]
NAME
recompute_select_qbz_spin
FUNCTION
Recompute select_qbz_spin table after the filtering on the k-points.
INPUTS
OUTPUT
SOURCE
1657 subroutine recompute_select_qbz_spin(gstore, qbz, qbz2ibz, qibz2bz, kbz, kibz, kbz2ibz, kibz2bz, & 1658 select_kbz_spin, select_qbz_spin) 1659 1660 !Arguments ------------------------------------ 1661 !scalars 1662 class(gstore_t),target,intent(inout) :: gstore 1663 !arrays 1664 real(dp),intent(in) :: qbz(3, gstore%nqbz) 1665 integer,intent(in) :: qbz2ibz(6,gstore%nqbz), qibz2bz(gstore%nqibz) 1666 integer,intent(in) :: kbz2ibz(6,gstore%nkbz), kibz2bz(gstore%nkibz) 1667 real(dp),target,intent(in) :: kibz(3, gstore%nkibz), kbz(3, gstore%nkbz) 1668 integer,intent(in) :: select_kbz_spin(gstore%nkbz, gstore%nsppol) 1669 integer,intent(out) :: select_qbz_spin(gstore%nqbz, gstore%nsppol) 1670 1671 !Local variables------------------------------- 1672 !scalars 1673 integer :: all_nproc, my_rank, ierr 1674 integer :: ii, iq_bz, iq_ibz, ikq_ibz, ikq_bz, len_kpts_ptr, ebands_timrev 1675 !arrays 1676 integer,allocatable :: map_kq(:,:) 1677 real(dp):: qpt(3) 1678 real(dp),contiguous, pointer :: kpts_ptr(:,:) 1679 1680 ! ************************************************************************* 1681 1682 ABI_UNUSED(kbz2ibz) 1683 ABI_UNUSED(qbz2ibz) 1684 1685 all_nproc = xmpi_comm_size(gstore%comm); my_rank = xmpi_comm_rank(gstore%comm) 1686 ebands_timrev = kpts_timrev_from_kptopt(gstore%ebands%kptopt) 1687 1688 select_qbz_spin = 0 1689 1690 if (gstore%kzone == "ibz") kpts_ptr => kibz 1691 if (gstore%kzone == "bz") kpts_ptr => kbz 1692 len_kpts_ptr = size(kpts_ptr, dim=2) 1693 ABI_MALLOC(map_kq, (6, len_kpts_ptr)) 1694 1695 select case (gstore%qzone) 1696 case ("ibz") 1697 do iq_ibz=1,gstore%nqibz 1698 if (mod(iq_ibz, all_nproc) /= my_rank) cycle ! MPI parallelism. 1699 qpt = gstore%qibz(:, iq_ibz) 1700 iq_bz = qibz2bz(iq_ibz) 1701 ! k + q_ibz --> k IBZ --> k BZ 1702 1703 if (kpts_map("symrel", ebands_timrev, gstore%cryst, gstore%krank_ibz, len_kpts_ptr, kpts_ptr, map_kq, qpt=qpt) /= 0) then 1704 ABI_ERROR("Cannot map k+q to IBZ!") 1705 end if 1706 1707 do ii=1,len_kpts_ptr 1708 ikq_ibz = map_kq(1, ii) 1709 ikq_bz = kibz2bz(ikq_ibz) 1710 select_qbz_spin(iq_bz, :) = select_qbz_spin(iq_bz, :) + select_kbz_spin(ikq_bz, :) 1711 end do 1712 end do ! iq_ibz 1713 1714 case ("bz") 1715 do iq_bz=1,gstore%nqbz 1716 if (mod(iq_bz, all_nproc) /= my_rank) cycle ! MPI parallelism. 1717 qpt = qbz(:, iq_bz) 1718 !iq_ibz = qbz2ibz(1, iq_bz) 1719 ! k + q_bz --> k IBZ --> k BZ 1720 1721 if (kpts_map("symrel", ebands_timrev, gstore%cryst, gstore%krank_ibz, len_kpts_ptr, kpts_ptr, map_kq, qpt=qpt) /= 0) then 1722 ABI_ERROR("Cannot map k+q to IBZ!") 1723 end if 1724 1725 do ii=1,len_kpts_ptr 1726 ikq_ibz = map_kq(1, ii) 1727 ikq_bz = kibz2bz(ikq_ibz) 1728 select_qbz_spin(iq_bz, :) = select_qbz_spin(iq_bz, :) + select_kbz_spin(ikq_bz, :) 1729 end do 1730 end do 1731 end select 1732 1733 call xmpi_sum(select_qbz_spin, gstore%comm, ierr) 1734 1735 ABI_FREE(map_kq) 1736 1737 end subroutine recompute_select_qbz_spin