Branch data 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, 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, contiguous, pointer :: nband(:) => null()
670 : : ! Number of bands for each k-point and polarization
671 : :
672 : : integer, contiguous, pointer :: npwarr(:) => null()
673 : : ! Number of plane waves on current process for each k-point
674 : :
675 : : integer, contiguous, pointer :: typat(:) => null()
676 : : ! Type of each atom
677 : :
678 : : real(dp), contiguous, pointer :: dmft_shiftself(:) => null()
679 : : ! Initial shift of the self-energy for each atom
680 : :
681 : : real(dp), 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), 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 : 5662 : subroutine init_sc_dmft(dtset,mpsang,paw_dmft,gprimd,kg,mpi_enreg,npwarr,occ,pawang, &
743 : 11324 : & 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 : 5712 : integer, allocatable :: ind_msml(:,:)
774 : 5712 : logical, allocatable :: lcycle(:),typcycle(:)
775 : 5712 : real(dp), allocatable :: rmax(:),kpg(:,:),kpg_norm(:)
776 : : character(len=500) :: dc_string,lda_string,message
777 : : !************************************************************************
778 : :
779 : 5712 : mband = dtset%mband
780 : 5712 : mkmem = dtset%mkmem
781 : 5712 : natom = dtset%natom
782 : 5712 : nkpt = dtset%nkpt
783 : 5712 : nspinor = dtset%nspinor
784 : 5712 : nsppol = dtset%nsppol
785 : 5712 : nsym = dtset%nsym
786 : 5712 : ntypat = dtset%ntypat
787 : :
788 : 5712 : dmftbandi = dtset%dmftbandi
789 : 5712 : dmftbandf = dtset%dmftbandf
790 : 5712 : dmft_dc = dtset%dmft_dc
791 : 5712 : dmft_solv = dtset%dmft_solv
792 : 5712 : off_diag = dtset%dmft_triqs_off_diag
793 : 5712 : use_dmft = abs(dtset%usedmft)
794 : 5712 : paw_dmft%use_dmft = use_dmft
795 : 5712 : paw_dmft%use_sc_dmft = 0
796 : :
797 : 5712 : paw_dmft%dmftbandf = dmftbandf
798 : 5712 : paw_dmft%dmftbandi = dmftbandi
799 : 5712 : paw_dmft%mband = mband
800 : 5712 : paw_dmft%mkmem = mkmem
801 : 5712 : paw_dmft%nkpt = nkpt
802 : 5712 : paw_dmft%nsym = nsym
803 : 5712 : paw_dmft%ntypat = ntypat
804 : :
805 : : ! Spin related variables
806 : 5712 : paw_dmft%nsppol = nsppol
807 : 5712 : paw_dmft%nspinor = nspinor
808 : 5712 : paw_dmft%idmftloop = 0
809 : 5712 : paw_dmft%mbandc = 0
810 : : !paw_dmft%nspden = nspden
811 : :
812 : 5712 : paw_dmft%dmft_read_occnd = dtset%dmft_read_occnd
813 : :
814 [ - + ]: 5712 : 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 [ - + + - : 34272 : ABI_MALLOC(paw_dmft%occnd,(2,mband,mband,nkpt,nsppol*use_dmft))
- + + - -
+ + - + +
+ - + - +
+ - + - +
- + ]
821 [ + + + - : 17136 : ABI_MALLOC(paw_dmft%band_in,(mband*use_dmft))
+ - + + -
+ - + -
+ ]
822 [ + + + - : 17136 : ABI_MALLOC(paw_dmft%include_bands,((dmftbandf-dmftbandi+1)*use_dmft))
+ - + + -
+ - + -
+ ]
823 [ + + + - : 11424 : ABI_MALLOC(paw_dmft%exclude_bands,(mband*use_dmft))
- + - + -
+ ]
824 : : endif
825 : :
826 [ + + ]: 5712 : 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 [ + + + + : 17036 : 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 : 5711 : 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 [ + - + - : 5711 : & .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 [ + - ]: 5711 : ABI_SFREE(paw_dmft%occnd)
2260 [ + - ]: 5711 : ABI_SFREE(paw_dmft%band_in)
2261 [ + - ]: 5711 : ABI_SFREE(paw_dmft%include_bands)
2262 [ + - ]: 5711 : ABI_SFREE(paw_dmft%exclude_bands)
2263 [ + + ]: 5711 : ABI_SFREE(paw_dmft%siz_proj)
2264 [ + + ]: 5711 : ABI_SFREE(paw_dmft%zarot)
2265 [ - + ]: 5711 : ABI_SFREE(paw_dmft%phimtphi)
2266 [ - + ]: 5711 : ABI_SFREE(paw_dmft%phimtphi_int)
2267 [ + + ]: 5711 : ABI_SFREE(paw_dmft%phi_int)
2268 [ - + ]: 5711 : ABI_SFREE(paw_dmft%ylm)
2269 [ - + ]: 5711 : ABI_SFREE(paw_dmft%dpro)
2270 [ - + ]: 5711 : ABI_SFREE(paw_dmft%bessel)
2271 [ - + ]: 5711 : ABI_SFREE(paw_dmft%bessel_int)
2272 [ + + ]: 5711 : ABI_SFREE(paw_dmft%lpawu)
2273 [ + + ]: 5711 : ABI_SFREE(paw_dmft%omega_lo)
2274 [ + + ]: 5711 : ABI_SFREE(paw_dmft%wgt_wlo)
2275 [ + + ]: 5711 : ABI_SFREE(paw_dmft%slm2ylm)
2276 [ + + ]: 5711 : ABI_SFREE(paw_dmft%jmj2ylm)
2277 : :
2278 : 5711 : paw_dmft%nband => null()
2279 : 5711 : paw_dmft%dmft_shiftself => null()
2280 : 5711 : paw_dmft%dmft_nominal => null()
2281 : 5711 : 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 : 5711 : paw_dmft%wtk => null()
2285 : 5711 : paw_dmft%typat => null()
2286 : 5711 : paw_dmft%int_meshsz => null()
2287 : :
2288 [ + + ]: 5711 : 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 : 5711 : call destroy_sc_dmft_paralkgb(paw_dmft)
2293 [ + + ]: 5711 : if (paw_dmft%use_dmft == 1) then
2294 : 43 : call destroy_paral_dmft(paw_dmft%distrib)
2295 : : end if
2296 : :
2297 : 5711 : 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 : 5711 : subroutine destroy_sc_dmft_paralkgb(paw_dmft)
2657 : :
2658 : : !Arguments ------------------------------------
2659 : : type(paw_dmft_type), intent(inout) :: paw_dmft
2660 : : ! *********************************************************************
2661 : :
2662 [ + + ]: 5711 : ABI_SFREE(paw_dmft%bandc_proc)
2663 [ + + ]: 5711 : ABI_SFREE(paw_dmft%use_bandc)
2664 : :
2665 : 5711 : 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 : : !!***
|
