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