Branch data 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 : : !!***
|
