Line data Source code
1 : !****m* ABINIT/m_paw_dmft
2 : !! NAME
3 : !! m_paw_dmft
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2006-2026 ABINIT group (BAmadon)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! INPUTS
14 : !!
15 : !! OUTPUT
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : MODULE m_paw_dmft
26 :
27 : use defs_basis
28 : use m_abicore
29 : use m_CtqmcInterface
30 : use m_data4entropyDMFT
31 : use m_dtset
32 : use m_errors
33 : use m_xmpi
34 :
35 : use defs_abitypes, only : MPI_type
36 : use m_crystal, only : crystal_t
37 : use m_fstrings, only : int2char4
38 : use m_geometry, only : symredcart
39 : use m_io_tools, only : open_file
40 : use m_mpinfo, only : proc_distrb_cycle
41 : use m_paw_numeric, only : paw_jbessel_4spline
42 : use m_pawang, only : pawang_type
43 : use m_pawrad, only : pawrad_free,pawrad_init,pawrad_type,simp_gen
44 : use m_pawtab, only : pawtab_type
45 :
46 : implicit none
47 :
48 : private
49 :
50 : public :: init_dmft
51 : public :: init_sc_dmft
52 : public :: construct_nwli_dmft
53 : public :: destroy_dmft
54 : public :: destroy_sc_dmft
55 : public :: print_dmft
56 : public :: print_sc_dmft
57 : public :: saveocc_dmft
58 : public :: readocc_dmft
59 :
60 : !!***
61 :
62 : !----------------------------------------------------------------------
63 :
64 : !!****t* m_paw_dmft/mpi_distrib_dmft_type
65 : !! NAME
66 : !! mpi_distrib_dmft_type
67 : !!
68 : !! FUNCTION
69 : !! This structured datatype contains the necessary data for the MPI
70 : !! parallelization over frequencies and kpts in DFT+DMFT.
71 : !!
72 : !! SOURCE
73 :
74 : type, public :: mpi_distrib_dmft_type
75 :
76 : ! Two types of parallelization
77 : ! Type 1: parallelization over kpt and then frequencies
78 : ! Type 2: parallelization over frequencies only
79 :
80 : integer :: comm_freq
81 : ! Frequency communicator (type 1)
82 :
83 : integer :: comm_kpt
84 : ! Kpt communicator (type 1)
85 :
86 : integer :: me_freq
87 : ! Rank in comm_freq (type 1)
88 :
89 : integer :: me_kpt
90 : ! Rank in comm_kpt (type 1)
91 :
92 : integer :: nw
93 : ! Number of frequencies (either nwlo or nwr, same for both types)
94 :
95 : integer :: shiftk
96 : ! Shift from kpt index on the current CPU to the physical index (type 1)
97 :
98 : integer, allocatable :: nkpt_mem(:)
99 : ! Number of kpt handled by each CPU of the kpt communicator (type 1)
100 :
101 : integer, allocatable :: nw_mem(:)
102 : ! Number of frequencies handled by each CPU of the global communicator (type 2)
103 :
104 : integer, allocatable :: nw_mem_kptparal(:)
105 : ! Number of frequencies handled by each CPU of the frequency communicator (type 1)
106 :
107 : integer, allocatable :: procb(:)
108 : ! Rank in comm_kpt of the CPU handling each kpt (type 1)
109 :
110 : integer, allocatable :: procf(:)
111 : ! Rank in the global communicator of the CPU handling each frequency (type 2)
112 :
113 : integer, allocatable :: proct(:)
114 : ! Rank in comm_freq of the CPU handling each frequency (type 1)
115 :
116 : end type mpi_distrib_dmft_type
117 : !!***
118 :
119 : !----------------------------------------------------------------------
120 :
121 : !!****t* m_paw_dmft/paw_dmft_type
122 : !! NAME
123 : !! paw_dmft_type
124 : !!
125 : !! FUNCTION
126 : !! This structured datatype contains the necessary data for the link
127 : !! between dmft and paw.
128 : !! occnd(non-diagonal band occupations for self-consistency), band_in
129 : !! (say which band are taken into account in the calculation), and the
130 : ! dimensions of these arrays.
131 : !!
132 : !! SOURCE
133 :
134 : type, public :: paw_dmft_type
135 :
136 : integer :: dmft_blockdiag
137 : ! Block diagonalize Hamiltonian in the local basis
138 :
139 : integer :: dmft_dc
140 : ! Type of double counting used in DMFT
141 :
142 : integer :: dmft_entropy
143 : ! = 0: do not compute entropy
144 : ! >= 1: compute entropy with an integration over self-consistent calculations
145 :
146 : integer :: dmft_hybri_limit
147 : ! = 1 : Use asymptotic limit to build hybridization function
148 :
149 : integer :: dmft_iter
150 : ! Nb of iterations for DMFT self-consistent cycle.
151 :
152 : integer :: dmft_kspectralfunc
153 : ! =0 Default
154 : ! =1 Activate calculation of k-resolved spectral function
155 :
156 : integer :: dmft_log_freq
157 : ! = 0: do not use log frequencies
158 : ! = 1: use log frequencies
159 :
160 : integer :: dmft_magnfield
161 : ! = 0: do nothing
162 : ! = 1: apply a magnetic field Bz via Zeeman Hamiltonian on Kohn-Sham energies
163 : ! = 2: apply a magnetic field Bz via Zeeman Hamiltonian on local impurity Hamiltonian
164 :
165 : integer :: dmft_nwli
166 : ! Physical index of the last imaginary frequency (/=dmft_nwlo when dmft_log_freq=1)
167 :
168 : integer :: dmft_nwlo
169 : ! Number of imaginary frequencies
170 :
171 : integer :: dmft_nwr
172 : ! Number of real frequencies
173 :
174 : integer :: dmft_prgn
175 : ! Specify the way of printing the green function.
176 : ! =1 print green
177 : ! =2 print self
178 :
179 : integer :: dmft_prt_maxent
180 : ! =1 to print Maxent files
181 :
182 : integer :: dmft_prtself
183 : ! =1 to keep self-energy files of all previous iterations
184 :
185 : integer :: dmft_prtwan
186 : ! =1 to print Wannier functions
187 :
188 : integer :: dmft_read_occnd
189 : ! Flag to read/write DMFT occupations
190 : ! =0 : Occupations are written but not read
191 : ! =1 : Occupations are read from I_DMFTOCCND, with I the root for input file
192 : ! =2 : Occupations are read from O_DMFTOCCND, with O the root for input file
193 :
194 : integer :: dmft_rslf
195 : ! Flag to read the self-energy at each iteration
196 : ! =-1 : Self-energy is set to 0
197 : ! =0 : Self-energy is set to double counting
198 : ! =1 : Self-energy is read from file
199 :
200 : integer :: dmft_solv
201 : ! Choice of solver for DMFT
202 :
203 : integer :: dmft_t2g
204 : ! Only use t2g orbitals
205 :
206 : integer :: dmft_triqs_chiloc
207 : ! Activate measurements for local spin-spin correlation function with triqs cthyb
208 : ! only relevant without spin-orbit coupling (nspinor = 1)
209 :
210 : integer :: dmft_triqs_chiloc_ins
211 : ! Option for local spin-spin correlation function measurements
212 : ! default value is 10
213 :
214 : integer :: dmft_triqs_compute_integral
215 : ! Only relevant when dmft_triqs_entropy=1.
216 : ! =1: Compute the impurity entropy by thermodynamic integration over interaction strength.
217 : ! =0: Do not compute the impurity entropy. All the other contributions to the free
218 : ! energy are still computed.
219 :
220 : integer :: dmft_triqs_det_init_size
221 : ! TRIQS CTQMC: Initial size of the hybridization matrix. If it is too low,
222 : ! the matrix will be resized very often, which can be slow.
223 :
224 : integer :: dmft_triqs_det_n_operations_before_check
225 : ! TRIQS CTQMC: Number of operations before check of the hybridization matrix.
226 : ! If it is low, the matrix will be checked too often, which can be slow.
227 :
228 : integer :: dmft_triqs_entropy
229 : ! TRIQS CTQMC: Compute the DMFT entropy.
230 :
231 : integer :: dmft_triqs_gaussorder
232 : ! Order of the Gauss-Legendre quadrature for each subdivision of the thermodynamic integration.
233 :
234 : integer :: dmft_triqs_loc_n_min
235 : ! TRIQS CTQMC: Only configurations with a number of electrons in
236 : ! [nlocmin,nlocmax] are taken into account.
237 :
238 : integer :: dmft_triqs_loc_n_max
239 : ! TRIQS CTQMC: Only configurations with a number of electrons in
240 : ! [nlocmin,nlocmax] are taken into account.
241 :
242 : integer :: dmft_triqs_n_cycles
243 : ! TRIQS CTQMC: Number of measurement cycles.
244 :
245 : integer :: dmft_triqs_nleg
246 : ! TRIQS CTQMC: Nb of Legendre polynomials used for the
247 : ! Green's function (Phys. Rev. B 84, 075145) [[cite:Boehnke2011]].
248 :
249 : integer :: dmft_triqs_nsubdivisions
250 : ! Number of regular subdivisions of the interval [0,U], each of which
251 : ! containing dmft_triqs_gaussorder points
252 :
253 : integer :: dmft_triqs_read_ctqmcdata
254 : ! TRIQS CTQMC: Read CTQMC data of the previous iteration
255 :
256 : integer :: dmft_triqs_seed_a
257 : ! TRIQS CTQMC: The CTQMC seed is seed_a + rank * seed_b.
258 :
259 : integer :: dmft_triqs_seed_b
260 : ! TRIQS CTQMC: The CTQMC seed is seed_a + rank * seed_b.
261 :
262 : integer :: dmft_triqs_therm_restart
263 : ! TRIQS CTQMC: Number of thermalization steps when we restart from a previous configuration.
264 :
265 : integer :: dmft_full_chipsi
266 : ! =0 do not use
267 : ! =1 build Wannier functions
268 :
269 : integer :: dmft_wanorthnorm
270 : ! =2 orthonormalization of Wannier functions for each k-point
271 : ! =3 orthonormalization over the sum over k-points
272 :
273 : integer :: dmft_x2my2d
274 : ! Only use x2my2d orbital
275 :
276 : integer :: dmftbandf
277 : ! Highest correlated band
278 :
279 : integer :: dmftbandi
280 : ! Lowest correlated band
281 :
282 : integer :: dmftcheck
283 : ! Check various part of the implementation
284 :
285 : integer :: dmftctqmc_basis
286 : ! Basis in which to perform the CTQMC calculation
287 : ! 0 : Slm basis, 1 : diagonalize local Hamiltonian, 2: diagonalize the density matrix
288 : ! Only for TRIQS: 3: Ylm, 4: JmJ
289 :
290 : integer :: dmftctqmc_check
291 : ! ABINIT CTQMC: perform a check on the impurity and/or bath operator
292 : ! only for debug
293 : ! 0 : nothing, 1 : impurity, 2 : bath, 3 : both
294 :
295 : integer :: dmftctqmc_correl
296 : ! ABINIT CTQMC: Gives analysis for CTQMC
297 : ! 0 : nothing, 1 : activated Correlations.dat
298 :
299 : integer :: dmftctqmc_gmove
300 : ! ABINIT CTQMC: add global move every dmftctqmc_gmove sweeps
301 : ! >= 0 ; done inside CT-QMC with warning
302 : ! == 0 ; no global moves
303 :
304 : integer :: dmftctqmc_grnns
305 : ! ABINIT CTQMC: compute green function noise for each imaginary time
306 : ! 0 : nothing, 1 : activated
307 :
308 : integer :: dmftctqmc_localprop
309 : ! ABINIT CTQMC: local properties calculations
310 : ! 0 : nothing, 1 : Histogram, 2 : magnetic susceptibility, 3 : charge susceptibility
311 :
312 : integer :: dmftctqmc_meas
313 : ! ABINIT/TRIQS CTQMC: measurements are done every dmftctqmc_meas step
314 :
315 : integer :: dmftctqmc_mov
316 : ! ABINIT CTQMC: Gives movie for CTQMC
317 : ! 0 : nothing, 1 : 1 file Movie_RANK.tex for each cpu
318 :
319 : integer :: dmftctqmc_chains
320 : ! ABINIT CTQMC: Number of chains per MPI-tasks
321 :
322 : integer :: dmftctqmc_mrka
323 : ! ABINIT CTQMC: Write a temporary file Spectra_RANK.dat with the sweep evolution of
324 : ! the number of electron for each flavor
325 : ! The measurement is done every dmftctqmc_meas*dmftctqmc_mrka sweep
326 : ! e.g. : meas=2 mrka=10 -> every 20 sweeps sum_i c+(ti)c(t'i) is measured
327 :
328 : integer :: dmftctqmc_order
329 : ! ABINIT CTQMC: Gives perturbation order of CTQMC solver
330 : ! 0 : nothing, >=1 max order evaluated in Perturbation.dat
331 :
332 : integer :: dmftqmc_l
333 : ! Number of points on the imaginary time grid for G(tau) and Delta(tau)
334 :
335 : ! integer :: dmft_mag
336 : ! ! 0 if non magnetic calculation, 1 if magnetic calculation
337 :
338 : integer :: dmftqmc_seed
339 : ! Seed for CTQMC (only for ABINIT)
340 :
341 : integer :: dmftqmc_therm
342 : ! Number of thermalization steps for CTQMC (only for ABINIT, and for TRIQS when we don't restart from a previous configuration)
343 :
344 : integer :: gpu_option
345 : ! Wether to use GPU implementation (expected values: ABI_GPU_DISABLED, ABI_GPU_OPENMP)
346 :
347 : integer :: idmftloop
348 : ! Current iteration in the DFT+DMFT loop
349 :
350 : integer :: ientropy
351 : ! activate evaluation of terms for alternative calculation of entropy in DMFT
352 :
353 : integer :: ireadctqmcdata
354 : ! Internal flag to indicate if an input CTQMC_DATA file must be read
355 :
356 : integer :: ireadself
357 : ! Internal flag to indicate if an input self file must be read
358 :
359 : integer :: ixc
360 : ! Index of the XC functional
361 :
362 : integer :: kptopt
363 : ! Option to generate kpts
364 :
365 : integer :: lchipsiortho
366 : ! Internal flag
367 : ! =0 <Chi|Psi> is not orthonormalized
368 : ! =1 <Chi|Psi> is orthonormalized
369 :
370 : integer :: maxlpawu
371 : ! Maximal correlated l over all atoms
372 :
373 : integer :: maxmeshsize
374 : ! Maximal size of the radial mesh over all atoms
375 :
376 : integer :: maxnproju
377 : ! Maximal number of correlated projectors over all atoms
378 :
379 : integer :: mband
380 : ! Total number of bands
381 :
382 : integer :: mbandc
383 : ! Total number of correlated bands
384 :
385 : integer :: mkmem
386 : ! Number of k-points handled by the current process within the DFT
387 : ! parallelization scheme
388 :
389 : integer :: myproc
390 : ! Rank in the global communicator
391 :
392 : integer :: natom
393 : ! Number of atoms
394 :
395 : integer :: natpawu
396 : ! Number of correlated atoms
397 :
398 : integer :: nkpt
399 : ! Number of k-points in the IBZ.
400 :
401 : !integer :: nspden
402 : ! Number of spin densities
403 :
404 : integer :: nproc
405 : ! Total number of MPI processes
406 :
407 : integer :: nspinor
408 : ! Number of spinor components
409 :
410 : integer :: nsppol
411 : ! Number of spin polarizations
412 :
413 : integer :: nsym
414 : ! Number of symmetries
415 :
416 : integer :: ntypat
417 : ! Number of atom types
418 :
419 : integer :: prtdos
420 : ! Print DOS when >=1
421 :
422 : integer :: prtvol
423 : ! Flag for different print options
424 :
425 : integer :: spacecomm
426 : ! MPI_COMM_WORLD
427 :
428 : integer :: unpaw
429 : ! File number for cprj
430 :
431 : integer :: use_dmft
432 : ! 1 if non diagonal occupations are used, else 0
433 :
434 : integer :: use_fixed_self
435 : ! Impose a fixed self-energy during the first use_fixed_self iterations
436 :
437 : integer :: use_sc_dmft
438 : ! 1 for charge-self consistent calculations
439 :
440 : logical :: dmft_triqs_leg_measure
441 : ! TRIQS CTQMC: Flag to activate Legendre measurement
442 :
443 : logical :: dmft_triqs_measure_density_matrix
444 : ! TRIQS CTQMC: Flag to activate the measurement of the density matrix
445 :
446 : logical :: dmft_triqs_move_double
447 : ! TRIQS CTQMC: Flag to activate the double moves
448 :
449 : logical :: dmft_triqs_move_shift
450 : ! TRIQS CTQMC: Flag to activate the shift move
451 :
452 : logical :: dmft_triqs_off_diag
453 : ! TRIQS CTQMC: Flag to sample the off-diagonal elements of the Green's function
454 :
455 : logical :: dmft_triqs_prt_entropy
456 : ! TRIQS CTQMC: Flag to print more information during entropy calculation
457 :
458 : logical :: dmft_triqs_time_invariance
459 : ! TRIQS CTQMC: Flag to activate the use of time invariance for the sampling
460 : ! of the density matrix
461 :
462 : logical :: dmft_triqs_use_norm_as_weight
463 : ! TRIQS CTQMC: Flag to activate the use of the norm of the matrix as weight
464 : ! instead of the trace
465 :
466 : real(dp) :: dmft_charge_prec
467 : ! Precision on charge required for determination of fermi level (fermi_green)
468 :
469 : real(dp) :: dmft_fermi_prec
470 : ! Required precision on Fermi level (fermi_green) during the DMFT SCF cycle, (=> ifermie_cv)
471 : ! used also for self (new_self) (=> iself_cv).
472 :
473 : real(dp) :: dmft_fermi_step
474 : ! When dmft_optim = 0, step increment to find the upper and lower bounds of the Fermi level
475 : ! When dmft_optim = 1, maximal step size in the Fermi level search
476 :
477 : real(dp) :: dmft_lcpr
478 : ! Required precision on local correlated charge in order to stop SCF
479 : ! DMFT cycle (integrate_green) => ichargeloc_cv
480 :
481 : real(dp) :: dmft_magnfield_b
482 : ! Value of the applied magnetic field in Tesla
483 :
484 : real(dp) :: dmft_mxsf
485 : ! Mixing coefficient for Self-Energy during the SCF DMFT cycle.
486 :
487 : real(dp) :: dmft_tolfreq
488 : ! Required precision on local correlated density matrix (depends on
489 : ! frequency mesh), used in m_dmft/dmft_solve
490 :
491 : real(dp) :: dmft_triqs_det_precision_error
492 : ! TRIQS CTQMC: Error threshold for the deviation of the determinant when a check is performed.
493 :
494 : real(dp) :: dmft_triqs_det_precision_warning
495 : ! TRIQS CTQMC: Warning threshold for the deviation of the determinant when a check is performed.
496 :
497 : real(dp) :: dmft_triqs_det_singular_threshold
498 : ! TRIQS CTQMC: Threshold when checking if the determinant is singular.
499 :
500 : real(dp) :: dmft_triqs_epsilon
501 : ! TRIQS CTQMC: Threshold for singular values of the kernel matrix for the DLR fit
502 :
503 : real(dp) :: dmft_triqs_imag_threshold
504 : ! TRIQS CTQMC: Threshold for the imaginary part of Delta(tau)
505 :
506 : real(dp) :: dmft_triqs_lambda
507 : ! TRIQS CTQMC: Cutoff for the real frequency grid for the DLR fit
508 :
509 : real(dp) :: dmft_triqs_pauli_prob
510 : ! TRIQS CTQMC: Probability for proposing Pauli-aware insert and remove
511 :
512 : real(dp) :: dmft_triqs_shift_mu
513 : ! TRIQS CTQMC: Shift of the chemical potential for thermodynamic integration
514 :
515 : real(dp) :: dmft_triqs_tol_block
516 : ! TRIQS CTQMC: Off-diagonal elements below this threshold are set to 0
517 :
518 : real(dp) :: dmft_wanrad
519 : ! Maximal radius for print of the Wannier functions
520 :
521 : real(dp) :: dmftqmc_n
522 : ! ABINIT CTQMC: Nb of sweeps
523 :
524 : real(dp) :: e_dc
525 : ! Double counting energy
526 :
527 : real(dp) :: e_hu
528 : ! Interaction energy
529 :
530 : real(dp) :: fermie
531 : ! DMFT Fermi level
532 :
533 : real(dp) :: fermie_dft
534 : ! DFT Fermi level
535 :
536 : real(dp) :: j_for_s
537 : ! Variable for evaluation of correlation energy for U=0 in the entropic
538 : ! calculation
539 :
540 : real(dp) :: nelectval
541 : ! Number of valence electrons
542 :
543 : real(dp) :: sdmft
544 : ! DFT+DMFT total entropy
545 :
546 : real(dp) :: simp
547 : ! DFT+DMFT entropy of the impurity electrons
548 :
549 : real(dp) :: temp
550 : ! Temperature (Ha)
551 :
552 : real(dp) :: u_for_s
553 : ! Variable for evaluation of correlation energy for U=0 in the entropic
554 : ! calculation
555 :
556 : character(len=fnlen) :: filapp
557 : ! Output file name
558 :
559 : character(len=fnlen) :: filctqmcdatain
560 : ! Input file name for CTQMC_DATA file
561 :
562 : character(len=fnlen) :: filnamei
563 : ! Input file name
564 :
565 : character(len=fnlen) :: filselfin
566 : ! Input file name for self file
567 :
568 : integer, allocatable :: bandc_proc(:)
569 : ! Proc index (on comm_band) for each correlated band in DMFT (for kgb paral)
570 :
571 : integer, allocatable :: exclude_bands(:)
572 : ! Gives the bands than are not in the DMFT calculations.
573 :
574 : integer, allocatable :: include_bands(:)
575 : ! For each bands included in the calculation (1..mbandc), include_bands
576 : ! gives the index in the full band index (1...mband)
577 :
578 : integer, allocatable :: lpawu(:)
579 : ! Correlated l for each atom (set to -1 if not correlated)
580 :
581 : integer, allocatable :: siz_proj(:)
582 : ! Size of the radial mesh for the DMFT orbital, for each atom type.
583 :
584 : logical, allocatable :: band_in(:)
585 : ! True for each band included in the calculation
586 :
587 : logical, allocatable :: use_bandc(:)
588 : ! True for each proc wich has at least one band involved in DMFT non diagonal
589 : ! occupations on band parallelism
590 :
591 : real(dp), allocatable :: edc(:)
592 : ! Double counting energy for each atom (only used as a temporary for dmft_dc=8)
593 :
594 : real(dp), allocatable :: edcdc(:)
595 : ! Integral of Vdc * rho for each atom (only used as a temporary for dmft_dc=8)
596 :
597 : real(dp), allocatable :: eigen_dft(:,:,:)
598 : ! DFT eigenvalues for each correlated band, k-point, polarization
599 :
600 : real(dp), allocatable :: occnd(:,:,:,:,:)
601 : ! Non diagonal band-occupation for each k-point, polarization.
602 :
603 : real(dp), allocatable :: omega_lo(:)
604 : ! Imaginary frequencies
605 :
606 : real(dp), allocatable :: omega_r(:)
607 : ! Real frequencies
608 :
609 : real(dp), allocatable :: phi_int(:,:)
610 : ! Integral of <Chi|Phi> for every correlated projector and atom type
611 :
612 : real(dp), allocatable :: phimtphi(:,:,:)
613 : ! Phi-Phi_tilde for every r,correlated projector and atom type
614 :
615 : real(dp), allocatable :: phimtphi_int(:,:)
616 : ! Integral of <Chi|Phi-Phi_tilde> for every correlated projector and atom type
617 :
618 : real(dp), allocatable :: symrec_cart(:,:,:)
619 : ! Symmetries in cartesian coordinates
620 :
621 : real(dp), allocatable :: wgt_wlo(:)
622 : ! Weight of the imaginary frequencies
623 :
624 : real(dp), allocatable :: ylm(:,:,:,:)
625 : ! Ylm(k+G) for each G,m,l,k
626 :
627 : ! real(dp), allocatable :: phi0phiiint(:)
628 : ! ! non diagonal band-occupation for each k-point, polarisation.
629 :
630 : complex(dp), allocatable :: bessel(:,:,:,:)
631 : ! 4*pi*(i**l)*jl(|k+G|r)*r/sqrt(ucvol) for each G,r,atom type and kpt
632 :
633 : complex(dp), allocatable :: bessel_int(:,:,:)
634 : ! Integral over r of bessel(ig,:,iat,ikpt)
635 :
636 : complex(dp), allocatable :: buf_psi(:)
637 : ! Temporary buffer for the computation of Wannier function
638 :
639 : complex(dp), allocatable :: chipsi(:,:,:,:,:)
640 : ! Hermitian product <Chi|Psi> for each flavor, correlated band,
641 : ! k-point, polarization and atom
642 :
643 : complex(dp), allocatable :: dpro(:,:,:)
644 : ! Exp(i(k+G).xred(iatom)) for each G,correlated atom and k
645 :
646 : complex(dp), allocatable :: jmj2ylm(:,:,:)
647 : ! Transformation matrix from JmJ to Ylm basis for each lpawu
648 :
649 : complex(dp), allocatable :: slm2ylm(:,:,:)
650 : ! Transformation matrix from real to complex harmonics for each lpawu
651 :
652 : complex(dp), allocatable :: wannier(:,:,:)
653 : ! Wannier functions for each r,flavor and atom
654 :
655 : complex(dp), allocatable :: zarot(:,:,:,:)
656 : ! Coeffs of the transformation of real spherical
657 : ! harmonics under the symmetry operations symrec.
658 :
659 : integer, ABI_CONTIGUOUS pointer :: dmft_nominal(:) => null()
660 : ! Only relevant when dmft_dc=7. Nominal occupancies for each atom.
661 :
662 : integer, pointer :: indsym(:,:) => null()
663 : ! Label of atom into which iatom is sent by the INVERSE of the
664 : ! symmetry operation symrel(isym)
665 :
666 : integer, pointer :: int_meshsz(:) => null()
667 : ! PAW integration radius for each atom type
668 :
669 : integer, ABI_CONTIGUOUS pointer :: nband(:) => null()
670 : ! Number of bands for each k-point and polarization
671 :
672 : integer, ABI_CONTIGUOUS pointer :: npwarr(:) => null()
673 : ! Number of plane waves on current process for each k-point
674 :
675 : integer, ABI_CONTIGUOUS pointer :: typat(:) => null()
676 : ! Type of each atom
677 :
678 : real(dp), ABI_CONTIGUOUS pointer :: dmft_shiftself(:) => null()
679 : ! Initial shift of the self-energy for each atom
680 :
681 : real(dp), ABI_CONTIGUOUS pointer :: eigen(:) => null()
682 : ! DFT eigenvalues
683 :
684 : real(dp), pointer :: fixed_self(:,:,:,:) => null()
685 : ! Fixed self-energy (only used when use_fixed_self > 0)
686 :
687 : real(dp), ABI_CONTIGUOUS pointer :: wtk(:) => null()
688 : ! Weights for each k-point
689 :
690 : type(CtqmcInterface), allocatable :: hybrid(:)
691 :
692 : type(data4entropyDMFT_t) :: forentropyDMFT
693 :
694 : type(pawrad_type), allocatable :: radgrid(:)
695 : ! Radial grid for each type of atom
696 :
697 : type(mpi_distrib_dmft_type) :: distrib
698 : ! MPI parallelization for imaginary frequencies
699 :
700 : type(mpi_distrib_dmft_type) :: distrib_r
701 : ! MPI parallelization for real frequencies
702 :
703 : end type paw_dmft_type
704 : !!***
705 :
706 : !----------------------------------------------------------------------
707 :
708 : CONTAINS !========================================================================================
709 : !!***
710 :
711 : !!****f* m_paw_dmft/init_sc_dmft
712 : !! NAME
713 : !! init_sc_dmft
714 : !!
715 : !! FUNCTION
716 : !! Allocate variables used in type paw_dmft_type.
717 : !!
718 : !! INPUTS
719 : !! dtset <type(dataset_type)>=all input variables for this dataset
720 : !! mpsang = highest angular momentum + 1
721 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
722 : !! gprimd(3,3) = dimensional reciprocal space primitive translations
723 : !! kg(3,mpw*mkmem) = reduced planewave coordinates.
724 : !! mpi_enreg = information about MPI parallelization
725 : !! npwarr(nkpt) = number of planewaves in basis at this k point
726 : !! occ = DFT occupations
727 : !! pawang <type(pawang)>=paw angular mesh and related data
728 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data
729 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
730 : !! rprimd(3,3) = dimensional primitive vectors
731 : !! ucvol = unit cell volume in bohr**3.
732 : !! unpaw = unit number for temporary PAW data
733 : !! use_sc_dmft = for charge self-consistent calculations
734 : !! xred(3,natom) = reduced dimensionless atomic coordinates
735 : !! ylm(mpw*mkmem,mpsang*mpsang*useylm) = real spherical harmonics for each G and k point
736 : !!
737 : !! OUTPUTS
738 : !! paw_dmft = datastructure for dmft
739 : !!
740 : !! SOURCE
741 :
742 5629 : subroutine init_sc_dmft(dtset,mpsang,paw_dmft,gprimd,kg,mpi_enreg,npwarr,occ,pawang, &
743 11258 : & pawrad,pawtab,rprimd,ucvol,unpaw,use_sc_dmft,xred,ylm)
744 :
745 : !Arguments ------------------------------------
746 : !scalars
747 : integer, intent(in) :: mpsang
748 : integer, optional, intent(in) :: unpaw,use_sc_dmft
749 : real(dp), optional, intent(in) :: ucvol
750 : !type
751 : type(paw_dmft_type), intent(inout) :: paw_dmft
752 : type(MPI_type), optional, intent(in) :: mpi_enreg
753 : type(dataset_type), target, intent(in) :: dtset
754 : type(pawtab_type), optional, intent(in) :: pawtab(dtset%ntypat)
755 : type(pawang_type), optional, intent(in) :: pawang
756 : type(pawrad_type), target, optional, intent(in) :: pawrad(dtset%ntypat)
757 : ! arrays
758 : real(dp), optional, intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
759 : integer, target, optional, intent(in) :: npwarr(dtset%nkpt)
760 : integer, optional, intent(in) :: kg(3,dtset%mpw*dtset%mkmem)
761 : real(dp), optional, intent(in) :: gprimd(3,3),rprimd(3,3),xred(3,dtset%natom)
762 : real(dp), optional, intent(in) :: ylm(dtset%mpw*dtset%mkmem,mpsang*mpsang)
763 : !Local variables ------------------------------------
764 : integer :: bdtot_index,dmft_dc,dmft_solv,dmftbandi,dmftbandf,fac,i,iatom
765 : integer :: iatom1,iband,icb,ig,ik,ikg,ikpt,im,im1,indproj,iproj,ir,isppol
766 : integer :: itypat,jc1,jj,jm,ll,lpawu,lpawu1,maxlpawu,mband,mbandc,mesh_size,mesh_type
767 : integer :: mkmem,ml1,mm,mpw,ms1,myproc,natom,nband_k,ndim,nkpt,nproc,nproju,npw
768 : integer :: nspinor,nsppol,nsym,ntypat,nwli,off_diag,siz_paw,siz_proj,siz_wan,use_dmft
769 : logical :: t2g,use_full_chipsi,verif,x2my2d
770 : real(dp) :: bes,besp,fac_bessel,invsqrt2lp1,lstep,norm,onem,rad,rint,rstep,sumwtk,xj,xmj
771 : complex(dpc) :: j_l
772 : integer, parameter :: mt2g(3) = (/1,2,4/)
773 5679 : integer, allocatable :: ind_msml(:,:)
774 5679 : logical, allocatable :: lcycle(:),typcycle(:)
775 5679 : real(dp), allocatable :: rmax(:),kpg(:,:),kpg_norm(:)
776 : character(len=500) :: dc_string,lda_string,message
777 : !************************************************************************
778 :
779 5679 : mband = dtset%mband
780 5679 : mkmem = dtset%mkmem
781 5679 : natom = dtset%natom
782 5679 : nkpt = dtset%nkpt
783 5679 : nspinor = dtset%nspinor
784 5679 : nsppol = dtset%nsppol
785 5679 : nsym = dtset%nsym
786 5679 : ntypat = dtset%ntypat
787 :
788 5679 : dmftbandi = dtset%dmftbandi
789 5679 : dmftbandf = dtset%dmftbandf
790 5679 : dmft_dc = dtset%dmft_dc
791 5679 : dmft_solv = dtset%dmft_solv
792 5679 : off_diag = dtset%dmft_triqs_off_diag
793 5679 : use_dmft = abs(dtset%usedmft)
794 5679 : paw_dmft%use_dmft = use_dmft
795 5679 : paw_dmft%use_sc_dmft = 0
796 :
797 5679 : paw_dmft%dmftbandf = dmftbandf
798 5679 : paw_dmft%dmftbandi = dmftbandi
799 5679 : paw_dmft%mband = mband
800 5679 : paw_dmft%mkmem = mkmem
801 5679 : paw_dmft%nkpt = nkpt
802 5679 : paw_dmft%nsym = nsym
803 5679 : paw_dmft%ntypat = ntypat
804 :
805 : ! Spin related variables
806 5679 : paw_dmft%nsppol = nsppol
807 5679 : paw_dmft%nspinor = nspinor
808 5679 : paw_dmft%idmftloop = 0
809 5679 : paw_dmft%mbandc = 0
810 : !paw_dmft%nspden = nspden
811 :
812 5679 : paw_dmft%dmft_read_occnd = dtset%dmft_read_occnd
813 :
814 5679 : if(use_dmft == 10) then
815 0 : ABI_MALLOC(paw_dmft%occnd,(2,mband,mband,nkpt,nsppol*1))
816 0 : ABI_MALLOC(paw_dmft%band_in,(mband*1))
817 0 : ABI_MALLOC(paw_dmft%include_bands,((dmftbandf-dmftbandi+1)*1))
818 0 : ABI_MALLOC(paw_dmft%exclude_bands,(mband*1))
819 : else
820 34074 : ABI_MALLOC(paw_dmft%occnd,(2,mband,mband,nkpt,nsppol*use_dmft))
821 17037 : ABI_MALLOC(paw_dmft%band_in,(mband*use_dmft))
822 17037 : ABI_MALLOC(paw_dmft%include_bands,((dmftbandf-dmftbandi+1)*use_dmft))
823 11358 : ABI_MALLOC(paw_dmft%exclude_bands,(mband*use_dmft))
824 : endif
825 :
826 5679 : if (use_dmft == 0) return
827 :
828 : ! In the case where no kpt is treated by the current CPU (ie sum(isppoltab)=0),
829 : ! dtset%mkmem is set to nkpt by convention in Abinit. We set it to its true value 0.
830 114 : if (sum(mpi_enreg%my_isppoltab(1:nsppol)) == 0) mkmem = 0
831 :
832 : ! Check processors for DMFT
833 : ! Initialize spaceComm, myproc, and nproc
834 : !spacecomm=mpi_enreg%comm_cell
835 : !myproc=mpi_enreg%me_cell
836 : !nproc=mpi_enreg%nproc_cell
837 : !spacecomm = mpi_enreg%comm_world
838 43 : myproc = mpi_enreg%me
839 43 : nproc = mpi_enreg%nproc
840 : !print *, " spacecomm,myproc,nproc",spacecomm,myproc,nproc
841 43 : paw_dmft%spacecomm = mpi_enreg%comm_world
842 43 : paw_dmft%myproc = myproc
843 43 : paw_dmft%nproc = nproc
844 :
845 43 : paw_dmft%unpaw = unpaw
846 :
847 43 : if (dtset%nbandkss == 0) paw_dmft%use_sc_dmft = use_sc_dmft
848 :
849 : ! Do not comment these lines: it guarantees the parallelism in DMFT/QMC will work.
850 43 : if (xmpi_comm_size(xmpi_world) /= xmpi_comm_size(mpi_enreg%comm_world)) &
851 0 : & ABI_ERROR("Someone changed the k-point parallelism again")
852 :
853 43 : if (dmft_solv == 0) then
854 3 : do itypat=1,ntypat
855 2 : if (pawtab(itypat)%lpawu == -1) cycle
856 2 : if ((pawtab(itypat)%upawu > tol5) .or. (pawtab(itypat)%jpawu > tol5)) then
857 0 : write(message,'(2a,i5,2a,2e15.6)') ch10,&
858 0 : & ' option dmft_solv=0 requires upaw=jpaw=0 for species',itypat,ch10,&
859 0 : & ' Value of upawu and jpawu are here',pawtab(itypat)%upawu,pawtab(itypat)%jpawu
860 0 : ABI_ERROR(message)
861 : end if
862 : end do ! itypat
863 : end if ! dmft_solv=0
864 :
865 : ! todo_ab: why upaw and jpawu are not zero (on bigmac) if lpawu==-1 ?
866 : ! if(paw_dmft%dmft_solv==0.and.&
867 : !& (maxval(abs(pawtab(:)%upawu))>tol5.or.maxval(abs(pawtab(:)%jpawu))>tol5))
868 : !then
869 : ! write(message, '(a,a,2f12.3)' )ch10,&
870 : !& ' option dmft_solv=0 requires
871 : !upaw=jpaw=0',maxval(abs(pawtab(:)%upawu)),maxval(abs(pawtab(:)%jpawu))
872 : ! ABI_WARNING(message)
873 : ! endif
874 :
875 43 : paw_dmft%dmftcheck = dtset%dmftcheck
876 :
877 43 : write(message,'(2a,i4)') ch10,'- ( number of procs used in dmft ) = ',nproc
878 129 : call wrtout([std_out,ab_out],message,'COLL')
879 43 : write(std_out_default,'(2a,i4)') ch10,' ( current proc is ) = ',myproc
880 : ! write(ab_out_default,'(2a,i3)') ch10,' ( current proc is ) =', myproc
881 43 : if (myproc == nproc-1) write(std_out_default,'(2a,i4)') ch10,' ( last proc ) = ',myproc
882 : ! write(ab_out_default,'(2a,i3)') ch10,' ( last proc ) =', myproc
883 :
884 : !#ifdef HAVE_MPI
885 : ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nb_procs,ie)
886 : ! write(6,*) "nprocs,nb_procs",nproc,nb_procs
887 : ! if(nb_procs/=nproc) then
888 : ! message = ' Number of procs used in DMFT is erroneously computed '
889 : ! ABI_ERROR(message)
890 : ! endif
891 : !#endif
892 :
893 : !=============================
894 : !== Associate pointers
895 : !=============================
896 :
897 43 : paw_dmft%nband => dtset%nband(:)
898 43 : paw_dmft%dmft_shiftself => dtset%dmft_shiftself(:)
899 43 : paw_dmft%dmft_nominal => dtset%dmft_nominal(:)
900 43 : paw_dmft%npwarr => npwarr(:)
901 :
902 : ! TODO: Make it work for usedmdft = -1 (interface with Wannier90 needs spinor
903 : ! generalization)
904 43 : if(nspinor==2.and.dtset%nspden==1.and.use_dmft==10) then
905 0 : message = ' nspinor==2 and nspden=1 and usedmft=10 is not implemented yet'
906 0 : ABI_ERROR(message)
907 : endif
908 :
909 1042 : paw_dmft%band_in(:) = .false.
910 1051880 : paw_dmft%occnd(:,:,:,:,:) = zero
911 43 : paw_dmft%use_dmft = use_dmft
912 :
913 : ! if (bandkss/=0) then
914 : ! paw_dmft%use_sc_dmft = 0
915 : ! else
916 : ! paw_dmft%use_sc_dmft = use_sc_dmft
917 : ! endif
918 : ! paw_dmft%dmft_read_occnd = dmft_read_occnd
919 : ! paw_dmft%idmftloop=0
920 : ! paw_dmft%mbandc = 0
921 :
922 43 : icb = 0
923 43 : mbandc = 0
924 1042 : do iband=1,mband
925 1042 : if (iband >= dmftbandi .and. iband <= dmftbandf) then
926 547 : paw_dmft%band_in(iband)=.true.
927 547 : mbandc = mbandc + 1
928 547 : paw_dmft%include_bands(mbandc) = iband
929 : else
930 452 : icb = icb + 1
931 452 : paw_dmft%exclude_bands(icb) = iband
932 : end if ! band>=bandi and band<=bandf
933 : end do ! iband
934 43 : paw_dmft%mbandc = mbandc
935 :
936 43 : bdtot_index = 1
937 114 : do isppol=1,nsppol
938 868 : do ikpt=1,nkpt
939 754 : nband_k = paw_dmft%nband(ikpt+(isppol-1)*nkpt)
940 16705 : do iband=1,nband_k
941 15880 : paw_dmft%occnd(1,iband,iband,ikpt,isppol) = occ(bdtot_index)
942 16634 : bdtot_index = bdtot_index + 1
943 : end do ! iband
944 : end do ! ikpt
945 : end do ! isppol
946 :
947 43 : if (paw_dmft%use_sc_dmft /= 0 .and. mpi_enreg%paral_kgb /= 0) then
948 4 : call init_sc_dmft_paralkgb(paw_dmft,mpi_enreg)
949 : end if
950 :
951 43 : if (mbandc /= dmftbandf-dmftbandi+1 .and. paw_dmft%use_dmft /= 10) then
952 0 : write(message,'(5a)') ' BUG init_sc_dmft',ch10,&
953 0 : & ' number of bands in dmft is not correctly computed ',ch10, &
954 0 : & ' Action : check the code'
955 0 : ABI_BUG(message)
956 : end if
957 :
958 43 : if (dmft_solv == 6 .or. dmft_solv == 7) then
959 :
960 0 : if (off_diag == 1) then
961 0 : write(message,'(3a)') "WARNING: You have activated the sampling of the off-diagonal elements ", &
962 0 : & "in TRIQS/CTHYB. Some features are not available, and you will not be ", &
963 0 : & "able to compute an energy."
964 0 : ABI_WARNING(message)
965 : end if
966 : #ifdef HAVE_TRIQS_COMPLEX
967 : if (off_diag == 0) then
968 : write(message,'(2a)') "WARNING: You have compiled with the complex version of TRIQS/CTHYB, yet you do not", &
969 : & " sample any off-diagonal element. This is a waste of computation time."
970 : ABI_WARNING(message)
971 : end if
972 : #else
973 : if (off_diag == 1) then
974 0 : write(message,'(3a)') "WARNING: You have compiled with the real version of TRIQS/CTHYB, yet you have", &
975 0 : & " activated the sampling of the off-diagonal elements. Thus their imaginary part will be", &
976 0 : & " neglected. You'll have to check that this is a valid approximation."
977 0 : ABI_WARNING(message)
978 : end if
979 : #endif
980 : end if ! dmft_solv=6 or 7
981 :
982 43 : write(message,'(7a)') ch10, &
983 43 : & ch10," ******************************************", &
984 43 : & ch10," DFT+DMFT Method is used", &
985 86 : & ch10," ******************************************"
986 129 : call wrtout([std_out,ab_out],message,'COLL')
987 :
988 43 : if (use_dmft /= 10) then
989 : if (dmft_solv == 0) then
990 1 : write(message,'(2a)') ch10,' DMFT check: no solver and U=J=0'
991 : else if (dmft_solv == 1) then
992 1 : write(message,'(2a)') ch10,' DMFT check: static solver'
993 : else if (dmft_solv == -1) then
994 1 : write(message,'(2a)') ch10,' DMFT check: static solver without renormalization of projectors: should recover DFT+U'
995 : else if (dmft_solv == 2) then
996 14 : write(message,'(2a)') ch10,' DMFT uses the Hubbard one solver'
997 : else if (dmft_solv == 4) then
998 0 : write(message,'(2a)') ch10,' DMFT uses the Hirsch Fye solver'
999 : else if (dmft_solv == 5) then
1000 26 : write(message,'(2a)') ch10,' DMFT uses the Continuous Time Quantum Monte Carlo solver of ABINIT'
1001 : else if (dmft_solv == 6) then
1002 0 : write(message,'(2a)') ch10,' DMFT uses the Continuous Time Quantum Monte Carlo solver of TRIQS &
1003 0 : &(with density density interactions)'
1004 : else if (dmft_solv == 7) then
1005 0 : write(message,'(2a)') ch10,' DMFT uses the Continuous Time Quantum Monte Carlo solver of TRIQS &
1006 0 : &(with rotationally invariant interactions)'
1007 : else if (dmft_solv == 8) then
1008 0 : write(message,'(2a)') ch10,' DMFT uses the Continuous Time Quantum Monte Carlo solver of ABINIT'
1009 : else if (dmft_solv == 9) then
1010 0 : write(message,'(2a)') ch10,' DMFT uses the python invocation of TRIQS, for which you need to &
1011 0 : & give your personal script'
1012 : else if (dmft_solv == 10) then
1013 0 : write(message,'(2a)') ch10,' DMFT uses the Complex Continuous Time Quantum Monte Carlo solver of ABINIT'
1014 : end if ! dmft_solv
1015 : else if(use_dmft == 10) then
1016 0 : write(message, '(a,a)') ch10,' DMFT uses the python invocation and orbitals constructed using Wannier90 '
1017 : endif
1018 129 : call wrtout([std_out,ab_out],message,'COLL')
1019 :
1020 : ! OG: What is all that? Something as moved? START
1021 43 : if (use_dmft /= 10) then
1022 : if (dmft_dc == 1) then
1023 26 : dc_string = "Magnetic FLL (Full Localized Limit)"
1024 : else if (dmft_dc == 2) then
1025 0 : dc_string = "Magnetic AMF (Around Mean Field)"
1026 : else if (dmft_dc == 5) then
1027 17 : dc_string = "Non-Magnetic FLL (Full Localized Limit)"
1028 : else if (dmft_dc == 6) then
1029 0 : dc_string = "Non-Magnetic AMF (Around Mean Field)"
1030 : else if (dmft_dc == 7) then
1031 0 : dc_string = "Non-Magnetic nominal"
1032 : else if (dmft_dc == 8) then
1033 0 : dc_string = "Non-Magnetic exact"
1034 : end if
1035 43 : dc_string = trim(dc_string) // " double counting"
1036 :
1037 43 : lda_string = "Magnetic DFT, with "
1038 43 : if (dtset%usepawu == 14) lda_string = "Non " // trim(adjustl(lda_string))
1039 43 : write(message,'(2(a,1x),a)') ch10,trim(adjustl(lda_string)),trim(adjustl(dc_string))
1040 129 : call wrtout([std_out,ab_out],message,'COLL')
1041 :
1042 43 : if (dtset%dmft_entropy == 0 .and. ((dmft_solv /= 6 .and. dmft_solv /= 7) &
1043 : & .or. (dtset%dmft_triqs_entropy == 0 .or. &
1044 : & dtset%dmft_triqs_compute_integral == 0))) then
1045 40 : write(message,'(a,1x,a)') ch10,"Entropy is not computed, only the internal energy is printed"
1046 120 : call wrtout([std_out,ab_out],message,'COLL')
1047 : end if
1048 :
1049 43 : if ((dmft_solv == 6 .or. dmft_solv == 7) .and. off_diag == 1) then
1050 : #ifndef HAVE_TRIQS_COMPLEX
1051 0 : write(message,'(a,1x,a)') ch10,"The imaginary part of the Green's function is neglected"
1052 0 : call wrtout([std_out,ab_out],message,'COLL')
1053 : #endif
1054 43 : else if (dmft_solv /= 6 .and. dmft_solv /= 7 .and. dmft_solv /= 10) then
1055 43 : write(message,'(a,1x,a)') ch10,"The imaginary part of the Green's function is neglected"
1056 129 : call wrtout([std_out,ab_out],message,'COLL')
1057 : end if
1058 43 : if (dmft_solv == 5 .or. ((dmft_solv == 6 .or. dmft_solv == 7) .and. off_diag == 0)) then
1059 26 : write(message,'(a,1x,a)') ch10,"The off-diagonal elements of the Green's function are neglected"
1060 78 : call wrtout([std_out,ab_out],message,'COLL')
1061 : end if
1062 :
1063 : !=============================
1064 : !== Define integers and reals
1065 : !=============================
1066 :
1067 43 : paw_dmft%nelectval = dble(dtset%nelect)
1068 :
1069 43 : if (dmft_solv /= 6 .and. dmft_solv /= 7) then
1070 43 : fac = merge(2,1,nsppol==1.and.nspinor==1)
1071 43 : paw_dmft%nelectval = dble(dtset%nelect-(dmftbandi-1)*nsppol*fac)
1072 : end if ! not use_all_bands
1073 :
1074 43 : paw_dmft%natpawu = dtset%natpawu
1075 43 : paw_dmft%natom = natom
1076 43 : paw_dmft%temp = dtset%tsmear!*unit_e
1077 43 : paw_dmft%dmft_hybri_limit = dtset%dmft_hybri_limit
1078 43 : paw_dmft%dmft_iter = dtset%dmft_iter
1079 43 : paw_dmft%dmft_entropy = dtset%dmft_entropy
1080 43 : paw_dmft%dmft_kspectralfunc = dtset%dmft_kspectralfunc
1081 43 : paw_dmft%dmft_magnfield = dtset%dmft_magnfield
1082 43 : paw_dmft%dmft_magnfield_b = dtset%dmft_magnfield_b
1083 43 : paw_dmft%dmft_dc = dmft_dc
1084 43 : paw_dmft%dmft_full_chipsi = dtset%dmft_full_chipsi
1085 43 : paw_dmft%dmft_wanorthnorm = dtset%dmft_wanorthnorm
1086 43 : paw_dmft%prtvol = dtset%prtvol
1087 43 : paw_dmft%prtdos = dtset%prtdos
1088 43 : paw_dmft%dmft_tolfreq = dtset%dmft_tolfreq
1089 43 : paw_dmft%dmft_lcpr = dtset%dmft_tollc
1090 43 : paw_dmft%dmft_charge_prec = dtset%dmft_charge_prec
1091 43 : paw_dmft%dmft_fermi_prec = dtset%dmft_charge_prec * ten
1092 43 : paw_dmft%dmft_fermi_step = dtset%dmft_fermi_step
1093 43 : paw_dmft%dmft_prt_maxent = dtset%dmft_prt_maxent
1094 43 : paw_dmft%dmft_prtself = dtset%dmft_prtself
1095 43 : paw_dmft%dmft_prtwan = dtset%dmft_prtwan
1096 43 : paw_dmft%dmft_wanrad = dtset%dmft_wanrad
1097 43 : paw_dmft%dmft_t2g = dtset%dmft_t2g
1098 43 : paw_dmft%dmft_x2my2d = dtset%dmft_x2my2d
1099 43 : paw_dmft%ixc = dtset%ixc
1100 :
1101 : ! for entropy (alternate external calculation)
1102 43 : paw_dmft%ientropy = 0
1103 43 : paw_dmft%u_for_s = 4.1_dp
1104 43 : paw_dmft%j_for_s = 0.5_dp
1105 :
1106 43 : paw_dmft%kptopt = dtset%kptopt
1107 :
1108 : !=======================
1109 : !== Choose solver
1110 : !=======================
1111 :
1112 43 : paw_dmft%dmft_solv = merge(2,dmft_solv,dmft_solv==-2)
1113 43 : paw_dmft%dmft_blockdiag = merge(1,0,dmft_solv==-2)
1114 :
1115 : ! 0: DFT, no solver
1116 : ! 1: DFT+U
1117 : ! -1: DFT+U but DFT values are not renormalized !
1118 : ! if((paw_dmft%dmft_solv==0.and.paw_dmft%prtvol>4).or.&
1119 : !& (paw_dmft%dmft_solv>=-1.and.paw_dmft%dmft_solv<=2)) then
1120 : ! call wrtout(std_out,message,'COLL')
1121 : ! call wrtout(ab_out,message,'COLL')
1122 : ! endif
1123 :
1124 : !=======================
1125 : !== Frequencies
1126 : !=======================
1127 :
1128 43 : paw_dmft%dmft_log_freq = merge(0,1,dmft_solv==6.or.dmft_solv==7.or.dmft_solv==9)
1129 43 : nwli = merge(dtset%dmft_triqs_n_iw,dtset%dmft_nwli,dmft_solv==6.or.dmft_solv==7)
1130 :
1131 43 : paw_dmft%dmft_nwli = nwli
1132 43 : paw_dmft%dmft_nwlo = merge(dtset%dmft_nwlo,nwli,paw_dmft%dmft_log_freq==1)
1133 43 : paw_dmft%dmft_nwr = 800
1134 :
1135 43 : paw_dmft%dmft_rslf = dtset%dmft_rslf
1136 43 : paw_dmft%dmft_mxsf = dtset%dmft_mxsf
1137 :
1138 : !=======================
1139 : !== CTQMC
1140 : !=======================
1141 :
1142 43 : paw_dmft%dmftqmc_l = merge(dtset%dmft_triqs_n_tau,dtset%dmftqmc_l,dmft_solv==6.or.dmft_solv==7)
1143 43 : paw_dmft%dmftqmc_n = dtset%dmftqmc_n
1144 43 : paw_dmft%dmftqmc_seed = dtset%dmftqmc_seed
1145 43 : paw_dmft%dmftqmc_therm = merge(dtset%dmft_triqs_n_warmup_cycles_init,dtset%dmftqmc_therm,dmft_solv==6.or.dmft_solv==7)
1146 :
1147 43 : paw_dmft%dmftctqmc_basis = merge(dtset%dmft_triqs_basis,dtset%dmftctqmc_basis,dmft_solv==6.or.dmft_solv==7)
1148 43 : paw_dmft%dmftctqmc_check = dtset%dmftctqmc_check
1149 43 : paw_dmft%dmftctqmc_correl = dtset%dmftctqmc_correl
1150 43 : paw_dmft%dmftctqmc_gmove = dtset%dmftctqmc_gmove
1151 43 : paw_dmft%dmftctqmc_grnns = dtset%dmftctqmc_grnns
1152 43 : paw_dmft%dmftctqmc_meas = merge(dtset%dmft_triqs_length_cycle,dtset%dmftctqmc_meas,dmft_solv==6.or.dmft_solv==7)
1153 43 : paw_dmft%dmftctqmc_mrka = dtset%dmftctqmc_mrka
1154 43 : paw_dmft%dmftctqmc_mov = dtset%dmftctqmc_mov
1155 43 : paw_dmft%dmftctqmc_order = dtset%dmftctqmc_order
1156 43 : paw_dmft%dmftctqmc_chains = dtset%dmftctqmc_chains
1157 43 : paw_dmft%dmftctqmc_localprop = dtset%dmftctqmc_localprop
1158 :
1159 43 : if (dmft_solv == 5 .or. dmft_solv >= 8) then
1160 26 : write(message,'(2a,i6)') ch10,&
1161 52 : & '=> Seed for CT-QMC inside DMFT is dmftqmc_seed=',paw_dmft%dmftqmc_seed
1162 26 : call wrtout(std_out,message,'COLL')
1163 : end if
1164 :
1165 : !=======================
1166 : !== TRIQS CTQMC
1167 : !=======================
1168 :
1169 43 : paw_dmft%dmft_triqs_nleg = dtset%dmft_triqs_n_l
1170 43 : paw_dmft%dmft_triqs_therm_restart = dtset%dmft_triqs_n_warmup_cycles_restart
1171 43 : paw_dmft%dmft_triqs_det_init_size = dtset%dmft_triqs_det_init_size
1172 43 : paw_dmft%dmft_triqs_det_n_operations_before_check = dtset%dmft_triqs_det_n_operations_before_check
1173 43 : paw_dmft%dmft_triqs_move_shift = (dtset%dmft_triqs_move_shift == 1)
1174 43 : paw_dmft%dmft_triqs_move_double = (dtset%dmft_triqs_move_double == 1)
1175 43 : paw_dmft%dmft_triqs_loc_n_min = dtset%dmft_triqs_loc_n_min
1176 43 : paw_dmft%dmft_triqs_loc_n_max = dtset%dmft_triqs_loc_n_max
1177 43 : paw_dmft%dmft_triqs_seed_a = dtset%dmft_triqs_random_seed_a
1178 43 : paw_dmft%dmft_triqs_seed_b = dtset%dmft_triqs_random_seed_b
1179 43 : paw_dmft%dmft_triqs_measure_density_matrix = (dtset%dmft_triqs_measure_density_matrix == 1)
1180 43 : paw_dmft%dmft_triqs_time_invariance = (dtset%dmft_triqs_time_invariance == 1)
1181 43 : paw_dmft%dmft_triqs_use_norm_as_weight = (dtset%dmft_triqs_use_norm_as_weight == 1)
1182 43 : paw_dmft%dmft_triqs_leg_measure = (dtset%dmft_triqs_measure_g_l == 1)
1183 43 : paw_dmft%dmft_triqs_off_diag = (off_diag == 1)
1184 43 : paw_dmft%dmft_triqs_imag_threshold = dtset%dmft_triqs_imag_threshold
1185 43 : paw_dmft%dmft_triqs_det_precision_warning = dtset%dmft_triqs_det_precision_warning
1186 43 : paw_dmft%dmft_triqs_det_precision_error = dtset%dmft_triqs_det_precision_error
1187 43 : paw_dmft%dmft_triqs_det_singular_threshold = dtset%dmft_triqs_det_singular_threshold
1188 43 : paw_dmft%dmft_triqs_epsilon = dtset%dmft_triqs_dlr_epsilon
1189 43 : paw_dmft%dmft_triqs_lambda = dtset%dmft_triqs_dlr_wmax / dtset%tsmear
1190 43 : paw_dmft%dmft_triqs_entropy = dtset%dmft_triqs_entropy
1191 43 : paw_dmft%dmft_triqs_compute_integral = dtset%dmft_triqs_compute_integral
1192 43 : paw_dmft%dmft_triqs_gaussorder = dtset%dmft_triqs_gaussorder
1193 43 : paw_dmft%dmft_triqs_nsubdivisions = dtset%dmft_triqs_nsubdivisions
1194 43 : paw_dmft%dmft_triqs_tol_block = dtset%dmft_triqs_tol_block
1195 43 : paw_dmft%dmft_triqs_read_ctqmcdata = dtset%dmft_triqs_read_ctqmcdata
1196 43 : paw_dmft%dmft_triqs_pauli_prob = dtset%dmft_triqs_pauli_prob
1197 43 : paw_dmft%dmft_triqs_n_cycles = dtset%dmft_triqs_n_cycles
1198 43 : paw_dmft%dmft_triqs_prt_entropy = (dtset%dmft_triqs_prt_entropy == 1)
1199 43 : paw_dmft%dmft_triqs_shift_mu = dtset%dmft_triqs_shift_mu
1200 43 : paw_dmft%dmft_triqs_chiloc = dtset%dmft_triqs_chiloc
1201 43 : paw_dmft%dmft_triqs_chiloc_ins = dtset%dmft_triqs_chiloc_ins
1202 :
1203 : !==============================
1204 : !== Variables for DMFT itself
1205 : !==============================
1206 :
1207 43 : paw_dmft%wtk => dtset%wtk(:)
1208 : ! In the case where we sample the full BZ, don't overwrite the wtk with 1/nkpt when we use TRIQS
1209 43 : if (dtset%iscf < 0 .and. (dtset%kptopt < 0 .or. &
1210 87 : & (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7))) paw_dmft%wtk(:) = one / dble(nkpt)
1211 483 : sumwtk = sum(paw_dmft%wtk(1:nkpt))
1212 43 : if (abs(sumwtk-one) > tol11) then
1213 0 : write(message,'(a,f15.11)') ' sum of k-point is incorrect',sumwtk
1214 0 : ABI_BUG(message)
1215 : end if
1216 :
1217 43 : t2g = (paw_dmft%dmft_t2g == 1)
1218 43 : x2my2d = (paw_dmft%dmft_x2my2d == 1)
1219 :
1220 43 : paw_dmft%typat => dtset%typat(:)
1221 :
1222 129 : ABI_MALLOC(paw_dmft%lpawu,(natom))
1223 :
1224 43 : maxlpawu = 0
1225 139 : do iatom=1,natom
1226 96 : lpawu = pawtab(paw_dmft%typat(iatom))%lpawu
1227 96 : if (t2g .and. lpawu /= -1) lpawu = 1
1228 96 : if (x2my2d .and. lpawu /= -1) lpawu = 0
1229 96 : if (lpawu > maxlpawu) maxlpawu = lpawu
1230 139 : paw_dmft%lpawu(iatom) = lpawu
1231 : end do ! iatom
1232 43 : paw_dmft%maxlpawu = maxlpawu
1233 43 : ndim = 2*maxlpawu + 1
1234 :
1235 43 : paw_dmft%maxnproju = 0
1236 129 : ABI_MALLOC(paw_dmft%siz_proj,(ntypat))
1237 129 : ABI_MALLOC(lcycle,(maxlpawu+1))
1238 :
1239 173 : lcycle(:) = .true.
1240 :
1241 111 : do itypat=1,ntypat
1242 68 : lpawu = pawtab(itypat)%lpawu
1243 68 : if (lpawu == -1) cycle
1244 43 : if (t2g) lpawu = 1
1245 43 : if (x2my2d) lpawu = 0
1246 43 : paw_dmft%siz_proj(itypat) = size(pawtab(itypat)%proj(:))
1247 43 : lcycle(lpawu+1) = .false.
1248 86 : if (pawtab(itypat)%nproju > paw_dmft%maxnproju) paw_dmft%maxnproju = pawtab(itypat)%nproju
1249 : end do ! itypat
1250 :
1251 215 : ABI_MALLOC(paw_dmft%slm2ylm,(ndim,ndim,maxlpawu+1))
1252 215 : ABI_MALLOC(paw_dmft%jmj2ylm,(2*ndim,2*ndim,maxlpawu+1))
1253 258 : ABI_MALLOC(paw_dmft%zarot,(ndim,ndim,nsym,maxlpawu+1))
1254 :
1255 4653 : paw_dmft%slm2ylm(:,:,:) = czero
1256 173 : do lpawu=0,maxlpawu
1257 130 : if (lcycle(lpawu+1)) cycle
1258 43 : ndim = 2*lpawu + 1
1259 303 : do im=1,ndim
1260 217 : mm = im - lpawu - 1 ; jm = - mm + lpawu + 1
1261 217 : onem = (-1)**mm
1262 217 : if (mm > 0) then
1263 87 : paw_dmft%slm2ylm(im,im,lpawu+1) = cmplx(onem/sqrt2,zero,kind=dp)
1264 87 : paw_dmft%slm2ylm(jm,im,lpawu+1) = cmplx(one/sqrt2,zero,kind=dp)
1265 : end if
1266 217 : if (mm == 0) paw_dmft%slm2ylm(im,im,lpawu+1) = cone
1267 347 : if (mm < 0) then
1268 87 : paw_dmft%slm2ylm(im,im,lpawu+1) = cmplx(zero,one/sqrt2,kind=dp)
1269 87 : paw_dmft%slm2ylm(jm,im,lpawu+1) = -cmplx(zero,onem/sqrt2,kind=dp)
1270 : end if
1271 : end do ! im
1272 : end do ! lpawu
1273 :
1274 16721 : paw_dmft%jmj2ylm(:,:,:) = czero
1275 130 : do ll=1,maxlpawu
1276 87 : if (lcycle(ll+1)) cycle
1277 129 : ABI_MALLOC(ind_msml,(2,-ll:ll))
1278 129 : jc1 = 0
1279 129 : do ms1=1,2
1280 563 : do ml1=-ll,ll
1281 434 : jc1 = jc1 + 1
1282 520 : ind_msml(ms1,ml1) = jc1
1283 : end do ! ml1
1284 : end do ! ms1
1285 43 : invsqrt2lp1 = one / sqrt(dble(2*ll+1))
1286 43 : jc1 = 0
1287 129 : do jj=ll,ll+1
1288 86 : xj = dble(jj) - half ! xj is in {ll-0.5,ll+0.5}
1289 563 : do jm=-jj,jj-1
1290 434 : xmj = dble(jm) + half ! xmj is in {-xj,xj}
1291 434 : jc1 = jc1 + 1 ! Global index for JMJ
1292 434 : if (nint(xj+half) == ll+1) then ! if xj=ll+0.5
1293 260 : if (nint(xmj+half) == ll+1) then
1294 43 : paw_dmft%jmj2ylm(ind_msml(1,ll),jc1,ll+1) = cone ! J=L+0.5 and m_J=L+0.5
1295 217 : else if (nint(xmj-half) == -ll-1) then
1296 43 : paw_dmft%jmj2ylm(ind_msml(2,-ll),jc1,ll+1) = cone ! J=L+0.5 and m_J=-L-0.5
1297 : else
1298 174 : paw_dmft%jmj2ylm(ind_msml(1,nint(xmj-half)),jc1,ll+1) = cmplx(invsqrt2lp1*(sqrt(dble(ll)+xmj+half)),zero,kind=dp)
1299 174 : paw_dmft%jmj2ylm(ind_msml(2,nint(xmj+half)),jc1,ll+1) = cmplx(invsqrt2lp1*(sqrt(dble(ll)-xmj+half)),zero,kind=dp)
1300 : end if
1301 : end if
1302 520 : if (nint(xj+half) == ll) then ! if xj=ll-0.5
1303 174 : paw_dmft%jmj2ylm(ind_msml(2,nint(xmj+half)),jc1,ll+1) = cmplx(invsqrt2lp1*(sqrt(dble(ll)+xmj+half)),zero,kind=dp)
1304 174 : paw_dmft%jmj2ylm(ind_msml(1,nint(xmj-half)),jc1,ll+1) = cmplx(-invsqrt2lp1*(sqrt(dble(ll)-xmj+half)),zero,kind=dp)
1305 : end if
1306 : end do ! jm
1307 : end do ! jj
1308 130 : ABI_FREE(ind_msml)
1309 : end do ! ll
1310 :
1311 173 : do lpawu=0,maxlpawu
1312 130 : if (lcycle(lpawu+1)) cycle
1313 43 : ndim = 2*lpawu + 1
1314 86 : if (t2g) then
1315 28 : do im1=1,ndim
1316 91 : do im=1,ndim
1317 3108 : paw_dmft%zarot(im,im1,:,lpawu+1) = cmplx(pawang%zarot(mt2g(im),mt2g(im1),3,1:nsym),zero,kind=dp)
1318 : end do ! im
1319 : end do ! im1
1320 36 : else if (x2my2d) then
1321 0 : paw_dmft%zarot(1,1,:,lpawu+1) = cmplx(pawang%zarot(5,5,3,1:nsym),zero,kind=dp)
1322 : else
1323 56727 : paw_dmft%zarot(1:ndim,1:ndim,:,lpawu+1) = cmplx(pawang%zarot(1:ndim,1:ndim,lpawu+1,1:nsym),zero,kind=dp)
1324 : end if
1325 : end do ! lpawu
1326 :
1327 : !=======================
1328 : ! Imaginary frequencies
1329 : !=======================
1330 : ! Set up log frequencies
1331 43 : if (dtset%ucrpa == 0 .and. paw_dmft%dmft_nwlo > 0) then
1332 42 : call construct_nwlo_dmft(paw_dmft)
1333 : end if
1334 :
1335 43 : if (paw_dmft%dmftcheck == 1 .and. dmft_solv < 4) paw_dmft%dmftqmc_l = 64
1336 :
1337 : !==============
1338 : ! Radial grid
1339 : !==============
1340 :
1341 43 : if (paw_dmft%dmft_prtwan == 1) then
1342 :
1343 : ! Initialize rmax (maximum radius for print of Wannier functions)
1344 0 : ABI_MALLOC(rmax,(maxlpawu+1))
1345 0 : rmax(:) = paw_dmft%dmft_wanrad
1346 :
1347 0 : if (paw_dmft%dmft_wanrad < 0) then ! default
1348 :
1349 : ! Set rmax to half the distance from the current atom
1350 : ! to the nearest atom with the same lpawu
1351 :
1352 : ! First take half * min(|Ri|)
1353 0 : rmax(:) = zero
1354 : verif = .true.
1355 0 : do i=1,3
1356 0 : norm = norm2(rprimd(1:3,i)) * half
1357 0 : if (verif .or. norm < rmax(1)) then
1358 0 : rmax(:) = norm
1359 : verif = .false.
1360 : end if
1361 : end do ! i
1362 :
1363 : ! Now look at the atoms in the same unit cell
1364 0 : do iatom=1,natom
1365 0 : lpawu = paw_dmft%lpawu(iatom)
1366 0 : if (lpawu == -1) cycle
1367 0 : do iatom1=1,natom
1368 0 : lpawu1 = paw_dmft%lpawu(iatom1)
1369 0 : if (lpawu1 == -1) cycle
1370 0 : if (lpawu /= lpawu1 .or. iatom == iatom1) cycle
1371 : norm = zero
1372 0 : do i=1,3
1373 0 : norm = norm + dot_product(xred(1:3,iatom)-xred(1:3,iatom1),rprimd(i,1:3))**2
1374 : end do
1375 0 : norm = sqrt(norm) * half
1376 0 : if (norm < rmax(lpawu+1)) rmax(lpawu+1) = norm
1377 : end do ! iatom1
1378 : end do ! iatom
1379 :
1380 : end if ! dmft_wanrad < 0
1381 :
1382 : end if ! prtwan=1
1383 :
1384 : ! Now build radial grid by extending the PAW mesh up to max(rmax,size(proj))
1385 : ! The mesh inside the PAW sphere is still exactly the same.
1386 : !use_full_chipsi = (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7)
1387 43 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) then
1388 : use_full_chipsi = .true.
1389 : else
1390 43 : use_full_chipsi = (paw_dmft%dmft_full_chipsi /= 0)
1391 : endif
1392 43 : paw_dmft%int_meshsz => pawrad(:)%int_meshsz
1393 :
1394 43 : if (use_full_chipsi) then
1395 0 : ABI_MALLOC(paw_dmft%phimtphi,(maxval(pawrad(1:ntypat)%int_meshsz),paw_dmft%maxnproju,ntypat))
1396 0 : ABI_MALLOC(paw_dmft%phimtphi_int,(paw_dmft%maxnproju,ntypat))
1397 : else
1398 172 : ABI_MALLOC(paw_dmft%phi_int,(paw_dmft%maxnproju,ntypat))
1399 : end if
1400 :
1401 197 : ABI_MALLOC(paw_dmft%radgrid,(ntypat))
1402 :
1403 43 : paw_dmft%maxmeshsize = 0
1404 111 : do itypat=1,ntypat
1405 68 : lpawu = pawtab(itypat)%lpawu
1406 68 : if (lpawu == -1) cycle
1407 43 : if (t2g) lpawu = 1
1408 43 : if (x2my2d) lpawu = 0
1409 43 : mesh_type = pawrad(itypat)%mesh_type
1410 43 : lstep = pawrad(itypat)%lstep
1411 43 : rstep = pawrad(itypat)%rstep
1412 43 : siz_proj = paw_dmft%siz_proj(itypat)
1413 43 : mesh_size = siz_proj
1414 43 : if (paw_dmft%dmft_prtwan == 1) then
1415 0 : if (mesh_type == 1) mesh_size = int(rmax(lpawu+1)/rstep) + 1
1416 0 : if (mesh_type == 2) mesh_size = int(log(rmax(lpawu+1)/rstep+1)/lstep) + 1
1417 0 : if (mesh_type == 3) mesh_size = int(log(rmax(lpawu+1)/rstep)/lstep) + 2
1418 0 : if (mesh_size < siz_proj) then
1419 0 : message = "Please set wanrad to a value greater than the radius of your DMFT orbital"
1420 0 : ABI_ERROR(message)
1421 : end if
1422 : end if ! prtwan=1
1423 : ! mesh_type > 3 cannot be extended outside the PAW sphere while keeping the
1424 : ! mesh inside the sphere unchanged.
1425 43 : if ((mesh_size /= pawrad(itypat)%mesh_size) .and. mesh_type > 3) then
1426 0 : message = "mesh_type > 3 is only compatible with dmft_orbital=1 and dmft_prtwan=0"
1427 0 : ABI_ERROR(message)
1428 : end if
1429 43 : if (mesh_size > pawrad(itypat)%int_meshsz .and. (.not. use_full_chipsi)) then
1430 0 : message = "You need to activate dmft_full_chipsi if you use an orbital that extends outside the PAW sphere"
1431 0 : ABI_ERROR(message)
1432 : end if
1433 43 : call pawrad_init(paw_dmft%radgrid(itypat),mesh_size,mesh_type,rstep,lstep)
1434 43 : if (mesh_size > paw_dmft%maxmeshsize) paw_dmft%maxmeshsize = mesh_size
1435 43 : siz_paw = min(mesh_size,pawrad(itypat)%int_meshsz)
1436 43 : siz_proj = min(siz_proj,pawrad(itypat)%int_meshsz)
1437 43 : rint = paw_dmft%radgrid(itypat)%rad(siz_proj)
1438 43 : nproju = pawtab(itypat)%nproju
1439 161 : do iproj=1,nproju
1440 75 : indproj = pawtab(itypat)%lnproju(iproj)
1441 143 : if (use_full_chipsi) then
1442 : ! Precompute <Chi|Phi-Phi_tilde>
1443 : paw_dmft%phimtphi(1:siz_paw,iproj,itypat) = pawtab(itypat)%phi(1:siz_paw,indproj) - &
1444 0 : & pawtab(itypat)%tphi(1:siz_paw,indproj)
1445 : call simp_gen(paw_dmft%phimtphi_int(iproj,itypat),pawtab(itypat)%proj(1:siz_proj)* &
1446 0 : & paw_dmft%phimtphi(1:siz_proj,iproj,itypat),paw_dmft%radgrid(itypat),r_for_intg=rint)
1447 : else
1448 : ! Precompute <Chi|Phi>
1449 : call simp_gen(paw_dmft%phi_int(iproj,itypat),pawtab(itypat)%proj(1:siz_proj)* &
1450 56739 : & pawtab(itypat)%phi(1:siz_proj,indproj),paw_dmft%radgrid(itypat),r_for_intg=rint)
1451 : end if ! use_full_chipsi
1452 : end do ! iproj
1453 : end do ! itypat
1454 :
1455 43 : if (paw_dmft%dmft_prtwan /= 1 .and. use_full_chipsi) then
1456 0 : ABI_FREE(paw_dmft%phimtphi)
1457 : end if
1458 43 : ABI_SFREE(rmax)
1459 :
1460 : !==============
1461 : ! Plane waves
1462 : !==============
1463 :
1464 43 : if (use_full_chipsi) then
1465 :
1466 : ! Compute ylm(k+G),exp(j*(k+G).R(iat)) and Bessel functions
1467 0 : mpw = dtset%mpw
1468 0 : ABI_MALLOC(paw_dmft%ylm,(mpw,2*maxlpawu+1,maxlpawu+1,mkmem))
1469 0 : ABI_MALLOC(paw_dmft%dpro,(mpw,natom,mkmem))
1470 0 : ABI_MALLOC(paw_dmft%bessel,(mpw,paw_dmft%maxmeshsize,ntypat,mkmem))
1471 0 : ABI_MALLOC(paw_dmft%bessel_int,(mpw,ntypat,mkmem))
1472 0 : ABI_MALLOC(typcycle,(ntypat))
1473 0 : ABI_MALLOC(kpg,(3,mpw))
1474 0 : ABI_MALLOC(kpg_norm,(mpw))
1475 :
1476 0 : ik = 0 ! kpt index on current CPU
1477 0 : ikg = 0
1478 0 : fac_bessel = four_pi / sqrt(ucvol)
1479 :
1480 0 : do ikpt=1,nkpt
1481 :
1482 0 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,1,-1,mpi_enreg%me_kpt)) cycle
1483 :
1484 0 : ik = ik + 1
1485 0 : npw = npwarr(ikpt)
1486 :
1487 0 : do ig=1,npw
1488 :
1489 0 : kpg(:,ig) = dtset%kptns(1:3,ikpt) + dble(kg(1:3,ikg+ig))
1490 :
1491 : norm = zero
1492 0 : do i=1,3
1493 0 : norm = norm + dot_product(kpg(:,ig),gprimd(i,1:3))**2
1494 : end do ! i
1495 :
1496 0 : kpg_norm(ig) = sqrt(norm)
1497 :
1498 : end do ! ig
1499 :
1500 0 : lcycle(:) = .false.
1501 0 : typcycle(:) = .false.
1502 :
1503 0 : do iatom=1,natom
1504 :
1505 0 : lpawu = paw_dmft%lpawu(iatom)
1506 0 : if (lpawu == -1) cycle
1507 0 : itypat = paw_dmft%typat(iatom)
1508 0 : ndim = 2*lpawu + 1
1509 :
1510 0 : if (.not. lcycle(lpawu+1)) then ! if this l has not been visited
1511 0 : if (t2g) then
1512 0 : do im=1,ndim
1513 0 : paw_dmft%ylm(1:npw,im,lpawu+1,ik) = ylm(ikg+1:ikg+npw,4+mt2g(im))
1514 : end do ! im
1515 0 : else if (x2my2d) then
1516 0 : paw_dmft%ylm(1:npw,1,lpawu+1,ik) = ylm(ikg+1:ikg+npw,9)
1517 : else
1518 0 : paw_dmft%ylm(1:npw,1:ndim,lpawu+1,ik) = ylm(ikg+1:ikg+npw,lpawu**2+1:lpawu**2+ndim)
1519 : end if
1520 0 : lcycle(lpawu+1) = .true.
1521 : end if ! not lcycle
1522 :
1523 0 : do ig=1,npw
1524 0 : paw_dmft%dpro(ig,iatom,ik) = exp(j_dpc*two_pi*dot_product(kpg(:,ig),xred(1:3,iatom)))
1525 : end do ! ig
1526 :
1527 0 : if (.not. typcycle(itypat)) then ! if this type has not been visited
1528 0 : lpawu1 = lpawu ! physical l
1529 0 : if (t2g .or. x2my2d) lpawu1 = 2
1530 0 : j_l = j_dpc**lpawu1
1531 0 : siz_proj = paw_dmft%siz_proj(itypat)
1532 0 : rint = paw_dmft%radgrid(itypat)%rad(siz_proj)
1533 0 : siz_wan = paw_dmft%radgrid(itypat)%mesh_size
1534 0 : do ir=1,siz_wan
1535 0 : rad = paw_dmft%radgrid(itypat)%rad(ir)
1536 0 : do ig=1,npw
1537 0 : call paw_jbessel_4spline(bes,besp,lpawu1,0,two_pi*kpg_norm(ig)*rad,tol3)
1538 : ! Multiply by r since we want to compute Psi(r) * r, for radial integration
1539 0 : paw_dmft%bessel(ig,ir,itypat,ik) = fac_bessel * bes * rad
1540 : end do ! ig
1541 : end do ! ir
1542 0 : do ig=1,npw
1543 : call simp_gen(bes,pawtab(itypat)%proj(1:siz_proj)*dble(paw_dmft%bessel(ig,1:siz_proj,itypat,ik)), &
1544 0 : & paw_dmft%radgrid(itypat),r_for_intg=rint)
1545 0 : paw_dmft%bessel_int(ig,itypat,ik) = bes * j_l ! CAREFUL: we multiply by j^l AFTER simp_gen since simp_gen doesn_t handle complex
1546 : end do ! ig
1547 0 : paw_dmft%bessel(1:npw,1:siz_wan,itypat,ik) = paw_dmft%bessel(1:npw,1:siz_wan,itypat,ik) * j_l
1548 0 : typcycle(itypat) = .true.
1549 : end if ! not typcycle
1550 :
1551 : end do ! iatom
1552 :
1553 0 : ikg = ikg + npw
1554 :
1555 : end do ! ikpt
1556 :
1557 0 : ABI_FREE(kpg)
1558 0 : ABI_FREE(kpg_norm)
1559 0 : ABI_FREE(typcycle)
1560 :
1561 0 : if (paw_dmft%dmft_prtwan /= 1) then
1562 0 : ABI_FREE(paw_dmft%bessel)
1563 : end if
1564 :
1565 : end if ! use_full_chipsi
1566 :
1567 43 : ABI_FREE(lcycle)
1568 :
1569 43 : call init_paral_dmft(paw_dmft,paw_dmft%distrib,paw_dmft%dmft_nwlo)
1570 :
1571 : ! OG: What is all that? Something as moved? START
1572 : endif
1573 :
1574 16937 : end subroutine init_sc_dmft
1575 : !!***
1576 :
1577 : !!****f* m_paw_dmft/init_dmft
1578 : !! NAME
1579 : !! init_dmft
1580 : !!
1581 : !! FUNCTION
1582 : !! Allocate variables and setup DFT hamiltonian and related data
1583 : !! (init_sc_dmft has to been called before)
1584 : !!
1585 : !! INPUTS
1586 : !! cryst_struc <type(crystal_t)> = crystal structure data
1587 : !! dmatpawu = fixed occupation matrix of correlated orbitals
1588 : !! dtset <type(dataset_type)> = all input variables for this dataset
1589 : !! fermie_dft = DFT Fermi level
1590 : !! filctqmcdatain = input file name for CTQMC_DATA file
1591 : !! filselfin = input file name for self file
1592 : !! fnamei = input file name
1593 : !! fnametmp_app = header for the output filename
1594 : !! ireadctqmcdata = flag to read CTQMC_DATA input file at first iteration
1595 : !! ireadself = flag to read self input file at first iteration
1596 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
1597 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1598 : !!
1599 : !! SOURCE
1600 : !!
1601 : !! NOTE
1602 : !! The part of the code which deals
1603 : !! with the use of logarithmic frequencies
1604 : !! is a modification of the GNU GPL
1605 : !! code available on http://dmft.rutgers.edu/ and
1606 : !! described in the RMP paper written by
1607 : !! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti.
1608 :
1609 97 : subroutine init_dmft(cryst_struc,dmatpawu,dtset,fermie_dft,filctqmcdatain,filselfin,fnamei,fnametmp_app,ireadctqmcdata,ireadself,paw_dmft,pawtab)
1610 :
1611 : !Arguments ------------------------------------
1612 : real(dp), intent(in) :: fermie_dft
1613 : type(dataset_type), intent(in) :: dtset
1614 : type(paw_dmft_type), intent(inout) :: paw_dmft
1615 : type(crystal_t), target, intent(in) :: cryst_struc
1616 : character(len=fnlen), intent(in) :: filctqmcdatain,filselfin,fnamei,fnametmp_app
1617 : integer, intent(in) :: ireadctqmcdata,ireadself
1618 : real(dp), target, intent(in) :: dmatpawu(:,:,:,:)
1619 : type(pawtab_type), intent(in) :: pawtab(dtset%ntypat)
1620 : !Local variables ------------------------------------
1621 : integer :: iatom,ierr,ifreq,ioerr,ir,irot,isym
1622 : integer :: itypat,lpawu,meshsz,nflavor,ngrid,nsym,unt
1623 : real(dp) :: int1,step
1624 : logical :: lexist
1625 : character(len=4) :: tag_at
1626 : character(len=500) :: message
1627 : character(len=fnlen) :: tmpfil
1628 : integer, parameter :: grid_unt = 2000
1629 : ! *********************************************************************
1630 :
1631 97 : if (dtset%ucrpa == 0) then
1632 97 : write(message,'(7a)') ch10,' ====================================', &
1633 97 : & ch10,' ===== Start of DMFT calculation', &
1634 194 : & ch10,' ====================================',ch10
1635 0 : else if (dtset%ucrpa > 0) then
1636 0 : write(message,'(6a)') ch10,' ============================================================', &
1637 0 : & ch10,' ===== Initialize construction of Wannier in DMFT routines',&
1638 0 : & ch10,' ============================================================'
1639 : end if ! ucrpa
1640 97 : call wrtout(std_out,message,'COLL')
1641 :
1642 97 : if (paw_dmft%dmft_dc == 8) then
1643 0 : ABI_MALLOC(paw_dmft%edc,(paw_dmft%natom))
1644 0 : ABI_MALLOC(paw_dmft%edcdc,(paw_dmft%natom))
1645 : end if
1646 :
1647 97 : nsym = paw_dmft%nsym
1648 :
1649 : !=======================
1650 : !== Check sym
1651 : !=======================
1652 :
1653 4417 : do isym=1,nsym
1654 4417 : if (dtset%symafm(isym) < 0) then
1655 0 : message = 'symafm negative is not implemented in DMFT'
1656 0 : ABI_ERROR(message)
1657 : end if
1658 : end do ! isym
1659 :
1660 97 : paw_dmft%nsym = cryst_struc%nsym ! very important to update it here
1661 97 : nsym = paw_dmft%nsym
1662 :
1663 : ! TODO: this really should be done in init_sc_dmft
1664 97 : paw_dmft%indsym => cryst_struc%indsym(4,1:nsym,1:paw_dmft%natom)
1665 97 : if (paw_dmft%nspinor == 2) then
1666 3 : ABI_MALLOC(paw_dmft%symrec_cart,(3,3,nsym))
1667 25 : do irot=1,nsym
1668 : call symredcart(cryst_struc%gprimd(:,:),cryst_struc%rprimd(:,:),&
1669 25 : & paw_dmft%symrec_cart(:,:,irot),cryst_struc%symrec(:,:,irot))
1670 : end do ! irot
1671 : end if ! nspinor=2
1672 :
1673 97 : paw_dmft%filapp = fnametmp_app
1674 97 : paw_dmft%filnamei = fnamei
1675 97 : paw_dmft%filselfin = filselfin
1676 97 : paw_dmft%filctqmcdatain = filctqmcdatain
1677 97 : paw_dmft%ireadctqmcdata = ireadctqmcdata
1678 97 : paw_dmft%ireadself = ireadself
1679 :
1680 : ! Write orbital on file
1681 97 : if (paw_dmft%myproc == 0) then
1682 130 : do itypat=1,paw_dmft%ntypat
1683 87 : lpawu = pawtab(itypat)%lpawu
1684 87 : if (lpawu == -1) cycle
1685 43 : meshsz = paw_dmft%siz_proj(itypat)
1686 :
1687 : call simp_gen(int1,pawtab(itypat)%proj(1:meshsz)**2,paw_dmft%radgrid(itypat), &
1688 38019 : & r_for_intg=paw_dmft%radgrid(itypat)%rad(meshsz))
1689 43 : int1 = sqrt(int1)
1690 :
1691 43 : call int2char4(itypat,tag_at)
1692 43 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1693 43 : if (open_file(trim(paw_dmft%filapp)//"_DMFTORBITAL_itypat"//tag_at//".dat",message,newunit=unt) /= 0) ABI_ERROR(message)
1694 :
1695 43 : write(unt,'(4a)') "# Correlated normalized radial orbital for DMFT. This", &
1696 43 : & " is not projected on any energy window (you need to use dmft_prtwan for that).",ch10, &
1697 86 : & "# Radius (Bohr) u_l(r) = R_l * r"
1698 :
1699 38019 : do ir=1,meshsz
1700 38019 : write(unt,*) paw_dmft%radgrid(itypat)%rad(ir),pawtab(itypat)%proj(ir)/int1
1701 : end do ! ir
1702 :
1703 130 : close(unt)
1704 : end do ! itypat
1705 : end if ! myproc=0
1706 :
1707 : !==================
1708 : ! Real frequencies
1709 : !==================
1710 :
1711 97 : if (dtset%iscf < 0 .and. paw_dmft%dmft_solv >= 5 .and. paw_dmft%dmft_solv <= 8) then
1712 4 : tmpfil = trim(paw_dmft%filapp)//'_spectralfunction_realgrid'
1713 4 : inquire(file=trim(tmpfil),exist=lexist)!,recl=nrecl)
1714 4 : if (.not. lexist) then
1715 0 : write(message,'(4x,a,i5,3a)') "File number",grid_unt," called ",trim(tmpfil)," does not exist"
1716 0 : call wrtout(std_out,message,'COLL')
1717 0 : message = "Cannot continue: the missing file coming from Maxent code is needed"
1718 0 : ABI_ERROR(message)
1719 : end if ! not lexist
1720 :
1721 4 : if (paw_dmft%myproc == 0) then
1722 : #ifdef FC_NAG
1723 : open(unit=grid_unt,file=trim(tmpfil),status='unknown',form='formatted',recl=ABI_RECL)
1724 : #else
1725 1 : open(unit=grid_unt,file=trim(tmpfil),status='unknown',form='formatted')
1726 : #endif
1727 1 : rewind(grid_unt)
1728 1 : write(message,'(3a)') ch10," == Read real frequency grid from file ",trim(tmpfil)
1729 1 : call wrtout(std_out,message,'COLL')
1730 1 : write(message,'(5x,3a,i4)') 'Opened file : ',trim(tmpfil),' on unit ',grid_unt
1731 1 : call wrtout(std_out,message,'COLL')
1732 1 : read(grid_unt,*,iostat=ioerr) ngrid
1733 3 : ABI_MALLOC(paw_dmft%omega_r,(ngrid))
1734 3001 : do ifreq=1,ngrid
1735 3001 : read(grid_unt,*,iostat=ioerr) paw_dmft%omega_r(ifreq)
1736 : end do ! ifreq
1737 1 : close(grid_unt)
1738 : end if ! myproc=0
1739 4 : call xmpi_bcast(ioerr,0,xmpi_world,ierr)
1740 4 : if (ioerr /= 0) ABI_ERROR("Error when reading grid file")
1741 4 : call xmpi_bcast(ngrid,0,xmpi_world,ierr)
1742 4 : if (paw_dmft%myproc /= 0) then
1743 9 : ABI_MALLOC(paw_dmft%omega_r,(ngrid))
1744 : end if
1745 4 : call xmpi_bcast(paw_dmft%omega_r(:),0,xmpi_world,ierr)
1746 : else
1747 279 : ABI_MALLOC(paw_dmft%omega_r,(2*paw_dmft%dmft_nwr))
1748 : ! Set up real frequencies for spectral function in Hubbard one.
1749 93 : step = 0.00005_dp
1750 93 : paw_dmft%omega_r(2*paw_dmft%dmft_nwr) = pi * step * (two*dble(paw_dmft%dmft_nwr-1)+one)
1751 148800 : do ifreq=1,2*paw_dmft%dmft_nwr-1
1752 148800 : paw_dmft%omega_r(ifreq) = pi*step*(two*dble(ifreq-1)+one) - paw_dmft%omega_r(2*paw_dmft%dmft_nwr)
1753 : ! write(std_out,*) ifreq,paw_dmft%omega_r(ifreq)
1754 : end do ! ifreq
1755 :
1756 : end if ! iscf<0 and dmft_solv>=5 and dmft_solv<=8
1757 :
1758 97 : call init_paral_dmft(paw_dmft,paw_dmft%distrib_r,size(paw_dmft%omega_r(:)))
1759 :
1760 : !unit_e=2_dp
1761 :
1762 : ! paw_dmft%dmft_mag=0
1763 : ! do iatom=1,dtset%natom
1764 : ! do ii=1,3
1765 : ! if ( dtset(ii,iatom) > 0.001 ) paw_dmft%dmft_mag=1
1766 : ! enddo
1767 : ! enddo
1768 :
1769 97 : paw_dmft%gpu_option = dtset%gpu_option
1770 97 : paw_dmft%fermie_dft = fermie_dft ! in Ha
1771 97 : paw_dmft%fermie = fermie_dft
1772 :
1773 : !========================
1774 : !== Fixed self as input
1775 : !========================
1776 97 : paw_dmft%use_fixed_self = dtset%usedmatpu
1777 97 : paw_dmft%fixed_self => dmatpawu(:,:,:,:)
1778 :
1779 97 : if (paw_dmft%dmftcheck == -1) then
1780 0 : message = ' init_dmft: dmftcheck=-1 should not happen here'
1781 0 : ABI_BUG(message)
1782 : end if
1783 :
1784 485 : ABI_MALLOC(paw_dmft%eigen_dft,(paw_dmft%mbandc,paw_dmft%nkpt,paw_dmft%nsppol))
1785 679 : ABI_MALLOC(paw_dmft%chipsi,(paw_dmft%nspinor*(2*paw_dmft%maxlpawu+1),paw_dmft%mbandc,paw_dmft%nkpt,paw_dmft%nsppol,paw_dmft%natom))
1786 :
1787 97 : paw_dmft%lchipsiortho = 0
1788 :
1789 : !=========================================================
1790 : !== if we use ctqmc impurity solver
1791 : !=========================================================
1792 : ! IMPORTANT : paw_dmft%hybrid is corrupted somewhere in DMFT routines on
1793 : ! tikal_psc and max2_open64. Use a local hybrid in qmc_prep even if not optimal.
1794 : ! Anyway initializing ctqmc here is not good and produce the same result for
1795 : ! dmft_iter=1 which speed up the convergence ...
1796 : ! FIXME : Move this to init_sc_dmft and find bug
1797 97 : if (paw_dmft%dmft_solv == 5) then ! CTQMC initialisation
1798 : ! write(message,'(a,2x,a,f13.5)') ch10," == Initializing CTQMC"
1799 : ! call wrtout(std_out,message,'COLL')
1800 :
1801 180 : ABI_MALLOC(paw_dmft%hybrid,(paw_dmft%natom))
1802 104 : do iatom=1,paw_dmft%natom
1803 66 : if (paw_dmft%lpawu(iatom) == -1) cycle
1804 38 : nflavor = 2 * (2*paw_dmft%lpawu(iatom)+1)
1805 : #ifdef HAVE_MPI
1806 : call CtqmcInterface_init(paw_dmft%hybrid(iatom),paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1807 : & paw_dmft%dmftqmc_therm,paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1808 38 : & std_out,paw_dmft%dmftctqmc_chains,paw_dmft%spacecomm,nspinor=paw_dmft%nspinor)
1809 : #else
1810 : call CtqmcInterface_init(paw_dmft%hybrid(iatom),paw_dmft%dmftqmc_seed,paw_dmft%dmftqmc_n, &
1811 : & paw_dmft%dmftqmc_therm,paw_dmft%dmftctqmc_meas,nflavor,paw_dmft%dmftqmc_l,one/paw_dmft%temp,zero,&
1812 : & std_out,paw_dmft%dmftctqmc_chains,nspinor=paw_dmft%nspinor)
1813 : #endif
1814 : call CtqmcInterface_setOpts(paw_dmft%hybrid(iatom),&
1815 : & opt_Fk = 1,&
1816 : & opt_order = paw_dmft%dmftctqmc_order, &
1817 : & opt_histo = paw_dmft%dmftctqmc_localprop,&
1818 : & opt_movie = paw_dmft%dmftctqmc_mov, &
1819 : & opt_analysis = paw_dmft%dmftctqmc_correl,&
1820 : & opt_check = paw_dmft%dmftctqmc_check, &
1821 : & opt_noise = paw_dmft%dmftctqmc_grnns, &
1822 : & opt_spectra = paw_dmft%dmftctqmc_mrka, &
1823 104 : & opt_gmove = paw_dmft%dmftctqmc_gmove )
1824 : end do ! iatom
1825 : ! write(message,'(a,2x,a,f13.5)') ch10,&
1826 : !& " == Initialization CTQMC done"
1827 : !call wrtout(std_out,message,'COLL')
1828 : end if ! dmft_solv=5
1829 :
1830 : !************************************************************************
1831 97 : end subroutine init_dmft
1832 : !!***
1833 :
1834 : !!****f* m_paw_dmft/construct_nwli_dmft
1835 : !! NAME
1836 : !! construct_nwli_dmft
1837 : !!
1838 : !! FUNCTION
1839 : !! Compute linear frequencies
1840 : !!
1841 : !! INPUTS
1842 : !! paw_dmft=structure for dmft
1843 : !! nwli=number of linear frequencies
1844 : !!
1845 : !! OUTPUTS
1846 : !! omegali(1:nwli)=computed frequencies
1847 : !!
1848 : !! SOURCE
1849 : !!
1850 :
1851 134 : subroutine construct_nwli_dmft(paw_dmft,nwli,omega_li)
1852 :
1853 : !Arguments ------------------------------------
1854 : type(paw_dmft_type), intent(in) :: paw_dmft
1855 : integer, intent(in) :: nwli
1856 : real(dp), intent(inout) :: omega_li(:)
1857 : !Local variables-------------------------------
1858 : !fortran2003 ?
1859 : !real(dp), allocatable, intent(inout) :: omega_li(:)
1860 : integer :: ifreq
1861 : real(dp) :: factor
1862 : character(len=100) :: message
1863 : ! *********************************************************************
1864 :
1865 : ! if (allocated(omega_li)) then
1866 134 : if (size(omega_li(:)) /= nwli) then
1867 0 : write(message,'(2a,i8,a,i8)') ch10,"Number of linear frequencies asked is", &
1868 0 : & nwli,"whereas dimension of array omega_li is",size(omega_li(:))
1869 0 : ABI_BUG(message)
1870 : ! ABI_FREE(omega_li)
1871 : ! ABI_MALLOC(omega_li,(nwli))
1872 : ! write(*,*) "RESIZE"
1873 : ! call flush(6)
1874 : end if
1875 : ! write(*,*) "NOTHING"
1876 : ! call flush(6)
1877 : ! else
1878 : ! write(*,*) "ALLOCATE"
1879 : ! call flush(6)
1880 : ! ABI_MALLOC(omega_li,(nwli))
1881 : ! endif
1882 :
1883 : ! Set up linear frequencies
1884 134 : factor = pi * paw_dmft%temp
1885 19189084 : do ifreq=1,nwli
1886 19189084 : omega_li(ifreq) = factor * dble(2*ifreq-1)
1887 : ! (2(ifreq-1)+1 = 2ifreq-1
1888 : end do ! ifreq
1889 :
1890 134 : end subroutine construct_nwli_dmft
1891 : !!***
1892 :
1893 : !!****f* m_paw_dmft/construct_nwlo_dmft
1894 : !! NAME
1895 : !! construct_nwlo_dmft
1896 : !!
1897 : !! FUNCTION
1898 : !! Allocate log frequencies if used and compute them as well as their weight
1899 : !!
1900 : !! INPUTS
1901 : !! paw_dmft=structure for dmft calculation
1902 : !!
1903 : !!
1904 : !! SOURCE
1905 : !!
1906 : !! NOTE
1907 : !! The part of the code which deals
1908 : !! with the use of logarithmic frequencies
1909 : !! is a modification of the GNU GPL
1910 : !! code available on http://dmft.rutgers.edu/ and
1911 : !! described in the RMP paper written by
1912 : !! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti.
1913 : !!
1914 :
1915 42 : subroutine construct_nwlo_dmft(paw_dmft)
1916 :
1917 : use m_splines
1918 :
1919 : !Arguments ------------------------------------
1920 : type(paw_dmft_type), intent(inout) :: paw_dmft
1921 : !Local variables-------------------------------
1922 : integer :: deltaw,ifreq,ifreq2,myproc,nlin,nproc,nwli
1923 : integer :: nwlo,omegaBegin,omegaEnd,residu,spacecomm
1924 : character(len=10) :: tag
1925 : character(len=500) :: message
1926 : real(dp) :: deltaomega,expfac,omegamaxmin,prefacexp,temp,wl
1927 : complex(dp):: ybcbeg,ybcend
1928 42 : integer, allocatable :: select_log(:)
1929 42 : real(dp), allocatable :: omega_li(:),omega_lo_tmp(:),wgt_wlo(:)
1930 42 : complex(dp), allocatable :: splined_li(:),tospline_lo(:),ysplin2_lo(:)
1931 : ! *********************************************************************
1932 :
1933 42 : nwlo = paw_dmft%dmft_nwlo
1934 42 : nwli = paw_dmft%dmft_nwli
1935 42 : temp = paw_dmft%temp
1936 :
1937 : !== Variables for DMFT related to frequencies
1938 : ! the part of the code which deals
1939 : ! with the use of logarithmic frequencies
1940 : ! is a modification of the GNU GPL
1941 : ! code available on http://dmft.rutgers.edu/ and
1942 : ! described in the RMP paper written by
1943 : ! G.Kotliar, S.Y.Savrasov, K.Haule, V.S.Oudovenko, O.Parcollet, C.A.Marianetti
1944 :
1945 : !========================================
1946 : !== construct log. freq.
1947 42 : if (paw_dmft%dmft_log_freq == 1) then
1948 : !=======================================
1949 :
1950 126 : ABI_MALLOC(omega_lo_tmp,(nwlo))
1951 84 : ABI_MALLOC(wgt_wlo,(nwlo))
1952 : !cubic_freq=0
1953 : !omegamaxmin=paw_dmft%omega_li(paw_dmft%dmft_nwli)-paw_dmft%omega_li(paw_dmft%dmftqmc_l+1)
1954 42 : nlin = paw_dmft%dmftqmc_l ! number of linear frequencies
1955 42 : if (paw_dmft%dmft_solv == 6 .or. paw_dmft%dmft_solv == 7) nlin = int(nwlo*half)
1956 42 : omegamaxmin = pi * temp * two * dble(nwli-nlin-1)
1957 :
1958 : !if(cubic_freq==1) then
1959 :
1960 : ! if (paw_dmft%dmft_solv .eq. 5 ) then
1961 : ! write(message, '(2a)') ch10, "Warning : Cubish Mesh not tested with CT-QMC"
1962 : ! ABI_WARNING(message)
1963 : ! end if
1964 : ! ------------ CUBIC MESH MESH
1965 : ! useless
1966 : !nlin=dble(paw_dmft%dmft_nwli)
1967 : !nlog=dble(paw_dmft%dmft_nwlo)
1968 : !AA=((nlin-one)/nlin/(nlog**2-one)-one/(three*nlin))/((nlog**3-one)/(nlog**2-one)-seven/three)
1969 : !BB=(one/nlin - seven*AA)/three
1970 : !CC=-AA-BB
1971 : ! AA=((nlin-one)/nlin/(nlog-one)-one/(nlin))/((nlog**2-one)/(nlog-one)-three)
1972 : ! BB=(one/nlin - three*AA)
1973 : ! CC=-AA-BB
1974 : !write(message, '(a,16x,2(2x,a))') ch10," Cubic Mesh Parameters are"
1975 : !call wrtout(std_out,message,'COLL')
1976 : !write(message, '(3x,a,3(2x,e13.5))') "AA,BB,CC",AA,BB,CC
1977 : !call wrtout(std_out,message,'COLL')
1978 : !do ifreq=1,paw_dmft%dmft_nwlo
1979 : ! t1=dble(ifreq)
1980 : !omega_lo_tmp(ifreq)=(AA*t1**3+BB*t1**2+CC)*omegamaxmin+paw_dmft%omega_li(1)
1981 : ! omega_lo_tmp(ifreq)=(AA*t1**3+BB*t1**2+CC)*omegamaxmin+paw_dmft%temp*pi
1982 : ! paw_dmft%omega_lo(ifreq)=(AA*t1**2+BB*t1+CC)*omegamaxmin+paw_dmft%omega_li(1)
1983 : ! write(69,*) paw_dmft%omega_lo(ifreq),0.5
1984 : !enddo
1985 : !else
1986 42 : if (paw_dmft%dmft_solv < 4) then
1987 16 : paw_dmft%dmftqmc_l = 0
1988 16 : nlin = 0
1989 : end if
1990 :
1991 : ! ------------ LOGARITHMIC MESH
1992 42 : deltaomega = half
1993 42 : expfac = log(omegamaxmin/deltaomega) / (dble(nwlo-nlin-1)*half)
1994 42 : prefacexp = omegamaxmin / (exp(expfac*dble(nwlo-nlin-1))-one)
1995 126 : ABI_MALLOC(select_log,(nwlo))
1996 :
1997 : ! ------------ IMPOSE LINEAR MESH for w < 2*w_n=(2*l-1)pi/beta
1998 : ! Check variables (Already done in chkinp if dmft_solv==5)
1999 42 : if (nlin > nwlo) then
2000 0 : write(message,'(2a,i6)') ch10, &
2001 0 : & ' ERROR: dmft_nwlo has to be at least equal to 2xdmftqmc_l :',2*paw_dmft%dmftqmc_l
2002 0 : ABI_ERROR(message)
2003 : end if
2004 : ! End Check
2005 :
2006 42 : call construct_nwli_dmft(paw_dmft,nlin,omega_lo_tmp(1:nlin))
2007 17942 : select_log(1:nlin) = (/ (ifreq,ifreq=1,nlin) /)
2008 :
2009 : !do ifreq=1,paw_dmft%dmftqmc_l
2010 : ! omega_lo_tmp(ifreq)=(two*DBLE(ifreq-1)+one)*pi*paw_dmft%temp
2011 : ! select_log(ifreq)=ifreq
2012 : !enddo
2013 :
2014 : ! ------------ COMPLETE FREQUENCIES WITH LOG MESH
2015 42 : wl = temp * pi * dble(2*nlin+1)
2016 26861 : do ifreq=1,nwlo-nlin
2017 : !omega_lo_tmp(ifreq+paw_dmft%dmftqmc_l)=prefacexp*(exp(expfac*float(ifreq-1))-one)+paw_dmft%omega_li(paw_dmft%dmftqmc_l+1)
2018 26819 : omega_lo_tmp(ifreq+nlin) = prefacexp*(exp(expfac*dble(ifreq-1))-one) + wl
2019 : ! -------- Impose that each frequency of the logarithmic mesh is on a Matsubara frequency
2020 : ! FIXME : This may be done for all solver, not only for QMCs
2021 26861 : if (paw_dmft%dmft_solv >= 4) then
2022 : ! Compute the index "n" of iwn
2023 25967 : ifreq2 = nint((omega_lo_tmp(ifreq+nlin)/(temp*pi)-one)*half)
2024 : ! Compute freq
2025 25967 : omega_lo_tmp(ifreq+nlin) = (dble(ifreq2)*two+one) * pi * temp
2026 :
2027 25967 : if ((ifreq2+1) > nwli) then
2028 0 : write(message,'(2a,i8)') ch10,&
2029 0 : & ' BUG: init_dmft, dimension of array select_log is about to be overflown',(ifreq2+1)
2030 0 : ABI_BUG(message)
2031 : end if
2032 25967 : select_log(nlin+ifreq) = ifreq2 + 1
2033 : end if ! dmft_solv>=4
2034 : end do ! ifreq
2035 :
2036 : ! -------- Suppress duplicate frequencies
2037 : ! FIXME : So this also should be done for all solver and remove useless
2038 : ! frequencies
2039 42 : if (paw_dmft%dmft_solv >= 4) then
2040 26 : ifreq2 = 1
2041 34891 : do ifreq=2,nwlo-1
2042 34865 : if (select_log(ifreq2) == select_log(ifreq)) cycle
2043 21813 : ifreq2 = ifreq2 + 1
2044 21813 : omega_lo_tmp(ifreq2) = omega_lo_tmp(ifreq)
2045 34891 : select_log(ifreq2) = select_log(ifreq)
2046 : end do ! ifreq
2047 26 : paw_dmft%dmft_nwlo = ifreq2 + 1
2048 26 : nwlo = paw_dmft%dmft_nwlo
2049 : end if ! dmft_solv>=4
2050 :
2051 42 : omega_lo_tmp(1) = temp * pi
2052 42 : omega_lo_tmp(nwlo) = temp * pi * dble(2*nwli-1)
2053 :
2054 : !==================================
2055 : !== Construct weight for log. freq.
2056 : !==================================
2057 :
2058 126 : ABI_MALLOC(tospline_lo,(nwlo))
2059 126 : ABI_MALLOC(splined_li,(nwli))
2060 84 : ABI_MALLOC(ysplin2_lo,(nwlo))
2061 126 : ABI_MALLOC(omega_li,(nwli))
2062 42 : call construct_nwli_dmft(paw_dmft,nwli,omega_li(:))
2063 :
2064 : !Parallelization over frequencies!
2065 : ! ============= Set up =============
2066 42 : myproc = paw_dmft%myproc
2067 42 : nproc = paw_dmft%nproc
2068 42 : spacecomm = paw_dmft%spacecomm
2069 42 : deltaw = nwlo / nproc
2070 42 : residu = nwlo - nproc*deltaw
2071 42 : if (myproc < nproc-residu) then
2072 37 : omegaBegin = 1 + myproc*deltaw
2073 37 : omegaEnd = (myproc+1) * deltaw
2074 : else
2075 5 : omegaBegin = 1 + myproc*(deltaw+1) - nproc + residu
2076 5 : omegaEnd = omegaBegin + deltaw
2077 : end if
2078 :
2079 22759 : wgt_wlo(1:nwlo) = zero ! very important for xmpi_sum
2080 42 : ybcbeg = czero
2081 42 : ybcend = czero
2082 : ! ============= END Set up =============
2083 :
2084 22759 : tospline_lo(:) = czero
2085 :
2086 7216 : do ifreq=omegaBegin,omegaEnd
2087 : ! do ifreq1=1,paw_dmft%dmft_nwlo
2088 7174 : tospline_lo(ifreq) = cone
2089 : ! tospline_lo(ifreq1)=ifreq1**2-ifreq1
2090 : ! enddo
2091 : ! ybcbeg=cmplx(one/tol16**2,zero)
2092 : ! ybcend=cmplx(one/tol16**2,zero)
2093 :
2094 : !== spline delta function
2095 : call spline_complex(omega_lo_tmp(:),tospline_lo(:),nwlo, &
2096 7174 : & ybcbeg,ybcend,ysplin2_lo(:))
2097 : ! do ifreq1=1,paw_dmft%dmft_nwlo
2098 : ! write(6588,*) paw_dmft%omega_lo(ifreq1),ysplin2_lo(ifreq1)
2099 : ! enddo
2100 :
2101 : call splint_complex(nwlo,omega_lo_tmp(:),tospline_lo(:), &
2102 7174 : & ysplin2_lo(:),nwli,omega_li(:),splined_li(:))
2103 :
2104 7174 : tospline_lo(ifreq) = czero
2105 :
2106 : !== accumulate weights
2107 6259027216 : wgt_wlo(ifreq) = sum(dble(splined_li(:)))
2108 : ! do ifreq1=1,paw_dmft%dmft_nwlo
2109 : ! write(6688,*) paw_dmft%omega_lo(ifreq1),tospline_lo(ifreq1)
2110 : ! enddo
2111 :
2112 : ! do ifreq1=1,paw_dmft%dmft_nwli
2113 : ! write(6788,*) paw_dmft%omega_li(ifreq1),splined_li(ifreq1)
2114 :
2115 : end do ! ifreq
2116 : ! ============= Gatherall =============
2117 42 : call xmpi_sum(wgt_wlo(1:nwlo),spacecomm,residu)
2118 : ! ============= END Gatherall ==========
2119 : ! end parallelisation over frequencies
2120 :
2121 42 : ABI_FREE(tospline_lo)
2122 42 : ABI_FREE(splined_li)
2123 42 : ABI_FREE(ysplin2_lo)
2124 : ! if(abs(dtset%pawprtvol)>=3) then
2125 42 : write(message,'(a,18x,2(2x,a21))') ch10,"Log. Freq","weight"
2126 42 : call wrtout(std_out,message,'COLL')
2127 22759 : do ifreq=1,nwlo
2128 22717 : write(message,'(3x,a9,i6,2(2x,e21.14))') "--ifreq--",ifreq,omega_lo_tmp(ifreq),wgt_wlo(ifreq)
2129 22759 : call wrtout(std_out,message,'COLL')
2130 : end do ! ifreq
2131 42 : write(message,'(3x,a,i6)') " Total number of log frequencies is",nwlo
2132 42 : call wrtout(std_out,message,'COLL')
2133 42 : ifreq2 = 1
2134 1302 : do ifreq=1,min(30,nwlo)
2135 1260 : write(message,'(3x,a9,i6,2(2x,e21.14))') "--ifreq--",ifreq,omega_li(ifreq)
2136 1260 : call wrtout(std_out,message,'COLL')
2137 1260 : if (select_log(ifreq2) == ifreq) then
2138 780 : write(message,'(3x,a,i4,2(2x,i5))') "--sel_log",1
2139 780 : ifreq2 = ifreq + 1
2140 : else
2141 480 : write(message,'(3x,a,i4,2(2x,i5))') "--sel_log",0
2142 : end if
2143 1302 : call wrtout(std_out,message,'COLL')
2144 : end do ! ifreq
2145 42 : write(message,'(3x,2a)') "--ifreq--","..."
2146 42 : call wrtout(std_out,message,'COLL')
2147 42 : write(message,'(3x,a,i6,2(2x,e13.5))') "--ifreq--",nwli,omega_li(nwli)
2148 42 : call wrtout(std_out,message,'COLL')
2149 : ! endif
2150 42 : ABI_FREE(select_log)
2151 42 : ABI_FREE(omega_li)
2152 126 : ABI_MALLOC(paw_dmft%omega_lo,(nwlo))
2153 84 : ABI_MALLOC(paw_dmft%wgt_wlo,(nwlo))
2154 22759 : paw_dmft%omega_lo(1:nwlo) = omega_lo_tmp(1:nwlo)
2155 22759 : paw_dmft%wgt_wlo(1:nwlo) = wgt_wlo(1:nwlo)
2156 42 : ABI_FREE(omega_lo_tmp)
2157 42 : ABI_FREE(wgt_wlo)
2158 :
2159 : !=========================================================
2160 : !== do not construct log. freq. and use linear frequencies
2161 : !=========================================================
2162 : else
2163 :
2164 0 : ABI_MALLOC(paw_dmft%omega_lo,(nwlo))
2165 0 : ABI_MALLOC(paw_dmft%wgt_wlo,(nwlo))
2166 0 : write(tag,'(i10)') nwlo
2167 0 : write(message,'(4a)') ch10," Use of ",trim(adjustl(tag))," linear Matsubara frequencies for DMFT calculation"
2168 0 : call wrtout(std_out,message,'COLL')
2169 0 : call construct_nwli_dmft(paw_dmft,nwli,paw_dmft%omega_lo(:))
2170 0 : paw_dmft%wgt_wlo(:) = one
2171 : end if ! dmft_log_freq
2172 :
2173 : ! Should be check but since type definition does not initialize pointer with
2174 : ! =>null() (fortran95 and later) it produces conditional jump in valgrind
2175 : !if ( associated(paw_dmft%omega_lo) ) then
2176 : ! ABI_FREE(paw_dmft%omega_lo)
2177 : !endif
2178 : !if ( associated(paw_dmft%wgt_wlo) ) then
2179 : ! ABI_FREE(paw_dmft%wgt_wlo)
2180 : !endif
2181 :
2182 42 : end subroutine construct_nwlo_dmft
2183 : !!***
2184 :
2185 : !!****f* m_paw_dmft/destroy_dmft
2186 : !! NAME
2187 : !! destroy_dmft
2188 : !!
2189 : !! FUNCTION
2190 : !! deallocate some variables related to paw_dmft
2191 : !!
2192 : !! INPUTS
2193 : !! paw_dmft=structure for dmft calculation
2194 : !!
2195 : !! OUTPUT
2196 : !!
2197 : !! SOURCE
2198 :
2199 97 : subroutine destroy_dmft(paw_dmft)
2200 :
2201 : !Arguments ------------------------------------
2202 : type(paw_dmft_type), intent(inout) :: paw_dmft
2203 : !Local variables-------------------------------
2204 : integer :: iatom
2205 : ! *********************************************************************
2206 :
2207 97 : if (paw_dmft%dmft_solv == 5 .and. allocated(paw_dmft%hybrid)) then
2208 104 : do iatom=1,size(paw_dmft%hybrid) !paw_dmft%natom
2209 : !if(paw_dmft%lpawu(iatom)/=-1) then
2210 104 : call ctqmcinterface_finalize(paw_dmft%hybrid(iatom))
2211 : !endif
2212 : end do ! iatom
2213 104 : ABI_FREE(paw_dmft%hybrid)
2214 : end if
2215 97 : ABI_SFREE(paw_dmft%chipsi)
2216 97 : ABI_SFREE(paw_dmft%edc)
2217 97 : ABI_SFREE(paw_dmft%edcdc)
2218 97 : ABI_SFREE(paw_dmft%eigen_dft)
2219 97 : ABI_SFREE(paw_dmft%omega_r)
2220 97 : ABI_SFREE(paw_dmft%symrec_cart)
2221 97 : paw_dmft%eigen => null()
2222 97 : paw_dmft%fixed_self => null()
2223 97 : paw_dmft%indsym => null()
2224 97 : call destroy_paral_dmft(paw_dmft%distrib_r)
2225 :
2226 97 : end subroutine destroy_dmft
2227 : !!***
2228 :
2229 : !!****f* m_paw_dmft/destroy_sc_dmft
2230 : !! NAME
2231 : !! destroy_sc_dmft
2232 : !!
2233 : !! FUNCTION
2234 : !! deallocate paw_dmft
2235 : !!
2236 : !! INPUTS
2237 : !! paw_dmft=structure for dmft calculation
2238 : !!
2239 : !! OUTPUT
2240 : !!
2241 : !! SOURCE
2242 :
2243 5678 : subroutine destroy_sc_dmft(paw_dmft)
2244 :
2245 : !Arguments ------------------------------------
2246 : type(paw_dmft_type), intent(inout) :: paw_dmft
2247 : !Local variables-------------------------------
2248 : character(len=500) :: message
2249 : ! *********************************************************************
2250 :
2251 : if ((.not. allocated(paw_dmft%occnd) .or. .not. allocated(paw_dmft%band_in) &
2252 : & .or. .not. allocated(paw_dmft%include_bands) .or. .not. allocated(paw_dmft%exclude_bands)) &
2253 5678 : & .and. paw_dmft%use_dmft == 1) then
2254 : write(message,'(3a)') &
2255 0 : & ' an array is not allocated and is not deallocated with use_dmft==1 ',ch10, &
2256 0 : & ' Action : check the code'
2257 0 : ABI_WARNING(message)
2258 : end if
2259 5678 : ABI_SFREE(paw_dmft%occnd)
2260 5678 : ABI_SFREE(paw_dmft%band_in)
2261 5678 : ABI_SFREE(paw_dmft%include_bands)
2262 5678 : ABI_SFREE(paw_dmft%exclude_bands)
2263 5678 : ABI_SFREE(paw_dmft%siz_proj)
2264 5678 : ABI_SFREE(paw_dmft%zarot)
2265 5678 : ABI_SFREE(paw_dmft%phimtphi)
2266 5678 : ABI_SFREE(paw_dmft%phimtphi_int)
2267 5678 : ABI_SFREE(paw_dmft%phi_int)
2268 5678 : ABI_SFREE(paw_dmft%ylm)
2269 5678 : ABI_SFREE(paw_dmft%dpro)
2270 5678 : ABI_SFREE(paw_dmft%bessel)
2271 5678 : ABI_SFREE(paw_dmft%bessel_int)
2272 5678 : ABI_SFREE(paw_dmft%lpawu)
2273 5678 : ABI_SFREE(paw_dmft%omega_lo)
2274 5678 : ABI_SFREE(paw_dmft%wgt_wlo)
2275 5678 : ABI_SFREE(paw_dmft%slm2ylm)
2276 5678 : ABI_SFREE(paw_dmft%jmj2ylm)
2277 :
2278 5678 : paw_dmft%nband => null()
2279 5678 : paw_dmft%dmft_shiftself => null()
2280 5678 : paw_dmft%dmft_nominal => null()
2281 5678 : paw_dmft%npwarr => null()
2282 : ! paw_dmft%wtk is only an explicit pointer =>dtset%wtk
2283 : ! if (associated(paw_dmft%wtk)) deallocate(paw_dmft%wtk)
2284 5678 : paw_dmft%wtk => null()
2285 5678 : paw_dmft%typat => null()
2286 5678 : paw_dmft%int_meshsz => null()
2287 :
2288 5678 : if (allocated(paw_dmft%radgrid)) then
2289 43 : call pawrad_free(paw_dmft%radgrid(:))
2290 111 : ABI_FREE(paw_dmft%radgrid)
2291 : end if
2292 5678 : call destroy_sc_dmft_paralkgb(paw_dmft)
2293 5678 : if (paw_dmft%use_dmft == 1) then
2294 43 : call destroy_paral_dmft(paw_dmft%distrib)
2295 : end if
2296 :
2297 5678 : end subroutine destroy_sc_dmft
2298 : !!***
2299 :
2300 : !!****f* m_paw_dmft/print_dmft
2301 : !! NAME
2302 : !! print_dmft
2303 : !!
2304 : !! FUNCTION
2305 : !! Print relevant data for DMFT cycle.
2306 : !!
2307 : !! INPUTS
2308 : !! paw_dmft=structure for dmft calculation
2309 : !! pawprtvol=flag for print
2310 : !!
2311 : !! OUTPUT
2312 : !!
2313 : !! SOURCE
2314 :
2315 97 : subroutine print_dmft(paw_dmft,pawprtvol)
2316 :
2317 : !Arguments ------------------------------------
2318 : integer, intent(in) :: pawprtvol
2319 : type(paw_dmft_type), intent(in) :: paw_dmft
2320 : !Local variables-------------------------------
2321 : integer :: iband,ifreq,ikpt,isppol
2322 : character(len=500) :: message
2323 : ! *********************************************************************
2324 :
2325 97 : if (abs(pawprtvol) >= 3 ) then
2326 : write(message,'(4a,3(a,2x,e21.14,a))') &
2327 73 : & " -------------------------------------------------",ch10,&
2328 73 : & " --- Data for DMFT ",ch10,&
2329 73 : & " --- fermie = ",paw_dmft%fermie ,ch10,&
2330 73 : & " --- fermie_dft = ",paw_dmft%fermie_dft,ch10,&
2331 146 : & " --- temp = ",paw_dmft%temp ,ch10
2332 73 : call wrtout(std_out,message,'COLL')
2333 73 : if (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7) then
2334 : write(message,'(7(a,15x,i8,a),a,2x,e21.14,3a)') &
2335 73 : & " --- natpawu = ",paw_dmft%natpawu ,ch10,&
2336 73 : & " --- dmft_iter = ",paw_dmft%dmft_iter ,ch10,&
2337 73 : & " --- dmft_solv = ",paw_dmft%dmft_solv ,ch10,&
2338 73 : & " --- dmft_nwlo = ",paw_dmft%dmft_nwlo ,ch10,&
2339 73 : & " --- dmft_nwli = ",paw_dmft%dmft_nwli ,ch10,&
2340 73 : & " --- dmft_dc = ",paw_dmft%dmft_dc ,ch10,&
2341 73 : & " --- dmftqmc_l = ",paw_dmft%dmftqmc_l ,ch10,&
2342 73 : & " --- dmftqmc_n = ",paw_dmft%dmftqmc_n ,ch10,&
2343 146 : & " -------------------------------------------------",ch10
2344 : else
2345 : write(message,'(7(a,15x,i8,a),2a)') &
2346 0 : & " --- natpawu = ",paw_dmft%natpawu ,ch10,&
2347 0 : & " --- dmft_iter = ",paw_dmft%dmft_iter ,ch10,&
2348 0 : & " --- dmft_solv = ",paw_dmft%dmft_solv ,ch10,&
2349 0 : & " --- n_iw = ",paw_dmft%dmft_nwli ,ch10,&
2350 0 : & " --- dmft_dc = ",paw_dmft%dmft_dc ,ch10,&
2351 0 : & " --- n_tau = ",paw_dmft%dmftqmc_l ,ch10,&
2352 0 : & " --- n_cycles = ",paw_dmft%dmft_triqs_n_cycles ,ch10,&
2353 0 : & " -------------------------------------------------",ch10
2354 : end if
2355 73 : call wrtout(std_out,message,'COLL')
2356 :
2357 : ! write(message,'(4a,3(a,2x,f8.3,a),8(a,2x,i8,a),a)') "-----------------------------------------------",ch10,&
2358 : !& "--- Data for DMFT ",ch10,&
2359 : !& "--- paw_dmft%fermie = ",paw_dmft%fermie ,ch10,&
2360 : !& "--- paw_dmft%fermie_dft = ",paw_dmft%fermie_dft,ch10,&
2361 : !& "--- paw_dmft%temp = ",paw_dmft%temp ,ch10,&
2362 : !& "--- paw_dmft%natpawu = ",paw_dmft%natpawu ,ch10,&
2363 : !& "--- paw_dmft%dmft_iter = ",paw_dmft%dmft_iter ,ch10,&
2364 : !& "--- paw_dmft%dmft_solv = ",paw_dmft%dmft_solv ,ch10,&
2365 : !& "--- paw_dmft%dmft_nwlo = ",paw_dmft%dmft_nwlo ,ch10,&
2366 : !& "--- paw_dmft%dmft_nwli = ",paw_dmft%dmft_nwli ,ch10,&
2367 : !& "--- paw_dmft%dmft_dc = ",paw_dmft%dmft_dc ,ch10,&
2368 : !& "--- paw_dmft%dmftqmc_l = ",paw_dmft%dmftqmc_l ,ch10,&
2369 : !& "--- paw_dmft%dmftqmc_n = ",paw_dmft%dmftqmc_n ,ch10,&
2370 : !& "-----------------------------------------------"
2371 73 : if (abs(pawprtvol) > 10) then
2372 0 : call wrtout(std_out,message,'COLL')
2373 0 : write(message,'(a)') " DFT Eigenvalues "
2374 0 : do isppol=1,paw_dmft%nsppol
2375 0 : write(message,'(a,i4)') "--isppol--",isppol
2376 0 : call wrtout(std_out,message,'COLL')
2377 0 : do ikpt=1,paw_dmft%nkpt
2378 0 : write(message,'(a,i4,2x,f14.5,a)') " -k-pt--",ikpt,paw_dmft%wtk(ikpt),"(<-weight(k-pt))"
2379 0 : call wrtout(std_out,message,'COLL')
2380 0 : do iband=1,paw_dmft%mbandc
2381 0 : write(message,'(a,i4,f10.5)') " -iband--",iband,paw_dmft%eigen_dft(iband,ikpt,isppol)
2382 0 : call wrtout(std_out,message,'COLL')
2383 : end do ! iband
2384 : end do ! ikpt
2385 : end do ! isppol
2386 0 : write(message,'(3x,a)') "Log. Freq"
2387 0 : call wrtout(std_out,message,'COLL')
2388 0 : do ifreq=1,paw_dmft%dmft_nwlo
2389 0 : write(message,'(3x,a,i4,2(2x,e13.5))') "--ifreq--",ifreq,paw_dmft%omega_lo(ifreq),paw_dmft%wgt_wlo(ifreq)
2390 0 : call wrtout(std_out,message,'COLL')
2391 : end do ! ifreq
2392 : end if ! abs(pawprtvol)>10
2393 : end if ! abs(pawprtvol)>=3
2394 :
2395 97 : end subroutine print_dmft
2396 : !!***
2397 :
2398 : !!****f* m_paw_dmft/print_sc_dmft
2399 : !! NAME
2400 : !! print_sc_dmft
2401 : !!
2402 : !! FUNCTION
2403 : !! Print relevant data for self-consistent
2404 : !! DFT+DMFT cycle.
2405 : !!
2406 : !! INPUTS
2407 : !! paw_dmft=structure for dmft calculation
2408 : !! pawprtvol=flag for print
2409 : !!
2410 : !! OUTPUT
2411 : !!
2412 : !! SOURCE
2413 :
2414 43 : subroutine print_sc_dmft(paw_dmft,pawprtvol)
2415 :
2416 : !Arguments ------------------------------------
2417 : integer, intent(in) :: pawprtvol
2418 : type(paw_dmft_type), intent(in) :: paw_dmft
2419 : !Local variables-------------------------------
2420 : integer :: iband
2421 : character(len=500) :: message
2422 : ! *********************************************************************
2423 :
2424 43 : if (abs(pawprtvol) >= 3) then
2425 29 : write(message,'(5a,7(a,2x,i5,a),a)') ch10, &
2426 29 : & " -----------------------------------------------",ch10,&
2427 29 : & " --- Data for self-consistent DFT+DMFT cycle",ch10,&
2428 29 : & " --- mband = ",paw_dmft%mband,ch10,&
2429 29 : & " --- dmftbandi = ",paw_dmft%dmftbandi,ch10,&
2430 29 : & " --- dmftbandf = ",paw_dmft%dmftbandf,ch10,&
2431 29 : & " --- nb of corr. bands = ",paw_dmft%mbandc,ch10,&
2432 29 : & " --- nkpt = ",paw_dmft%nkpt,ch10,&
2433 29 : & " --- nsppol = ",paw_dmft%nsppol,ch10,&
2434 29 : & " --- usedmft = ",paw_dmft%use_dmft,ch10,&
2435 : !& " --- use_sc_dmft = ",paw_dmft%use_sc_dmft,ch10,&
2436 58 : & " -----------------------------------------------"
2437 29 : call wrtout(std_out,message,'COLL')
2438 29 : write(message,'(2a)') ch10," Indicating correlated bands"
2439 29 : call wrtout(std_out,message,'COLL')
2440 670 : write(message,'(100i5)') (iband,iband=1,min(paw_dmft%mband,100))
2441 29 : call wrtout(std_out,message,'COLL')
2442 670 : write(message,'(100L5)') (paw_dmft%band_in(iband),iband=1,min(paw_dmft%mband,100))
2443 29 : call wrtout(std_out,message,'COLL')
2444 29 : write(message,'(2a)') ch10," Correlated index Band index"
2445 29 : call wrtout(std_out,message,'COLL')
2446 310 : do iband=1,paw_dmft%mbandc
2447 281 : write(message,'(5x,i5,10x,i5)') iband,paw_dmft%include_bands(iband)
2448 310 : call wrtout(std_out,message,'COLL')
2449 : end do ! iband
2450 29 : if (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7) then
2451 29 : write(message,'(2a,i4,a)') ch10,&
2452 29 : & 'The',paw_dmft%mband-paw_dmft%dmftbandf+paw_dmft%dmftbandi-1,&
2453 58 : & ' Following bands are excluded from the DMFT calculation'
2454 29 : call wrtout(std_out,message,'COLL')
2455 389 : write(message,'(100i5)') (paw_dmft%exclude_bands(iband),iband=1,min(paw_dmft%mband-paw_dmft%dmftbandf+paw_dmft%dmftbandi-1,100))
2456 29 : call wrtout(std_out,message,'COLL')
2457 : end if
2458 29 : write(message,*)
2459 29 : call wrtout(std_out,message,'COLL')
2460 : end if ! abs(pawprtvol)>=3
2461 :
2462 43 : end subroutine print_sc_dmft
2463 : !!***
2464 :
2465 : !!****f* m_paw_dmft/saveocc_dmft
2466 : !! NAME
2467 : !! saveocc_dmft
2468 : !!
2469 : !! FUNCTION
2470 : !! save occnd on disk
2471 : !!
2472 : !! INPUTS
2473 : !! paw_dmft=structure for dmft calculation
2474 : !!
2475 : !! OUTPUT
2476 : !!
2477 : !! SOURCE
2478 :
2479 42 : subroutine saveocc_dmft(paw_dmft)
2480 :
2481 : !Arguments ------------------------------------
2482 : type(paw_dmft_type), intent(in) :: paw_dmft
2483 : !Local variables-------------------------------
2484 : integer :: ib,ib1,ikpt,is,nband_k,nkpt,unitsaveocc
2485 : character(len=500) :: message
2486 : character(len=fnlen) :: tmpfil
2487 : ! ********************************************************************
2488 :
2489 42 : tmpfil = trim(paw_dmft%filapp)//'_DMFTOCCND'
2490 42 : if (open_file(tmpfil,message,newunit=unitsaveocc,status='unknown',form='formatted') /= 0) ABI_ERROR(message)
2491 :
2492 42 : rewind(unitsaveocc)
2493 42 : write(message,'(2a)') ch10," == Print DFT+DMFT non diagonal occupations on disk"
2494 42 : call wrtout(std_out,message,'COLL')
2495 42 : write(message,'(5a,2x,4i5,2a)') "# DFT+DMFT off-diagonal occupations f_{ib,ib1} = <Psi^{dagger}_{ib1}|Psi_{ib}>", &
2496 42 : & ch10,"# natom,nsppol,mbandc,nkpt",ch10, &
2497 42 : & "####",paw_dmft%natom,paw_dmft%nsppol,paw_dmft%mbandc,paw_dmft%nkpt,ch10, &
2498 84 : & "# isppol ikpt ib ib1 Re Imag"
2499 :
2500 42 : call wrtout(unitsaveocc,message,'COLL')
2501 42 : nkpt = paw_dmft%nkpt
2502 107 : do is=1,paw_dmft%nsppol
2503 579 : do ikpt=1,nkpt
2504 472 : nband_k = paw_dmft%nband(ikpt+(is-1)*nkpt)
2505 11333 : do ib=1,nband_k
2506 10796 : if ((.not. paw_dmft%band_in(ib)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
2507 106592 : do ib1=1,nband_k
2508 101944 : if ((.not. paw_dmft%band_in(ib1)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
2509 56568 : write(unitsaveocc,*) is,ikpt,ib,ib1,paw_dmft%occnd(1,ib,ib1,ikpt,is),&
2510 169308 : & paw_dmft%occnd(2,ib,ib1,ikpt,is)
2511 : end do ! ib1
2512 : end do ! ib
2513 : end do ! ikpt
2514 : end do ! is
2515 42 : write(message,'(3a)') "# end of record",ch10,"#### 1234 "
2516 42 : call wrtout(unitsaveocc,message,'COLL')
2517 42 : close(unitsaveocc)
2518 :
2519 42 : end subroutine saveocc_dmft
2520 : !!***
2521 :
2522 : !!****f* m_paw_dmft/readocc_dmft
2523 : !! NAME
2524 : !! readocc_dmft
2525 : !!
2526 : !! FUNCTION
2527 : !! read occnd on disk
2528 : !!
2529 : !! INPUTS
2530 : !! paw_dmft=structure for dmft calculation
2531 : !! filnam_ds3 = root for filname to read (input)
2532 : !! filnam_ds4 = root for filname to read (output)
2533 : !!
2534 : !! OUTPUT
2535 : !! paw_dmft: occnd
2536 : !!
2537 : !! SOURCE
2538 :
2539 22 : subroutine readocc_dmft(paw_dmft,filnam_ds3,filnam_ds4)
2540 :
2541 : !Arguments ------------------------------------
2542 : type(paw_dmft_type), intent(inout) :: paw_dmft
2543 : character(len=fnlen), intent(in) :: filnam_ds3,filnam_ds4
2544 : !Local variables-------------------------------
2545 : integer :: dum1,dum2,dum3,dum4,ib,ib1,ikpt,ioerr,is,nband_k,nkpt,unitsaveocc
2546 : logical :: lexist
2547 : character(len=4) :: chtemp
2548 : character(len=500) :: message
2549 : character(len=fnlen) :: tmpfil
2550 : ! *********************************************************************
2551 :
2552 22 : if (paw_dmft%dmft_read_occnd == 0) return
2553 2 : if (paw_dmft%dmft_read_occnd == 1) tmpfil = trim(filnam_ds3)//'_DMFTOCCND'
2554 2 : if (paw_dmft%dmft_read_occnd == 2) tmpfil = trim(filnam_ds4)//'_DMFTOCCND'
2555 2 : inquire(file=trim(tmpfil),exist=lexist)!,recl=nrecl)
2556 2 : unitsaveocc = 679
2557 2 : if (lexist) then
2558 0 : if (open_file(tmpfil,message,unit=unitsaveocc,status='unknown',form='formatted') /= 0) ABI_ERROR(message)
2559 0 : rewind(unitsaveocc)
2560 0 : write(message,'(3a)') ch10," == Read DMFT non diagonal occupations on disk"
2561 0 : call wrtout(std_out,message,'COLL')
2562 0 : read(unitsaveocc,*)
2563 0 : read(unitsaveocc,*)
2564 0 : read(unitsaveocc,*,iostat=ioerr) chtemp,dum1,dum2,dum3,dum4
2565 0 : read(unitsaveocc,*)
2566 0 : if (ioerr < 0) write(std_out,*) "read",dum1,dum2,dum3,dum4
2567 0 : write(message,'(2a,4i4)') ch10," == natom, nsppol, nbandc, nkpt read are",dum1,dum2,dum3,dum4
2568 0 : call wrtout(std_out,message,'COLL')
2569 0 : nkpt = paw_dmft%nkpt
2570 0 : do is=1,paw_dmft%nsppol
2571 0 : do ikpt=1,nkpt
2572 0 : nband_k = paw_dmft%nband(ikpt+(is-1)*nkpt)
2573 0 : do ib=1,nband_k
2574 0 : if ((.not. paw_dmft%band_in(ib)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
2575 0 : do ib1=1,nband_k
2576 0 : if ((.not. paw_dmft%band_in(ib1)) .and. (paw_dmft%dmft_solv /= 6 .and. paw_dmft%dmft_solv /= 7)) cycle
2577 0 : read(unitsaveocc,*) dum1,dum2,dum3,dum4,&
2578 0 : & paw_dmft%occnd(1,ib,ib1,ikpt,is),paw_dmft%occnd(2,ib,ib1,ikpt,is)
2579 : end do ! ib1
2580 : end do ! ib
2581 : end do ! ikpt
2582 : end do ! is
2583 : ! write(read,'(3a)') "# end of record",ch10&
2584 : !& ,"#### 1234 "
2585 : ! call wrtout(unitsaveocc,message,'COLL')
2586 : else
2587 2 : write(message,'(2a,2x,2a)') ch10," File",trim(tmpfil),"is not available"
2588 2 : call wrtout(std_out,message,'COLL')
2589 2 : write(message,'(4a)') ch10," ==> DMFT Occupations not available for restart", ch10, &
2590 4 : & " -> The calculation is started with Fermi Dirac scheme for occupations"
2591 2 : call wrtout(std_out,message,'COLL')
2592 : end if ! lexist
2593 :
2594 : end subroutine readocc_dmft
2595 : !!***
2596 :
2597 : !!****f* m_paw_dmft/init_sc_dmft_paralkgb
2598 : !! NAME
2599 : !! init_sc_dmft_paralkgb
2600 : !!
2601 : !! FUNCTION
2602 : !! Init some values used with KGB parallelism in self consistent DMFT
2603 : !! calculation.
2604 : !!
2605 : !! INPUTS
2606 : !! paw_dmft=structure for dmft calculation
2607 : !! mpi_enreg = information about MPI parallelization
2608 : !!
2609 : !! OUTPUT
2610 : !! paw_dmft: bandc_proc, use_bandc
2611 : !!
2612 : !! SOURCE
2613 :
2614 4 : subroutine init_sc_dmft_paralkgb(paw_dmft,mpi_enreg)
2615 :
2616 : !Arguments ------------------------------------
2617 : type(paw_dmft_type), intent(inout) :: paw_dmft
2618 : type(MPI_type), intent(in) :: mpi_enreg
2619 : !Local variables-------------------------------
2620 : integer :: ib,ibc,mbandc,nproc,proc
2621 : ! *********************************************************************
2622 :
2623 4 : mbandc = paw_dmft%mbandc
2624 4 : nproc = mpi_enreg%nproc_band
2625 :
2626 12 : ABI_MALLOC(paw_dmft%bandc_proc,(mbandc))
2627 12 : ABI_MALLOC(paw_dmft%use_bandc,(nproc))
2628 52 : paw_dmft%bandc_proc(:) = 0
2629 20 : paw_dmft%use_bandc(:) = .false.
2630 :
2631 52 : do ibc=1,mbandc
2632 48 : ib = paw_dmft%include_bands(ibc)
2633 48 : proc = mod((ib-1)/mpi_enreg%bandpp,nproc)
2634 48 : paw_dmft%bandc_proc(ibc) = proc
2635 52 : paw_dmft%use_bandc(proc+1) = .true.
2636 : end do ! ibc
2637 :
2638 4 : end subroutine init_sc_dmft_paralkgb
2639 :
2640 : !!***
2641 :
2642 : !!****f* m_paw_dmft/destroy_sc_dmft_paralkgb
2643 : !! NAME
2644 : !! destroy_sc_dmft_paralkgb
2645 : !!
2646 : !! FUNCTION
2647 : !! deallocate bandc_proc and use_bandc
2648 : !!
2649 : !! INPUTS
2650 : !! paw_dmft=structure for dmft calculation
2651 : !!
2652 : !! OUTPUT
2653 : !!
2654 : !! SOURCE
2655 :
2656 5678 : subroutine destroy_sc_dmft_paralkgb(paw_dmft)
2657 :
2658 : !Arguments ------------------------------------
2659 : type(paw_dmft_type), intent(inout) :: paw_dmft
2660 : ! *********************************************************************
2661 :
2662 5678 : ABI_SFREE(paw_dmft%bandc_proc)
2663 5678 : ABI_SFREE(paw_dmft%use_bandc)
2664 :
2665 5678 : end subroutine destroy_sc_dmft_paralkgb
2666 : !!***
2667 :
2668 : !!****f* m_paw_dmft/init_paral_dmft
2669 : !! NAME
2670 : !! init_paral_dmft
2671 : !!
2672 : !! FUNCTION
2673 : !! Initialize MPI distribution
2674 : !!
2675 : !! INPUTS
2676 : !! paw_dmft=structure for dmft calculation
2677 : !! distrib <type(mpi_distrib_dmft_type)> = mpi related data
2678 : !! nfreq=number of frequencies over which to parallelize
2679 : !!
2680 : !! OUTPUT
2681 : !!
2682 : !! SOURCE
2683 :
2684 140 : subroutine init_paral_dmft(paw_dmft,distrib,nfreq)
2685 :
2686 : !Arguments ------------------------------------
2687 : type(paw_dmft_type), intent(in) :: paw_dmft
2688 : type(mpi_distrib_dmft_type), intent(inout) :: distrib
2689 : integer, intent(in) :: nfreq
2690 : !Local variables-------------------------------
2691 : integer :: deltakpt,deltaw,ierr,ifreq,ikpt,irank,myproc,nfreq_proc
2692 : integer :: nkpt,nkpt_proc,nproc,nproc_freq,nproc_kpt,residu,spacecomm
2693 : ! *********************************************************************
2694 :
2695 140 : myproc = paw_dmft%myproc
2696 140 : nkpt = paw_dmft%nkpt
2697 140 : nproc = paw_dmft%nproc
2698 140 : spacecomm = paw_dmft%spacecomm
2699 :
2700 140 : distrib%nw = nfreq
2701 140 : distrib%shiftk = 0
2702 :
2703 420 : ABI_MALLOC(distrib%nw_mem,(nproc))
2704 420 : ABI_MALLOC(distrib%procb,(nkpt))
2705 420 : ABI_MALLOC(distrib%procf,(nfreq))
2706 280 : ABI_MALLOC(distrib%proct,(nfreq))
2707 :
2708 : ! First initialize parallelization over frequencies only
2709 :
2710 140 : deltaw = nfreq / nproc
2711 140 : residu = nfreq - deltaw*nproc
2712 :
2713 140 : nproc_freq = min(nfreq,nproc)
2714 :
2715 141 : if (nproc_freq < nproc) distrib%nw_mem(nproc_freq+1:nproc) = 0
2716 : ifreq = 1
2717 579 : do irank=0,nproc_freq-1
2718 439 : nfreq_proc = merge(deltaw+1,deltaw,irank<residu)
2719 439 : distrib%nw_mem(irank+1) = nfreq_proc
2720 183956 : distrib%procf(ifreq:ifreq+nfreq_proc-1) = irank
2721 140 : ifreq = ifreq + nfreq_proc
2722 : end do ! irank
2723 :
2724 : ! Next initialize parallelization on kpt and then frequencies
2725 :
2726 140 : nproc_kpt = min(nkpt,nproc)
2727 140 : nproc_freq = nproc / nkpt
2728 :
2729 280 : ABI_MALLOC(distrib%nkpt_mem,(nproc))
2730 420 : ABI_MALLOC(distrib%nw_mem_kptparal,(nproc_freq+1))
2731 :
2732 152 : if (nproc > nproc_kpt) distrib%nkpt_mem(nproc_kpt+1:nproc) = 0
2733 140 : distrib%nw_mem_kptparal(nproc_freq+1) = 0
2734 :
2735 140 : if (nproc_freq <= 1) then ! parallelization on kpt only
2736 :
2737 140 : deltakpt = nkpt / nproc
2738 140 : residu = nkpt - deltakpt*nproc
2739 :
2740 140 : distrib%nw_mem_kptparal(1) = nfreq
2741 183657 : distrib%proct(:) = 0
2742 :
2743 : ikpt = 1
2744 568 : do irank=0,nproc_kpt-1
2745 428 : nkpt_proc = merge(deltakpt+1,deltakpt,irank<residu)
2746 428 : distrib%nkpt_mem(irank+1) = nkpt_proc
2747 1622 : distrib%procb(ikpt:ikpt+nkpt_proc-1) = irank
2748 428 : if (myproc == irank) distrib%shiftk = ikpt - 1
2749 140 : ikpt = ikpt + nkpt_proc
2750 : end do ! irank
2751 :
2752 140 : distrib%comm_kpt = spacecomm
2753 140 : distrib%me_kpt = mod(myproc,nproc_kpt)
2754 140 : distrib%me_freq = myproc / nproc_kpt
2755 :
2756 140 : call xmpi_comm_split(spacecomm,distrib%me_kpt,distrib%me_freq,distrib%comm_freq,ierr)
2757 :
2758 140 : distrib%me_kpt = myproc
2759 :
2760 : else ! parallelization on both kpt and frequencies
2761 :
2762 0 : deltaw = nfreq / nproc_freq
2763 0 : residu = nfreq - deltaw*nproc_freq
2764 :
2765 0 : distrib%nkpt_mem(1:nkpt) = 1
2766 :
2767 0 : do ikpt=1,nkpt
2768 0 : distrib%procb(ikpt) = ikpt - 1
2769 : end do ! ikpt
2770 :
2771 : ifreq = 1
2772 0 : do irank=0,nproc_freq-1
2773 0 : nfreq_proc = merge(deltaw+1,deltaw,irank<residu)
2774 0 : if (nfreq_proc > 0) distrib%proct(ifreq:ifreq+nfreq_proc-1) = irank
2775 0 : ifreq = ifreq + nfreq_proc
2776 0 : distrib%nw_mem_kptparal(irank+1) = nfreq_proc
2777 : end do ! irank
2778 :
2779 0 : distrib%me_kpt = myproc / nproc_freq
2780 0 : distrib%me_freq = mod(myproc,nproc_freq)
2781 0 : distrib%shiftk = distrib%me_kpt
2782 :
2783 0 : call xmpi_comm_split(spacecomm,distrib%me_freq,distrib%me_kpt,distrib%comm_kpt,ierr)
2784 :
2785 0 : if (myproc >= nkpt*nproc_freq) then
2786 0 : distrib%me_kpt = myproc - nkpt*nproc_freq
2787 0 : distrib%me_freq = nproc_freq
2788 : end if
2789 :
2790 0 : call xmpi_comm_split(spacecomm,distrib%me_kpt,distrib%me_freq,distrib%comm_freq,ierr)
2791 :
2792 0 : if (myproc >= nkpt*nproc_freq) distrib%me_kpt = myproc / nproc_freq
2793 :
2794 : end if ! nproc_freq<=1
2795 :
2796 140 : end subroutine init_paral_dmft
2797 : !!***
2798 :
2799 : !!****f* m_paw_dmft/destroy_paral_dmft
2800 : !! NAME
2801 : !! destroy_paral_dmft
2802 : !!
2803 : !! FUNCTION
2804 : !! Deallocate MPI distribution
2805 : !!
2806 : !! INPUTS
2807 : !! paw_dmft=structure for dmft calculation
2808 : !! distrib <type(mpi_distrib_dmft_type)> = mpi related data
2809 : !!
2810 : !! OUTPUT
2811 : !!
2812 : !! SOURCE
2813 :
2814 140 : subroutine destroy_paral_dmft(distrib)
2815 :
2816 : !Arguments ------------------------------------
2817 : type(mpi_distrib_dmft_type), intent(inout) :: distrib
2818 : !Local variables-------------------------------
2819 : ! *********************************************************************
2820 :
2821 140 : ABI_FREE(distrib%nkpt_mem)
2822 140 : ABI_FREE(distrib%nw_mem)
2823 140 : ABI_FREE(distrib%nw_mem_kptparal)
2824 140 : ABI_FREE(distrib%procb)
2825 140 : ABI_FREE(distrib%procf)
2826 140 : ABI_FREE(distrib%proct)
2827 :
2828 140 : call xmpi_comm_free(distrib%comm_freq)
2829 140 : call xmpi_comm_free(distrib%comm_kpt)
2830 :
2831 140 : end subroutine destroy_paral_dmft
2832 : !!***
2833 :
2834 0 : END MODULE m_paw_dmft
2835 : !!***
|
