Branch data Line data Source code
1 : : !!****m* ABINIT/m_cumulant
2 : : !! NAME
3 : : !! m_cumulant
4 : : !!
5 : : !! FUNCTION
6 : : !! Module to compute the e-ph spectral functions within the cumulant formalism
7 : : !! and, optionally, transport properties within the Kubo formalism.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2008-2026 ABINIT group (JCA, 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 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_cumulant
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_xmpi
28 : : use m_errors
29 : : use m_ebands
30 : : use m_sigmaph
31 : : use m_dtset
32 : : use m_dtfil
33 : : use netcdf
34 : : use m_nctk
35 : :
36 : : !use m_ebands, only : ebands_t
37 : : use defs_abitypes, only : MPI_type
38 : : use m_io_tools, only : open_file, file_exists, is_open
39 : : use m_time, only : cwtime, cwtime_report
40 : : use m_crystal, only : crystal_t
41 : : use m_numeric_tools, only : simpson_cplx, arth, c2r, simpson, safe_div, simpson_int, ctrap, linfit, linspace
42 : : use m_fstrings, only : strcat, sjoin, itoa, ltoa, stoa, ftoa
43 : : !use m_kpts, only : kpts_timrev_from_kptopt
44 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
45 : : use m_fft, only : fourdp
46 : : use m_fftcore, only : ngfft_seq,fftalg_isavailable
47 : : use m_occ, only : occ_fd, occ_dfde
48 : :
49 : : implicit none
50 : :
51 : : private
52 : : !!****
53 : :
54 : : public :: cumulant_driver ! Main entry point for cumulant computations
55 : : !!****
56 : :
57 : : !----------------------------------------------------------------------
58 : :
59 : : !!****t* m_cumulant/cumulant_t
60 : : !! NAME
61 : : !! cumulant_t
62 : : !!
63 : : !! FUNCTION
64 : : !! Container for cumulant and transport quantities
65 : : !!
66 : : !! SOURCE
67 : :
68 : : type,public :: cumulant_t
69 : :
70 : : integer :: my_nspins
71 : : ! Number of spins treated by this MPI rank
72 : :
73 : : integer :: nkcalc
74 : : ! Total number of k-points computed. global variable
75 : :
76 : : integer :: my_nkcalc
77 : : ! Number of k-points treated by this MPI rank i.e. local variable.
78 : :
79 : : integer :: max_nbcalc
80 : : ! Maximum number of bands computed (max over nkcalc and spin). global variable
81 : :
82 : : integer :: nsppol
83 : : ! Number of independent spin polarizations.
84 : :
85 : : integer :: nspinor
86 : : ! Number of spinor components.
87 : :
88 : : integer :: nqbz
89 : : ! Number of q-points in the (dense) BZ for sigma integration
90 : :
91 : : integer :: nqibz
92 : : ! Number of q-points in the (dense) IBZ for sigma integration
93 : :
94 : : integer :: nbsum
95 : : ! Total number of bands used in sum over states without taking into account MPI distribution.
96 : :
97 : : integer :: natom3
98 : : ! 3 * natom.
99 : :
100 : : integer :: nwr
101 : : ! Number of frequency points along the real axis for Sigma(w) and spectral function A(w)
102 : : ! Read from SIGEPH.nc
103 : : ! Odd number so that the mesh is centered on the KS energy.
104 : : ! The spectral function is computed only if nwr > 0 (taken from dtset%nfreqsp)
105 : :
106 : : integer :: nwr_ce
107 : : ! Number of frequency points along the real axis for Sigma(w) and spectral function A(w)
108 : : ! Read from SIGEPH.nc
109 : : ! Odd number so that the mesh is centered on the KS energy.
110 : : ! The spectral function is computed only if nwr > 0 (taken from dtset%nfreqsp)
111 : :
112 : : integer :: ntemp
113 : : ! Number of temperatures.
114 : :
115 : : integer :: comm
116 : : ! MPI communicator
117 : :
118 : : !type(xcomm_t) :: all_comm
119 : :
120 : : type(xcomm_t) :: spin_comm
121 : : ! MPI communicator for parallelism over spins (high-level)
122 : :
123 : : type(xcomm_t) :: kcalc_comm
124 : : ! MPI communicator for parallelism over k-points
125 : :
126 : : type(xcomm_t) :: wt_comm
127 : : ! MPI communicator for time/frequency domain (low-level)
128 : : ! NB: memory is not distributed
129 : :
130 : : type(xcomm_t) :: ncwrite_comm
131 : : ! MPI communicator for parallel netcdf IO used to write results for the different k-points/spins
132 : :
133 : : !integer :: coords(5)
134 : : ! Cartesian coordinates of this processor in the Cartesian grid.
135 : :
136 : : ! real(dp) :: wr_step
137 : : ! Step of the linear mesh along the real axis (Ha units).
138 : :
139 : : integer :: ngqpt(3)
140 : : ! Number of divisions in the Q mesh in the BZ.
141 : :
142 : : integer,allocatable :: kcalc2ebands(:,:)
143 : : ! Mapping ikcalc --> ebands IBZ
144 : : ! Note that this array is not necessarily equation to kcalc2ibz computed in sigmaph
145 : : ! because we may have used sigma_nkpt to downsample the initial nkpt mesh.
146 : : ! This array is computed in get_ebands and is equal to kcalc2ibz if sigma_nkpt == ngkpt
147 : :
148 : : real(dp),allocatable :: linewidths(:,:,:,:,:)
149 : : ! (ntemp, bmin:bmax, nkpt, nsppol, nrta)
150 : : ! Linewidth in the IBZ computed in the SERTA/MRTA.
151 : : ! Non-zero only for the kcalc k-points.
152 : :
153 : : real(dp),allocatable :: vbks(:,:,:,:)
154 : : ! (3, bmin:bmax, nkpt, nsppol))
155 : : ! band velocity in Cartesian coordinates in the IBZ
156 : : ! Non-zero only for the kcalc k-points.
157 : :
158 : : integer :: bmin, bmax, bsize
159 : : ! Only bands between bmin and bmax are considered in the integrals
160 : : ! as we don't compute linewidths for all states.
161 : : ! bmin = minval(%bstart_ks); bmax = maxval(%bstop_ks)
162 : : ! bisze = bmax - bmin + 1
163 : :
164 : : type(ebands_t) :: ebands
165 : : ! bandstructure object used to compute the transport properties
166 : : ! Allocate using only the relevant bands for transport
167 : : ! including valence states to allow to compute different doping
168 : :
169 : : complex(dp) :: ieta
170 : : ! Used to shift the poles in the complex plane (Ha units)
171 : : ! Corresponds to `i eta` term in equations.
172 : :
173 : : real(dp) :: tolcum
174 : :
175 : : integer :: debug
176 : :
177 : : integer,allocatable :: bstart_ks(:,:)
178 : : ! bstart_ks(nkcalc, nsppol)
179 : : ! Initial KS band index included in self-energy matrix elements for each k-point in kcalc.
180 : : ! Depends on spin because all denerate states should be included when symmetries are used.
181 : :
182 : : !integer,allocatable :: bstop_ks(:,:)
183 : : ! bstop_ks(nkcalc, nsppol)
184 : :
185 : : integer,allocatable :: nbcalc_ks(:,:)
186 : : ! nbcalc_ks(nkcalc, nsppol)
187 : : ! Number of bands included in self-energy matrix elements for each k-point in kcalc.
188 : : ! Depends on spin because all denerate states should be included when symmetries are used.
189 : :
190 : : integer,allocatable :: coords_kws(:,:)
191 : : ! (2, self%nsppol))
192 : : ! Cartesian coordinates of this processor in the Cartesian grid.
193 : :
194 : : integer,allocatable :: my_spins(:)
195 : : ! my_spins(my_nspins)
196 : : ! Indirect table giving the spin indices treated by this rank.
197 : : ! Used only the collinear case with nspinor == 1
198 : :
199 : : integer,allocatable :: my_ikcalc(:)
200 : : ! my_ikcalc(my_nkcalc)
201 : : ! List of ikcalc indices treated by this pool if k-point parallelism is activated.
202 : :
203 : : real(dp),allocatable :: time_mesh(:,:,:,:,:)
204 : : ! time_mesh(nwr, max_nbcalc, my_nkcalc, nsppol)
205 : : ! Time mesh in atomic units
206 : :
207 : : real(dp),allocatable :: kcalc(:,:)
208 : : ! kcalc(3, nkcalc)
209 : : ! List of k-points where the self-energy is computed.
210 : : ! This array is not MPI-distributed.
211 : :
212 : : real(dp),allocatable :: kTmesh(:)
213 : : ! kTmesh(ntemp)
214 : : ! List of temperatures (kT units).
215 : :
216 : : real(dp),allocatable :: mu_e(:)
217 : : ! mu_e(ntemp)
218 : : ! chemical potential of electrons for the different temperatures.
219 : :
220 : : real(dp),allocatable :: l0(:,:,:,:,:), l1(:,:,:,:,:), l2(:,:,:,:,:)
221 : : real(dp),allocatable :: l0_dm(:,:,:,:,:), l1_dm(:,:,:,:,:), l2_dm(:,:,:,:,:)
222 : : ! (3, 3, 2, nsppol, ntemp)
223 : : ! Onsager coeficients in Cartesian coordinates
224 : : ! _dm is version with Dyson Migdal
225 : :
226 : : integer,allocatable :: kcalc2ibz(:,:)
227 : : !kcalc2ibz(nkcalc, 6))
228 : : ! Mapping ikcalc --> IBZ as reported by listkk.
229 : :
230 : : real(dp),allocatable :: dw_vals(:,:,:,:)
231 : : ! dw_vals(ntemp, max_nbcalc, my_nkcalc, nsppol)
232 : : ! Debye-Waller term (static).
233 : :
234 : : real(dp),allocatable :: e0vals(:,:,:)
235 : : ! (max_nbcalc, my_nkcalc, nsppol)
236 : : ! KS energies at the calculated states retrieved from EPH output file
237 : :
238 : : real(dp),allocatable :: wrmesh_b(:,:,:,:)
239 : : real(dp),allocatable :: wrmesh_ce(:,:,:,:)
240 : : ! wrmesh_b(nwr, max_nbcalc, my_nkcalc, nsppol))
241 : : ! Frequency mesh along the real axis (Ha units) used for the different bands
242 : : ! Each mesh is **centered** on the corresponding KS energy.
243 : :
244 : : complex(dp),allocatable :: vals_e0ks(:,:,:,:)
245 : : ! vals_e0ks(ntemp, max_nbcalc, my_nkcalc, nsppol))
246 : : ! Sigma_eph(omega=eKS, kT, band, ikcalc, spin).
247 : : ! Fan-Migdal + Debye-Waller
248 : :
249 : : complex(dp),allocatable :: vals_wr(:,:,:,:,:)
250 : : ! vals_wr(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
251 : : ! Sigma_eph(omega, kT, band, ikcalc, spin).
252 : : ! enk_KS corresponds to nwr/2 + 1.
253 : :
254 : : complex(dp),allocatable :: ct_vals(:,:,:,:,:)
255 : : ! ct_vals(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
256 : : ! Cumulant function (time, kT, band, ikcalc, spin).
257 : :
258 : : complex(dp),allocatable :: c1(:,:,:,:,:)
259 : : ! FIXME ct_vals(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
260 : : ! Cumulant function (time, kT, band, ikcalc, spin).
261 : :
262 : : complex(dp),allocatable :: c2(:,:,:,:,:)
263 : : ! FIXME ct_vals(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
264 : : ! Cumulant function (time, kT, band, ikcalc, spin).
265 : :
266 : : complex(dp),allocatable :: c3(:,:,:,:,:)
267 : : ! FIXME ct_vals(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
268 : : ! Cumulant function (time, kT, band, ikcalc, spin).
269 : :
270 : : complex(dp),allocatable :: gt_vals(:,:,:,:,:)
271 : : ! FIXME vals_wr(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
272 : : ! Green's function in time domain (time, kT, band, ikcalc, spin).
273 : :
274 : : complex(dp),allocatable :: gw_vals(:,:,:,:,:)
275 : : ! gw_vals(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
276 : : ! Green's function in frequency domain(omega, kT, band) for given (ikcalc, spin).
277 : :
278 : : real(dp),allocatable :: spfunc_dm_wr(:,:,:,:,:)
279 : : real(dp),allocatable :: spfunc_wr(:,:,:,:,:)
280 : : ! ce_spfunc_wr(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
281 : : ! spectral function for dyson migdal case (omega, kT, band, ikcalc, spin).
282 : :
283 : : ! real(dp),allocatable :: ce_spfunc_wr(:,:,:,:,:)
284 : : ! ce_spfunc_wr(nwr, ntemp, max_nbcalc, my_nkcalc, nsppol)
285 : : ! Absorption spectrum (omega, kT, band, ikcalc, spin).
286 : :
287 : : ! for Cumulant and for Dyson Migdal (_dm)
288 : : real(dp),allocatable :: seebeck(:,:,:,:,:)
289 : : real(dp),allocatable :: seebeck_dm(:,:,:,:,:)
290 : : real(dp),allocatable :: kappa(:,:,:,:,:)
291 : : real(dp),allocatable :: kappa_dm(:,:,:,:,:)
292 : : ! (3, 3, 2, nsppol, ntemp)
293 : : ! Transport coefficients in Cartesian coordinates, for Cumulant and for Dyson Migdal (_dm)
294 : :
295 : : real(dp),allocatable :: conductivity_mu(:,:,:,:,:)
296 : : real(dp),allocatable :: conductivity_mu_dm(:,:,:,:,:)
297 : : ! (3, 3, 2, nsppol, ntemp)
298 : : ! Conductivity in Siemens * cm-1
299 : : ! computed by summing over k-points rather that by performing an energy integration).
300 : : ! for Cumulant and for Dyson Migdal (_dm)
301 : :
302 : : real(dp),allocatable :: print_dfdw(:,:)
303 : :
304 : : real(dp),allocatable :: mobility_mu(:,:,:,:,:)
305 : : real(dp),allocatable :: mobility_mu_dm(:,:,:,:,:)
306 : : ! (3, 3, 2, nsppol, ntemp)
307 : : ! mobility for electrons and holes (third dimension) at transport_mu_e(ntemp)
308 : : ! Third dimension is for electron/hole
309 : : ! for Cumulant and for Dyson Migdal (_dm)
310 : :
311 : : real(dp),allocatable :: transport_mu_e(:)
312 : : ! (%ntemp)
313 : : ! Chemical potential at this carrier concentration and temperature
314 : :
315 : : real(dp) :: transport_fermie
316 : : ! Fermi level specified in the input file when computing the SIGEPH file.
317 : :
318 : : real(dp) :: transport_extrael
319 : : ! Extra electrons per unit cell specified in the input file when computing the SIGEPH file.
320 : :
321 : : real(dp),allocatable :: n_ehst(:,:,:)
322 : : ! (2, %nsppol, %ntemp)
323 : : ! Number of electrons (e) and holes (h) per unit cell
324 : : ! The first dimension is for electrons/holes.
325 : : ! If nsppol == 2, the second dimension is the number of e/h for spin else the total number of e/h summed over spins.
326 : :
327 : : real(dp) :: eph_extrael
328 : : ! Extra electrons per unit cell used to compute SERTA lifetimes in sigmaph.
329 : :
330 : : real(dp) :: eph_fermie
331 : : ! Fermi level specified in the input file when computing the SIGEPH file.
332 : :
333 : : integer :: ce_ngfft(18)
334 : : integer :: ce_ngfft_g(18)
335 : : type(mpi_type) :: ce_mpi_enreg
336 : :
337 : : contains
338 : :
339 : : procedure :: init => cumulant_init
340 : : procedure :: compute => cumulant_compute
341 : : procedure :: kubo_transport => cumulant_kubo_transport
342 : : procedure :: sigmaph_ncread => cumulant_sigmaph_ncread
343 : : procedure :: ncwrite => cumulant_ncwrite
344 : : !procedure :: print_txt_files => cumulant_print_txt_files
345 : : procedure :: free => cumulant_free
346 : :
347 : : end type cumulant_t
348 : : !!***
349 : :
350 : : !----------------------------------------------------------------------
351 : :
352 : : contains !=====================================================
353 : : !!***
354 : :
355 : : !----------------------------------------------------------------------
356 : :
357 : : !!****f* m_cumulant/cumulant_driver
358 : : !! NAME
359 : : !! cumulant_driver
360 : : !!
361 : : !! FUNCTION
362 : : !! General driver to compute transport properties
363 : : !!
364 : : !! INPUTS
365 : : !! dtfil<datafiles_type>=variables related to files.
366 : : !! dtset<dataset_type>=All input variables for this dataset.
367 : : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
368 : : !! cryst<crystal_t>=Crystalline structure
369 : : !! comm=Initial MPI communicator with all procs provided by user.
370 : : !!
371 : : !! SOURCE
372 : :
373 : 2 : subroutine cumulant_driver(dtfil, dtset, ebands, cryst, comm)
374 : :
375 : : !Arguments ------------------------------------
376 : : !scalars
377 : : integer, intent(in) :: comm
378 : : type(datafiles_type),intent(in) :: dtfil
379 : : type(dataset_type),intent(in) :: dtset
380 : : type(crystal_t),intent(in) :: cryst
381 : : type(ebands_t),intent(in) :: ebands
382 : :
383 : : !Local variables ------------------------------
384 : : integer,parameter :: master = 0
385 : : integer :: my_rank
386 : : character(len=fnlen) :: path
387 : 2 : type(cumulant_t) :: cumulant
388 : 2 : type(sigmaph_t) :: sigmaph
389 : : !arrays
390 : : integer :: unts(2)
391 : :
392 : : ! *************************************************************************
393 : :
394 : 2 : my_rank = xmpi_comm_rank(comm)
395 [ + + ]: 6 : unts = [std_out, ab_out]
396 : :
397 : : !sigeph_filepath = dtfil%filsigephin
398 : 2 : call wrtout(unts, ch10//' Entering cumulant expansion computation driver.')
399 : :
400 : : ! Initialize cumulant object
401 : 2 : call cumulant%init(dtset, dtfil, cryst, ebands, comm, sigmaph)
402 : : !if (cumulant%all_comm%me == -1) goto 100
403 : :
404 : : ! Compute C(t), G(t), G(w) and A(w) using cumulant expansion
405 : 2 : call cumulant%compute()
406 : :
407 : : ! Compute transport properties within the Kubo formalism.
408 [ + - + + : 3 : if (any(abs(dtset%sigma_erange) > zero)) call cumulant%kubo_transport(dtset, cryst)!path, ncid)
+ - ]
409 : :
410 : : ! Output netcdf file with cumulant results e.g. A_nk(w).
411 : : ! Use EPH prefix because we may also have EE and PHE cumulants.
412 [ + - ]: 2 : path = strcat(dtfil%filnam_ds(4), "_EPH_CUMULANT.nc")
413 : : !call cumulant%ncwrite(path, cryst, ebands, dtset)
414 : 2 : call cumulant%ncwrite(path, cryst, dtset)
415 : :
416 : : !if (my_rank == master) then
417 : : ! ! Print cumulant expansion results to stdout and other external txt files (for the test suite)
418 : : ! !call cumulant%print_txt_files(cryst, dtset, dtfil)
419 : : !end if
420 : :
421 : : ! Free memory
422 : : !100
423 : 2 : call cumulant%free()
424 : 2 : call sigmaph%free()
425 : :
426 [ - + - + : 2 : end subroutine cumulant_driver
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
+ - - - -
- + - - -
- - - - -
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + ]
427 : : !!***
428 : :
429 : : !----------------------------------------------------------------------
430 : :
431 : : !!****f* m_cumulant/cumulant_init
432 : : !! NAME
433 : : !! cumulant_init
434 : : !!
435 : : !! FUNCTION
436 : : !! Initialization of the cumulant, ex: time mesh
437 : : !!
438 : : !! INPUTS
439 : : !! dtset<dataset_type>=All input variables for this dataset.
440 : : !! dtfil<datafiles_type>=variables related to files.
441 : : !! comm=Initial MPI communicator with all procs provided by user.
442 : : !!
443 : : !! SOURCE
444 : :
445 [ - + - + : 2 : subroutine cumulant_init(self, dtset, dtfil, cryst, ebands, comm, sigmaph )
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - -
- - - + -
- - - - +
- - - - -
- - - ]
446 : :
447 : : !Arguments --------------------------------------
448 : : integer,intent(in) :: comm
449 : : type(dataset_type),intent(in) :: dtset
450 : : type(datafiles_type),intent(in) :: dtfil
451 : : type(crystal_t),intent(in) :: cryst
452 : : type(ebands_t),intent(in) :: ebands
453 : : type(sigmaph_t), intent(out) :: sigmaph
454 : : class(cumulant_t),intent(inout) :: self
455 : :
456 : : !Local variables --------------------------------
457 : : integer, parameter :: master = 0
458 : : integer :: ncerr, ncid, my_rank, nprocs, ierr, spin, ikcalc, ib, color !my_spin, my_kcalc, cnt, nbands,
459 : : integer :: fftalg, fftalga
460 : : real(dp) :: cpu, wall, gflops, rsize
461 : : logical :: is_prime
462 : : character(len=500) :: msg
463 : : character(len=fnlen) :: sigeph_filepath
464 : : integer :: unts(2), facts(2)
465 : 2 : real(dp), allocatable :: rtmp_vals_wr(:,:,:,:,:,:), rtmp_vals_e0ks(:,:,:,:,:)
466 : : !type(ebands_t) :: tmp_ebands
467 : : real(dp) :: extrael_fermie(2),sigma_erange(2)
468 : : integer :: sigma_ngkpt(3)
469 : : #ifdef HAVE_MPI
470 : : integer :: ndims, comm_cart, me_cart
471 : : logical :: reorder
472 : 2 : integer,allocatable :: dims(:)
473 [ - + ]: 2 : logical,allocatable :: periods(:), keepdim(:)
474 : : #endif
475 : :
476 : : !************************************************************************
477 : :
478 : 2 : call cwtime(cpu, wall, gflops, "start")
479 : :
480 : 2 : self%comm = comm
481 : 2 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
482 : :
483 : :
484 [ + + ]: 6 : unts = [std_out, ab_out]
485 : :
486 : : ! Read data from out_SIGEPH.out file
487 : 2 : sigeph_filepath = dtfil%filsigephin
488 : 2 : call wrtout(unts, sjoin("- Reading Sigma results from:", sigeph_filepath), newlines=1, do_flush=.True.)
489 : 2 : call self%sigmaph_ncread(sigeph_filepath, ncid, comm)
490 : : sigmaph = sigmaph_read(sigeph_filepath, dtset, comm, msg, ierr, keep_open = .false., &
491 [ - + - + : 2 : extrael_fermie=extrael_fermie, sigma_ngkpt=sigma_ngkpt, sigma_erange=sigma_erange)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- - - - -
+ - - - -
- + - - -
- - - -
- ]
492 : :
493 [ + + + + : 30 : self%bmin = minval(sigmaph%bstart_ks); self%bmax = maxval(sigmaph%bstop_ks)
+ + + + ]
494 : : !ABI_CALLOC(self%vbks, (3, self%bmin:self%bmax, ebands%nkpt, sigmaph%nsppol))
495 : :
496 : : !if (any(abs(dtset%sigma_erange) > zero)) then
497 : : self%ebands = sigmaph%get_ebands(cryst, ebands, [self%bmin, self%bmax], &
498 [ + + - + : 6 : self%kcalc2ebands, self%linewidths, self%vbks, xmpi_comm_self)
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ ]
499 : : !end if
500 : :
501 : 2 : self%eph_extrael = extrael_fermie(1)
502 : 2 : self%eph_fermie = extrael_fermie(2)
503 : :
504 [ - + + - : 6 : ABI_MALLOC(self%mu_e, (self%ntemp))
+ - + - -
+ - + -
+ ]
505 : : ! Possibility to increase nwr in case of interpolation in cumulant to add extra points
506 : : ! TODO: not working yet
507 : 2 : self%nwr_ce = self%nwr!*4 -1 ! Odd
508 : :
509 : 2 : self%ieta = j_dpc * sigmaph%ieta
510 : : !self%ebands = tmp_ebands
511 [ + - - + : 6 : self%mu_e = sigmaph%mu_e
- - - - -
- + + ]
512 : :
513 : : ! Setting variables to launch 1d FFT calculations later on
514 [ + + ]: 8 : call ngfft_seq(self%ce_ngfft, [self%nwr, 1, 1])
515 [ + + ]: 8 : self%ce_ngfft(4:6) = self%ce_ngfft(1:3)
516 : :
517 : 2 : fftalg = self%ce_ngfft(7); fftalga = fftalg/100
518 [ - + ]: 2 : if (fftalga == FFT_SG) then
519 : 0 : self%ce_ngfft(7)= 102
520 : 0 : ABI_WARNING("Setting fftalg to 102. Please link with FFTW3 or DFTI for better performance!")
521 : : end if
522 : : !self%ce_ngfft(7)= 102
523 : :
524 : 2 : call initmpi_seq(self%ce_mpi_enreg)
525 : 2 : call self%ce_mpi_enreg%distribfft%init_seq('c', self%ce_ngfft(2), self%ce_ngfft(3), 'all')
526 : 2 : call self%ce_mpi_enreg%distribfft%init_seq('f', self%ce_ngfft(2), self%ce_ngfft(3), 'all')
527 : :
528 [ + + ]: 8 : call ngfft_seq(self%ce_ngfft_g, [self%nwr_ce, 1, 1])
529 [ + + ]: 8 : self%ce_ngfft_g(4:6) = self%ce_ngfft_g(1:3)
530 : :
531 : 2 : fftalg = self%ce_ngfft_g(7); fftalga = fftalg/100
532 [ - + ]: 2 : if (fftalga == FFT_SG) then
533 : 0 : self%ce_ngfft_g(7)= 102
534 : 0 : ABI_WARNING("Setting fftalg to 102. Please link with FFTW3 or DFTI for better performance!")
535 : : end if
536 : : !self%ce_ngfft_g(7)= 102
537 : :
538 : : ! Setting debugging ( higher verbosity )
539 : 2 : self%tolcum = dtset%tolcum
540 [ + - ]: 2 : if (self%tolcum > zero) then
541 : 2 : self%debug = 0
542 : : else
543 : 0 : self%debug = 1
544 : : endif
545 : 2 : self%tolcum = abs(self%tolcum)
546 : :
547 : 2 : call wrtout(unts, " Cumulant parameters:")
548 : 2 : call wrtout(unts, sjoin(" Number of spins: ", itoa(self%nsppol)))
549 : 2 : call wrtout(unts, sjoin(" Number of spinor components: ", itoa(self%nspinor)))
550 : 2 : call wrtout(unts, sjoin(" Number of k-points computed: ", itoa(self%nkcalc)))
551 : 2 : call wrtout(unts, sjoin(" Maximum number of bands computed: ", itoa(self%max_nbcalc)))
552 : 2 : call wrtout(unts, sjoin(" Number of frequencies in Sigma(w): ", itoa(self%nwr)))
553 : 2 : call wrtout(unts, sjoin(" Number of frequencies in G(w): ", itoa(self%nwr_ce)))
554 : 2 : call wrtout(unts, sjoin(" Number of Temperatures: ", itoa(self%ntemp)))
555 : :
556 : : ! ========================
557 : : ! === MPI DISTRIBUTION ===
558 : : ! ========================
559 : : ! At present, there are two levels of MPI parallelism: kcalc and time/frequencies.
560 : : ! Only the kcalc parallelism leads to memory distribution (my_nkcalc) whereas the secon level is only
561 : : ! used to distribute the loop over frequencies/time points.
562 : : ! Note the following:
563 : : !
564 : : ! 1) We distribute contigous blocks of kcalc k-points so that we can read/write data in mykcalc chuncks.
565 : :
566 : : ! TODO:
567 : : ! *) Spin and nkcalc(nsppol)
568 : : ! *) Decide whther it makes sense to create a new comm if nprocs > nkcalc and mod(nprocs, nkcalc) /= 0
569 : : ! *) No need anymore for parallelism in time/frequency
570 : :
571 : :
572 : 2 : self%kcalc_comm%nproc = nprocs
573 : : ! self%wt_comm%nproc = nprocs / self%kcalc_comm%nproc
574 : :
575 [ - + ]: 2 : if (nprocs > self%nkcalc) then
576 : 0 : call ifact2(nprocs, facts, is_prime)
577 [ # # ]: 0 : if (is_prime) then
578 : : ! TODO
579 : 0 : ABI_ERROR("prime numbers are evil!")
580 : : end if
581 [ # # ]: 0 : if (facts(1) > self%nkcalc) then
582 : : ! Avoid procs with no k-points.
583 [ # # # # ]: 0 : facts(1:2) = facts(2:1:-1)
584 : : end if
585 : 0 : self%kcalc_comm%nproc = facts(1)
586 : : ! self%wt_comm%nproc = facts(2)
587 : : !do cnt=self%nkcalc,1,-1
588 : : ! if (mod(nprocs, cnt) == 0 .and. mod(self%nkcalc, cnt) == 0) then
589 : : ! self%kcalc_comm%nproc = cnt; self%my_nkcalc = self%nkcalc / cnt; exit
590 : : ! end if
591 : : !end do
592 : : end if
593 : :
594 [ - + + - : 12 : ABI_ICALLOC(self%coords_kws, (2, self%nsppol))
+ - + - -
+ - + - +
+ + + + ]
595 : : #ifdef HAVE_MPI
596 : : ! Create 2d cartesian communicator: kpoints in Sigma_k, (time/frequency)
597 : : ! FIXME: Fix spin
598 : 2 : spin = 1
599 : 2 : ndims = 1
600 [ - + ]: 2 : ABI_MALLOC(dims, (ndims))
601 [ - + ]: 2 : ABI_MALLOC(periods, (ndims))
602 [ - + ]: 2 : ABI_MALLOC(keepdim, (ndims))
603 [ + + ]: 4 : periods(:) = .False.; reorder = .False.
604 : : ! dims = [self%kcalc_comm%nproc, self%wt_comm%nproc]
605 [ - + + + ]: 6 : dims = [self%kcalc_comm%nproc]
606 : :
607 : 2 : call MPI_CART_CREATE(comm, ndims, dims, periods, reorder, comm_cart, ierr)
608 : : ! Find the index and coordinates of the current processor
609 : 2 : call MPI_COMM_RANK(comm_cart, me_cart, ierr)
610 : 2 : call MPI_CART_COORDS(comm_cart, me_cart, ndims, self%coords_kws(:, spin), ierr)
611 : :
612 : : ! Create communicator for kpoints.
613 : 2 : keepdim = .False.; keepdim(1) = .True.
614 : 2 : call MPI_CART_SUB(comm_cart, keepdim, self%kcalc_comm%value, ierr); self%kcalc_comm%me = xmpi_comm_rank(self%kcalc_comm%value)
615 : :
616 : : ! Create communicator for w/t.
617 : : ! keepdim = .False.; keepdim(2) = .True.
618 : : ! call MPI_CART_SUB(comm_cart, keepdim, self%wt_comm%value, ierr); self%wt_comm%me = xmpi_comm_rank(self%wt_comm%value)
619 : :
620 : : ! Create communicator for spins.
621 : : !keepdim = .False.; keepdim(5) = .True.
622 : : !call MPI_CART_SUB(comm_cart, keepdim, self%spin_comm%value, ierr); self%spin_comm%me = xmpi_comm_rank(self%spin_comm%value)
623 : :
624 : : ! Create communicator for parallel IO
625 : : ! If parallelism over time/frequencies is ON, we select only a column of procs (kcalc dimension)
626 : : ! to avoid writing the same data multiple times as omega/time are not MPI distributed
627 : : ! Obviously I'm assuming HDF5 + MPI-IO
628 [ + - ]: 2 : color = xmpi_undefined; if (self%coords_kws(2, spin) == 0) color = 1
629 : 2 : call xmpi_comm_split(comm, color, my_rank, self%ncwrite_comm%value, ierr)
630 [ + - ]: 2 : if (color == 1) then
631 : 2 : self%ncwrite_comm%me = xmpi_comm_rank(self%ncwrite_comm%value)
632 : 2 : self%ncwrite_comm%nproc = xmpi_comm_size(self%ncwrite_comm%value)
633 : : else
634 : 0 : call self%ncwrite_comm%set_to_null()
635 : : end if
636 : :
637 [ - + ]: 2 : ABI_FREE(dims)
638 [ - + ]: 2 : ABI_FREE(periods)
639 : 2 : ABI_FREE(keepdim)
640 : 2 : call xmpi_comm_free(comm_cart)
641 : : #endif
642 : :
643 : : call wrtout(unts, &
644 : 2 : sjoin("- Using: ", itoa(self%kcalc_comm%nproc), "MPI procs to distribute:", itoa(self%nkcalc), "k-points."))
645 : : !call wrtout(unts, &
646 : : ! sjoin("- Using: ", itoa(self%wt_comm%nproc), "MPI procs to parallelize:", itoa(self%nwr_ce), "frequency/time loops."))
647 : 2 : call wrtout(std_out, sjoin("- Performing parallel IO with:", itoa(self%ncwrite_comm%nproc), "procs"))
648 : :
649 [ - + - - ]: 2 : if (self%ncwrite_comm%nproc /= 1 .and. .not. nctk_has_mpiio) then
650 : 0 : ABI_WARNING("Cumulant requires netcdf with MPI-IO support! Continuing anyway but runtime error is expected!")
651 : : end if
652 : :
653 : : ! Consistency check.
654 : : !if (self%kcalc_comm%nproc * self%wt_comm%nproc /= nprocs) then
655 [ - + ]: 2 : if (self%kcalc_comm%nproc /= nprocs) then
656 : : write(msg, "(a,i0,3a, 3(a,1x,i0))") &
657 : 0 : "Cannot create 2d Cartesian grid with total nprocs: ", nprocs, ch10, &
658 : 0 : "Idle processes are not supported. The product of comm%nproc should be equal to nprocs.", ch10, &
659 : 0 : "kcalc_nproc (", self%kcalc_comm%nproc, ") != ", & ! x wt_nproc (", self%wt_comm%nproc, ") != ", &
660 : 0 : self%kcalc_comm%nproc ! * self%wt_comm%nproc
661 : 0 : ABI_ERROR(msg)
662 : : end if
663 : :
664 : : ! Here we distribute the k-points inside kcalc_comm.
665 [ - + ]: 2 : call xmpi_split_block(self%nkcalc, self%kcalc_comm%value, self%my_nkcalc, self%my_ikcalc)
666 [ - + ]: 2 : ABI_CHECK(self%my_nkcalc > 0, sjoin("nkcalc (", itoa(self%nkcalc), ") > nprocs (", itoa(nprocs), ")"))
667 : :
668 : 2 : ib = mod(self%nkcalc, self%my_nkcalc)
669 : 2 : call xmpi_sum(ib, self%comm, ierr)
670 [ - + ]: 2 : if (ib /= 0) then
671 : 0 : ABI_COMMENT("The number of MPI procs should be divisible by nkcalc to distribute memory equally!")
672 : : end if
673 : :
674 : : ! Distribute spins and create mapping to spin index.
675 [ - + ]: 2 : if (self%nsppol == 2) then
676 : 0 : ABI_ERROR("cumulant with nsppol 2 is not supported.")
677 [ # # ]: 0 : call xmpi_split_block(self%nsppol, comm, self%my_nspins, self%my_spins)
678 : : !ABI_CHECK(self%my_nspins > 0, sjoin("nsppol (", itoa(self%nsppol), ") > spin_comm_nproc (", itoa(comm), ")"))
679 : : else
680 : : ! No nsppol parallelism DOH!
681 : 2 : self%my_nspins = 1
682 [ - + - + ]: 2 : ABI_MALLOC(self%my_spins, (self%my_nspins))
683 [ + + ]: 4 : self%my_spins = 1
684 : : end if
685 : :
686 : : !ABI_CALLOC(self%ce_spfunc_wr, (self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
687 [ - + + - : 5508 : ABI_CALLOC(self%spfunc_dm_wr, (self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ + + + +
+ + + + +
+ ]
688 [ - + + - : 5508 : ABI_CALLOC(self%spfunc_wr, (self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ + + + +
+ + + + +
+ ]
689 [ - + + - : 5508 : ABI_CALLOC(self%gw_vals, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ + + + +
+ + + + +
+ ]
690 : :
691 [ - + ]: 2 : if (self%debug == 1) then
692 [ # # # # : 0 : ABI_CALLOC(self%time_mesh, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
693 [ # # # # : 0 : ABI_CALLOC(self%ct_vals, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
694 [ # # # # : 0 : ABI_CALLOC(self%c1, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
695 [ # # # # : 0 : ABI_CALLOC(self%c2, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
696 [ # # # # : 0 : ABI_CALLOC(self%c3, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
697 [ # # # # : 0 : ABI_CALLOC(self%gt_vals, (self%nwr_ce, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
698 : : endif
699 : :
700 [ - + + - : 14 : ABI_MALLOC(rtmp_vals_wr, (2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ ]
701 [ - + + - : 12 : ABI_MALLOC(rtmp_vals_e0ks, (2, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
702 [ - + + - : 14 : ABI_MALLOC(self%vals_wr, (self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + -
+ ]
703 [ - + + - : 12 : ABI_MALLOC(self%vals_e0ks, (self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
704 [ - + + - : 12 : ABI_MALLOC(self%wrmesh_b, (self%nwr, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
705 [ - + + - : 12 : ABI_MALLOC(self%wrmesh_ce, (self%nwr_ce, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
706 : : ! ABI_MALLOC(self%wrmesh_ce, (self%nwr_ce, self%max_nbcalc, self%my_nkcalc, self%nsppol))
707 [ - + + - : 10 : ABI_MALLOC(self%e0vals, (self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - + -
+ - - + -
+ - + ]
708 [ - + + - : 12 : ABI_MALLOC(self%dw_vals, (self%ntemp, self%max_nbcalc, self%my_nkcalc, self%nsppol))
- + + - -
+ + - - +
+ - + - +
- - + - +
- + ]
709 : :
710 : : ! Estimate memory (only nwr arrays)
711 [ - + ]: 2 : rsize = five; if (self%debug == 1) rsize = rsize + 3
712 : 2 : rsize = rsize * self%nwr_ce * self%ntemp * self%max_nbcalc * self%my_nkcalc * self%nsppol
713 : 2 : call xmpi_max(rsize, self%comm, ierr)
714 : 2 : write(msg,'(a,f8.1,a)')' Memory needed for cumulant arrays per MPI proc: ',dp * rsize * b2Mb, ' [Mb] <<< MEM'
715 : 2 : call wrtout(std_out, msg)
716 : :
717 : : ! Each MPI proc reads a chunk of my_nkcalc entries starting at (ikcalc, spin) from the SIGEPH.nc file.
718 : : ! so that memory is MPI-distributed.
719 : : ! On disk, we have the following netcdf arrays:
720 : : !
721 : : ! nctkarr_t("wrmesh_b", "dp", "nwr, max_nbcalc, nkcalc, nsppol"), &
722 : : ! nctkarr_t("vals_wr", "dp", "two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
723 : : ! nctkarr_t("vals_e0ks", "dp", "two, ntemp, max_nbcalc, nkcalc, nsppol"), &
724 : : ! nctkarr_t("ks_enes", "dp", "max_nbcalc, nkcalc, nsppol"), &
725 : : !
726 : : ! Note the nkcalc global dimension instead of my_nkcalc.
727 : : ! Results are in a.u.
728 : :
729 : : ! Offse of this proc
730 : 2 : spin = self%my_spins(1)
731 : 2 : ikcalc = self%my_ikcalc(1)
732 : :
733 : : ncerr = nf90_get_var(ncid, vid("vals_wr"), rtmp_vals_wr, start=[1,1,1,1,ikcalc,spin], &
734 [ + + + + ]: 26 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
735 [ - + ]: 2 : NCF_CHECK(ncerr)
736 : :
737 : : ncerr = nf90_get_var(ncid, vid("wrmesh_b"), self%wrmesh_b, start=[1,1,ikcalc,spin], &
738 [ + + + + ]: 18 : count=[self%nwr, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
739 [ - + ]: 2 : NCF_CHECK(ncerr)
740 : :
741 : : ncerr = nf90_get_var(ncid, vid("vals_e0ks"), rtmp_vals_e0ks, start=[1,1,1,ikcalc,spin], &
742 [ + + + + ]: 22 : count=[2, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
743 [ - + ]: 2 : NCF_CHECK(ncerr)
744 : :
745 : : ncerr = nf90_get_var(ncid, vid("ks_enes"), self%e0vals, start=[1,ikcalc,spin], &
746 [ + + + + ]: 14 : count=[self%max_nbcalc, self%my_nkcalc, self%my_nspins])
747 [ - + ]: 2 : NCF_CHECK(ncerr)
748 : :
749 : : ncerr = nf90_get_var(ncid, vid("dw_vals"), self%dw_vals, start=[1,1,ikcalc,spin], &
750 [ + + + + ]: 18 : count=[self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
751 [ - + ]: 2 : NCF_CHECK(ncerr)
752 : :
753 : : ! Close the file here but the Kubo equation needs to read v_nk from file.
754 : : ! We will have to rationalize this part.
755 [ - + ]: 2 : NCF_CHECK(nf90_close(ncid))
756 : :
757 : : ! Convert from real arrays to Fortran complex
758 [ + - + - : 5498 : self%vals_wr = (rtmp_vals_wr(1,:,:,:,:,:) + j_dpc * rtmp_vals_wr(2,:,:,:,:,:))!*Ha_eV
+ - + - +
- - + - -
+ + + + +
+ + + +
+ ]
759 [ + - + - : 98 : self%vals_e0ks = (rtmp_vals_e0ks(1,:,:,:,:) + j_dpc * rtmp_vals_e0ks(2,:,:,:,:))!*Ha_eV
+ - + - -
+ - - + +
+ + + + +
+ ]
760 : :
761 : : ! self%wrmesh_b = self%wrmesh_b*Ha_eV
762 : : ! self%e0vals = self%e0vals*Ha_eV
763 : : ! self%max_nbcalc=self%max_nbcalc*Ha_eV
764 : :
765 : : ! self%e0vals = self%e0vals * Ha_eV
766 : : ! self%wrmesh_b = self%wrmesh_b * Ha_eV
767 : :
768 [ + - ]: 2 : ABI_SFREE(rtmp_vals_wr)
769 [ + - ]: 2 : ABI_SFREE(rtmp_vals_e0ks)
770 : :
771 [ - + - + ]: 12 : call cwtime_report(" cumulant init", cpu, wall, gflops)
772 : :
773 : : contains
774 : 10 : integer function vid(vname)
775 : : character(len=*),intent(in) :: vname
776 : 10 : vid = nctk_idname(ncid, vname)
777 : : end function vid
778 : :
779 : : end subroutine cumulant_init
780 : : !!***
781 : :
782 : : !----------------------------------------------------------------------
783 : :
784 : : !!****f* m_cumulant/cumulant_compute
785 : : !! NAME
786 : : !! cumulant_compute
787 : : !!
788 : : !! FUNCTION
789 : : !! Compute cumulant.
790 : : !!
791 : : !! INPUTS
792 : : !!
793 : : !! SOURCE
794 : :
795 : 2 : subroutine cumulant_compute(self)
796 : :
797 : : !Arguments ------------------------------------
798 : : class(cumulant_t),intent(inout) :: self
799 : :
800 : : !Local variables ------------------------------
801 : : integer,parameter :: master = 0
802 : : integer :: nbands, ib, ikcalc, it, iw, spin, itemp, comm
803 : : integer :: nwr, nwr_ce
804 : : integer :: my_rank, nprocs, my_spin, my_ik
805 : : real(dp) :: time, omega, init_t, time_step, time_step_init, time_max!, wmesh_step_init
806 : : real(dp) :: cpu, wall, gflops, cpu_kloop, wall_kloop, gflops_kloop
807 : : real(dp) :: wr_step, Ha_fs ! wr_step_ce
808 : : character(len=500) :: msg
809 : : !arrays
810 : 2 : real(dp),allocatable :: temp_g(:,:,:), temp_r(:,:), temp_r_cplx(:,:), temp_g_ce(:,:,:)
811 : 2 : real(dp),allocatable :: betaoverw2(:) !, dfft(:)
812 : : !complex(dp),allocatable :: temp_reflex(:) ! betaoverw2c(:),
813 : 2 : real(dp),allocatable :: wrmesh_shifted(:), wrmesh_shifted_ce(:), beta(:), c3(:)
814 : 2 : real(dp),allocatable :: time_mesh(:), time_mesh_temp(:)
815 : : real(dp) :: output_c3!, output_test2r, output_test2i
816 : : real(dp) :: m_fit_re, b_fit_re, m_fit_im, b_fit_im, res_re, res_im
817 : 2 : complex(dp),allocatable :: c1(:), ct_temp(:), c_temp(:)
818 : 2 : complex(dp),allocatable :: c2(:), ct(:), gt(:), gw(:), g1(:)
819 : : integer :: fftalg, fftalga
820 : : logical :: use_fft
821 : :
822 : : !************************************************************************
823 : :
824 : : ! Initialization of parallel variables
825 : 2 : comm = self%comm
826 : 2 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
827 : :
828 : 2 : call wrtout(std_out, " Computing cumulant. This may take some time depending on the dims of the problem...")
829 : 2 : call cwtime(cpu_kloop, wall_kloop, gflops_kloop, "start")
830 : :
831 : : ! Initialization of variables and allocation of arrays
832 : 2 : nwr = self%nwr
833 : 2 : nwr_ce = self%nwr_ce
834 [ - + + - : 276 : ABI_CALLOC(c1, (nwr))
+ - + - -
+ - + +
+ ]
835 [ - + + - : 276 : ABI_CALLOC(c_temp, (nwr_ce))
+ - + - -
+ - + +
+ ]
836 : : !ABI_CALLOC(temp_reflex, (nwr/2))
837 [ - + + - : 274 : ABI_CALLOC(c2, (nwr))
- + - + +
+ ]
838 [ - + + - : 276 : ABI_CALLOC(c3, (nwr))
+ - + - -
+ - + +
+ ]
839 [ - + + - : 274 : ABI_CALLOC(g1, (nwr_ce))
- + - + +
+ ]
840 [ - + + - : 274 : ABI_CALLOC(gw, (nwr_ce))
- + - + +
+ ]
841 : : !ABI_CALLOC(dfft, (nwr_ce))
842 [ - + + - : 276 : ABI_CALLOC(time_mesh, (nwr_ce))
+ - + - -
+ - + +
+ ]
843 [ - + + - : 274 : ABI_CALLOC(time_mesh_temp, (nwr))
- + - + +
+ ]
844 [ - + + - : 274 : ABI_CALLOC(ct, (nwr_ce))
- + - + +
+ ]
845 [ - + + - : 274 : ABI_CALLOC(ct_temp, (nwr))
- + - + +
+ ]
846 [ - + + - : 274 : ABI_CALLOC(gt, (nwr_ce))
- + - + +
+ ]
847 [ - + + - : 818 : ABI_CALLOC(temp_g, (2,nwr,1))
+ - + - -
+ - + + +
+ + + + ]
848 [ - + + - : 818 : ABI_CALLOC(temp_g_ce, (2,nwr_ce,1))
+ - + - -
+ - + + +
+ + + + ]
849 : :
850 [ - + + - : 276 : ABI_CALLOC(temp_r, (nwr,1))
- + - + +
+ + + ]
851 [ + - - + : 548 : ABI_CALLOC(temp_r_cplx, (2*nwr_ce,1))
+ - + - +
- - + - +
+ + + + ]
852 : : !ABI_CALLOC(inv_wrmesh_shifted_sq, (nwr_ce))
853 [ - + + - : 274 : ABI_CALLOC(wrmesh_shifted, (nwr))
- + - + +
+ ]
854 [ - + + - : 274 : ABI_CALLOC(wrmesh_shifted_ce, (nwr_ce))
- + - + +
+ ]
855 [ - + + - : 274 : ABI_CALLOC(beta, (nwr))
- + - + +
+ ]
856 [ - + + - : 274 : ABI_CALLOC(betaoverw2, (nwr))
- + - + -
+ + + ]
857 : : !ABI_CALLOC(betaoverw2c, (nwr))
858 : 2 : Ha_fs = 8.955433106
859 : :
860 : : ! Setting direct domain after fft
861 : : ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= gc
862 : : ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index ig
863 : :
864 : : !do iw=1,nwr_ce
865 : : ! dfft(iw) = ig2gfft(iw,nwr)
866 : : !end do
867 : : !ABI_FREE(dfft)
868 : :
869 : 2 : fftalg = self%ce_ngfft(7)
870 : 2 : fftalga = fftalg/100
871 : : ! Loops are MPI-parallelized over k-points and spin.
872 [ + + ]: 4 : do my_spin=1,self%my_nspins
873 : 2 : spin = self%my_spins(my_spin)
874 [ + + ]: 16 : do my_ik=1,self%my_nkcalc
875 : 12 : ikcalc= self%my_ikcalc(my_ik)
876 : : !if (ikcalc /= 1) cycle ! For debugging purpose
877 : 12 : nbands = self%nbcalc_ks(ikcalc, spin)
878 : 12 : call cwtime(cpu, wall, gflops, "start")
879 : :
880 [ + + ]: 32 : do ib=1,nbands
881 : : !if (ib > 1) cycle ! MG HACK To be able to run tests quickly.
882 : :
883 : : ! Shifting w-mesh to center at KS-energy and remove the value w=0.0 from integrations ( Principal Value )
884 : 20 : wr_step = self%wrmesh_b(2,ib,my_ik,spin) - self%wrmesh_b(1,ib,my_ik,spin)
885 : : wrmesh_shifted(:) = self%wrmesh_b(:,ib,my_ik,spin) - self%e0vals(ib,my_ik,spin) &
886 [ + + ]: 2720 : - 0.5 * wr_step
887 : :
888 : :
889 [ + + ]: 52 : do itemp=1,self%ntemp
890 : : ! if (itemp > 1) cycle ! MG HACK To be able to run tests quickly.
891 : :
892 : : ! Defining time mesh
893 : : ! all variables with nwr_ce size are to use after interpolation of the
894 : : ! cumulant function in case we want to add more points
895 [ + - ]: 20 : if (abs(aimag(self%vals_e0ks(itemp, ib, my_ik, spin))) > 1.0e-16) then
896 : 20 : time_max = log(self%tolcum) / (-abs(aimag(self%vals_e0ks(itemp, ib, my_ik, spin))))
897 : : else
898 : 0 : time_max = log(self%tolcum) / (-1.0e-16)
899 : 0 : ABI_WARNING("Imaginary part of the self energy is 0 in cumulant calculation. Reset to 1.e-16")
900 : : endif
901 : 20 : init_t = 0.0
902 : 20 : time_step_init = 1.0/(nwr - 1.0)
903 : 20 : time_mesh_temp(:) = arth(init_t, time_step_init, nwr)
904 [ + + ]: 2720 : time_mesh_temp = time_mesh_temp * 2*PI/wr_step
905 : 20 : time_step = time_mesh_temp(2) - time_mesh_temp(1)
906 : :
907 [ + + ]: 2720 : time_mesh(:) = time_mesh_temp(:) !! arth(init_t, time_step_init, nwr_ce)
908 [ + + ]: 2720 : time_mesh = time_mesh * 2*PI/wr_step
909 : :
910 : : ! Defining frequency mesh for after interpolation
911 : :
912 : : !! TODO: Fix problems from when adding more points to time mesh by
913 : : !interpolation
914 : :
915 : : !!wmesh_step_init = 1.0/(nwr_ce/2.0 -1 )
916 : : !!init_w = -1.0
917 : : !!end_w = 0.0
918 : : !!wrmesh_shifted_ce(:nwr_ce/2+1) = linspace(init_w,end_w, nwr_ce/2+1)
919 : : !!init_w = 0.0
920 : : !!end_w = 1.0
921 : : !!wrmesh_shifted_ce(nwr_ce/2+1:) = linspace(init_w,end_w, nwr_ce/2+1)
922 : :
923 : : !!wr_step_ce = wrmesh_shifted_ce(2) - wrmesh_shifted_ce(1)
924 : : !!wrmesh_shifted_ce(:) = wrmesh_shifted_ce(:) * PI/time_step - 0.5 *wr_step_ce
925 : : !!self%wrmesh_ce(:,ib,my_ik,spin) = wrmesh_shifted_ce(:) + self%e0vals(ib,my_ik,spin) + 0.5 *wr_step_ce
926 : : !!wrmesh_shifted_ce = wrmesh_shifted
927 : :
928 [ - + - - ]: 20 : if (self%debug == 1) self%time_mesh(:,itemp,ib,my_ik,spin) = time_mesh_temp(:)
929 : :
930 [ - + ]: 20 : if (time_mesh_temp(nwr) < time_max) ABI_WARNING(sjoin("KBT",itoa(my_ik),itoa(ib),itoa(itemp)))
931 : : msg = sjoin( &
932 : : "Time mesh smaller than needed to reach tolerance value. Actual value of", &
933 [ + - ]: 20 : ftoa(time_mesh_temp(nwr)), ", desirable ",ftoa(time_max),". Increase nfreqsp.")
934 [ - + ]: 20 : if (time_mesh_temp(nwr) < time_max) ABI_WARNING(msg)
935 [ + + ]: 2720 : beta(:) = abs( aimag( self%vals_wr(:, itemp, ib, my_ik, spin ) ) ) / pi
936 : :
937 : : ! Calculation of the different terms of the cumulant ( c1, c2 and c3 ),
938 [ + + ]: 2720 : betaoverw2(:) = beta(:) / (wrmesh_shifted(:) ** 2)
939 : :
940 [ + + ]: 2720 : c2(:) = -j_dpc * real( self%vals_e0ks( itemp, ib, my_ik, spin ) ) * time_mesh_temp(:)
941 : :
942 : 20 : output_c3 = simpson(wr_step, betaoverw2)
943 [ + + ]: 2720 : c3(:) = -1.0 * output_c3
944 : :
945 : :
946 [ + + ]: 2720 : temp_r(:,1) = betaoverw2(:)
947 : :
948 : : !@Joao: I got different results when I set use_fft to .False.
949 : : !Can you recheck this part, if use_fft = .False. still needed?
950 : :
951 : 20 : use_fft = .True.
952 : : !use_fft = .not. fftalga == FFT_SG
953 : : if (use_fft) then
954 : 20 : call fourdp(1, temp_g, temp_r, -1, self%ce_mpi_enreg, nwr, 1, self%ce_ngfft , 0 )
955 [ + + ]: 2720 : c1(:) = temp_g(1, :, 1) + j_dpc* temp_g(2, :, 1)
956 [ + + ]: 2720 : c1 = c1 * wr_step * nwr
957 [ + + ]: 1360 : c1(2::2) = -1.0 * c1(2::2)
958 : : else
959 : : msg = sjoin("FFT not available, using DFT instead.", &
960 : : " Slower but reliable.", &
961 : : " Parallelism over frequency for integration of C(t)")
962 : : ABI_COMMENT(msg)
963 : : do it=1, nwr
964 : : ! if (mod(it, self%wt_comm%nproc) /= self%wt_comm%me) cycle ! MPI parallelism over time
965 : : time = time_mesh_temp(it)
966 : :
967 : : ! Discrete Fourier Transform of C1 ( cumulant function without the +iwt and -1 parts
968 : : c_temp(:) = betaoverw2(:) * exp( - j_dpc * time * wrmesh_shifted(:) )
969 : : c1(it) = simpson_cplx( nwr, wr_step, c_temp)
970 : : enddo
971 : :
972 : : endif
973 : :
974 : :
975 : : ! The cumulant function as sum of the different components
976 [ + + ]: 2720 : ct_temp(:) = c1 + c2 + c3
977 : : ! Fitting the cumulant function
978 [ + - + + : 2720 : res_re = linfit(nwr,time_mesh_temp(:), real(ct_temp(:)), m_fit_re, b_fit_re)
- + ]
979 [ + + - + ]: 2720 : res_im = linfit(nwr,time_mesh_temp(:), aimag(ct_temp(:)), m_fit_im, b_fit_im)
980 : : ! call xmpi_sum(ct, self%wt_comm%value, ierr)
981 [ - + ]: 20 : if (self%debug == 1) then ! .and. self%wt_comm%nproc > 1) then
982 [ # # ]: 0 : self%ct_vals(:, itemp, ib, my_ik, spin) = ct(:)
983 [ # # ]: 0 : self%c1(:, itemp, ib, my_ik, spin) = c1(:)
984 [ # # ]: 0 : self%c2(:, itemp, ib, my_ik, spin) = c2(:)
985 [ # # ]: 0 : self%c3(:, itemp, ib, my_ik, spin) = c3(:)
986 : : endif
987 : :
988 : : ! Adding extra interpolated points to the cumulant function
989 : : !!ct(:nwr/2) = ct_temp(:nwr/2)
990 : : !!ct(nwr/2:) = time_mesh(nwr/2:) * m_fit_re + b_fit_re + j_dpc * ( time_mesh(nwr/2:) * m_fit_im + b_fit_im )
991 : :
992 : : !!do iw=1, nwr/2
993 : : !! temp_reflex(nwr/2-iw) = - ct(iw) + time_mesh(iw) * m_fit_re + b_fit_re + j_dpc * ( time_mesh(iw) * m_fit_im + b_fit_im )
994 : : !!enddo
995 : : !!ct(nwr_ce - nwr/2:) = ct(nwr_ce - nwr/2:) + temp_reflex(:)
996 [ + + ]: 2720 : ct(:) = ct_temp(:)
997 : : ! Retarded Green's function in time domain
998 [ - + + + ]: 2740 : gt = - j_dpc * exp( ct )
999 : :
1000 : : ! Collect data if wt parallelism.
1001 : : !call xmpi_sum(gt, self%wt_comm%value, ierr)
1002 : :
1003 [ - + ]: 20 : if (self%debug == 1) then ! .and. self%wt_comm%nproc > 1) then
1004 [ # # ]: 0 : self%gt_vals(:, itemp, ib, my_ik, spin) = gt(:)
1005 : : end if
1006 : :
1007 : 20 : use_fft = .True.
1008 : : !use_fft = .not. fftalga == FFT_SG
1009 : 20 : if (use_fft) then
1010 : :
1011 : : ! Fast Fourier Transform to obtain the Green's function in frequency
1012 : : ! domain
1013 [ + + ]: 2720 : temp_g_ce(1,:,1) = real(gt(:))
1014 [ + + ]: 2720 : temp_g_ce(2,:,1) = aimag(gt(:))
1015 : 20 : call fourdp(2, temp_g_ce, temp_r_cplx, 1, self%ce_mpi_enreg, nwr_ce, 1, self%ce_ngfft_g , 0 )
1016 [ + + ]: 2720 : gw(:) = temp_r_cplx(1::2,1) + j_dpc* temp_r_cplx(2::2,1)
1017 : :
1018 : :
1019 : : ! TODO use ig2gfft from 52_fft_mpi_noabirule/m_fftcore.F90 instead of the two following lines
1020 [ - + ]: 20 : if ( mod(nwr_ce,2) == 0 ) then
1021 [ # # ]: 0 : self%gw_vals(1:int(nwr_ce/2.0), itemp, ib, my_ik, spin) = gw(int(nwr_ce/2.0):nwr_ce)
1022 [ # # ]: 0 : self%gw_vals(int(nwr_ce/2.0)+1:nwr_ce, itemp, ib, my_ik, spin) = gw(1:int(nwr_ce/2.0)-1)
1023 : : else
1024 [ + + ]: 1380 : self%gw_vals(1:int(nwr_ce/2.0)+1, itemp, ib, my_ik, spin) = gw(int(nwr_ce/2.0)+1:nwr_ce)
1025 [ + + ]: 1360 : self%gw_vals(int(nwr_ce/2.0)+2:nwr_ce, itemp, ib, my_ik, spin) = gw(1:int(nwr_ce/2.0))
1026 : : endif
1027 : :
1028 : :
1029 [ + + ]: 2720 : self%gw_vals(:, itemp, ib, my_ik, spin) = self%gw_vals(:, itemp, ib, my_ik, spin) * time_step
1030 : :
1031 : : ! FFT is different from integration methods and the two extreme points
1032 : : ! need to be compensated
1033 : : self%gw_vals(:, itemp, ib, my_ik, spin) = self%gw_vals(:, itemp, ib, my_ik, spin) &
1034 [ + + ]: 2720 : - 0.5* gt(1)*time_step - 0.5 * gt(nwr_ce)*exp(j_dpc*wrmesh_shifted_ce(:)*time_mesh(nwr_ce)) * time_step
1035 : :
1036 : :
1037 : : else
1038 : : msg = sjoin("FFT not available, using DFT instead.", &
1039 : : " Slower but reliable.", &
1040 : : " Parallelism over time for integration of G(t)")
1041 : : ABI_COMMENT(msg)
1042 : :
1043 : : do iw=1, nwr_ce
1044 : : ! if (mod(iw, self%wt_comm%nproc) /= self%wt_comm%me) cycle ! MPI parallelism over freqs
1045 : : omega = wrmesh_shifted_ce(iw)
1046 : :
1047 : : ! Discrete Fourier Transform of the Green's function
1048 : : g1(:) = exp( j_dpc * omega * time_mesh(:) ) * gt(:)
1049 : :
1050 : : ! Retarded Green's function in frequency domain
1051 : : self%gw_vals(iw, itemp, ib, my_ik, spin) = simpson_cplx(nwr_ce, time_step, g1)
1052 : :
1053 : : !if (my_rank == 0) write(ab_out, *)"gw_vals", self%gw_vals(iw, itemp, ib, my_ik, spin)
1054 : : end do ! iw
1055 : : end if
1056 : :
1057 : : ! Collect data if wt parallelism.
1058 : : ! call xmpi_sum(self%gw_vals(:, itemp, ib, my_ik, spin) , self%wt_comm%value, ierr)
1059 : :
1060 : : end do ! itemp
1061 : : end do ! ib
1062 : :
1063 : 12 : write(msg,'(4(a,i0),a,f8.2)') " k-point [", ikcalc, "/", self%nkcalc, "]"
1064 : 26 : call cwtime_report(msg, cpu, wall, gflops)
1065 : : end do ! my_ik
1066 : : end do ! my_spin
1067 : :
1068 [ + - ]: 2 : ABI_SFREE(c1)
1069 [ - + ]: 2 : ABI_FREE(ct_temp)
1070 [ - + ]: 2 : ABI_FREE(c_temp)
1071 [ + - ]: 2 : ABI_SFREE(c2)
1072 [ + - ]: 2 : ABI_SFREE(c3)
1073 [ + - ]: 2 : ABI_SFREE(ct)
1074 [ + - ]: 2 : ABI_SFREE(gt)
1075 [ + - ]: 2 : ABI_SFREE(g1)
1076 [ + - ]: 2 : ABI_SFREE(gw)
1077 [ + - ]: 2 : ABI_SFREE(beta)
1078 [ + - ]: 2 : ABI_SFREE(temp_g)
1079 [ + - ]: 2 : ABI_SFREE(temp_g_ce)
1080 [ + - ]: 2 : ABI_SFREE(temp_r)
1081 [ + - ]: 2 : ABI_SFREE(temp_r_cplx)
1082 [ + - ]: 2 : ABI_SFREE(time_mesh)
1083 [ + - ]: 2 : ABI_SFREE(wrmesh_shifted)
1084 [ + - ]: 2 : ABI_SFREE(wrmesh_shifted_ce)
1085 : : !ABI_SFREE(inv_wrmesh_shifted_sq)
1086 [ + - ]: 2 : ABI_SFREE(betaoverw2)
1087 [ - + ]: 2 : ABI_FREE(time_mesh_temp)
1088 : :
1089 : 2 : call cwtime_report(" cumulant_compute", cpu_kloop, wall_kloop, gflops_kloop)
1090 : :
1091 : : ! contains
1092 : :
1093 : : ! real function lreg(x,y, n, nsteps) result(m,b)
1094 : : ! !
1095 : : ! ! Determines a linear regression from x and y values
1096 : : ! ! What are the best m and b values to produce y= m*x + b ?
1097 : : ! !
1098 : : ! !
1099 : : ! ! Cost function ( Root Mean Squared Error ): J = 1/n sum_i^n (pred_i - y_i)^2
1100 : : ! ! where pred is the predicted value and y the true value
1101 : : ! !
1102 : : ! ! Goal: minimize J
1103 : : ! ! How? Using Gradient Descent
1104 : : ! ! - Learning rate is the step that the new value will be ( too small, takes longer; too large, it can be instable )
1105 : : ! ! - Initial m or b are chosen randomly
1106 : : ! ! - n is the number of points
1107 : : ! !
1108 : : ! ! new b = old b - 2*(learning rate)/n sum_i^n (pred(x_i) - y) * x_i
1109 : : ! ! new m = old m - 2*(learning rate)/n sum_i^n (pred(x_i) - y)
1110 : : ! !
1111 : : !
1112 : : ! integer, intent(in) :: n, nsteps
1113 : : ! real(dp), intent(in) :: x(n), y(n)
1114 : : ! !real(dp) :: m,b
1115 : : ! integer :: istep, i
1116 : : ! real(dp) :: lrate, cost, acc
1117 : : ! real(dp) :: pred(n)
1118 : : !
1119 : : ! m = 0 ! Initial guesses
1120 : : ! b = 0
1121 : : ! do istep=1,nsteps
1122 : : !
1123 : : ! pred(:) = m * x(:) + b ! Prediction with the new coefficients m, b
1124 : : !
1125 : : ! ! Check accuracy comparing the linear regression and the data
1126 : : ! do i=1, n
1127 : : ! if (abs(y(i) ) < 1e-4 ) cycle
1128 : : ! acc = acc + sum( abs(pred(i) - y(i))/y(i) )
1129 : : ! end do i
1130 : : ! acc = 1 - acc
1131 : : ! print *, "Accuracy: ",istep, acc
1132 : : !
1133 : : ! cost = 1.0/n * sum(pred(:) - y(:))**2
1134 : : !
1135 : : ! print *, "Cost: ", istep, cost
1136 : : !
1137 : : ! ! Update coefficients
1138 : : !
1139 : : ! m = m - 2.0*lrate/n * sum(pred(:) - y(:))
1140 : : ! b = b - 2.0*lrate/n * sum((pred(:) - y(:))*x(:))
1141 : : ! enddo ! istep
1142 : : !
1143 : : ! end function lreg
1144 : :
1145 : : ! complex function trapz(f_size, f_step, f)
1146 : : !
1147 : : ! integer,intent(in) :: f_size
1148 : : ! real(dp),intent(in) :: f_step
1149 : : ! complex(dp),intent(in) :: f(f_size)
1150 : : !
1151 : : ! trapz = ( sum(f) - 0.5* f(1) - 0.5* f(f_size) )* f_step
1152 : : !
1153 : : ! end function trapz
1154 : :
1155 [ - + ]: 4 : end subroutine cumulant_compute
1156 : : !!***
1157 : :
1158 : : !----------------------------------------------------------------------
1159 : :
1160 : : !!****f* m_cumulant/cumulant_kubo_transport
1161 : : !! NAME
1162 : : !! cumulant_kubo_transport
1163 : : !!
1164 : : !! FUNCTION
1165 : : !! Compute conductivity/mobility within Kubo Greenwood formula, using Cumulant self energy spectral function
1166 : : !! NB: also modified to use Dyson Migdal for comparison
1167 : : !!
1168 : : !! INPUTS
1169 : : !!
1170 : : !! SOURCE
1171 : :
1172 : 2 : subroutine cumulant_kubo_transport(self, dtset, cryst)
1173 : :
1174 : : !Arguments ------------------------------------
1175 : : class(cumulant_t),intent(inout) :: self
1176 : : type(dataset_type),intent(in) :: dtset
1177 : : type(crystal_t),intent(in) :: cryst
1178 : :
1179 : : !Local variables ------------------------------
1180 : : integer,parameter :: master = 0
1181 : : integer :: cnt
1182 : : integer :: nbands, ib, ikcalc, iw, spin, itemp, comm, ik_ibz
1183 : : integer :: ieh, ib_eph
1184 : : integer :: ii, jj, time_opt, isym_k, trev_k
1185 : : integer :: my_rank, nprocs, my_spin, my_ik
1186 : : real(dp) :: omega!, time_step!, time_max
1187 : : real(dp) :: cpu, wall, gflops, cpu_kloop, wall_kloop, gflops_kloop
1188 : : real(dp) :: eig_nk, sp_func, wr_step
1189 : : real(dp) :: sp_func_dm
1190 : : real(dp) :: integration
1191 : : real(dp) :: integration_dm
1192 : : real(dp) :: mu_e, max_occ, fact0, fact
1193 : : character(len=500) :: msg
1194 : : !arrays
1195 : 2 : real(dp),allocatable :: kernel(:), dfdw_acc(:),Aw(:),Aw_l0(:),dfdw_l0(:)!test_Aw(:), test_dfdw(:)
1196 : 2 : real(dp),allocatable :: kernel_dm(:), Aw_dm(:),Aw_l0_dm(:)
1197 : 2 : complex(dp),allocatable :: gdm_vals(:)
1198 : : real(dp) :: int_Aw
1199 : : real(dp) :: int_Aw_dm
1200 : : real(dp) :: int_dfdw, dfdw
1201 : : real(dp) :: vr(3), vv_tens(3,3), S(3,3), wtk!, spfunc
1202 : : ! real(dp), allocatable :: onsager_coeff(:,:)
1203 : : real(dp) :: work_33(3,3)
1204 : : real(dp) :: l0inv_33nw(3,3,2)
1205 : : real(dp) :: l0inv_33nw_dm(3,3,2)
1206 : : real(dp) :: Tkelv
1207 : :
1208 : : !************************************************************************
1209 : :
1210 : 2 : comm = self%comm
1211 : 2 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
1212 : :
1213 : 2 : call wrtout(std_out, " Computing conductivity using Kubo-Greenwood method.")
1214 : 2 : call cwtime(cpu_kloop, wall_kloop, gflops_kloop, "start")
1215 : :
1216 : :
1217 [ - + + - : 6 : ABI_MALLOC(kernel, (self%nwr))
+ - + - -
+ - + -
+ ]
1218 [ - + + - : 4 : ABI_MALLOC(kernel_dm, (self%nwr))
- + - + -
+ ]
1219 : : ! ABI_MALLOC(test_Aw, (self%ntemp))
1220 : : ! ABI_MALLOC(test_dfdw, (self%ntemp))
1221 [ - + + - : 4 : ABI_MALLOC(Aw, (self%nwr))
- + - + ]
1222 [ - + + - : 4 : ABI_MALLOC(Aw_dm, (self%nwr))
- + - + ]
1223 [ - + + - : 4 : ABI_MALLOC(dfdw_acc, (self%nwr))
- + - + ]
1224 [ - + + - : 6 : ABI_MALLOC(Aw_l0, (self%ntemp))
+ - + - -
+ - + ]
1225 [ - + + - : 4 : ABI_MALLOC(Aw_l0_dm, (self%ntemp))
- + - + ]
1226 [ - + + - : 4 : ABI_MALLOC(dfdw_l0, (self%ntemp))
- + - + ]
1227 : :
1228 [ - + + - : 6 : ABI_MALLOC(gdm_vals, (self%nwr))
+ - + - -
+ - + ]
1229 : :
1230 [ - + + - : 64 : ABI_CALLOC(self%l0, (3, 3, 2, self%nsppol, self%ntemp))
- + + - +
- + - - +
- + - + +
+ + + + +
+ + + + ]
1231 [ - + + - : 62 : ABI_CALLOC(self%l1, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1232 [ - + + - : 62 : ABI_CALLOC(self%l2, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1233 [ - + + - : 62 : ABI_CALLOC(self%seebeck, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1234 [ - + + - : 62 : ABI_CALLOC(self%kappa, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1235 : :
1236 [ - + + - : 62 : ABI_CALLOC(self%l0_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1237 [ - + + - : 62 : ABI_CALLOC(self%l1_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1238 [ - + + - : 62 : ABI_CALLOC(self%l2_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1239 [ - + + - : 62 : ABI_CALLOC(self%seebeck_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1240 [ - + + - : 62 : ABI_CALLOC(self%kappa_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1241 : :
1242 : :
1243 [ - + + - : 62 : ABI_CALLOC(self%mobility_mu, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1244 [ - + + - : 62 : ABI_CALLOC(self%conductivity_mu, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1245 : :
1246 [ - + + - : 62 : ABI_CALLOC(self%mobility_mu_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1247 [ - + + - : 62 : ABI_CALLOC(self%conductivity_mu_dm, (3, 3, 2, self%nsppol, self%ntemp))
- + + - -
+ - + - +
+ + + + +
+ + + +
+ ]
1248 : :
1249 [ - + + - : 280 : ABI_CALLOC(self%print_dfdw, (self%nwr, self%ntemp))
- + + - +
- + - - +
- + - + +
+ + + ]
1250 : :
1251 [ - + + - : 4 : ABI_MALLOC(self%transport_mu_e, (self%ntemp))
- + - + -
+ ]
1252 : :
1253 [ - + + - : 16 : ABI_CALLOC(self%n_ehst, (2, self%nsppol, self%ntemp))
- + + - +
- + - - +
- + - + +
+ + + +
+ ]
1254 : :
1255 : 2 : self%transport_fermie = dtset%eph_fermie
1256 : 2 : self%transport_extrael = dtset%eph_extrael
1257 [ + - - + : 6 : self%transport_mu_e = self%mu_e
- - - - +
+ ]
1258 [ - + - - ]: 2 : if (self%transport_fermie /= zero) self%transport_mu_e = self%transport_fermie
1259 : :
1260 [ + - - + ]: 2 : if (self%transport_fermie == zero .and. self%transport_extrael /= self%eph_extrael) then
1261 : :
1262 : : if (self%transport_extrael /= self%eph_extrael) then
1263 : : write(msg,'(2(a,e18.8),3a)') &
1264 : 0 : ' extrael from SIGEPH: ',self%transport_extrael, ' and input file: ',self%eph_extrael, "differ", ch10, &
1265 : 0 : ' Will recompute the chemical potential'
1266 : 0 : call wrtout(std_out, msg)
1267 : : end if
1268 : :
1269 : : ! Compute Fermi level for different T values.
1270 : 0 : call self%ebands%get_muT_with_fd(self%ntemp, self%kTmesh, dtset%spinmagntarget, dtset%prtvol, self%transport_mu_e, comm)
1271 : : end if
1272 : :
1273 : 2 : call self%ebands%get_carriers(self%ntemp, self%kTmesh, self%transport_mu_e, self%n_ehst)
1274 : :
1275 : :
1276 : 2 : time_opt =0
1277 : 2 : cnt = 0
1278 [ + + ]: 4 : do my_spin=1,self%my_nspins
1279 : 2 : spin = self%my_spins(my_spin)
1280 [ + + ]: 16 : do my_ik=1,self%my_nkcalc
1281 : 12 : ikcalc= self%my_ikcalc(my_ik)
1282 : : !if (ikcalc > 1) cycle
1283 : 12 : ik_ibz = self%kcalc2ibz(ikcalc, 1)
1284 : 12 : isym_k = self%kcalc2ibz(ikcalc, 2)
1285 : 12 : trev_k = self%kcalc2ibz(ikcalc, 6)
1286 : :
1287 : 12 : wtk = self%ebands%wtk(ik_ibz)
1288 : : ! TODO: this S is not used below - remove
1289 : : S = transpose(cryst%symrel_cart(:,:,isym_k))
1290 : :
1291 : 12 : nbands = self%nbcalc_ks(ikcalc, spin)
1292 : 12 : call cwtime(cpu, wall, gflops, "start")
1293 : : !do ib=self%bmin,self%bmax
1294 [ + + ]: 34 : do ib=self%bstart_ks(ikcalc, spin),self%bstart_ks(ikcalc, spin)+self%nbcalc_ks(ikcalc, spin)-1
1295 : : !if (ib > self%bmin) cycle ! MG HACK To be able to run tests quickly.
1296 : : !ib_eph = ib - self%bmin + 1
1297 : 20 : ib_eph = ib - self%bstart_ks(ikcalc, spin) + 1
1298 : 20 : eig_nk = self%ebands%eig(ib, ik_ibz, spin)
1299 : 20 : wr_step = self%wrmesh_b(2,ib_eph,my_ik,spin) - self%wrmesh_b(1,ib_eph,my_ik,spin)
1300 [ + + ]: 80 : vr(:) = self%vbks(:, ib, ik_ibz, spin)
1301 : : ! Store outer product (v_bks x v_bks) in vv_tens. This part does not depend on T and irta.
1302 [ + + ]: 80 : do ii=1,3
1303 [ + + ]: 260 : do jj=1,3
1304 : 240 : vv_tens(ii, jj) = vr(ii) * vr(jj)
1305 : : end do
1306 : : end do
1307 : : ! Calculation of the velocity tensor
1308 [ + + + + ]: 260 : vv_tens = cryst%symmetrize_cart_tens33(vv_tens, time_opt)
1309 : :
1310 : : ! temperature loop in Self Energy dependency
1311 [ + + ]: 52 : do itemp=1,self%ntemp
1312 : : !if (itemp > 1) cycle
1313 : :
1314 : : ! calculate Dyson Migdal Green's function as well, for comparison
1315 : : gdm_vals(:) = one / (self%wrmesh_b(:,ib_eph,my_ik,spin) &
1316 : : & - self%e0vals(ib_eph,my_ik,spin) - 0.5 * wr_step &
1317 [ + + ]: 2720 : & - self%vals_wr(:, itemp, ib_eph, my_ik, my_spin))
1318 : :
1319 : 20 : Tkelv = self%kTmesh(itemp) / kb_HaK; if (Tkelv < one) Tkelv = one
1320 [ + + ]: 2720 : do iw=1, self%nwr
1321 : : ! if (mod(iw, self%wt_comm%nproc) /= self%wt_comm%me) cycle ! MPI parallelism over freqs
1322 : :
1323 : : ! Preparing all elements needed for conductivity
1324 : 2700 : omega = self%wrmesh_b(iw,ib_eph,my_ik,my_spin)
1325 : :
1326 : : ! this retrieves the Cumulant spectral function
1327 : : ! TODO: add the Dyson Migdal as well, to compare properly the transport with the same KG equation
1328 : 2700 : sp_func = -aimag (self%gw_vals(iw, itemp, ib_eph, my_ik, my_spin) ) / pi
1329 : 2700 : sp_func_dm = -aimag (gdm_vals(iw)) / pi
1330 : :
1331 : :
1332 : 2700 : self%spfunc_dm_wr(iw, itemp, ib_eph, my_ik, my_spin) = sp_func_dm
1333 : :
1334 : 2700 : self%spfunc_wr(iw, itemp, ib_eph, my_ik, my_spin) = sp_func
1335 : : ! test_Aw(itemp) = test_Aw(itemp) + sp_func
1336 : 2700 : dfdw = occ_dfde(omega, self%kTmesh(itemp), self%mu_e(itemp))
1337 : 2700 : self%print_dfdw(iw,itemp) = dfdw
1338 : : ! test_dfdw(itemp) = test_dfdw(itemp) + dfdw
1339 : 2700 : kernel(iw) = - dfdw * sp_func**2
1340 : 2700 : kernel_dm(iw) = - dfdw * sp_func_dm**2
1341 : 2700 : Aw(iw) = sp_func**2
1342 : 2700 : Aw_dm(iw) = sp_func_dm**2
1343 : 2720 : dfdw_acc(iw) = dfdw
1344 : :
1345 : : end do !iw
1346 : 20 : mu_e = self%transport_mu_e(itemp)
1347 [ + + ]: 20 : ieh = 2; if (eig_nk >= mu_e) ieh = 1
1348 : 20 : integration = simpson( wr_step, kernel)
1349 : 20 : int_Aw = simpson(wr_step,Aw)
1350 : 20 : integration_dm = simpson( wr_step, kernel_dm)
1351 : 20 : int_Aw_dm = simpson(wr_step,Aw_dm)
1352 : 20 : int_dfdw = simpson(wr_step,dfdw_acc)
1353 : : ! Calculation of the conductivity
1354 : : self%l0( :, :, ieh, spin, itemp ) = self%l0( :, :, ieh, spin, itemp ) &
1355 [ + + + + ]: 260 : & + integration*vv_tens(:,:)*wtk
1356 : : self%l0_dm( :, :, ieh, spin, itemp ) = self%l0_dm( :, :, ieh, spin, itemp ) &
1357 [ + + + + ]: 260 : & + integration_dm*vv_tens(:,:)*wtk
1358 : :
1359 : : ! TODO: document this bit, why only use the xx component of vv_tens?? Aw_l0 is never used or output
1360 : 20 : Aw_l0(itemp) = Aw_l0(itemp) + int_Aw*wtk*vv_tens(1,1)
1361 : 20 : Aw_l0_dm(itemp) = Aw_l0_dm(itemp) + int_Aw_dm*wtk*vv_tens(1,1)
1362 : 20 : dfdw_l0(itemp) = dfdw_l0(itemp) + int_dfdw*wtk*vv_tens(1,1)
1363 : : ! END TODO
1364 : :
1365 : : ! TODO : check these equations for l1 l2, they should not depend on the band indices, as we are inside the ib loop!!
1366 : : self%l1( :, :, ieh, spin, itemp ) = self%l1( :, :, ieh, spin, itemp ) &
1367 [ + + + + ]: 260 : & + integration*vv_tens(:,:)*wtk*(eig_nk - self%mu_e(itemp))
1368 : : self%l2( :, :, ieh, spin, itemp ) = self%l2( :, :, ieh, spin, itemp ) &
1369 [ + + + + ]: 260 : & + integration*vv_tens(:,:)*wtk*(eig_nk - self%mu_e(itemp))**2
1370 : : self%l1_dm( :, :, ieh, spin, itemp ) = self%l1_dm( :, :, ieh, spin, itemp ) &
1371 [ + + + + ]: 260 : & + integration_dm*vv_tens(:,:)*wtk*(eig_nk - self%mu_e(itemp))
1372 : : self%l2_dm( :, :, ieh, spin, itemp ) = self%l2_dm( :, :, ieh, spin, itemp ) &
1373 [ + + + + ]: 280 : & + integration_dm*vv_tens(:,:)*wtk*(eig_nk - self%mu_e(itemp))**2
1374 : :
1375 : : end do ! itemp
1376 : :
1377 : : end do !ib
1378 : :
1379 : : end do ! my_ik
1380 : :
1381 : : end do !my_spin
1382 : :
1383 : : ! Collect data if k-points parallelism. TODO: if this is reactivated, mpi_sum the l0 l1 l2 arrays instead
1384 : : !call xmpi_sum(self%l0, self%kcalc_comm%value, ierr)
1385 : : !call xmpi_sum(self%l1, self%kcalc_comm%value, ierr)
1386 : : !call xmpi_sum(self%l2, self%kcalc_comm%value, ierr)
1387 : :
1388 [ + + ]: 4 : do itemp = 1, self%ntemp
1389 : 2 : Tkelv = self%kTmesh(itemp) / kb_HaK; if (Tkelv < one) Tkelv = one
1390 [ + + ]: 6 : do my_spin=1,self%my_nspins
1391 : 2 : spin = self%my_spins(my_spin)
1392 [ + + ]: 8 : do ieh = 1, 2
1393 : : ! calculate the transport coefficients from the l0 l1 l2
1394 : 4 : call inv33(self%l0(:, :, ieh, spin, itemp), l0inv_33nw(:,:,ieh))
1395 [ + - + + : 580 : self%seebeck(:,:,ieh,spin,itemp) = matmul(l0inv_33nw(:,:,ieh), self%l1(:,:,ieh,spin,itemp)) / Tkelv
+ + + + +
+ + + + -
+ - + - +
- - - - -
+ + + + ]
1396 : :
1397 : 4 : call inv33(self%l0_dm(:, :, ieh, spin, itemp), l0inv_33nw_dm(:,:,ieh))
1398 [ + - + + : 580 : self%seebeck_dm(:,:,ieh,spin,itemp) = matmul(l0inv_33nw_dm(:,:,ieh), self%l1_dm(:,:,ieh,spin,itemp)) / Tkelv
+ + + + +
+ + + + -
+ - + - +
- - - - -
+ + + + ]
1399 : :
1400 [ + + + + ]: 52 : work_33 = self%l1(:, :, ieh, spin, itemp)
1401 : : ! TODO: check if one of the work_33 below is incorrect or needs a transpose
1402 [ + + + + : 364 : work_33 = self%l2(:, :, ieh, spin, itemp) - matmul(work_33, matmul(l0inv_33nw(:, :, ieh), work_33))
+ + + + +
+ + + + +
+ + ]
1403 [ + + + + ]: 52 : self%kappa(:,:,ieh,spin,itemp) = work_33 / Tkelv
1404 : :
1405 [ + + + + ]: 52 : work_33 = self%l1_dm(:, :, ieh, spin, itemp)
1406 [ + + + + : 364 : work_33 = self%l2_dm(:, :, ieh, spin, itemp) - matmul(work_33, matmul(l0inv_33nw_dm(:, :, ieh), work_33))
+ + + + +
+ + + + +
+ + ]
1407 [ + + + + ]: 54 : self%kappa_dm(:,:,ieh,spin,itemp) = work_33 / Tkelv
1408 : :
1409 : : !TODO: leftover below - can this be removed?
1410 : : !self%conductivity_mu( :, :, ieh, spin, itemp ) = self%conductivity_mu( :, :, ieh, spin, itemp ) + integration*vv_tens(:,:)*wtk
1411 : : end do ! ieh
1412 : : end do !my_spin
1413 : : end do ! itemp
1414 : :
1415 : 2 : max_occ = two / (self%nspinor * self%nsppol)
1416 : 2 : fact0 = max_occ * (siemens_SI / Bohr_meter / cryst%ucvol) / 100
1417 [ + - + - : 60 : self%conductivity_mu = fact0 * self%l0 ! siemens cm^-1
+ - + - +
- - + - -
+ + + + +
+ + + +
+ ]
1418 [ + + + + : 58 : self%seebeck = - volt_SI * max_occ * self%seebeck
+ + + + +
+ ]
1419 [ + + + + : 58 : self%kappa = + volt_SI**2 * fact0 * self%kappa
+ + + + +
+ ]
1420 : :
1421 [ + - + - : 60 : self%conductivity_mu_dm = fact0 * self%l0_dm ! siemens cm^-1
+ - + - +
- - + - -
+ + + + +
+ + + +
+ ]
1422 [ + + + + : 58 : self%seebeck_dm = - volt_SI * max_occ * self%seebeck_dm
+ + + + +
+ ]
1423 [ + + + + : 58 : self%kappa_dm = + volt_SI**2 * fact0 * self%kappa_dm
+ + + + +
+ ]
1424 : :
1425 : : ! Scale by the carrier concentration
1426 : 2 : fact = 100**3 / e_Cb
1427 [ + + ]: 4 : do my_spin=1,self%my_nspins
1428 : 2 : spin = self%my_spins(my_spin)
1429 [ + + ]: 6 : do itemp=1,self%ntemp
1430 [ + + ]: 8 : do ieh=1,2 ! e/h
1431 : : call safe_div(fact * self%conductivity_mu(:,:,ieh,spin,itemp), &
1432 : : self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, &
1433 [ + + + + ]: 52 : self%mobility_mu(:,:,ieh,spin,itemp))
1434 : : call safe_div(fact * self%conductivity_mu_dm(:,:,ieh,spin,itemp), &
1435 : : self%n_ehst(ieh, spin, itemp) / cryst%ucvol / Bohr_meter**3, zero, &
1436 [ + + + + ]: 54 : self%mobility_mu_dm(:,:,ieh,spin,itemp))
1437 : : end do
1438 : : end do
1439 : : end do
1440 : :
1441 : :
1442 : 2 : call cwtime_report(" cumulant_kubo_transport", cpu_kloop, wall_kloop, gflops_kloop)
1443 : :
1444 [ + - ]: 2 : ABI_SFREE(kernel)
1445 [ + - ]: 2 : ABI_SFREE(kernel_dm)
1446 : : ! ABI_SFREE(test_Aw)
1447 : : ! ABI_SFREE(test_dfdw)
1448 [ + - ]: 2 : ABI_SFREE(Aw)
1449 [ + - ]: 2 : ABI_SFREE(Aw_dm)
1450 [ + - ]: 2 : ABI_SFREE(dfdw_acc)
1451 [ + - ]: 2 : ABI_SFREE(Aw_l0)
1452 [ + - ]: 2 : ABI_SFREE(Aw_l0_dm)
1453 [ + - ]: 2 : ABI_SFREE(dfdw_l0)
1454 [ + - ]: 2 : ABI_SFREE(gdm_vals)
1455 : :
1456 [ - + - + : 6 : end subroutine cumulant_kubo_transport
- + - + -
+ ]
1457 : :
1458 : : !!***
1459 : :
1460 : : !----------------------------------------------------------------------
1461 : :
1462 : : !!****f* m_cumulant/cumulant_sigmaph_ncread
1463 : : !! NAME
1464 : : !! cumulant_sigmaph_ncread
1465 : : !!
1466 : : !! FUNCTION
1467 : : !! read out_SIGPH.nc file
1468 : : !!
1469 : : !! INPUTS
1470 : : !! cryst<crystal_t>=Crystalline structure
1471 : : !! dtset<dataset_type>=All input variables for this dataset.
1472 : : !! ncid=Netcdf file handle.
1473 : : !!
1474 : : !! SOURCE
1475 : :
1476 : 2 : subroutine cumulant_sigmaph_ncread(self, path, ncid, comm)
1477 : :
1478 : : !Arguments --------------------------------------
1479 : : class(cumulant_t),intent(inout) :: self
1480 : : character(len=fnlen),intent(in) :: path
1481 : : integer,intent(in) :: comm
1482 : : integer,intent(out) :: ncid
1483 : :
1484 : : !Local variables --------------------------------
1485 : : integer :: ierr
1486 : : real(dp) :: cpu, wall, gflops
1487 : : character(len=1000) :: msg
1488 : : !real(dp), allocatable :: vals_wr(:,:,:,:,:,:),vals_e0ks(:,:,:,:,:)
1489 : :
1490 : : !************************************************************************
1491 : :
1492 : 2 : call cwtime(cpu, wall, gflops, "start")
1493 : :
1494 : : ! Open netcdf file
1495 : 2 : ierr = 0
1496 : :
1497 [ - + ]: 2 : if (.not. file_exists(path)) then
1498 : 0 : msg = sjoin("Cannot find file", path)
1499 : 0 : ierr = 1; return
1500 : : end if
1501 : :
1502 : 2 : call cwtime(cpu, wall, gflops, "start")
1503 [ - + ]: 2 : NCF_CHECK(nctk_open_read(ncid, path, comm))
1504 : :
1505 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "nkcalc", self%nkcalc))
1506 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "max_nbcalc", self%max_nbcalc))
1507 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "nsppol", self%nsppol))
1508 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "number_of_spinor_components", self%nspinor))
1509 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "ntemp", self%ntemp))
1510 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "nwr", self%nwr))
1511 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "nqbz", self%nqbz))
1512 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "nqibz", self%nqibz))
1513 [ - + ]: 2 : NCF_CHECK(nctk_get_dim(ncid, "natom3", self%natom3))
1514 : :
1515 [ - + + - : 6 : ABI_MALLOC(self%kcalc, (3, self%nkcalc))
+ - + - -
+ - + -
+ ]
1516 [ - + + - : 8 : ABI_MALLOC(self%nbcalc_ks, (self%nkcalc, self%nsppol))
- + + - +
- + - - +
- + - + ]
1517 [ - + + - : 6 : ABI_MALLOC(self%bstart_ks, (self%nkcalc, self%nsppol))
- + + - -
+ - + -
+ ]
1518 [ - + + - : 6 : ABI_MALLOC(self%kcalc2ibz, (self%nkcalc, 6))
+ - + - -
+ - + -
+ ]
1519 [ - + + - : 6 : ABI_MALLOC(self%kTmesh, (self%ntemp))
+ - + - -
+ - + -
+ ]
1520 : : !ABI_MALLOC(self%wrmesh_b, (self%nwr, self%max_nbcalc, self%nkcalc, self%nsppol))
1521 : : !ABI_MALLOC(self%vals_wr, ( self%nwr, self%ntemp, self%max_nbcalc, self%nkcalc, self%nsppol))
1522 : : !ABI_MALLOC(self%vals_e0ks, ( self%ntemp, self%max_nbcalc, self%nkcalc, self%nsppol))
1523 : : !ABI_MALLOC(self%e0vals, (self%max_nbcalc, self%nkcalc, self%nsppol))
1524 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("ngqpt"), self%ngqpt))
1525 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("kcalc"), self%kcalc))
1526 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("kTmesh"), self%kTmesh))
1527 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("nbcalc_ks"), self%nbcalc_ks))
1528 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("bstart_ks"), self%bstart_ks))
1529 [ - + ]: 2 : NCF_CHECK(nf90_get_var(ncid, vid("kcalc2ibz"), self%kcalc2ibz))
1530 : :
1531 : 6 : call cwtime_report(" sigmaph_ncread", cpu, wall, gflops)
1532 : :
1533 : : contains
1534 : 12 : integer function vid(vname)
1535 : : character(len=*),intent(in) :: vname
1536 : 12 : vid = nctk_idname(ncid, vname)
1537 : : end function vid
1538 : :
1539 : : end subroutine cumulant_sigmaph_ncread
1540 : : !!***
1541 : :
1542 : : !----------------------------------------------------------------------
1543 : :
1544 : : !!****f* m_cumulant/cumulant_ncwrite
1545 : : !! NAME
1546 : : !! cumulant_ncwrite
1547 : : !!
1548 : : !! FUNCTION
1549 : : !!
1550 : : !! INPUTS
1551 : : !! path=Filenae of output netcdf file.
1552 : : !! cryst<crystal_t>=Crystalline structure
1553 : : !! dtset<dataset_type>=All input variables for this dataset.
1554 : : !! ncid=Netcdf file handle.
1555 : : !!
1556 : : !! SOURCE
1557 : :
1558 : : !subroutine cumulant_ncwrite(self, path, cryst, ebands, dtset)
1559 : 6 : subroutine cumulant_ncwrite(self, path, cryst, dtset)
1560 : :
1561 : : !Arguments --------------------------------------
1562 : : class(cumulant_t),intent(in) :: self
1563 : : type(crystal_t),intent(in) :: cryst
1564 : : ! type(ebands_t),intent(in) :: ebands
1565 : : type(dataset_type),intent(in) :: dtset
1566 : : character(len=*),intent(in) :: path
1567 : :
1568 : : !Local variables --------------------------------
1569 : : integer,parameter :: master = 0
1570 : : integer :: ncerr, ncid, ikcalc, spin, my_ik, ib, itemp, ntemp
1571 : : integer :: ii, ieh
1572 : : real(dp) :: cpu, wall, gflops
1573 : : character(len=1000) :: msg
1574 : :
1575 : : !************************************************************************
1576 : :
1577 : : !comm = self%comm my_rank =
1578 [ + + ]: 6 : call wrtout([std_out, ab_out], ch10//sjoin("- Writing cumulant results to:", path))
1579 : 2 : call cwtime(cpu, wall, gflops, "start")
1580 : :
1581 : : ! Only one proc create the file, write structure and define basic dimensions.
1582 : : ! Then we reopen the file in MPI-IO mode.
1583 : :
1584 [ + - ]: 2 : if (xmpi_comm_rank(self%comm) == master) then
1585 : :
1586 [ - + ]: 2 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
1587 : :
1588 : : ! Write to netcdf file
1589 [ - + ]: 2 : NCF_CHECK(cryst%ncwrite(ncid))
1590 : : ! FIXME: Cannot write ebands because it crashes in
1591 : : ! k_dependent = "no"; if (any(ebands%nband(1) /= ebands%nband)) k_dependent = "yes"
1592 : : ! Should understand why!
1593 : : !NCF_CHECK(ebands_ncwrite(ebands, ncid))
1594 : :
1595 : : ! Add cumulant dimensions.
1596 : : ncerr = nctk_def_dims(ncid, [ &
1597 : : nctkdim_t("nkcalc", self%nkcalc), nctkdim_t("max_nbcalc", self%max_nbcalc), &
1598 : : nctkdim_t("nsppol", self%nsppol), nctkdim_t("ntemp", self%ntemp), &
1599 : : nctkdim_t("nqbz", self%nsppol), nctkdim_t("nqibz", self%ntemp), &
1600 : : nctkdim_t("nwr", self%nwr)], &
1601 [ + + ]: 16 : defmode=.True.)
1602 [ - + ]: 2 : NCF_CHECK(ncerr)
1603 : : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
1604 : : "eph_task", "nbsum", "bsum_start", "bsum_stop", "symdynmat", &
1605 : : "ph_intmeth", "eph_intmeth", "qint_method", "eph_transport", &
1606 [ + + ]: 30 : "imag_only", "symv1scf", "dvdb_add_lr", "mrta", "ibte_prep"])
1607 [ - + ]: 2 : NCF_CHECK(ncerr)
1608 : : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
1609 [ + + ]: 20 : "eta", "wr_step", "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear", "eph_phwinfact"])
1610 [ - + ]: 2 : NCF_CHECK(ncerr)
1611 : :
1612 : :
1613 : : ! Define arrays. Note nkcalc instead of my_nkcalc
1614 : : ncerr = nctk_def_arrays(ncid, [ &
1615 : : nctkarr_t("bstart_ks", "int", "nkcalc, nsppol"), &
1616 : : nctkarr_t("nbcalc_ks", "int", "nkcalc, nsppol"), &
1617 : : nctkarr_t("ngqpt", "dp", "three"), &
1618 : : nctkarr_t("kcalc", "dp", "three, nkcalc"), &
1619 : : nctkarr_t("kcalc2ibz", "dp", " nkcalc, six"), &
1620 : : nctkarr_t("kTmesh", "dp", "ntemp"), &
1621 : : !nctkarr_t("mu_e", "dp", "ntemp"), &
1622 : : nctkarr_t("wrmesh_b", "dp", "nwr, max_nbcalc, nkcalc, nsppol"), &
1623 : : !nctkarr_t("vals_wr", "dp", "two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1624 : : nctkarr_t("gw_vals", "dp", "two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1625 : : nctkarr_t("ks_enes", "dp", "max_nbcalc, nkcalc, nsppol"), &
1626 : : nctkarr_t("dw_vals", "dp", "ntemp, max_nbcalc, nkcalc, nsppol"), &
1627 : : nctkarr_t("spfunc_dm_wr", "dp", "nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1628 : : nctkarr_t("spfunc_wr", "dp", "nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1629 : : nctkarr_t('dfdw',"dp", "nwr, ntemp"), &
1630 : : nctkarr_t('conductivity_mu',"dp", "three, three, two, nsppol, ntemp"), &
1631 : : nctkarr_t('conductivity_mu_dm',"dp", "three, three, two, nsppol, ntemp"), &
1632 : : nctkarr_t('mobility_mu', "dp", "three, three, two, nsppol, ntemp"), &
1633 : : nctkarr_t('mobility_mu_dm', "dp", "three, three, two, nsppol, ntemp"), &
1634 : : nctkarr_t('seebeck',"dp", "three, three, two, nsppol, ntemp"), &
1635 : : nctkarr_t('seebeck_dm',"dp", "three, three, two, nsppol, ntemp"), &
1636 : : nctkarr_t('kappa',"dp", "three, three, two, nsppol, ntemp"), &
1637 : : nctkarr_t('kappa_dm',"dp", "three, three, two, nsppol, ntemp") &
1638 [ + + ]: 44 : ])
1639 [ - + ]: 2 : NCF_CHECK(ncerr)
1640 : :
1641 [ - + ]: 2 : if (self%debug == 1) then
1642 : : ncerr = nctk_def_arrays(ncid, [ &
1643 : : nctkarr_t("time_mesh", "dp", "nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1644 : : nctkarr_t("ct_vals", "dp","two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1645 : : nctkarr_t("c1", "dp","two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1646 : : nctkarr_t("c2", "dp","two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1647 : : nctkarr_t("c3", "dp","two, nwr, ntemp, max_nbcalc, nkcalc, nsppol"), &
1648 [ # # ]: 0 : nctkarr_t("gt_vals", "dp","two, nwr, ntemp, max_nbcalc, nkcalc, nsppol") ] )
1649 [ # # ]: 0 : NCF_CHECK(ncerr)
1650 : : endif
1651 : :
1652 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eta"), self%ieta))
1653 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nbsum"), self%nbsum))
1654 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), self%ngqpt))
1655 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc"), self%kcalc))
1656 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kTmesh"), self%kTmesh))
1657 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nbcalc_ks"), self%nbcalc_ks))
1658 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "bstart_ks"), self%bstart_ks))
1659 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kcalc2ibz"), self%kcalc2ibz))
1660 [ - + ]: 2 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "dfdw"), self%print_dfdw))
1661 : : ! FIXME This part is wrong since these arrays are MPI distributed
1662 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "conductivity_mu"), self%conductivity_mu))
1663 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "mobility_mu"), self%mobility_mu))
1664 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "seebeck"), self%seebeck))
1665 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kappa"), self%kappa))
1666 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "wrmesh_b"), self%wrmesh_b))
1667 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vals_wr"), c2r(self%vals_wr)))
1668 : : !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ks_enes"), self%e0vals))
1669 : :
1670 [ - + ]: 2 : NCF_CHECK(nf90_close(ncid))
1671 : : end if ! master
1672 : :
1673 : : ! Barrier to avoid race conditions.
1674 : : !call wrtout(std_out, "before barrier")
1675 : 2 : call xmpi_barrier(self%comm)
1676 : :
1677 : : ! Only the procs in the ncwrite_comm communicator write to disk.
1678 [ + - ]: 2 : if (self%ncwrite_comm%value == xmpi_comm_null) goto 100
1679 : :
1680 : : ! open file for parallel-IO mode inside comm. All procs in ncwrite_comm enter this part.
1681 : 2 : call wrtout(std_out, sjoin(" Performing parallel IO with:", itoa(self%ncwrite_comm%nproc), "procs"))
1682 [ - + ]: 2 : NCF_CHECK(nctk_open_modify(ncid, path, self%ncwrite_comm%value))
1683 : :
1684 [ - + ]: 2 : NCF_CHECK(nctk_set_datamode(ncid))
1685 : :
1686 : : ! Activate collectve IO.
1687 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "gw_vals"), nf90_collective)
1688 [ - + ]: 2 : NCF_CHECK(ncerr)
1689 : :
1690 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "wrmesh_b"), nf90_collective)
1691 [ - + ]: 2 : NCF_CHECK(ncerr)
1692 : :
1693 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "ks_enes"), nf90_collective)
1694 [ - + ]: 2 : NCF_CHECK(ncerr)
1695 : :
1696 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "dw_vals"), nf90_collective)
1697 [ - + ]: 2 : NCF_CHECK(ncerr)
1698 : :
1699 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "spfunc_dm_wr"), nf90_collective)
1700 [ - + ]: 2 : NCF_CHECK(ncerr)
1701 : :
1702 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "spfunc_wr"), nf90_collective)
1703 [ - + ]: 2 : NCF_CHECK(ncerr)
1704 : :
1705 [ + - + + : 3 : if (any(abs(dtset%sigma_erange) > zero)) then
+ - ]
1706 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "dfdw"), nf90_collective)
1707 [ - + ]: 2 : NCF_CHECK(ncerr)
1708 : :
1709 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "conductivity_mu"), nf90_collective)
1710 [ - + ]: 2 : NCF_CHECK(ncerr)
1711 : :
1712 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "mobility_mu"), nf90_collective)
1713 [ - + ]: 2 : NCF_CHECK(ncerr)
1714 : :
1715 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "seebeck"), nf90_collective)
1716 [ - + ]: 2 : NCF_CHECK(ncerr)
1717 : :
1718 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "kappa"), nf90_collective)
1719 [ - + ]: 2 : NCF_CHECK(ncerr)
1720 : :
1721 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "conductivity_mu_dm"), nf90_collective)
1722 [ - + ]: 2 : NCF_CHECK(ncerr)
1723 : :
1724 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "mobility_mu_dm"), nf90_collective)
1725 [ - + ]: 2 : NCF_CHECK(ncerr)
1726 : :
1727 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "seebeck_dm"), nf90_collective)
1728 [ - + ]: 2 : NCF_CHECK(ncerr)
1729 : :
1730 : 2 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "kappa_dm"), nf90_collective)
1731 [ - + ]: 2 : NCF_CHECK(ncerr)
1732 : :
1733 : :
1734 : : end if
1735 : :
1736 [ - + ]: 2 : if (self%debug == 1) then
1737 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "time_mesh"), nf90_collective)
1738 [ # # ]: 0 : NCF_CHECK(ncerr)
1739 : :
1740 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "ct_vals"), nf90_collective)
1741 [ # # ]: 0 : NCF_CHECK(ncerr)
1742 : :
1743 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "c1"), nf90_collective)
1744 [ # # ]: 0 : NCF_CHECK(ncerr)
1745 : :
1746 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "c2"), nf90_collective)
1747 [ # # ]: 0 : NCF_CHECK(ncerr)
1748 : :
1749 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "c3"), nf90_collective)
1750 [ # # ]: 0 : NCF_CHECK(ncerr)
1751 : :
1752 : 0 : ncerr = nf90_var_par_access(ncid, nctk_idname(ncid, "gt_vals"), nf90_collective)
1753 [ # # ]: 0 : NCF_CHECK(ncerr)
1754 : : endif
1755 : 2 : spin = self%my_spins(1)
1756 : 2 : ikcalc = self%my_ikcalc(1) ! index of the first kcalc treated by this rank.
1757 : :
1758 : : ! Start to write my **contiguous block** of kpoints from this **global** location
1759 : : ! Each MPI proc writes my_nkcalc entries.
1760 : :
1761 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "gw_vals"), c2r(self%gw_vals), &
1762 : : start=[1,1,1,1,ikcalc,spin], &
1763 [ + - + + : 26 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
+ + ]
1764 [ - + ]: 2 : NCF_CHECK(ncerr)
1765 : :
1766 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "wrmesh_b"), self%wrmesh_b, &
1767 : : start=[1,1,ikcalc,spin], &
1768 [ + + + + ]: 18 : count=[self%nwr, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1769 [ - + ]: 2 : NCF_CHECK(ncerr)
1770 : :
1771 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "ks_enes"), self%e0vals, &
1772 : : start=[1,ikcalc,spin], &
1773 [ + + + + ]: 14 : count=[self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1774 [ - + ]: 2 : NCF_CHECK(ncerr)
1775 : :
1776 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "dw_vals"), self%dw_vals, &
1777 : : start=[1,1,ikcalc,spin], &
1778 [ + + + + ]: 18 : count=[self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1779 [ - + ]: 2 : NCF_CHECK(ncerr)
1780 : :
1781 : : !DEBUG STUFF MJV
1782 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "spfunc_dm_wr"), self%spfunc_dm_wr, &
1783 : : start=[1,1,1,ikcalc,spin], &
1784 [ + + + + ]: 22 : count=[self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1785 [ - + ]: 2 : NCF_CHECK(ncerr)
1786 : :
1787 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "spfunc_wr"), self%spfunc_wr, &
1788 : : start=[1,1,1,ikcalc,spin], &
1789 [ + + + + ]: 22 : count=[self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1790 [ - + ]: 2 : NCF_CHECK(ncerr)
1791 : :
1792 : :
1793 [ + - + + : 3 : if (any(abs(dtset%sigma_erange) > zero)) then
+ - ]
1794 : :
1795 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "seebeck"), self%seebeck, &
1796 : : start=[1,1,1,spin,1], &
1797 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1798 [ - + ]: 2 : NCF_CHECK(ncerr)
1799 : :
1800 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "kappa"), self%kappa, &
1801 : : start=[1,1,1,spin,1], &
1802 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1803 [ - + ]: 2 : NCF_CHECK(ncerr)
1804 : :
1805 : :
1806 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "mobility_mu"), self%mobility_mu, &
1807 : : start=[1,1,1,spin,1], &
1808 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1809 [ - + ]: 2 : NCF_CHECK(ncerr)
1810 : :
1811 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "conductivity_mu"), self%conductivity_mu, &
1812 : : start=[1,1,1,spin,1], &
1813 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1814 [ - + ]: 2 : NCF_CHECK(ncerr)
1815 : :
1816 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "seebeck_dm"), self%seebeck_dm, &
1817 : : start=[1,1,1,spin,1], &
1818 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1819 [ - + ]: 2 : NCF_CHECK(ncerr)
1820 : :
1821 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "kappa_dm"), self%kappa_dm, &
1822 : : start=[1,1,1,spin,1], &
1823 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1824 [ - + ]: 2 : NCF_CHECK(ncerr)
1825 : :
1826 : :
1827 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "mobility_mu_dm"), self%mobility_mu_dm, &
1828 : : start=[1,1,1,spin,1], &
1829 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1830 [ - + ]: 2 : NCF_CHECK(ncerr)
1831 : :
1832 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "conductivity_mu_dm"), self%conductivity_mu_dm, &
1833 : : start=[1,1,1,spin,1], &
1834 [ + + + + ]: 22 : count=[3, 3, 2, self%my_nspins , self%ntemp])
1835 [ - + ]: 2 : NCF_CHECK(ncerr)
1836 : :
1837 : : end if
1838 : :
1839 : :
1840 [ - + ]: 2 : if (self%debug == 1) then
1841 : :
1842 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "time_mesh"), self%time_mesh, &
1843 : : start=[1,1,1,ikcalc,spin], &
1844 [ # # # # ]: 0 : count=[self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
1845 [ # # ]: 0 : NCF_CHECK(ncerr)
1846 : :
1847 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "ct_vals"), c2r(self%ct_vals), &
1848 : : start=[1,1,1,1,ikcalc,spin], &
1849 [ # # # # : 0 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
# # ]
1850 [ # # ]: 0 : NCF_CHECK(ncerr)
1851 : :
1852 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "c1"), c2r(self%c1), &
1853 : : start=[1,1,1,1,ikcalc,spin], &
1854 [ # # # # : 0 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
# # ]
1855 [ # # ]: 0 : NCF_CHECK(ncerr)
1856 : :
1857 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "c2"), c2r(self%c2), &
1858 : : start=[1,1,1,1,ikcalc,spin], &
1859 [ # # # # : 0 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
# # ]
1860 [ # # ]: 0 : NCF_CHECK(ncerr)
1861 : :
1862 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "c3"), c2r(self%c3), &
1863 : : start=[1,1,1,1,ikcalc,spin], &
1864 [ # # # # : 0 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
# # ]
1865 [ # # ]: 0 : NCF_CHECK(ncerr)
1866 : :
1867 : : ncerr = nf90_put_var(ncid, nctk_idname(ncid, "gt_vals"), c2r(self%gt_vals), &
1868 : : start=[1,1,1,1,ikcalc,spin], &
1869 [ # # # # : 0 : count=[2, self%nwr, self%ntemp, self%max_nbcalc, self%my_nkcalc, self%my_nspins])
# # ]
1870 [ # # ]: 0 : NCF_CHECK(ncerr)
1871 : : endif
1872 : :
1873 [ - + ]: 2 : NCF_CHECK(nf90_close(ncid))
1874 : :
1875 : : ! Write to ab_out for automatic testing.
1876 [ + - + - ]: 2 : if (xmpi_comm_rank(self%comm) == master .and. is_open(ab_out)) then
1877 : 2 : write(ab_out, "(/,a)")" Print first 10 frequencies in gw_vals array (re-im) for testing purposes:"
1878 : 2 : write(ab_out, "(2(a, i0))")" spin: ", spin, ", ikcalc: ", ikcalc
1879 : 2 : my_ik = 1
1880 [ + + ]: 7 : do ib=1,self%nbcalc_ks(ikcalc, spin)
1881 [ + + ]: 12 : do itemp=1,self%ntemp
1882 : 5 : write(ab_out, "(2(a,i0))")" gw_vals for itemp:", itemp, "ib: ", ib
1883 [ + + ]: 55 : write(ab_out, "(*(es13.5))")dble(self%gw_vals(1:min(10, self%nwr), itemp, ib, my_ik, spin))
1884 [ + + ]: 60 : write(ab_out, "(*(es13.5))")aimag(self%gw_vals(1:min(10, self%nwr), itemp, ib, my_ik, spin))
1885 : : end do
1886 : : end do
1887 : : end if
1888 [ + - + + : 3 : if (xmpi_comm_rank(self%comm) == master .and. is_open(ab_out) .and. any(abs(dtset%sigma_erange) > zero)) then
+ - + - +
- ]
1889 : : msg = sjoin(" Print first 5 temperatures of diagonal mobility_mu", &
1890 [ + - ]: 2 : " > 1e-6 (with ieh as electrons or holes) for testing purposes:")
1891 : 2 : write(ab_out, "(/,a)") trim(msg)
1892 : 2 : write(ab_out, "(2(a, i0))")" spin: ", spin
1893 : 2 : if (self%ntemp > 5) then
1894 : : ntemp = 5
1895 : : else
1896 : : ntemp = self%ntemp
1897 : : end if
1898 [ + + ]: 6 : do ieh=1,2
1899 [ + + ]: 18 : do ii=1,3
1900 [ + + ]: 28 : do itemp=1,ntemp
1901 [ + + ]: 24 : if (dble(self%mobility_mu(ii, ii, ieh, spin, itemp)) > 1e-6) then
1902 : 6 : write(ab_out, "(3(a,i0))")" mobility_mu for itemp:", itemp, " ieh: ", ieh, " xyz: ", ii
1903 : 6 : write(ab_out, "(*(es13.5))")dble(self%mobility_mu(ii, ii, ieh, spin, itemp))
1904 : : end if
1905 : : end do
1906 : : end do
1907 : : end do
1908 : : end if
1909 : :
1910 : :
1911 : 2 : 100 call cwtime_report(" cumulant_ncwrite", cpu, wall, gflops)
1912 : :
1913 : 2 : end subroutine cumulant_ncwrite
1914 : : !!***
1915 : :
1916 : : !----------------------------------------------------------------------
1917 : :
1918 : : !!****f* m_cumulant/cumulant_free
1919 : : !! NAME
1920 : : !! cumulant_free
1921 : : !!
1922 : : !! FUNCTION
1923 : : !! Free dynamic memory.
1924 : : !!
1925 : : !! INPUTS
1926 : : !!
1927 : : !! SOURCE
1928 : :
1929 : 2 : subroutine cumulant_free(self)
1930 : :
1931 : : !Arguments --------------------------------------
1932 : : class(cumulant_t),intent(inout) :: self
1933 : :
1934 : : !************************************************************************
1935 : :
1936 [ + - ]: 2 : ABI_SFREE(self%kcalc2ebands)
1937 [ + - ]: 2 : ABI_SFREE(self%linewidths)
1938 [ + - ]: 2 : ABI_SFREE(self%vbks)
1939 [ + - ]: 2 : ABI_SFREE(self%nbcalc_ks)
1940 [ + - ]: 2 : ABI_SFREE(self%bstart_ks)
1941 : : !ABI_SFREE(self%bstop_ks)
1942 [ + - ]: 2 : ABI_SFREE(self%kcalc2ibz)
1943 [ + - ]: 2 : ABI_SFREE(self%coords_kws)
1944 [ + - ]: 2 : ABI_SFREE(self%my_spins)
1945 [ + - ]: 2 : ABI_SFREE(self%my_ikcalc)
1946 [ + - ]: 2 : ABI_SFREE(self%e0vals)
1947 [ + - ]: 2 : ABI_SFREE(self%dw_vals)
1948 [ + - ]: 2 : ABI_SFREE(self%gw_vals)
1949 : : !ABI_SFREE(self%ce_spfunc_wr)
1950 [ + - ]: 2 : ABI_SFREE(self%spfunc_dm_wr)
1951 [ + - ]: 2 : ABI_SFREE(self%spfunc_wr)
1952 [ + - ]: 2 : ABI_SFREE(self%conductivity_mu)
1953 [ + - ]: 2 : ABI_SFREE(self%conductivity_mu_dm)
1954 [ + - ]: 2 : ABI_SFREE(self%mobility_mu)
1955 [ + - ]: 2 : ABI_SFREE(self%mobility_mu_dm)
1956 [ + - ]: 2 : ABI_SFREE(self%transport_mu_e)
1957 [ + - ]: 2 : ABI_SFREE(self%print_dfdw)
1958 [ + - ]: 2 : ABI_SFREE(self%seebeck)
1959 [ + - ]: 2 : ABI_SFREE(self%seebeck_dm)
1960 [ + - ]: 2 : ABI_SFREE(self%kappa)
1961 [ + - ]: 2 : ABI_SFREE(self%kappa_dm)
1962 [ + - ]: 2 : ABI_SFREE(self%l0)
1963 [ + - ]: 2 : ABI_SFREE(self%l0_dm)
1964 [ + - ]: 2 : ABI_SFREE(self%l1)
1965 [ + - ]: 2 : ABI_SFREE(self%l1_dm)
1966 [ + - ]: 2 : ABI_SFREE(self%l2)
1967 [ + - ]: 2 : ABI_SFREE(self%l2_dm)
1968 [ - + ]: 2 : ABI_SFREE(self%time_mesh)
1969 [ - + ]: 2 : ABI_SFREE(self%ct_vals)
1970 [ - + ]: 2 : ABI_SFREE(self%c1)
1971 [ - + ]: 2 : ABI_SFREE(self%c2)
1972 [ - + ]: 2 : ABI_SFREE(self%c3)
1973 [ - + ]: 2 : ABI_SFREE(self%gt_vals)
1974 [ + - ]: 2 : ABI_SFREE(self%wrmesh_b)
1975 [ + - ]: 2 : ABI_SFREE(self%wrmesh_ce)
1976 [ + - ]: 2 : ABI_SFREE(self%vals_e0ks)
1977 [ + - ]: 2 : ABI_SFREE(self%vals_wr)
1978 [ + - ]: 2 : ABI_SFREE(self%kcalc)
1979 [ + - ]: 2 : ABI_SFREE(self%kTmesh)
1980 [ + - ]: 2 : ABI_SFREE(self%mu_e)
1981 [ + - ]: 2 : ABI_SFREE(self%n_ehst)
1982 : :
1983 : 2 : call destroy_mpi_enreg(self%ce_mpi_enreg)
1984 : 2 : call self%ebands%free()
1985 : 2 : call self%spin_comm%free()
1986 : 2 : call self%kcalc_comm%free()
1987 : 2 : call self%wt_comm%free()
1988 : 2 : call self%ncwrite_comm%free()
1989 : :
1990 : 2 : end subroutine cumulant_free
1991 : : !!***
1992 : :
1993 : :
1994 : : integer function gcd(m, n) result(answer)
1995 : : ! The greatest common divisor (GCD) of two nonzero integers
1996 : : ! i.e. the largest positive integer that divides each of the integers.
1997 : : ! gcd(a, 0) = gcd(0, a) = |a|
1998 : : ! gcd(0, 0) is commonly defined as 0.
1999 : :
2000 : : integer,intent(in) :: m, n
2001 : : integer :: irest,ifirst
2002 : :
2003 : : ifirst = iabs(m)
2004 : : answer = iabs(n)
2005 : : if (answer == 0) then
2006 : : answer = ifirst
2007 : : else
2008 : : do
2009 : : irest = mod(ifirst,answer)
2010 : : if (irest == 0) exit
2011 : : ifirst = answer
2012 : : answer = irest
2013 : : end do
2014 : : answer= iabs(answer)
2015 : : end if
2016 : :
2017 : : end function gcd
2018 : :
2019 : : integer function lcm(a, b)
2020 : : integer,intent(in) :: a, b
2021 : : lcm = iabs(a * b) / gcd(a,b)
2022 : : end function lcm
2023 : : !
2024 : : !integer function gcd(a,b)
2025 : : ! integer :: a,b,t
2026 : : ! do while (b/=0)
2027 : : ! t = b
2028 : : ! b = mod(a,b)
2029 : : ! a = t
2030 : : ! end do
2031 : : ! gcd = abs(a)
2032 : : !end function gcd
2033 : :
2034 : : ! Two Factor decomposotionns of positive integer
2035 : : ! The first entry in facts is the largest factor unless is_print is True.
2036 : : ! 1 facts (1, 1) is_prime: True
2037 : : ! 2 facts (2, 1) is_prime: False
2038 : : ! 12 facts (6, 2) is_prime: False
2039 : : ! 11 facts (1, 11) is_prime: True
2040 : : ! 17 facts (1, 17) is_prime: True
2041 : : ! 33 facts (11, 3) is_prime: False
2042 : : ! 35 facts (7, 5) is_prime: False
2043 : :
2044 : 0 : subroutine ifact2(nn, facts, is_prime)
2045 : :
2046 : : integer,intent(in) :: nn
2047 : : integer,intent(out) :: facts(2)
2048 : : logical,intent(out) :: is_prime
2049 : :
2050 : : !Local variables ------------------------------
2051 : : integer :: start, ii
2052 : :
2053 : : ! *************************************************************************
2054 : :
2055 [ # # ]: 0 : ABI_CHECK(nn > 0, sjoin("invalid nn:", itoa(nn)))
2056 : :
2057 : 0 : start = nn / 2 + 1
2058 [ # # ]: 0 : do ii=start, 1, -1
2059 [ # # ]: 0 : if (mod(nn, ii) == 0) exit
2060 : : end do
2061 : :
2062 [ # # ]: 0 : facts = [ii, nn / ii]
2063 : 0 : is_prime = facts(1) == 1
2064 : :
2065 : 0 : end subroutine ifact2
2066 : : !!***
2067 : :
2068 : : ! Invert 3x3 matrix, copied from matr3inv
2069 : 8 : pure subroutine inv33(aa, ait)
2070 : :
2071 : : !Arguments ------------------------------------
2072 : : !arrays
2073 : : real(dp),intent(in) :: aa(3,3)
2074 : : real(dp),intent(out) :: ait(3,3)
2075 : :
2076 : : !Local variables-------------------------------
2077 : : !scalars
2078 : : real(dp) :: dd,det,t1,t2,t3
2079 : :
2080 : : ! *************************************************************************
2081 : :
2082 : 8 : t1 = aa(2,2) * aa(3,3) - aa(3,2) * aa(2,3)
2083 : 8 : t2 = aa(3,2) * aa(1,3) - aa(1,2) * aa(3,3)
2084 : 8 : t3 = aa(1,2) * aa(2,3) - aa(2,2) * aa(1,3)
2085 : 8 : det = aa(1,1) * t1 + aa(2,1) * t2 + aa(3,1) * t3
2086 : :
2087 : : ! Make sure matrix is not singular
2088 [ + + ]: 8 : if (abs(det) > 100 * tiny(one)) then
2089 : 6 : dd = one / det
2090 : 6 : ait(1,1) = t1 * dd
2091 : 6 : ait(2,1) = t2 * dd
2092 : 6 : ait(3,1) = t3 * dd
2093 : 6 : ait(1,2) = (aa(3,1)*aa(2,3)-aa(2,1)*aa(3,3)) * dd
2094 : 6 : ait(2,2) = (aa(1,1)*aa(3,3)-aa(3,1)*aa(1,3)) * dd
2095 : 6 : ait(3,2) = (aa(2,1)*aa(1,3)-aa(1,1)*aa(2,3)) * dd
2096 : 6 : ait(1,3) = (aa(2,1)*aa(3,2)-aa(3,1)*aa(2,2)) * dd
2097 : 6 : ait(2,3) = (aa(3,1)*aa(1,2)-aa(1,1)*aa(3,2)) * dd
2098 : 6 : ait(3,3) = (aa(1,1)*aa(2,2)-aa(2,1)*aa(1,2)) * dd
2099 [ + + + + : 156 : ait = transpose(ait)
+ + + + ]
2100 : : else
2101 : 2 : ait = zero
2102 : : end if
2103 : :
2104 : 8 : end subroutine inv33
2105 : :
2106 : :
2107 [ # # # # : 0 : end module m_cumulant
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ]
2108 : : !!***
|
