Line data Source code
1 : !!****m* ABINIT/m_eph_path
2 : !! NAME
3 : !! m_eph_path
4 : !!
5 : !! FUNCTION
6 : !! Compute e-ph matrix elements g(k,q) along an arbitrary path either in k- or q-space.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_eph_path
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_xmpi
27 : use m_mpinfo
28 : use m_errors
29 : use m_copy
30 : use m_hamiltonian
31 : use m_pawcprj
32 : use m_ephtk
33 : use netcdf
34 : use m_nctk
35 :
36 : use defs_abitypes, only : mpi_type
37 : use defs_datatypes, only : pseudopotential_type
38 : use m_time, only : cwtime, cwtime_report, timab, sec2str
39 : use m_fstrings, only : itoa, ftoa, sjoin, ktoa, ltoa, strcat
40 : use m_cgtools, only : cg_zdotc
41 : use m_crystal, only : crystal_t
42 : use m_dtset, only : dataset_type
43 : use m_dtfil, only : datafiles_type
44 : use m_ebands, only : ebands_t
45 : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack
46 : use m_pawang, only : pawang_type
47 : use m_pawrad, only : pawrad_type
48 : use m_pawtab, only : pawtab_type
49 : use m_pawfgr, only : pawfgr_type
50 : use m_cgwf, only : nscf_t
51 : use m_bz_mesh, only : kpath_t
52 : use m_wfd, only : u0_cache_t
53 : use m_ifc, only : ifc_type
54 : use m_dvdb, only : dvdb_t
55 : use m_mlwfovlp, only : wan_t
56 :
57 : implicit none
58 :
59 : private
60 : !!***
61 :
62 : public :: eph_path_run
63 :
64 : contains !=====================================================
65 : !!***
66 :
67 : !----------------------------------------------------------------------
68 :
69 : !!****f* m_eph_path/eph_path_run
70 : !! NAME
71 : !! eph_path_run
72 : !!
73 : !! FUNCTION
74 : !! Compute e-ph matrix elements g(k,q) along an arbitrary path either in k- or q-space.
75 : !! Wavefunctions at k and k+q are computed non-self-consistently by invoking the CG eigensolver
76 : !! starting from the GS potential read from file.
77 : !! The DFPT potentials at q are usually obtained via Fourier interpolation, but it is also possible
78 : !! to use fully ab-initio potentials provided the DVDB file contains all the q-points along the path.
79 : !! This requires performing DFPT calculations for all the q-points, and then merging
80 : !! all the POT1 files with the mrgdv utility.
81 : !!
82 : !! INPUTS
83 : !! dtfil<datafiles_type>=Variables related to files.
84 : !! dtset<dataset_type>=All input variables for this dataset.
85 : !! cryst<crystal_t>=crystal structure parameters
86 : !! wfk_ebands: electron bands from the input WFK file (used to propagate nelect and fermie to GPATH.nc)
87 : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
88 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
89 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
90 : !! pawang<pawang_type)>=PAW angular mesh and related data.
91 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
92 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
93 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
94 : !! comm=MPI communicator.
95 : !!
96 : !! OUTPUT
97 : !! Only writing.
98 : !!
99 : !! SOURCE
100 :
101 2 : subroutine eph_path_run(dtfil, dtset, cryst, wfk_ebands, dvdb, ifc, pawfgr, pawang, pawrad, pawtab, psps, comm)
102 :
103 : !Arguments ------------------------------------
104 : !scalars
105 : type(datafiles_type),intent(in) :: dtfil
106 : type(dataset_type),intent(in) :: dtset
107 : type(crystal_t),intent(in) :: cryst
108 : type(ebands_t),intent(in) :: wfk_ebands
109 : type(dvdb_t),intent(inout) :: dvdb
110 : type(pawang_type),intent(in) :: pawang
111 : type(pseudopotential_type),intent(in) :: psps
112 : type(pawfgr_type),intent(in) :: pawfgr
113 : type(ifc_type),intent(in) :: ifc
114 : integer,intent(in) :: comm
115 : !arrays
116 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
117 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
118 :
119 : !Local variables ------------------------------
120 : !scalars
121 : integer,parameter :: istwfk_1 = 1, tim_getgh1c = 1, berryopt0 = 0, master = 0, ndims = 3, ndat1 = 1
122 : integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1, nu
123 : integer :: spin, iq, ik, nk_path, nq_path, ierr, npw_k, npw_kq, my_rank, nprocs, n1, n2, n3, n4, n5, n6, cplex
124 : integer :: natom, natom3, nsppol, nspden, nspinor, qptopt, comm_cart, me_cart
125 : integer :: nfft,nfftf,mgfft,mgfftf, my_npert, my_ip, idir, ipert, ipc, ncerr, ncid, my_nkpath, my_nqpath
126 : integer :: in_k, im_kq, my_is, my_ik, my_iq, nband, nb_in_g, ii, band_n, band_m, bstart, bstop, my_nspins, np, tot_nscf_ierr
127 : integer :: my_pert_start_glob, nwan_glob, iglob
128 : real(dp) :: cpu_all,wall_all,gflops_all, eig0nk, eshift
129 : logical :: qq_is_gamma, need_ftinterp, gen_eigenpb, use_cg_k, use_cg_kq, use_cache, has_gwan
130 2 : type(gs_hamiltonian_type) :: gs_ham_k, gs_ham_kq
131 2 : type(rf_hamiltonian_type) :: rf_ham_kq
132 2 : type(nscf_t) :: nscf
133 2 : type(kpath_t) :: qpath, kpath
134 : type(xcomm_t),target :: kpt_comm, qpt_comm, pert_comm
135 14 : type(u0_cache_t) :: ucache_kq, ucache_k
136 14 : type(wan_t) :: wan
137 : character(len=fnlen) :: gpath_path
138 : character(len=5000) :: msg
139 : character(len=10) :: priority
140 : !!arrays
141 : integer :: units(2), ngfft(18),ngfftf(18), coords_spin(ndims), dims(ndims)
142 2 : integer,allocatable :: kg_k(:,:), kg_kq(:,:), qmap_symrec(:,:), my_ik_inds(:), my_iq_inds(:), my_spins(:), my_iperts(:)
143 2 : integer,allocatable :: pert_table(:,:), my_pinfo(:,:)
144 4 : real(dp) :: kk(3), qq(3), kq(3), phfreqs(3*cryst%natom), phfreqs_ev(3*cryst%natom), fake_path(3,2)
145 2 : real(dp),allocatable :: grad_berry(:,:), kinpw_k(:), kinpw_kq(:)
146 2 : real(dp),allocatable :: cg_k(:,:,:), cg_kq(:,:,:), gsc_k(:,:,:), gsc_kq(:,:,:),eig_k(:), eig_kq(:)
147 2 : real(dp),allocatable :: eig_k_wan(:), eig_kq_wan(:,:)
148 4 : real(dp),allocatable :: v1scf(:,:,:,:), vlocal1(:,:,:,:), vlocal(:,:,:,:), gkq_atm(:,:,:,:), gkq_nu(:,:,:,:), gkq2_nu(:,:,:)
149 2 : real(dp),allocatable :: gkq_atm_wan(:,:,:,:), gkq_nu_wan(:,:,:,:), gkq2_nu_wan(:,:,:)
150 2 : complex(dp),allocatable :: g_atm_wan_local(:,:,:,:)
151 2 : real(dp),allocatable :: gvnlx1(:,:), gs1c(:,:), h1_kets_kq(:,:,:), displ_cart(:,:,:,:),displ_red_qq(:,:,:,:)
152 4 : real(dp),allocatable :: kpg_k(:,:), ph3d_k(:,:,:), ffnl_k(:,:,:,:), vlocal_k(:,:,:,:)
153 2 : real(dp),allocatable :: kpg_kq(:,:), ph3d_kq(:,:,:), ffnl_kq(:,:,:,:), vlocal_kq(:,:,:,:), real_vec(:)
154 : logical :: reorder, periods(ndims), keepdim(ndims)
155 2 : type(pawcprj_type),allocatable :: cwaveprj0(:,:)
156 2 : type(xcomm_t),allocatable :: comm_my_is(:)
157 : !************************************************************************
158 :
159 2 : if (psps%usepaw == 1) then
160 0 : ABI_ERROR("PAW not implemented")
161 : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
162 : end if
163 :
164 2 : call cwtime(cpu_all, wall_all, gflops_all, "start")
165 2 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
166 6 : units = [std_out, ab_out]
167 :
168 : ! Copy important dimensions.
169 2 : natom = cryst%natom; natom3 = 3 * natom; nsppol = dtset%nsppol; nspinor = dtset%nspinor; nspden = dtset%nspden
170 :
171 : ! Build (k/q)-path. NB: Input variables have been already checked for consistency in chkinp.
172 1 : select case (dtset%eph_fix_korq)
173 : case ("k")
174 1 : call qpath%init(dtset%ph_qpath(:,1:dtset%ph_nqpath), cryst%gprimd, dtset%ph_ndivsm)
175 1 : nq_path = qpath%npts
176 1 : call qpath%print(units, header=sjoin("q-point path for g(k,q) with fixed k:", ktoa(dtset%eph_fix_wavevec)), prtvol=dtset%prtvol)
177 7 : fake_path(:,1) = dtset%eph_fix_wavevec; fake_path(:,2) = dtset%eph_fix_wavevec + one
178 1 : call kpath%init(fake_path, cryst%gprimd, 0)
179 1 : nk_path = 1
180 :
181 : case ("q")
182 1 : call kpath%init(dtset%kptbounds(:,1:dtset%nkpath), cryst%gprimd, dtset%ndivsm)
183 1 : nk_path = kpath%npts
184 1 : call kpath%print(units, header=sjoin("k-point path for g(k,q) with fixed q:", ktoa(dtset%eph_fix_wavevec)), prtvol=dtset%prtvol)
185 7 : fake_path(:,1) = dtset%eph_fix_wavevec; fake_path(:,2) = dtset%eph_fix_wavevec + one
186 1 : call qpath%init(fake_path, cryst%gprimd, 0)
187 1 : nq_path = 1
188 :
189 : case default
190 2 : ABI_ERROR(sjoin("Invalid value of eph_fix_korq:", dtset%eph_fix_korq))
191 : end select
192 :
193 : ! Define band range and nb_in_g from eph_path_brange.
194 2 : nband = dtset%mband; bstart = dtset%eph_path_brange(1); bstop = dtset%eph_path_brange(2)
195 :
196 : ! Check whether e-ph matrix elements should also be computed via Wannier interpolation from
197 : ! GWAN.nc, for comparison against the ab-initio results computed below.
198 2 : has_gwan = (dtfil%filgwanin /= ABI_NOFILE)
199 2 : nwan_glob = 0
200 2 : if (has_gwan) then
201 0 : ABI_CHECK(dtfil%filabiwanin /= ABI_NOFILE, "getgwan_filepath requires getabiwan_filepath to be given as well when eph_task=18.")
202 0 : ABI_CHECK(dtset%nsppol == 1, "has_gwan with nsppol == 2 is not yet supported in eph_path_run.")
203 :
204 : ! Cheap, non-collective read of ABIWAN.nc metadata, purely to validate/auto-set eph_path_brange
205 : ! (mirrors gstore_init's has_abiwan branch, m_gstore.F90).
206 0 : call wan%from_abiwan(dtfil%filabiwanin, 1, dtset%nsppol, .False., "", xmpi_comm_self)
207 0 : if (bstart <= 0 .and. bstop <= 0) then
208 0 : bstart = wan%bmin; bstop = wan%bmax
209 0 : call wrtout(units, sjoin(" has_gwan: auto-setting eph_path_brange to Wannier band window:", ltoa([bstart,bstop])))
210 : else
211 0 : if (bstart <= 0) bstart = 1
212 0 : if (bstop <= 0) bstop = nband
213 : msg = sjoin("eph_path_brange:", ltoa([bstart,bstop]), "must equal [wan%bmin, wan%bmax]:", ltoa([wan%bmin, wan%bmax]), &
214 0 : "when getgwan_filepath is given (the gauge-invariant trace comparison needs the same electronic subspace).")
215 0 : ABI_CHECK(bstart == wan%bmin .and. bstop == wan%bmax, msg)
216 : end if
217 0 : nwan_glob = wan%nwan
218 0 : call wan%free()
219 : else
220 2 : if (bstart <= 0) bstart = 1
221 2 : if (bstop <= 0) bstop = nband
222 : end if
223 :
224 2 : nb_in_g = bstop - bstart + 1
225 :
226 : ! The values of eph_path_brange must be validated at this level!
227 2 : ABI_CHECK_IRANGE(bstart, 1, nband, "Wrong eph_path_brange(1)")
228 2 : ABI_CHECK_IRANGE(bstop, 1, nband, "Wrong eph_path_brange(2)")
229 2 : ABI_CHECK_IGEQ(bstop, bstart, "eph_path_brange(2) < eph_path_brange(1)")
230 6 : call wrtout(units, sjoin(" Computing g(k,q) with eph_path_brange:", ltoa([bstart, bstop])))
231 :
232 : ! Distribute spins inside input comm.
233 2 : call xmpi_split_nsppol(comm, nsppol, my_nspins, my_spins, comm_my_is)
234 :
235 : ! ==================
236 : ! MPI cartesian grid
237 : ! ==================
238 2 : my_npert = natom3
239 4 : do my_is=1,my_nspins
240 2 : spin = my_spins(my_is)
241 2 : np = comm_my_is(my_is)%nproc; priority = "12"
242 :
243 12 : if (any(dtset%eph_np_pqbks /= 0)) then
244 : ! Take MPI grid from input.
245 0 : pert_comm%nproc = dtset%eph_np_pqbks(1)
246 0 : qpt_comm%nproc = dtset%eph_np_pqbks(2)
247 0 : ABI_CHECK_IEQ(dtset%eph_np_pqbks(3), 1, "Band parallelism not implemented in eph_path")
248 0 : kpt_comm%nproc = dtset%eph_np_pqbks(4)
249 :
250 : else
251 : ! Automatic MPI grid generation.
252 2 : if (nk_path == 1) then
253 1 : call xmpi_distrib_2d(np, priority, nq_path, natom3, qpt_comm%nproc, pert_comm%nproc, ierr)
254 : else
255 1 : call xmpi_distrib_2d(np, priority, nk_path, natom3, kpt_comm%nproc, pert_comm%nproc, ierr)
256 : end if
257 2 : ABI_CHECK(ierr == 0, sjoin("Cannot distribute nprocs:", itoa(np), " with priority: ", priority, ". Decrease MPI nprocs"))
258 : end if
259 :
260 : ! Consistency check
261 2 : write(msg, "(a,i2,a,3(i0,1x))")"P Cartesian grid for spin", spin, ": (pert_comm%nproc, qpt_comm%nproc, kpt_comm%nproc) = ", &
262 4 : pert_comm%nproc, qpt_comm%nproc, kpt_comm%nproc
263 2 : call wrtout(units, msg)
264 :
265 2 : if (pert_comm%nproc * qpt_comm%nproc * kpt_comm%nproc /= np) then
266 : write(msg, "(a,i0,3a, 4(a,1x,i0))") &
267 0 : "Cannot create Cartesian grid with nproc(spin): ", np, ch10, &
268 0 : "Idle processes are not supported. The product of the `nproc_*` vars should be equal to nproc.", ch10, &
269 0 : "qpt_nproc (", qpt_comm%nproc, ") x kpt_nproc (", kpt_comm%nproc, ") x pert_nproc", pert_comm%nproc, &
270 0 : ") != ", np
271 0 : ABI_ERROR(msg)
272 : end if
273 :
274 : ! For each spin treated by this rank, create MPI cartesian communicator of rank ndims.
275 2 : periods(:) = .False.; reorder = .False.
276 8 : dims = [pert_comm%nproc, qpt_comm%nproc, kpt_comm%nproc]
277 :
278 : #ifdef HAVE_MPI
279 2 : call MPI_CART_CREATE(comm_my_is(my_is)%value, ndims, dims, periods, reorder, comm_cart, ierr)
280 : ! Find the index and coordinates of the current processor
281 2 : call MPI_COMM_RANK(comm_cart, me_cart, ierr)
282 2 : call MPI_CART_COORDS(comm_cart, me_cart, ndims, coords_spin, ierr)
283 :
284 : ! Communicator for q-points in g(k,q)
285 2 : keepdim = .False.; keepdim(1) = .True.; call pert_comm%from_cart_sub(comm_cart, keepdim)
286 2 : keepdim = .False.; keepdim(2) = .True.; call qpt_comm%from_cart_sub(comm_cart, keepdim)
287 2 : keepdim = .False.; keepdim(3) = .True.; call kpt_comm%from_cart_sub(comm_cart, keepdim)
288 4 : call xmpi_comm_free(comm_cart)
289 : #endif
290 : end do ! my_is
291 :
292 : ! Distribute k-points (q-points) inside kpt_comm (qpt_comm) using block distribution.
293 2 : call xmpi_split_block(nk_path, kpt_comm%value, my_nkpath, my_ik_inds)
294 2 : call xmpi_split_block(nq_path, qpt_comm%value, my_nqpath, my_iq_inds)
295 2 : call xmpi_split_block(natom3, pert_comm%value, my_npert, my_iperts)
296 2 : my_pert_start_glob = my_iperts(1)
297 2 : ABI_FREE(my_iperts)
298 :
299 : ! Idle processors are not supported (tested).
300 2 : ABI_CHECK_IGEQ(my_nkpath, 1, "Too many procs for k-point parallelism.")
301 2 : ABI_CHECK_IGEQ(my_nqpath, 1, "Too many procs for q-point parallelism.")
302 2 : ABI_CHECK_IGEQ(my_npert, 1, "Too many procs for perturbation parallelism.")
303 :
304 2 : if (pert_comm%nproc > 1) then
305 : ! Build table with list of perturbations treated by this CPU inside pert_comm
306 0 : call ephtk_set_pertables(cryst%natom, my_npert, pert_table, my_pinfo, pert_comm%value)
307 : ! Activate parallelism over perturbations
308 0 : call dvdb%set_pert_distrib(my_npert, natom3, my_pinfo, pert_table, pert_comm%value)
309 0 : ABI_FREE(my_pinfo)
310 0 : ABI_FREE(pert_table)
311 0 : ABI_WARNING("Parallelism over perturbations should be tested!")
312 : end if
313 :
314 : ! Load KS potential from file.
315 2 : call nscf%init(dtset, dtfil, cryst, comm)
316 :
317 : ! FFT meshes (taken from the GS POT file)
318 74 : ngfft = nscf%ngfft; ngfftf = nscf%ngfftf
319 8 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
320 8 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
321 2 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3); n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6)
322 :
323 : ! Open the DVDB file and make sure we have POT1 files.
324 2 : call dvdb%open_read(ngfftf, xmpi_comm_self)
325 2 : ABI_CHECK(dvdb%has_fields("pot1", msg), msg)
326 :
327 : ! Check if all the q-points are present in the DVDB.
328 2 : qptopt = dtset%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt
329 2 : call dvdb%need_ftinterp(nq_path, qpath%points, qptopt, qmap_symrec, need_ftinterp)
330 :
331 2 : if (.not. need_ftinterp .and. dtset%eph_use_ftinterp /= 0) then
332 0 : ABI_COMMENT("Enforcing FT interpolation for q-points even if it's not strictly needed.")
333 0 : need_ftinterp = .True.
334 : end if
335 :
336 2 : if (need_ftinterp) then
337 2 : call wrtout(units, " Cannot find all q-points in the DVDB --> Activating Fourier interpolation.")
338 2 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt, 1, dtset%ddb_shiftq, nfftf, ngfftf, xmpi_comm_self)
339 : else
340 0 : call wrtout(units, " DVDB file contains all q-points along the path --> Reading DFPT potentials from file.")
341 : end if
342 :
343 : ! Prepare call to getgh1c
344 2 : usevnl = 0
345 2 : optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
346 2 : optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
347 2 : opt_gvnlx1 = 0 ! gvnlx1 is output
348 2 : usecprj = 0
349 :
350 2 : ABI_MALLOC(gvnlx1, (2, usevnl))
351 2 : ABI_MALLOC(grad_berry, (2, nspinor*(berryopt0/4)))
352 4 : ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
353 8 : ABI_MALLOC(displ_cart, (2, 3, cryst%natom, natom3))
354 6 : ABI_MALLOC(displ_red_qq, (2, 3, cryst%natom, natom3))
355 10 : ABI_MALLOC(gkq_atm, (2, nb_in_g, nb_in_g, natom3))
356 8 : ABI_MALLOC(gkq_nu, (2, nb_in_g, nb_in_g, natom3))
357 10 : ABI_MALLOC(gkq2_nu, (nb_in_g, nb_in_g, natom3))
358 :
359 2 : if (has_gwan) then
360 0 : ABI_MALLOC(gkq_atm_wan, (2, nwan_glob, nwan_glob, natom3))
361 0 : ABI_MALLOC(gkq_nu_wan, (2, nwan_glob, nwan_glob, natom3))
362 0 : ABI_MALLOC(gkq2_nu_wan, (nwan_glob, nwan_glob, natom3))
363 0 : ABI_MALLOC(eig_k_wan, (nwan_glob))
364 0 : ABI_MALLOC(eig_kq_wan, (nwan_glob, 1))
365 : end if
366 :
367 : ! Master writes metadata to GPATH file.
368 2 : gpath_path = strcat(dtfil%filnam_ds(4), "_GPATH.nc")
369 :
370 2 : if (my_rank == master) then
371 2 : NCF_CHECK(nctk_open_create(ncid, gpath_path, xmpi_comm_self))
372 : ! Add crystalline structure.
373 2 : NCF_CHECK(cryst%ncwrite(ncid))
374 :
375 : ! Write dimensions.
376 : ncerr = nctk_def_dims(ncid, [ &
377 : nctkdim_t("nspinor", nspinor), nctkdim_t("nspden", nspden), nctkdim_t("nsppol", nsppol), &
378 : nctkdim_t("nband", nband), nctkdim_t("nb_in_g", nb_in_g), &
379 : nctkdim_t("nq_path", nq_path), nctkdim_t("nk_path", nk_path), &
380 : nctkdim_t("natom", cryst%natom), nctkdim_t("natom3", natom3), nctkdim_t("number_of_phonon_modes", natom3) &
381 22 : ], defmode=.True.)
382 2 : NCF_CHECK(ncerr)
383 :
384 : ! integer scalars
385 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
386 : "bstart", "bstop", "dvdb_add_lr", "used_ftinterp" &
387 10 : ])
388 2 : NCF_CHECK(ncerr)
389 :
390 : ! double precision scalars
391 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
392 : "nelect", "fermie" &
393 6 : ])
394 2 : NCF_CHECK(ncerr)
395 :
396 : ncerr = nctk_def_arrays(ncid, [ &
397 : nctkarr_t("kpoints", "dp", "three, nk_path"), &
398 : nctkarr_t("qpoints", "dp", "three, nq_path"), &
399 : nctkarr_t("qweights", "dp", "nq_path"), &
400 : nctkarr_t("all_eigens_k", "dp", "nband, nk_path, nsppol"), &
401 : nctkarr_t("all_eigens_kq", "dp", "nband, nq_path, nsppol"), &
402 : nctkarr_t("gkq2_nu", "dp", "nb_in_g, nb_in_g, natom3, nq_path, nk_path, nsppol"), &
403 : nctkarr_t("phfreqs", "dp", "natom3, nq_path"), &
404 : nctkarr_t("phdispl_cart", "dp", "two, three, natom, natom3, nq_path"), &
405 : nctkarr_t("eph_fix_korq", "c", "one"), &
406 : nctkarr_t("eph_fix_wavevec", "dp", "three") &
407 22 : ])
408 2 : NCF_CHECK(ncerr)
409 :
410 2 : NCF_CHECK(nf90_def_var_fill(ncid, vid("gkq2_nu"), NF90_FILL, -one))
411 2 : NCF_CHECK(nf90_def_var_fill(ncid, vid("phfreqs"), NF90_FILL, -one))
412 :
413 2 : if (has_gwan) then
414 0 : ncerr = nctk_def_dims(ncid, [nctkdim_t("nwan", nwan_glob)], defmode=.True.)
415 0 : NCF_CHECK(ncerr)
416 : ncerr = nctk_def_arrays(ncid, [ &
417 : ! Use the same path convention as the ab-initio eigenvalues: k varies when q
418 : ! is fixed, whereas k+q varies with q when k is fixed.
419 : nctkarr_t("all_eigens_wan_k", "dp", "nwan, nk_path, nsppol"), &
420 : nctkarr_t("all_eigens_wan_kq", "dp", "nwan, nq_path, nsppol"), &
421 : nctkarr_t("gkq2_nu_wan", "dp", "nwan, nwan, natom3, nq_path, nk_path, nsppol") &
422 0 : ])
423 0 : NCF_CHECK(ncerr)
424 0 : NCF_CHECK(nf90_def_var_fill(ncid, vid("gkq2_nu_wan"), NF90_FILL, -one))
425 : end if
426 :
427 : ! Write data.
428 2 : NCF_CHECK(nctk_set_datamode(ncid))
429 2 : ii = merge(1, 0, need_ftinterp)
430 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
431 : "bstart", "bstop", "dvdb_add_lr", "used_ftinterp"], &
432 : [bstart, bstop, dtset%dvdb_add_lr, ii &
433 18 : ])
434 2 : NCF_CHECK(ncerr)
435 :
436 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
437 : "nelect", "fermie"], &
438 : [wfk_ebands%nelect, wfk_ebands%fermie &
439 10 : ])
440 2 : NCF_CHECK(ncerr)
441 :
442 : ! arrays
443 2 : NCF_CHECK(nf90_put_var(ncid, vid("kpoints"), kpath%points(:,1:nk_path)))
444 2 : NCF_CHECK(nf90_put_var(ncid, vid("qpoints"), qpath%points(:,1:nq_path)))
445 6 : ABI_MALLOC(real_vec, (nq_path))
446 45 : real_vec = one
447 2 : NCF_CHECK(nf90_put_var(ncid, vid("qweights"), real_vec))
448 2 : ABI_FREE(real_vec)
449 2 : NCF_CHECK(nf90_put_var(ncid, vid("eph_fix_korq"), dtset%eph_fix_korq))
450 2 : NCF_CHECK(nf90_put_var(ncid, vid("eph_fix_wavevec"), dtset%eph_fix_wavevec))
451 :
452 : ! Compute non-analytic phonons for q--> 0 in polar materials.
453 2 : if (nq_path > 1 .and. (any(ifc%zeff /= zero))) call ifc%calcnwrite_nana_terms_qpath(qpath, cryst, ncid, units)
454 2 : NCF_CHECK(nf90_close(ncid))
455 : end if ! master
456 :
457 : ! Make sure the netcdf file has been written by master before continuing.
458 2 : call xmpi_barrier(comm)
459 :
460 : ! All procs open the GPATH file here.
461 : ! FIXME
462 : !NCF_CHECK(nctk_open_modify(ncid, gpath_path, comm))
463 2 : NCF_CHECK(nctk_open_modify(ncid, gpath_path, xmpi_comm_self))
464 :
465 : ! The cache allows one the reuse the wavefunctions of the previous k/q to init the NSCF cycle
466 : ! It usually reduces the number of iterations by 3-4 but it requires more memory.
467 2 : tot_nscf_ierr = 0; use_cache = .True. !; use_cache = .False.
468 3 : call ucache_k%init(use_cache .and. my_nkpath > 1, ngfft)
469 3 : call ucache_kq%init(use_cache .and. my_nqpath > 1, ngfft)
470 :
471 : ! Loop over spins (MPI parallelized)
472 4 : do my_is=1,my_nspins
473 2 : spin = my_spins(my_is)
474 :
475 2 : if (has_gwan) then
476 : ! Build wan for this spin from the ABIWAN.nc file, set the perturbation distribution
477 : ! from THIS routine's own pert_comm (same natom3 block-distribution used for the
478 : ! ab-initio gkq_atm below), then load g(R_e, R_p) from GWAN.nc.
479 0 : call wan%from_abiwan(dtfil%filabiwanin, spin, dtset%nsppol, .False., "", comm_my_is(my_is)%value)
480 0 : wan%my_pert_start = my_pert_start_glob; wan%my_npert = my_npert; wan%pert_comm => pert_comm
481 0 : call wan%load_gwan(dtfil%filgwanin, cryst, spin, dtset%nsppol, comm_my_is(my_is))
482 : end if
483 :
484 : ! Loop over k-points in k-path (MPI parallelized).
485 45 : do my_ik=1,my_nkpath
486 172 : ik = my_ik_inds(my_ik); kk = kpath%points(:, ik)
487 : !print *, "ik, kk", ik, kk
488 :
489 : ! Prepare NSCF run at k.
490 : ! gs_ham_k has pointers to the *_k arrays in output so we cannot deallocate them till the end.
491 : ! This is the reason why we use vlocal_k although this term does not depend on k
492 : call nscf%setup_kpt(spin, kk, istwfk_1, nband, cryst, dtset, psps, pawtab, pawfgr, & ! in
493 43 : npw_k, kg_k, kpg_k, ph3d_k, kinpw_k, ffnl_k, vlocal_k, cg_k, gsc_k, gs_ham_k) ! out
494 :
495 : ! Cache to initialize u_{nk}(g).
496 43 : use_cg_k = (my_ik > 1 .and. ucache_k%use_cache)
497 43 : if (use_cg_k) call ucache_k%get_kpt(kk, istwfk_1, npw_k, nspinor, nband, kg_k, cg_k)
498 :
499 : ! Compute u_{nk}(g)
500 : call nscf%solve_kpt(spin, kk, istwfk_1, nband, cryst, dtset, dtfil, gs_ham_k, &
501 43 : use_cg_k, npw_k, cg_k, gsc_k, eig_k, msg, ierr)
502 :
503 43 : ABI_WARNING_IF(ierr /= 0, msg)
504 43 : tot_nscf_ierr = tot_nscf_ierr + ierr
505 :
506 43 : call ucache_k%store_kpt(kk, istwfk_1, npw_k, nspinor, nband, kg_k, cg_k)
507 :
508 : !if (pert_comm%me == master) then
509 172 : NCF_CHECK(nf90_put_var(ncid, vid("all_eigens_k"), eig_k, start=[1,ik,spin]))
510 : !end if
511 :
512 : ! Make sure all procs in pert_comm have the same wavefunctions at k.
513 43 : call xmpi_bcast(cg_k, master, pert_comm%value, ierr)
514 43 : if (psps%usepaw == 1) call xmpi_bcast(gsc_k, master, pert_comm%value, ierr)
515 :
516 : ! Allocate vlocal. Note nvloc
517 258 : ABI_MALLOC(vlocal, (n4, n5, n6, gs_ham_k%nvloc))
518 :
519 : ! Loop over q-points in q-path (MPI parallelized).
520 : ! All procs in pert_comm enter this loop with the same ik/iq indices.
521 127 : do my_iq=1,my_nqpath
522 588 : iq = my_iq_inds(my_iq); qq = qpath%points(:,iq); qq_is_gamma = sum(qq**2) < tol14
523 336 : kq = kk + qq
524 : !print *, "iq, kq", iq, kq
525 :
526 : ! Prepare NSCF run at k+q.
527 : ! gs_ham_kq has pointers to the *_kq arrays in output so we cannot deallocate them till the end.
528 : ! This is the reason why we use vlocal_kq although this term does not depend on k+q.
529 : call nscf%setup_kpt(spin, kq, istwfk_1, nband, cryst, dtset, psps, pawtab, pawfgr, & ! in
530 84 : npw_kq, kg_kq, kpg_kq, ph3d_kq, kinpw_kq, ffnl_kq, vlocal_kq, cg_kq, gsc_kq, gs_ham_kq) ! out
531 :
532 : ! Cache for u_{m k+q}(g).
533 84 : use_cg_kq = (my_iq > 1 .and. ucache_kq%use_cache)
534 84 : if (use_cg_kq) call ucache_kq%get_kpt(kq, istwfk_1, npw_kq, nspinor, nband, kg_kq, cg_kq)
535 :
536 : ! We can use cg_k as input for the NSCF for a very quick return.
537 84 : if (qq_is_gamma) then
538 134284 : use_cg_kq = .True.; cg_kq = cg_k
539 : end if
540 :
541 : ! Compute u_{m k+q}(g)
542 : call nscf%solve_kpt(spin, kq, istwfk_1, nband, cryst, dtset, dtfil, gs_ham_kq, &
543 84 : use_cg_kq, npw_kq, cg_kq, gsc_kq, eig_kq, msg, ierr)
544 :
545 84 : ABI_WARNING_IF(ierr /= 0, msg)
546 84 : tot_nscf_ierr = tot_nscf_ierr + ierr
547 :
548 84 : call ucache_kq%store_kpt(kq, istwfk_1, npw_kq, nspinor, nband, kg_kq, cg_kq)
549 :
550 : ! This to have the same gauge when qq = 0.
551 134366 : if (qq_is_gamma) cg_kq = cg_k
552 :
553 : ! Make sure all procs in pert_comm have the same wavefunctions at k+q.
554 84 : call xmpi_bcast(cg_kq, master, pert_comm%value, ierr)
555 84 : if (psps%usepaw == 1) call xmpi_bcast(gsc_kq, master, pert_comm%value, ierr)
556 :
557 : ! Get phonons for this q-point.
558 84 : call ifc%fourq(cryst, qq, phfreqs, displ_cart, out_displ_red=displ_red_qq)
559 588 : phfreqs_eV = phfreqs * Ha_eV
560 :
561 : ! Only the global first k-point contributes to q-path-only variables.
562 84 : if (ik == 1) then
563 172 : NCF_CHECK(nf90_put_var(ncid, vid("all_eigens_kq"), eig_kq, start=[1,iq,spin]))
564 : ! Write phonons for this qpt.
565 43 : if (spin == 1) then
566 129 : NCF_CHECK(nf90_put_var(ncid, vid("phfreqs"), phfreqs_ev, start=[1,iq]))
567 258 : NCF_CHECK(nf90_put_var(ncid, vid("phdispl_cart"), displ_cart, start=[1,1,1,1,iq]))
568 : end if
569 : end if
570 :
571 : ! For PAW, one has to solve a generalized eigenproblem.
572 84 : gen_eigenpb = psps%usepaw == 1; sij_opt = 0; if (gen_eigenpb) sij_opt = 1
573 252 : ABI_MALLOC(gs1c, (2, npw_kq*nspinor*((sij_opt+1)/2)))
574 336 : ABI_MALLOC(h1_kets_kq, (2, npw_kq*nspinor, nb_in_g))
575 :
576 : ! ====================================
577 : ! Get DFPT potentials for this q-point
578 : ! ====================================
579 : ! After this branch we have allocated v1scf(cplex, nfftf, nspden, my_npert)).
580 84 : if (need_ftinterp) then
581 84 : call dvdb%get_ftqbz(qq, cplex, nfftf, ngfftf, v1scf, pert_comm%value)
582 : else
583 : ! Read and reconstruct the dvscf potentials for qq and my_npert perturbations.
584 0 : call dvdb%readsym_qbz(cryst, qq, qmap_symrec(:,iq), cplex, nfftf, ngfftf, v1scf, pert_comm%value)
585 : end if
586 :
587 : ! Allocate vlocal1 with correct cplex. Note nvloc.
588 504 : ABI_MALLOC(vlocal1, (cplex*n4, n5, n6, gs_ham_kq%nvloc))
589 :
590 : ! Load the k/k+q dependent parts of the Hamiltonian.
591 : ! NB: In this routine we have to use gs_ham_k to have {k+q}_H0_k.
592 : ! Using gs_ham_kq would be wrong as it would lead to {k+q}_H0_{k+q}.
593 :
594 : call gs_ham_k%load_kprime(kpt_kp=kq, npw_kp=npw_kq, istwf_kp=istwfk_1, kg_kp=kg_kq, kpg_kp=kpg_kq, kinpw_kp=kinpw_kq, &
595 84 : ph3d_kp=ph3d_kq, ffnl_kp=ffnl_kq, compute_ph3d=.true., compute_gbound=.true.)
596 : !call gs_ham_k%print([std_out], "gs_ham_k after load", dtset%prtvol)
597 :
598 : ! Loop over my atomic perturbations: apply H1_{kappa, alpha} and compute gkq_atm.
599 15708 : gkq_atm = zero
600 588 : do my_ip=1,my_npert
601 504 : idir = dvdb%my_pinfo(1, my_ip); ipert = dvdb%my_pinfo(2, my_ip); ipc = dvdb%my_pinfo(3, my_ip)
602 :
603 : ! Set up local potential vlocal1 with proper dimensioning from vtrial1 taking into account the spin.
604 : call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, gs_ham_kq%nvloc,&
605 504 : pawfgr, nscf%mpi_enreg, nscf%vtrial, v1scf(:,:,:,my_ip), vlocal, vlocal1)
606 :
607 : ! Prepare application of the NL part.
608 504 : call rf_ham_kq%init(cplex, gs_ham_k, ipert, has_e1kbsc=.true.)
609 504 : call rf_ham_kq%load_spin(spin, vlocal1=vlocal1, with_nonlocal=.true.)
610 : ! Load k-dependent part in the 1st-order Hamiltonian datastructure
611 : !call rf_ham_kq%load_k(npw_k=npw_k)
612 :
613 : ! Calculate dvscf * psi_k, results stored in h1_kets_kq on the k+q sphere.
614 : ! Compute H(1) applied to GS wavefunction Psi(0).
615 2016 : do in_k=1,nb_in_g
616 1512 : band_n = in_k + bstart - 1
617 1512 : eig0nk = eig_k(band_n)
618 : ! Use scissor shift on 0-order eigenvalue.
619 1512 : eshift = eig0nk - dtset%dfpt_sciss
620 :
621 : call getgh1c(berryopt0, cg_k(:,:,band_n), cwaveprj0, h1_kets_kq(:,:,in_k), &
622 : grad_berry, gs1c, gs_ham_k, gvnlx1, idir, ipert, [eshift], nscf%mpi_enreg, ndat1, optlocal, &
623 3528 : optnl, opt_gvnlx1, rf_ham_kq, sij_opt, tim_getgh1c, usevnl)
624 : !print *, "maxval(abs(h1_kets_kq(:,:,in_k))): ", maxval(abs(h1_kets_kq(:,:,in_k)))
625 : end do ! in_k
626 :
627 : ! Calculate <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation. No need to handle istwf_kq because it's always 1.
628 : !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(band_m)
629 2100 : do in_k=1,nb_in_g
630 6552 : do im_kq=1,nb_in_g
631 4536 : band_m = im_kq + bstart - 1
632 6048 : gkq_atm(:, im_kq, in_k, ipc) = cg_zdotc(npw_kq*nspinor, cg_kq(:,:,band_m), h1_kets_kq(:,:,in_k))
633 : end do
634 : end do
635 :
636 : end do ! my_ip
637 :
638 : ! Collect gkq_atm inside pert_comm so that all procs can operate on the data.
639 84 : if (pert_comm%nproc > 1) call xmpi_sum(gkq_atm, pert_comm%value, ierr)
640 :
641 : ! From atom to phonon mode representation. Results stored in gkq_nu.
642 84 : call ephtk_gkknu_from_atm(nb_in_g, nb_in_g, 1, natom, gkq_atm, phfreqs, displ_red_qq, gkq_nu)
643 : !print *, "gkq_atm:", gkq_atm; print *, "displ_red_qq:", displ_red_qq; print *, "gkq_nu:", gkq_nu
644 :
645 : ! Write |g|^2 for this q.
646 : !if (pert_comm%me == master) then
647 6720 : gkq2_nu = gkq_nu(1,:,:,:)**2 + gkq_nu(2,:,:,:)**2
648 588 : NCF_CHECK(nf90_put_var(ncid, vid("gkq2_nu"), gkq2_nu, start=[1,1,1,iq,ik,spin]))
649 : !end if
650 :
651 84 : if (has_gwan) then
652 : ! Interpolate e-ph matrix elements from GWAN.nc at this same (k,q) for comparison.
653 0 : ABI_MALLOC(g_atm_wan_local, (nwan_glob, nwan_glob, my_npert, 1))
654 : ! The routine already diagonalizes H^W(k) and H^W(k+q) to rotate g to the
655 : ! interpolated eigenstate basis, hence it also returns these eigenvalues.
656 : call wan%interp_eph_manyq(cryst, 1, qq, kk, g_atm_wan_local, &
657 0 : out_eigens_k=eig_k_wan, out_eigens_kq=eig_kq_wan)
658 0 : if (iq == 1) then
659 0 : NCF_CHECK(nf90_put_var(ncid, vid("all_eigens_wan_k"), eig_k_wan, start=[1,ik,spin]))
660 : end if
661 0 : if (ik == 1) then
662 0 : NCF_CHECK(nf90_put_var(ncid, vid("all_eigens_wan_kq"), eig_kq_wan(:,1), start=[1,iq,spin]))
663 : end if
664 :
665 0 : gkq_atm_wan = zero
666 0 : do my_ip=1,my_npert
667 0 : iglob = wan%my_pert_start + my_ip - 1
668 0 : gkq_atm_wan(1,:,:,iglob) = real(g_atm_wan_local(:,:,my_ip,1), kind=dp)
669 0 : gkq_atm_wan(2,:,:,iglob) = aimag(g_atm_wan_local(:,:,my_ip,1))
670 : end do
671 0 : if (pert_comm%nproc > 1) call xmpi_sum(gkq_atm_wan, pert_comm%value, ierr)
672 :
673 : ! Reuses phfreqs/displ_red_qq already computed above (ifc%fourq) for this same q --
674 : ! no new gauge issue, same literal-q call already used for the ab-initio side.
675 0 : call ephtk_gkknu_from_atm(nwan_glob, nwan_glob, 1, natom, gkq_atm_wan, phfreqs, displ_red_qq, gkq_nu_wan)
676 0 : gkq2_nu_wan = gkq_nu_wan(1,:,:,:)**2 + gkq_nu_wan(2,:,:,:)**2
677 0 : NCF_CHECK(nf90_put_var(ncid, vid("gkq2_nu_wan"), gkq2_nu_wan, start=[1,1,1,iq,ik,spin]))
678 :
679 0 : ABI_FREE(g_atm_wan_local)
680 : end if
681 :
682 84 : ABI_FREE(gs1c)
683 84 : ABI_FREE(vlocal1)
684 84 : ABI_FREE(v1scf)
685 84 : ABI_FREE(vlocal_kq)
686 84 : ABI_FREE(ph3d_kq)
687 84 : ABI_FREE(kpg_kq)
688 84 : ABI_FREE(kinpw_kq)
689 84 : ABI_FREE(ffnl_kq)
690 84 : ABI_FREE(kg_kq)
691 84 : ABI_FREE(eig_kq)
692 84 : ABI_FREE(cg_kq)
693 84 : ABI_FREE(gsc_kq)
694 84 : ABI_FREE(h1_kets_kq)
695 127 : call gs_ham_kq%free(); call rf_ham_kq%free()
696 : end do ! my_iq
697 :
698 43 : ABI_FREE(vlocal)
699 43 : ABI_FREE(vlocal_k)
700 43 : ABI_FREE(ph3d_k)
701 43 : ABI_FREE(kpg_k)
702 43 : ABI_FREE(kinpw_k)
703 43 : ABI_FREE(ffnl_k)
704 43 : ABI_FREE(kg_k)
705 43 : ABI_FREE(eig_k)
706 43 : ABI_FREE(cg_k)
707 43 : ABI_FREE(gsc_k)
708 45 : call gs_ham_k%free()
709 : end do ! my_ik
710 :
711 4 : if (has_gwan) call wan%free()
712 : end do ! my_is
713 :
714 2 : NCF_CHECK(nf90_close(ncid))
715 2 : call xmpi_barrier(comm)
716 :
717 : ! ===========================================
718 : ! Write results to ab_out for automatic tests
719 : ! ===========================================
720 2 : call xmpi_sum(tot_nscf_ierr, comm, ierr)
721 2 : tot_nscf_ierr = int(tot_nscf_ierr / dble(pert_comm%nproc))
722 :
723 2 : if (my_rank == master) then
724 2 : if (tot_nscf_ierr == 0) then
725 : call wrtout(units, &
726 2 : sjoin("Computation of g(k,q) completed. All NSCF runs converged within tolwfr: ", ftoa(dtset%tolwfr)), pre_newlines=1)
727 : else
728 0 : msg = sjoin("WARNING:", itoa(tot_nscf_ierr), "NSCF runs did not converge within tolwfr: ", ftoa(dtset%tolwfr))
729 0 : msg = sjoin(msg, ". Use nbdbuf and/or increase nstep!")
730 0 : call wrtout(ab_out, msg)
731 0 : ABI_WARNING(msg)
732 : end if
733 :
734 2 : NCF_CHECK(nctk_open_read(ncid, gpath_path, xmpi_comm_self))
735 :
736 : ! Write k/q wavevectors.
737 2 : call wrtout(units, "kpoints:")
738 45 : do ik=1, nk_path
739 45 : call wrtout(units, sjoin(char(9), itoa(ik), ktoa(kpath%points(:,ik))))
740 : end do
741 2 : call wrtout(units, "qpoints:")
742 45 : do iq=1, nq_path
743 45 : call wrtout(units, sjoin(char(9), itoa(iq), ktoa(qpath%points(:,iq))))
744 : end do
745 :
746 : ! Write KS eigenvalues.
747 2 : if (nq_path > 1) then
748 3 : ABI_MALLOC(eig_kq, (nband))
749 2 : do spin=1,nsppol
750 1 : call wrtout(units, sjoin(" Energies_kq in eV for spin:", itoa(spin)))
751 44 : do iq=1, nq_path
752 42 : if (all(iq /= [1, 2, nq_path-1, nq_path])) cycle
753 16 : NCF_CHECK(nf90_get_var(ncid, vid("all_eigens_kq"), eig_kq, start=[1,iq,spin]))
754 17 : do ii=0,(nband-1)/8
755 84 : write(msg, '(a, 8es16.6)' )' ene:',(eig_kq(band_m) * Ha_eV, band_m=1+ii*8,min(nband,8+ii*8))
756 54 : call wrtout(units, msg)
757 : end do
758 : end do
759 : end do
760 1 : ABI_FREE(eig_kq)
761 :
762 : else
763 : ! nk_path > 1
764 3 : ABI_MALLOC(eig_k, (nband))
765 2 : do spin=1,nsppol
766 1 : call wrtout(units, sjoin(" Energies_k in eV for spin:", itoa(spin)))
767 44 : do ik=1, nk_path
768 168 : NCF_CHECK(nf90_get_var(ncid, vid("all_eigens_k"), eig_k, start=[1,ik,spin]))
769 42 : if (all(ik /= [1, 2, nk_path-1, nk_path])) cycle
770 17 : do ii=0,(nband-1)/8
771 84 : write(msg, '(a, 8es16.6)' )' ene:',(eig_k(band_n) * Ha_eV, band_n=1+ii*8,min(nband,8+ii*8))
772 54 : call wrtout(units, msg)
773 : end do
774 : end do
775 : end do
776 1 : ABI_FREE(eig_k)
777 : end if
778 :
779 : ! Save g^2 in nc format without any average. This operation will be performed by AbiPy (need ph freqs and eigenergies)
780 2 : call wrtout(units, " Writing sqrt(1/N_b^2 \sum_{mn} |g_{mn,nu}(k, q)|^2) in meV for testing purpose.", pre_newlines=2)
781 2 : write(msg, "(1x,4(a5,1x),a16)") "nu","iq", "ik", "spin", "|g| in meV"
782 2 : call wrtout(units, msg)
783 :
784 4 : do spin=1,nsppol
785 47 : do ik=1, nk_path
786 43 : if (all(ik /= [1, 2, nk_path-1, nk_path])) cycle
787 53 : do iq=1, nq_path
788 46 : if (all(iq /= [1, 2, nq_path-1, nq_path])) cycle
789 56 : NCF_CHECK(nf90_get_var(ncid, vid("gkq2_nu"), gkq2_nu, start=[1,1,1,iq,ik,spin]))
790 99 : do nu=1,natom3
791 624 : write(msg, "(1x,4(i5,1x),es16.6)") nu, iq, ik, spin, sqrt(sum(gkq2_nu(:,:, nu)) / nb_in_g**2) * Ha_meV
792 94 : call wrtout(units, msg)
793 : end do
794 : end do ! iq
795 : end do ! ik
796 : end do ! spin
797 :
798 2 : if (has_gwan) then
799 : call wrtout(units, &
800 0 : " Writing sqrt(1/N^2 sum |g|^2): ab-initio vs Wannier-interpolated (gauge-invariant trace), meV.", pre_newlines=2)
801 0 : write(msg, "(1x,4(a5,1x),2(a18,1x))") "nu","iq","ik","spin", "|g|_abinitio(meV)", "|g|_wannier(meV)"
802 0 : call wrtout(units, msg)
803 0 : do spin=1,nsppol
804 0 : do ik=1, nk_path
805 0 : if (all(ik /= [1, 2, nk_path-1, nk_path])) cycle
806 0 : do iq=1, nq_path
807 0 : if (all(iq /= [1, 2, nq_path-1, nq_path])) cycle
808 0 : NCF_CHECK(nf90_get_var(ncid, vid("gkq2_nu"), gkq2_nu, start=[1,1,1,iq,ik,spin]))
809 0 : NCF_CHECK(nf90_get_var(ncid, vid("gkq2_nu_wan"), gkq2_nu_wan, start=[1,1,1,iq,ik,spin]))
810 0 : do nu=1,natom3
811 0 : write(msg, "(1x,4(i5,1x),2(es18.6,1x))") nu, iq, ik, spin, &
812 0 : sqrt(sum(gkq2_nu(:,:,nu)) / nb_in_g**2) * Ha_meV, sqrt(sum(gkq2_nu_wan(:,:,nu)) / nwan_glob**2) * Ha_meV
813 0 : call wrtout(units, msg)
814 : end do
815 : end do ! iq
816 : end do ! ik
817 : end do ! spin
818 : end if
819 :
820 2 : NCF_CHECK(nf90_close(ncid))
821 : end if ! master
822 :
823 : ! Free memory.
824 2 : ABI_FREE(my_ik_inds)
825 2 : ABI_FREE(my_iq_inds)
826 2 : ABI_FREE(gvnlx1)
827 2 : ABI_FREE(grad_berry)
828 2 : ABI_FREE(qmap_symrec)
829 2 : ABI_FREE(gkq_atm)
830 2 : ABI_FREE(gkq_nu)
831 2 : ABI_FREE(gkq2_nu)
832 2 : if (has_gwan) then
833 0 : ABI_FREE(gkq_atm_wan)
834 0 : ABI_FREE(gkq_nu_wan)
835 0 : ABI_FREE(gkq2_nu_wan)
836 0 : ABI_FREE(eig_k_wan)
837 0 : ABI_FREE(eig_kq_wan)
838 : end if
839 2 : ABI_FREE(displ_cart)
840 2 : ABI_FREE(displ_red_qq)
841 2 : ABI_FREE(my_spins)
842 :
843 2 : call pawcprj_free(cwaveprj0)
844 2 : ABI_FREE(cwaveprj0)
845 4 : do my_is=1,my_nspins
846 4 : call comm_my_is(my_is)%free()
847 : end do
848 2 : ABI_FREE(comm_my_is)
849 :
850 2 : call qpath%free(); call kpath%free(); call ucache_k%free(); call ucache_kq%free()
851 2 : call qpt_comm%free(); call kpt_comm%free(); call pert_comm%free(); call nscf%free()
852 :
853 8 : call cwtime_report(" eph_path: MPI barrier before returning.", cpu_all, wall_all, gflops_all, end_str=ch10, comm=comm)
854 : !stop
855 :
856 : contains
857 324 : integer function vid(var_name)
858 : character(len=*),intent(in) :: var_name
859 324 : vid = nctk_idname(ncid, var_name)
860 : end function vid
861 :
862 : end subroutine eph_path_run
863 : !!***
864 :
865 : end module m_eph_path
866 : !!***
|
