Line data Source code
1 : !!****m* ABINIT/m_phgamma
2 : !! NAME
3 : !!
4 : !! FUNCTION
5 : !! Computation of phonon linewidths, isotropic superconducting properties,
6 : !! and transport properties in metals within the LOVA approximation
7 : !! to the linearized Boltzmann equation.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! TODO
16 : !! 1) Implement restart capabilities (eph_restart).
17 : !! 2) Sum over the IBZ(q) on the FS instead of full BZ.
18 : !! 3) Gaussian adaptive smearing for double delta (tetra version from libtetrabz really slow!).
19 : !! 4) Automatic detections of energy window, improve filtering techniques.
20 : !! 5) Interface with KERANGE trick
21 : !! 6) More examples and tutorials (using precomputed Netcdf files).
22 : !! 7) Implement SKW interpolation for ph linewidths and/or linear interpolation (I don't trust plain Fourier interpolation).
23 : !! 8) Perform more benchmarks with dense meshes to detect hotspots and memory bottlenecks.
24 : !! 9) Test spin and SOC.
25 : !!
26 : !! SOURCE
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : module m_phgamma
35 :
36 : use, intrinsic :: iso_c_binding
37 : use defs_basis
38 : use m_abicore
39 : use m_xmpi
40 : use m_errors
41 : use m_krank
42 : use m_htetra
43 : use m_tetrahedron
44 : use libtetrabz
45 : use m_ifc
46 : use m_ebands
47 : use m_fstab
48 : use m_nctk
49 : use netcdf
50 : use m_wfk
51 : use m_ddb
52 : use m_ddk
53 : use m_dvdb
54 : use m_crystal
55 : use m_fft
56 : use m_hamiltonian
57 : use m_pawcprj
58 : use m_dtset
59 : use m_dtfil
60 : use m_ephtk
61 : use m_mkffnl
62 :
63 : use defs_abitypes, only : mpi_type
64 : use defs_datatypes, only : pseudopotential_type
65 : use m_time, only : cwtime, cwtime_report
66 : use m_fstrings, only : toupper, itoa, sjoin, ktoa, ltoa, strcat
67 : use m_numeric_tools, only : arth, wrap2_pmhalf, simpson_int, simpson, mkherm, get_diag, isdiagmat
68 : use m_io_tools, only : open_file, iomode_from_fname
69 : use m_symtk, only : littlegroup_q
70 : use m_geometry, only : normv
71 : use m_special_funcs, only : gaussian
72 : use m_fftcore, only : ngfft_seq, get_kg
73 : use m_cgtools, only : cg_zdotc
74 : use m_kg, only : getph, mkkpg
75 : use m_dynmat, only : symdyma, ftgam_init, ftgam, asrif9
76 : use m_bz_mesh, only : kpath_t
77 : use m_special_funcs, only : fermi_dirac
78 : use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt, kpts_timrev_from_kptopt, kpts_map
79 : use defs_elphon, only : complete_gamma !, complete_gamma_tr
80 : use m_getgh1c, only : getgh1c, rf_transgrid_and_pack
81 : use m_pawang, only : pawang_type
82 : use m_pawrad, only : pawrad_type
83 : use m_pawtab, only : pawtab_type
84 : use m_pawfgr, only : pawfgr_type
85 : use m_wfd, only : wfd_t
86 : use m_pstat, only : pstat_proc
87 : use m_lgroup, only : lgroup_t
88 : use m_abi_linalg, only : abi_gpu_xgemm_d
89 :
90 : implicit none
91 :
92 : private
93 : !!***
94 :
95 : public :: eph_phgamma
96 : public :: find_ewin
97 :
98 : !----------------------------------------------------------------------
99 :
100 : !!****t* m_phgamma/phgamma_t
101 : !! NAME
102 : !! phgamma_t
103 : !!
104 : !! FUNCTION
105 : !! Provides methods for computing phonon linewidths, interpolating the results
106 : !! in q-space and evaluate superconducting properties within the isotropic formalism.
107 : !!
108 : !! SOURCE
109 :
110 : type,public :: phgamma_t
111 :
112 : integer :: natom = -1
113 : ! Number of atoms per unit cell.
114 :
115 : integer :: natom3 = -1
116 : ! Number of phonon branches i.e. 3*natom.
117 :
118 : integer :: nsppol = -1
119 : ! Number of independent spin polarizations.
120 :
121 : integer :: nspinor = -1
122 : ! Number of spinorial components.
123 :
124 : integer :: nqibz = -1
125 : ! Number of q-points in the IBZ.
126 :
127 : integer :: my_nqibz = -1
128 : ! Number of q-points in the IBZ treated by the current MPI processor
129 :
130 : integer :: nqbz = -1
131 : ! Number of q-points in the BZ.
132 :
133 : integer :: eph_scalprod = 0
134 : ! This to call anaddb routines. Note that eph_scalprod 1 is not supported in eph.
135 :
136 : integer :: bcorr = 0
137 : ! 1 to include Blochl correction in the tetrahedron method else 0.
138 :
139 : integer :: prteliash = 0
140 : ! This flag activates the computation of the Eliashberg function.
141 :
142 : integer :: nrpt = -1
143 : ! Number of R-points in the real space representation of the gamma matrices.
144 :
145 : integer :: symgamma
146 : ! 1 if gamma matrices should be symmetrized by symdyma when using Fourier interpolation
147 :
148 : integer :: asr
149 : ! If the "Acoustic rule" at Gamma should be enforced.
150 :
151 : integer :: ndir_transp = 0
152 : ! 0 if no transport, otherwise 3.
153 :
154 : integer :: ngqpt(3)
155 : ! Number of divisions in the q-mesh.
156 :
157 : integer :: nene = -1
158 : ! Number of chemical potential values used for inelastic integration
159 :
160 : !integer :: my_nqpt
161 : !integer,allocatable :: my_iqpt(:)
162 :
163 : integer :: my_nfsk_q = -1
164 : ! Number of k-points in the FS treated by this MPI processor for a given q.
165 : ! Computed in phgamma_setup_qpoint.
166 :
167 : integer,allocatable :: my_ifsk_q(:)
168 : ! Index of the FS k-points treated by this processor for a given q
169 : ! Computed in phgamma_setup_qpoint.
170 :
171 : integer :: my_nspins = -1
172 : ! Number of spins treated by the MPI rank.
173 :
174 : integer,allocatable :: my_spins(:)
175 : ! my_spins(my_nspins)
176 : ! Indirect table giving the spin indices treated by this rank.
177 : ! Used only the collinear case with nspinor == 1
178 :
179 : !integer :: my_npert
180 : ! Number of atomic perturbations or phonon modes treated by this MPI rank.
181 :
182 : !integer(i1b),allocatable :: itreat_qibz(:)
183 : ! itreat_qibz(nqibz)
184 : ! Table used to distribute potentials over q-points in the IBZ.
185 : ! The loop over qpts in the IBZ(k) is MPI distributed inside qpt_comm accordinging to this table.
186 : ! 0 if this IBZ point is not treated by this proc.
187 : ! 1 if this IBZ is treated.
188 :
189 : !integer,allocatable :: my_pinfo(:,:)
190 : ! my_pinfo(3, my_npert)
191 : ! my_pinfo(1, ip) gives the `idir` index of the ip-th perturbation.
192 : ! my_pinfo(2, ip) gives the `ipert` index of the ip-th perturbation.
193 : ! my_pinfo(3, ip) gives `pertcase`=idir + (ipert-1)*3
194 :
195 : !integer,allocatable :: pert_table(:,:)
196 : ! pert_table(2, natom3)
197 : ! pert_table(1, npert): rank of the processor treating this atomic perturbation.
198 : ! pert_table(2, npert): my_ip index in my_pinfo table, -1 if this rank is not treating ipert.
199 :
200 : integer, allocatable :: my_iqibz(:)
201 : ! indices of ibz iq in local array. -1 if iq does not belong to current proc
202 :
203 : integer, allocatable :: my_iqbz(:)
204 : ! List of iq_ibz indices treated by current proc
205 :
206 : real(dp) :: enemin
207 : ! Minimal chemical potential value used for inelastic integration Copied from fstab
208 :
209 : real(dp) :: deltaene
210 : ! Chemical potential increment for inelastic integration Copied from fstab
211 : ! for simplicity could be made equal to phonon frequency step
212 :
213 : real(dp) :: gprim(3,3)
214 : ! Needed for Fourier interpolation.
215 : ! NOTE: gprim (not gprimd) is used for all FT interpolations,
216 : ! to be consistent with the dimensions of the rpt, which come from anaddb.
217 :
218 : real(dp),allocatable :: n0(:)
219 : ! (%nsppol)
220 : ! Density of states at the Fermi level per spin in a.u.
221 :
222 : real(dp),allocatable :: qibz(:,:)
223 : ! qibz(3,nqibz)
224 : ! Reduced coordinates of the q-points in the IBZ.
225 :
226 : real(dp),allocatable :: wtq(:)
227 : ! wtq(nqibz)
228 : ! Weights of the q-points in the IBZ (normalized to one)
229 :
230 : real(dp),allocatable :: qbz(:,:)
231 : ! qbz(3,nqbz)
232 : ! Reduced coordinates of the q-points in the BZ.
233 :
234 : real(dp),allocatable :: rpt(:,:)
235 : ! rpt(3,nrpt)
236 : ! Reduced coordinates ***in terms of rprim*** of the lattice points used
237 : ! for the Fourier transform of the phonon linewidths.
238 :
239 : real(dp),allocatable :: wghatm(:,:,:)
240 : ! wghatm(natom,natom,nrpt)
241 : ! Weights used in the FT of the phonon linewidths.
242 :
243 : real(dp),allocatable :: vals_qibz(:,:,:,:,:)
244 : ! vals_qibz(2,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ.
245 : ! vals_qibz {\tau'\alpha',\tau\alpha} = sum over k, im_kq, in_k
246 : ! <psi_{k+q,im_kq} | H(1)_{\tau'\alpha'} | psi_{k,in_k}>* \cdot
247 : ! <psi_{k+q,im_kq} | H(1)_{\tau \alpha } | psi_{k,in_k}>
248 :
249 : !NOTE: choice to put nsppol before or after nqbz is a bit arbitrary
250 : ! abinit uses nband,nkpt,nsppol, but here for convenience nkpt_phon,nsppol,nqbz interpolation is on qpt
251 : !MG: I think that nsppol should be the last dimension set call to ftgam.
252 :
253 : real(dp),allocatable :: vals_rpt(:,:,:,:)
254 : ! vals_rpt(2,natom3**2,nrpt,nsppol)
255 : ! tgamma matrices in real space in reduced coordinates. Used for the Fourier interpolation.
256 :
257 : ! transport stuff with velocity factors
258 : real(dp),allocatable :: vals_in_qibz(:,:,:,:,:,:)
259 : real(dp),allocatable :: vals_out_qibz(:,:,:,:,:,:)
260 : ! vals_XX_qibz(2,9,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ.
261 :
262 : real(dp),allocatable :: vals_in_rpt(:,:,:,:,:)
263 : real(dp),allocatable :: vals_out_rpt(:,:,:,:,:)
264 : ! vals_XX_rpt(2,9,natom3**2,nrpt,nsppol)
265 : ! tgamma matrices in real space in reduced coordinates.
266 :
267 : ! gamma matrices keeping full electron energy dependency
268 : real(dp),allocatable :: vals_ee(:,:,:,:,:,:,:)
269 : ! vals_eew(2, nene, nene, natom3, natom3, nqibz, nsppol)
270 :
271 : contains
272 :
273 : procedure :: free => phgamma_free
274 : ! Free memory.
275 :
276 : procedure :: interp => phgamma_interp
277 : ! Interpolates the phonon linewidths.
278 :
279 : procedure :: eval_qibz => phgamma_eval_qibz
280 : ! Evaluate phonon linewidths without Fourier interpolation.
281 :
282 : procedure :: interp_setup => phgamma_interp_setup
283 : ! Compute the tables used for the interpolation in q-space.
284 :
285 : procedure :: linwid => phgamma_linwid
286 : ! Interpolate linewidths along an arbitrary q-path.
287 :
288 : end type phgamma_t
289 :
290 : public :: phgamma_init ! Creation method.
291 : !!***
292 :
293 : !----------------------------------------------------------------------
294 :
295 : !!****t* m_phgamma/a2fw_t
296 : !! NAME
297 : !! a2fw_t
298 : !!
299 : !! FUNCTION
300 : !! Store the Eliashberg function a2F(w).
301 : !!
302 : !! SOURCE
303 :
304 : type,public :: a2fw_t
305 :
306 : integer :: nomega
307 : ! Number of frequency points in a2f(w).
308 :
309 : integer :: nsppol
310 : ! Number of independent spin polarizations.
311 :
312 : integer :: natom3
313 : ! Number of phonon modes.
314 :
315 : integer :: nene
316 : ! Number of chemical potential values used for inelastic integration
317 :
318 : real(dp) :: enemin
319 : ! Minimal chemical potential value used for inelastic integration Copied from fstab
320 :
321 : real(dp) :: deltaene
322 : ! Chemical potential increment for inelastic integration Copied from fstab
323 : ! for simplicity could be made equal to phonon frequency step
324 :
325 : real(dp) :: omega_min, omega_max
326 : ! min and Max frequency (Ha) in the linear mesh.
327 :
328 : real(dp) :: wstep
329 : ! Step of the linear mesh
330 :
331 : real(dp) :: smear
332 : ! Gaussian broadening used to approximated the Dirac distribution.
333 :
334 : integer :: nqshift
335 : ! Number of shifts in the q-mesh
336 :
337 : integer :: ngqpt(3)
338 : ! The q-mesh used for calculating vals(w).
339 :
340 : real(dp),allocatable :: qshift(:,:)
341 : ! qshift(3,nqshift)
342 : ! The shifts used to generate the q-mesh.
343 :
344 : real(dp),allocatable :: n0(:)
345 : ! n0(nsppol)
346 : ! Electronic DOS at the Fermi level.
347 :
348 : real(dp),allocatable :: omega(:)
349 : ! omega(nomega)
350 : ! Frequency mesh in Hartree (linear).
351 :
352 : real(dp),allocatable :: vals(:,:,:)
353 : ! vals(nomega, 0:natom3, nsppol)
354 : ! Eliashberg function
355 : ! vals(w, 1:natom3, 1:nsppol): a2f(w) decomposed per phonon branch and spin
356 : ! vals(w, 0 , 1:nsppol): a2f(w) summed over phonons modes, decomposed in spin
357 :
358 : real(dp),allocatable :: vals_ee(:,:,:,:)
359 : ! vals_ee(nene,nene,nomega,nsppol)
360 : ! Eliashberg function
361 : ! vals(e,e',w,0,1:nsppol): a2f(e,e',w) summed over phonons modes, decomposed in spin
362 :
363 : real(dp),allocatable :: lambdaw(:,:,:)
364 : ! lambda(nomega,0:natom3,nsppol)
365 :
366 : contains
367 :
368 : procedure :: init => a2fw_init !
369 : ! Calculates the FS averaged alpha^2F(w) function.
370 :
371 : procedure :: free => a2fw_free
372 : ! Free the memory allocated in the structure.
373 :
374 : procedure :: write => a2fw_write
375 : ! Write alpha^2F(w) to an external file in text/netcdf format
376 :
377 : procedure :: get_moment => a2fw_get_moment
378 : ! Compute moments of alpha^2F(w)/w .
379 : end type a2fw_t
380 : !!***
381 :
382 : !!****t* m_phgamma/a2fw_tr_t
383 : !! NAME
384 : !! a2fw_tr_t
385 : !!
386 : !! FUNCTION
387 : !! Store the Eliashberg transport spectral functions:
388 : !! a2F_trin(w, x, x')
389 : !! a2F_trout(w, x, x')
390 : !! a2F_tr(w, x, x') = in - out
391 : !! a2F_tr_gen(e, e', w, x, x')
392 : !!
393 : !! SOURCE
394 :
395 : type,public :: a2fw_tr_t
396 :
397 : integer :: nomega
398 : ! Number of frequency points in a2f_tr(w).
399 :
400 : integer :: nene
401 : ! Number of chemical potential values used for inelastic integration
402 : ! Number of electron points in a2f_tr_gen(e,e',w).
403 :
404 : integer :: nsppol
405 : ! Number of independent spin polarizations.
406 :
407 : integer :: natom3
408 : ! Number of phonon modes.
409 :
410 : real(dp) :: enemin
411 : ! Minimal chemical potential value used for inelastic integration
412 :
413 : real(dp) :: deltaene
414 : ! Chemical potential increment for inelastic integration
415 :
416 : real(dp) :: omega_min,omega_max
417 : ! min and Max frequency (Ha) in the linear mesh.
418 :
419 : real(dp) :: wstep
420 : ! Step of the linear mesh
421 :
422 : real(dp) :: smear
423 : ! Gaussian broadening used to approximated the Dirac distribution.
424 :
425 : integer :: nqshift
426 : ! Number of shifts in the q-mesh
427 :
428 : integer :: ngqpt(3)
429 : ! The q-mesh used for calculating vals(w).
430 :
431 : real(dp),allocatable :: qshift(:,:)
432 : ! qshift(3,nqshift)
433 : ! The shifts used to generate the q-mesh.
434 :
435 : real(dp),allocatable :: n0(:)
436 : ! n0(nsppol)
437 : ! Electronic DOS at the Fermi level.
438 :
439 : real(dp),allocatable :: omega(:)
440 : ! omega(nomega)
441 : ! Frequency mesh in Hartree (linear).
442 :
443 : real(dp),allocatable :: vals_in(:,:,:,:,:)
444 : real(dp),allocatable :: vals_out(:,:,:,:,:)
445 : ! vals_in(nomega,3,3,0:natom3,nsppol)
446 : ! Eliashberg transport functions for in and out scattering
447 : ! vals_in(w,3,3,1:natom3,1:nsppol): a2f_tr(w) decomposed per phonon branch and spin
448 : ! vals_in(w,3,3,0,1:nsppol): a2f_tr(w) summed over phonons modes, decomposed in spin
449 :
450 : real(dp),allocatable :: vals_tr(:,:,:,:,:)
451 : ! vals_tr(nomega,3,3,0:natom3,nsppol)
452 : ! transport spectral function = in-out
453 :
454 : real(dp),allocatable :: vals_tr_gen(:,:,:,:,:,:,:)
455 : ! vals(nene,nene,nomega,3,3,nsppol)
456 : ! generalized transport spectral function from PB Allen Phys. Rev. Lett. 59, 1460 (1987) [[cite:Allen1987]]
457 :
458 : real(dp),allocatable :: lambdaw_tr(:,:,:,:,:)
459 : ! lambda(nomega,3,3,0:natom3,nsppol)
460 :
461 : contains
462 :
463 : procedure :: init => a2fw_tr_init
464 : ! Calculates the FS averaged alpha^2F_tr,in,out(w, x, x') functions.
465 :
466 : procedure :: free => a2fw_tr_free
467 : ! Free the memory allocated in the structure.
468 :
469 : procedure :: write => a2fw_tr_write
470 : ! Write alpha^2F(w) to an external file in text/netcdf format
471 :
472 : end type a2fw_tr_t
473 : !!***
474 :
475 : real(dp),private,parameter :: EPH_Q0TOL = 0.01_dp
476 :
477 : contains !=====================================================
478 : !!***
479 :
480 : !----------------------------------------------------------------------
481 :
482 : !!****f* m_phgamma/phgamma_free
483 : !! NAME
484 : !! phgamma_free
485 : !!
486 : !! FUNCTION
487 : !! Free the dynamic memory
488 : !!
489 : !! SOURCE
490 :
491 5 : subroutine phgamma_free(gams)
492 :
493 : !Arguments ------------------------------------
494 : class(phgamma_t),intent(inout) :: gams
495 : ! *************************************************************************
496 :
497 : !real
498 5 : ABI_SFREE(gams%n0)
499 5 : ABI_SFREE(gams%qibz)
500 5 : ABI_SFREE(gams%wtq)
501 5 : ABI_SFREE(gams%qbz)
502 5 : ABI_SFREE(gams%rpt)
503 5 : ABI_SFREE(gams%wghatm)
504 5 : ABI_SFREE(gams%vals_qibz)
505 5 : ABI_SFREE(gams%vals_rpt)
506 5 : ABI_SFREE(gams%vals_in_qibz)
507 5 : ABI_SFREE(gams%vals_in_rpt)
508 5 : ABI_SFREE(gams%vals_out_qibz)
509 5 : ABI_SFREE(gams%vals_out_rpt)
510 5 : ABI_SFREE(gams%vals_ee)
511 5 : ABI_SFREE(gams%my_iqibz)
512 5 : ABI_SFREE(gams%my_ifsk_q)
513 5 : ABI_SFREE(gams%my_spins)
514 :
515 5 : end subroutine phgamma_free
516 : !!***
517 :
518 : !----------------------------------------------------------------------
519 :
520 : !!****f* m_phgamma/phgamma_init
521 : !! NAME
522 : !! phgamma_init
523 : !!
524 : !! FUNCTION
525 : !! Creation method for the phgamma_t datatype.
526 : !!
527 : !! INPUTS
528 : !! cryst<crystal_t>
529 : !! ifc<ifc_type>=Interatomic force constants.
530 : !! symdynmat=1 to activa symmetrization of gamma matrices.
531 : !! ngqpt(3)=Q-mesh divisions
532 : !! n0(dtset%nsppol)=Density of states at the Fermi level per spin.
533 : !!
534 : !! OUTPUT
535 : !! gams<phgamma_t>
536 : !!
537 : !! SOURCE
538 :
539 5 : subroutine phgamma_init(gams, cryst, ifc, ebands, fstab, dtset, eph_scalprod, ngqpt, n0, comm)
540 :
541 : !Arguments ------------------------------------
542 : !scalars
543 : class(phgamma_t),intent(out) :: gams
544 : integer,intent(in) :: eph_scalprod, comm
545 : type(crystal_t),intent(in) :: cryst
546 : type(ifc_type),intent(in) :: ifc
547 : type(ebands_t),intent(in) :: ebands
548 : type(fstab_t), intent(in) :: fstab
549 : type(dataset_type),intent(in) :: dtset
550 : !arrays
551 : integer,intent(in) :: ngqpt(3)
552 : real(dp),intent(in) :: n0(dtset%nsppol)
553 :
554 : !Local variables-------------------------------
555 : !scalars
556 : integer :: my_rank, nproc, ierr, nsppol, natom3, qptopt, qtimrev
557 : !arrays
558 : integer :: qptrlatt(3,3)
559 : ! *************************************************************************
560 :
561 5 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
562 :
563 : ! Set basic dimensions.
564 5 : nsppol = dtset%nsppol; gams%nsppol = nsppol; gams%nspinor = dtset%nspinor
565 :
566 5 : gams%natom = cryst%natom; gams%natom3 = 3*cryst%natom; natom3 = gams%natom3
567 :
568 5 : gams%symgamma = dtset%symdynmat; gams%eph_scalprod = eph_scalprod
569 : !gams%asr = ifc%asr
570 5 : gams%asr = 0
571 5 : gams%prteliash = dtset%prteliash
572 :
573 5 : gams%ndir_transp = 0; if (dtset%eph_transport > 0) gams%ndir_transp = 3
574 :
575 : ! Copy DOS.
576 15 : ABI_MALLOC(gams%n0, (nsppol))
577 15 : gams%n0 = n0
578 :
579 5 : gams%nene = fstab%nene
580 5 : gams%enemin = fstab%enemin
581 5 : gams%deltaene = fstab%deltaene
582 :
583 : ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion
584 20 : gams%ngqpt = ngqpt
585 5 : qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
586 5 : qptopt = ebands%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt
587 : qtimrev = kpts_timrev_from_kptopt(qptopt)
588 :
589 : call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt, 1, [zero, zero, zero], &
590 5 : gams%nqibz, gams%qibz, gams%wtq, gams%nqbz, gams%qbz)
591 :
592 : ! Allocate matrices in the IBZ.
593 30 : ABI_MALLOC_OR_DIE(gams%vals_qibz, (2, natom3, natom3, gams%nqibz, nsppol), ierr)
594 1901 : gams%vals_qibz = zero
595 :
596 5 : if (dtset%eph_transport > 0) then
597 30 : ABI_MALLOC_OR_DIE(gams%vals_in_qibz, (2, 9, natom3, natom3, gams%nqibz, nsppol), ierr)
598 15626 : gams%vals_in_qibz = zero
599 30 : ABI_MALLOC_OR_DIE(gams%vals_out_qibz, (2, 9, natom3, natom3, gams%nqibz, nsppol), ierr)
600 15626 : gams%vals_out_qibz = zero
601 : end if
602 :
603 5 : if (gams%prteliash == 3) then
604 0 : ABI_MALLOC_OR_DIE(gams%vals_ee, (2, gams%nene, gams%nene, natom3, natom3, gams%nqibz, nsppol), ierr)
605 0 : gams%vals_ee = zero
606 : end if
607 :
608 : ! Prepare Fourier interpolation.
609 65 : gams%gprim = ifc%gprim
610 5 : gams%nrpt = ifc%nrpt
611 15 : ABI_MALLOC(gams%rpt, (3, gams%nrpt))
612 1870 : gams%rpt = ifc%rpt
613 25 : ABI_MALLOC(gams%wghatm, (gams%natom, gams%natom, gams%nrpt))
614 1405 : gams%wghatm = ifc%wghatm
615 :
616 5 : end subroutine phgamma_init
617 : !!***
618 :
619 : !----------------------------------------------------------------------
620 :
621 : !!****f* m_phgamma/phgamma_ncwrite
622 : !! NAME
623 : !! phgamma_ncwrite
624 : !!
625 : !! FUNCTION
626 : !! Write the results stored in the phgamma_t datatype to netcdf file.
627 : !!
628 : !! INPUTS
629 : !! cryst<crystal_t>=Crystalline structure.
630 : !! ifc<ifc_type>=Interatomic force constants.
631 : !! ncid=Netcdf file handler (already open in the caller).
632 : !!
633 : !! SOURCE
634 :
635 5 : subroutine phgamma_ncwrite(gams, cryst, ifc, ncid)
636 :
637 : !Arguments ------------------------------------
638 : !scalars
639 : class(phgamma_t),intent(inout) :: gams
640 : type(crystal_t),intent(in) :: cryst
641 : type(ifc_type),intent(in) :: ifc
642 : integer,intent(in) :: ncid
643 :
644 : !Local variables-------------------------------
645 : !scalars
646 : integer :: max_nqpt = 50
647 : integer :: iq_ibz,spin,mu, units(2)
648 : real(dp) :: lambda_tot, cpu, wall, gflops
649 : character(len=500) :: msg
650 : !arrays
651 10 : real(dp) :: phfrq(3*cryst%natom), gamma_ph(3*cryst%natom), lambda_ph(3*cryst%natom)
652 10 : real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
653 : ! *************************************************************************
654 :
655 15 : units = [std_out, ab_out]
656 5 : call cwtime(cpu, wall, gflops, "start")
657 :
658 : ! Write data to files for each q point, also compute total lambda.
659 5 : lambda_tot = zero
660 10 : do spin=1,gams%nsppol
661 66 : do iq_ibz=1,gams%nqibz
662 :
663 : ! Get phonon frequencies, gamma(q,nu) and lambda(q,nu).
664 : ! Quantities are already summed over (collinear) spin channels if nsppol == 1.
665 61 : call gams%eval_qibz(cryst, ifc, iq_ibz, spin, phfrq, gamma_ph, lambda_ph, displ_cart)
666 :
667 244 : do mu=1,gams%natom3
668 244 : lambda_tot = lambda_tot + lambda_ph(mu) * gams%wtq(iq_ibz)
669 : end do
670 :
671 : ! Write data to netcdf file
672 61 : if (ncid /= nctk_noid) then
673 61 : if (spin == 1) then
674 183 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qibz"), phfrq, start=[1, iq_ibz]))
675 305 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart_qibz'), displ_cart, start=[1, 1, 1, iq_ibz]))
676 : end if
677 244 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phgamma_qibz'), gamma_ph, start=[1, iq_ibz, spin]))
678 244 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phlambda_qibz'), lambda_ph, start=[1, iq_ibz, spin]))
679 : end if
680 :
681 : ! Output to the main output file
682 66 : if (iq_ibz <= max_nqpt) then
683 61 : if (gams%nsppol == 2) then
684 0 : write(msg,'(2a,3es16.6,a,i1,a,a)')ch10,&
685 0 : ' q-point =',gams%qibz(:, iq_ibz),' spin = ',spin,ch10,&
686 0 : ' Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n)'
687 : else
688 61 : write(msg,'(2a,3es16.6,a,a)')ch10,&
689 61 : ' q-point =',gams%qibz(:, iq_ibz),ch10,&
690 122 : ' Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n)'
691 : end if
692 61 : call wrtout(units, msg)
693 244 : do mu=1,gams%natom3
694 183 : write(msg,'(i5,es20.6,2es16.6)')mu, phfrq(mu), gamma_ph(mu), lambda_ph(mu)
695 244 : call wrtout(units, msg)
696 : end do
697 : else
698 0 : if (iq_ibz == max_nqpt + 1) then
699 0 : call wrtout(units, sjoin(" No more that:", itoa(max_nqpt), " q-points are reported in the main output file"))
700 0 : call wrtout(units, " Use A2f.nc and AbiPy to analyze the results.")
701 : end if
702 : end if
703 :
704 : end do
705 : ! Add blank lines to output files between spins
706 10 : call wrtout(units, "", do_flush=.True.)
707 : end do
708 :
709 5 : write(ab_out,"(a,f8.4)")" lambda= ",lambda_tot
710 : !write(ab_out,"(a,f8.4)")" omega_log= ",omega_log
711 :
712 5 : call cwtime_report(" phgamma_ncwrite:", cpu, wall, gflops)
713 :
714 5 : end subroutine phgamma_ncwrite
715 : !!***
716 :
717 : !----------------------------------------------------------------------
718 :
719 : !!****f* m_phgamma/tgamma_symm
720 : !! NAME
721 : !! tgamma_symm
722 : !!
723 : !! FUNCTION
724 : !! Symmetrize the tgamma matrix
725 : !!
726 : !! INPUTS
727 : !! qpt(3)=phonon wavevector in reduced coordinates.
728 : !! cryst<crystal_t>=Crystalline structure.
729 : !!
730 : !! SOURCE
731 :
732 14809 : subroutine tgamma_symm(cryst, qpt, tgamma)
733 :
734 : !Arguments ------------------------------------
735 : !scalars
736 : type(crystal_t),intent(in) :: cryst
737 : !arrays
738 : real(dp),intent(in) :: qpt(3)
739 : real(dp),intent(inout) :: tgamma(2,3*cryst%natom,3*cryst%natom)
740 :
741 : !Local variables-------------------------------
742 : !scalars
743 : integer :: ii,natom3,k
744 : !arrays
745 29618 : real(dp) :: tgcart(2,3*cryst%natom,3*cryst%natom)
746 14809 : real(dp) :: umat(2,3*cryst%natom,3*cryst%natom),tmp_mat(2,3*cryst%natom,3*cryst%natom)
747 : ! *********************************************************************
748 :
749 : ! Build U matrix.
750 459079 : umat = zero; k = 1
751 29618 : do ii=1,cryst%natom
752 192517 : umat(1,k:k+2, k:k+2) = cryst%gprimd
753 29618 : k = k + 3
754 : end do
755 :
756 14809 : natom3 = 3 * cryst%natom
757 :
758 : ! Reduced --> Cartesian
759 14809 : call zgemm('N', 'N', natom3, natom3, natom3, cone, tgamma, natom3, umat, natom3, czero, tmp_mat, natom3)
760 14809 : call zgemm('T', 'N', natom3, natom3, natom3, cone, umat, natom3, tmp_mat, natom3, czero, tgcart, natom3)
761 :
762 : ! Make the matrix hermitian
763 14809 : call mkherm(tgcart, 3*cryst%natom)
764 :
765 : ! Symmetrize tgamma matrix.
766 14809 : call symdyma(tgcart, cryst%indsym, cryst%natom, cryst%nsym, qpt, cryst%rprimd, cryst%symrel, cryst%symafm)
767 :
768 459079 : umat = zero; k = 1
769 29618 : do ii=0,cryst%natom-1
770 192517 : umat(1,k:k+2, k:k+2) = cryst%rprimd
771 29618 : k = k + 3
772 : end do
773 :
774 : ! Cartesian --> Reduced
775 14809 : call zgemm('N', 'N', natom3, natom3, natom3, cone, tgcart, natom3, umat, natom3, czero, tmp_mat, natom3)
776 14809 : call zgemm('T', 'N', natom3, natom3, natom3, cone, umat, natom3, tmp_mat, natom3, czero, tgamma, natom3)
777 :
778 14809 : end subroutine tgamma_symm
779 : !!***
780 :
781 : !----------------------------------------------------------------------
782 :
783 : !!****f* m_phgamma/phgamma_eval_qibz
784 : !! NAME
785 : !! phgamma_eval_qibz
786 : !!
787 : !! FUNCTION
788 : !! Compute the phonon linewidths for q-points in the IBZ without performing the interpolation.
789 : !!
790 : !! INPUTS
791 : !! cryst<crystal_t>=Crystal structure.
792 : !! ifc<ifc_type>=Interatomic force constants.
793 : !! iq_ibz=Index of the q-point in the IBZ array.
794 : !! spin=Spin index. See notes below.
795 : !!
796 : !! OUTPUT
797 : !! phfrq(gams%natom3)=Phonon frequencies
798 : !! gamma_ph(gams%natom3)=Phonon linewidths.
799 : !! lambda_ph(gams%natom3)=coupling strength coefficients
800 : !! displ_cart(2,3,cry%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates.
801 : !! [gamma_ph_ee]
802 : !!
803 : !! NOTES
804 : !! If nsppol == 1 and nspinor == 1, lambda and gamma are already summed over the two equivalent spin channels.
805 : !! If nsppol == 2, lambda and gamma are the particular contributions given by the input spin index.
806 : !! Client code is responsible for assembling the final observables by summing over spins.
807 : !!
808 : !! SOURCE
809 :
810 0 : subroutine phgamma_eval_qibz(gams, cryst, ifc, iq_ibz, spin, phfrq, gamma_ph, lambda_ph, displ_cart, gamma_ph_ee)
811 :
812 : !Arguments ------------------------------------
813 : !scalars
814 : class(phgamma_t),intent(inout) :: gams
815 : integer,intent(in) :: iq_ibz,spin
816 : type(crystal_t),intent(in) :: cryst
817 : type(ifc_type),intent(in) :: ifc
818 : !arrays
819 : real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3)
820 : real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
821 : real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3)
822 :
823 : !Local variables-------------------------------
824 : !scalars
825 : integer :: natom3, nu1, iene, jene
826 : real(dp) :: spinfact
827 : !character(len=500) :: msg
828 : !arrays
829 244 : real(dp) :: displ_red(2,gams%natom3,gams%natom3), work_qnu(gams%natom3), gam_atm(2,gams%natom3,gams%natom3)
830 : ! *************************************************************************
831 :
832 122 : natom3 = gams%natom3
833 :
834 : ! Get phonon frequencies and eigenvectors.
835 122 : call ifc%fourq(cryst, gams%qibz(:,iq_ibz), phfrq, displ_cart, out_displ_red=displ_red)
836 :
837 : ! Scalar product with the displ_red vectors.
838 : ! Note that the factor 1 / (2 * omega_qnu) coming from |g|^2 is absorbed in the expressions below.
839 488 : gam_atm = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2, natom3, natom3])
840 122 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gamma_ph)
841 :
842 122 : if (present(gamma_ph_ee) .and. gams%my_iqibz(iq_ibz) /= -1) then
843 0 : do iene = 1, gams%nene
844 0 : do jene = 1, gams%nene
845 0 : gam_atm = reshape(gams%vals_ee(:,jene,iene,:,:,iq_ibz,spin), [2, natom3, natom3])
846 0 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu)
847 0 : gamma_ph_ee(jene, iene, :) = work_qnu
848 : end do
849 : end do
850 : end if
851 :
852 : ! Compute lambda
853 : ! TODO: check this - looks like a factor of 2 wrt the inline documentation!
854 : ! NB: one factor of 2 comes from the phonon propagator and BE factor,
855 : ! then you have to be careful with the convention for the Fermi level DOS
856 : !
857 : ! spinfact should be 1 for a normal non sppol calculation without spinorbit
858 : ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
859 : ! for nsppol 2 it should be 0.5 as we have 2 spin channels to sum
860 122 : spinfact = two / (gams%nsppol * gams%nspinor)
861 :
862 488 : do nu1=1,gams%natom3
863 366 : gamma_ph(nu1) = gamma_ph(nu1) * pi * spinfact
864 366 : lambda_ph(nu1) = zero
865 366 : if (abs(phfrq(nu1)) > EPHTK_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
866 : !if (abs(phfrq(nu1)) > EPHTK_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * sum(gams%n0(:)) * phfrq(nu1)**2)
867 488 : if (present(gamma_ph_ee)) gamma_ph_ee(:,:,nu1) = gamma_ph_ee(:,:,nu1) * pi * spinfact
868 : end do
869 :
870 : ! This to avoid spurious results for the acoustic modes as gamma(q) --> 0 and lambda(q) --> 0 for q --> 0.
871 : ! Here we set everything to zero when we are inside a sphere of radius EPH_Q0TOL as this acoustic rule
872 : ! is not fulfilled due to numerical inaccuracies.
873 122 : if (normv(gams%qibz(:,iq_ibz), cryst%gmet, "G") < EPH_Q0TOL) then
874 40 : gamma_ph(1:3) = zero
875 40 : lambda_ph(1:3) = zero
876 : end if
877 :
878 122 : end subroutine phgamma_eval_qibz
879 : !!***
880 :
881 : !----------------------------------------------------------------------
882 :
883 : !!****f* m_phgamma/phgamma_interp
884 : !! NAME
885 : !! phgamma_interp
886 : !!
887 : !! FUNCTION
888 : !! Interpolate the phonon linewidths at a given q-point.
889 : !!
890 : !! INPUTS
891 : !! cryst<crystal_t>=crystalline structure.
892 : !! ifc<ifc_type>=Interatomic force constants.
893 : !! spin=Spin index
894 : !! qpt(3)=q-point in reduced coordinates
895 : !! gamma_ph(3*natom)=Phonon linewidths
896 : !!
897 : !! OUTPUT
898 : !! gamma_ph(gams%natom3)=Interpolated Phonon linewidths.
899 : !! lamda_ph(3*natom)=Lambda coefficients for the different phonon modes.
900 : !! phfrq(3*natom)=phonon frequencies at current q
901 : !! displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates
902 : !!
903 : !! SOURCE
904 :
905 0 : subroutine phgamma_interp(gams, cryst, ifc, spin, qpt, phfrq, gamma_ph, lambda_ph, displ_cart, gamma_ph_ee)
906 :
907 : !Arguments ------------------------------------
908 : !scalars
909 : class(phgamma_t),intent(inout) :: gams
910 : integer,intent(in) :: spin
911 : type(crystal_t),intent(in) :: cryst
912 : type(ifc_type),intent(in) :: ifc
913 : !arrays
914 : real(dp),intent(in) :: qpt(3)
915 : real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3)
916 : real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
917 : real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3)
918 :
919 : !Local variables-------------------------------
920 : !scalars
921 : integer,parameter :: qtor0 = 0
922 : integer, save :: icall=0
923 : integer :: natom3,nu1
924 : real(dp) :: spinfact
925 : character(len=500) :: msg
926 : !arrays
927 1860 : real(dp) :: displ_red(2,gams%natom3,gams%natom3), gam_now(2,gams%natom3**2), gam_atm(2,gams%natom3,gams%natom3)
928 930 : real(dp),allocatable :: coskr(:,:),sinkr(:,:)
929 : ! *************************************************************************
930 :
931 : ! Compute internal tables used for Fourier interpolation.
932 930 : if (.not. allocated(gams%vals_rpt)) call gams%interp_setup(cryst)
933 :
934 930 : if (present(gamma_ph_ee) .and. icall == 0) then
935 0 : gamma_ph_ee = zero
936 : write (msg,'(3a)')&
937 0 : " For the moment gams_ee matrix elements are not FT interpolated wrt q,",ch10,&
938 0 : " only evaluated on the electron k grid. The resulting a2feew will be 0"
939 0 : ABI_WARNING(msg)
940 0 : icall = 1
941 : end if
942 :
943 930 : natom3 = gams%natom3
944 :
945 : ! Taken from mkph_linwid
946 : ! This reduced version of ftgkk supposes the kpoints have been integrated
947 : ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt.
948 2790 : ABI_MALLOC(coskr, (1, gams%nrpt))
949 1860 : ABI_MALLOC(sinkr, (1, gams%nrpt))
950 : ! TODO: This is not optimal
951 930 : call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr)
952 :
953 930 : call ftgam(gams%wghatm, gam_now, gams%vals_rpt(:,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr)
954 :
955 : ! This call is not executed in elphon!
956 930 : if (gams%symgamma == 1) call tgamma_symm(cryst, qpt, gam_now)
957 :
958 930 : ABI_FREE(coskr)
959 930 : ABI_FREE(sinkr)
960 :
961 : ! Get phonon frequencies and eigenvectors.
962 930 : call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
963 :
964 : ! Scalar product with the displ_red
965 3720 : gam_atm = reshape(gam_now, [2, natom3, natom3])
966 930 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gamma_ph)
967 :
968 : ! Compute lambda
969 : ! spinfact should be 1 for a normal non sppol calculation without spinorbit
970 : ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
971 : ! for nsppol 2 it should be 0.5 as we have 2 spin channels to sum
972 930 : spinfact = two / (gams%nsppol * gams%nspinor)
973 :
974 : ! Compute lambda
975 3720 : do nu1=1,gams%natom3
976 2790 : gamma_ph(nu1) = gamma_ph(nu1) * pi * spinfact
977 2790 : lambda_ph(nu1) = zero
978 3720 : if (abs(phfrq(nu1)) > EPHTK_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
979 : end do
980 :
981 : ! This to avoid spurious results for the acoustic modes.
982 : ! In principle, gamma(q) --> 0 and lambda(q) --> 0 for q --> 0
983 : ! but the Fourier interpolated gammas do not fulfill this property so we set everything
984 : ! to zero when we are inside a sphere or radius
985 930 : if (normv(qpt, cryst%gmet, "G") < EPH_Q0TOL) then
986 : !write(std_out,*)"Setting values to zero."
987 40 : gamma_ph(1:3) = zero
988 40 : lambda_ph(1:3) = zero
989 : end if
990 :
991 930 : end subroutine phgamma_interp
992 : !!***
993 :
994 : !----------------------------------------------------------------------
995 :
996 : !!****f* m_phgamma/phgamma_interp_setup
997 : !! NAME
998 : !! phgamma_interp_setup
999 : !!
1000 : !! FUNCTION
1001 : !! This routines prepares the internal tables used to interpolate the linewidths in q-space
1002 : !!
1003 : !! INPUTS
1004 : !!
1005 : !! SOURCE
1006 :
1007 5 : subroutine phgamma_interp_setup(gams, cryst)
1008 :
1009 : !Arguments ------------------------------------
1010 : !scalars
1011 : class(phgamma_t),intent(inout) :: gams
1012 : type(crystal_t),intent(in) :: cryst
1013 :
1014 : !Local variables-------------------------------
1015 : !scalars
1016 : integer,parameter :: qtor1 = 1
1017 : integer :: iq_bz,iq_ibz,spin,ierr,ii
1018 : !character(len=500) :: msg
1019 : !arrays
1020 5 : integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:)
1021 5 : real(dp),allocatable :: coskr(:,:),sinkr(:,:), gamma_qpt(:,:,:,:),atmfrc(:,:), vals_bz(:,:,:,:)
1022 : ! *************************************************************************
1023 :
1024 25 : ABI_MALLOC_OR_DIE(vals_bz, (2, gams%natom3**2, gams%nqbz, gams%nsppol), ierr)
1025 21514 : vals_bz = zero
1026 :
1027 : ! Build tables needed by complete_gamma.
1028 5 : call ephtk_mkqtabs(cryst, gams%nqibz, gams%qibz, gams%nqbz, gams%qbz, qirredtofull, qpttoqpt)
1029 :
1030 : ! Fill BZ array with IBZ data.
1031 10 : do spin=1,gams%nsppol
1032 71 : do iq_ibz=1,gams%nqibz
1033 61 : iq_bz = qirredtofull(iq_ibz)
1034 188 : vals_bz(:,:,iq_bz,spin) = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2, gams%natom3**2])
1035 : end do
1036 : end do
1037 :
1038 : ! Complete vals_bz in the full BZ.
1039 : ! FIXME: Change complete_gamma API to pass (..., nsppol)
1040 25 : ABI_MALLOC(gamma_qpt, (2, gams%natom3**2, gams%nsppol, gams%nqbz))
1041 10 : do spin=1,gams%nsppol
1042 21514 : gamma_qpt(:, :, spin, :) = vals_bz(:, :, :, spin)
1043 : end do
1044 :
1045 : call complete_gamma(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz, &
1046 5 : gams%eph_scalprod, qirredtofull, qpttoqpt, gamma_qpt)
1047 :
1048 10 : do spin=1,gams%nsppol
1049 21514 : vals_bz(:, :, :, spin) = gamma_qpt(:, :, spin, :)
1050 : end do
1051 :
1052 5 : ABI_FREE(gamma_qpt)
1053 5 : ABI_FREE(qirredtofull)
1054 5 : ABI_FREE(qpttoqpt)
1055 :
1056 : ! This call is not executed in elphon!
1057 5 : if (gams%symgamma == 1) then
1058 10 : do spin=1,gams%nsppol
1059 778 : do iq_bz=1,gams%nqbz
1060 773 : call tgamma_symm(cryst, gams%qbz(:,iq_bz), vals_bz(:,:,iq_bz,spin))
1061 : end do
1062 : end do
1063 : end if
1064 :
1065 : ! Now FT to real space too
1066 : ! NOTE: gprim (not gprimd) is used for all FT interpolations,
1067 : ! to be consistent with the dimensions of the rpt, which come from anaddb.
1068 25 : ABI_MALLOC_OR_DIE(gams%vals_rpt, (2, gams%natom3**2, gams%nrpt, gams%nsppol), ierr)
1069 13030 : gams%vals_rpt = zero
1070 :
1071 : ! q --> r
1072 20 : ABI_MALLOC(coskr, (gams%nqbz, gams%nrpt))
1073 15 : ABI_MALLOC(sinkr, (gams%nqbz, gams%nrpt))
1074 5 : call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr)
1075 :
1076 10 : do spin=1,gams%nsppol
1077 : call ftgam(gams%wghatm, vals_bz(:,:,:,spin), gams%vals_rpt(:,:,:,spin), gams%natom, gams%nqbz,&
1078 5 : gams%nrpt, qtor1, coskr, sinkr)
1079 :
1080 : ! Enforce "acoustic" rule on vals_rpt. NB: This call is not executed in elphon!
1081 10 : if (gams%asr /= 0) then
1082 0 : ABI_MALLOC(atmfrc, (3*gams%natom*3*gams%natom, gams%nrpt))
1083 0 : do ii=1,2
1084 0 : atmfrc = gams%vals_rpt(ii,:,:,spin)
1085 : !gals%vals_rpt(2,:,,spin) = zero
1086 0 : call asrif9(gams%asr, atmfrc, gams%natom, gams%nrpt, gams%rpt, gams%wghatm)
1087 0 : gams%vals_rpt(ii,:,:,spin) = atmfrc
1088 : end do
1089 0 : ABI_FREE(atmfrc)
1090 : end if
1091 : end do
1092 :
1093 5 : ABI_FREE(vals_bz)
1094 5 : ABI_FREE(coskr)
1095 5 : ABI_FREE(sinkr)
1096 :
1097 5 : end subroutine phgamma_interp_setup
1098 : !!***
1099 :
1100 : !----------------------------------------------------------------------
1101 :
1102 : !!****f* m_phgamma/phgamma_vv_eval_qibz
1103 : !! NAME
1104 : !! phgamma_vv_eval_qibz
1105 : !!
1106 : !! FUNCTION
1107 : !! Compute the phonon linewidths times velocity squared, for q-points in the IBZ without performing interpolation.
1108 : !!
1109 : !! INPUTS
1110 : !! gams<phgamma_t>
1111 : !! cryst<crystal_t>=Crystal structure.
1112 : !! ifc<ifc_type>=Interatomic force constants.
1113 : !! iq_ibz=Index of the q-point in the IBZ array.
1114 : !! spin=Spin index
1115 : !!
1116 : !! OUTPUT
1117 : !! phfrq(gams%natom3)=Phonon frequencies
1118 : !! gamma_in_ph(gams%natom3)=Phonon linewidths.
1119 : !! gamma_out_ph(gams%natom3)=Phonon linewidths.
1120 : !! lambda_in_ph(gams%natom3)=Phonon linewidths.
1121 : !! lambda_out_ph(gams%natom3)=Phonon linewidths.
1122 : !!
1123 : !! SOURCE
1124 :
1125 61 : subroutine phgamma_vv_eval_qibz(gams, cryst, ifc, iq_ibz, spin, phfrq, gamma_in_ph, gamma_out_ph, lambda_in_ph, lambda_out_ph)
1126 :
1127 : !Arguments ------------------------------------
1128 : !scalars
1129 : class(phgamma_t),intent(inout) :: gams
1130 : integer,intent(in) :: iq_ibz,spin
1131 : type(crystal_t),intent(in) :: cryst
1132 : type(ifc_type),intent(in) :: ifc
1133 : !arrays
1134 : real(dp),intent(out) :: phfrq(gams%natom3)
1135 : real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3)
1136 : real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3)
1137 :
1138 : !Local variables-------------------------------
1139 : !scalars
1140 : integer,parameter :: qtor0 = 0
1141 : integer :: natom3, nu1, idir, jdir, ii
1142 : real(dp) :: spinfact
1143 : !character(len=500) :: msg
1144 : !arrays
1145 122 : real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom), displ_red(2,gams%natom3,gams%natom3)
1146 122 : real(dp) :: work_qnu(gams%natom3), gam_atm(2,gams%natom3,gams%natom3)
1147 : ! *************************************************************************
1148 :
1149 61 : natom3 = gams%natom3
1150 :
1151 : ! Get phonon frequencies and eigenvectors.
1152 61 : call ifc%fourq(cryst, gams%qibz(:,iq_ibz), phfrq, displ_cart, out_displ_red=displ_red)
1153 :
1154 244 : do jdir=1,3
1155 793 : do idir=1,3
1156 549 : ii = idir + gams%ndir_transp * (jdir - 1)
1157 :
1158 : ! Scalar product with the displ_red vectors.
1159 2196 : gam_atm = reshape(gams%vals_in_qibz(:,ii,:,:,iq_ibz,spin), [2, natom3, natom3])
1160 549 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu)
1161 2196 : gamma_in_ph(idir,jdir,:) = work_qnu
1162 :
1163 2196 : gam_atm = reshape(gams%vals_out_qibz(:,ii,:,:,iq_ibz,spin), [2, natom3, natom3])
1164 549 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu)
1165 2379 : gamma_out_ph(idir,jdir,:) = work_qnu
1166 : end do ! idir
1167 : end do ! jdir
1168 :
1169 : ! TODO : check this - looks like a factor of 2 wrt the inline documentation!
1170 : !spinfact should be 1 for a normal non sppol calculation without spinorbit
1171 : !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1172 : !for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1173 61 : spinfact = two / (gams%nsppol*gams%nspinor)
1174 :
1175 : ! Compute lambda transport.
1176 244 : do jdir =1,3
1177 793 : do idir =1,3
1178 2379 : do nu1=1,gams%natom3
1179 1647 : gamma_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) * pi * spinfact
1180 1647 : gamma_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) * pi * spinfact
1181 1647 : lambda_in_ph(idir,jdir,nu1) = zero
1182 1647 : lambda_out_ph(idir,jdir,nu1) = zero
1183 2196 : if (abs(phfrq(nu1)) > EPHTK_WTOL) then
1184 1512 : lambda_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1185 1512 : lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1186 : end if
1187 : end do
1188 : end do ! idir
1189 : end do ! jdir
1190 :
1191 61 : end subroutine phgamma_vv_eval_qibz
1192 : !!***
1193 :
1194 : !----------------------------------------------------------------------
1195 :
1196 : !!****f* m_phgamma/phgamma_vv_interp
1197 : !! NAME
1198 : !! phgamma_vv_interp
1199 : !!
1200 : !! FUNCTION
1201 : !! Interpolate the linewidths at a given q-point.
1202 : !!
1203 : !! INPUTS
1204 : !! cryst<crystal_t>=crystalline structure.
1205 : !! ifc<ifc_type>=Interatomic force constants.
1206 : !! spin=Spin index
1207 : !! qpt(3)=q-point in reduced coordinates
1208 : !! gamma_ph(3*natom)=Phonon linewidths
1209 : !! displ_cart(2,3,cryst%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates.
1210 : !!
1211 : !! OUTPUT
1212 : !! phfrq(3*natom)=phonon frequencies at current q
1213 : !! gamma_out_ph(gams%natom3)=Interpolated Phonon linewidths.
1214 : !! lambda_in_ph(3*natom)=Lambda coefficients for the different phonon modes.
1215 : !! lambda_out_ph(3*natom)=Lambda coefficients for the different phonon modes.
1216 : !!
1217 : !! SOURCE
1218 :
1219 725 : subroutine phgamma_vv_interp(gams, cryst, ifc, spin, qpt, phfrq, gamma_in_ph, gamma_out_ph, lambda_in_ph, lambda_out_ph)
1220 :
1221 : !Arguments ------------------------------------
1222 : !scalars
1223 : class(phgamma_t),intent(inout) :: gams
1224 : type(crystal_t),intent(in) :: cryst
1225 : type(ifc_type),intent(in) :: ifc
1226 : integer,intent(in) :: spin
1227 : !arrays
1228 : real(dp),intent(in) :: qpt(3)
1229 : real(dp),intent(out) :: phfrq(gams%natom3)
1230 : real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3)
1231 : real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3)
1232 :
1233 : !Local variables-------------------------------
1234 : !scalars
1235 : integer,parameter :: qtor0 = 0
1236 : integer :: natom3, nu1, idir,jdir,ii
1237 : real(dp) :: spinfact
1238 : !character(len=500) :: msg
1239 : !arrays
1240 1450 : real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom), displ_red(2,gams%natom3,gams%natom3),work_qnu(gams%natom3)
1241 1450 : real(dp) :: gam_in_now(2,3,3,gams%natom3**2), gam_out_now(2,3,3,gams%natom3**2), gam_atm(2,gams%natom3,gams%natom3)
1242 725 : real(dp),allocatable :: coskr(:,:),sinkr(:,:)
1243 : ! *************************************************************************
1244 :
1245 : ! Compute internal tables used for Fourier interpolation.
1246 725 : if (.not.allocated(gams%vals_in_rpt)) call phgamma_vv_interp_setup(gams, cryst)
1247 :
1248 725 : natom3 = gams%natom3
1249 :
1250 : ! Get phonon frequencies and eigenvectors.
1251 725 : call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
1252 :
1253 : ! Taken from mkph_linwid
1254 : ! This reduced version of ftgkk supposes the kpoints have been integrated
1255 : ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt.
1256 2175 : ABI_MALLOC(coskr, (1,gams%nrpt))
1257 1450 : ABI_MALLOC(sinkr, (1,gams%nrpt))
1258 : ! TODO: This is not optimal
1259 725 : call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr)
1260 :
1261 2900 : do idir=1,3
1262 9425 : do jdir=1,3
1263 6525 : ii = idir+gams%ndir_transp*(jdir-1)
1264 :
1265 33988725 : call ftgam(gams%wghatm, gam_in_now, gams%vals_in_rpt(:,ii,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr)
1266 33988725 : call ftgam(gams%wghatm, gam_out_now, gams%vals_out_rpt(:,ii,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr)
1267 :
1268 : ! This call is not executed in elphon!
1269 : ! TODO: needs to take into account the matrix nature of _in_ and _out_ gammas
1270 6525 : if (gams%symgamma == 1) then
1271 6525 : call tgamma_symm(cryst, qpt, gam_in_now)
1272 6525 : call tgamma_symm(cryst, qpt, gam_out_now)
1273 : end if
1274 :
1275 : ! Scalar product with the displ_red
1276 26100 : gam_atm = reshape(gam_in_now, [2, natom3, natom3])
1277 6525 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu)
1278 26100 : gamma_in_ph(idir,jdir,:) = work_qnu
1279 :
1280 26100 : gam_atm = reshape(gam_out_now, [2, natom3, natom3])
1281 6525 : call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu)
1282 28275 : gamma_out_ph(idir,jdir,:) = work_qnu
1283 : end do
1284 : end do
1285 :
1286 725 : ABI_FREE(coskr)
1287 725 : ABI_FREE(sinkr)
1288 :
1289 : ! Compute lambda
1290 : ! spinfact should be 1 for a normal non sppol calculation without spinorbit
1291 : ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2
1292 : ! for sppol 2 it should be 0.5 as we have 2 spin channels to sum
1293 725 : spinfact = two / (gams%nsppol * gams%nspinor)
1294 :
1295 : ! Compute lambda
1296 2900 : do nu1=1,gams%natom3
1297 9425 : do idir=1,3
1298 28275 : do jdir=1,3
1299 19575 : gamma_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) * pi * spinfact
1300 19575 : gamma_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) * pi * spinfact
1301 19575 : lambda_in_ph(idir,jdir,nu1) = zero
1302 19575 : lambda_out_ph(idir,jdir,nu1) = zero
1303 26100 : if (abs(phfrq(nu1)) > EPHTK_WTOL) then
1304 19440 : lambda_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1305 19440 : lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2)
1306 : end if
1307 : end do
1308 : end do
1309 : end do
1310 :
1311 725 : end subroutine phgamma_vv_interp
1312 : !!***
1313 :
1314 : !----------------------------------------------------------------------
1315 :
1316 : !!****f* m_phgamma/phgamma_vv_interp_setup
1317 : !! NAME
1318 : !! phgamma_vv_interp_setup
1319 : !!
1320 : !! FUNCTION
1321 : !! This routines prepare the internal tables used to interpolate the vv_linewidths in q-space
1322 : !!
1323 : !! INPUTS
1324 : !! action =
1325 : !! "INIT" to allocate and compute the internal tables (default)
1326 : !! "FREE" to deallocate the internal tables.
1327 : !!
1328 : !! SOURCE
1329 :
1330 5 : subroutine phgamma_vv_interp_setup(gams, cryst)
1331 :
1332 : !Arguments ------------------------------------
1333 : class(phgamma_t),intent(inout) :: gams
1334 : type(crystal_t),intent(in) :: cryst
1335 :
1336 : !Local variables-------------------------------
1337 : !scalars
1338 : integer,parameter :: qtor1 = 1
1339 : integer :: iq_bz,iq_ibz,spin,ierr, ii, idir, jdir
1340 : !character(len=500) :: msg
1341 : !arrays
1342 5 : integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:)
1343 5 : real(dp),allocatable :: coskr(:,:),sinkr(:,:), vals_in_bz(:,:,:,:,:), vals_out_bz(:,:,:,:,:)
1344 : ! *************************************************************************
1345 :
1346 25 : ABI_MALLOC_OR_DIE(vals_in_bz,(2, 9, gams%natom3**2, gams%nqbz, gams%nsppol), ierr)
1347 25 : ABI_MALLOC_OR_DIE(vals_out_bz,(2, 9, gams%natom3**2, gams%nqbz, gams%nsppol), ierr)
1348 194314 : vals_in_bz = zero
1349 194314 : vals_out_bz = zero
1350 :
1351 : ! Build tables needed by complete_gamma.
1352 5 : call ephtk_mkqtabs(cryst, gams%nqibz, gams%qibz, gams%nqbz, gams%qbz, qirredtofull, qpttoqpt)
1353 :
1354 : ! Fill BZ array with IBZ data.
1355 10 : do spin=1,gams%nsppol
1356 71 : do iq_ibz=1,gams%nqibz
1357 61 : iq_bz = qirredtofull(iq_ibz)
1358 244 : vals_in_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_in_qibz(:,:,:,:,iq_ibz,spin), [2, 9, gams%natom3**2])
1359 249 : vals_out_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_out_qibz(:,:,:,:,iq_ibz,spin), [2, 9, gams%natom3**2])
1360 : end do
1361 : end do
1362 :
1363 : ! Complete vals_bz in the full BZ.
1364 : !TODO!!! rotate the vv in and out matrices, according to the symmetry operation, instead of just copying them
1365 : !call complete_gamma_vv(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz,&
1366 : ! gams%eph_scalprod, qirredtofull, qpttoqpt, vals_in_bz(:,:,:,:,spin))
1367 : !call complete_gamma_vv(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz,&
1368 : ! gams%eph_scalprod, qirredtofull, qpttoqpt, vals_out_bz(:,:,:,:,spin))
1369 :
1370 : ! TODO: replace the above with these calls from anaddb
1371 : !call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,&
1372 : !& elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_trout,elph_ds%qirredtofull,qpttoqpt)
1373 :
1374 : ! TODO: idem for vv_vals 3x3 matrices
1375 :
1376 5 : ABI_FREE(qirredtofull)
1377 5 : ABI_FREE(qpttoqpt)
1378 :
1379 : !! This call is not executed in elphon!
1380 : !if (gams%symgamma == 1) then
1381 : ! do spin=1,gams%nsppol
1382 : ! do iq_bz=1,gams%nqbz
1383 : ! call tgamma_symm_vv(cryst, gams%qbz(:,iq_bz), vals_in_bz(:,:,:,iq_bz,spin))
1384 : ! call tgamma_symm_vv(cryst, gams%qbz(:,iq_bz), vals_out_bz(:,:,:,iq_bz,spin))
1385 : ! end do
1386 : ! end do
1387 : !end if
1388 :
1389 : ! Now FT to real space
1390 : ! NOTE: gprim (not gprimd) is used for all FT interpolations,
1391 : ! to be consistent with the dimensions of the rpt, which come from anaddb.
1392 : ! TODO: this is needed only if FT is used, not when the linear interpolation is employed.
1393 : !if (.not. allocated(gams%vals_in_rpt)) then
1394 25 : ABI_MALLOC_OR_DIE(gams%vals_in_rpt, (2, 9, gams%natom3**2, gams%nrpt, gams%nsppol), ierr)
1395 25 : ABI_MALLOC_OR_DIE(gams%vals_out_rpt, (2, 9, gams%natom3**2, gams%nrpt, gams%nsppol), ierr)
1396 117655 : gams%vals_in_rpt = zero
1397 117655 : gams%vals_out_rpt = zero
1398 :
1399 : ! q --> r
1400 20 : ABI_MALLOC(coskr, (gams%nqbz,gams%nrpt))
1401 15 : ABI_MALLOC(sinkr, (gams%nqbz,gams%nrpt))
1402 5 : call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr)
1403 :
1404 10 : do spin=1,gams%nsppol
1405 25 : do idir=1,3
1406 65 : do jdir=1,3
1407 45 : ii = idir+gams%ndir_transp*(jdir-1)
1408 : ! TODO: this is no contiguous in memory and will be slow. Make adapted ftgam?
1409 : call ftgam(gams%wghatm, vals_in_bz(:,ii,:,:,spin), gams%vals_in_rpt(:,ii,:,:,spin), gams%natom, gams%nqbz,&
1410 621477 : gams%nrpt, qtor1, coskr, sinkr)
1411 : call ftgam(gams%wghatm, vals_out_bz(:,ii,:,:,spin), gams%vals_out_rpt(:,ii,:,:,spin), gams%natom, gams%nqbz,&
1412 621492 : gams%nrpt, qtor1, coskr, sinkr)
1413 : end do
1414 : end do
1415 : end do
1416 :
1417 5 : ABI_FREE(vals_in_bz)
1418 5 : ABI_FREE(vals_out_bz)
1419 5 : ABI_FREE(coskr)
1420 5 : ABI_FREE(sinkr)
1421 :
1422 5 : end subroutine phgamma_vv_interp_setup
1423 : !!***
1424 : !----------------------------------------------------------------------
1425 :
1426 : !!****f* m_phgamma/phgamma_linwid
1427 : !! NAME
1428 : !! phgamma_linwid
1429 : !!
1430 : !! FUNCTION
1431 : !! Interpolate the phonon linewidths along an arbitrary q-path (use Fourier interpolation).
1432 : !!
1433 : !! INPUTS
1434 : !! cryst<crystal_t>=Info on the unit cell and symmetries.
1435 : !! ifc<ifc_type>=Interatomic force constants.
1436 : !! ndivsm=Number of points used to sample the smallest segment
1437 : !! nvert = Number of extrema in qverts
1438 : !! qverts(3,nvert) = vertices of reciprocal space trajectory
1439 : !! basename=name used to create the different output files (text format).
1440 : !! ncid=Netcdf file handler (already open in the caller).
1441 : !! comm=MPI communicator
1442 : !!
1443 : !! OUTPUT
1444 : !! wminmax=Minimum and max phonon frequency obtained on the path (Hartree units)
1445 : !!
1446 : !! SOURCE
1447 :
1448 5 : subroutine phgamma_linwid(gams, cryst, ifc, ndivsm, nvert, qverts, basename, ncid, wminmax, comm)
1449 :
1450 : !Arguments ------------------------------------
1451 : !scalars
1452 : class(phgamma_t),intent(inout) :: gams
1453 : integer,intent(in) :: nvert,ndivsm,comm,ncid
1454 : type(crystal_t),intent(in) :: cryst
1455 : type(ifc_type),intent(in) :: ifc
1456 : character(len=*),intent(in) :: basename
1457 : !arrays
1458 : real(dp),intent(in) :: qverts(3,nvert)
1459 : real(dp),intent(out) :: wminmax(2)
1460 :
1461 : !Local variables-------------------------------
1462 : !scalars
1463 : integer,parameter :: master = 0
1464 : integer :: natom,ii,mu,iqpt,natom3,nsppol,ierr, spin,unt,nqpt,nrpt,cnt,nproc,my_rank, ncerr
1465 : real(dp) :: omega_min,omega_max,wtmp,omega
1466 : character(len=500) :: msg
1467 5 : type(kpath_t) :: qpath
1468 : !arrays
1469 10 : real(dp) :: gamma_spin(gams%nsppol),lambda_spin(gams%nsppol)
1470 10 : real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom)
1471 10 : real(dp) :: phfrq(3*cryst%natom),gamma_ph(3*cryst%natom),lambda_ph(3*cryst%natom), qpt(3),shift(3)
1472 5 : real(dp),allocatable :: all_phfreq(:,:),all_gammaq(:,:,:),all_lambdaq(:,:,:),all_displ_cart(:,:,:,:)
1473 : ! *********************************************************************
1474 :
1475 5 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1476 5 : natom = cryst%natom; natom3 = gams%natom3; nsppol = gams%nsppol; nrpt = gams%nrpt
1477 :
1478 : ! Define the q-path along which phonon linwid will be interpolated.
1479 5 : call qpath%init(qverts, cryst%gprimd, ndivsm)
1480 5 : nqpt = qpath%npts
1481 :
1482 : ! Allocate workspace arrays for MPI.
1483 840 : ABI_CALLOC(all_phfreq, (natom3, nqpt))
1484 850 : ABI_CALLOC(all_gammaq, (natom3, nqpt, nsppol))
1485 845 : ABI_CALLOC(all_lambdaq, (natom3, nqpt, nsppol))
1486 6380 : ABI_CALLOC(all_displ_cart, (2, natom3, natom3, nqpt))
1487 :
1488 : ! initialize the minimum and maximum phonon frequency
1489 5 : omega_min = huge(one); omega_max = -huge(one)
1490 :
1491 : ! Interpolation along specified path in q space (keep spin dep.)
1492 5 : cnt = 0
1493 10 : do spin=1,nsppol
1494 215 : do iqpt=1,nqpt
1495 205 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
1496 820 : call wrap2_pmhalf(qpath%points(:,iqpt), qpt, shift)
1497 :
1498 : ! Get phgamma
1499 205 : call gams%interp(cryst, ifc, spin, qpt, phfrq, gamma_ph, lambda_ph, displ_cart)
1500 820 : all_gammaq(:, iqpt, spin) = gamma_ph
1501 820 : all_lambdaq(:, iqpt, spin) = lambda_ph
1502 205 : if (spin == 1) then
1503 820 : all_phfreq(:, iqpt) = phfrq
1504 6355 : all_displ_cart(:, :, :, iqpt) = displ_cart
1505 : end if
1506 :
1507 : ! Find max/min phonon frequency along path chosen
1508 : ! presumed to be representative of full BZ to within 10 percent
1509 205 : omega_min = min(omega_min, phfrq(1))
1510 210 : omega_max = max(omega_max, phfrq(natom3))
1511 : end do ! end iqpt do
1512 : end do ! spin
1513 :
1514 : ! Collect results
1515 5 : if (omega_min > tol12) omega_min = zero
1516 5 : wtmp = omega_min; call xmpi_min(wtmp, omega_min, comm, ierr)
1517 5 : wtmp = omega_max; call xmpi_max(wtmp, omega_max, comm, ierr)
1518 15 : wminmax = [omega_min, omega_max]
1519 :
1520 5 : call xmpi_sum_master(all_gammaq, master, comm, ierr)
1521 5 : call xmpi_sum_master(all_lambdaq, master, comm, ierr)
1522 5 : call xmpi_sum_master(all_phfreq, master, comm, ierr)
1523 5 : call xmpi_sum_master(all_displ_cart, master, comm, ierr)
1524 :
1525 : ! Master writes text file with final results.
1526 : ! 3 * natom blocks, one block for each phonon mode.
1527 : ! Each block contains:
1528 : ! iqpt omega_(q) gamma(q) lambda(q) nesting(q) ....
1529 :
1530 5 : if (xmpi_comm_rank(comm) == master) then
1531 5 : if (open_file(strcat(basename, '_PHGAMMA'), msg, newunit=unt, action="write", form="formatted", status="unknown") /= 0) then
1532 5 : ABI_ERROR(msg)
1533 : end if
1534 :
1535 5 : write(unt,'(a)') '#'
1536 5 : write(unt,'(a)') '# ABINIT package: E-PH band structure file. Hartree units'
1537 5 : write(unt,'(a)') '#'
1538 5 : write(unt,'(a,i0,a)')'# Phonon frequencies, ph linewidths and lambda calculated on ',nqpt,' q-points'
1539 10 : call qpath%print([unt], header="Description of the q-path:", pre="#")
1540 15 : do ii=1,2; write(unt,'(a)') "# "; end do
1541 :
1542 10 : write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(gams%n0)
1543 10 : do spin=1,nsppol
1544 10 : write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",gams%n0(spin)
1545 : end do
1546 :
1547 20 : do mu=1,natom3
1548 15 : write(unt,'(a)')"#"
1549 15 : if (nsppol == 1) write(unt,'(a,i0,a)')"# phonon mode [",mu,"] q-index omega gamma lambda"
1550 15 : if (nsppol == 2) write(unt,'(a,i0,a)')&
1551 0 : "# phonon mode [",mu,"] q-index omega gamma_tot lambda_tot gamma[spin=1] lambda[spin=1] gamma[2] lambda[2]"
1552 15 : write(unt,'(a)')"#"
1553 635 : do iqpt=1,nqpt
1554 615 : omega = all_phfreq(mu, iqpt)
1555 1230 : gamma_spin = all_gammaq(mu, iqpt, :)
1556 1230 : lambda_spin = all_lambdaq(mu, iqpt, :)
1557 630 : if (nsppol == 1) then
1558 615 : write(unt,'(i8,3es16.6)' )iqpt,omega,gamma_spin(1),lambda_spin(1)
1559 : else
1560 0 : write(unt,'(i8,es20.6,6es16.6)' )iqpt,omega,&
1561 0 : sum(gamma_spin),sum(lambda_spin),gamma_spin(1),lambda_spin(1),gamma_spin(2),lambda_spin(2)
1562 : end if
1563 : end do
1564 : end do
1565 :
1566 5 : close(unt)
1567 :
1568 : ! Write data to netcdf file
1569 5 : if (ncid /= nctk_noid) then
1570 : ncerr = nctk_def_dims(ncid, [&
1571 : nctkdim_t("natom3", 3*natom), nctkdim_t("nqpath", nqpt), nctkdim_t("number_of_spins", nsppol) &
1572 20 : ], defmode=.True.)
1573 5 : NCF_CHECK(ncerr)
1574 :
1575 : ncerr = nctk_def_arrays(ncid, [&
1576 : nctkarr_t('qpath', "dp", "number_of_reduced_dimensions, nqpath"), &
1577 : nctkarr_t('phfreq_qpath', "dp", "natom3, nqpath"), &
1578 : nctkarr_t('phdispl_cart_qpath', "dp", "two, natom3, natom3, nqpath"), &
1579 : nctkarr_t('phgamma_qpath', "dp", "natom3, nqpath, number_of_spins"), &
1580 30 : nctkarr_t('phlambda_qpath', "dp", "natom3, nqpath, number_of_spins")])
1581 5 : NCF_CHECK(ncerr)
1582 :
1583 5 : NCF_CHECK(nctk_set_datamode(ncid))
1584 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpath"), qpath%points))
1585 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qpath"), all_phfreq))
1586 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phdispl_cart_qpath"), all_displ_cart))
1587 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phgamma_qpath"), all_gammaq))
1588 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phlambda_qpath"), all_lambdaq))
1589 : end if
1590 : end if ! master
1591 :
1592 5 : ABI_FREE(all_phfreq)
1593 5 : ABI_FREE(all_gammaq)
1594 5 : ABI_FREE(all_lambdaq)
1595 5 : ABI_FREE(all_displ_cart)
1596 :
1597 5 : call qpath%free()
1598 :
1599 10 : end subroutine phgamma_linwid
1600 : !!***
1601 :
1602 : !----------------------------------------------------------------------
1603 :
1604 : !!****f* m_phgamma/a2fw_free
1605 : !! NAME
1606 : !! a2fw_free
1607 : !!
1608 : !! FUNCTION
1609 : !! Free the memory allocated in a2f
1610 : !!
1611 : !! SOURCE
1612 :
1613 10 : subroutine a2fw_free(a2f)
1614 :
1615 : !Arguments ------------------------------------
1616 : class(a2fw_t),intent(inout) :: a2f
1617 : ! *********************************************************************
1618 :
1619 : ! integer
1620 10 : ABI_SFREE(a2f%qshift)
1621 :
1622 : ! real
1623 10 : ABI_SFREE(a2f%n0)
1624 10 : ABI_SFREE(a2f%omega)
1625 10 : ABI_SFREE(a2f%vals)
1626 10 : ABI_SFREE(a2f%vals_ee)
1627 10 : ABI_SFREE(a2f%lambdaw)
1628 :
1629 10 : end subroutine a2fw_free
1630 : !!***
1631 :
1632 : !----------------------------------------------------------------------
1633 :
1634 : !!****f* m_phgamma/a2fw_init
1635 : !! NAME
1636 : !! a2fw_init
1637 : !!
1638 : !! FUNCTION
1639 : !! Calculates the FS averaged alpha^2F(w) function
1640 : !!
1641 : !! INPUTS
1642 : !! cryst<crystal_t>=Info on the unit cell.
1643 : !! ifc<ifc_type>=Interatomic force constants.
1644 : !! gams<phgamma_t>=Structure storing the phonon linewidths.
1645 : !! wstep=Step for linear frequency mesh in Ha.
1646 : !! wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w).
1647 : !! ph_intmeth=Integration method for phonons: 1 for gaussian, 2 for tetrahedra
1648 : !! smear=Gaussian broadening used to approximate the Dirac delta.
1649 : !! ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift).
1650 : !! nqshift=Number of shifts used to generated the Q-mesh.
1651 : !! qshift(3,nqshift)=The shifts.
1652 : !! comm=MPI communicator
1653 : !! [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop are ignored and
1654 : !! A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation.
1655 : !! [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account
1656 : !! the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1
1657 : !! Other values of qptopt can be used for debugging purpose.
1658 : !!
1659 : !! OUTPUT
1660 : !! a2f<a2fw_t>=Structure storing the Eliashberg function a2F(w).
1661 : !!
1662 : !! SOURCE
1663 :
1664 10 : subroutine a2fw_init(a2f, gams, cryst, ifc, ph_intmeth, wstep, wminmax, smear, ngqpt, nqshift, qshift, comm, &
1665 : qintp, qptopt) ! optional
1666 :
1667 : !Arguments ------------------------------------
1668 : !scalars
1669 : class(a2fw_t),target,intent(out) :: a2f
1670 : integer,intent(in) :: ph_intmeth,nqshift,comm
1671 : integer,intent(in),optional :: qptopt
1672 : real(dp),intent(in) :: wstep,smear
1673 : logical,optional,intent(in) :: qintp
1674 : type(phgamma_t),intent(inout) :: gams
1675 : type(ifc_type),intent(in) :: ifc
1676 : type(crystal_t),intent(in) :: cryst
1677 : !arrays
1678 : integer,intent(in) :: ngqpt(3)
1679 : real(dp),intent(in) :: wminmax(2),qshift(3,nqshift)
1680 :
1681 : !Local variables -------------------------
1682 : !scalars
1683 : integer,parameter :: master = 0
1684 : integer :: my_qptopt,iq_ibz,nqibz,ount,my_rank,nproc,cnt
1685 : integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz
1686 : integer :: iene, jene, itemp, ntemp, jene_jump, spinfactor
1687 : real(dp) :: cpu,wall,gflops
1688 : real(dp) :: lambda_iso,omega,omega_log,xx,omega_min,omega_max,ww,mustar,tc_macmill
1689 : real(dp) :: temp_el, min_temp, delta_temp, chempot, ene1, ene2, G0
1690 : logical :: do_qintp
1691 : character(len=500) :: msg
1692 10 : type(htetra_t) :: tetra_q
1693 : !arrays
1694 : integer :: qptrlatt(3,3),new_qptrlatt(3,3)
1695 20 : real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom)
1696 40 : real(dp) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3), invphfrq(gams%natom3)
1697 10 : real(dp),allocatable :: my_qshift(:,:), gamma_ph_ee(:,:,:,:), tmp_a2f(:)
1698 10 : real(dp),contiguous, pointer :: a2f_1d(:)
1699 10 : real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:), a2f_1mom(:),a2flogmom(:),a2flogmom_int(:),wdt(:,:)
1700 10 : real(dp),allocatable :: lambda_tetra(:,:,:),phfreq_tetra(:,:), tmp_gaussian(:,:)
1701 10 : real(dp),allocatable :: a2feew_partial(:), a2feew_partial_int(:), a2feew_w(:), a2feew_w_int(:)
1702 : ! *********************************************************************
1703 :
1704 10 : call cwtime(cpu, wall, gflops, "start")
1705 10 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
1706 :
1707 10 : my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt
1708 10 : do_qintp = .True.; if (present(qintp)) do_qintp = qintp
1709 10 : nsppol = gams%nsppol; natom3 = gams%natom3
1710 :
1711 10 : if (gams%prteliash == 3) then
1712 0 : ABI_MALLOC(gamma_ph_ee, (gams%nene, gams%nene, gams%natom3, gams%nsppol))
1713 : end if
1714 :
1715 10 : if (do_qintp) then
1716 : ! Generate fine q-mesh, find the IBZ and the corresponding weights.
1717 5 : qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
1718 :
1719 : call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, &
1720 5 : new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift)
1721 5 : ABI_FREE(qbz)
1722 :
1723 : ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
1724 20 : a2f%ngqpt = ngqpt; a2f%nqshift = size(my_qshift, dim=2)
1725 15 : ABI_MALLOC(a2f%qshift, (3, a2f%nqshift))
1726 30 : a2f%qshift = my_qshift
1727 5 : ABI_FREE(my_qshift)
1728 :
1729 : else
1730 : ! No interpolation. Use q-mesh parameters from gams%
1731 20 : a2f%ngqpt = gams%ngqpt; a2f%nqshift = 1
1732 5 : nqibz = gams%nqibz
1733 15 : ABI_MALLOC(qibz, (3, nqibz))
1734 15 : ABI_MALLOC(wtq, (nqibz))
1735 5 : ABI_MALLOC(a2f%qshift, (3, a2f%nqshift))
1736 320 : qibz = gams%qibz; wtq = gams%wtq
1737 25 : a2f%qshift = zero ! Note: assuming q-mesh centered on Gamma.
1738 : end if
1739 :
1740 10 : call cwtime_report(" a2fw_init, q-setup:", cpu, wall, gflops)
1741 :
1742 : ! Define Min and max frequency for the mesh (enlarge it a bit)
1743 10 : omega_min = wminmax(1); omega_max = wminmax(2)
1744 10 : omega_min = omega_min - 0.1*abs(omega_min)
1745 10 : if (omega_min >= zero) omega_min = one/Ha_meV
1746 10 : omega_max = omega_max + 0.1*abs(omega_max)
1747 :
1748 10 : a2f%nsppol = nsppol; a2f%natom3 = gams%natom3; a2f%smear = smear
1749 10 : a2f%omega_min = omega_min; a2f%omega_max = omega_max
1750 10 : nomega = int((omega_max - omega_min) / wstep); a2f%nomega = nomega; a2f%wstep = wstep
1751 10 : a2f%nene = gams%nene
1752 10 : a2f%enemin = gams%enemin
1753 10 : a2f%deltaene = gams%deltaene
1754 :
1755 30 : ABI_MALLOC(a2f%n0, (nsppol))
1756 30 : a2f%n0 = gams%n0
1757 : ! Build linear mesh.
1758 30 : ABI_MALLOC(a2f%omega, (nomega))
1759 4750 : a2f%omega = arth(omega_min, wstep, nomega)
1760 :
1761 19020 : ABI_CALLOC(a2f%vals, (nomega,0:natom3, nsppol))
1762 19010 : ABI_CALLOC(a2f%lambdaw, (nomega,0:natom3, nsppol))
1763 :
1764 10 : if (gams%prteliash == 3) then
1765 0 : ABI_CALLOC(a2f%vals_ee, (gams%nene, gams%nene, nomega, nsppol))
1766 : end if
1767 :
1768 20 : ABI_MALLOC(tmp_a2f, (nomega))
1769 :
1770 10 : if (ph_intmeth == 2) then
1771 : ! Prepare tetrahedron integration.
1772 8 : call cwtime(cpu, wall, gflops, "start")
1773 :
1774 8 : qptrlatt = 0; qptrlatt(1, 1) = a2f%ngqpt(1); qptrlatt(2, 2) = a2f%ngqpt(2); qptrlatt(3, 3) = a2f%ngqpt(3)
1775 8 : tetra_q = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f%nqshift, a2f%qshift, nqibz, qibz, comm, msg, ierr)
1776 8 : if (ierr /= 0) ABI_ERROR(msg)
1777 :
1778 1971 : ABI_CALLOC(lambda_tetra, (nqibz, natom3, nsppol))
1779 1955 : ABI_CALLOC(phfreq_tetra, (nqibz, natom3))
1780 8 : cnt = 0
1781 641 : do iq_ibz = 1, nqibz
1782 633 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
1783 633 : call ifc%fourq(cryst, qibz(:,iq_ibz), phfrq, displ_cart)
1784 : ! save for tetrahedron interpolation
1785 2540 : phfreq_tetra(iq_ibz,:) = phfrq(:)
1786 : end do
1787 8 : call xmpi_sum(phfreq_tetra, comm, ierr)
1788 16 : call cwtime_report(" a2fw_init%tetra", cpu, wall, gflops)
1789 : end if
1790 :
1791 10 : call cwtime(cpu, wall, gflops, "start")
1792 :
1793 : ! DEV_MJV
1794 : !open(unit=900, file="a2fvals_ee.dat"); write (900,*) '# do_qintp ', do_qintp
1795 :
1796 : ! Loop over spins and q-points in the IBZ. For the moment parallelize over iq_ibz
1797 20 : do spin=1,nsppol
1798 10 : cnt = 0
1799 806 : do iq_ibz=1,nqibz
1800 : ! TODO: for the moment the memory is not distributed, only the calculation exception is vals_ee
1801 786 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
1802 :
1803 : ! Interpolate or evaluate gamma directly.
1804 786 : if (do_qintp) then
1805 725 : if (gams%prteliash == 3) then
1806 0 : call gams%interp(cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee=gamma_ph_ee(:,:,:,spin))
1807 : else
1808 725 : call gams%interp(cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart)
1809 : end if
1810 : else
1811 61 : if (gams%prteliash == 3) then
1812 0 : call gams%eval_qibz(cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee=gamma_ph_ee(:,:,:,spin))
1813 : else
1814 61 : call gams%eval_qibz(cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph, displ_cart)
1815 : end if
1816 : end if
1817 :
1818 10 : select case (ph_intmeth)
1819 : case (1)
1820 : ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f
1821 : ! (note that unstable modes are included).
1822 612 : ABI_MALLOC(tmp_gaussian, (nomega, natom3))
1823 612 : do mu=1,natom3
1824 217566 : tmp_a2f = zero
1825 217566 : do iw=1,nomega
1826 217107 : xx = a2f%omega(iw) - phfrq(mu)
1827 217107 : tmp_gaussian(iw, mu) = gaussian(xx, smear)
1828 217566 : tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * lambda_ph(mu) * abs(phfrq(mu))
1829 : end do
1830 217719 : a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + tmp_a2f * wtq(iq_ibz)
1831 : end do
1832 :
1833 153 : if (gams%prteliash == 3) then
1834 : ! reset phfrq for low freq modes
1835 0 : invphfrq = zero
1836 0 : do mu=1,natom3
1837 0 : if (abs(phfrq(mu)) > EPHTK_WTOL) invphfrq(mu) = one / abs(phfrq(mu))
1838 : end do
1839 :
1840 0 : do iene= 1, gams%nene
1841 0 : do jene= 1, gams%nene
1842 0 : tmp_a2f = zero
1843 : ! TODO: following block is just a GEMM
1844 0 : do mu=1,natom3
1845 0 : do iw=1,nomega
1846 0 : tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * gamma_ph_ee(jene,iene,mu,spin) * invphfrq(mu)
1847 : end do
1848 : end do
1849 :
1850 0 : a2f%vals_ee(jene, iene, :, spin) = a2f%vals_ee(jene, iene, :, spin) + tmp_a2f(:) * wtq(iq_ibz)
1851 : !if (iene == gams%nene/2 .and. jene == gams%nene/2) then
1852 : ! write (900, '(a,E20.10,2x,2x,I6,3E20.10)') '#', wtq(iq_ibz), iq_ibz, invphfrq(1:3)
1853 : ! do iw=1,nomega
1854 : ! write (900, '(i6,2x,E20.10,2x,3E20.10,2x,3E20.10)') iw, a2f%vals_ee(jene, iene, iw, spin), &
1855 : ! gamma_ph_ee(jene,iene,:,spin), tmp_gaussian(iw,1:3)
1856 : ! end do
1857 : ! write (900,*)
1858 : !end if
1859 : end do
1860 : end do
1861 : end if
1862 :
1863 153 : ABI_FREE(tmp_gaussian)
1864 :
1865 : case (2)
1866 : ! Tetra: store data.
1867 2532 : do mu=1,natom3
1868 2532 : lambda_tetra(iq_ibz, mu, spin) = lambda_ph(mu) * abs(phfrq(mu))
1869 : end do
1870 :
1871 633 : if (gams%prteliash == 3) then
1872 0 : ABI_ERROR("Eliashberg function with tetra not coded")
1873 : end if
1874 :
1875 : case default
1876 786 : ABI_ERROR(sjoin("Wrong ph_intmeth:", itoa(ph_intmeth)))
1877 : end select
1878 :
1879 : end do ! iq_ibz
1880 : end do ! spin
1881 :
1882 : !DEV_MJV
1883 : !close(900)
1884 :
1885 10 : if (ph_intmeth == 2) then
1886 : ! workspace for tetra.
1887 24 : ABI_MALLOC(wdt, (nomega, 2))
1888 :
1889 : ! For each mode get its contribution
1890 16 : do spin=1,nsppol
1891 40 : do mu=1,natom3
1892 24 : cnt = 0
1893 1931 : do iq_ibz=1,nqibz
1894 : ! NB: if we are interpolating the gamma, nqibz > gams%nqibz
1895 1899 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
1896 :
1897 1899 : call tetra_q%get_onewk(iq_ibz, gams%bcorr, nomega, nqibz, phfreq_tetra(:,mu), omega_min, omega_max, one, wdt)
1898 1802151 : wdt = wdt * wtq(iq_ibz)
1899 :
1900 : ! Accumulate (Integral of a2F is computed afterwards)
1901 900150 : a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + wdt(:,1) * lambda_tetra(iq_ibz,mu,spin)
1902 : !a2f%lambdaw(:,mu,spin) = a2f%lambdaw(:,mu,spin) + wdt(:,2) * lambda_tetra(iq_ibz, mu, spin)
1903 : end do
1904 : end do
1905 : end do
1906 :
1907 : ! Free memory allocated for tetra.
1908 8 : ABI_FREE(wdt)
1909 8 : ABI_FREE(lambda_tetra)
1910 8 : ABI_FREE(phfreq_tetra)
1911 8 : call tetra_q%free()
1912 : end if
1913 :
1914 : ! Collect final results on each node
1915 10 : call xmpi_sum(a2f%vals, comm, ierr)
1916 20 : do spin=1,nsppol
1917 18940 : a2f%vals(:,0,spin) = sum(a2f%vals(:,1:natom3,spin), dim=2)
1918 : ! previously would divide by g(eF, spin)
1919 : !a2f%vals(:,:,spin) = a2f%vals(:,:,spin) / (two_pi*a2f%n0(spin))
1920 : end do
1921 :
1922 10 : if (gams%prteliash == 3) then
1923 : ! For the moment vals_ee only works with gaussians
1924 0 : call xmpi_sum(a2f%vals_ee, comm, ierr)
1925 0 : do spin=1,nsppol
1926 0 : a2f%vals_ee(:,:,:,spin) = a2f%vals_ee(:,:,:,spin) / (two * pi * gams%n0(spin))
1927 : end do
1928 : end if
1929 :
1930 10 : call cwtime_report(" a2fw_init, a2f_eval:", cpu, wall, gflops)
1931 :
1932 20 : ABI_MALLOC(a2f_1mom, (nomega))
1933 20 : ABI_MALLOC(a2flogmom, (nomega))
1934 20 : ABI_MALLOC(a2flogmom_int, (nomega))
1935 :
1936 : !call a2fw_print_info()
1937 :
1938 : ! Compute lambda(w) resolved in spin and phonon mode.
1939 : ! lambda(w) = int a2F(w')/w' dw'
1940 20 : do spin=1,nsppol
1941 60 : do mu=0,natom3
1942 :
1943 18960 : do iw=1,nomega
1944 18920 : ww = a2f%omega(iw)
1945 18960 : if (abs(ww) > EPHTK_WTOL) then
1946 18920 : a2f_1mom(iw) = a2f%vals(iw, mu, spin) / abs(ww)
1947 : else
1948 0 : a2f_1mom(iw) = zero
1949 : end if
1950 : end do
1951 :
1952 50 : call simpson_int(nomega, wstep, a2f_1mom, a2f%lambdaw(:, mu, spin))
1953 : end do
1954 : end do
1955 :
1956 : ! print log moments of the alpha^2 F function (for each spin, note degeneracy factor).
1957 10 : spinfactor = int(2 / nsppol / gams%nspinor)
1958 20 : do spin=1,nsppol
1959 10 : a2f_1d => a2f%vals(:,0,spin)
1960 :
1961 10 : lambda_iso = a2f%get_moment(0, spin)
1962 :
1963 : ! Get log moment of alpha^2F.
1964 4740 : a2flogmom = zero
1965 4740 : do iw=1,nomega
1966 4730 : omega = a2f%omega(iw)
1967 4740 : if (abs(omega) > EPHTK_WTOL) then
1968 : !a2flogmom(iw) = (two / lambda_iso) * a2f_1d(iw) * log(abs(omega)) / abs(omega)
1969 : ! I (WHO??? GMATTEO?) think this is the correct expression.
1970 4730 : a2flogmom(iw) = (spinfactor / lambda_iso) * a2f_1d(iw) * log(abs(omega)) / abs(omega)
1971 : end if
1972 : end do
1973 10 : call simpson_int(nomega, wstep, a2flogmom, a2flogmom_int)
1974 10 : omega_log = exp(a2flogmom_int(nomega))
1975 :
1976 10 : mustar = 0.12
1977 10 : tc_macmill = omega_log/1.2_dp * exp((-1.04_dp*(one+lambda_iso)) / (lambda_iso-mustar*(one+0.62_dp*lambda_iso)))
1978 :
1979 20 : if (my_rank == master) then
1980 10 : ount = std_out
1981 10 : if (nsppol > 1) then
1982 0 : write(msg,'(3a)') ch10,' Warning: some of the following quantities should be integrated over spin', ch10
1983 0 : call wrtout(ount, msg)
1984 : end if
1985 :
1986 10 : if (do_qintp) then
1987 5 : write(ount,'(a,i3)') ' Superconductivity: isotropic evaluation of parameters from electron-phonon coupling (interpolated) for spin ', spin
1988 : else
1989 5 : write(ount,'(a,i3)') ' Superconductivity: isotropic evaluation of parameters from electron-phonon coupling (coarse grid) for spin ', spin
1990 : endif
1991 10 : write(ount,'(a,es16.6)') ' isotropic lambda = ',lambda_iso
1992 10 : write(ount,'(a,es16.6,a,es16.6,a)' )' omegalog = ',omega_log,' (Ha) ', omega_log * Ha_K, ' (Kelvin) '
1993 10 : write(ount,'(a,es16.6,a,es16.6,a)') ' MacMillan Tc = ',tc_macmill,' (Ha) ', tc_macmill * Ha_K, ' (Kelvin) '
1994 10 : write(ount,"(a)") ' Positive moments of alpha2F:'
1995 10 : write(ount,'(a,es16.6)' ) ' lambda <omega^2> = ',a2f%get_moment(2, spin)
1996 10 : write(ount,'(a,es16.6)' ) ' lambda <omega^3> = ',a2f%get_moment(3, spin)
1997 10 : write(ount,'(a,es16.6)' ) ' lambda <omega^4> = ',a2f%get_moment(4, spin)
1998 10 : write(ount,'(a,es16.6)' ) ' lambda <omega^5> = ',a2f%get_moment(5, spin)
1999 : end if
2000 : end do
2001 :
2002 : ! Calculate the temperature dependence of the a2f(e,e',w) integrals (G_0(T_e)
2003 : ! as in PRL 110 016405 (2013) [[cite:Arnaud2013]])
2004 10 : if (gams%prteliash == 3 .and. my_rank == master) then
2005 0 : ntemp = 100
2006 0 : min_temp = zero
2007 0 : delta_temp = 40._dp ! Kelvin
2008 0 : ABI_MALLOC(a2feew_partial, (a2f%nene))
2009 0 : ABI_MALLOC(a2feew_partial_int, (a2f%nene))
2010 0 : ABI_MALLOC(a2feew_w, (nomega))
2011 0 : ABI_MALLOC(a2feew_w_int, (nomega))
2012 :
2013 0 : if (open_file("EPC_strength_aafo_T.dat", msg, newunit=ount, form="formatted", action="write", status="unknown") /= 0) then
2014 0 : ABI_ERROR(msg)
2015 : end if
2016 :
2017 0 : write(ount, "(a)")"# temp_el, G_0(T_e) in W/m^3/K, spin"
2018 0 : do spin=1,nsppol
2019 0 : do itemp = 1, ntemp
2020 0 : temp_el = min_temp + (itemp-1)*delta_temp
2021 :
2022 : ! TODO: need to evolve the chemical potential with T, but I do not have access to the full DOS here!!
2023 : ! possible fix using local information on DOS variation near E_F...
2024 0 : chempot = a2f%enemin + half * a2f%nene * a2f%deltaene
2025 :
2026 0 : do iw=1,nomega
2027 0 : omega = a2f%omega(iw)
2028 0 : jene_jump = nint(omega/a2f%deltaene)
2029 0 : a2feew_partial = zero
2030 0 : do iene=1,a2f%nene
2031 0 : ene1 = a2f%enemin + (iene-1)*a2f%deltaene
2032 0 : ene2 = ene1 + omega
2033 : a2feew_partial(iene) = a2f%vals_ee(min(a2f%nene,iene+jene_jump), iene, iw, spin) * &
2034 0 : (fermi_dirac(ene1, chempot, temp_el/Ha_K) - fermi_dirac(ene2, chempot, temp_el/Ha_K))
2035 : end do
2036 0 : call simpson_int(a2f%nene, a2f%deltaene, a2feew_partial, a2feew_partial_int)
2037 0 : a2feew_w(iw) = a2feew_partial_int(a2f%nene)
2038 : end do
2039 0 : call simpson_int(nomega, wstep, a2feew_w, a2feew_w_int)
2040 0 : G0 = a2feew_w_int(nomega) * two_pi * a2f%n0(spin) / cryst%ucvol
2041 : ! conversion factor for G0 to SI units = Ha_J / Time_Sec / (Bohr_meter)**3 ~ 1.2163049915755545e+30
2042 0 : write(ount, "(2(e20.10,2x),i5)") temp_el, G0 * kb_HaK / Time_Sec / (Bohr_meter)**3, spin !* Ha_J???
2043 : end do
2044 : end do
2045 0 : close(ount)
2046 :
2047 0 : ABI_FREE(a2feew_partial)
2048 0 : ABI_FREE(a2feew_partial_int)
2049 0 : ABI_FREE(a2feew_w)
2050 0 : ABI_FREE(a2feew_w_int)
2051 0 : ABI_FREE(gamma_ph_ee)
2052 : end if
2053 :
2054 10 : ABI_FREE(tmp_a2f)
2055 10 : ABI_FREE(a2f_1mom)
2056 10 : ABI_FREE(a2flogmom)
2057 10 : ABI_FREE(a2flogmom_int)
2058 10 : ABI_FREE(qibz)
2059 10 : ABI_FREE(wtq)
2060 :
2061 70 : end subroutine a2fw_init
2062 : !!***
2063 :
2064 : !----------------------------------------------------------------------
2065 :
2066 : !!****f* m_phgamma/a2fw_get_moment
2067 : !! NAME
2068 : !! a2fw_get_moment
2069 : !!
2070 : !! FUNCTION
2071 : !! Compute \int dw [a2F(w)/w] w^n. From Allen PRL 59 1460 [[cite:Allen1987]].
2072 : !! See also [[cite:Grimvall1981]], Eq 6.72 page 175)
2073 : !!
2074 : !! INPUTS
2075 : !! nn=Value of n
2076 : !! spin=The spin component. 0 to sum over spins.
2077 : !!
2078 : !! OUTPUT
2079 : !! a2fw_get_moment = \int dw [a2F(w)/w] w^n
2080 : !! [out_int(x)] = \int^{x} dw [a2F(w)/w] w^n
2081 : !!
2082 : !! SOURCE
2083 :
2084 50 : real(dp) function a2fw_get_moment(a2f, nn, spin, out_int)
2085 :
2086 : !Arguments ------------------------------------
2087 : class(a2fw_t),intent(in) :: a2f
2088 : integer,intent(in) :: spin,nn
2089 : !arrays
2090 : real(dp),intent(out),optional :: out_int(a2f%nomega)
2091 :
2092 : !Local variables -------------------------
2093 : !scalars
2094 : integer :: iw
2095 : real(dp) :: omg, omg_nm1
2096 : !arrays
2097 50 : real(dp) :: ff(a2f%nomega), int_ff(a2f%nomega), values(a2f%nomega)
2098 : ! *********************************************************************
2099 :
2100 : ! Construct the integrand function. [a2F(w)/w] w^n
2101 47350 : ff = zero; int_ff = zero
2102 :
2103 50 : if (spin == 0) then
2104 0 : if (a2f%nsppol == 1) then
2105 0 : values = a2f%vals(:, 0, 1)
2106 : !if (a2f%nspinor == 1) values = two * values
2107 0 : else if (a2f%nsppol == 2) then
2108 0 : values = sum(a2f%vals(:, 0, :), dim=2)
2109 : end if
2110 :
2111 : else
2112 23700 : values = a2f%vals(:, 0, spin)
2113 : end if
2114 :
2115 50 : if (nn - 1 >= 0) then
2116 18960 : do iw=1,a2f%nomega
2117 18920 : omg = a2f%omega(iw)
2118 18920 : omg_nm1 = omg ** (nn - 1)
2119 18960 : ff(iw) = values(iw) * omg_nm1
2120 : end do
2121 : else
2122 4740 : do iw=1,a2f%nomega
2123 4730 : omg = a2f%omega(iw)
2124 4730 : omg_nm1 = zero; if (abs(omg) > EPHTK_WTOL) omg_nm1 = omg**(nn-1)
2125 4740 : ff(iw) = values(iw) * omg_nm1
2126 : end do
2127 : end if
2128 :
2129 : ! Integration with simpson rule on a linear mesh.
2130 50 : call simpson_int(a2f%nomega, a2f%wstep, ff, int_ff)
2131 :
2132 50 : a2fw_get_moment = int_ff(a2f%nomega)
2133 50 : if (present(out_int)) out_int = int_ff
2134 :
2135 50 : end function a2fw_get_moment
2136 : !!***
2137 :
2138 : !----------------------------------------------------------------------
2139 :
2140 : !!****f* m_phgamma/a2fw_tr_moment
2141 : !! NAME
2142 : !! a2fw_tr_moment
2143 : !!
2144 : !! FUNCTION
2145 : !! Compute \int dw [a2F_tr(w)/w] w^n
2146 : !! From Allen PRL 59 1460 [[cite:Allen1987]] and later PRB papers.
2147 : !! See also [[cite:Grimvall1981]] book.
2148 : !!
2149 : !! INPUTS
2150 : !! nn=Value of n
2151 : !! spin=The spin component
2152 : !!
2153 : !! OUTPUT
2154 : !! a2fw_tr_moment = \int dw [a2F_tr(w)/w] w^n
2155 : !! [out_int(x)] = \int^{x} dw [a2F_tr(w)/w] w^n
2156 : !!
2157 : !! SOURCE
2158 :
2159 50 : function a2fw_tr_moment(a2f_tr, nn, spin, out_int)
2160 :
2161 : !Arguments ------------------------------------
2162 : !scalars
2163 : class(a2fw_tr_t),intent(in) :: a2f_tr
2164 : integer,intent(in) :: spin,nn
2165 : !arrays
2166 : real(dp),intent(out),optional :: out_int(a2f_tr%nomega,3,3)
2167 : real(dp) :: a2fw_tr_moment(3,3)
2168 :
2169 : !Local variables -------------------------
2170 : !scalars
2171 : integer :: iw, idir, jdir
2172 : real(dp) :: omg,omg_nm1
2173 : !arrays
2174 50 : real(dp) :: ff(a2f_tr%nomega),int_ff(a2f_tr%nomega)
2175 : ! *********************************************************************
2176 :
2177 200 : do jdir = 1, 3
2178 650 : do idir = 1, 3
2179 : ! Construct the integrand function. [a2F_tr(w)/w] w^n
2180 426150 : ff = zero; int_ff = zero
2181 :
2182 450 : if (nn-1 >= 0) then
2183 170640 : do iw=1,a2f_tr%nomega
2184 170280 : omg = a2f_tr%omega(iw)
2185 170280 : omg_nm1 = omg ** (nn-1)
2186 170640 : ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1
2187 : end do
2188 : else
2189 42660 : do iw=1,a2f_tr%nomega
2190 42570 : omg = a2f_tr%omega(iw)
2191 42570 : omg_nm1 = zero; if (abs(omg) > EPHTK_WTOL) omg_nm1 = omg**(nn-1)
2192 42660 : ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1
2193 : end do
2194 : end if
2195 :
2196 : ! Integration with simpson rule on a linear mesh.
2197 450 : call simpson_int(a2f_tr%nomega, a2f_tr%wstep, ff, int_ff)
2198 :
2199 450 : a2fw_tr_moment(idir, jdir) = int_ff(a2f_tr%nomega)
2200 600 : if (present(out_int)) out_int(:,idir,jdir) = int_ff(:)
2201 : end do
2202 : end do
2203 :
2204 : end function a2fw_tr_moment
2205 : !!***
2206 :
2207 : !----------------------------------------------------------------------
2208 :
2209 : !!****f* m_phgamma/a2fw_write
2210 : !! NAME
2211 : !! a2fw_write
2212 : !!
2213 : !! FUNCTION
2214 : !! Write alpha^2F(w) to an external file in text form
2215 : !!
2216 : !! INPUTS
2217 : !! basename=Filename for output.
2218 : !! post=String appended to netcdf variables e.g. _qcoarse, _qintp
2219 : !! ncid=Netcdf file handler. Set it to nctk_noid to disable output.
2220 : !!
2221 : !! OUTPUT
2222 : !! Output is written to file. This routine should be called by one MPI proc.
2223 : !!
2224 : !! SOURCE
2225 :
2226 10 : subroutine a2fw_write(a2f, basename, post, ncid)
2227 :
2228 : !Arguments ------------------------------------
2229 : !scalars
2230 : class(a2fw_t),intent(in) :: a2f
2231 : character(len=*),intent(in) :: basename, post
2232 : integer,intent(in) :: ncid
2233 :
2234 : !Local variables -------------------------
2235 : !scalars
2236 : integer :: iw,spin,unt,ii,mu, ncerr
2237 : character(len=500) :: dim1_name, msg
2238 : character(len=fnlen) :: path
2239 : ! *********************************************************************
2240 :
2241 : ! Write spin-resolved a2F(w)
2242 10 : path = strcat(basename, "_A2FW")
2243 10 : if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then
2244 0 : ABI_ERROR(msg)
2245 : end if
2246 :
2247 10 : call write_a2fw_header()
2248 :
2249 10 : if (a2f%nsppol == 1) then
2250 10 : write(unt,'(a)')"# Frequency, a2F(w), lambda(w)"
2251 4740 : do iw=1,a2f%nomega
2252 4740 : write(unt,*) a2f%omega(iw), a2f%vals(iw,0,1), a2f%lambdaw(iw,0,1)
2253 : end do
2254 :
2255 : else
2256 0 : write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..."
2257 0 : do iw=1,a2f%nomega
2258 0 : write(unt,'(E20.10,3x,3(2E20.10, 2x))') a2f%omega(iw), &
2259 0 : sum(a2f%vals(iw,0,:)) , sum(a2f%lambdaw(iw,0,:)), & ! TOT
2260 0 : a2f%vals(iw,0,1) , a2f%lambdaw(iw,0,1), & ! UP
2261 0 : a2f%vals(iw,0,2) , a2f%lambdaw(iw,0,2) ! DOWN
2262 : end do
2263 : end if
2264 :
2265 10 : close(unt)
2266 :
2267 : ! Write phonon contributions to a2F(w)
2268 10 : path = strcat(basename, "_PH_A2FW")
2269 10 : if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then
2270 0 : ABI_ERROR(msg)
2271 : end if
2272 10 : call write_a2fw_header()
2273 :
2274 10 : if (a2f%nsppol == 1) then
2275 50 : do mu=0,a2f%natom3
2276 40 : write(unt,'(a,i0)')"# Phonon mode ",mu
2277 40 : write(unt,'(a)')"# Frequency, a2F(w), lambda(w)"
2278 18960 : do iw=1,a2f%nomega
2279 18960 : write(unt,*) a2f%omega(iw), a2f%vals(iw,mu,1), a2f%lambdaw(iw,mu,1)
2280 : end do
2281 130 : do ii=1,2; write(unt,'(a)')""; end do
2282 : end do
2283 :
2284 : else
2285 0 : do mu=0,a2f%natom3
2286 0 : write(unt,'(a,i0)')"# Phonon mode ",mu
2287 0 : write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..."
2288 0 : do iw=1,a2f%nomega
2289 0 : write(unt,'(E20.10,3x,3(2E20.10, 2x))') a2f%omega(iw), &
2290 0 : sum(a2f%vals(iw,mu,:)), sum(a2f%lambdaw(iw,mu,:)), & ! TOT
2291 0 : a2f%vals(iw,mu,1) , a2f%lambdaw(iw,mu,1), & ! UP
2292 0 : a2f%vals(iw,mu,2) , a2f%lambdaw(iw,mu,2) ! DOWN
2293 : end do
2294 : end do
2295 : end if
2296 :
2297 10 : close(unt)
2298 :
2299 10 : if (ncid /= nctk_noid) then
2300 : ! Define dimensions.
2301 10 : dim1_name = strcat("a2f_nomega", post)
2302 : ncerr = nctk_def_dims(ncid, [ &
2303 : nctkdim_t(dim1_name, a2f%nomega), &
2304 : nctkdim_t("natom3p1", a2f%natom3 + 1), nctkdim_t("number_of_spins", a2f%nsppol) &
2305 40 : ], defmode=.True.)
2306 10 : NCF_CHECK(ncerr)
2307 :
2308 : ncerr = nctk_def_arrays(ncid, [ &
2309 : nctkarr_t(strcat('a2f_mesh', post), "dp", dim1_name), &
2310 : nctkarr_t(strcat('a2f_values', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")), &
2311 : nctkarr_t(strcat('a2f_lambdaw', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")) &
2312 40 : ])
2313 10 : NCF_CHECK(ncerr)
2314 :
2315 10 : NCF_CHECK(nctk_set_datamode(ncid))
2316 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_mesh", post)), a2f%omega))
2317 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_values", post)), a2f%vals))
2318 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_lambdaw", post)), a2f%lambdaw))
2319 : end if
2320 :
2321 : contains
2322 :
2323 20 : subroutine write_a2fw_header()
2324 :
2325 : ! Output the header.
2326 20 : write(unt,'(a)') '#'
2327 20 : write(unt,'(a)') '# ABINIT package: a2F(w) file'
2328 20 : write(unt,'(a)') '#'
2329 20 : write(unt,'(a)') '# a2F(w) function integrated over the FS. omega in a.u.'
2330 20 : write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f%ngqpt))
2331 20 : write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min,&
2332 40 : ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep
2333 20 : write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f%smear
2334 40 : write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0)
2335 40 : do spin=1,a2f%nsppol
2336 40 : write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin)
2337 : end do
2338 60 : do ii=1,2; write(unt,'(a)') "# "; end do
2339 :
2340 20 : end subroutine write_a2fw_header
2341 :
2342 : end subroutine a2fw_write
2343 : !!***
2344 :
2345 : !----------------------------------------------------------------------
2346 : !!****f* m_phgamma/a2fw_ee_write
2347 : !! NAME
2348 : !! a2fw_ee_write
2349 : !!
2350 : !! FUNCTION
2351 : !! Write alpha^2F(e,e',w) to an external file in text form
2352 : !!
2353 : !! INPUTS
2354 : !! basename=Filename for output.
2355 : !!
2356 : !! OUTPUT
2357 : !! Output is written to file. This routine should be called by one MPI proc.
2358 : !!
2359 : !! SOURCE
2360 :
2361 0 : subroutine a2fw_ee_write(a2f, basename)
2362 :
2363 : !Arguments ------------------------------------
2364 : class(a2fw_t),intent(in) :: a2f
2365 : character(len=*),intent(in) :: basename
2366 :
2367 : !Local variables -------------------------
2368 : !scalars
2369 : integer :: iw, spin, unt, ii, iene, jene
2370 : real(dp) :: ene1, ene2
2371 : character(len=500) :: msg
2372 : character(len=fnlen) :: path
2373 : ! *********************************************************************
2374 :
2375 : ! Write spin-resolved a2F(e,e',w)
2376 0 : path = strcat(basename, "_A2FEEW")
2377 0 : if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then
2378 0 : ABI_ERROR(msg)
2379 : end if
2380 :
2381 0 : call write_a2fw_header()
2382 :
2383 0 : write(unt,'(a)')"# en2, en1, Frequency, a2F_tot(w) (presently summed over spin)"
2384 0 : do iw=1,a2f%nomega
2385 0 : do iene=1,a2f%nene
2386 0 : ene1 = a2f%enemin + (iene-1)*a2f%deltaene
2387 0 : do jene=1,a2f%nene
2388 0 : ene2 = a2f%enemin + (jene-1)*a2f%deltaene
2389 0 : write(unt,'(3E20.10,2x,E20.10)') ene2, ene1, a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:)) ! TOT
2390 : end do
2391 : end do
2392 : end do
2393 :
2394 0 : close(unt)
2395 :
2396 : ! Write spin-resolved a2F(ef, ef, w)
2397 0 : path = strcat(basename, "_A2FW_reference")
2398 0 : if (open_file(path, msg, newunit=unt, form="formatted", status="unknown") /= 0) then
2399 0 : ABI_ERROR(msg)
2400 : end if
2401 :
2402 0 : call write_a2fw_header()
2403 :
2404 0 : write(unt,'(a)')"# Frequency, a2F_tot(ef,ef,w) for comparison with normal a2F(w) (presently summed over spin)"
2405 0 : iene = int(a2f%nene/2)
2406 0 : jene = int(a2f%nene/2)
2407 0 : do iw=1,a2f%nomega
2408 0 : write(unt,'(E20.10,2x,E20.10)') a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:)) ! TOT
2409 : end do
2410 :
2411 0 : close(unt)
2412 :
2413 : contains
2414 :
2415 0 : subroutine write_a2fw_header()
2416 :
2417 : ! Output the header.
2418 0 : write(unt,'(a)') '#'
2419 0 : write(unt,'(a)') '# ABINIT package: a2F(e,eprime,w) file'
2420 0 : write(unt,'(a)') '#'
2421 0 : write(unt,'(a)') '# a2F(e,eprime,w) function integrated over the FS. omega in a.u.'
2422 0 : write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f%ngqpt))
2423 0 : write(unt,'(a,i0,2(a,e16.6))')'# number of energies: ',a2f%nene," from enemin: ",a2f%enemin,' Ha with step ', a2f%deltaene
2424 0 : write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min, &
2425 0 : ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep
2426 0 : write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f%smear
2427 0 : write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0)
2428 0 : do spin=1,a2f%nsppol
2429 0 : write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin)
2430 : end do
2431 0 : do ii=1,2
2432 0 : write(unt,'(a)') "# "
2433 : end do
2434 :
2435 0 : end subroutine write_a2fw_header
2436 :
2437 : end subroutine a2fw_ee_write
2438 : !!***
2439 :
2440 : !----------------------------------------------------------------------
2441 :
2442 : !!****f* m_phgamma/a2fw_tr_free
2443 : !! NAME
2444 : !! a2fw_tr_free
2445 : !!
2446 : !! FUNCTION
2447 : !! Free the memory allocated in a2fw_tr_t
2448 : !!
2449 : !! SOURCE
2450 :
2451 10 : subroutine a2fw_tr_free(a2f_tr)
2452 :
2453 : !Arguments ------------------------------------
2454 : class(a2fw_tr_t),intent(inout) :: a2f_tr
2455 : ! *********************************************************************
2456 :
2457 : ! integer
2458 10 : ABI_SFREE(a2f_tr%qshift)
2459 :
2460 : ! real
2461 10 : ABI_SFREE(a2f_tr%n0)
2462 10 : ABI_SFREE(a2f_tr%omega)
2463 10 : ABI_SFREE(a2f_tr%vals_in)
2464 10 : ABI_SFREE(a2f_tr%vals_out)
2465 10 : ABI_SFREE(a2f_tr%vals_tr)
2466 10 : ABI_SFREE(a2f_tr%vals_tr_gen)
2467 10 : ABI_SFREE(a2f_tr%lambdaw_tr)
2468 :
2469 10 : end subroutine a2fw_tr_free
2470 : !!***
2471 :
2472 : !----------------------------------------------------------------------
2473 :
2474 : !!****f* m_phgamma/a2fw_tr_init
2475 : !! NAME
2476 : !! a2fw_tr_init
2477 : !!
2478 : !! FUNCTION
2479 : !! Calculates the FS averaged alpha^2F_tr,in,out(w) functions
2480 : !!
2481 : !! INPUTS
2482 : !! cryst<crystal_t>=Info on the unit cell.
2483 : !! ifc<ifc_type>=Interatomic force constants.
2484 : !! gams<phgamma_t>=Structure storing the phonon linewidths.
2485 : !! wstep=Step for linear frequency mesh in Ha.
2486 : !! wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w).
2487 : !! ph_intmeth=Integration method for phonons: 1 for gaussian, 2 for tetrahedra
2488 : !! smear=Gaussian broadening used to approximate the Dirac delta.
2489 : !! ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift).
2490 : !! nqshift=Number of shifts used to generated the Q-mesh.
2491 : !! qshift(3,nqshift)=The shifts.
2492 : !! comm=MPI communicator
2493 : !! [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop are ignored and
2494 : !! A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation.
2495 : !! [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account
2496 : !! the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1
2497 : !! Other values of qptopt can be used for debugging purpose.
2498 : !!
2499 : !! OUTPUT
2500 : !! a2f_tr<a2fw_tr_t>=Structure storing the Eliashberg transport function a2F_tr(w).
2501 : !!
2502 : !! SOURCE
2503 :
2504 10 : subroutine a2fw_tr_init(a2f_tr, gams, cryst, ifc, ph_intmeth, wstep, wminmax, smear, ngqpt, nqshift, qshift, comm,&
2505 : qintp, qptopt) ! optional
2506 :
2507 : !Arguments ------------------------------------
2508 : !scalars
2509 : class(a2fw_tr_t),target,intent(out) :: a2f_tr
2510 : integer,intent(in) :: ph_intmeth,nqshift,comm
2511 : integer,intent(in),optional :: qptopt
2512 : real(dp),intent(in) :: wstep,smear
2513 : logical,optional,intent(in) :: qintp
2514 : type(phgamma_t),intent(inout) :: gams
2515 : type(ifc_type),intent(in) :: ifc
2516 : type(crystal_t),intent(in) :: cryst
2517 : !arrays
2518 : integer,intent(in) :: ngqpt(3)
2519 : real(dp),intent(in) :: wminmax(2),qshift(3,nqshift)
2520 :
2521 : !Local variables -------------------------
2522 : !scalars
2523 : integer,parameter :: master = 0
2524 : integer :: my_qptopt,iq_ibz,nqibz,ount,my_rank,nproc,cnt
2525 : integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz, idir, jdir
2526 : real(dp) :: omega,xx,omega_min,omega_max,ww, cpu, wall, gflops
2527 : logical :: do_qintp
2528 : character(len=500) :: msg
2529 10 : type(htetra_t) :: tetra_q
2530 : !arrays
2531 : integer :: qptrlatt(3,3),new_qptrlatt(3,3)
2532 20 : real(dp) :: lambda_iso(3,3), omega_log(3,3), phfrq(gams%natom3)
2533 20 : real(dp) :: gamma_in_ph(3,3,gams%natom3), gamma_out_ph(3,3,gams%natom3)
2534 30 : real(dp) :: lambda_in_ph(3,3,gams%natom3), lambda_out_ph(3,3,gams%natom3)
2535 10 : real(dp),allocatable :: tmp_a2f_in(:,:,:), tmp_a2f_out(:,:,:), my_qshift(:,:)
2536 10 : real(dp), contiguous, pointer :: a2f_tr_1d(:)
2537 10 : real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:)
2538 10 : real(dp),allocatable :: a2f_tr_1mom(:),a2f_tr_logmom(:),a2f_tr_logmom_int(:),wdt(:,:)
2539 10 : real(dp),allocatable :: lambda_in_tetra(:,:,:,:,:),phfreq_tetra(:,:,:), lambda_out_tetra(:,:,:,:,:)
2540 : ! *********************************************************************
2541 :
2542 10 : my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt
2543 10 : do_qintp = .True.; if (present(qintp)) do_qintp = qintp
2544 :
2545 10 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
2546 10 : nsppol = gams%nsppol; natom3 = gams%natom3
2547 :
2548 10 : call cwtime(cpu, wall, gflops, "start")
2549 :
2550 10 : if (do_qintp) then
2551 : ! Generate the q-mesh by finding the IBZ and the corresponding weights.
2552 5 : qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3)
2553 :
2554 : call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, &
2555 5 : new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift)
2556 5 : ABI_FREE(qbz)
2557 :
2558 : ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
2559 20 : a2f_tr%ngqpt = ngqpt; a2f_tr%nqshift = size(my_qshift, dim=2)
2560 15 : ABI_MALLOC(a2f_tr%qshift, (3, a2f_tr%nqshift))
2561 30 : a2f_tr%qshift = my_qshift
2562 5 : ABI_FREE(my_qshift)
2563 :
2564 : else
2565 : ! No interpolation. Use q-mesh parameters from gams%
2566 20 : a2f_tr%ngqpt = gams%ngqpt; a2f_tr%nqshift = 1
2567 5 : nqibz = gams%nqibz
2568 15 : ABI_MALLOC(qibz, (3,nqibz))
2569 15 : ABI_MALLOC(wtq, (nqibz))
2570 5 : ABI_MALLOC(a2f_tr%qshift, (3, a2f_tr%nqshift))
2571 320 : qibz = gams%qibz; wtq = gams%wtq
2572 25 : a2f_tr%qshift = zero ! Note: assuming q-mesh centered on Gamma.
2573 : end if
2574 :
2575 10 : call cwtime_report(" a2fw_tr_init, q-setup", cpu, wall, gflops)
2576 :
2577 : ! Define Min and max frequency for the mesh (enlarge it a bit)
2578 10 : omega_min = wminmax(1); omega_max = wminmax(2)
2579 10 : omega_min = omega_min - 0.1*abs(omega_min)
2580 10 : if (omega_min >= zero) omega_min = one/Ha_meV
2581 10 : omega_max = omega_max + 0.1*abs(omega_max)
2582 :
2583 10 : a2f_tr%nsppol = nsppol; a2f_tr%natom3 = gams%natom3; a2f_tr%smear = smear
2584 10 : a2f_tr%omega_min = omega_min; a2f_tr%omega_max = omega_max
2585 10 : nomega = int((omega_max - omega_min) / wstep); a2f_tr%nomega = nomega; a2f_tr%wstep = wstep
2586 :
2587 30 : ABI_MALLOC(a2f_tr%n0, (nsppol))
2588 30 : a2f_tr%n0 = gams%n0
2589 : ! Build linear mesh.
2590 30 : ABI_MALLOC(a2f_tr%omega, (nomega))
2591 4750 : a2f_tr%omega = arth(omega_min, wstep, nomega)
2592 :
2593 170870 : ABI_CALLOC(a2f_tr%vals_in, (nomega, 3, 3, 0:natom3, nsppol))
2594 170850 : ABI_CALLOC(a2f_tr%vals_out, (nomega,3, 3, 0:natom3, nsppol))
2595 170850 : ABI_CALLOC(a2f_tr%vals_tr, (nomega,3, 3, 0:natom3, nsppol))
2596 170850 : ABI_CALLOC(a2f_tr%lambdaw_tr, (nomega, 3, 3, 0:natom3, nsppol))
2597 30 : ABI_MALLOC(tmp_a2f_in, (nomega, 3, 3))
2598 20 : ABI_MALLOC(tmp_a2f_out, (nomega, 3, 3))
2599 :
2600 10 : if (ph_intmeth == 2) then
2601 8 : call cwtime(cpu, wall, gflops, "start")
2602 :
2603 : ! Prepare tetrahedron integration.
2604 8 : qptrlatt = 0; qptrlatt(1, 1) = a2f_tr%ngqpt(1); qptrlatt(2, 2) = a2f_tr%ngqpt(2); qptrlatt(3, 3) = a2f_tr%ngqpt(3)
2605 8 : tetra_q = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f_tr%nqshift, a2f_tr%qshift, nqibz, qibz, comm, msg, ierr)
2606 8 : if (ierr/=0) ABI_ERROR(msg)
2607 :
2608 48 : ABI_MALLOC_OR_DIE(lambda_in_tetra, (nqibz, 3, 3, natom3, nsppol), ierr)
2609 32 : ABI_MALLOC_OR_DIE(lambda_out_tetra, (nqibz, 3, 3, natom3, nsppol), ierr)
2610 17419 : lambda_in_tetra = zero
2611 17419 : lambda_out_tetra = zero
2612 :
2613 40 : ABI_MALLOC_OR_DIE(phfreq_tetra, (nqibz, natom3, nsppol), ierr)
2614 1939 : phfreq_tetra = zero
2615 :
2616 16 : call cwtime_report(" a2fw_tr_init%tetra,", cpu, wall, gflops)
2617 : end if
2618 :
2619 10 : call cwtime(cpu, wall, gflops, "start")
2620 :
2621 : ! Loop over spins and qpoints in the IBZ
2622 10 : cnt = 0
2623 20 : do spin=1,nsppol
2624 806 : do iq_ibz=1,nqibz
2625 786 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle
2626 :
2627 : ! Interpolate or evaluate lambda transport coefficients.
2628 786 : if (do_qintp) then
2629 725 : call phgamma_vv_interp(gams,cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
2630 : else
2631 61 : call phgamma_vv_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph)
2632 : end if
2633 :
2634 10 : select case (ph_intmeth)
2635 : case (1)
2636 : ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f_tr
2637 : ! (note that unstable modes are included).
2638 612 : do mu=1,natom3
2639 1959930 : tmp_a2f_in = zero
2640 1959930 : tmp_a2f_out = zero
2641 217566 : do iw=1,nomega
2642 217107 : xx = a2f_tr%omega(iw) - phfrq(mu)
2643 2822391 : tmp_a2f_in(iw,:,:) = tmp_a2f_in(iw,:,:) + gaussian(xx, smear) * lambda_in_ph(:,:,mu) * abs(phfrq(mu))
2644 2822850 : tmp_a2f_out(iw,:,:) = tmp_a2f_out(iw,:,:) + gaussian(xx, smear) * lambda_out_ph(:,:,mu) * abs(phfrq(mu))
2645 : end do
2646 1959930 : a2f_tr%vals_in(:,:,:,mu,spin) = a2f_tr%vals_in(:,:,:,mu,spin) + tmp_a2f_in(:,:,:) * wtq(iq_ibz)
2647 1960083 : a2f_tr%vals_out(:,:,:,mu,spin) = a2f_tr%vals_out(:,:,:,mu,spin) + tmp_a2f_out(:,:,:) * wtq(iq_ibz)
2648 : end do
2649 :
2650 : case (2)
2651 : ! Tetra: store data.
2652 2685 : do mu=1,natom3
2653 24687 : lambda_in_tetra(iq_ibz, :,:, mu, spin) = lambda_in_ph(:,:,mu) * abs(phfrq(mu))
2654 24687 : lambda_out_tetra(iq_ibz, :,:, mu, spin) = lambda_out_ph(:,:,mu) * abs(phfrq(mu))
2655 2532 : phfreq_tetra(iq_ibz, mu, spin) = phfrq(mu)
2656 : end do
2657 : end select
2658 :
2659 : end do ! iq_ibz
2660 : end do ! spin
2661 :
2662 10 : if (ph_intmeth == 2) then
2663 : ! Collect results on each node.
2664 8 : call xmpi_sum(lambda_in_tetra, comm, ierr)
2665 8 : call xmpi_sum(lambda_out_tetra, comm, ierr)
2666 8 : call xmpi_sum(phfreq_tetra, comm, ierr)
2667 :
2668 : ! workspace for tetra.
2669 24 : ABI_MALLOC(wdt, (nomega, 2))
2670 :
2671 : ! TODO: with the tetra_get_onewk call we can integrate this above
2672 : ! and avoid allocating all of lambda_in_tetra and phfreq_tetra!!!
2673 : ! For each mode get its contribution
2674 8 : cnt = 0
2675 16 : do spin=1,nsppol
2676 40 : do mu=1,natom3
2677 1931 : do iq_ibz=1,nqibz
2678 1899 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! mpi-parallelism
2679 :
2680 1899 : call tetra_q%get_onewk(iq_ibz, gams%bcorr, nomega, nqibz, phfreq_tetra(:,mu,spin), omega_min, omega_max, one, wdt)
2681 1802151 : wdt = wdt * wtq(iq_ibz)
2682 :
2683 : ! Accumulate (Integral of a2F_tr is computed afterwards)
2684 7620 : do idir=1,3
2685 24687 : do jdir=1,3
2686 : a2f_tr%vals_in(:,idir,jdir,mu,spin) = a2f_tr%vals_in(:,idir,jdir,mu,spin) &
2687 8101134 : + wdt(:,1) * lambda_in_tetra(iq_ibz,idir,jdir,mu,spin)
2688 : a2f_tr%vals_out(:,idir,jdir,mu,spin) = a2f_tr%vals_out(:,idir,jdir,mu,spin) &
2689 8106831 : + wdt(:,1) * lambda_out_tetra(iq_ibz,idir,jdir,mu,spin)
2690 : !a2f_tr%lambdaw(:,mu,spin) = a2f_tr%lambdaw(:,mu,spin) + wdt(:,2) * lambda_XX_tetra(iq_ibz, mu, spin)
2691 : end do
2692 : end do
2693 : end do
2694 : end do
2695 : end do
2696 :
2697 8 : ABI_FREE(wdt)
2698 8 : ABI_FREE(lambda_in_tetra)
2699 8 : ABI_FREE(lambda_out_tetra)
2700 8 : ABI_FREE(phfreq_tetra)
2701 16 : call tetra_q%free()
2702 : end if
2703 :
2704 : ! Collect final results on each node
2705 10 : call xmpi_sum(a2f_tr%vals_in, comm, ierr)
2706 10 : call xmpi_sum(a2f_tr%vals_out, comm, ierr)
2707 :
2708 : ! TODO: divide by g(eF, spin)
2709 : ! check normalization with N(0) and <v^2> from Savrasov eq 20
2710 20 : do spin=1,nsppol
2711 170410 : a2f_tr%vals_in(:,:,:,0,spin) = sum(a2f_tr%vals_in(:,:,:,1:natom3,spin), dim=4)
2712 170420 : a2f_tr%vals_out(:,:,:,0,spin) = sum(a2f_tr%vals_out(:,:,:,1:natom3,spin), dim=4)
2713 : !a2f_tr%vals(:,:,spin) = a2f_tr%vals(:,:,spin) / (two_pi * a2f_tr%n0(spin))
2714 : end do
2715 :
2716 170830 : a2f_tr%vals_tr = a2f_tr%vals_out - a2f_tr%vals_in
2717 :
2718 10 : call cwtime_report(" a2fw_tr_init, a2f_eval", cpu, wall, gflops)
2719 :
2720 : ! NOTE: for the moment calculate the log moms etc on just the trace of the a2f functions
2721 20 : ABI_MALLOC(a2f_tr_1mom, (nomega))
2722 20 : ABI_MALLOC(a2f_tr_logmom, (nomega))
2723 20 : ABI_MALLOC(a2f_tr_logmom_int, (nomega))
2724 :
2725 : !call a2fw_tr_print_info()
2726 :
2727 : ! Compute lambda_tr(w) (resolved in spin mode and directions)
2728 20 : do spin=1,nsppol
2729 60 : do mu=0,natom3
2730 170 : do idir=1,3
2731 520 : do jdir=1,3
2732 :
2733 170640 : do iw=1,nomega
2734 170280 : ww = a2f_tr%omega(iw)
2735 170640 : if (abs(ww) > EPHTK_WTOL) then
2736 170280 : a2f_tr_1mom(iw) = a2f_tr%vals_tr(iw,idir,jdir,mu,spin) / abs(ww)
2737 : else
2738 0 : a2f_tr_1mom(iw) = zero
2739 : end if
2740 : ! TODO: Strange that the first value of int(f) is not set to zero!
2741 : ! FIXME: The frequency integration overestimates lambda if there are negative phonon frequencies!
2742 : end do
2743 :
2744 480 : call simpson_int(nomega, wstep, a2f_tr_1mom, a2f_tr%lambdaw_tr(:,idir,jdir,mu,spin))
2745 : end do ! jdir
2746 : end do ! idir
2747 : end do ! mu
2748 : end do ! spin
2749 :
2750 10 : if (my_rank == master) then
2751 10 : ount = std_out
2752 10 : if (nsppol > 1) then
2753 0 : write(msg,'(3a)') ch10,'Comment: some of the following quantities should be integrated over spin', ch10
2754 0 : call wrtout(ount, msg)
2755 : end if
2756 : end if
2757 :
2758 20 : do spin=1,nsppol
2759 10 : lambda_iso = a2fw_tr_moment(a2f_tr, 0, spin)
2760 :
2761 40 : do idir=1,3
2762 130 : do jdir=1,3
2763 90 : a2f_tr_1d => a2f_tr%vals_tr(:,idir,jdir,0,spin)
2764 :
2765 : ! Get log moment of alpha^2F.
2766 42660 : a2f_tr_logmom = zero
2767 42660 : do iw=1,nomega
2768 42570 : omega = a2f_tr%omega(iw)
2769 42660 : if (abs(omega) > EPHTK_WTOL .and. abs(lambda_iso(idir,jdir)) > EPHTK_WTOL) then
2770 42570 : a2f_tr_logmom(iw) = (two/lambda_iso(idir,jdir)) * a2f_tr_1d(iw)*log(abs(omega))/abs(omega)
2771 : end if
2772 : end do
2773 90 : call simpson_int(nomega, wstep, a2f_tr_logmom, a2f_tr_logmom_int)
2774 : !write(std_out,*)' iw,nomega,greatest_real,a2f_tr_logmom_int(nomega)=',& iw,nomega,greatest_real,a2f_tr_logmom_int(nomega)
2775 120 : if(abs(a2f_tr_logmom_int(nomega)) < log(greatest_real*tol6)) then
2776 90 : omega_log(idir,jdir) = exp(a2f_tr_logmom_int(nomega))
2777 : else
2778 0 : omega_log(idir,jdir)=greatest_real*tol6
2779 : endif
2780 : end do
2781 : end do
2782 :
2783 : ! TODO: make output only for irred values xx yy zz and top half of matrix
2784 20 : if (my_rank == master) then
2785 10 : if (do_qintp) then
2786 5 : write(ount,'(a)')' Evaluation of parameters analogous to electron-phonon coupling for 3x3 directions (interpolated) '
2787 : else
2788 5 : write(ount,'(a)')' Evaluation of parameters analogous to electron-phonon coupling for 3x3 directions (coarse grid) '
2789 : endif
2790 10 : write(ount,'(a,3(3es10.3,2x))') ' lambda = ',lambda_iso
2791 10 : write(ount,'(a,3(3es10.3,2x),a)' )' omegalog = ',omega_log,' (Ha) '
2792 130 : write(ount,'(a,3(3es10.3,2x),a)' )' ',omega_log*Ha_K, ' (Kelvin) '
2793 10 : write(ount,"(a)")' positive moments of alpha2Ftr:'
2794 130 : write(ount,'(a,9(es10.3,2x))' )' lambda <omega^2> = ',a2fw_tr_moment(a2f_tr, 2, spin)
2795 130 : write(ount,'(a,9(es10.3,2x))' )' lambda <omega^3> = ',a2fw_tr_moment(a2f_tr, 3, spin)
2796 130 : write(ount,'(a,9(es10.3,2x))' )' lambda <omega^4> = ',a2fw_tr_moment(a2f_tr, 4, spin)
2797 130 : write(ount,'(a,9(es10.3,2x))' )' lambda <omega^5> = ',a2fw_tr_moment(a2f_tr, 5, spin)
2798 : end if
2799 : end do
2800 :
2801 10 : ABI_FREE(tmp_a2f_in)
2802 10 : ABI_FREE(tmp_a2f_out)
2803 10 : ABI_FREE(a2f_tr_1mom)
2804 10 : ABI_FREE(a2f_tr_logmom)
2805 10 : ABI_FREE(a2f_tr_logmom_int)
2806 10 : ABI_FREE(qibz)
2807 10 : ABI_FREE(wtq)
2808 :
2809 70 : end subroutine a2fw_tr_init
2810 : !!***
2811 :
2812 : !----------------------------------------------------------------------
2813 :
2814 : !!****f* m_phgamma/a2fw_tr_write
2815 : !! NAME
2816 : !! a2fw_tr_write
2817 : !!
2818 : !! FUNCTION
2819 : !! Write alpha^2F_tr(w) to an external file in text form
2820 : !!
2821 : !! INPUTS
2822 : !! basename=Filename for output.
2823 : !! post=String appended to netcdf variables e.g. _qcoarse, _qintp
2824 : !! ncid=Netcdf file handler. Set it to nctk_noid to disable output.
2825 : !!
2826 : !! OUTPUT
2827 : !! Output is written to file. This routine should be called by one MPI proc.
2828 : !!
2829 : !! SOURCE
2830 :
2831 10 : subroutine a2fw_tr_write(a2f_tr, basename, post, ncid)
2832 :
2833 : !Arguments ------------------------------------
2834 : !scalars
2835 : class(a2fw_tr_t),intent(in) :: a2f_tr
2836 : integer,intent(in) :: ncid
2837 : character(len=*),intent(in) :: basename, post
2838 :
2839 : !Local variables -------------------------
2840 : !scalars
2841 : integer :: iw,spin,unt,ii,mu,idir, jdir, ncerr
2842 : character(len=500) :: dim1_name, msg
2843 : character(len=fnlen) :: path
2844 : ! *********************************************************************
2845 :
2846 : ! Write spin-resolved a2F_tr(w)
2847 10 : path = strcat(basename, "_A2FW_tr")
2848 10 : if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then
2849 0 : ABI_ERROR(msg)
2850 : end if
2851 :
2852 10 : call write_a2fw_tr_header()
2853 :
2854 : ! Reminder: vals_tr(nomega,3,3,0:natom3,nsppol)
2855 :
2856 10 : if (a2f_tr%nsppol == 1) then
2857 10 : write(unt,'(a)')"# Frequency, a2F_tr(w,i,j)"
2858 4740 : do iw=1,a2f_tr%nomega
2859 4740 : write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,0,1)
2860 : end do
2861 10 : write(unt,'(3a)')ch10, ch10, "# Frequency, lambda(w,i,j)"
2862 4740 : do iw=1,a2f_tr%nomega
2863 4740 : write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,0,1)
2864 : end do
2865 :
2866 : else
2867 0 : write(unt,'(a)')"# Frequency, a2F_tr_tot(w,i,j), a2F_tr_spin1(w,i,j) ..."
2868 0 : do iw=1,a2f_tr%nomega
2869 0 : write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), &
2870 0 : ((sum(a2f_tr%vals_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , &
2871 0 : a2f_tr%vals_tr(iw,:,:,0,1), & ! UP
2872 0 : a2f_tr%vals_tr(iw,:,:,0,2) ! DOWN
2873 : end do
2874 0 : write(unt,'(3a)') ch10, ch10, "# Frequency, lambda_tr_tot(w,i,j)dw, lambda_tr_spin1(w,i,j) ..."
2875 0 : do iw=1,a2f_tr%nomega
2876 0 : write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), &
2877 0 : ((sum(a2f_tr%lambdaw_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , & ! TOT
2878 0 : a2f_tr%lambdaw_tr(iw,:,:,0,1), & ! UP
2879 0 : a2f_tr%lambdaw_tr(iw,:,:,0,2) ! DOWN
2880 : end do
2881 : end if
2882 :
2883 10 : close(unt)
2884 :
2885 : ! Write phonon mode contributions to a2F_tr(w,i,j)
2886 10 : path = strcat(basename, "_PH_A2FW_tr")
2887 10 : if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then
2888 0 : ABI_ERROR(msg)
2889 : end if
2890 10 : call write_a2fw_tr_header()
2891 :
2892 : ! Reminder: vals_tr(nomega,3,3,0:natom3,nsppol)
2893 10 : if (a2f_tr%nsppol == 1) then
2894 50 : do mu=0,a2f_tr%natom3
2895 40 : write(unt,'(a,i0)')"# Phonon mode ",mu
2896 40 : write(unt,'(a)')"# Frequency, a2F_tr(w)"
2897 18960 : do iw=1,a2f_tr%nomega
2898 18960 : write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,mu,1)
2899 : end do
2900 120 : do ii=1,2; write(unt,'(a)')""; end do
2901 40 : write(unt,'(a)')"# Frequency, lambda_tr(w)"
2902 18960 : do iw=1,a2f_tr%nomega
2903 18960 : write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,mu,1)
2904 : end do
2905 130 : do ii=1,2; write(unt,'(a)')""; end do
2906 : end do
2907 :
2908 : else
2909 0 : do mu=0,a2f_tr%natom3
2910 0 : write(unt,'(a,i0)')"# Phonon mode ",mu
2911 0 : write(unt,'(a)')"# Frequency, a2F_tot(w), a2F_spin1(w) ..."
2912 0 : do iw=1,a2f_tr%nomega
2913 0 : write(unt,*) a2f_tr%omega(iw), &
2914 0 : sum(a2f_tr%vals_tr(iw,:,:,mu,:)), & ! TOT
2915 0 : a2f_tr%vals_tr(iw,:,:,mu,1) , & ! UP
2916 0 : a2f_tr%vals_tr(iw,:,:,mu,2) ! DOWN
2917 : end do
2918 0 : write(unt,'(3a)')ch10,ch10,"# Frequency, lambda_tot(w)dw, lambda_spin1(w) ..."
2919 0 : do iw=1,a2f_tr%nomega
2920 0 : write(unt,*) a2f_tr%omega(iw), &
2921 0 : sum(a2f_tr%lambdaw_tr(iw,:,:,mu,:)), & ! TOT
2922 0 : a2f_tr%lambdaw_tr(iw,:,:,mu,1), & ! UP
2923 0 : a2f_tr%lambdaw_tr(iw,:,:,mu,2) ! DOWN
2924 : end do
2925 : end do
2926 : end if
2927 :
2928 10 : close(unt)
2929 :
2930 10 : if (ncid /= nctk_noid) then
2931 : ! Define dimensions.
2932 10 : dim1_name = strcat("a2ftr_nomega", post)
2933 : ncerr = nctk_def_dims(ncid, [ &
2934 : nctkdim_t(dim1_name, a2f_tr%nomega), &
2935 : nctkdim_t("natom3p1", a2f_tr%natom3 + 1), nctkdim_t("number_of_spins", a2f_tr%nsppol) &
2936 40 : ], defmode=.True.)
2937 10 : NCF_CHECK(ncerr)
2938 :
2939 : ncerr = nctk_def_arrays(ncid, [ &
2940 : nctkarr_t(strcat('a2ftr_mesh', post), "dp", dim1_name), &
2941 : nctkarr_t(strcat('a2ftr_values', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")), &
2942 : nctkarr_t(strcat('a2ftr_lambdaw', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")) &
2943 40 : ])
2944 10 : NCF_CHECK(ncerr)
2945 :
2946 10 : NCF_CHECK(nctk_set_datamode(ncid))
2947 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_mesh", post)), a2f_tr%omega))
2948 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_values", post)), a2f_tr%vals_tr))
2949 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_lambdaw", post)), a2f_tr%lambdaw_tr))
2950 : end if
2951 :
2952 : contains
2953 :
2954 20 : subroutine write_a2fw_tr_header()
2955 :
2956 : ! Output the header.
2957 20 : write(unt,'(a)') '#'
2958 20 : write(unt,'(a)') '# ABINIT package: a2F_tr(w) file'
2959 20 : write(unt,'(a)') '#'
2960 20 : write(unt,'(a)') '# a2F_tr(w) function integrated over the FS. omega in a.u.'
2961 20 : write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f_tr%ngqpt))
2962 20 : write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f_tr%nomega," between omega_min: ",a2f_tr%omega_min,&
2963 40 : ' Ha and omega_max: ',a2f_tr%omega_max,' Ha with step:',a2f_tr%wstep
2964 20 : write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f_tr%smear
2965 20 : if (a2f_tr%nsppol == 2) then
2966 0 : write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f_tr%n0)
2967 : end if
2968 40 : do spin=1,a2f_tr%nsppol
2969 40 : write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f_tr%n0(spin)
2970 : end do
2971 60 : do ii=1,2; write(unt,'(a)') "# "; end do
2972 :
2973 20 : end subroutine write_a2fw_tr_header
2974 :
2975 : end subroutine a2fw_tr_write
2976 : !!***
2977 :
2978 : !----------------------------------------------------------------------
2979 :
2980 : !!****f* m_phgamma/eph_phgamma
2981 : !! NAME
2982 : !! eph_phgamma
2983 : !!
2984 : !! FUNCTION
2985 : !! Compute phonon linewidths in metals.
2986 : !!
2987 : !! INPUTS
2988 : !! wk0_path=String with the path to the GS unperturbed WFK file.
2989 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
2990 : !! dtset<dataset_type>=All input variables for this dataset.
2991 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
2992 : !! dvdb<dbdb_type>=Database with the DFPT SCF potentials.
2993 : !! ifc<ifc_type>=interatomic force constants and corresponding real space grid info.
2994 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
2995 : !! pawang<pawang_type)>=PAW angular mesh and related data.
2996 : !! pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data.
2997 : !! pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data.
2998 : !! psps<pseudopotential_type>=Variables related to pseudopotentials.
2999 : !! comm=MPI communicator.
3000 : !!
3001 : !! OUTPUT
3002 : !!
3003 : !! SOURCE
3004 :
3005 5 : subroutine eph_phgamma(wfk0_path, dtfil, ngfft, ngfftf, dtset, cryst, ebands, dvdb, ifc, &
3006 5 : pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm)
3007 :
3008 : !Arguments ------------------------------------
3009 : !scalars
3010 : character(len=*),intent(in) :: wfk0_path
3011 : integer,intent(in) :: comm
3012 : type(datafiles_type),intent(in) :: dtfil
3013 : type(dataset_type),intent(in) :: dtset
3014 : type(crystal_t),intent(in) :: cryst
3015 : type(ebands_t),intent(in) :: ebands
3016 : type(dvdb_t),intent(inout) :: dvdb
3017 : type(pawang_type),intent(in) :: pawang
3018 : type(pseudopotential_type),intent(in) :: psps
3019 : type(pawfgr_type),intent(in) :: pawfgr
3020 : type(ifc_type),intent(in) :: ifc
3021 : type(mpi_type),intent(in) :: mpi_enreg
3022 : !arrays
3023 : integer,intent(in) :: ngfft(18),ngfftf(18)
3024 : type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw)
3025 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
3026 :
3027 : !Local variables ------------------------------
3028 : !scalars
3029 : integer,parameter :: tim_getgh1c = 1, berryopt0 = 0, qptopt1 = 1, master = 0, ndat1 = 1, eph_scalprod0 = 0
3030 : integer :: my_rank,nproc,mband,nsppol,nkibz,idir,ipert,iq_ibz
3031 : integer :: cplex,db_iqpt,natom,natom3,ipc,ipc1,ipc2,nspinor,onpw
3032 : integer :: bstart_k,bstart_kq,nb_k,nb_kq,band_k, band_kq, in_k, im_kq
3033 : integer :: ik_ibz,ik_bz,ikq_bz,ikq_ibz,isym_k,isym_kq,trev_k,trev_kq !,timrev_q
3034 : integer :: ik_fs, my_ik, my_is, spin, istwf_k, istwf_kq, npw_k, npw_kq
3035 : integer :: ii,jj,ipw,mpw,my_mpw,mnb,ierr,cnt,ncid
3036 : integer :: n1, n2, n3, n4, n5, n6, nspden, ltetra
3037 : integer :: sij_opt, usecprj, usevnl, optlocal, optnl, opt_gvnlx1
3038 : integer :: nfft,nfftf,mgfft,mgfftf,nkpg
3039 : integer :: jene, iene, comm_rpt, nesting, my_npert, my_ip, my_iq, ncerr, edos_intmeth
3040 : real(dp) :: cpu, wall, gflops, cpu_q, wall_q, gflops_q, cpu_k, wall_k, gflops_k, cpu_all, wall_all, gflops_all
3041 : real(dp) :: edos_step, edos_broad, sigma, ecut, eig0nk ! eshift,
3042 : logical :: gen_eigenpb, need_velocities, isirr_k, isirr_kq, print_time_k, need_ftinterp
3043 5 : type(wfd_t) :: wfd
3044 : type(fstab_t),pointer :: fs
3045 5 : type(gs_hamiltonian_type) :: gs_hamkq
3046 5 : type(rf_hamiltonian_type) :: rf_hamkq
3047 5 : type(edos_t) :: edos
3048 5 : type(phgamma_t) :: gams
3049 5 : type(a2fw_t) :: a2fw
3050 5 : type(a2fw_tr_t) :: a2fw_tr
3051 50 : type(ddkop_t) :: ddkop
3052 : type(xcomm_t) :: pert_comm, qs_comm, qpt_comm, bsum_comm, kpt_comm, spin_comm, pkb_comm !, ncwrite_comm
3053 5 : type(krank_t) :: krank
3054 5 : type(htetra_t) :: tetra
3055 : !type(lgroup_t) :: lg_q
3056 : character(len=500) :: msg
3057 : character(len=fnlen) :: path
3058 : !arrays
3059 : integer :: g0_k(3),g0bz_kq(3),g0_kq(3),units(2)
3060 : integer :: indkk_kq(6,1), work_ngfft(18),gmax(3),my_gmax(3),gamma_ngqpt(3) !g0ibz_kq(3),
3061 5 : integer,allocatable :: kg_k(:,:),kg_kq(:,:),gtmp(:,:),nband(:,:),wfd_istwfk(:), my_pinfo(:,:), pert_table(:,:) !, qibz_done(:)
3062 10 : real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3),qpt(3), lf(2),rg(2),res(2), vk(3), vkq(3), wminmax(2), n0(ebands%nsppol)
3063 10 : real(dp) :: resvv_in(2,9), resvv_out(2,9), phfrq(3*cryst%natom), abc(3)
3064 5 : real(dp),allocatable :: kinpw_k(:), kinpw_kq(:), displ_cart(:,:,:,:), displ_red(:,:,:,:)
3065 5 : real(dp),allocatable :: grad_berry(:,:), kpg_kq(:,:), kpg_k(:,:)
3066 5 : real(dp),allocatable :: ffnl_k(:,:,:,:), ffnl_kq(:,:,:,:), ph3d_k(:,:,:), ph3d_kq(:,:,:)
3067 5 : real(dp),allocatable :: v1scf(:,:,:,:), tgam(:,:,:), gkq_atm(:,:,:,:), gkq_atm_ipc(:,:,:), lambda(:)
3068 10 : real(dp),allocatable :: bras_kq(:,:,:), kets_k(:,:,:), h1_kets_kq(:,:,:), cg_work(:,:)
3069 5 : real(dp),allocatable :: ph1d(:,:), vlocal(:,:,:,:), vlocal1(:,:,:,:,:)
3070 5 : real(dp),allocatable :: dummy_vtrial(:,:), gvnlx1(:,:,:), work(:,:,:,:)
3071 5 : real(dp),allocatable :: gs1c_kq(:,:), v1_work(:,:,:,:), vcart_ibz(:,:,:,:)
3072 10 : real(dp),allocatable :: wt_ek(:,:), wt_ekq(:,:), dbldelta_wts(:,:)
3073 5 : real(dp),allocatable :: tgamvv_in(:,:,:,:), vv_kk(:,:,:), tgamvv_out(:,:,:,:), vv_kkq(:,:,:), tmp_vals_ee(:,:,:,:,:), emesh(:)
3074 10 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
3075 5 : type(fstab_t),target,allocatable :: fstab(:)
3076 5 : type(pawcprj_type),allocatable :: cwaveprj0(:,:)
3077 : #ifdef HAVE_MPI
3078 : integer :: ndims, comm_cart, me_cart, coords(5)
3079 : logical :: reorder
3080 5 : integer,allocatable :: dims(:)
3081 5 : logical,allocatable :: periods(:), keepdim(:)
3082 : #endif
3083 : !************************************************************************
3084 :
3085 5 : if (psps%usepaw == 1) then
3086 0 : ABI_ERROR("PAW not implemented")
3087 : ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/))
3088 : end if
3089 :
3090 5 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
3091 15 : units = [std_out, ab_out]
3092 5 : call wrtout(std_out, " Computing ph linewidths in metals.")
3093 5 : call cwtime(cpu_all, wall_all, gflops_all, "start")
3094 5 : call pstat_proc%print(_PSTAT_ARGS_)
3095 :
3096 : ! Copy important dimensions
3097 5 : natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden
3098 5 : nkibz = ebands%nkpt; mband = ebands%mband
3099 :
3100 : ! FFT meshes
3101 20 : nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3))
3102 35 : nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3))
3103 5 : n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3)
3104 5 : n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6)
3105 :
3106 : ! Compute electron DOS.
3107 : !call dtset%get_edos_params(edos_intmeth, edos_step, edos_broad)
3108 5 : edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos
3109 : edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear
3110 5 : edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha
3111 5 : edos = ebands%get_edos(cryst, edos_intmeth, edos_step, edos_broad, comm)
3112 :
3113 : ! Store DOS per spin channel.
3114 10 : n0(:) = edos%gef(1:edos%nsppol)
3115 5 : if (my_rank == master) then
3116 5 : call edos%print(units)
3117 5 : call edos%write(strcat(dtfil%filnam_ds(4), "_EDOS"))
3118 : end if
3119 :
3120 : ! Find Fermi surface k-points
3121 : ! TODO: support kptopt, change setup of k-points if tetra: fist tetra weights then k-points on the Fermi surface!
3122 20 : ABI_MALLOC(fstab, (nsppol))
3123 5 : call fstab_init(fstab, ebands, cryst, dtset, tetra, comm)
3124 5 : call tetra%free()
3125 :
3126 : ! Define q-mesh. eph_ngqpt_fine activates the Fourier interpolation of the DFPT potentials.
3127 26 : gamma_ngqpt = ifc%ngqpt; if (all(dtset%eph_ngqpt_fine /= 0)) gamma_ngqpt = dtset%eph_ngqpt_fine
3128 :
3129 5 : call phgamma_init(gams, cryst, ifc, ebands, fstab(1), dtset, eph_scalprod0, gamma_ngqpt, n0, comm)
3130 :
3131 5 : call wrtout(std_out, sjoin(" q-mesh for phonon linewidths:", ltoa(gamma_ngqpt)))
3132 5 : call wrtout(std_out, sjoin(" Will compute:", itoa(gams%nqibz), "q-points in the IBZ"))
3133 5 : call pstat_proc%print(_PSTAT_ARGS_)
3134 :
3135 : ! Select option for double delta with tetra.
3136 : ! 2 for the optimized tetrahedron method.
3137 : ! -2 for the linear tetrahedron method.
3138 5 : ltetra = 0
3139 5 : if (dtset%eph_intmeth == 2) ltetra = 2
3140 5 : if (dtset%eph_intmeth == -2) ltetra = 1
3141 :
3142 : ! Compute optimal energy window for tetra.
3143 : !call find_ewin(gams%nqibz, gams%qibz, cryst, ebands, ltetra, sigma, comm)
3144 :
3145 : ! ========================
3146 : ! === MPI DISTRIBUTION ===
3147 : ! ========================
3148 : ! Init for sequential execution.
3149 5 : my_npert = natom3
3150 :
3151 30 : if (any(dtset%eph_np_pqbks /= 0)) then
3152 : ! Use parameters from input file.
3153 0 : pert_comm%nproc = dtset%eph_np_pqbks(1)
3154 0 : qpt_comm%nproc = dtset%eph_np_pqbks(2)
3155 0 : bsum_comm%nproc = dtset%eph_np_pqbks(3)
3156 0 : kpt_comm%nproc = dtset%eph_np_pqbks(4)
3157 0 : spin_comm%nproc = dtset%eph_np_pqbks(5)
3158 0 : my_npert = natom3 / pert_comm%nproc
3159 0 : ABI_CHECK(my_npert > 0, "pert_comm_nproc cannot be greater than 3 * natom.")
3160 0 : ABI_CHECK(mod(natom3, pert_comm%nproc) == 0, "pert_comm_nproc must divide 3 * natom.")
3161 : else
3162 : ! Automatic grid generation
3163 : ! Not the most efficient distribution if large number of MPI procs.
3164 : ! TODO: Add parallelism over perturbations although it's less efficient than the parallelism over k
3165 : ! It starts to be interesting if we implement symmetries in the k-integration though.
3166 :
3167 : ! TODO: Spin
3168 : ! Automatic grid generation over q-points and spins.
3169 : !if (new%nsppol == 2 .and. mod(nprocs, 2) == 0) then
3170 : ! spin_comm%nproc = 2
3171 : ! new%qpt_comm%nproc = nprocs / 2
3172 : !else
3173 : ! new%qpt_comm%nproc = nprocs
3174 : !end if
3175 :
3176 : ! Handle parallelism over perturbations first.
3177 : ! Use MPI communicator to distribute the 3 * natom perturbations to reduce memory requirements for DFPT potentials.
3178 : ! Ideally, perturbations are equally distributed --> total number of CPUs should be divisible by 3 * natom.
3179 : ! or at least, divisible by one integer i for i in [2, 3 * natom - 1].
3180 :
3181 5 : if (nsppol == 2 .and. mod(nproc, 2) == 0) then
3182 0 : spin_comm%nproc = 2
3183 0 : kpt_comm%nproc = nproc / 2
3184 : else
3185 : ! Try to have 3 perts per proc first because the q-point parallelism is more efficient.
3186 : ! The memory for W(R,r,ipert) will increase though.
3187 : !do cnt=natom,2,-1
3188 : ! if (mod(nprocs, cnt) == 0 .and. mod(natom3, cnt) == 0) then
3189 : ! pert_comm%nproc = cnt; new%my_npert = natom3 / cnt; exit
3190 : ! end if
3191 : !end do
3192 :
3193 : !if (pert_comm%nproc == 1) then
3194 : ! ! Try again with more procs.
3195 : ! do cnt=natom3,2,-1
3196 : ! if (mod(nprocs, cnt) == 0 .and. mod(natom3, cnt) == 0) then
3197 : ! pert_comm%nproc = cnt; new%my_npert = natom3 / cnt; exit
3198 : ! end if
3199 : ! end do
3200 : !end if
3201 :
3202 : !if (new%my_npert == natom3 .and. nprocs > 1) then
3203 : ! ABI_WARNING("The number of MPI procs should be divisible by 3*natom to reduce memory requirements!")
3204 : !end if
3205 :
3206 5 : kpt_comm%nproc = nproc
3207 : end if
3208 : end if
3209 :
3210 : ! Consistency check.
3211 5 : if (pert_comm%nproc * qpt_comm%nproc * bsum_comm%nproc * kpt_comm%nproc * spin_comm%nproc /= nproc) then
3212 : write(msg, "(a,i0,3a, 6(a,1x,i0))") &
3213 0 : "Cannot create 5d Cartesian grid with total nproc: ", nproc, ch10, &
3214 0 : "Idle processes are not supported. The product of the `nproc_*` vars should be equal to nproc.", ch10, &
3215 0 : "pert_nproc (", pert_comm%nproc, ") x qpt_nproc (", qpt_comm%nproc, ") x bsum_nproc (", bsum_comm%nproc, &
3216 0 : ") x kcalc_nproc (", kpt_comm%nproc, ") x spin_nproc (", spin_comm%nproc, ") != ", &
3217 0 : pert_comm%nproc * qpt_comm%nproc * bsum_comm%nproc * kpt_comm%nproc * spin_comm%nproc
3218 0 : ABI_ERROR(msg)
3219 : end if
3220 :
3221 5 : ABI_CHECK(bsum_comm%nproc == 1, "Band parallelism not implemented in m_phgamma")
3222 :
3223 : #ifdef HAVE_MPI
3224 : ! Create 5d cartesian communicator: 3*natom perturbations, q-points in IBZ, bands in FS, kpoints in FS, spins
3225 5 : ndims = 5
3226 5 : ABI_MALLOC(dims, (ndims))
3227 5 : ABI_MALLOC(periods, (ndims))
3228 5 : ABI_MALLOC(keepdim, (ndims))
3229 30 : periods(:) = .False.; reorder = .False.
3230 35 : dims = [pert_comm%nproc, qpt_comm%nproc, bsum_comm%nproc, kpt_comm%nproc, spin_comm%nproc]
3231 :
3232 5 : call MPI_CART_CREATE(comm, ndims, dims, periods, reorder, comm_cart, ierr)
3233 : ! Find the index and coordinates of the current processor
3234 5 : call MPI_COMM_RANK(comm_cart, me_cart, ierr)
3235 5 : call MPI_CART_COORDS(comm_cart, me_cart, ndims, coords, ierr)
3236 :
3237 : ! Create communicator to distribute natom3 perturbations.
3238 30 : keepdim = .False.; keepdim(1) = .True.
3239 5 : call MPI_CART_SUB(comm_cart, keepdim, pert_comm%value, ierr); pert_comm%me = xmpi_comm_rank(pert_comm%value)
3240 : ! Create communicator for qpoints in self-energy integration.
3241 30 : keepdim = .False.; keepdim(2) = .True.
3242 5 : call MPI_CART_SUB(comm_cart, keepdim, qpt_comm%value, ierr); qpt_comm%me = xmpi_comm_rank(qpt_comm%value)
3243 : ! Create communicator for bands for band summation
3244 30 : keepdim = .False.; keepdim(3) = .True.
3245 5 : call MPI_CART_SUB(comm_cart, keepdim, bsum_comm%value, ierr); bsum_comm%me = xmpi_comm_rank(bsum_comm%value)
3246 : ! Create communicator for kpoints.
3247 30 : keepdim = .False.; keepdim(4) = .True.
3248 5 : call MPI_CART_SUB(comm_cart, keepdim, kpt_comm%value, ierr); kpt_comm%me = xmpi_comm_rank(kpt_comm%value)
3249 : ! Create communicator for spins.
3250 30 : keepdim = .False.; keepdim(5) = .True.
3251 5 : call MPI_CART_SUB(comm_cart, keepdim, spin_comm%value, ierr); spin_comm%me = xmpi_comm_rank(spin_comm%value)
3252 : ! Create communicator for the (qpoint, spin) loops
3253 30 : keepdim = .False.; keepdim(2) = .True.; keepdim(5) = .True.
3254 5 : call MPI_CART_SUB(comm_cart, keepdim, qs_comm%value, ierr); qs_comm%me = xmpi_comm_rank(qs_comm%value)
3255 : ! Create communicator for the (perturbation, k-point, band_sum)
3256 40 : keepdim = .False.; keepdim(1) = .True.; keepdim(3:4) = .True.
3257 5 : call MPI_CART_SUB(comm_cart, keepdim, pkb_comm%value, ierr); pkb_comm%me = xmpi_comm_rank(pkb_comm%value)
3258 :
3259 5 : ABI_FREE(dims)
3260 5 : ABI_FREE(periods)
3261 5 : ABI_FREE(keepdim)
3262 5 : call xmpi_comm_free(comm_cart)
3263 : #endif
3264 :
3265 : ! Distribute spins and q-points in IBZ.
3266 5 : call xmpi_split_cyclic(nsppol, spin_comm%value, gams%my_nspins, gams%my_spins)
3267 5 : ABI_CHECK(gams%my_nspins > 0, sjoin("nsppol (", itoa(nsppol), ") < spin_comm_nproc (", itoa(spin_comm%nproc), ")"))
3268 :
3269 5 : call xmpi_split_cyclic(gams%nqibz, qpt_comm%value, gams%my_nqibz, gams%my_iqibz)
3270 5 : ABI_CHECK(gams%my_nqibz > 0, sjoin("nqibz (", itoa(gams%nqibz), ") < qpt_comm_nproc (", itoa(qpt_comm%nproc), ")"))
3271 :
3272 : !ABI_CALLOC(qibz_done, (gams%nqibz))
3273 5 : path = strcat(dtfil%filnam_ds(4), "_A2F.nc")
3274 5 : ncid = nctk_noid
3275 :
3276 5 : if (my_rank == master) then
3277 5 : call cwtime(cpu, wall, gflops, "start")
3278 :
3279 5 : write(std_out, "(/,a)")" === MPI parallelism ==="
3280 : !write(std_out, "(2(a,i0))")"P Allocating and summing bands from my_bsum_start: ", self%my_bsum_start, &
3281 : ! " up to my_bsum_stop: ", self%my_bsum_stop
3282 5 : write(std_out, "(a,i0)")"P Number of CPUs for parallelism over perturbations: ", pert_comm%nproc
3283 5 : write(std_out, "(a,i0)")"P Number of perturbations treated by this CPU: ", my_npert
3284 5 : write(std_out, "(a,i0)")"P Number of CPUs for parallelism over q-points: ", qpt_comm%nproc
3285 : !write(std_out, "(2(a,i0))")"P Number of q-points in the IBZ treated by this proc: " , &
3286 : ! count(self%itreat_qibz == 1), " of ", self%nqibz
3287 5 : write(std_out, "(a,i0)")"P Number of CPUs for parallelism over bands: ", bsum_comm%nproc
3288 5 : write(std_out, "(a,i0)")"P Number of CPUs for parallelism over spins: ", spin_comm%nproc
3289 5 : write(std_out, "(a,i0)")"P Number of CPUs for parallelism over k-points: ", kpt_comm%nproc
3290 : !write(std_out, "(2(a,i0),/)")"P Number of k-point in Sigma_nk treated by this proc: ", self%my_nkcalc, " of ", self%nkcalc
3291 :
3292 : ! Master creates the netcdf file used to store the results of the calculation.
3293 5 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
3294 5 : NCF_CHECK(cryst%ncwrite(ncid))
3295 5 : NCF_CHECK(ebands%ncwrite(ncid))
3296 5 : NCF_CHECK(edos%ncwrite(ncid))
3297 :
3298 : ! Add dimensions.
3299 15 : ncerr = nctk_def_dims(ncid, [nctkdim_t("nqibz", gams%nqibz), nctkdim_t("natom3", natom3)], defmode=.True.)
3300 5 : NCF_CHECK(ncerr)
3301 20 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat"])
3302 5 : NCF_CHECK(ncerr)
3303 10 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "ph_intmeth"])
3304 5 : NCF_CHECK(ncerr)
3305 25 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie"])
3306 5 : NCF_CHECK(ncerr)
3307 15 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ph_wstep", "ph_smear"])
3308 5 : NCF_CHECK(ncerr)
3309 :
3310 : ! Define arrays for results.
3311 : ncerr = nctk_def_arrays(ncid, [ &
3312 : nctkarr_t("ngqpt", "int", "three"), &
3313 : nctkarr_t("eph_ngqpt_fine", "int", "three"), &
3314 : nctkarr_t("ddb_ngqpt", "int", "three"), &
3315 : nctkarr_t("ph_ngqpt", "int", "three"), &
3316 : nctkarr_t('qibz', "dp", "number_of_reduced_dimensions, nqibz"), &
3317 : nctkarr_t('wtq', "dp", "nqibz"), &
3318 : !nctkarr_t('qibz_done', "int", "nqibz"), &
3319 : nctkarr_t('phfreq_qibz', "dp", "natom3, nqibz"), &
3320 : nctkarr_t('phdispl_cart_qibz', "dp", "two, natom3, natom3, nqibz"), &
3321 : nctkarr_t('phgamma_qibz', "dp", "natom3, nqibz, number_of_spins"), &
3322 : nctkarr_t('phlambda_qibz', "dp", "natom3, nqibz, number_of_spins") &
3323 55 : ])
3324 5 : NCF_CHECK(ncerr)
3325 :
3326 : ! ======================================================
3327 : ! Write data that do not depend on the (kpt, spin) loop.
3328 : ! ======================================================
3329 5 : NCF_CHECK(nctk_set_datamode(ncid))
3330 :
3331 : ncerr = nctk_write_iscalars(ncid, &
3332 : [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat", "ph_intmeth"], &
3333 45 : [dtset%eph_intmeth, dtset%eph_transport, dtset%symdynmat, dtset%ph_intmeth])
3334 5 : NCF_CHECK(ncerr)
3335 : ncerr = nctk_write_dpscalars(ncid, &
3336 : [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear"], &
3337 65 : [dtset%eph_fsewin, dtset%eph_fsmear, dtset%eph_extrael, dtset%eph_fermie, dtset%ph_wstep, dtset%ph_smear])
3338 5 : NCF_CHECK(ncerr)
3339 :
3340 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qibz'), gams%qibz))
3341 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'wtq'), gams%wtq))
3342 : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qibz_done'), qibz_done))
3343 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), gamma_ngqpt))
3344 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine))
3345 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ph_ngqpt"), dtset%ph_ngqpt))
3346 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt))
3347 : !NCF_CHECK(nf90_close(ncid))
3348 5 : call cwtime_report(" Creation of A2F.nc", cpu, wall, gflops, end_str=ch10)
3349 : end if
3350 : !call xmpi_barrier(comm)
3351 :
3352 : ! Now reopen the file inside ncwrite_comm to perform parallel-IO (required for q-point parallelism).
3353 : !if (self%ncwrite_comm%value /= xmpi_comm_null) then
3354 : ! NCF_CHECK(nctk_open_modify(self%ncid, path, self%ncwrite_comm%value))
3355 : ! NCF_CHECK(nctk_set_datamode(self%ncid))
3356 : !end if
3357 5 : call edos%free()
3358 :
3359 : ! Open the DVDB file
3360 5 : call dvdb%open_read(ngfftf, xmpi_comm_self)
3361 5 : ABI_CHECK(dvdb%has_fields("pot1", msg), msg)
3362 :
3363 5 : if (pert_comm%nproc > 1) then
3364 : ! Activate parallelism over perturbations
3365 : ! Build table with list of perturbations treated by this CPU inside pert_comm
3366 0 : call ephtk_set_pertables(cryst%natom, my_npert, pert_table, my_pinfo, pert_comm%value)
3367 0 : call dvdb%set_pert_distrib(my_npert, natom3, my_pinfo, pert_table, pert_comm%value)
3368 0 : ABI_FREE(my_pinfo)
3369 0 : ABI_FREE(pert_table)
3370 : end if
3371 :
3372 : ! Check whether all q-points are available in the input DVDB file.
3373 5 : need_ftinterp = .False.
3374 66 : do iq_ibz=1,gams%nqibz
3375 66 : if (dvdb%findq(gams%qibz(:,iq_ibz) ) == -1) then
3376 0 : need_ftinterp = .True.; exit
3377 : end if
3378 : end do
3379 5 : call wrtout(std_out, " ", do_flush=.True.)
3380 :
3381 : ! Find correspondence IBZ --> set of q-points in DVDB.
3382 : ! Activate FT interpolation automatically if required q-points in the IBZ are not found in the DVDB.
3383 : !need_ftinterp = .False.
3384 : !ABI_MALLOC(qibz2dvdb, (sigma%nqibz))
3385 : !if (dvdb%find_qpts(sigma%nqibz, sigma%qibz, qibz2dvdb, comm) /= 0) then
3386 : ! call wrtout(units, " Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.")
3387 : ! need_ftinterp = .True.
3388 : !else
3389 : ! call wrtout(units, " DVDB file contains all q-points in the IBZ --> Reading DFPT potentials from file.")
3390 : ! need_ftinterp = .False.
3391 : !end if
3392 :
3393 : !if (.not. need_ftinterp .and. dtset%eph_use_ftinterp /= 0) then
3394 : ! ABI_WARNING("Enforcing FT interpolation for q-points even if it's not strictly needed.")
3395 : ! need_ftinterp = .True.
3396 : !end if
3397 :
3398 5 : if (need_ftinterp) then
3399 0 : call wrtout(units, " Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.")
3400 : ! Prepare Fourier interpolation of DFPT potentials.
3401 0 : comm_rpt = xmpi_comm_self
3402 : !comm_rpt = bqs_comm%value
3403 0 : call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt)
3404 : end if
3405 :
3406 : ! Initialize the wave function descriptor.
3407 : ! Only wavefunctions on the FS are stored in wfd.
3408 : ! Need all k-points on the FS because of k+q, spin is not distributed for the time being.
3409 : ! It would be possible to reduce the memory allocated per MPI-rank via OpenMP.
3410 :
3411 20 : ABI_MALLOC(nband, (nkibz, nsppol))
3412 25 : ABI_MALLOC(bks_mask, (mband, nkibz, nsppol))
3413 20 : ABI_MALLOC(keep_ur, (mband, nkibz ,nsppol))
3414 1905 : nband = mband; bks_mask = .False.; keep_ur = .False.
3415 :
3416 10 : do my_is=1,gams%my_nspins
3417 5 : spin = gams%my_spins(my_is); fs => fstab(spin)
3418 700 : do ik_bz=1,fs%nkfs
3419 690 : ik_ibz = fs%indkk_fs(1, ik_bz)
3420 690 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nb_k = fs%bstart_cnt_ibz(2, ik_ibz)
3421 1385 : bks_mask(bstart_k:bstart_k+nb_k-1, ik_ibz, spin) = .True.
3422 : end do
3423 : end do
3424 : !bks_mask(1:mband,:,:) = .True. ! no memory distribution, each node has the full set of states.
3425 :
3426 : ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts)
3427 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
3428 15 : ABI_MALLOC(wfd_istwfk, (nkibz))
3429 150 : wfd_istwfk = 1
3430 :
3431 5 : ecut = dtset%ecut ! dtset%dilatmx
3432 : call wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkibz, nsppol, bks_mask,&
3433 : nspden, nspinor, ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,&
3434 5 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
3435 :
3436 10 : call wfd%print([std_out], header="Wavefunctions on the Fermi surface")
3437 :
3438 5 : ABI_FREE(nband)
3439 5 : ABI_FREE(keep_ur)
3440 5 : ABI_FREE(wfd_istwfk)
3441 :
3442 : ! Read wavefunctions.
3443 5 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
3444 :
3445 : ! one-dimensional structure factor information on the coarse grid.
3446 15 : ABI_MALLOC(ph1d, (2, 3*(2*mgfft+1)*natom))
3447 5 : call getph(cryst%atindx, natom, n1, n2, n3, ph1d, cryst%xred)
3448 :
3449 : ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ.
3450 : ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work
3451 : ! that will be used to symmetrize the wavefunctions in G-space.
3452 5 : call cwtime(cpu, wall, gflops, "start")
3453 5 : call wrtout(std_out, " Computing mpw. This may take some time for dense k/q meshes...")
3454 : !call fstab_get_mpw_gmax(nsppol, fstab, mpw, gmax, comm)
3455 5 : mpw = 0; gmax = 0; cnt = 0
3456 :
3457 : ! FIXME: This is an hotspot due to the loop over nkfs.
3458 : ! Should use a geometric approach to compute mpw gmax.
3459 10 : do spin=1,nsppol
3460 5 : fs => fstab(spin)
3461 700 : do ik_bz=1,fs%nkfs
3462 690 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism inside comm
3463 2760 : kk = fs%kpts(:, ik_bz)
3464 : ! Compute G sphere, returning npw. Note istwfk == 1.
3465 690 : call get_kg(kk, 1, ecut, cryst%gmet, onpw, gtmp)
3466 690 : mpw = max(mpw, onpw)
3467 52350 : do ipw=1,onpw
3468 207330 : do ii=1,3
3469 206640 : gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw)))
3470 : end do
3471 : end do
3472 690 : ABI_FREE(gtmp)
3473 :
3474 9113 : do iq_ibz=1,gams%nqibz
3475 33672 : qpt = gams%qibz(:,iq_ibz)
3476 : ! TODO: g0 umklapp here can enter into play!
3477 : ! fstab should contains the max of the umlapp G-vectors.
3478 : ! gmax could not be large enough!
3479 33672 : kq = kk + qpt
3480 8418 : call get_kg(kq, 1, ecut, cryst%gmet, onpw, gtmp)
3481 8418 : mpw = max(mpw, onpw)
3482 633595 : do ipw=1,onpw
3483 2509126 : do ii=1,3
3484 2500708 : gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw)))
3485 : end do
3486 : end do
3487 9108 : ABI_FREE(gtmp)
3488 : end do
3489 :
3490 : end do ! ik_bz
3491 : end do ! spin
3492 :
3493 5 : my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr)
3494 5 : my_gmax = gmax; call xmpi_max(my_gmax, gmax, comm, ierr)
3495 5 : call wrtout(std_out, sjoin(' Optimal value of mpw: ', itoa(mpw)))
3496 5 : call cwtime_report(" Computation of gmax and mpw:", cpu, wall, gflops)
3497 :
3498 : ! Init work_ngfft
3499 20 : gmax = gmax + 4 ! FIXME: this is to account for umklapp
3500 20 : gmax = 2*gmax + 1
3501 5 : call ngfft_seq(work_ngfft, gmax)
3502 : !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
3503 25 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
3504 :
3505 : ! Allow PW-arrays dimensioned with mpw
3506 15 : ABI_MALLOC(kg_k, (3, mpw))
3507 10 : ABI_MALLOC(kg_kq, (3, mpw))
3508 :
3509 : ! TODO FOR PAW
3510 5 : usecprj = 0
3511 10 : ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj))
3512 :
3513 : ! Prepare call to getgh1c
3514 5 : usevnl = 0
3515 5 : optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C>
3516 5 : optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
3517 5 : opt_gvnlx1 = 0 ! gvnlx1 is output
3518 :
3519 5 : ABI_MALLOC(grad_berry, (2, nspinor*(berryopt0/4)))
3520 :
3521 : ! This part is taken from dfpt_vtorho
3522 : !==== Initialize most of the Hamiltonian (and derivative) ====
3523 : !1) Allocate all arrays and initialize quantities that do not depend on k and spin.
3524 : !2) Perform the setup needed for the non-local factors:
3525 : ! Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk.
3526 : ! PAW: Initialize the overlap coefficients and allocate the Dij coefficients.
3527 :
3528 : call gs_hamkq%init(psps, pawtab, nspinor, nsppol, nspden, natom, &
3529 : dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg, &
3530 : comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, &
3531 5 : usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option)
3532 :
3533 : ! Allocate vlocal. Note nvloc
3534 : ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data)
3535 30 : ABI_MALLOC(vlocal, (n4, n5, n6, gs_hamkq%nvloc))
3536 10990 : vlocal = huge(one)
3537 :
3538 : ! Allocate work space arrays.
3539 20 : ABI_MALLOC(tgam, (2, natom3, natom3))
3540 20 : ABI_MALLOC(displ_cart, (2, 3, cryst%natom, natom3))
3541 15 : ABI_MALLOC(displ_red, (2, 3, cryst%natom, natom3))
3542 8665 : ABI_CALLOC(dummy_vtrial, (nfftf, nspden))
3543 :
3544 : ! TODO: Save data to netcdf file for each q in IBZ and then read data to build a2Fw once all big
3545 : ! datastructures (wfd, dvdb) have been deallocated.
3546 : ! As a side effect, one can also implement restart over q-points
3547 :
3548 : ! Create ddkop object to compute group velocities if needed.
3549 : !
3550 : ! 1) precompute group velocities in the IBZ and the ihave_ikibz_spin file (common to all procs)
3551 : ! 2) Use symmetries to reconstruct v_kq from vcart_ibz
3552 : !
3553 : ! NB: All procs store in memory the same set of Bloch states inside the energy window.
3554 :
3555 5 : call ddkop%init(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft)
3556 :
3557 5 : call cwtime(cpu, wall, gflops, "start", msg=" Computing v_nk matrix elements for all states on the FS...")
3558 5 : ii = huge(1); jj = -1
3559 10 : do spin=1,nsppol
3560 5 : ii = min(ii, fstab(spin)%bmin)
3561 10 : jj = max(jj, fstab(spin)%bmax)
3562 : end do
3563 1335 : ABI_CALLOC(vcart_ibz, (3, ii:jj, nkibz, nsppol))
3564 15 : ABI_MALLOC(cg_work, (2, mpw * wfd%nspinor))
3565 :
3566 5 : cnt = 0
3567 10 : do spin=1,nsppol
3568 5 : fs => fstab(spin)
3569 155 : do ik_ibz=1,ebands%nkpt
3570 580 : kk = ebands%kptns(:, ik_ibz); npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz)
3571 : ! NB: The two checks below are global --> all procs will cycle.
3572 760 : if (all(bks_mask(:, ik_ibz, spin) .eqv. .False.)) cycle
3573 30 : if (npw_k == 1) cycle
3574 30 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
3575 :
3576 30 : call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk, istwf_k, npw_k, wfd%kdata(ik_ibz)%kg_k)
3577 :
3578 95 : do band_k=fs%bmin,fs%bmax
3579 60 : if (.not. bks_mask(band_k, ik_ibz, spin)) cycle
3580 30 : call wfd%copy_cg(band_k, ik_ibz, spin, cg_work)
3581 30 : eig0nk = ebands%eig(band_k, ik_ibz, spin)
3582 30 : vk = ddkop%get_vdiag(eig0nk, istwf_k, npw_k, wfd%nspinor, cg_work, cwaveprj0)
3583 265 : vcart_ibz(:, band_k, ik_ibz, spin) = vk
3584 :
3585 : ! TODO: Use ebands_get_edos_matrix_elements
3586 : ! reald(dp) :: vv_fs(3,3,nsppol)
3587 : !vv_fs = zero
3588 : !sigma = fs%eph_fsmear
3589 : !if (fs%eph_fsmear < zero) then
3590 : ! sigma = max(maxval([(abs(dot_product(fs%vk(:, ib2), fs%kmesh_cartvec(:,ii))), ii=1,3)]), fs%min_smear)
3591 : !end if
3592 : !fs_wtk = gaussian(emesh - ebands%eig(band_k, ik_ibz, spin), sigma) / (one * fs%nktot)
3593 : !do ii=1,3
3594 : ! do jj=1,3
3595 : ! vv_fs(ii,jj,spin) = vv_fs(ii,jj,spin) + vk(ii) * vk(jj) * fs_wtk
3596 : ! end do
3597 : !end do
3598 : !call xmpi_sum(vv_fs, comm, ierr)
3599 : end do
3600 : end do
3601 : end do ! spin
3602 :
3603 5 : call xmpi_sum(vcart_ibz, comm, ierr)
3604 5 : ABI_FREE(cg_work)
3605 5 : call cwtime_report(" Computation of group velocities", cpu, wall, gflops)
3606 :
3607 5 : ABI_FREE(bks_mask)
3608 :
3609 : !Write v_nk to disk.
3610 : !if (my_rank == master) then
3611 : ! NCF_CHECK(nf90_put_var(sigma%ncid, nctk_idname(sigma%ncid, "vcart_ibz"), vcart_ibz))
3612 : !end if
3613 :
3614 5 : if (dtset%eph_transport > 0) then
3615 20 : ABI_MALLOC(tgamvv_in, (2, 9, natom3, natom3))
3616 15 : ABI_MALLOC(tgamvv_out, (2, 9, natom3, natom3))
3617 : end if
3618 :
3619 : ! Build krank object to find k-points
3620 5 : call krank%from_kptrlatt(ebands%nkpt, ebands%kptns, ebands%kptrlatt, compute_invrank=.False.)
3621 5 : call pstat_proc%print(_PSTAT_ARGS_)
3622 :
3623 : ! Loop over my q-points in the IBZ.
3624 66 : do my_iq=1,gams%my_nqibz
3625 61 : iq_ibz = gams%my_iqibz(my_iq)
3626 :
3627 : ! TODO: Restart capabilities
3628 : ! Check if this (kpoint, spin) was already calculated
3629 : !if (all(qibz_done(iq_ibz, :) == 1)) cycle
3630 :
3631 61 : call cwtime(cpu_q, wall_q, gflops_q, "start")
3632 :
3633 244 : qpt = gams%qibz(:, iq_ibz)
3634 61 : msg = sjoin("[", itoa(iq_ibz), "/", itoa(gams%nqibz), "]")
3635 61 : call wrtout(std_out, sjoin(" Computing ph linewidths for IBZ q-point:", ktoa(qpt), msg))
3636 :
3637 1891 : tgam = zero
3638 61 : if (dtset%eph_transport > 0) then
3639 31171 : tgamvv_in = zero; tgamvv_out = zero
3640 : end if
3641 :
3642 61 : if (.not. need_ftinterp) then
3643 : ! No interpolation --> find the index of the q-point in the DVDB.
3644 61 : db_iqpt = dvdb%findq(qpt)
3645 :
3646 61 : if (db_iqpt /= -1) then
3647 : ! Read and reconstruct the dvscf potentials for all 3*natom perturbations.
3648 : ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom))
3649 61 : call dvdb%readsym_allv1(db_iqpt, cplex, nfftf, ngfftf, v1scf, pkb_comm%value)
3650 :
3651 61 : if (dvdb%my_npert /= natom3) then
3652 : ! Extract my npert from v1scf
3653 0 : ABI_MALLOC(v1_work, (cplex, nfftf, nspden, dvdb%my_npert))
3654 0 : do my_ip=1,dvdb%my_npert
3655 0 : v1_work(:,:,:,my_ip) = v1scf(:,:,:,dvdb%my_pinfo(3, my_ip))
3656 : end do
3657 0 : ABI_MOVE_ALLOC(v1_work, v1scf)
3658 : end if
3659 :
3660 : else
3661 0 : ABI_ERROR(sjoin("Cannot find q-point:", ktoa(qpt), "in the DVDB file."))
3662 : end if
3663 :
3664 : else
3665 : ! Use Fourier interpolation of DFPT potentials to get my_npert potentials.
3666 0 : cplex = 2
3667 0 : ABI_MALLOC(v1scf, (cplex, nfft, nspden, dvdb%my_npert))
3668 0 : call dvdb%ftinterp_qpt(qpt, nfftf, ngfftf, v1scf, dvdb%comm_rpt)
3669 : end if
3670 :
3671 : ! Compute the little group of the q-point so that we can compute g(k,q) only for k in the IBZ_q
3672 : !if (dtset%symsigma /= 0) then
3673 : ! timrev_q = kpts_timrev_from_kptopt(my_qptopt)
3674 : ! call lg_q%init(cryst, qpt, timrev_q, gstore%nkbz, gstore%kbz, gstore%nkibz, gstore%kibz, xmpi_comm_self)
3675 : !end if
3676 :
3677 : !if (dtset%symsigma /= 0) then
3678 : ! ii = lg_q%findq_ibzk(kk)
3679 : ! if (ii == -1) then
3680 : ! call wrtout(std_out, sjoin(" my_ik:", itoa(my_ik), kk_string, " not in IBZ_q --> skipping iteration"))
3681 : ! cycle
3682 : ! ! TODO: Check fillvalue (should be zero)
3683 : ! end if
3684 : ! weight_k = lg_q%weights(ii)
3685 : !end if
3686 : !call lg_q%free()
3687 :
3688 : ! Get phonon frequencies and eigenvectors for this q-point.
3689 61 : call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red)
3690 61 : wminmax(1) = min(wminmax(1), phfrq(1))
3691 61 : wminmax(2) = max(wminmax(2), phfrq(3*cryst%natom))
3692 :
3693 : ! Allocate vlocal1 with correct cplex. Note nvloc and my_npert.
3694 427 : ABI_MALLOC(vlocal1, (cplex*n4, n5, n6, gs_hamkq%nvloc, my_npert))
3695 :
3696 : ! Loop over my spins.
3697 122 : do my_is=1,gams%my_nspins
3698 61 : spin = gams%my_spins(my_is); fs => fstab(spin)
3699 :
3700 61 : if (dtset%prteliash == 3) then
3701 0 : ABI_CALLOC_OR_DIE(tmp_vals_ee, (2, gams%nene, gams%nene, natom3, natom3), ierr)
3702 : ! Energy mesh for electrons in a2F(e,e',w) NB: It depends on the spin through enemin
3703 0 : ABI_MALLOC(emesh, (gams%nene))
3704 0 : emesh = arth(fs%enemin, fs%deltaene, gams%nene)
3705 : end if
3706 :
3707 : ! Set up local potential vlocal1 with proper dimensioning from vtrial1 taking into account the spin.
3708 244 : do my_ip=1,my_npert
3709 : call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, gs_hamkq%nvloc,&
3710 244 : pawfgr, mpi_enreg, dummy_vtrial, v1scf(:,:,:,my_ip), vlocal, vlocal1(:,:,:,:,my_ip))
3711 : end do
3712 :
3713 : ! Continue to initialize the GS Hamiltonian
3714 61 : call gs_hamkq%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.)
3715 :
3716 : ! Allocate workspace for wavefunctions. Make npw larger than expected.
3717 : ! maxnb is the maximum number of bands crossing the FS, used to dimension arrays.
3718 61 : mnb = fs%maxnb
3719 : ! The weights for the integration of the double-delta.
3720 244 : ABI_MALLOC(dbldelta_wts, (mnb, mnb))
3721 :
3722 61 : if (dtset%eph_transport > 0) then
3723 244 : ABI_MALLOC(vv_kk, (9, mnb, mnb))
3724 183 : ABI_MALLOC(vv_kkq, (9, mnb, mnb))
3725 : end if
3726 :
3727 61 : if (dtset%prteliash == 3) then
3728 0 : ABI_MALLOC(wt_ek, (gams%nene, mnb))
3729 0 : ABI_MALLOC(wt_ekq, (gams%nene, mnb))
3730 : end if
3731 :
3732 : ! =====================================
3733 : ! Integration over the FS for this spin
3734 : ! =====================================
3735 : ! Compute integration weights and distribute k-points (gams%my_nfsk_q) inside kpt_comm
3736 61 : call phgamma_setup_qpoint(gams, fs, cryst, ebands, spin, ltetra, qpt, nesting, kpt_comm%value)
3737 :
3738 6265 : do my_ik=1,gams%my_nfsk_q
3739 6204 : print_time_k = (my_ik <= 2 .or. (fs%nkfs > 100 .and. mod(my_ik, 200) == 0))
3740 122 : if (print_time_k) call cwtime(cpu_k, wall_k, gflops_k, "start")
3741 :
3742 : ! The k-point and the symmetries relating the BZ k-point to the IBZ.
3743 6204 : ik_fs = gams%my_ifsk_q(my_ik)
3744 24816 : kk = fs%kpts(:, ik_fs)
3745 6204 : ik_ibz = fs%indkk_fs(1, ik_fs); isym_k = fs%indkk_fs(2, ik_fs)
3746 24816 : trev_k = fs%indkk_fs(6, ik_fs); g0_k = fs%indkk_fs(3:5,ik_fs)
3747 19906 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
3748 24816 : kk_ibz = ebands%kptns(:,ik_ibz)
3749 :
3750 : ! Number of bands crossing the Fermi level at k.
3751 6204 : bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nb_k = fs%bstart_cnt_ibz(2, ik_ibz)
3752 :
3753 : ! Find k+q in the extended zone and extract symmetry info. cycle if k+q not in FS.
3754 : ! Be careful here because there are two umklapp vectors to be considered:
3755 : !
3756 : ! k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz
3757 : !
3758 24816 : kq = kk + qpt; ikq_bz = fs%findkg0(kq, g0bz_kq)
3759 :
3760 : ! Skip this point if kq does not belong to the FS window.
3761 6204 : if (ikq_bz == -1) cycle
3762 :
3763 2864 : if (kpts_map("symrel", ebands%kptopt, cryst, krank, 1, kq, indkk_kq) /= 0) then
3764 : write(msg, '(9a)' ) &
3765 0 : "The WFK file cannot be used to compute phonon linewidths.",ch10, &
3766 0 : "At least one of the k-points on the FS could not be generated from a symmetrical one.", ch10, &
3767 0 : "q-mesh: ", trim(ltoa(gamma_ngqpt)), ", k-mesh (from kptrlatt) ", trim(ltoa(get_diag(ebands%kptrlatt))), &
3768 0 : 'Action: check your WFK file and the (k, q) point input variables.'
3769 0 : ABI_ERROR(msg)
3770 : end if
3771 :
3772 : ! The k+q point and the symmetries relating it to the IBZ.
3773 2864 : ikq_ibz = indkk_kq(1, 1); isym_kq = indkk_kq(2, 1)
3774 11456 : trev_kq = indkk_kq(6, 1); g0_kq = indkk_kq(3:5, 1)
3775 6934 : isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0))
3776 11456 : kq_ibz = ebands%kptns(:, ikq_ibz)
3777 :
3778 : ! If we have used the KERANGE trick, we may have k or k+q points with just one G component set to zero
3779 : ! so we skip this transition immediately. This should happen only if fsewin > sigma_erange.
3780 2864 : if (wfd%npwarr(ik_ibz) == 1 .or. wfd%npwarr(ikq_ibz) == 1) cycle
3781 :
3782 : ! Number of bands crossing the Fermi level at k+q.
3783 2864 : bstart_kq = fs%bstart_cnt_ibz(1, ikq_ibz); nb_kq = fs%bstart_cnt_ibz(2, ikq_ibz)
3784 2864 : ABI_CHECK(nb_k <= mnb .and. nb_kq <= mnb, "wrong nband")
3785 :
3786 : ! Get npw_k, kg_k and symmetrize wavefunctions from IBZ (if needed).
3787 : call wfd%sym_ug_kg_npw(ecut, kk, kk_ibz, bstart_k, nb_k, spin, fs%indkk_fs(:,ik_fs), cryst, &
3788 2864 : work_ngfft, work, istwf_k, npw_k, kg_k, kets_k)
3789 :
3790 : ! Get npw_kq, kg_kq and symmetrize wavefunctions from IBZ (if needed).
3791 : call wfd%sym_ug_kg_npw(ecut, kq, kq_ibz, bstart_kq, nb_kq, spin, indkk_kq(:,1), cryst, &
3792 2864 : work_ngfft, work, istwf_kq, npw_kq, kg_kq, bras_kq)
3793 :
3794 : call gs_hamkq%eph_setup_k("k" , kk, istwf_k, npw_k, kg_k, dtset, cryst, psps, & ! in
3795 2864 : nkpg, kpg_k, ffnl_k, kinpw_k, ph3d_k, pert_comm%value) ! out
3796 :
3797 : call gs_hamkq%eph_setup_k("kq", kq, istwf_k, npw_kq, kg_kq, dtset, cryst, psps, & ! in
3798 2864 : nkpg, kpg_kq, ffnl_kq, kinpw_kq, ph3d_kq, pert_comm%value) ! out
3799 :
3800 : ! If PAW, one has to solve a generalized eigenproblem
3801 2864 : gen_eigenpb = psps%usepaw == 1; sij_opt = 0; if (gen_eigenpb) sij_opt = 1
3802 :
3803 8592 : ABI_MALLOC(lambda, (nb_k))
3804 14320 : ABI_MALLOC(gkq_atm, (2, nb_kq, nb_k, natom3))
3805 11456 : ABI_MALLOC(h1_kets_kq, (2, npw_kq*nspinor, nb_k))
3806 8592 : ABI_MALLOC(gvnlx1, (2, npw_kq*nspinor, nb_k))
3807 8592 : ABI_MALLOC(gs1c_kq, (2, npw_kq*nspinor*nb_k*((sij_opt+1)/2)))
3808 11456 : ABI_MALLOC(gkq_atm_ipc, (2, nb_kq, nb_k))
3809 : #ifdef HAVE_OPENMP_OFFLOAD
3810 : !$OMP TARGET ENTER DATA MAP(alloc:gkq_atm_ipc) IF (dtset%gpu_option == ABI_GPU_OPENMP)
3811 : !$OMP TARGET ENTER DATA MAP(to:kets_k, bras_kq) IF (dtset%gpu_option == ABI_GPU_OPENMP)
3812 : !$OMP TARGET ENTER DATA MAP(alloc:h1_kets_kq, gvnlx1, kets_k, bras_kq) IF (dtset%gpu_option == ABI_GPU_OPENMP)
3813 : !$OMP TARGET ENTER DATA MAP(alloc:gs1c_kq) IF (dtset%gpu_option == ABI_GPU_OPENMP .and. sij_opt /= 0)
3814 : #endif
3815 :
3816 : ! Loop over all my atomic perturbations and compute gkq_atm.
3817 11456 : do my_ip=1,my_npert
3818 8592 : idir = dvdb%my_pinfo(1, my_ip); ipert = dvdb%my_pinfo(2, my_ip); ipc = dvdb%my_pinfo(3, my_ip)
3819 :
3820 : ! Prepare application of the NL part.
3821 8592 : call rf_hamkq%init(cplex, gs_hamkq, ipert, has_e1kbsc=.true.)
3822 8592 : call rf_hamkq%load_spin(spin, vlocal1=vlocal1(:,:,:,:,my_ip), with_nonlocal=.true.)
3823 :
3824 : #if 1
3825 : ! Calculate dvscf * psi_k, results stored in h1_kets_kq on the k+q sphere.
3826 : ! Compute H(1) applied to GS wavefunction Psi(0)
3827 17184 : do in_k=1,nb_k
3828 : ! Use scissor shift on 0-order eigenvalue
3829 8592 : band_k = in_k + bstart_k - 1
3830 17184 : lambda(in_k) = ebands%eig(band_k, ik_ibz, spin) - dtset%dfpt_sciss
3831 : end do
3832 : !call ebands%get_dfpt_eshifted(gqk%bstart_k, nb_k, ik_ibz, spin, dtset%dfpt_sciss, lambda)
3833 :
3834 : call getgh1c(berryopt0, kets_k, cwaveprj0, h1_kets_kq, &
3835 : grad_berry, gs1c_kq, gs_hamkq, gvnlx1, idir, ipert, lambda, mpi_enreg, nb_k, optlocal, &
3836 8592 : optnl, opt_gvnlx1, rf_hamkq, sij_opt, tim_getgh1c, usevnl)
3837 :
3838 : #else
3839 : do in_k=1,nb_k
3840 : ! Use scissor shift on 0-order eigenvalue
3841 : band_k = in_k + bstart_k - 1
3842 : eshift = ebands%eig(band_k, ik_ibz, spin) - dtset%dfpt_sciss
3843 :
3844 : call getgh1c(berryopt0, kets_k(:,:,in_k), cwaveprj0, h1_kets_kq(:,:,in_k), &
3845 : grad_berry, gs1c_kq, gs_hamkq, gvnlx1, idir, ipert, [eshift], mpi_enreg, ndat1, optlocal, &
3846 : optnl, opt_gvnlx1, rf_hamkq, sij_opt, tim_getgh1c, usevnl)
3847 : end do ! in_k
3848 : #endif
3849 :
3850 8592 : call rf_hamkq%free()
3851 :
3852 : ! Calculate <psi_{k+q,m}|dvscf_q|psi_{k,n}> for this perturbation.
3853 : ! No need to handle istwf_kq because it's always 1.
3854 8592 : if (dtset%gpu_option == ABI_GPU_OPENMP) then
3855 : call abi_gpu_xgemm_d(2, 'C', 'N', nb_kq, nb_k, npw_kq*nspinor, cone, bras_kq, npw_kq*nspinor, &
3856 0 : h1_kets_kq, npw_kq*nspinor, czero, gkq_atm_ipc, nb_kq)
3857 : #ifdef HAVE_OPENMP_OFFLOAD
3858 : !$OMP TARGET UPDATE FROM(gkq_atm_ipc)
3859 : #endif
3860 :
3861 : else
3862 : call ZGEMM('C', 'N', nb_kq, nb_k, npw_kq*nspinor, cone, bras_kq, npw_kq*nspinor, &
3863 8592 : h1_kets_kq, npw_kq*nspinor, czero, gkq_atm_ipc, nb_kq)
3864 : end if
3865 :
3866 : ! Transfer data
3867 45824 : gkq_atm(:,:,:,ipc) = gkq_atm_ipc
3868 : end do ! my_ip (loop over my_npert atomic perturbations)
3869 :
3870 : ! Collect gkq_atm inside pert_comm so that all procs can operate on the data.
3871 2864 : if (pert_comm%nproc > 1) call xmpi_sum(gkq_atm, pert_comm%value, ierr)
3872 :
3873 2864 : ABI_FREE(ffnl_k)
3874 2864 : ABI_FREE(ffnl_kq)
3875 2864 : ABI_FREE(kpg_kq)
3876 2864 : ABI_FREE(kpg_k)
3877 2864 : ABI_FREE(kinpw_k)
3878 2864 : ABI_FREE(kinpw_kq)
3879 2864 : ABI_FREE(ph3d_k)
3880 2864 : ABI_FREE(ph3d_kq)
3881 :
3882 : #ifdef HAVE_OPENMP_OFFLOAD
3883 : !$OMP TARGET EXIT DATA MAP(delete:h1_kets_kq, gvnlx1, kets_k, bras_kq) IF (dtset%gpu_option == ABI_GPU_OPENMP)
3884 : !$OMP TARGET EXIT DATA MAP(delete:gs1c_kq) IF (dtset%gpu_option == ABI_GPU_OPENMP .and. sij_opt /= 0)
3885 : #endif
3886 2864 : ABI_FREE(kets_k)
3887 2864 : ABI_FREE(bras_kq)
3888 2864 : ABI_FREE(h1_kets_kq)
3889 2864 : ABI_FREE(gs1c_kq)
3890 2864 : ABI_FREE(gvnlx1)
3891 :
3892 : ! Compute group velocities if we are in transport mode or adaptive gaussian or
3893 : ! tetrahedron with libtetrabz returning nesting condition.
3894 2864 : need_velocities = dtset%eph_transport > 0 .or. fs%eph_fsmear < zero .or. nesting /= 0
3895 :
3896 : if (need_velocities) then
3897 : ! Compute diagonal matrix elements of velocity operator with DFPT routines
3898 : ! Velocities are in Cartesian coordinates.
3899 : !
3900 : ! If k+q is not in the IBZ, we need to recostruct the value by symmetry using v(Sq) = S v(q).
3901 : ! Use transpose(R) because we are using the tables for the wavefunctions
3902 : ! In this case listkk has been called with symrec and use_symrec=False
3903 : ! so q_bz = S^T q_ibz where S is the isym_kq symmetry
3904 :
3905 : !call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk, istwf_k, npw_k, kg_k)
3906 5728 : do in_k=1,nb_k
3907 2864 : band_k = in_k + bstart_k - 1
3908 11456 : vk = vcart_ibz(:, band_k, ik_ibz, spin)
3909 2864 : if (.not. isirr_k) then
3910 2754 : vk = matmul(transpose(cryst%symrel_cart(:,:,isym_k)), vk)
3911 2754 : if (trev_k /= 0) vk = -vk
3912 : end if
3913 : !vk = ddkop%get_vdiag(ebands%eig(band_k, ik_ibz, spin), &
3914 : ! istwf_k, npw_k, wfd%nspinor, kets_k(:,:,in_k), cwaveprj0)
3915 14320 : fs%vk(:,in_k) = vk
3916 : end do
3917 :
3918 : !call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kq, istwf_kq, npw_kq, kg_kq)
3919 5728 : do im_kq=1,nb_kq
3920 2864 : band_kq = im_kq + bstart_kq - 1
3921 11456 : vkq = vcart_ibz(:, band_kq, ikq_ibz, spin)
3922 2864 : if (.not. isirr_kq) then
3923 2833 : vkq = matmul(transpose(cryst%symrel_cart(:,:,isym_kq)), vkq)
3924 2833 : if (trev_kq /= 0) vkq = -vkq
3925 : end if
3926 : !vkq = ddkop%get_vdiag(ebands%eig(band_kq, ikq_ibz, spin), &
3927 : ! istwf_kq, npw_kq, wfd%nspinor, bras_kq(:,:,im_kq), cwaveprj0)
3928 14320 : fs%vkq(:,im_kq) = vkq
3929 : end do
3930 : end if
3931 :
3932 : ! Compute weights for double delta integration at the Fermi level.
3933 : ! Note that we have to call the routine here after the initialization of fs%vk and fs%vkq
3934 : ! TODO: Could precompute weights before entering the loop, apply filter and compute MPI-distribution
3935 2864 : call fs%get_dbldelta_weights(ebands, ik_fs, ik_ibz, ikq_ibz, spin, nesting, dbldelta_wts)
3936 : ! Multiply by the weight of the q-point if we are summing over the IBZ(q).
3937 : !dbldelta_wts = dbldelta_wts * wtk_lgq
3938 :
3939 : ! Accumulate results in tgam (sum over FS k-points and bands for this spin).
3940 11456 : do ipc2=1,natom3
3941 37232 : do ipc1=1,natom3
3942 60144 : do in_k=1,nb_k
3943 77328 : do im_kq=1,nb_kq
3944 77328 : lf = gkq_atm(:, im_kq, in_k, ipc1)
3945 77328 : rg = gkq_atm(:, im_kq, in_k, ipc2)
3946 25776 : res(1) = lf(1) * rg(1) + lf(2) * rg(2)
3947 25776 : res(2) = lf(1) * rg(2) - lf(2) * rg(1)
3948 103104 : tgam(:,ipc1,ipc2) = tgam(:,ipc1,ipc2) + res(:) * dbldelta_wts(im_kq, in_k)
3949 : end do
3950 : end do
3951 : end do
3952 : end do
3953 :
3954 2864 : if (dtset%eph_transport > 0) then
3955 : ! TODO: could almost make this a BLAS call plus a reshape...
3956 5728 : do in_k = 1,nb_k
3957 8592 : do im_kq = 1,nb_kq
3958 14320 : do ipc2 = 1,3
3959 37232 : do ipc1 = 1,3
3960 : ! vk x vk
3961 25776 : vv_kk(ipc1 + (ipc2-1)*3, im_kq, in_k) = fs%vk(ipc1, im_kq) * fs%vk(ipc2, in_k)
3962 : ! vk x vk+q
3963 34368 : vv_kkq(ipc1 + (ipc2-1)*3, im_kq, in_k) = fs%vkq(ipc1, im_kq) * fs%vk(ipc2, in_k)
3964 : end do
3965 : end do
3966 : end do
3967 : end do
3968 :
3969 : ! Accumulate results (sum over FS and bands).
3970 11456 : do ipc2=1,natom3
3971 37232 : do ipc1=1,natom3
3972 60144 : do in_k=1,nb_k
3973 77328 : do im_kq=1,nb_kq
3974 77328 : lf = gkq_atm(:, im_kq, in_k, ipc1)
3975 77328 : rg = gkq_atm(:, im_kq, in_k, ipc2)
3976 : ! res was missing in MJV version!
3977 25776 : res(1) = lf(1) * rg(1) + lf(2) * rg(2)
3978 25776 : res(2) = lf(1) * rg(2) - lf(2) * rg(1)
3979 257760 : resvv_in(1,:) = res(1) * vv_kkq(:,im_kq, in_k)
3980 257760 : resvv_in(2,:) = res(2) * vv_kkq(:,im_kq, in_k)
3981 257760 : resvv_out(1,:) = res(1) * vv_kk(:,im_kq, in_k)
3982 257760 : resvv_out(2,:) = res(2) * vv_kk(:,im_kq, in_k)
3983 721728 : tgamvv_in(:,:,ipc1, ipc2) = tgamvv_in(:,:,ipc1,ipc2) + resvv_in * dbldelta_wts(im_kq, in_k)
3984 747504 : tgamvv_out(:,:,ipc1, ipc2) = tgamvv_out(:,:,ipc1,ipc2) + resvv_out * dbldelta_wts(im_kq, in_k)
3985 : end do
3986 : end do
3987 : end do
3988 : end do
3989 : end if ! add transport things
3990 :
3991 2864 : if (dtset%prteliash == 3) then
3992 : ! Precompute deltas with gaussian (tetra is not supported here
3993 : ! also because one should allocate (nene, nene, mband, mband)
3994 : !call get_dbldelta_weights_emesh()
3995 0 : do in_k=1,nb_k
3996 0 : band_k = in_k + bstart_k - 1
3997 0 : sigma = fs%eph_fsmear
3998 0 : if (fs%eph_fsmear < zero) then
3999 0 : do ii=1,3
4000 0 : abc(ii) = abs(dot_product(fs%vk(:, in_k), fs%kmesh_cartvec(:,ii)))
4001 : end do
4002 0 : sigma = max(maxval(abc), fs%min_smear)
4003 : end if
4004 0 : wt_ek(:, in_k) = gaussian(emesh - ebands%eig(band_k, ik_ibz, spin), sigma) / sqrt(one * fs%nktot)
4005 : end do
4006 :
4007 0 : do im_kq=1,nb_kq
4008 0 : band_kq = im_kq + bstart_kq - 1
4009 0 : sigma = fs%eph_fsmear
4010 0 : if (fs%eph_fsmear < zero) then
4011 0 : do ii=1,3
4012 0 : abc(ii) = abs(dot_product(fs%vkq(:, im_kq), fs%kmesh_cartvec(:,ii)))
4013 : end do
4014 0 : sigma = max(maxval(abc), fs%min_smear)
4015 : end if
4016 0 : wt_ekq(:, im_kq) = gaussian(emesh - ebands%eig(band_kq, ikq_ibz, spin), sigma) / sqrt(one * fs%nktot)
4017 : end do
4018 :
4019 0 : do in_k=1,nb_k
4020 0 : do im_kq=1,nb_kq
4021 0 : do ipc2=1,natom3
4022 0 : do ipc1=1,natom3
4023 0 : lf = gkq_atm(:, im_kq, in_k, ipc1)
4024 0 : rg = gkq_atm(:, im_kq, in_k, ipc2)
4025 0 : res(1) = lf(1) * rg(1) + lf(2) * rg(2)
4026 0 : res(2) = lf(1) * rg(2) - lf(2) * rg(1)
4027 0 : do jene = 1, gams%nene
4028 0 : do iene = 1, gams%nene
4029 : ! TODO: distribute this in procs over q. Make a temp array here for 1 q
4030 : ! then mpi sync it and save it only on 1 processor below after mpisum over k
4031 : tmp_vals_ee(:,iene,jene,ipc1,ipc2) = tmp_vals_ee(:,iene,jene,ipc1,ipc2) + &
4032 0 : res(:) * wt_ekq(iene, im_kq) * wt_ek(jene, in_k)
4033 : end do
4034 : end do
4035 : end do
4036 : end do
4037 : end do
4038 : end do
4039 : end if ! prteliash == 3
4040 :
4041 2864 : ABI_FREE(lambda)
4042 2864 : ABI_FREE(gkq_atm)
4043 : #ifdef HAVE_OPENMP_OFFLOAD
4044 : !$OMP TARGET EXIT DATA MAP(delete:gkq_atm_ipc) IF (dtset%gpu_option == ABI_GPU_OPENMP)
4045 : #endif
4046 2864 : ABI_FREE(gkq_atm_ipc)
4047 :
4048 2925 : if (print_time_k) then
4049 66 : write(msg,'(5x,2(a,i0),a)')"k-point [", my_ik, "/", gams%my_nfsk_q, "]"
4050 66 : call cwtime_report(msg, cpu_k, wall_k, gflops_k)
4051 : end if
4052 : end do ! imyk: sum over k-points on the BZ FS for this spin.
4053 :
4054 61 : ABI_FREE(dbldelta_wts)
4055 61 : ABI_SFREE(wt_ek)
4056 61 : ABI_SFREE(wt_ekq)
4057 :
4058 : ! Collect tgam values inside kpt comm and divide by the total number of k-points in the full mesh.
4059 61 : call xmpi_sum(tgam, kpt_comm%value, ierr)
4060 :
4061 : ! Symmetrize gamma matrices (default). Note that this call is not executed in elphon!
4062 61 : if (gams%symgamma == 1) call tgamma_symm(cryst, gams%qibz(:,iq_ibz), tgam)
4063 :
4064 : ! Save results in gams.
4065 1891 : gams%vals_qibz(:,:,:,iq_ibz, spin) = tgam
4066 :
4067 61 : if (dtset%eph_transport > 0) then
4068 61 : ABI_FREE(vv_kk)
4069 61 : ABI_FREE(vv_kkq)
4070 61 : call xmpi_sum(tgamvv_in, kpt_comm%value, ierr)
4071 61 : call xmpi_sum(tgamvv_out, kpt_comm%value, ierr)
4072 15616 : gams%vals_in_qibz(:,:,:,:,iq_ibz,spin) = tgamvv_in
4073 15616 : gams%vals_out_qibz(:,:,:,:,iq_ibz,spin) = tgamvv_out
4074 : end if
4075 :
4076 183 : if (dtset%prteliash == 3) then
4077 0 : call xmpi_sum(tmp_vals_ee, kpt_comm%value, ierr)
4078 : !TODO: write arrays to ncfile.
4079 0 : gams%vals_ee(:,:,:,:,:, iq_ibz, spin) = tmp_vals_ee
4080 0 : ABI_FREE(tmp_vals_ee)
4081 0 : ABI_FREE(emesh)
4082 : end if
4083 : end do ! spin
4084 :
4085 61 : ABI_FREE(v1scf)
4086 61 : ABI_FREE(vlocal1)
4087 :
4088 61 : write(msg,'(2(a,i0),a)')" Computation of q-point [", iq_ibz, "/", gams%nqibz, "]"
4089 127 : call cwtime_report(msg, cpu_q, wall_q, gflops_q, end_str=ch10)
4090 : end do ! iq_ibz
4091 :
4092 5 : call cwtime_report(" phonon linewidths q-loop", cpu_all, wall_all, gflops_all, end_str=ch10)
4093 :
4094 : ! Free memory
4095 5 : ABI_FREE(grad_berry)
4096 5 : ABI_FREE(dummy_vtrial)
4097 5 : ABI_FREE(work)
4098 5 : ABI_FREE(ph1d)
4099 5 : ABI_FREE(vlocal)
4100 5 : ABI_FREE(kg_k)
4101 5 : ABI_FREE(kg_kq)
4102 5 : ABI_FREE(tgam)
4103 5 : ABI_FREE(displ_cart)
4104 5 : ABI_FREE(displ_red)
4105 : !ABI_FREE(qibz_done)
4106 5 : ABI_SFREE(vcart_ibz)
4107 5 : call krank%free()
4108 :
4109 5 : ABI_SFREE(tgamvv_in)
4110 5 : ABI_SFREE(tgamvv_out)
4111 5 : call ddkop%free(); call gs_hamkq%free(); call wfd%free()
4112 10 : do spin=1,ebands%nsppol
4113 10 : call fstab(spin)%free()
4114 : end do
4115 10 : ABI_FREE(fstab)
4116 5 : call pawcprj_free(cwaveprj0)
4117 5 : ABI_FREE(cwaveprj0)
4118 :
4119 : ! Collect results on each node
4120 5 : call cwtime(cpu_all, wall_all, gflops_all, "start")
4121 5 : call xmpi_sum(gams%vals_qibz, qs_comm%value, ierr)
4122 5 : if (dtset%eph_transport > 0) then
4123 5 : call xmpi_sum(gams%vals_in_qibz, qs_comm%value, ierr)
4124 5 : call xmpi_sum(gams%vals_out_qibz, qs_comm%value, ierr)
4125 : end if
4126 5 : if (dtset%prteliash == 3) call xmpi_sum(gams%vals_ee, qs_comm%value, ierr)
4127 5 : call cwtime_report(" xmpi_sum", cpu_all, wall_all, gflops_all)
4128 :
4129 : ! Deallocate MPI communicators.
4130 5 : call pert_comm%free(); call qpt_comm%free(); call bsum_comm%free(); call qs_comm%free()
4131 5 : call kpt_comm%free(); call spin_comm%free(); call pkb_comm%free()
4132 : !call ncwrite_comm%free()
4133 :
4134 : ! Print gamma(IBZ) to ab_out and ncid.
4135 5 : if (my_rank == master) call phgamma_ncwrite(gams, cryst, ifc, ncid)
4136 :
4137 : ! Interpolate linewidths along the q-path.
4138 5 : if (dtset%ph_nqpath <= 0) then
4139 : write(msg, '(7a,es16.6,4a)' ) &
4140 0 : 'You have not specified a path for the linewidth calculation - no interpolation or output will be done ',ch10,&
4141 0 : 'Action: check your input variables ph_nqpath and ph_qpath'
4142 0 : ABI_WARNING(msg)
4143 : else
4144 5 : call gams%linwid(cryst, ifc, dtset%ph_ndivsm, dtset%ph_nqpath, dtset%ph_qpath, dtfil%filnam_ds(4), ncid, wminmax, comm)
4145 : end if
4146 :
4147 : ! Compute a2Fw using the ab-initio q-points (no interpolation here)
4148 : call a2fw%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear, &
4149 5 : dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qintp=.False., qptopt=1)
4150 :
4151 5 : if (my_rank == master) then
4152 5 : call a2fw%write(strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid)
4153 5 : if (dtset%prteliash == 3) call a2fw_ee_write(a2fw, strcat(dtfil%filnam_ds(4), "_NOINTP"))
4154 : end if
4155 5 : call a2fw%free()
4156 :
4157 : ! Compute a2Fw using Fourier interpolation (R -> q) and ph_ngqpt grid.
4158 : call a2fw%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear, &
4159 5 : dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=1)
4160 :
4161 5 : if (my_rank == master) then
4162 5 : call a2fw%write(dtfil%filnam_ds(4), "_qintp", ncid)
4163 5 : if (dtset%prteliash == 3) call a2fw_ee_write(a2fw, dtfil%filnam_ds(4))
4164 : end if
4165 :
4166 5 : call a2fw%free()
4167 :
4168 : ! Compute A2fw using Fourier interpolation and full BZ for debugging purposes.
4169 : !call a2fw%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,&
4170 : ! dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=3)
4171 : !if (my_rank == master) call a2fw%write(strcat(dtfil%filnam_ds(4), "_A2FW_QPTOPT3"), "fake", nctk_noid)
4172 : !call a2fw%free()
4173 :
4174 5 : if (dtset%eph_transport == 1) then
4175 : ! Calculate and output transport quantities
4176 :
4177 : ! Compute a2Fw_tr using ab-initio q-points (no interpolation)
4178 : call a2fw_tr%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,&
4179 5 : dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qintp=.False., qptopt=1)
4180 :
4181 5 : if (my_rank == master) call a2fw_tr%write(strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid)
4182 5 : call a2fw_tr%free()
4183 :
4184 : ! Compute a2Fw_tr using Fourier interpolation (R --> q) and ph_ngqpt grid
4185 : call a2fw_tr%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,&
4186 5 : dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=1)
4187 :
4188 5 : if (my_rank == master) call a2fw_tr%write(dtfil%filnam_ds(4), "_qintp", ncid)
4189 5 : call a2fw_tr%free()
4190 :
4191 : ! Compute A2fw_tr using Fourier interpolation and full BZ for debugging purposes.
4192 : !call a2fw_tr%init(gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,&
4193 : ! dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=3)
4194 : !if (my_rank == master) call a2fw_tr%write(strcat(dtfil%filnam_ds(4), "_A2FWTR_QPTOPT3"), "fake", nctk_noid)
4195 : !call a2fw_tr%free()
4196 : end if
4197 :
4198 5 : if (my_rank == master) then
4199 5 : NCF_CHECK(nf90_close(ncid))
4200 : end if
4201 :
4202 5 : call gams%free()
4203 :
4204 115 : end subroutine eph_phgamma
4205 : !!***
4206 :
4207 : !----------------------------------------------------------------------
4208 :
4209 : !!****f* m_phgamma/phgamma_setup_qpoint
4210 : !! NAME
4211 : !! phgamma_setup_qpoint
4212 : !!
4213 : !! FUNCTION
4214 : !!
4215 : !! SOURCE
4216 :
4217 61 : subroutine phgamma_setup_qpoint(gams, fs, cryst, ebands, spin, ltetra, qpt, nesting, comm)
4218 :
4219 : !Arguments ------------------------------------
4220 : class(phgamma_t),intent(inout) :: gams
4221 : type(fstab_t),intent(inout) :: fs
4222 : type(crystal_t),intent(in) :: cryst
4223 : type(ebands_t),intent(in) :: ebands
4224 : integer,intent(in) :: spin, ltetra, comm
4225 : integer,intent(out) :: nesting
4226 : !arrays
4227 : real(dp),intent(in) :: qpt(3)
4228 :
4229 : !Local variables-------------------------------
4230 : !scalars
4231 : integer,parameter :: enough = 5
4232 : integer :: nkbz, ierr, nb, ik_bz, ik_ibz, ikq_ibz, ikq_fs, ik_fs, i1, i2, i3, nkfs_q, nene, ib1, ib2 ! band_k, band_kq
4233 : real(dp),parameter :: max_occ1 = one
4234 : real(dp) :: cpu, wall, gflops, enemin, enemax
4235 : character(len=500) :: msg
4236 61 : type(krank_t) :: ibz_krank
4237 61 : type(t_tetrahedron) :: tetra
4238 : character(len=80) :: errorstring
4239 : !arrays
4240 : integer :: nge(3), ngw(3), g0bz_kq(3)
4241 61 : integer,allocatable :: select_ikfs(:), indkpt(:), kbz2fs(:) !, symrecfm(:,:,:)
4242 : real(dp) :: kk(3), kq(3)
4243 61 : real(dp),allocatable :: eig_k(:,:), eig_kq(:,:), wght_bz(:,:,:), kbz(:,:)
4244 61 : real(dp),allocatable :: work_k(:), work_kq(:), dtweightde(:,:,:), tweight(:,:,:)
4245 : ! *************************************************************************
4246 :
4247 61 : call cwtime(cpu, wall, gflops, "start")
4248 :
4249 : ! The double delta with tetra is ill-defined for q == 0. Set nesting to 1 and return
4250 61 : nesting = 0
4251 76 : if (abs(fs%eph_intmeth) == 2 .and. all(abs(qpt) < tol12)) then
4252 2 : ABI_COMMENT("Tetrahedron for double grid with q = 0 is ill-defined. Using adaptive gaussian.")
4253 2 : nesting = 1
4254 : end if
4255 :
4256 61 : if (fs%eph_intmeth == 1 .or. nesting == 1) then
4257 : ! Gaussian method:
4258 : ! Distribute k-points within the FS window inside comm.
4259 : ! 1) Select k-points such that k+q is still inside the FS window
4260 : ! 2) Distribute effective k-points assuming all procs in comm have all FS k-points (no filtering).
4261 78 : ABI_MALLOC(select_ikfs, (fs%nkfs))
4262 26 : nkfs_q = 0
4263 3614 : do ik_fs=1,fs%nkfs
4264 14352 : kq = fs%kpts(:, ik_fs) + qpt
4265 3588 : ikq_fs = fs%findkg0(kq, g0bz_kq); if (ikq_fs == -1) cycle
4266 1536 : nkfs_q = nkfs_q + 1
4267 3614 : select_ikfs(nkfs_q) = ik_fs
4268 : end do
4269 26 : ABI_SFREE(gams%my_ifsk_q)
4270 26 : call xmpi_split_list(nkfs_q, select_ikfs, comm, gams%my_nfsk_q, gams%my_ifsk_q)
4271 26 : ABI_FREE(select_ikfs)
4272 26 : call cwtime_report(" phgamma_setup_qpoint", cpu, wall, gflops)
4273 26 : return
4274 : end if
4275 :
4276 : ! Tetrahedron method:
4277 : ! 1) Compute weights for double delta integration.
4278 : ! 2) Filter k-points according to the weights (k must be in FS window with non-zero weight)
4279 : ! 3) Distribute effective k-points assuming all procs have all FS k-points.
4280 :
4281 33845 : fs%dbldelta_tetra_weights_kfs = zero
4282 :
4283 35 : ABI_CHECK(isdiagmat(ebands%kptrlatt), "kptrlatt must be diagonal when tetra is used.")
4284 35 : ABI_CHECK(ebands%nshiftk == 1, "nshiftk must be 1 when tetra is used")
4285 35 : nge = get_diag(ebands%kptrlatt); ngw = nge
4286 140 : nkbz = product(nge(1:3))
4287 :
4288 : ! TODO: Handle symmetries in a cleaner way. Change API of krank_new to pass symafm and kptopt
4289 : call ibz_krank%init(ebands%nkpt, ebands%kptns, nsym=cryst%nsym, symrec=cryst%symrec, &
4290 35 : time_reversal=kpts_timrev_from_kptopt(ebands%kptopt) == 1)
4291 :
4292 : ! Compute eig_k and eig_kq in full BZ for the relevant bands around Ef.
4293 35 : nb = fs%maxnb
4294 105 : ABI_MALLOC(kbz, (3, nkbz))
4295 105 : ABI_MALLOC(indkpt, (nkbz))
4296 140 : ABI_MALLOC(eig_k, (nb, nkbz))
4297 105 : ABI_MALLOC(eig_kq, (nb, nkbz))
4298 70 : ABI_MALLOC(kbz2fs, (nkbz))
4299 :
4300 17955 : kbz2fs = -1; ierr = 0; ik_bz = 0
4301 315 : do i3=0,nge(3) - 1
4302 2555 : do i2=0,nge(2) - 1
4303 20440 : do i1=0,nge(1) - 1
4304 17920 : ik_bz = ik_bz + 1
4305 17920 : indkpt(ik_bz) = ik_bz
4306 :
4307 : ! Find correspondence between libtetra mesh and the IBZ.
4308 71680 : kk = ([i1, i2, i3] + ebands%shiftk(:, 1)) / nge(:)
4309 71680 : kbz(:, ik_bz) = kk
4310 :
4311 17920 : ik_fs = fs%krank%get_index(kk)
4312 17920 : if (ik_fs /= -1) then
4313 4830 : kbz2fs(ik_bz) = ik_fs
4314 : else
4315 : !ABI_ERROR(sjoin('kpt:', trim(ktoa(kk)), 'is not in FS!!'))
4316 : end if
4317 :
4318 17920 : ik_ibz = ibz_krank%get_index(kk)
4319 17920 : if (ik_ibz < 1) then
4320 0 : if (ierr <= enough) then
4321 0 : ABI_WARNING(sjoin('kpt:', trim(ktoa(kk)), 'has no symmetric among the k-points!'))
4322 : end if
4323 0 : ierr = ierr + 1; cycle
4324 : end if
4325 :
4326 53760 : eig_k(:, ik_bz) = ebands%eig(fs%bmin:fs%bmax, ik_ibz, spin)
4327 :
4328 : ! Find correspondence between k+q in the BZ and the IBZ image.
4329 71680 : kq = kk + qpt
4330 17920 : ikq_ibz = ibz_krank%get_index(kq)
4331 :
4332 17920 : if (ikq_ibz < 1) then
4333 0 : if (ierr <= enough) then
4334 0 : ABI_WARNING(sjoin('kpt + qpt:', trim(ktoa(kq)), 'has no symmetric among the k-points!'))
4335 : end if
4336 0 : ierr = ierr + 1; cycle
4337 : end if
4338 :
4339 56000 : eig_kq(:, ik_bz) = ebands%eig(fs%bmin:fs%bmax, ikq_ibz, spin)
4340 : end do
4341 : end do
4342 : end do
4343 :
4344 35 : ABI_CHECK(ierr == 0, "See above warnings")
4345 35 : call ibz_krank%free()
4346 :
4347 70 : select case (ltetra)
4348 : case (1, 2)
4349 : ! Use libtetra routines.
4350 : ! Compute weights for double delta integration. Note that libtetra assumes Ef set to zero.
4351 : ! TODO: Average weights over degenerate states?
4352 35 : if (ltetra == 1) call wrtout(std_out, " Using linear tetrahedron method from libtetrabz (ltetra: 1)")
4353 35 : if (ltetra == 2) call wrtout(std_out, " Using optimized tetrahedron method from libtetrabz (ltetra: 2)")
4354 107555 : eig_k = eig_k - ebands%fermie; eig_kq = eig_kq - ebands%fermie
4355 175 : ABI_MALLOC(wght_bz, (nb, nb, nkbz))
4356 35 : call libtetrabz_dbldelta(ltetra, cryst%gprimd, nb, nge, eig_k, eig_kq, ngw, wght_bz, comm=comm)
4357 : ! Revert changes to eig arrays
4358 53795 : eig_k = eig_k + ebands%fermie !; eig_kq = eig_kq + ebands%fermie
4359 :
4360 : ! Reindex from full BZ to fs% kpoints.
4361 17955 : do ik_bz=1,nkbz
4362 17920 : ik_fs = kbz2fs(ik_bz)
4363 17955 : if (ik_fs /= -1) then
4364 33810 : fs%dbldelta_tetra_weights_kfs(:,:,ik_fs) = wght_bz(:,:,ik_bz)
4365 : else
4366 : !write(std_out,*)"should be zero :", wght_bz(:,:,ik_bz)
4367 : end if
4368 : end do
4369 35 : ABI_FREE(wght_bz)
4370 :
4371 : case (3)
4372 : ! Tetrahedron method with Allen's approach for double delta.
4373 : ! WARNING: wrong results maybe bug somewhere.
4374 0 : call wrtout(std_out, " Calling Allen's version for q-point")
4375 0 : nene = 3
4376 0 : enemin = ebands%fermie - tol6
4377 0 : enemax = ebands%fermie + tol6
4378 :
4379 0 : ABI_MALLOC(dtweightde, (nkbz, nene, nene))
4380 0 : ABI_MALLOC(tweight, (nkbz, nene, nene))
4381 0 : ABI_MALLOC(work_k, (nkbz))
4382 0 : ABI_MALLOC(work_kq, (nkbz))
4383 :
4384 0 : call init_tetra(indkpt, cryst%gprimd, fs%klatt, kbz, nkbz, tetra, ierr, errorstring, comm)
4385 0 : ABI_CHECK(ierr == 0, "init_tetra returned ierr /= 0")
4386 :
4387 0 : fs%dbldelta_tetra_weights_kfs = zero
4388 0 : do ib2=1,nb
4389 0 : work_k = eig_k(ib2, :)
4390 0 : do ib1=1,nb
4391 0 : work_kq = eig_kq(ib1, :)
4392 : ! TODO: Average weights over degenerate states?
4393 : call get_dbl_tetra_weight(work_k, work_kq, enemin, enemax, enemin, enemax, &
4394 0 : max_occ1, nene, nene, nkbz, tetra, tweight, dtweightde, ierr)
4395 0 : ABI_CHECK(ierr == 0, "get_dbldelta_weights returned ierr /= 0")
4396 :
4397 : ! Reindex from full BZ to fs% kpoints.
4398 0 : do ik_bz=1,nkbz
4399 0 : ik_fs = kbz2fs(ik_bz)
4400 0 : if (ik_fs /= -1) then
4401 0 : fs%dbldelta_tetra_weights_kfs(ib1,ib2,ik_fs) = dtweightde(ik_bz, 2, 2) ! / nkbz
4402 : else
4403 : !write(std_out,*)"should be zero :", wght_bz(:,:,ik_bz)
4404 : end if
4405 : end do
4406 : end do
4407 : end do
4408 :
4409 0 : call destroy_tetra(tetra)
4410 0 : ABI_FREE(work_k)
4411 0 : ABI_FREE(work_kq)
4412 0 : ABI_FREE(dtweightde)
4413 0 : ABI_FREE(tweight)
4414 :
4415 : case default
4416 35 : ABI_ERROR(sjoin("Invalid value of ltetra:", itoa(ltetra)))
4417 : end select
4418 :
4419 : ! Now we can filter the k-points according to the tetra weights and distribute them inside comm.
4420 : ! Assuming all procs in comm have all k-points in the IBZ.
4421 : ! nkfs_q is the total number of BZ k-points on the FS contributing for this q-point.
4422 105 : ABI_MALLOC(select_ikfs, (fs%nkfs))
4423 35 : nkfs_q = 0
4424 4865 : do ik_fs=1,fs%nkfs
4425 6079 : if (any(abs(fs%dbldelta_tetra_weights_kfs(:,:,ik_fs)) > zero)) then
4426 4668 : nkfs_q = nkfs_q + 1; select_ikfs(nkfs_q) = ik_fs
4427 : end if
4428 : end do
4429 :
4430 : ! Here we compute my_ifsk_q.
4431 35 : ABI_SFREE(gams%my_ifsk_q)
4432 35 : call xmpi_split_list(nkfs_q, select_ikfs, comm, gams%my_nfsk_q, gams%my_ifsk_q)
4433 35 : ABI_FREE(select_ikfs)
4434 :
4435 : !write(std_out,"(2(a,i0),/)")" Treating ", gams%my_nfsk_q, " my k-points in the FS window over total nkfs: ", fs%nkfs
4436 :
4437 : write(msg, "(2(a,i0),a)") &
4438 35 : " Number of k-points in the FS window treated by me: ", gams%my_nfsk_q, " over: ", nkfs_q, ch10
4439 : !" Number of MPI procs in kpt_comm: ", gams%kpt_comm%nproc
4440 35 : call wrtout(std_out, msg)
4441 :
4442 35 : ABI_FREE(kbz2fs)
4443 35 : ABI_FREE(indkpt)
4444 35 : ABI_FREE(kbz)
4445 35 : ABI_FREE(eig_k)
4446 35 : ABI_FREE(eig_kq)
4447 :
4448 35 : call cwtime_report(sjoin(" phgamma_setup_qpoint for spin:", itoa(spin)), cpu, wall, gflops)
4449 :
4450 122 : end subroutine phgamma_setup_qpoint
4451 : !!***
4452 :
4453 : !----------------------------------------------------------------------
4454 :
4455 : !!****f* m_phgamma/find_ewin
4456 : !! NAME
4457 : !! find_ewin
4458 : !!
4459 : !! FUNCTION
4460 : !! Use bisection to find the optimal energy window around the Fermi level
4461 : !!
4462 : !! SOURCE
4463 :
4464 0 : subroutine find_ewin(nqibz, qibz, cryst, ebands, ltetra, fs_ewin, comm)
4465 :
4466 : !Arguments ------------------------------------
4467 : integer,intent(in) :: nqibz
4468 : type(crystal_t),intent(in) :: cryst
4469 : type(ebands_t),intent(in) :: ebands
4470 : integer,intent(in) :: ltetra, comm
4471 : real(dp),intent(out) :: fs_ewin
4472 : !arrays
4473 : real(dp),intent(in) :: qibz(3, nqibz)
4474 :
4475 : !Local variables-------------------------------
4476 : !scalars
4477 : integer,parameter :: master = 0
4478 : integer :: my_rank, iq_ibz, spin, ii, unt
4479 : real(dp) :: cpu, wall, gflops
4480 : character(len=500) :: msg
4481 : !arrays
4482 : integer :: bstarts(3), bstops(3)
4483 : real(dp) :: elows(3), ehighs(3), ewins(3), qsums(3)
4484 0 : real(dp), allocatable :: wtqs(:,:,:)
4485 : ! *************************************************************************
4486 :
4487 0 : call cwtime(cpu, wall, gflops, "start")
4488 0 : my_rank = xmpi_comm_rank(comm)
4489 :
4490 0 : unt = -1
4491 0 : if (xmpi_comm_rank(comm) == master) then
4492 0 : if (open_file("tetra.dat", msg, newunit=unt, action="write", form="formatted", status="unknown") /= 0) then
4493 0 : ABI_ERROR(msg)
4494 : end if
4495 : end if
4496 :
4497 : ! 1 is low, 2 is high, 3 is for the workspace
4498 0 : ABI_MALLOC(wtqs, (nqibz, ebands%nsppol, 3))
4499 0 : ewins(1) = half * eV_Ha; elows(1) = ebands%fermie - ewins(1); ehighs(1) = ebands%fermie + ewins(1)
4500 0 : call ebands%get_bands_from_erange(elows(1), ehighs(1), bstarts(1), bstops(1))
4501 0 : call calc_dbldelta(cryst, ebands, ltetra, bstarts(1), bstops(1), nqibz, qibz, wtqs(:,:,1), comm)
4502 0 : if (unt /= -1) write(unt, *) wtqs(:,:,1)
4503 :
4504 : !ewins(2) = five * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2)
4505 : ewins(2) = two * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2)
4506 0 : ewins(2) = one * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2)
4507 0 : call ebands%get_bands_from_erange(elows(2), ehighs(2), bstarts(2), bstops(2))
4508 0 : call calc_dbldelta(cryst, ebands, ltetra, bstarts(2), bstops(2), nqibz, qibz, wtqs(:,:,2), comm)
4509 0 : if (unt /= -1) write(unt, *) wtqs(:,:,2)
4510 :
4511 0 : if (abs(sum(abs(wtqs(:,:,1)) - sum(abs(wtqs(:,:,2))))) / sum(abs(wtqs(:,:,2))) < tol2) then
4512 0 : fs_ewin = ewins(1)
4513 0 : call print_weights_index(1); goto 100
4514 : end if
4515 :
4516 : ! Bisection part.
4517 0 : do
4518 0 : ewins(3) = (ewins(1) + ewins(2)) / two; elows(3) = ebands%fermie - ewins(3); ehighs(3) = ebands%fermie + ewins(3)
4519 0 : call ebands%get_bands_from_erange(elows(3), ehighs(3), bstarts(3), bstops(3))
4520 0 : call calc_dbldelta(cryst, ebands, ltetra, bstarts(3), bstops(3), nqibz, qibz, wtqs(:,:,3), comm)
4521 0 : if (unt /= -1) write(unt, *) wtqs(:,:,3)
4522 :
4523 0 : do ii=1,3
4524 0 : qsums(ii) = sum(abs(wtqs(:,:,ii))) / nqibz
4525 : end do
4526 :
4527 0 : write(std_out,*)qsums(1), ", for ewins(1)_eV", ewins(1) * Ha_eV, " bstart(1)", bstarts(1), " bstopt(1)", bstops(1)
4528 0 : write(std_out,*)qsums(2), ", for ewins(2)_eV", ewins(2) * Ha_eV, " bstart(1)", bstarts(2), " bstopt(1)", bstops(2)
4529 0 : write(std_out,*)qsums(3), ", for ewins(3)_eV", ewins(3) * Ha_eV, " bstart(1)", bstarts(3), " bstopt(1)", bstops(3)
4530 :
4531 : !if (abs(sum(abs(wtqs(:,:,3)) - sum(abs(wtqs(:,:,2))))) / sum(abs(wtqs(:,:,2))) < tol3) then
4532 : ! ! If 3 is close to 2, try to reduce the window by moving towards 1
4533 : ! ewins(2) = ewins(3); elows(2) = elows(3); ehighs(2) = ehighs(3)
4534 :
4535 0 : if (abs(qsums(3) - qsums(2)) / qsums(2) < tol2) then
4536 0 : fs_ewin = ewins(3)
4537 0 : call print_weights_index(3)
4538 0 : exit
4539 :
4540 : else
4541 0 : ewins(1) = ewins(3); elows(1) = elows(3); ehighs(1) = ehighs(3)
4542 : end if
4543 :
4544 : end do
4545 :
4546 : 100 continue
4547 0 : ABI_FREE(wtqs)
4548 0 : if (unt /= -1) close(unt)
4549 :
4550 0 : call cwtime_report(" find_ewin", cpu, wall, gflops)
4551 :
4552 : contains
4553 :
4554 0 : subroutine print_weights_index(ind)
4555 :
4556 : integer,intent(in) :: ind
4557 :
4558 0 : if (my_rank == master) then
4559 0 : write(std_out, "(a, f8.3, a)")" Optimal FS energy window: ", fs_ewin * Ha_eV, " (eV)"
4560 0 : do iq_ibz=1,nqibz
4561 0 : write(std_out, "(2a)")" For q-point:", trim(ktoa(qibz(:, iq_ibz)))
4562 0 : do spin=1,ebands%nsppol
4563 0 : write(std_out, *) wtqs(iq_ibz, spin, ind)
4564 : end do
4565 : end do
4566 : end if
4567 :
4568 0 : end subroutine print_weights_index
4569 :
4570 : end subroutine find_ewin
4571 : !!***
4572 :
4573 : !----------------------------------------------------------------------
4574 :
4575 : !!****f* m_phgamma/calc_dbldelta
4576 : !! NAME
4577 : !! calc_dbldelta
4578 : !!
4579 : !! FUNCTION
4580 : !! Compute Tetrahedron weights for the double delta.
4581 : !!
4582 : !! SOURCE
4583 :
4584 0 : subroutine calc_dbldelta(cryst, ebands, ltetra, bstart, bstop, nqibz, qibz, wtqs, comm)
4585 :
4586 : !Arguments ------------------------------------
4587 : type(crystal_t),intent(in) :: cryst
4588 : type(ebands_t),intent(in) :: ebands
4589 : integer,intent(in) :: ltetra, bstart, bstop, nqibz, comm
4590 : !arrays
4591 : real(dp),intent(in) :: qibz(3,nqibz)
4592 : real(dp),intent(out) :: wtqs(nqibz, ebands%nsppol)
4593 :
4594 : !Local variables-------------------------------
4595 : !scalars
4596 : integer,parameter :: enough = 50
4597 : integer :: nkbz, ierr, nb, ik_bz, ik_ibz, ikq_ibz, i1, i2, i3, nene, spin, iq_ibz
4598 : integer :: ib1, ib2, my_rank, nproc, cnt
4599 : real(dp),parameter :: max_occ1 = one
4600 : real(dp) :: enemin, enemax !cpu, wall, gflops,
4601 : !character(len=500) :: msg
4602 0 : type(krank_t) :: ibz_krank
4603 0 : type(t_tetrahedron) :: tetra
4604 : !character(len=80) :: errorstring
4605 : !arrays
4606 : integer :: nge(3), ngw(3)
4607 0 : integer,allocatable :: indkpt(:) !, symrecfm(:,:,:)
4608 : real(dp) :: qpt(3), kk(3), kq(3)
4609 0 : real(dp),allocatable :: eig_k(:,:), eig_kq(:,:), wght_bz(:,:,:), kbz(:,:)
4610 0 : real(dp),allocatable :: work_k(:), work_kq(:), dtweightde(:,:,:), tweight(:,:,:)
4611 : ! *************************************************************************
4612 :
4613 0 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
4614 : !call cwtime(cpu, wall, gflops, "start")
4615 :
4616 0 : ABI_CHECK(nqibz > 1, "Need more that 1 q-point")
4617 0 : ABI_CHECK(isdiagmat(ebands%kptrlatt), "kptrlatt must be diagonal when tetra is used.")
4618 0 : ABI_CHECK(ebands%nshiftk == 1, "nshiftk must be 1 when tetra is used")
4619 0 : nge = get_diag(ebands%kptrlatt); ngw = nge
4620 0 : nkbz = product(nge(1:3))
4621 :
4622 : ! TODO: Handle symmetries in a cleaner way. Change API of krank_new to pass symafm and kptopt
4623 : call ibz_krank%init(ebands%nkpt, ebands%kptns, nsym=cryst%nsym, symrec=cryst%symrec, &
4624 0 : time_reversal=kpts_timrev_from_kptopt(ebands%kptopt) == 1)
4625 :
4626 0 : nb = bstop - bstart + 1
4627 0 : ABI_MALLOC(kbz, (3, nkbz))
4628 0 : ABI_MALLOC(indkpt, (nkbz))
4629 0 : ABI_MALLOC(eig_k, (nb, nkbz))
4630 0 : ABI_MALLOC(eig_kq, (nb, nkbz))
4631 :
4632 0 : wtqs = zero; cnt = 0
4633 0 : do spin=1,ebands%nsppol
4634 0 : do iq_ibz=1,nqibz
4635 0 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism.
4636 0 : qpt = qibz(:, iq_ibz)
4637 :
4638 : ! The double delta is ill-defined for q == 0
4639 0 : if (all(abs(qpt) < tol12)) cycle
4640 :
4641 0 : ierr = 0; ik_bz = 0
4642 0 : do i3=0,nge(3) - 1
4643 0 : do i2=0,nge(2) - 1
4644 0 : do i1=0,nge(1) - 1
4645 0 : ik_bz = ik_bz + 1
4646 0 : indkpt(ik_bz) = ik_bz
4647 :
4648 : ! Find correspondence between the grid and the IBZ
4649 0 : kk = ([i1, i2, i3] + ebands%shiftk(:, 1)) / nge(:)
4650 0 : kbz(:, ik_bz) = kk
4651 :
4652 0 : ik_ibz = ibz_krank%get_index(kk)
4653 0 : if (ik_ibz < 1) then
4654 0 : if (ierr <= enough) then
4655 0 : ABI_WARNING(sjoin('kpt:', trim(ktoa(kk)), 'has no symmetric among the k-points!'))
4656 : end if
4657 0 : ierr = ierr + 1; cycle
4658 : end if
4659 :
4660 0 : eig_k(:, ik_bz) = ebands%eig(bstart:bstop, ik_ibz, spin)
4661 :
4662 : ! Find correspondence between the BZ grid and the IBZ.
4663 0 : kq = kk + qpt
4664 0 : ikq_ibz = ibz_krank%get_index(kq)
4665 :
4666 0 : if (ikq_ibz < 1) then
4667 0 : if (ierr <= enough) then
4668 0 : ABI_WARNING(sjoin('kpt + qpt:', trim(ktoa(kq)), 'has no symmetric among the k-points!'))
4669 : end if
4670 0 : ierr = ierr + 1; cycle
4671 : end if
4672 :
4673 0 : eig_kq(:, ik_bz) = ebands%eig(bstart:bstop, ikq_ibz, spin)
4674 : end do
4675 : end do
4676 : end do
4677 0 : ABI_CHECK(ierr == 0, "See above warnings")
4678 :
4679 0 : if (any(ltetra == [1, 2])) then
4680 : ! Compute weights for double delta integration. Note that libtetra assumes Ef set to zero.
4681 : ! TODO: Average weights over degenerate states?
4682 : !write(std_out,"(a,i0,2a)")" Calling libtetrabz_dbldelta with ltetra: ", ltetra, " for q-point:", trim(ktoa(qpt))
4683 0 : eig_k = eig_k - ebands%fermie; eig_kq = eig_kq - ebands%fermie
4684 0 : ABI_MALLOC(wght_bz, (nb, nb, nkbz))
4685 0 : call libtetrabz_dbldelta(ltetra, cryst%gprimd, nb, nge, eig_k, eig_kq, ngw, wght_bz, comm=comm)
4686 0 : eig_k = eig_k + ebands%fermie !; eig_kq = eig_kq + ebands%fermie
4687 :
4688 0 : wtqs(iq_ibz, spin) = sum(abs(wght_bz))
4689 0 : ABI_FREE(wght_bz)
4690 :
4691 0 : else if (ltetra == 3) then
4692 : ! Tetrahedron method with Allen's approach for double delta.
4693 : !write(std_out,"(2a)")" Calling Allen's version for q-point: ", trim(ktoa(qpt))
4694 0 : nene = 3
4695 0 : enemin = ebands%fermie - tol6
4696 0 : enemax = ebands%fermie + tol6
4697 :
4698 0 : ABI_MALLOC(dtweightde, (nkbz, nene, nene))
4699 0 : ABI_MALLOC(tweight, (nkbz, nene, nene))
4700 0 : ABI_MALLOC(work_k, (nkbz))
4701 0 : ABI_MALLOC(work_kq, (nkbz))
4702 :
4703 : ! TODO
4704 0 : ABI_ERROR("Not Implemented")
4705 : !call init_tetra(indkpt, cryst%gprimd, fs%klatt, kbz, nkbz, tetra, ierr, errorstring, comm)
4706 0 : ABI_CHECK(ierr == 0, "init_tetra returned ierr /= 0")
4707 :
4708 : !fs%dbldelta_tetra_weights_kfs = zero
4709 0 : do ib2=1,nb
4710 0 : work_k = eig_k(ib2, :)
4711 0 : do ib1=1,nb
4712 0 : work_kq = eig_kq(ib1, :)
4713 : ! TODO: Average weights over degenerate states?
4714 : call get_dbl_tetra_weight(work_k, work_kq, enemin, enemax, enemin, enemax, &
4715 0 : max_occ1, nene, nene, nkbz, tetra, tweight, dtweightde, ierr)
4716 0 : ABI_CHECK(ierr == 0, "get_dbldelta_weights returned ierr /= 0")
4717 :
4718 0 : wtqs(iq_ibz, spin) = wtqs(iq_ibz, spin) + sum(abs(dtweightde(:, 2, 2)))
4719 : end do
4720 : end do
4721 :
4722 0 : call destroy_tetra(tetra)
4723 0 : ABI_FREE(work_k)
4724 0 : ABI_FREE(work_kq)
4725 0 : ABI_FREE(dtweightde)
4726 0 : ABI_FREE(tweight)
4727 : else
4728 0 : ABI_ERROR(sjoin("Invalid value of ltetra:", itoa(ltetra)))
4729 : end if
4730 : end do ! iq_ibz
4731 : end do ! spin
4732 :
4733 0 : ABI_FREE(indkpt)
4734 0 : ABI_FREE(kbz)
4735 0 : ABI_FREE(eig_k)
4736 0 : ABI_FREE(eig_kq)
4737 :
4738 0 : call ibz_krank%free()
4739 0 : call xmpi_sum(wtqs, comm, ierr)
4740 : !call cwtime_report(" calc_dbldelta", cpu, wall, gflops)
4741 :
4742 0 : end subroutine calc_dbldelta
4743 : !!***
4744 :
4745 75 : end module m_phgamma
4746 : !!***
|
