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