Line data Source code
1 : !! NAME
2 : !! m_varpeq
3 : !!
4 : !! FUNCTION
5 : !! Description
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2023-2026 ABINIT group (VV, MG)
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 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_varpeq
22 :
23 : use, intrinsic :: iso_c_binding
24 : use defs_basis
25 : use m_abicore
26 : use m_dtset
27 : use m_dtfil
28 : use m_crystal
29 : use m_ebands
30 : use m_errors
31 : use m_krank
32 : use netcdf
33 : use m_nctk
34 : use m_xmpi
35 : use m_ifc
36 : use m_wfd
37 : use m_sigtk
38 :
39 : use defs_datatypes, only : pseudopotential_type
40 : use m_numeric_tools, only : interpolate_ur
41 : use m_fstrings, only : sjoin, ktoa, ftoa, strcat, ltoa, itoa, yesno
42 : use m_time, only : cwtime_report, cwtime
43 : use m_io_tools, only : file_exists, iomode_from_fname, open_file
44 : use m_pptools, only : write_xsf
45 : use m_geometry, only : xcart2xred, phdispl_cart2red_nmodes
46 : use m_kpts, only : kpts_map, kpts_timrev_from_kptopt, bzlint_t, kptrlatt_from_ngkpt
47 : use m_fft_mesh, only : calc_ceikr
48 : use m_pawtab, only : pawtab_type
49 : use m_gstore, only : gstore_t, gqk_t
50 : use m_supercell, only : supercell_type
51 : use m_paw_sphharm, only : ylm_angular_mesh
52 : use m_fftcore, only : ngfft_seq
53 : use m_ephtk, only : ephtk_get_mpw_gmax, EPHTK_WTOL
54 : use m_dynmat, only : phdispl_from_eigvec
55 : use m_phonons, only : pheigvec_rotate
56 : use m_splines, only : spline_complex, splint_complex
57 : use m_occ, only : occ_be, occ_fd
58 :
59 :
60 : implicit none
61 :
62 : private
63 : !!***
64 :
65 : !----------------------------------------------------------------------
66 :
67 : !!****t* m_varpeq/polstate_t
68 : !! NAME
69 : !! polstate_t
70 : !!
71 : !! FUNCTION
72 : !! Datatype defining MPI-distributed parameters of polaronic states for a
73 : !! given spin index (if collinear magnetism i.e. nsppol 2). Local variables and
74 : !! arrays start with `my_`, global have `*glob*` in their names. The MPI grid is
75 : !! inherited from a gstore%gqk object at initialization.
76 : !!
77 : !! SOURCE
78 :
79 : type, public :: polstate_t
80 :
81 : character(len=abi_slen) :: aseed = " "
82 : ! Specifies the type of initial seed for charge localization A_nk
83 : ! Possible values: "gau_energy", "gau_length", "random", "even", "localize"
84 :
85 : logical :: translate = .false.
86 : ! Flag controlling treatment of polaronic solution invariant by primitive
87 : ! translations inside supercell
88 : ! if .true. and np > 1, the next solution is found to be orthogonal to all
89 : ! previous states + their translated replicas
90 :
91 : integer :: spin = -1
92 : ! Spin index
93 :
94 : integer :: ip = -1
95 : ! Current polaronic state
96 :
97 : integer :: np = -1
98 : ! Number of polaronic states
99 :
100 : integer :: nkbz = -1
101 : ! Number of k-points in full BZ
102 :
103 : integer :: nqbz = -1
104 : ! Number of q-points in full BZ
105 :
106 : integer :: psign = 1
107 : ! Sign of the charge carrier
108 : ! +1 -- electron addition (electron polaron)
109 : ! -1 -- electron removal (hole polaron)
110 :
111 : real(dp) :: e_frohl
112 : ! Long-range divergence correction of polaron binding energy due to g(0) avg
113 :
114 : real(dp) :: efilter
115 : ! Energy filter
116 :
117 : integer :: ngkpt(3)
118 : ! Number of points in the uniform k-grid defining the electronic subspace
119 :
120 : real(dp) :: gpr_energy(2)
121 : ! Gaussian parameters for energy-based initialization
122 :
123 : real(dp) :: gpr_length(3)
124 : ! Gaussian parameters for localization length-based initialization
125 :
126 : logical, allocatable :: has_prev_grad(:)
127 : ! (np)
128 : ! Flag indicating if an electronic gradient has been computed at previous step
129 :
130 : real(dp), allocatable :: gradres(:)
131 : ! (np)
132 : ! L^2-norm of the electronic gradient for each state
133 :
134 : real(dp), allocatable :: enterms(:,:)
135 : ! (4, np)
136 : ! Polaronic energy terms for each state
137 : ! 1 -> electronic term
138 : ! 2 -> vibrational term
139 : ! 3 -> electron-phonon term
140 : ! 4 -> polaron energy level
141 :
142 : real(dp), allocatable :: eig(:,:)
143 : ! (gqk%nb_k, gstore%ebands%nkpt)
144 : ! Electronic eigenstates participating in the polaron formation
145 : ! Correspond either to valence or conduction manifold, not both
146 : ! Band edge is set to 0; valence states are also inverted
147 :
148 : real(dp), allocatable :: my_g0(:)
149 : ! (gqk%my_npert)
150 : ! Long-range correction for electron-phonon matrix elements at q=\Gamma
151 :
152 : real(dp), allocatable :: my_qpts(:,:)
153 : ! (3, gqk%my_nq)
154 : ! q-points treated by this MPI proc
155 :
156 : real(dp), allocatable :: displ(:,:,:)
157 : ! (3, self%scell%natom, np)
158 : ! Atomic displacements \Delta \tau_{\alpha \kappa*p} for each state
159 :
160 : real(dp), pointer :: my_kpts(:,:) => null()
161 : ! (3, gqk%my_nk)
162 : ! k-points treated by this MPI proc, points to gqk%my_kpts(:,:)
163 :
164 : complex(dp), allocatable :: my_a(:,:,:)
165 : ! (gqk%nb_k, gqk%my_nk, np)
166 : ! Electronic coefficients A_nk for each state treated by this MPI proc
167 :
168 : complex(dp), allocatable :: a_glob(:,:,:)
169 : ! (gqk%nb_k, gqk%glob_nk, np)
170 : ! Global array of electronic coefficients A_nk at each state
171 :
172 : complex(dp), allocatable :: my_b(:,:,:)
173 : ! (gqk%my_npert, gqk%my_nq, np)
174 : ! Vibrational coefficients B_q\nu for each state treated by this MPI proc
175 :
176 : complex(dp), allocatable :: my_prev_b(:,:)
177 : ! (gqk%my_npert, gqk%my_nq)
178 : ! Previous vibrational coefficients B_q\nu for each state treated by this MPI proc
179 :
180 : complex(dp), allocatable :: my_prev_b_hop(:,:,:)
181 : ! (gqk%my_npert, gqk%my_nq)
182 : ! Previous vibrational coefficients B_q\nu for each state (hopping optimization)
183 :
184 : complex(dp), allocatable :: my_pc(:,:)
185 : ! (gqk%nb_k, gqk%my_nk)
186 : ! Preconditioner at current state treated by this MPI proc
187 :
188 : complex(dp), allocatable :: my_grad(:,:)
189 : ! (gqk%nb_k, gqk%my_nk)
190 : ! Electronic gradient D_nk at current state treated by this MPI proc
191 : ! orthogonal to states
192 :
193 : complex(dp), allocatable :: my_prev_grad(:,:)
194 : ! (gqk%nb_k, gqk%my_nk)
195 : ! Previous electronic gradient D_nk at current state, orthogonal to states
196 :
197 : complex(dp), allocatable :: my_pcgrad(:,:)
198 : ! (gqk%nb_k, gqk%my_nk)
199 : ! Preconditioned gradient at current state, orthogonal to states
200 :
201 : complex(dp), allocatable :: my_prev_pcgrad(:,:)
202 : ! (gqk%nb_k, gqk%my_nk)
203 : ! Previous preconditioned gradient at current state, orthogonal to states
204 :
205 : complex(dp), allocatable :: my_pcjgrad(:,:)
206 : ! (gqk%nb_k, gqk%my_nk)
207 : ! Preconditioned conjugate gradient at current state treated by this MPI proc
208 : ! orthogonal to current state & normalized
209 :
210 : complex(dp), allocatable :: my_prev_pcjgrad(:,:)
211 : ! (gqk%nb_k, gqk%my_nk)
212 : ! Previous preconditioned conjugate gradient at current state
213 :
214 : complex(dp), allocatable :: pcjgrad_glob(:,:)
215 : ! (gqk%nb_k, gqk%glob_nk)
216 : ! Global preconditioned conjugate gradient at current state
217 : ! orthogonal to current state & normalized
218 :
219 : complex(dp), allocatable :: my_phgrad(:,:,:)
220 : ! (gqk%my_npert, gqk%my_nq, np)
221 : ! Phonon gradient \nabla E_{B_q\nu} for each state treated by this MPI proc
222 :
223 : complex(dp), allocatable :: my_eff_phforce(:,:,:)
224 : ! (gqk%my_npert, gqk%my_nq, np)
225 : ! Effective phonon force F_q\nu at current state treated by this MPI proc
226 :
227 : real(dp), allocatable :: hop_ts(:)
228 : ! (np)
229 : ! Line minimization timestep for hopping optimization for each state
230 :
231 : class(gqk_t), pointer :: gqk => null()
232 : ! Datastructure storing e-ph matrix elements treated by this MPI proc
233 :
234 : type(krank_t) :: krank_kpts
235 : ! Object used to find k-points in BZ
236 :
237 : type(krank_t) :: krank_qpts
238 : ! Object used to find q-points in BZ
239 :
240 : type(crystal_t) :: cryst
241 : ! Object storing information on crystal structure & symmetries
242 :
243 : contains
244 :
245 : procedure :: setup => polstate_setup
246 : ! Set up optimization process by specifying initial electronic vector A_nk
247 :
248 : procedure :: localize => polstate_localize
249 : ! Localize polaron at current state. From A_nk calculate:
250 : ! vibrational coefficients B_q\nu, energy terms, polaron level
251 :
252 : procedure :: get_enel => polstate_get_enel
253 : ! Calculate and return electronic energy term
254 :
255 : procedure :: get_enph => polstate_get_enph
256 : ! Calculate and return vibrational energy term
257 :
258 : procedure :: get_enelph => polstate_get_enelph
259 : ! Calculate and return electron-phonon energy term
260 :
261 : procedure :: get_lm_theta => polstate_get_lm_theta
262 : ! Calculate and return line minimization parameter \theta
263 :
264 : procedure :: calc_grad => polstate_calc_grad
265 : ! Calculate steepest descent vector (electronic gradient)
266 :
267 : procedure :: calc_pcjgrad => polstate_calc_pcjgrad
268 : ! Calculate preconditioned conjugate gradient direction
269 :
270 : procedure :: update_pc => polstate_update_pc
271 : ! Update preconditioner
272 :
273 : procedure :: update_a => polstate_update_a
274 : ! Update array of electronic coefficients
275 :
276 : procedure :: ort_to_states => polstate_ort_to_states
277 : ! Orthogonalize a given vector to a set of polaronic states
278 :
279 : procedure :: calc_b_from_a => polstate_calc_b_from_a
280 : ! Calculate vibrational coefficients B_q\nu from a known set of electronic
281 : ! coefficients A_nk
282 :
283 : procedure :: calc_b_from_displ => polstate_calc_b_from_displ
284 : ! Calculate vibrational coefficients B_q\nu from a known set of displacements
285 : ! \Delta \tau_{\alpha \kappa p}
286 :
287 : procedure :: calc_displ_from_b => polstate_calc_displ_from_b
288 : ! Calculate polaron-induced displacements \Delta \tau_{\alpha \kappa p} from
289 : ! a known set of vibrational coefficients B_q\nu
290 :
291 : procedure :: calc_clb_displ => polstate_calc_clb_displ
292 : ! Calculate displacements induced by the Coulomb forces for a unit charge
293 : ! localized on a supercell atom
294 :
295 : procedure :: calc_phgrad => polstate_calc_phgrad
296 : ! Calculate steepest descent vector (phonon gradient)
297 :
298 : procedure :: calc_hpol => polstate_calc_hpol
299 : ! Construct and calculate the polaron Hamiltonian, <A_bra|H(B)|A_ket>
300 :
301 : procedure :: seed_a => polstate_seed_a
302 : ! Seed an initial vector of electronic coefficients A_nk
303 :
304 : procedure :: load_a => polstate_load_a
305 : ! Initialize A_nk from an existing vector of electronic coefficients
306 :
307 : procedure :: load_b => polstate_load_b
308 : ! Initialize B_q\nu from an existing vector of electronic coefficients
309 :
310 : procedure :: linterp_b => polstate_linterp_b
311 : ! Linear interpolation of B_\qnu coefficients between initial and final state
312 :
313 : procedure :: redistr_b => polstate_redistr_b
314 : ! Redistribute B_\qnu coefficients with spline interpolation
315 :
316 : procedure :: calc_eff_phforce => polstate_calc_eff_phforce
317 : ! Calculate effective phonon force for hopping optimization
318 :
319 : procedure :: calc_hop_timestep => polstate_calc_hop_timestep
320 : ! Calculate timestep from effective phonon force for this state
321 :
322 : procedure :: get_sqnorm => polstate_get_sqnorm
323 : ! Helper function to compute squared L^2-norm of MPI-distributed array
324 :
325 : procedure :: gather => polstate_gather
326 : ! Helper function to gather MPI-distributed array into a global one
327 :
328 : procedure :: filter => polstate_filter
329 : ! Helper function to filter MPI-distributed array
330 :
331 : procedure :: get_krank_glob => polstate_get_krank_glob
332 : ! Helper function to calculate global krank objects
333 :
334 : procedure :: free => polstate_free
335 : ! Free memory
336 :
337 : end type polstate_t
338 : !!***
339 :
340 : !----------------------------------------------------------------------
341 :
342 : !!****t* m_varpeq/varpeq_t
343 : !! NAME
344 : !! varpeq_t
345 : !!
346 : !! FUNCTION
347 : !! Variational Polaron Equations datatype. Stores variational optimization
348 : !! parameters, polaronic states and provides higher-level procedures required
349 : !! for variational optimization and related input and output.
350 : !!
351 : !! SOURCE
352 :
353 : type, public :: varpeq_t
354 :
355 : character(len=abi_slen) :: mode = " "
356 : ! Specifies the calculation mode
357 : ! Possible values: "polaron", "hopping"
358 :
359 : character(len=abi_slen) :: pkind = " "
360 : ! Specifies the kind of polaron
361 : ! Possible values: "hole", "electron"
362 :
363 : character(len=abi_slen) :: aseed = " "
364 : ! Specifies the type of initial seed for charge localization A_nk
365 : ! Possible values: "gau_energy", "gau_length", "random", "even", "localize"
366 :
367 : logical :: fix_displ = .false.
368 : ! Flag indicating if vibrational coefficients need to be fixed between steps.
369 : ! True for hopping calculations.
370 :
371 : logical :: use_filter = .false.
372 : ! Flag indicating if the energy filtering for electronic states was used
373 :
374 : logical :: is_complete = .false.
375 : ! Flag indicating if the datatype is completely or partially initialized
376 : ! Required to distinguish between newly created and loaded-from-disk datatype
377 :
378 : logical :: restart = .false.
379 : ! Flag to check if a restart from a *VPQ.nc file is needed
380 :
381 : logical :: interp = .false.
382 : ! Flag to check if an interpolation from a *VPQ.nc file is needed
383 :
384 : logical :: ld_flag = .false.
385 : ! Flag indicating if internal variables have been loaded from source
386 :
387 : logical :: g0_flag = .true.
388 : ! Flag indicating if the average of el-ph matrix elements at Gamma is computed
389 :
390 : logical :: translate = .false.
391 : ! Flag controlling the translational invariance of polaronic solutions
392 :
393 : integer :: ncid = nctk_noid
394 : ! Netcdf file handle used to save results
395 :
396 : integer :: atloc = -1
397 : ! At which atom the polaron is localized
398 : ! from 1 to cryst%natom if aseed = "localize"
399 :
400 : integer :: hop_nstep = -1
401 : ! Maximum number of iterations for hopping transport optimization
402 :
403 : integer :: nstep = -1
404 : ! Maximum number of iterations for optimization of a single polaronic state
405 :
406 : integer :: nstep_ort = -1
407 : ! Maximum number of iterations for which orthogonalization to all previous
408 : ! states is performed
409 :
410 : integer :: nsppol = -1
411 : ! Number of independent spin polarizations
412 :
413 : integer :: nstates = -1
414 : ! Number of polaronic states to be optimized for each spin polarization
415 :
416 : integer :: natom3 = -1
417 : ! 3*gstore%cryst%natom
418 : ! Number of atomic perturbations, mainly used to dimensionalize arrays
419 :
420 : integer :: max_nk = -1
421 : ! Maximum number of k-points among between independent spin polarizations
422 :
423 : integer :: max_nq = -1
424 : ! Maximum number of q-points among between independent spin polarizations
425 :
426 : integer :: max_nb = -1
427 : ! Maximum number of bands among between independent spin polarizations
428 :
429 : integer :: frohl_ntheta = -1
430 : ! Number of angular mesh division for spherical integration of long-range
431 : ! divergence of electron-phonon matrix elements
432 :
433 : real(dp) :: e_frohl
434 : ! Long-range divergence correction of polaron binding energy
435 :
436 : real(dp) :: mixing_factor
437 : ! Mixing factor to be used in the solver for vibrational coefficients
438 :
439 : real(dp) :: tolgrs
440 : ! L^2 gradient norm tolerance
441 :
442 : real(dp) :: hop_tolgrs
443 : ! L^2 force norm tolerance for hopping optimization
444 :
445 : real(dp) :: hop_ts
446 : ! Timestep for hopping
447 :
448 : integer :: ngkpt(3)
449 : ! Number of points in the uniform k-grid defining the electronic subspace
450 :
451 : integer, allocatable :: nk_spin(:)
452 : ! (nsppol)
453 : ! Number of k-points for each spin polarization
454 :
455 : integer, allocatable :: nq_spin(:)
456 : ! (nsppol)
457 : ! Number of q-points for each spin polarization
458 :
459 : integer, allocatable :: nb_spin(:)
460 : ! (nsppol)
461 : ! Number of bands for each spin polarization
462 :
463 : integer, allocatable :: brange_spin(:,:)
464 : ! (2, nsppol)
465 : ! Number of bands for each spin polarization
466 :
467 : integer, allocatable :: cvflag_spin(:,:,:)
468 : ! (nstates, hop_nstep, nsppol)
469 : ! Convergence flags at each state for each spin:
470 : ! 0 --> calculation is not converged
471 : ! 1 --> calculation is converged
472 :
473 : integer, allocatable :: hop_nstep2cv_spin(:)
474 : ! (nsppol)
475 : ! number of steps to convergence the hopping optimization for each spin
476 :
477 : integer, allocatable :: nstep2cv_spin(:,:,:)
478 : ! (nstates, hop_nstep, nsppol)
479 : ! number of steps to convergence at each state for each spin
480 :
481 : real(dp), allocatable :: erange_spin(:)
482 : ! (nsppol)
483 : ! energy window wrt to vbm/cbm for hole/electron polaron for each spin
484 :
485 : integer, allocatable :: k2ibz_spin(:,:)
486 : ! (max_nk, nsppol)
487 : ! bz->ibz index table for kpoints (related to sell%gstore%kibz)
488 :
489 : integer, allocatable :: q2ibz_spin(:,:)
490 : ! (max_nq, nsppol)
491 : ! bz->ibz index table for qpoints (related to sell%gstore%qibz)
492 :
493 : real(dp), allocatable :: chrgat(:)
494 : ! (natom)
495 : ! charge of the cations
496 :
497 : real(dp), allocatable :: hop_hist_spin(:,:,:,:)
498 : ! (5, nstates, hop_nstep, nsppol)
499 : ! hopping optimization history at each state (image) for each spin
500 : ! 1 - bare phonon force norms, i.e. |dE_pol/dB| at each image
501 : ! 2 - string forces, i.e. |B^n - B^(n-1)|/timestep for step n
502 : ! 3 - diabatic energy curve for initial state: <A_in|B_i|A_in>
503 : ! 4 - diabatic energy curve for final state: <A_end|B_i|A_end>
504 : ! 5 - overlap integral <A_in|B_i|A_end>
505 :
506 : real(dp), allocatable :: scf_hist_spin(:,:,:,:,:)
507 : ! (6, nstep, nstates, hop_nstep, nsppol)
508 : ! scf optimization history at each state for each spin
509 : ! 1 - binding energy, 2 - electronic energy, 3 - phonon energy
510 : ! 4 - electron-phonon energy, 5 - localization energy
511 : ! 6 - gradient norm
512 :
513 : real(dp), allocatable :: kpts_spin(:,:,:)
514 : ! (3, max_nk, nsppol)
515 : ! k-points for each spin
516 :
517 : real(dp), allocatable :: qpts_spin(:,:,:)
518 : ! (3, max_nq, nsppol)
519 : ! q-points for each spin
520 :
521 : complex(dp), allocatable :: a_spin(:,:,:,:)
522 : ! (max_nb, max_nk, nstates, nsppol)
523 : ! optimized electronic coefficients a_nk at each state for each spin
524 :
525 : complex(dp), allocatable :: b_spin(:,:,:,:)
526 : ! (natom3, max_nq, nstates, nsppol)
527 : ! optimized vibrational coefficients b_q\nu at each state for each spin
528 :
529 : class(gstore_t), pointer :: gstore => null()
530 : ! object storing el-ph matrix elements and other related quantities
531 :
532 : type(crystal_t) :: cryst
533 : ! Object storing information on crystal structure & symmetries
534 :
535 : type(gaps_t) :: gaps
536 : ! Object used to get information on bandgap
537 :
538 : type(supercell_type) :: scell
539 : ! Object used to get information on supercell
540 :
541 : type(polstate_t), allocatable :: polstate(:)
542 : ! (nsppol)
543 : ! Datatype providing data and lower-level methods for polaronic states
544 : ! at each spin polarization
545 :
546 : contains
547 :
548 : procedure :: init => varpeq_init
549 : ! Initialize object
550 :
551 : procedure :: load => varpeq_load
552 : ! Load the initial electronic vector from a *VPQ.nc netcdf file
553 :
554 : procedure :: solve => varpeq_solve
555 : ! Solve variational polaron equations at each polaronic state for each spin
556 :
557 : procedure :: record => varpeq_record
558 : ! Record the current SCF optimization data
559 :
560 : procedure :: collect => varpeq_collect
561 : ! Collect SCF optimization results from each spin
562 :
563 : procedure :: print_scf_results => varpeq_print_scf_results
564 : ! Output SCF optimization final results
565 :
566 : procedure :: print_hop_results => varpeq_print_hop_results
567 : ! Output hopping optimization final results
568 :
569 : procedure :: print_metadata => varpeq_print_metadata
570 : ! Output parameters defining varpeq calculation
571 :
572 : procedure :: ncwrite => varpeq_ncwrite
573 : ! Save results to a *VPQ.nc netcdf file
574 :
575 : procedure :: ncread => varpeq_ncread
576 : ! Initialize an incomplete object from a *VPQ.nc netcdf file
577 :
578 : procedure :: compare => varpeq_compare
579 : ! Compare basic dimensions with another instance of varpeq_t datatype
580 :
581 : procedure :: calc_fravg => varpeq_calc_fravg
582 : ! Calculate average Fr\"ohlich (long-range) contribution to the polaron
583 : ! binding energy & average of electron-phonon matrix elements at q=\Gamma
584 :
585 : procedure :: hop_setup => varpeq_hop_setup
586 : ! Set up the data structure for hopping transport calculations
587 :
588 : procedure :: hop => varpeq_hop
589 : ! Polaron hopping optimization for self%nstates images
590 :
591 : procedure :: free => varpeq_free
592 : ! Free memory
593 :
594 : end type varpeq_t
595 : !!***
596 :
597 : public :: varpeq_run
598 : ! Main entry point
599 :
600 : public :: varpeq_plot
601 : ! Compute polaron wavefunctions and atomic displacements in the supercell and write results to XSF files
602 :
603 : contains !=====================================================================
604 :
605 :
606 : !!****f* m_varpeq/varpeq_run
607 : !! NAME
608 : !! varpeq_run
609 : !!
610 : !! FUNCTION
611 : !! Higher-level subroutine that solves variational polaron equations, produces
612 : !! necessary output and writes results to a *VPQ.nc file.
613 : !!
614 : !! INPUTS
615 : !! gstore<gstore_t>=Electron-phonon matrix elements and related quantities.
616 : !! dtset<dataset_type>=All input variables for this dataset.
617 : !! dtfil<datafiles_types>=Variables related to files.
618 : !!
619 : !! OUTPUT
620 : !!
621 : !! SOURCE
622 :
623 9 : subroutine varpeq_run(gstore, dtset, dtfil)
624 :
625 : !Arguments ------------------------------------
626 : class(gstore_t), intent(in) :: gstore
627 : type(dataset_type), intent(in) :: dtset
628 : type(datafiles_type), intent(in) :: dtfil
629 : !arrays
630 : integer :: units(2)
631 :
632 : !Local variables-------------------------------
633 468 : type(varpeq_t) :: vpq
634 : !----------------------------------------------------------------------
635 :
636 27 : units = [std_out, ab_out]
637 9 : call wrtout(units, sjoin(ch10, " === Variational Polaron Equations ==="))
638 :
639 9 : call vpq%init(gstore, dtset)
640 9 : if (vpq%frohl_ntheta > 0) call vpq%calc_fravg(avg_g0=vpq%g0_flag)
641 9 : if (vpq%interp .or. vpq%restart) call vpq%load(dtfil, dtset%vpq_select)
642 :
643 9 : call vpq%print_metadata()
644 :
645 8 : select case(vpq%mode)
646 : case ("polaron")
647 8 : call vpq%solve(verbose=.true.)
648 :
649 : case ("hopping")
650 1 : ABI_CHECK(vpq%nstates > 1, "varpeq_run: vpq_nstates must be > 1 for hopping")
651 1 : call vpq%hop_setup(dtset)
652 1 : call vpq%hop()
653 :
654 : case default
655 9 : ABI_ERROR(sjoin("varpeq_run, unsupported varpeq mode: ", vpq%mode))
656 : end select
657 :
658 9 : call vpq%collect()
659 :
660 8 : select case(vpq%mode)
661 : case ("polaron")
662 8 : call vpq%print_scf_results(ihop=1)
663 : case ("hopping")
664 9 : call vpq%print_hop_results()
665 : end select
666 :
667 9 : call vpq%ncwrite(dtset, dtfil)
668 9 : call vpq%free()
669 :
670 9 : end subroutine varpeq_run
671 : !!***
672 :
673 : !----------------------------------------------------------------------
674 :
675 : !!****f* m_varpeq/varpeq_free
676 : !! NAME
677 : !! varpeq_free
678 : !!
679 : !! FUNCTION
680 : !! Free dynamic memory
681 : !!
682 : !! SOURCE
683 :
684 13 : subroutine varpeq_free(self)
685 :
686 : !Arguments ------------------------------------
687 : class(varpeq_t), target, intent(inout) :: self
688 :
689 : !Local variables-------------------------------
690 : integer :: my_is
691 : !----------------------------------------------------------------------
692 :
693 : ! Free allocatable arrays
694 : ! integer
695 13 : ABI_SFREE(self%nk_spin)
696 13 : ABI_SFREE(self%nq_spin)
697 13 : ABI_SFREE(self%nb_spin)
698 13 : ABI_SFREE(self%brange_spin)
699 13 : ABI_SFREE(self%cvflag_spin)
700 13 : ABI_SFREE(self%hop_nstep2cv_spin)
701 13 : ABI_SFREE(self%nstep2cv_spin)
702 : ! real
703 13 : ABI_SFREE(self%chrgat)
704 13 : ABI_SFREE(self%erange_spin)
705 13 : ABI_SFREE(self%hop_hist_spin)
706 13 : ABI_SFREE(self%scf_hist_spin)
707 13 : ABI_SFREE(self%k2ibz_spin)
708 13 : ABI_SFREE(self%q2ibz_spin)
709 13 : ABI_SFREE(self%kpts_spin)
710 13 : ABI_SFREE(self%qpts_spin)
711 :
712 : ! complex
713 13 : ABI_SFREE(self%a_spin)
714 13 : ABI_SFREE(self%b_spin)
715 :
716 : ! Free local datatypes
717 13 : call self%cryst%free()
718 :
719 : ! Close netcdf file
720 13 : if (self%ncid /= nctk_noid) then
721 9 : NCF_CHECK(nf90_close(self%ncid))
722 9 : self%ncid = nctk_noid
723 : end if
724 :
725 : ! If entry is completely initialized (e.g. from self%init call), free remaining
726 : ! datatypes and nullify pointers
727 13 : if (self%is_complete) then
728 9 : call self%gaps%free()
729 18 : do my_is=1,self%gstore%my_nspins
730 18 : call self%polstate(my_is)%free()
731 : enddo
732 9 : call self%scell%free()
733 18 : ABI_SFREE(self%polstate)
734 9 : self%gstore => null()
735 : endif
736 :
737 13 : end subroutine varpeq_free
738 : !!***
739 :
740 : !----------------------------------------------------------------------
741 :
742 : !!****f* m_varpeq/varpeq_compare
743 : !! NAME
744 : !! varpeq_compare
745 : !!
746 : !! FUNCTION
747 : !! Compares the current instance of varpeq datatype with another one (possibly
748 : !! incomplete).
749 : !!
750 : !! INPUTS
751 : !! other<varpeq_t>=Varpeq datatype to compare with.
752 : !! bz_mismatch [optional]=if .true. mismatch between BZ sampling is allowed
753 : !! (required for comparison prior to an interpolation)
754 : !!
755 : !! OUTPUT
756 : !!
757 : !! SOURCE
758 :
759 4 : subroutine varpeq_compare(self, other, bz_mismatch)
760 :
761 : !Arguments ------------------------------------
762 : class(varpeq_t), target, intent(in) :: self, other
763 : logical, optional, intent(in) :: bz_mismatch
764 :
765 : !Local variables-------------------------------
766 : integer :: ierr
767 : real(dp) :: cpu, wall, gflops
768 : !----------------------------------------------------------------------
769 :
770 : ! TODO: provide more flexibility for nstates
771 :
772 2 : call cwtime(cpu, wall, gflops, "start")
773 :
774 2 : ierr = 0
775 :
776 : ! Compare basic dimensions
777 2 : call check_(self%pkind == other%pkind, "Difference found in pkind.")
778 2 : call check_(self%nsppol == other%nsppol, "Difference found in nsppol.")
779 : !call check_(self%nstates == other%nstates, "Difference found in nstates.")
780 : call check_(self%cryst%compare(other%cryst) == 0, &
781 2 : "Difference found in cryst.")
782 : call check_(all(self%brange_spin == other%brange_spin), &
783 8 : "Difference found in brange_spin.")
784 :
785 : ! If bz_mismatch is not allowed, also compare k/q-grids
786 2 : if (present(bz_mismatch)) then
787 2 : if (.not. bz_mismatch) then
788 : call check_(all(abs(self%kpts_spin - other%kpts_spin) < tol6), &
789 1004 : "Difference found in kpts_spin.")
790 : call check_(all(abs(self%qpts_spin - other%qpts_spin) < tol6), &
791 1004 : "Difference found in nq_spin.")
792 : end if
793 : endif
794 :
795 2 : ABI_CHECK(ierr == 0, "Error in varpeq_compare, see previous messages!")
796 :
797 2 : call cwtime_report(" varpeq_compare", cpu, wall, gflops)
798 :
799 : contains
800 12 : subroutine check_(cond, msg)
801 : logical, intent(in) :: cond
802 : character(len=*), intent(in) :: msg
803 12 : ABI_CHECK_NOSTOP(cond, msg, ierr)
804 12 : end subroutine check_
805 :
806 : end subroutine varpeq_compare
807 : !!***
808 :
809 : !----------------------------------------------------------------------
810 :
811 : !!****f* m_varpeq/varpeq_ncread
812 : !! NAME
813 : !! varpeq_ncread
814 : !!
815 : !! FUNCTION
816 : !! Reads basic dimensions of varpeq_t datatype from a *VPQ.nc netcdf file.
817 : !!
818 : !! INPUTS
819 : !! path=Path a *VPQ.nc file to be read.
820 : !! comm=MPI communicator.
821 : !! keep_open [optional]=if .true. keep the nc file handle open for further
822 : !! reading. Default: .false.
823 : !!
824 : !! OUTPUT
825 : !!
826 : !! SOURCE
827 :
828 4 : subroutine varpeq_ncread(self, path, comm, keep_open)
829 :
830 : !Arguments ------------------------------------
831 : class(varpeq_t), target, intent(inout) :: self
832 : character(len=fnlen), intent(in) :: path
833 : integer, intent(in) :: comm
834 : logical, optional, intent(in) :: keep_open
835 :
836 : !Local variables-------------------------------
837 : integer :: ncid, nsppol, natom3, nstates
838 : real(dp) :: cpu, wall, gflops
839 4 : real(dp), contiguous, pointer :: rpt_d5(:,:,:,:,:)
840 : !----------------------------------------------------------------------
841 :
842 4 : call cwtime(cpu, wall, gflops, "start")
843 :
844 4 : ABI_CHECK(file_exists(path), sjoin(" varpeq_ncread: cannot find *VPQ.nc file", path))
845 :
846 4 : NCF_CHECK(nctk_open_read(ncid, path, comm))
847 :
848 : ! Read crystal structure
849 4 : call self%cryst%ncread(ncid)
850 :
851 : ! Read varpeq dimensions
852 4 : NCF_CHECK(nctk_get_dim(ncid, "nsppol", self%nsppol))
853 4 : NCF_CHECK(nctk_get_dim(ncid, "nstates", self%nstates))
854 4 : NCF_CHECK(nctk_get_dim(ncid, "natom3", self%natom3))
855 4 : NCF_CHECK(nctk_get_dim(ncid, "max_nk", self%max_nk))
856 4 : NCF_CHECK(nctk_get_dim(ncid, "max_nq", self%max_nq))
857 4 : NCF_CHECK(nctk_get_dim(ncid, "max_nb", self%max_nb))
858 4 : nsppol = self%nsppol
859 4 : nstates = self%nstates
860 4 : natom3 = self%natom3
861 :
862 : ! Read data
863 : ! Static arrays
864 4 : NCF_CHECK(nf90_get_var(ncid, vid("vpq_pkind"), self%pkind))
865 4 : NCF_CHECK(nf90_get_var(ncid, vid("ngkpt"), self%ngkpt))
866 :
867 : ! Allocatable arrays
868 12 : ABI_MALLOC(self%nk_spin, (nsppol))
869 8 : ABI_MALLOC(self%nq_spin, (nsppol))
870 8 : ABI_MALLOC(self%nb_spin, (nsppol))
871 12 : ABI_MALLOC(self%brange_spin, (2, nsppol))
872 16 : ABI_MALLOC(self%kpts_spin, (3, self%max_nk, nsppol))
873 16 : ABI_MALLOC(self%qpts_spin, (3, self%max_nq, nsppol))
874 24 : ABI_MALLOC(self%a_spin, (self%max_nb, self%max_nk, nstates, nsppol))
875 24 : ABI_MALLOC(self%b_spin, (natom3, self%max_nq, nstates, nsppol))
876 :
877 : ! integer
878 4 : NCF_CHECK(nf90_get_var(ncid, vid("nk_spin"), self%nk_spin))
879 4 : NCF_CHECK(nf90_get_var(ncid, vid("nq_spin"), self%nq_spin))
880 4 : NCF_CHECK(nf90_get_var(ncid, vid("nb_spin"), self%nb_spin))
881 4 : NCF_CHECK(nf90_get_var(ncid, vid("brange_spin"), self%brange_spin))
882 : ! real
883 4 : NCF_CHECK(nf90_get_var(ncid, vid("kpts_spin"), self%kpts_spin))
884 4 : NCF_CHECK(nf90_get_var(ncid, vid("qpts_spin"), self%qpts_spin))
885 : ! complex
886 : ! TODO: is it possible to encapsulate this trick as some abstraction?
887 : call c_f_pointer(c_loc(self%a_spin), rpt_d5, &
888 24 : [2, self%max_nb, self%max_nk, nstates, nsppol])
889 4 : NCF_CHECK(nf90_get_var(ncid, vid("a_spin"), rpt_d5))
890 :
891 : call c_f_pointer(c_loc(self%b_spin), rpt_d5, &
892 24 : [2, natom3, self%max_nq, nstates, nsppol])
893 4 : NCF_CHECK(nf90_get_var(ncid, vid("b_spin"), rpt_d5))
894 :
895 4 : if (present(keep_open)) then
896 2 : if (keep_open) self%ncid = ncid
897 : else
898 2 : NCF_CHECK(nf90_close(ncid))
899 2 : self%ncid = nctk_noid
900 : end if
901 :
902 4 : self%is_complete = .false.
903 :
904 8 : call cwtime_report(" varpeq_ncread", cpu, wall, gflops)
905 :
906 : !----------------------------------------------------------------------
907 :
908 : contains
909 40 : integer function vid(var_name)
910 : character(len=*),intent(in) :: var_name
911 40 : vid = nctk_idname(ncid, var_name)
912 : end function vid
913 :
914 : end subroutine varpeq_ncread
915 : !!***
916 :
917 : !----------------------------------------------------------------------
918 :
919 : !!****f* m_varpeq/varpeq_ncwrite
920 : !! NAME
921 : !! varpeq_ncwrite
922 : !!
923 : !! FUNCTION
924 : !! Dump varpeq variables in a newly created *VPQ.nc netcdf file.
925 : !!
926 : !! INPUTS
927 : !! dtset<dataset_type>=All input variables for this dataset.
928 : !! dtfil<datafiles_types>=Variables related to files.
929 : !!
930 : !! OUTPUT
931 : !!
932 : !! SOURCE
933 :
934 27 : subroutine varpeq_ncwrite(self, dtset, dtfil)
935 :
936 : !Arguments ------------------------------------
937 : class(varpeq_t), target, intent(inout) :: self
938 : type(dataset_type), intent(in) :: dtset
939 : type(datafiles_type),intent(in) :: dtfil
940 :
941 : !Local variables-------------------------------
942 : character(len=fnlen) :: path
943 : integer, parameter :: master = 0
944 : integer :: my_rank, ncid, ncerr
945 : real(dp) :: cpu, wall, gflops
946 9 : real(dp), contiguous, pointer :: rpt_d5(:,:,:,:,:)
947 : integer :: units(2)
948 : !----------------------------------------------------------------------
949 :
950 27 : units = [std_out, ab_out]
951 :
952 9 : call cwtime(cpu, wall, gflops, "start")
953 :
954 9 : my_rank = xmpi_comm_rank(self%gstore%comm)
955 :
956 : ! Create netcdf file (only master works, HDF5 + MPI-IO can be handled after
957 : ! by reopening the file inside ncwrite_comm)
958 9 : path = strcat(dtfil%filnam_ds(4), "_VPQ.nc")
959 9 : if (my_rank == master) then
960 :
961 9 : call wrtout(units, sjoin(ch10, sjoin("Saving results to:", path)))
962 :
963 : ! Master creates the netcdf file used to store the data.
964 9 : NCF_CHECK(nctk_open_create(self%ncid, path, xmpi_comm_self))
965 9 : ncid = self%ncid
966 :
967 : ! Write the crystal (TR & invsersion symmetry only) & ebands dataset_type
968 9 : NCF_CHECK(self%cryst%ncwrite(ncid))
969 9 : NCF_CHECK(self%gstore%ebands%ncwrite(ncid))
970 :
971 : ! Add varpeq dimensions.
972 : ncerr = nctk_def_dims(ncid, [ &
973 : nctkdim_t("hop_nstep", self%hop_nstep), &
974 : nctkdim_t("nstep", self%nstep), nctkdim_t("nsppol", self%nsppol), &
975 : nctkdim_t("nstates", self%nstates), nctkdim_t("natom3", self%natom3), &
976 : nctkdim_t("max_nk", self%max_nk), nctkdim_t("max_nq", self%max_nq), &
977 : nctkdim_t("max_nb", self%max_nb), nctkdim_t("nkibz", self%gstore%nkibz), &
978 : nctkdim_t("nqibz", self%gstore%nqibz)], &
979 99 : defmode=.true.)
980 9 : NCF_CHECK(ncerr)
981 :
982 : ! Define scalars
983 : ! integers
984 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
985 : "eph_task", "nkbz", "nqbz", "frohl_ntheta", "vpq_avg_g", "vpq_translate", &
986 : "vpq_interp", "vpq_nstates", "vpq_nstep_ort", "vpq_select", "vpq_mesh_fact", &
987 117 : "vpq_hop_nstep"])
988 9 : NCF_CHECK(ncerr)
989 : ! real
990 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: &
991 36 : "vpq_tolgrs", "e_frohl", "vpq_mix_fact"])
992 9 : NCF_CHECK(ncerr)
993 :
994 : ! Define arrays with results
995 : ncerr = nctk_def_arrays(ncid, [ &
996 : nctkarr_t("vpq_mode", "c", "character_string_length"), &
997 : nctkarr_t("vpq_pkind", "c", "character_string_length"), &
998 : nctkarr_t("vpq_aseed", "c", "character_string_length"), &
999 : nctkarr_t("ngkpt", "int", "three"), &
1000 : nctkarr_t("gstore_ngqpt", "int", "three"), &
1001 : nctkarr_t("nk_spin", "int", "nsppol"), &
1002 : nctkarr_t("nq_spin", "int", "nsppol"), &
1003 : nctkarr_t("nb_spin", "int", "nsppol"), &
1004 : nctkarr_t("brange_spin", "int", "two, nsppol"), &
1005 : nctkarr_t("cvflag_spin", "int", "nstates, hop_nstep, nsppol"), &
1006 : nctkarr_t("hop_nstep2cv_spin", "int", "nsppol"), &
1007 : nctkarr_t("nstep2cv_spin", "int", "nstates, hop_nstep, nsppol"), &
1008 : nctkarr_t("vpq_trvec", "int", "three"), &
1009 : nctkarr_t("k2ibz_spin", "int", "max_nk, nsppol"), &
1010 : nctkarr_t("q2ibz_spin", "int", "max_nq, nsppol"), &
1011 : nctkarr_t("erange_spin", "dp", "nsppol"), &
1012 : nctkarr_t("hop_hist_spin", "dp", "five, nstates, hop_nstep, nsppol"), &
1013 : nctkarr_t("scf_hist_spin", "dp", "six, nstep, nstates, hop_nstep, nsppol"), &
1014 : nctkarr_t("kibz", "dp", "three, nkibz"), &
1015 : nctkarr_t("qibz", "dp", "three, nqibz"), &
1016 : nctkarr_t("kpts_spin", "dp", "three, max_nk, nsppol"), &
1017 : nctkarr_t("qpts_spin", "dp", "three, max_nq, nsppol"), &
1018 : nctkarr_t("cb_min_spin", "dp", "nsppol"), &
1019 : nctkarr_t("vb_max_spin", "dp", "nsppol"), &
1020 : nctkarr_t("vpq_gpr_energy", "dp", "two"), &
1021 : nctkarr_t("vpq_gpr_length", "dp", "three"), &
1022 : nctkarr_t("a_spin", "dp", "two, max_nb, max_nk, nstates, nsppol"), &
1023 : nctkarr_t("b_spin", "dp", "two, natom3, max_nq, nstates, nsppol") &
1024 261 : ])
1025 9 : NCF_CHECK(ncerr)
1026 :
1027 : ! Write data
1028 9 : NCF_CHECK(nctk_set_datamode(ncid))
1029 : ! Scalars
1030 : ! integer
1031 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
1032 : "eph_task", "nkbz", "nqbz", "frohl_ntheta", "vpq_avg_g", "vpq_translate", &
1033 : "vpq_interp", "vpq_nstates", "vpq_nstep_ort", "vpq_select", "vpq_mesh_fact", &
1034 : "vpq_hop_nstep"], &
1035 : [dtset%eph_task, self%gstore%nkbz, self%gstore%nqbz, self%frohl_ntheta, &
1036 : dtset%vpq_avg_g, dtset%vpq_translate, dtset%vpq_interp, dtset%vpq_nstates, &
1037 225 : dtset%vpq_nstep_ort, dtset%vpq_select, dtset%vpq_mesh_fact, dtset%vpq_hop_nstep])
1038 9 : NCF_CHECK(ncerr)
1039 : ! real
1040 : ncerr = nctk_write_dpscalars(ncid, [character(len=nctk_slen) :: &
1041 : "vpq_tolgrs", "e_frohl", "vpq_mix_fact"], &
1042 63 : [self%tolgrs, self%e_frohl, dtset%vpq_mix_fact])
1043 9 : NCF_CHECK(ncerr)
1044 :
1045 : ! Arrays
1046 : ! character
1047 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_mode"), self%mode))
1048 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_pkind"), self%pkind))
1049 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_aseed"), self%aseed))
1050 : ! integer
1051 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngkpt"), self%ngkpt))
1052 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "gstore_ngqpt"), self%gstore%ngqpt))
1053 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nk_spin"), self%nk_spin))
1054 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nq_spin"), self%nq_spin))
1055 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nb_spin"), self%nb_spin))
1056 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "brange_spin"), self%brange_spin))
1057 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cvflag_spin"), self%cvflag_spin))
1058 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "hop_nstep2cv_spin"), self%hop_nstep2cv_spin))
1059 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "nstep2cv_spin"), self%nstep2cv_spin))
1060 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_trvec"), dtset%vpq_trvec))
1061 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "k2ibz_spin"), self%k2ibz_spin))
1062 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "q2ibz_spin"), self%q2ibz_spin))
1063 : ! real
1064 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "hop_hist_spin"), self%hop_hist_spin))
1065 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "scf_hist_spin"), self%scf_hist_spin))
1066 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kibz"), self%gstore%kibz))
1067 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qibz"), self%gstore%qibz))
1068 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "kpts_spin"), self%kpts_spin))
1069 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpts_spin"), self%qpts_spin))
1070 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "cb_min_spin"), self%gaps%cb_min))
1071 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vb_max_spin"), self%gaps%vb_max))
1072 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_gpr_energy"), dtset%vpq_gpr_energy))
1073 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "vpq_gpr_length"), dtset%vpq_gpr_length))
1074 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "erange_spin"), self%erange_spin))
1075 : ! complex
1076 : call c_f_pointer(c_loc(self%a_spin), rpt_d5, &
1077 54 : [2, self%max_nb, self%max_nk, self%nstates, self%nsppol])
1078 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "a_spin"), rpt_d5))
1079 :
1080 : call c_f_pointer(c_loc(self%b_spin), rpt_d5, &
1081 54 : [2, self%natom3, self%max_nq, self%nstates, self%nsppol])
1082 9 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "b_spin"), rpt_d5))
1083 : end if ! master
1084 :
1085 9 : call xmpi_barrier(self%gstore%comm)
1086 9 : call cwtime_report(" varpeq: ncwrite", cpu, wall, gflops)
1087 :
1088 9 : end subroutine varpeq_ncwrite
1089 : !!***
1090 :
1091 : !!----------------------------------------------------------------------
1092 :
1093 : !!****f* m_varpeq/varpeq_print_metadata
1094 : !! NAME
1095 : !! varpeq_print_metadata
1096 : !!
1097 : !! FUNCTION
1098 : !! Output parameters defining varpeq calculation
1099 : !!
1100 : !! INPUTS
1101 : !!
1102 : !! OUTPUT
1103 : !!
1104 : !! SOURCE
1105 :
1106 : !subroutine varpeq_print_metadata(self, dtset)
1107 9 : subroutine varpeq_print_metadata(self)
1108 :
1109 : !Arguments ------------------------------------
1110 : class(varpeq_t), target, intent(inout) :: self
1111 : !type(dataset_type), intent(in) :: dtset
1112 :
1113 : !Local variables-------------------------------
1114 : !scalars
1115 : character(len=abi_slen) :: msg
1116 : !character(len=abi_slen) :: use_frohl = "no"
1117 : integer :: spin, my_rank
1118 : integer, parameter :: master = 0
1119 : !arrays
1120 : integer :: units(2)
1121 :
1122 : !----------------------------------------------------------------------
1123 :
1124 9 : my_rank = xmpi_comm_rank(self%gstore%comm)
1125 :
1126 27 : units = [std_out, ab_out]
1127 9 : if (my_rank == master) then
1128 :
1129 9 : call title_("Polaron:")
1130 :
1131 9 : call entry_("Polaron kind", self%pkind)
1132 9 : write(msg, '(i0,a,i0,a,i0)') self%ngkpt(1), "x", self%ngkpt(2), "x", self%ngkpt(3)
1133 9 : call entry_("BvK supercell", msg)
1134 9 : call entry_("Number of independent spin polarizations", itoa(self%nsppol))
1135 9 : call entry_("Number of polaronic states", itoa(self%nstates))
1136 9 : call entry_("Filtering of electronic states", yesno(self%use_filter))
1137 9 : if (self%use_filter) then
1138 2 : do spin=1,self%nsppol
1139 1 : write(msg, '(a,i0,a,i0)') "Energy filter wrt band edge for spin ", &
1140 2 : spin, "/", self%nsppol
1141 2 : call entry_(msg, sjoin(ftoa(self%erange_spin(spin)*Ha_eV, "es8.2"), "eV"))
1142 : enddo
1143 : endif
1144 :
1145 9 : call title_("Long-range corrections:")
1146 9 : call entry_("Frohlich correction", yesno(self%frohl_ntheta > 0))
1147 9 : if (self%frohl_ntheta > 0) then
1148 9 : call entry_("Frohlich correction value", sjoin(ftoa(self%e_frohl*Ha_eV), "eV"))
1149 9 : call entry_("Frohlich correction included in matrix elements", yesno(self%g0_flag))
1150 : endif
1151 :
1152 9 : call title_("Optimization parameters:")
1153 :
1154 : ! TODO change this once a support for starting from B is added
1155 9 : call entry_("Initial seed", "charge localization A_nk")
1156 9 : call entry_("Initial seed type", strseed())
1157 :
1158 9 : call entry_("Tolerance on the gradient norm", ftoa(self%tolgrs, "es8.2"))
1159 9 : call entry_("Maximum number of iterations per state", itoa(self%nstep))
1160 9 : if (self%nstates > 1) then
1161 2 : call wrtout(units, " (for pstate > 1)")
1162 2 : call entry_("Number of orthogonaliztion steps", itoa(self%nstep_ort))
1163 2 : call entry_("Translational invariance", yesno(self%translate))
1164 : endif
1165 :
1166 : endif
1167 :
1168 : contains
1169 27 : subroutine title_(str)
1170 : character(len=*), intent(in) :: str
1171 27 : call wrtout(units, sjoin(ch10, str))
1172 27 : end subroutine title_
1173 :
1174 113 : subroutine entry_(name, val)
1175 : character(len=*), intent(in) :: name
1176 : character(len=*), intent(in) :: val
1177 : character(len=abi_slen) :: output
1178 113 : write(output, '(a5,a,a2,a)') "* ", trim(name), ": ", trim(val)
1179 113 : call wrtout(units, output)
1180 113 : end subroutine entry_
1181 :
1182 9 : character(len=abi_slen) function strseed()
1183 18 : select case(self%aseed)
1184 : case ("gau_energy")
1185 9 : write(strseed, *) "Gaussian, based on the electroic energies"
1186 : case ("gau_length")
1187 0 : write(strseed, *) "Gaussian, based on the localizaiton length"
1188 : case ("random")
1189 0 : write(strseed, *) "random"
1190 : case ("even")
1191 0 : write(strseed, *) "even"
1192 : case ("localize")
1193 0 : write(strseed, *) "manual loclization"
1194 : case default
1195 9 : write(strseed, *) "undefined"
1196 : end select
1197 :
1198 9 : if (self%interp .or. self%restart) then
1199 0 : write(strseed, *) "loaded from file"
1200 : endif
1201 :
1202 9 : strseed = adjustl(strseed)
1203 9 : end function strseed
1204 :
1205 : end subroutine varpeq_print_metadata
1206 : !!***
1207 :
1208 : !!----------------------------------------------------------------------
1209 :
1210 : !!****f* m_varpeq/varpeq_print_scf_results
1211 : !! NAME
1212 : !! varpeq_print_scf_results
1213 : !!
1214 : !! FUNCTION
1215 : !! Output SCF optimization results
1216 : !!
1217 : !! INPUTS
1218 : !! ineb [optional]=Hopping iteration. Relevant only if self%vpq_mode="hopping".
1219 : !!
1220 : !! OUTPUT
1221 : !!
1222 : !! SOURCE
1223 :
1224 8 : subroutine varpeq_print_scf_results(self, ihop)
1225 :
1226 : !Arguments ------------------------------------
1227 : class(varpeq_t), target, intent(inout) :: self
1228 : integer, optional, intent(in) :: ihop
1229 :
1230 : !Local variables-------------------------------
1231 : !scalars
1232 : character(len=5000) :: msg
1233 : integer, parameter :: master = 0
1234 : integer :: my_rank, spin, ip, ii, ihop_
1235 : !arrays
1236 : integer :: units(2)
1237 : !----------------------------------------------------------------------
1238 :
1239 8 : my_rank = xmpi_comm_rank(self%gstore%comm)
1240 :
1241 8 : ihop_ = 1
1242 8 : if (present(ihop)) ihop_ = ihop
1243 :
1244 24 : units = [std_out, ab_out]
1245 8 : if (my_rank == master) then
1246 16 : do spin=1,self%gstore%nsppol
1247 26 : do ip=1,self%nstates
1248 10 : call header_(spin, ip)
1249 :
1250 334 : do ii=1,self%nstep2cv_spin(ip, ihop_, spin)
1251 326 : call report_(spin, ip, ii, ihop_, self%nstep2cv_spin(ip, ihop_, spin))
1252 :
1253 : enddo
1254 : enddo
1255 : enddo
1256 : endif
1257 :
1258 : contains
1259 10 : subroutine header_(spin, state)
1260 : integer, intent(in) :: spin, state
1261 : character(len=5000) :: sep
1262 :
1263 10 : call wrtout(units, "")
1264 10 : if (state == 1) then
1265 8 : call wrtout(units, " Printing the optimization logs")
1266 : endif
1267 :
1268 10 : write(sep, "(a3,a)") "", repeat('-', 86)
1269 : write(msg, '(a5,a,i0,a,i0,a,i0,a,i0)') &
1270 10 : "* ", "spin ", spin, "/", self%nsppol, ", pstate ", state, "/", self%nstates
1271 :
1272 10 : call wrtout(units, sep)
1273 10 : call wrtout(units, msg)
1274 10 : write(msg, '(a3,a)') "", "* values in the optimization log are in (a.u.)"
1275 10 : call wrtout(units, msg)
1276 :
1277 10 : if ((state > 1) .and. (self%nstep_ort > 1)) then
1278 2 : write(msg, "(a3,a,i0)") "", "(o) - orthogonal all pstates < ", state
1279 2 : call wrtout(units, msg)
1280 : endif
1281 :
1282 10 : call wrtout(units, sep)
1283 : !write(msg, '(a3,a4,7a13,a5)') "", "Step", "E_pol", "E_el", "E_ph", "E_elph", "E2_elph", &
1284 : ! "epsilon", "||grad||", ""
1285 10 : write(msg, '(a3,a4,6a13,a5)') "", "Step", "E_pol", "E_el", "E_ph", "E_elph", &
1286 20 : "epsilon", "||grad||", ""
1287 10 : call wrtout(units, msg)
1288 10 : end subroutine header_
1289 :
1290 316 : subroutine report_(spin, state, step, hop_step, step2conv)
1291 : integer, intent(in) :: spin, state, step, hop_step, step2conv
1292 : character(len=5000) :: sep
1293 : character(len=abi_slen) :: ort_flag
1294 : logical :: is_conv
1295 : real(dp) :: enpol, enel, enph, enelph, eps, grs
1296 :
1297 316 : enpol = self%scf_hist_spin(1, step, state, hop_step, spin)
1298 316 : enel = self%scf_hist_spin(2, step, state, hop_step, spin)
1299 316 : enph = self%scf_hist_spin(3, step, state, hop_step, spin)
1300 316 : enelph = self%scf_hist_spin(4, step, state, hop_step, spin)
1301 316 : eps = self%scf_hist_spin(5, step, state, hop_step, spin)
1302 316 : grs = self%scf_hist_spin(6, step, state, hop_step, spin)
1303 :
1304 316 : write(sep, "(a3,a)") "", repeat('-', 86)
1305 :
1306 : ! for states > 1 add flag if orthogonalization is used
1307 316 : write(ort_flag, '(a)') ""
1308 316 : if ((ip > 1) .and. (step <= self%nstep_ort)) then
1309 60 : write(ort_flag, '(a)') "(o)"
1310 : endif
1311 :
1312 : write(msg,'(a3,i4,6es13.4,a5)') &
1313 316 : "", step, enpol, enel, enph, enelph, eps, grs, trim(ort_flag)
1314 316 : call wrtout(units, msg)
1315 :
1316 316 : if (step == step2conv) then
1317 :
1318 10 : is_conv = (self%cvflag_spin(state, hop_step, spin) == 1)
1319 10 : if (is_conv) then
1320 9 : write(msg, '(a3,a,es11.4,a,es11.4)') "", "Converged: ||grad||=", grs, &
1321 18 : " < vpq_tolgrs=", self%tolgrs
1322 : else
1323 1 : write(msg, '(a3,a,es11.4,a,es11.4)') "", "Unconverged: ||grad||=", grs, &
1324 2 : " > vpq_tolgrs=", self%tolgrs
1325 : endif
1326 :
1327 10 : call wrtout(units, sep)
1328 10 : call wrtout(units, msg)
1329 :
1330 10 : write(msg, '(a3,a,es17.8)') "", "E_pol (eV):", enpol*Ha_eV
1331 10 : call wrtout(units, msg)
1332 10 : write(msg, '(a3,a,es17.8)') "", " eps (eV):", eps*Ha_eV
1333 10 : call wrtout(units, msg)
1334 :
1335 10 : call wrtout(units, sep)
1336 :
1337 : endif
1338 316 : end subroutine report_
1339 :
1340 : end subroutine varpeq_print_scf_results
1341 : !!***
1342 :
1343 : !!----------------------------------------------------------------------
1344 :
1345 : !!****f* m_varpeq/varpeq_print_hop_results
1346 : !! NAME
1347 : !! varpeq_print_hop_results
1348 : !!
1349 : !! FUNCTION
1350 : !! Output hopping optimization results
1351 : !!
1352 : !! INPUTS
1353 : !!
1354 : !! OUTPUT
1355 : !!
1356 : !! SOURCE
1357 :
1358 1 : subroutine varpeq_print_hop_results(self)
1359 :
1360 : !Arguments ------------------------------------
1361 : class(varpeq_t), target, intent(inout) :: self
1362 :
1363 : !Local variables-------------------------------
1364 : !scalars
1365 : character(len=5000) :: msg
1366 : integer, parameter :: master = 0
1367 : integer :: my_rank, spin, ip, ihop
1368 : integer :: hop_nstep2cv, scf_nstep2cv
1369 : !arrays
1370 : integer :: units(2)
1371 :
1372 : !----------------------------------------------------------------------
1373 :
1374 1 : my_rank = xmpi_comm_rank(self%gstore%comm)
1375 :
1376 :
1377 3 : units = [std_out, ab_out]
1378 :
1379 1 : if (my_rank == master) then
1380 2 : do spin=1,self%gstore%nsppol
1381 1 : hop_nstep2cv = self%hop_nstep2cv_spin(spin)
1382 :
1383 1 : call hop_header_(spin)
1384 16 : do ihop=1,self%hop_nstep2cv_spin(spin)
1385 16 : call hop_report_(spin, ihop)
1386 : enddo
1387 :
1388 1 : call header_(spin)
1389 11 : do ip=1,self%nstates
1390 9 : scf_nstep2cv = self%nstep2cv_spin(ip, hop_nstep2cv, spin)
1391 10 : call report_(spin, ip, scf_nstep2cv, hop_nstep2cv)
1392 : enddo
1393 :
1394 : enddo
1395 : endif
1396 :
1397 : contains
1398 1 : subroutine header_(spin)
1399 : integer, intent(in) :: spin
1400 : character(len=5000) :: sep
1401 :
1402 1 : call wrtout(units, "")
1403 1 : call wrtout(units, " Printing the minimal energy path")
1404 :
1405 1 : write(sep, "(a3,a)") "", repeat('-', 86)
1406 : write(msg, '(a5,a,i0,a,i0)') &
1407 1 : "* ", "spin ", spin, "/", self%nsppol
1408 :
1409 1 : call wrtout(units, sep)
1410 1 : call wrtout(units, msg)
1411 1 : write(msg, '(a3,a)') "", "* values are in (a.u.)"
1412 1 : call wrtout(units, msg)
1413 :
1414 1 : call wrtout(units, sep)
1415 1 : write(msg, '(a3,a4,6a13)') "", "Image", "E_pol", "E_el", "E_ph", "E_elph", &
1416 2 : "epsilon", "||el_grad||"
1417 1 : call wrtout(units, msg)
1418 1 : end subroutine header_
1419 :
1420 1 : subroutine hop_header_(spin)
1421 : integer, intent(in) :: spin
1422 : character(len=5000) :: sep
1423 :
1424 1 : call wrtout(units, "")
1425 1 : call wrtout(units, " Printing the hopping optimization log")
1426 :
1427 1 : write(sep, "(a3,a)") "", repeat('-', 43)
1428 : write(msg, '(a5,a,i0,a,i0)') &
1429 1 : "* ", "spin ", spin, "/", self%nsppol
1430 :
1431 1 : call wrtout(units, sep)
1432 1 : call wrtout(units, msg)
1433 1 : write(msg, '(a3,a)') "", "* values are in (a.u.)"
1434 1 : call wrtout(units, msg)
1435 :
1436 1 : call wrtout(units, sep)
1437 1 : write(msg, '(a3,a4,a17,a18)') "", "Step", &
1438 2 : "max||ph_grad||", "max||hop_grad||"
1439 1 : call wrtout(units, msg)
1440 1 : end subroutine hop_header_
1441 :
1442 9 : subroutine report_(spin, state, step, hop_step)
1443 : integer, intent(in) :: spin, state, step, hop_step
1444 : character(len=5000) :: sep
1445 : real(dp) :: enpol, enel, enph, enelph, eps, grs
1446 : !real(dp) :: min_epol, max_epol, ehop
1447 :
1448 9 : enpol = self%scf_hist_spin(1, step, state, hop_step, spin)
1449 9 : enel = self%scf_hist_spin(2, step, state, hop_step, spin)
1450 9 : enph = self%scf_hist_spin(3, step, state, hop_step, spin)
1451 9 : enelph = self%scf_hist_spin(4, step, state, hop_step, spin)
1452 9 : eps = self%scf_hist_spin(5, step, state, hop_step, spin)
1453 9 : grs = self%scf_hist_spin(6, step, state, hop_step, spin)
1454 :
1455 9 : write(sep, "(a3,a)") "", repeat('-', 86)
1456 :
1457 : write(msg,'(a3,i4,6es13.4)') &
1458 9 : "", state, enpol, enel, enph, enelph, eps, grs
1459 9 : call wrtout(units, msg)
1460 :
1461 9 : if (state == self%nstates) then
1462 1 : call wrtout(units, sep)
1463 : endif
1464 :
1465 9 : end subroutine report_
1466 :
1467 15 : subroutine hop_report_(spin, hop_step)
1468 : integer, intent(in) :: spin, hop_step
1469 : character(len=5000) :: sep
1470 : real(dp) :: ph_max_grs, hop_max_grs
1471 :
1472 165 : ph_max_grs = maxval(self%hop_hist_spin(1, :, hop_step, spin))
1473 165 : hop_max_grs = maxval(self%hop_hist_spin(2, :, hop_step, spin))
1474 :
1475 15 : write(sep, "(a3,a)") "", repeat('-', 43)
1476 :
1477 : write(msg,'(a3,i4,es17.4,es18.4)') &
1478 15 : "", hop_step, ph_max_grs, hop_max_grs
1479 15 : call wrtout(units, msg)
1480 :
1481 15 : if (hop_step == self%hop_nstep2cv_spin(spin)) then
1482 1 : call wrtout(units, sep)
1483 : endif
1484 :
1485 15 : end subroutine hop_report_
1486 :
1487 : end subroutine varpeq_print_hop_results
1488 : !!***
1489 :
1490 : !!----------------------------------------------------------------------
1491 :
1492 : !!****f* m_varpeq/varpeq_collect
1493 : !! NAME
1494 : !! varpeq_collect
1495 : !!
1496 : !! FUNCTION
1497 : !! Collect SCF optimization results for each spin
1498 : !!
1499 : !! INPUTS
1500 : !!
1501 : !! OUTPUT
1502 : !!
1503 : !! SOURCE
1504 :
1505 198 : subroutine varpeq_collect(self)
1506 :
1507 : !Arguments ------------------------------------
1508 : class(varpeq_t), target, intent(inout) :: self
1509 :
1510 : !Local variables-------------------------------
1511 : class(gqk_t), pointer :: gqk
1512 : class(polstate_t), pointer :: polstate
1513 : integer :: ierr
1514 : integer :: my_is, spin, my_ik, ik_glob, ik_ibz, ip
1515 : integer :: my_iq, iq_glob, iq_ibz, my_pert, pert_glob
1516 : integer :: oc_scf, oc_a, oc_b, oc_k, oc_q
1517 : !----------------------------------------------------------------------
1518 :
1519 : ! Gather the SCF process evolution data
1520 18 : call xmpi_sum(self%cvflag_spin, self%gstore%comm, ierr)
1521 18 : call xmpi_sum(self%scf_hist_spin, self%gstore%comm, ierr)
1522 18 : call xmpi_sum(self%nstep2cv_spin, self%gstore%comm, ierr)
1523 :
1524 18 : call xmpi_sum(self%hop_hist_spin, self%gstore%comm, ierr)
1525 18 : call xmpi_sum(self%hop_nstep2cv_spin, self%gstore%comm, ierr)
1526 :
1527 : ! Gather electron/phonon vectors and k/q points
1528 13534 : self%a_spin(:,:,:,:) = zero
1529 31364 : self%b_spin(:,:,:,:) = zero
1530 1382 : self%k2ibz_spin(:,:) = zero
1531 2398 : self%q2ibz_spin(:,:) = zero
1532 9484 : self%qpts_spin(:,:,:) = zero
1533 5420 : self%kpts_spin(:,:,:) = zero
1534 36 : do my_is=1,self%gstore%my_nspins
1535 18 : spin = self%gstore%my_spins(my_is)
1536 18 : gqk => self%gstore%gqk(my_is)
1537 18 : polstate => self%polstate(spin)
1538 :
1539 56 : do ip=1,self%nstates
1540 :
1541 : ! electronic vector
1542 3492 : do my_ik=1,gqk%my_nk
1543 3454 : ik_glob = gqk%my_kstart + my_ik - 1
1544 13498 : self%a_spin(:, ik_glob, ip, spin) = polstate%my_a(:, my_ik, ip)
1545 : enddo
1546 :
1547 : ! phonon vector
1548 4526 : do my_iq=1,gqk%my_nq
1549 4470 : iq_glob = gqk%my_qstart + my_iq - 1
1550 31328 : do my_pert=1,gqk%my_npert
1551 26820 : pert_glob = gqk%my_pert_start + my_pert - 1
1552 168054 : self%b_spin(pert_glob, iq_glob, :, spin) = polstate%my_b(my_pert, my_iq, :)
1553 : enddo
1554 : enddo
1555 :
1556 : enddo
1557 :
1558 : ! k-points
1559 1364 : do my_ik=1,gqk%my_nk
1560 1346 : ik_ibz = gqk%my_k2ibz(1, my_ik)
1561 1346 : ik_glob = gqk%my_kstart + my_ik - 1
1562 5384 : self%kpts_spin(:, ik_glob, spin) = polstate%my_kpts(:, my_ik)
1563 1364 : self%k2ibz_spin(ik_glob, spin) = ik_ibz
1564 : enddo
1565 :
1566 : ! q-points
1567 2398 : do my_iq=1,gqk%my_nq
1568 2362 : iq_ibz = gqk%my_q2ibz(1, my_iq)
1569 2362 : iq_glob = gqk%my_qstart + my_iq - 1
1570 9448 : self%qpts_spin(:, iq_glob, spin) = polstate%my_qpts(:, my_iq)
1571 2380 : self%q2ibz_spin(iq_glob, spin) = iq_ibz
1572 : enddo
1573 :
1574 : enddo
1575 :
1576 18 : call xmpi_sum(self%a_spin, self%gstore%comm, ierr)
1577 18 : call xmpi_sum(self%b_spin, self%gstore%comm, ierr)
1578 18 : call xmpi_sum(self%k2ibz_spin, self%gstore%comm, ierr)
1579 18 : call xmpi_sum(self%q2ibz_spin, self%gstore%comm, ierr)
1580 18 : call xmpi_sum(self%kpts_spin, self%gstore%comm, ierr)
1581 18 : call xmpi_sum(self%qpts_spin, self%gstore%comm, ierr)
1582 :
1583 : ! Hack to mimic the summation over a non-existing spin communicator
1584 : ! Divide by the number of times we overcount, as we use the global communicator
1585 36 : do my_is=1,self%gstore%my_nspins
1586 18 : spin = self%gstore%my_spins(my_is)
1587 18 : gqk => self%gstore%gqk(my_is)
1588 18 : polstate => self%polstate(spin)
1589 :
1590 18 : oc_scf = gqk%comm%nproc
1591 18 : oc_a = gqk%qpt_pert_comm%nproc
1592 18 : oc_b = gqk%kpt_comm%nproc
1593 18 : oc_k = oc_a
1594 18 : oc_q = oc_b * gqk%pert_comm%nproc
1595 :
1596 354 : self%cvflag_spin(:,:,spin) = self%cvflag_spin(:,:,spin) / oc_scf
1597 197264 : self%scf_hist_spin(:,:,:,:,spin) = self%scf_hist_spin(:,:,:,:,spin) / oc_scf
1598 354 : self%nstep2cv_spin(:,:,spin) = self%nstep2cv_spin(:,:,spin) / oc_scf
1599 1804 : self%hop_hist_spin(:,:,:,spin) = self%hop_hist_spin(:,:,:,spin) / oc_scf
1600 18 : self%hop_nstep2cv_spin(spin) = self%hop_nstep2cv_spin(spin) / oc_scf
1601 :
1602 13516 : self%a_spin(:,:,:,spin) = self%a_spin(:,:,:,spin) / oc_a
1603 31346 : self%b_spin(:,:,:,spin) = self%b_spin(:,:,:,spin) / oc_b
1604 1364 : self%k2ibz_spin(:,spin) = self%k2ibz_spin(:,spin) / oc_k
1605 2380 : self%q2ibz_spin(:,spin) = self%q2ibz_spin(:,spin) / oc_q
1606 5402 : self%kpts_spin(:,:,spin) = self%kpts_spin(:,:,spin) / oc_k
1607 9484 : self%qpts_spin(:,:,spin) = self%qpts_spin(:,:,spin) / oc_q
1608 : enddo
1609 :
1610 18 : end subroutine varpeq_collect
1611 : !!***
1612 :
1613 : !----------------------------------------------------------------------
1614 :
1615 : !!****f* m_varpeq/varpeq_load
1616 : !! NAME
1617 : !! varpeq_load
1618 : !!
1619 : !! FUNCTION
1620 : !! Load and (optionally) interpolate the initial electronic vector A_nk from
1621 : !! a *VPQ.nc netcdf file. Store result in the self%a_spin variable.
1622 : !!
1623 : !! INPUTS
1624 : !! dtfil<datafiles_types>=Variables related to files.
1625 : !! pselect=Which state to select for reload/interpolation. Non-positive value
1626 : !! selects all states.
1627 : !!
1628 : !! OUTPUT
1629 : !!
1630 : !! SOURCE
1631 :
1632 0 : subroutine varpeq_load(self, dtfil, pselect)
1633 :
1634 : !Arguments ------------------------------------
1635 : class(varpeq_t), target, intent(inout) :: self
1636 : type(datafiles_type), intent(in) :: dtfil
1637 : integer,intent(in) :: pselect
1638 :
1639 : !Local variables-------------------------------
1640 : !scalars
1641 0 : type(varpeq_t) :: vpq_ld
1642 0 : type(bzlint_t) :: bzlint
1643 : logical :: single_state
1644 : integer, parameter :: master = 0
1645 : integer :: my_rank, comm, ierr
1646 : integer :: spin, ip, nk, nb, ik, ib
1647 : real(dp) :: cpu, wall, gflops
1648 : !arrays
1649 : integer :: units(2)
1650 : real(dp) :: kpt(3)
1651 0 : real(dp), allocatable :: ak(:), kpts_ld(:,:)
1652 0 : real(dp), contiguous, pointer :: rpt_d2(:,:)
1653 0 : complex(dp), allocatable, target :: a_ld(:,:)
1654 : !----------------------------------------------------------------------
1655 :
1656 0 : call cwtime(cpu, wall, gflops, "start")
1657 :
1658 0 : units = [std_out, ab_out]
1659 0 : comm = self%gstore%comm; my_rank = xmpi_comm_rank(comm)
1660 :
1661 : ! Read A_nk from file. Only the master processor reads, then broadcasts the data
1662 0 : if (my_rank == master) then
1663 0 : call vpq_ld%ncread(dtfil%filvpqin, xmpi_comm_self, keep_open=.false.)
1664 :
1665 : ! Consitency check
1666 0 : call self%compare(vpq_ld, bz_mismatch=self%interp)
1667 0 : if (pselect > 0) then
1668 0 : ABI_CHECK(self%nstates == 1, "vpq_pstates must be 1 if vpq_select > 0.")
1669 0 : ABI_CHECK(pselect <= vpq_ld%nstates, "vpq_select must be <= loaded nstates.")
1670 : single_state = .true.
1671 : else
1672 0 : ABI_CHECK(self%nstates == vpq_ld%nstates, "Diefference found in nstates.")
1673 : single_state = .false.
1674 : endif
1675 :
1676 0 : self%a_spin(:,:,:,:) = zero
1677 :
1678 0 : if (self%interp) then ! Interpolation
1679 0 : call wrtout(units, " - interpolating previous A_nk")
1680 0 : do spin=1,self%nsppol
1681 : ! Setting basic dimensions & arrays
1682 0 : nk = vpq_ld%nk_spin(spin); nb = vpq_ld%nb_spin(spin)
1683 0 : ABI_MALLOC(kpts_ld, (3, nk))
1684 0 : ABI_MALLOC(a_ld, (nb, nk))
1685 0 : ABI_MALLOC(ak, (2*nb))
1686 0 : kpts_ld(:,:) = vpq_ld%kpts_spin(:, 1:nk, spin)
1687 :
1688 : ! Here, interpolation is performed
1689 0 : do ip=1,self%nstates
1690 :
1691 0 : if (single_state) then
1692 0 : a_ld(:,:) = vpq_ld%a_spin(:, 1:nk, pselect, spin)
1693 : else
1694 0 : a_ld(:,:) = vpq_ld%a_spin(:, 1:nk, ip, spin)
1695 : endif
1696 :
1697 0 : call c_f_pointer(c_loc(a_ld), rpt_d2, [2*nb, nk])
1698 :
1699 0 : call bzlint%init(vpq_ld%ngkpt, 2*nb, nk, kpts_ld, rpt_d2)
1700 :
1701 0 : do ik=1,self%nk_spin(spin)
1702 0 : kpt(:) = self%kpts_spin(:, ik, spin)
1703 0 : call bzlint%interp(kpt, ak)
1704 :
1705 0 : do ib=1,nb
1706 0 : self%a_spin(ib, ik, ip, spin) = ak(2*ib-1) + j_dpc*ak(2*ib)
1707 : enddo
1708 : enddo
1709 0 : call bzlint%free()
1710 :
1711 : enddo
1712 0 : ABI_FREE(kpts_ld)
1713 0 : ABI_FREE(a_ld)
1714 0 : ABI_FREE(ak)
1715 : enddo
1716 :
1717 : else ! Restart
1718 0 : call wrtout(units, " - restarting from previous A_nk")
1719 0 : if (single_state) then
1720 0 : self%a_spin(:,:,:,:) = vpq_ld%a_spin(:,:,pselect:pselect,:)
1721 : else
1722 0 : self%a_spin(:,:,:,:) = vpq_ld%a_spin(:,:,:,:)
1723 : endif
1724 : endif
1725 :
1726 0 : call vpq_ld%free()
1727 : endif
1728 :
1729 0 : call xmpi_bcast(self%a_spin, master, comm, ierr)
1730 0 : self%ld_flag = .true.
1731 :
1732 0 : call cwtime_report(" varpeq: load", cpu, wall, gflops)
1733 :
1734 0 : end subroutine varpeq_load
1735 : !!***
1736 :
1737 : !----------------------------------------------------------------------
1738 :
1739 : !!****f* m_varpeq/varpeq_hop
1740 : !! NAME
1741 : !! varpeq_hop
1742 : !!
1743 : !! FUNCTION
1744 : !! Solve the Variational Polaron Equations for polaronic hopping between
1745 : !! initial and final state.
1746 : !! self%nstates act as number of polarnoc images, including the initial
1747 : !! and final one.
1748 : !! Optimization is performed via string method:
1749 : !! [Weinan, J. Chem. Phys. 126, 164103 (2007)]
1750 : !!
1751 : !! INPUTS
1752 : !!
1753 : !! OUTPUT
1754 : !!
1755 : !! SOURCE
1756 :
1757 1 : subroutine varpeq_hop(self)
1758 :
1759 : !Arguments ------------------------------------
1760 : class(varpeq_t), target, intent(inout) :: self
1761 :
1762 : !Local variables-------------------------------
1763 : !scalars
1764 : class(polstate_t), pointer :: polstate
1765 : character(len=5000) :: msg
1766 : integer :: ierr
1767 : integer :: my_is, spin, ip, ihop
1768 : real(dp) :: cpu, wall, gflops
1769 : real(dp) :: ts, dist2
1770 : real(dp) :: phforce_grs, phforce_hop_grs
1771 : !arrays
1772 2 : logical :: is_conv_spin(self%nsppol)
1773 : integer :: units(2)
1774 2 : real(dp) :: b_mesh(self%nstates)
1775 2 : real(dp) :: even_mesh(self%nstates)
1776 :
1777 : !----------------------------------------------------------------------
1778 :
1779 3 : units = [std_out, ab_out]
1780 :
1781 : call wrtout(units, &
1782 1 : sjoin(ch10, "Solving for the polaron hopping..."))
1783 :
1784 1 : call cwtime(cpu, wall, gflops, "start")
1785 :
1786 10 : do ip=1,self%nstates
1787 10 : even_mesh(ip) = (ip - one) / (self%nstates - one)
1788 : enddo
1789 :
1790 2 : is_conv_spin(:) = .false.
1791 16 : do ihop=1,self%hop_nstep
1792 15 : write(msg, '(a5,a,i0,a,i0)') "* ", "Hopping: step ", ihop, "/", self%hop_nstep
1793 15 : call wrtout(units, msg)
1794 :
1795 : ! Solve the variational polaron equations
1796 15 : call self%solve(verbose=.false., ihop=ihop)
1797 :
1798 30 : do my_is=1,self%gstore%my_nspins
1799 15 : spin = self%gstore%my_spins(my_is)
1800 15 : polstate => self%polstate(my_is)
1801 :
1802 7899 : if (ihop == 1) polstate%my_prev_b_hop(:,:,:) = zero
1803 15 : self%hop_nstep2cv_spin(spin) = ihop
1804 :
1805 : ! Calculate the forces
1806 150 : do ip=1,self%nstates
1807 135 : call polstate%calc_eff_phforce("fw_euler", ip)
1808 :
1809 : ! Record the bare phonon force norm
1810 135 : phforce_grs = sqrt(polstate%get_sqnorm("phgrad", ip))
1811 150 : self%hop_hist_spin(1, ip, ihop, spin) = phforce_grs
1812 :
1813 : enddo
1814 :
1815 : ! Calculate optimal timestep
1816 15 : if (ihop == 1) then
1817 8 : do ip=2,self%nstates-1
1818 8 : call polstate%calc_hop_timestep(ip)
1819 : enddo
1820 9 : ts = minval(abs(polstate%hop_ts(2:self%nstates-1)))
1821 :
1822 1 : if (self%hop_ts > zero) ts = self%hop_ts
1823 :
1824 1 : write(msg, '(a5,a,es8.2)') "", "line minimiation step t = ", ts
1825 1 : call wrtout(units, msg)
1826 : endif
1827 :
1828 : ! Evolve B_\qnu coefficients
1829 : polstate%my_b(:,:,:) = &
1830 118275 : polstate%my_b(:,:,:) + ts * polstate%my_eff_phforce(:,:,:)
1831 :
1832 : ! Reparametrize the string:
1833 : ! Get current parametrization
1834 15 : b_mesh(1) = zero
1835 135 : do ip=2,self%nstates
1836 105120 : dist2 = sum(abs(polstate%my_b(:,:,ip) - polstate%my_b(:,:,ip-1))**2)
1837 120 : call xmpi_sum(dist2, polstate%gqk%qpt_pert_comm%value, ierr)
1838 135 : b_mesh(ip) = b_mesh(ip-1) + sqrt(dist2)
1839 : enddo
1840 150 : b_mesh(:) = b_mesh(:) / b_mesh(self%nstates)
1841 :
1842 : ! Redistribute B_q\nu
1843 15 : call polstate%redistr_b(b_mesh, even_mesh, self%nstates)
1844 :
1845 : ! Record phonon force norm, computed from optimization
1846 : ! F_q\nu = 1/ts * |B_q\nu^(n) - B_q\nu^(n-1)|
1847 150 : do ip=1,self%nstates
1848 : phforce_hop_grs = &
1849 118260 : sum(abs(polstate%my_b(:,:,ip) - polstate%my_prev_b_hop(:,:,ip))**2)
1850 135 : call xmpi_sum(phforce_hop_grs, polstate%gqk%qpt_pert_comm%value, ierr)
1851 :
1852 150 : self%hop_hist_spin(2, ip, ihop, spin) = sqrt(phforce_hop_grs) / ts
1853 : enddo
1854 :
1855 165 : if (maxval(self%hop_hist_spin(2, :, ihop, spin)) < self%hop_tolgrs) then
1856 0 : is_conv_spin(spin) = .true.
1857 : endif
1858 :
1859 118290 : polstate%my_prev_b_hop(:,:,:) = polstate%my_b(:,:,:)
1860 : enddo
1861 :
1862 15 : if (all(is_conv_spin) .or. (ihop == self%hop_nstep)) then
1863 10 : do ip=1,self%nstates
1864 : ! Compute the overlaps
1865 : ! <A_in|B_i|A_in>
1866 : self%hop_hist_spin(3, ip, ihop, spin) = &
1867 : polstate%calc_hpol(polstate%my_a(:,:,1), polstate%a_glob(:,:,1), &
1868 9 : polstate%my_a(:,:,1), polstate%my_b(:,:,ip))
1869 : ! <A_end|B_i|A_end>
1870 : self%hop_hist_spin(4, ip, ihop, spin) = &
1871 : polstate%calc_hpol(polstate%my_a(:,:,self%nstates), polstate%a_glob(:,:,self%nstates), &
1872 9 : polstate%my_a(:,:,self%nstates), polstate%my_b(:,:,ip))
1873 : ! <A_in|B_i|A_end>
1874 : self%hop_hist_spin(5, ip, ihop, spin) = &
1875 : polstate%calc_hpol(polstate%my_a(:,:,1), polstate%a_glob(:,:,1), &
1876 24 : polstate%my_a(:,:,self%nstates), polstate%my_b(:,:,ip))
1877 : enddo
1878 : endif
1879 :
1880 16 : if (all(is_conv_spin)) exit
1881 :
1882 : enddo
1883 :
1884 1 : call cwtime_report(" varpeq: hop", cpu, wall, gflops)
1885 :
1886 1 : end subroutine varpeq_hop
1887 : !!***
1888 :
1889 : !----------------------------------------------------------------------
1890 :
1891 : !!****f* m_varpeq/varpeq_solve
1892 : !! NAME
1893 : !! varpeq_solve
1894 : !!
1895 : !! FUNCTION
1896 : !! Solve the Variational Polaron Equations for each spin for self%nstates
1897 : !! polaronic states.
1898 : !!
1899 : !! INPUTS
1900 : !! ineb [optional]=Hopping iteration. Relevant only if self%vpq_mode="hopping".
1901 : !!
1902 : !! OUTPUT
1903 : !!
1904 : !! SOURCE
1905 :
1906 69 : subroutine varpeq_solve(self, verbose, ihop)
1907 :
1908 : !Arguments ------------------------------------
1909 : class(varpeq_t), target, intent(inout) :: self
1910 : logical, intent(in) :: verbose
1911 : integer, optional, intent(in) :: ihop
1912 :
1913 : !Local variables-------------------------------
1914 : class(polstate_t), pointer :: polstate
1915 : character(len=5000) :: msg
1916 : integer :: my_is, spin, ip, ii, ihop_
1917 : real(dp) :: grad_sqnorm
1918 : real(dp) :: cpu, wall, gflops
1919 : integer :: units(2)
1920 : !----------------------------------------------------------------------
1921 :
1922 69 : units = [std_out, ab_out]
1923 :
1924 23 : ihop_ = 1
1925 23 : if (present(ihop)) ihop_ = ihop
1926 :
1927 23 : if (verbose) then
1928 : call wrtout(units, &
1929 8 : sjoin(ch10, "Solving the variational polaron equations for each state..."))
1930 : endif
1931 :
1932 23 : call cwtime(cpu, wall, gflops, "start")
1933 :
1934 46 : do my_is=1,self%gstore%my_nspins
1935 23 : spin = self%gstore%my_spins(my_is)
1936 23 : polstate => self%polstate(my_is)
1937 :
1938 191 : do ip=1,self%nstates
1939 :
1940 145 : if (verbose) then
1941 10 : write(msg, '(a5,a,i0,a,i0,a,i0,a,i0,a)') "* ", "spin ", spin, "/", &
1942 20 : self%nsppol, ", pstate ", ip, "/", self%nstates, "..."
1943 10 : call wrtout(units, msg)
1944 : endif
1945 :
1946 : ! initialize A_nk at this state, orthogonalize to the previous ones
1947 : ! and normalize
1948 145 : if (ihop_ == 1) then
1949 : call polstate%setup(ip, a_src=self%a_spin(:,:,ip,spin), load=self%ld_flag, &
1950 : atloc=self%atloc, chrgat=self%chrgat, scell=self%scell, cryst=self%cryst, &
1951 19 : nstep_loc=self%nstep)
1952 : endif
1953 :
1954 2286 : do ii=1,self%nstep
1955 : ! gather A, get B_qnu, get energies
1956 2285 : call polstate%localize(ip, self%mixing_factor, fix_displ=self%fix_displ)
1957 :
1958 : ! get bare gradient
1959 2285 : call polstate%calc_grad(ip)
1960 :
1961 : ! calculate and save the L^2 gradient norm
1962 2285 : grad_sqnorm = polstate%get_sqnorm("grad", ip)
1963 2285 : polstate%gradres(ip) = sqrt(grad_sqnorm)
1964 :
1965 : ! record the energies & gradient norm to varepq datatype
1966 2285 : call self%record(ii, ip, ihop_, my_is)
1967 :
1968 : ! check if gradient norm is lower than convergence threshold
1969 2285 : if (polstate%gradres(ip) < self%tolgrs) then
1970 144 : self%cvflag_spin(ip, ihop_, spin) = 1
1971 144 : exit
1972 : endif
1973 :
1974 : ! calculate the preconditioner
1975 2141 : call polstate%update_pc(ip)
1976 :
1977 : ! get preconditioned conjugate gradient direction
1978 2141 : call polstate%calc_pcjgrad(ip, ii, self%nstep_ort)
1979 :
1980 : ! update a based on line minimization and pcj direction
1981 2142 : call polstate%update_a(ip)
1982 :
1983 : enddo
1984 :
1985 168 : if (verbose) then
1986 10 : call wrtout(units, " Done")
1987 : endif
1988 :
1989 : enddo
1990 : enddo
1991 :
1992 23 : call cwtime_report(" varpeq: solve", cpu, wall, gflops)
1993 :
1994 23 : end subroutine varpeq_solve
1995 : !!***
1996 :
1997 : !----------------------------------------------------------------------
1998 :
1999 : !!****f* m_varpeq/varpeq_record
2000 : !! NAME
2001 : !! varpeq_record
2002 : !!
2003 : !! FUNCTION
2004 : !! Record variational polaron equations results for current iteration.
2005 : !! Used to propagate the SCF cycle history from self%polstate to the
2006 : !! varpeq datatype itself.
2007 : !!
2008 : !! INPUTS
2009 : !! iter=Current iteration.
2010 : !! ihop=Current iteration (hopping).
2011 : !! ip=Index of a polaronic state.
2012 : !! my_is=Spin polarization treated by this MPI proc.
2013 : !!
2014 : !! OUTPUT
2015 : !!
2016 : !! SOURCE
2017 :
2018 2285 : subroutine varpeq_record(self, iter, ip, ihop, my_is)
2019 :
2020 : !Arguments ------------------------------------
2021 : class(varpeq_t), target, intent(inout) :: self
2022 : integer, intent(in) :: iter, ip, ihop, my_is
2023 :
2024 : !Local variables-------------------------------
2025 : class(polstate_t), pointer :: polstate
2026 : integer :: spin, psign
2027 : real(dp) :: enel, enph, enelph, eps
2028 : !----------------------------------------------------------------------
2029 :
2030 2285 : spin = self%gstore%my_spins(my_is)
2031 2285 : polstate => self%polstate(my_is)
2032 :
2033 2285 : psign = 1
2034 2285 : if (self%pkind == "hole") psign = -1
2035 :
2036 2285 : enel = polstate%enterms(1, ip); enph = polstate%enterms(2, ip)
2037 2285 : enelph = polstate%enterms(3, ip); eps = polstate%enterms(4, ip)
2038 :
2039 2285 : self%scf_hist_spin(1, iter, ip, ihop, spin) = (enel + enph + enelph)
2040 2285 : self%scf_hist_spin(2, iter, ip, ihop, spin) = enel
2041 2285 : self%scf_hist_spin(3, iter, ip, ihop, spin) = enph
2042 2285 : self%scf_hist_spin(4, iter, ip, ihop, spin) = enelph
2043 2285 : self%scf_hist_spin(5, iter, ip, ihop, spin) = psign*eps
2044 2285 : self%scf_hist_spin(6, iter, ip, ihop, spin) = polstate%gradres(ip)
2045 2285 : self%nstep2cv_spin(ip, ihop, spin) = iter
2046 :
2047 2285 : end subroutine varpeq_record
2048 : !!***
2049 :
2050 : !----------------------------------------------------------------------
2051 :
2052 : !!****f* m_varpeq/varpeq_hop_setup
2053 : !! NAME
2054 : !! varpeq_hop_setup
2055 : !!
2056 : !! FUNCTION
2057 : !! Setup the hopping optimization.
2058 : !! This routine specfifies the set of B_q\nu images via linear
2059 : !! interpolation between an initial and final images.
2060 : !!
2061 : !! INPUTS
2062 : !! dtset<dataset_type>=All input variables for this dataset.
2063 : !!
2064 : !! OUTPUT
2065 : !!
2066 : !! SOURCE
2067 :
2068 1 : subroutine varpeq_hop_setup(self, dtset)
2069 :
2070 : !Arguments ------------------------------------
2071 : class(varpeq_t), target, intent(inout) :: self
2072 : type(dataset_type), intent(in) :: dtset
2073 :
2074 : !Local variables-------------------------------
2075 : !scalars
2076 103 : type(varpeq_t) :: vpq_from, vpq_to
2077 : class(polstate_t), pointer :: polstate
2078 : class(gqk_t), pointer :: gqk
2079 : integer, parameter :: master = 0
2080 : integer :: my_rank, comm, ierr
2081 : integer :: my_is, spin
2082 : real(dp) :: cpu, wall, gflops
2083 : !arrays
2084 : integer :: center(3)
2085 2 : complex(dp) :: b_spin_from(self%natom3, self%max_nq, self%nsppol)
2086 1 : complex(dp) :: b_spin_to(self%natom3, self%max_nq, self%nsppol)
2087 :
2088 : !----------------------------------------------------------------------
2089 :
2090 1 : call cwtime(cpu, wall, gflops, "start")
2091 :
2092 1 : comm = self%gstore%comm; my_rank = xmpi_comm_rank(comm)
2093 :
2094 1 : if (my_rank == master) then
2095 1 : call vpq_from%ncread(dtset%vpq_hop_from_filepath, xmpi_comm_self)
2096 :
2097 1 : if (dtset%vpq_hop_to_filepath /= ABI_NOFILE) then
2098 1 : call vpq_to%ncread(dtset%vpq_hop_to_filepath, xmpi_comm_self)
2099 : else
2100 0 : call vpq_to%ncread(dtset%vpq_hop_from_filepath, xmpi_comm_self)
2101 : endif
2102 :
2103 :
2104 : ! TODO: allow BZ mismatch by adding a BZ interpolation step?
2105 1 : call self%compare(vpq_from, bz_mismatch=.false.)
2106 1 : call self%compare(vpq_to, bz_mismatch=.false.)
2107 :
2108 877 : b_spin_from(:,:,:) = vpq_from%b_spin(:,:, dtset%vpq_hop_from_ip, :)
2109 877 : b_spin_to(:,:,:) = vpq_to%b_spin(:,:, dtset%vpq_hop_to_ip, :)
2110 :
2111 1 : call vpq_from%free()
2112 1 : call vpq_to%free()
2113 : endif
2114 :
2115 1 : call xmpi_bcast(b_spin_from, master, comm, ierr)
2116 1 : call xmpi_bcast(b_spin_to, master, comm, ierr)
2117 :
2118 4 : center(:) = self%ngkpt(:) / 2
2119 :
2120 2 : do my_is=1,self%gstore%my_nspins
2121 1 : spin = self%gstore%my_spins(my_is)
2122 1 : gqk => self%gstore%gqk(my_is)
2123 1 : polstate => self%polstate(spin)
2124 :
2125 : call polstate%load_b(b_spin_from(:,:,spin), ip=1, &
2126 4 : trvec=(center(:) - dtset%vpq_hop_from_site(:)))
2127 :
2128 1 : if (dtset%vpq_hop_to_filepath /= ABI_NOFILE) then
2129 : call polstate%load_b(b_spin_to(:,:,spin), ip=self%nstates, &
2130 4 : trvec=(center(:) - dtset%vpq_hop_to_site(:) + dtset%vpq_hop_vec(:)))
2131 : else
2132 : call polstate%load_b(b_spin_to(:,:,spin), ip=self%nstates, &
2133 0 : trvec=(center(:) - dtset%vpq_hop_from_site(:) + dtset%vpq_hop_vec(:)))
2134 : endif
2135 :
2136 2 : call polstate%linterp_b()
2137 : enddo
2138 :
2139 1 : self%fix_displ = .true.
2140 :
2141 1 : call cwtime_report(" varpeq: hop_setup", cpu, wall, gflops)
2142 :
2143 1 : end subroutine varpeq_hop_setup
2144 : !!***
2145 :
2146 : !----------------------------------------------------------------------
2147 :
2148 : !!****f* m_varpeq/varpeq_init
2149 : !! NAME
2150 : !! varpeq_init
2151 : !!
2152 : !! FUNCTION
2153 : !! Initialize the oboject by setting basic variables, allocate dynamic arrays.
2154 : !!
2155 : !! INPUTS
2156 : !! gstore<gstore_t>=Electron-phonon matrix elements and related quantities.
2157 : !! dtset<dataset_type>=All input variables for this dataset.
2158 : !!
2159 : !! OUTPUT
2160 : !!
2161 : !! SOURCE
2162 :
2163 9 : subroutine varpeq_init(self, gstore, dtset)
2164 :
2165 : !Arguments ------------------------------------
2166 : class(varpeq_t), target, intent(inout) :: self
2167 : type(gstore_t), target, intent(in) :: gstore
2168 : type(dataset_type), intent(in) :: dtset
2169 :
2170 : !Local variables-------------------------------
2171 : !scalars
2172 : character(len=5000) :: msg
2173 : class(gqk_t), pointer :: gqk
2174 : class(crystal_t), pointer :: cryst
2175 : class(polstate_t), pointer :: polstate
2176 : integer :: ierr, my_is, spin, bstart, bend, my_iq
2177 : real(dp) :: wtq, cpu, wall, gflops
2178 : !----------------------------------------------------------------------
2179 :
2180 9 : call cwtime(cpu, wall, gflops, "start")
2181 :
2182 : ! Consistency check
2183 9 : ierr = 0
2184 9 : if (gstore%check_cplex_qkzone_gmode(2, "bz", "bz", "phonon") /= 0) then
2185 0 : ABI_ERROR_NOSTOP("The gstore object is inconsistent with varpeq. See messages above.", ierr)
2186 : end if
2187 9 : if (gstore%has_used_lgk /= 0) then
2188 0 : ABI_ERROR_NOSTOP("The varpeq formalism does not support use_lgk /=0 .", ierr)
2189 : end if
2190 9 : if (gstore%has_used_lgq /= 0) then
2191 0 : ABI_ERROR_NOSTOP("The varpeq formalism does not support use_lgq /=0 .", ierr)
2192 : end if
2193 9 : if (ierr > 1) then
2194 : write(msg,'(a,i0,5a)')&
2195 0 : 'Checking consistency of input data against itself gave ',ierr,' inconsistencies.',ch10,&
2196 0 : 'The details of the problems can be FOUND ABOVE (or in output or log file), in an earlier WARNING.',ch10,&
2197 0 : 'In parallel, the details might not even be printed there. Then, try running in sequential to see the details.'
2198 0 : ABI_ERROR(msg)
2199 : end if
2200 :
2201 9 : cryst => gstore%cryst
2202 :
2203 : ! Scalars
2204 : ! character
2205 9 : self%mode = dtset%vpq_mode
2206 9 : self%pkind = dtset%vpq_pkind
2207 9 : self%aseed = dtset%vpq_aseed
2208 : ! logical
2209 9 : self%restart = (dtset%eph_restart /= 0)
2210 9 : self%interp = (dtset%vpq_interp /= 0)
2211 9 : self%g0_flag = (dtset%vpq_avg_g /= 0)
2212 9 : self%translate = (dtset%vpq_translate /= 0)
2213 : ! integer
2214 9 : self%atloc = dtset%vpq_atloc
2215 9 : self%hop_nstep = dtset%vpq_hop_nstep
2216 9 : self%nstep = dtset%vpq_nstep
2217 9 : self%nstep_ort = dtset%vpq_nstep_ort
2218 9 : self%nsppol = gstore%nsppol
2219 9 : self%nstates = dtset%vpq_nstates
2220 9 : self%natom3 = cryst%natom*3
2221 18 : self%max_nk = maxval(gstore%glob_nk_spin)
2222 18 : self%max_nq = maxval(gstore%glob_nq_spin)
2223 18 : self%max_nb = maxval(gstore%brange_k_spin(2,:) - gstore%brange_k_spin(1,:)) + 1
2224 9 : self%frohl_ntheta = dtset%eph_frohl_ntheta
2225 : ! real
2226 9 : self%hop_tolgrs = dtset%vpq_hop_tolgrs
2227 9 : self%hop_ts = dtset%vpq_hop_ts
2228 9 : self%tolgrs = dtset%vpq_tolgrs
2229 9 : self%mixing_factor = dtset%vpq_mix_fact
2230 :
2231 : ! Static arrays
2232 : ! integer
2233 36 : self%ngkpt(:) = dtset%ngkpt(:)
2234 :
2235 : ! Dynamic arrays
2236 : ! integer
2237 27 : ABI_MALLOC(self%nk_spin, (gstore%nsppol))
2238 27 : ABI_MALLOC(self%nq_spin, (gstore%nsppol))
2239 27 : ABI_MALLOC(self%nb_spin, (gstore%nsppol))
2240 27 : ABI_MALLOC(self%brange_spin, (2, gstore%nsppol))
2241 18 : self%nk_spin(:) = gstore%glob_nk_spin(:)
2242 18 : self%nq_spin(:) = gstore%glob_nq_spin(:)
2243 18 : self%nb_spin(:) = gstore%brange_k_spin(2,:) - gstore%brange_k_spin(1,:) + 1
2244 36 : self%brange_spin(:,:) = gstore%brange_k_spin(:,:)
2245 :
2246 45 : ABI_MALLOC(self%cvflag_spin, (self%nstates, self%hop_nstep, gstore%nsppol))
2247 45 : ABI_MALLOC(self%nstep2cv_spin, (self%nstates, self%hop_nstep, gstore%nsppol))
2248 27 : ABI_MALLOC(self%hop_nstep2cv_spin, (gstore%nsppol))
2249 186 : self%cvflag_spin(:,:,:) = zero
2250 186 : self%nstep2cv_spin(:,:,:) = zero
2251 18 : self%hop_nstep2cv_spin(:) = zero
2252 :
2253 36 : ABI_MALLOC(self%k2ibz_spin, (self%max_nk, gstore%nsppol))
2254 36 : ABI_MALLOC(self%q2ibz_spin, (self%max_nq, gstore%nsppol))
2255 :
2256 : ! real
2257 27 : ABI_MALLOC(self%chrgat, (cryst%natom))
2258 45 : ABI_MALLOC(self%hop_hist_spin, (5, self%nstates, self%hop_nstep, gstore%nsppol))
2259 54 : ABI_MALLOC(self%scf_hist_spin, (6, self%nstep, self%nstates, self%hop_nstep, gstore%nsppol))
2260 27 : self%chrgat(:) = dtset%chrgat(:)
2261 911 : self%hop_hist_spin(:,:,:,:) = zero
2262 98641 : self%scf_hist_spin(:,:,:,:,:) = zero
2263 :
2264 36 : ABI_MALLOC(self%kpts_spin, (3, self%max_nk, gstore%nsppol))
2265 36 : ABI_MALLOC(self%qpts_spin, (3, self%max_nq, gstore%nsppol))
2266 27 : ABI_MALLOC(self%erange_spin, (gstore%nsppol))
2267 18 : self%erange_spin(:) = zero
2268 9 : if (gstore%kfilter == "erange") then
2269 1 : self%use_filter = .true.
2270 1 : if (dtset%vpq_pkind == "hole") then
2271 0 : self%erange_spin(:) = gstore%erange_spin(1,:)
2272 : else
2273 2 : self%erange_spin(:) = gstore%erange_spin(2,:)
2274 : endif
2275 : endif
2276 :
2277 : ! complex
2278 54 : ABI_MALLOC(self%a_spin, (self%max_nb, self%max_nk, self%nstates, gstore%nsppol))
2279 54 : ABI_MALLOC(self%b_spin, (self%natom3, self%max_nq, self%nstates, gstore%nsppol))
2280 :
2281 : ! Datatypes and pointers
2282 9 : self%gstore => gstore
2283 9 : call cryst%copy(self%cryst)
2284 9 : self%gaps = gstore%ebands%get_gaps(ierr)
2285 :
2286 : call self%scell%init(cryst%natom, gstore%ebands%kptrlatt, cryst%rprimd, &
2287 9 : cryst%typat, cryst%xcart, cryst%znucl, xyz_order="xyz")
2288 9 : ABI_CHECK(gstore%same_nbands(msg), sjoin("VarPEq requires nb_k == nb_kq.", msg))
2289 :
2290 : ! Initialize polaronic states for each spin
2291 495 : ABI_MALLOC(self%polstate, (gstore%my_nspins))
2292 18 : do my_is=1,gstore%my_nspins
2293 9 : spin = gstore%my_spins(my_is)
2294 9 : gqk => gstore%gqk(my_is)
2295 9 : polstate => self%polstate(my_is)
2296 :
2297 : ! Scalars
2298 : ! character
2299 9 : polstate%aseed = dtset%vpq_aseed
2300 : ! logical
2301 9 : polstate%translate = (dtset%vpq_translate /= 0)
2302 : ! integers
2303 9 : polstate%np = dtset%vpq_nstates
2304 9 : polstate%nkbz = gstore%nkbz
2305 9 : polstate%nqbz = gstore%nqbz
2306 : ! real
2307 9 : polstate%e_frohl = zero
2308 9 : polstate%efilter = dtset%vpq_efilter
2309 :
2310 : ! Static arrays
2311 : ! integer
2312 36 : polstate%ngkpt(:) = dtset%ngkpt(:)
2313 : ! real
2314 27 : polstate%gpr_energy(:) = dtset%vpq_gpr_energy(:)
2315 36 : polstate%gpr_length(:) = dtset%vpq_gpr_length(:)
2316 :
2317 : ! Dynamic arrays
2318 : ! logical
2319 27 : ABI_MALLOC(polstate%has_prev_grad, (dtset%vpq_nstates))
2320 28 : polstate%has_prev_grad(:) = .false.
2321 : ! real
2322 27 : ABI_MALLOC(polstate%gradres, (dtset%vpq_nstates))
2323 27 : ABI_MALLOC(polstate%enterms, (4, dtset%vpq_nstates))
2324 :
2325 36 : ABI_MALLOC(polstate%eig, (gqk%nb_k, gstore%ebands%nkpt))
2326 : msg = sjoin(self%gaps%errmsg_spin(spin), &
2327 9 : "VarPEq is incompatible with metals and requires band gap.")
2328 9 : ABI_CHECK(self%gaps%ierr(spin) == 0, msg)
2329 :
2330 9 : bstart = gstore%brange_k_spin(1, spin)
2331 9 : bend = bstart + gqk%nb_k - 1
2332 1 : select case(dtset%vpq_pkind)
2333 : case ("electron")
2334 : polstate%eig = &
2335 170 : gstore%ebands%eig(bstart:bend, :, spin) - self%gaps%cb_min(spin)
2336 1 : polstate%psign = 1
2337 : case ("hole")
2338 : polstate%eig = &
2339 248 : -(gstore%ebands%eig(bstart:bend, :, spin) - self%gaps%vb_max(spin))
2340 17 : polstate%psign = -1
2341 : end select
2342 :
2343 27 : ABI_MALLOC(polstate%my_g0, (gqk%my_npert))
2344 63 : polstate%my_g0(:) = zero
2345 :
2346 27 : ABI_MALLOC(polstate%my_qpts, (3, gqk%my_nq))
2347 1190 : do my_iq=1,gqk%my_nq
2348 1190 : call gqk%myqpt(my_iq, gstore, wtq, polstate%my_qpts(:, my_iq))
2349 : enddo
2350 :
2351 36 : ABI_MALLOC(polstate%displ, (3, self%scell%natom, dtset%vpq_nstates))
2352 :
2353 : ! complex
2354 45 : ABI_MALLOC(polstate%my_a, (gqk%nb_k, gqk%my_nk, dtset%vpq_nstates))
2355 45 : ABI_MALLOC(polstate%a_glob, (gqk%nb_k, gqk%glob_nk, dtset%vpq_nstates))
2356 45 : ABI_MALLOC(polstate%my_b, (gqk%my_npert, gqk%my_nq, dtset%vpq_nstates))
2357 36 : ABI_MALLOC(polstate%my_prev_b, (gqk%my_npert, gqk%my_nq))
2358 45 : ABI_MALLOC(polstate%my_prev_b_hop, (gqk%my_npert, gqk%my_nq, dtset%vpq_nstates))
2359 36 : ABI_MALLOC(polstate%my_pc, (gqk%nb_k, gqk%my_nk))
2360 36 : ABI_MALLOC(polstate%my_grad, (gqk%nb_k, gqk%my_nk))
2361 36 : ABI_MALLOC(polstate%my_prev_grad, (gqk%nb_k, gqk%my_nk))
2362 36 : ABI_MALLOC(polstate%my_pcgrad, (gqk%nb_k, gqk%my_nk))
2363 36 : ABI_MALLOC(polstate%my_prev_pcgrad, (gqk%nb_k, gqk%my_nk))
2364 36 : ABI_MALLOC(polstate%my_pcjgrad, (gqk%nb_k, gqk%my_nk))
2365 36 : ABI_MALLOC(polstate%my_prev_pcjgrad, (gqk%nb_k, gqk%my_nk))
2366 36 : ABI_MALLOC(polstate%pcjgrad_glob, (gqk%nb_k, gqk%glob_nk))
2367 45 : ABI_MALLOC(polstate%my_phgrad, (gqk%my_npert, gqk%my_nq, dtset%vpq_nstates))
2368 45 : ABI_MALLOC(polstate%my_eff_phforce, (gqk%my_npert, gqk%my_nq, dtset%vpq_nstates))
2369 27 : ABI_MALLOC(polstate%hop_ts, (dtset%vpq_nstates))
2370 :
2371 : ! Datatypes ans pointers
2372 9 : polstate%gqk => gqk
2373 9 : polstate%my_kpts => gqk%my_kpts(:,:)
2374 9 : polstate%krank_kpts = polstate%get_krank_glob("k", gstore%ebands%kptrlatt)
2375 9 : polstate%krank_qpts = polstate%get_krank_glob("q", gstore%ebands%kptrlatt)
2376 :
2377 18 : call cryst%copy(polstate%cryst)
2378 :
2379 : enddo
2380 :
2381 : ! This is needed to fill kpts_spin/qpts_spin before the calculation, as they may
2382 : ! be required by the varpeq_load subroutine
2383 9 : call xmpi_barrier(gstore%comm)
2384 9 : call self%collect()
2385 :
2386 9 : self%is_complete = .true.
2387 :
2388 9 : call cwtime_report(" varpeq: init", cpu, wall, gflops)
2389 :
2390 18 : end subroutine varpeq_init
2391 : !!***
2392 :
2393 : !----------------------------------------------------------------------
2394 :
2395 : !!****f* m_varpeq/varpeq_calc_fravg
2396 : !! NAME
2397 : !! varpeq_calc_fravg
2398 : !!
2399 : !! FUNCTION
2400 : !! Calculate the average Fr\"ohlich long-range contribution to the polaron
2401 : !! binding energy at Gamma using spherical integration in the spherical region
2402 : !! arond Gamma-point.
2403 : !!
2404 : !! INPUTS
2405 : !! avg_g0 [optional]=If .true., average electron-phonon matrix elements at
2406 : !! Gamma-point. Default: .true.
2407 : !!
2408 : !! OUTPUT
2409 : !!
2410 : !! SOURCE
2411 :
2412 9 : subroutine varpeq_calc_fravg(self, avg_g0)
2413 :
2414 : !Arguments ------------------------------------
2415 : class(varpeq_t), target, intent(inout) :: self
2416 : logical, optional, intent(in) :: avg_g0
2417 :
2418 : !Local variables-------------------------------
2419 : !scalars
2420 : type(ifc_type), pointer :: ifc
2421 : class(polstate_t), pointer :: polstate
2422 : class(gqk_t), pointer :: gqk
2423 : integer :: comm, my_rank, nproc, ierr
2424 : integer :: ntheta, nphi, angl_size
2425 : integer :: iang, iatom
2426 : integer :: nu, my_is, my_pert, pert_glob
2427 : real(dp) :: inv_qepsq, wqnu, prefactor
2428 : real(dp) :: cpu, wall, gflops
2429 : complex(dp) :: cnum
2430 : !arrays
2431 : real(dp) :: qpt_cart(3)
2432 9 : real(dp), allocatable :: phfreq(:), displ_cart(:, :, :, :)
2433 9 : real(dp), allocatable :: qvers_cart(:, :)
2434 9 : real(dp), allocatable :: angweight(:)
2435 : real(dp), allocatable :: e_frohl_mode(:)
2436 : complex(dp) :: cp3(3)
2437 : !----------------------------------------------------------------------
2438 :
2439 9 : call cwtime(cpu, wall, gflops, "start")
2440 :
2441 9 : comm = self%gstore%comm
2442 9 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2443 :
2444 : ! Integrate electron-phonon energy in the spherical region around Gamma
2445 27 : ABI_MALLOC(phfreq, (self%natom3))
2446 18 : ABI_MALLOC(e_frohl_mode, (self%natom3))
2447 36 : ABI_MALLOC(displ_cart, (2, 3, self%cryst%natom, self%natom3))
2448 :
2449 9 : ifc => self%gstore%ifc
2450 :
2451 : ! Create the mesh for spherical integrtion
2452 9 : ntheta = self%frohl_ntheta
2453 9 : nphi = 2*ntheta
2454 9 : call ylm_angular_mesh(ntheta, nphi, angl_size, qvers_cart, angweight)
2455 :
2456 : ! Integrate the contribution from each phonon mode
2457 : ! This is similar to what is done in the src/78_eph/m_sigmaph.f90 module
2458 63 : e_frohl_mode(:) = zero
2459 18441 : do iang=1,angl_size
2460 18432 : if (mod(iang, nproc) /= my_rank) cycle
2461 :
2462 73728 : qpt_cart = qvers_cart(:, iang)
2463 294912 : inv_qepsq = one / dot_product(qpt_cart, matmul(ifc%dielt, qpt_cart))
2464 18432 : call ifc%fourq(self%cryst, qpt_cart, phfreq, displ_cart, nanaqdir="cart")
2465 :
2466 73737 : do nu=4,self%natom3
2467 55296 : wqnu = phfreq(nu)
2468 :
2469 55296 : cp3(:) = zero
2470 165888 : do iatom=1,self%cryst%natom
2471 331776 : cp3(:) = cp3(:) + matmul(ifc%zeff(:, :, iatom), &
2472 4036608 : cmplx(displ_cart(1,:,iatom,nu), displ_cart(2,:,iatom,nu), kind=dp))
2473 : enddo
2474 221184 : cnum = dot_product(qpt_cart, cp3)
2475 55296 : if (abs(cnum) < tol12) cycle
2476 :
2477 : e_frohl_mode(nu) = e_frohl_mode(nu) + &
2478 73728 : angweight(iang)*(abs(cnum)*inv_qepsq/wqnu)**2
2479 : enddo
2480 : enddo
2481 9 : call xmpi_sum(e_frohl_mode, comm, ierr)
2482 :
2483 : prefactor = (eight * pi / self%cryst%ucvol) * &
2484 9 : (three / (four_pi * self%cryst%ucvol * self%gstore%nqbz))**third
2485 :
2486 63 : e_frohl_mode(:) = prefactor*e_frohl_mode(:)
2487 63 : self%e_frohl = sum(e_frohl_mode)
2488 : ! For an electron polaron, the correction has to be negative
2489 9 : if (self%pkind == "electron") self%e_frohl = -self%e_frohl
2490 :
2491 : ! If corresponding flag is provided, average the matrix elements at Gamma
2492 9 : if (present(avg_g0) .and. avg_g0) then
2493 12 : do my_is=1,self%gstore%my_nspins
2494 6 : polstate => self%polstate(my_is)
2495 6 : gqk => polstate%gqk
2496 :
2497 6 : polstate%e_frohl = self%e_frohl
2498 48 : do my_pert=1,gqk%my_npert
2499 36 : pert_glob = gqk%my_pert_start + my_pert - 1
2500 36 : wqnu = phfreq(pert_glob)
2501 : polstate%my_g0(my_pert) = &
2502 42 : sqrt(polstate%nqbz * wqnu * e_frohl_mode(pert_glob) / two)
2503 :
2504 : enddo
2505 : enddo
2506 : endif
2507 :
2508 9 : ABI_FREE(phfreq)
2509 9 : ABI_FREE(e_frohl_mode)
2510 9 : ABI_FREE(displ_cart)
2511 9 : ABI_FREE(qvers_cart)
2512 9 : ABI_FREE(angweight)
2513 :
2514 9 : call cwtime_report(" varpeq: calc_fravg", cpu, wall, gflops)
2515 :
2516 27 : end subroutine varpeq_calc_fravg
2517 : !!***
2518 :
2519 : !----------------------------------------------------------------------
2520 :
2521 : !!****f* m_varpeq/polstate_free
2522 : !! NAME
2523 : !! polstate_free
2524 : !!
2525 : !! FUNCTION
2526 : !! Free dynamic memory
2527 : !!
2528 : !! OUTPUT
2529 : !!
2530 : !! SOURCE
2531 :
2532 9 : subroutine polstate_free(self)
2533 :
2534 : !Arguments ------------------------------------
2535 : class(polstate_t), intent(inout) :: self
2536 : !----------------------------------------------------------------------
2537 :
2538 : ! Free allocatable arrays
2539 :
2540 : ! logical
2541 9 : ABI_SFREE(self%has_prev_grad)
2542 : ! real
2543 9 : ABI_SFREE(self%gradres)
2544 9 : ABI_SFREE(self%enterms)
2545 9 : ABI_SFREE(self%eig)
2546 9 : ABI_SFREE(self%my_g0)
2547 9 : ABI_SFREE(self%my_qpts)
2548 9 : ABI_SFREE(self%displ)
2549 :
2550 : ! complex
2551 9 : ABI_SFREE(self%my_a)
2552 9 : ABI_SFREE(self%a_glob)
2553 9 : ABI_SFREE(self%my_b)
2554 9 : ABI_SFREE(self%my_prev_b)
2555 9 : ABI_SFREE(self%my_prev_b_hop)
2556 9 : ABI_SFREE(self%my_pc)
2557 9 : ABI_SFREE(self%my_grad)
2558 9 : ABI_SFREE(self%my_pcgrad)
2559 9 : ABI_SFREE(self%my_prev_pcgrad)
2560 9 : ABI_SFREE(self%my_prev_grad)
2561 9 : ABI_SFREE(self%my_pcjgrad)
2562 9 : ABI_SFREE(self%my_prev_pcjgrad)
2563 9 : ABI_SFREE(self%pcjgrad_glob)
2564 9 : ABI_SFREE(self%my_phgrad)
2565 9 : ABI_SFREE(self%my_eff_phforce)
2566 9 : ABI_SFREE(self%hop_ts)
2567 :
2568 : ! Free local datatypes & nullify pointers
2569 9 : self%my_kpts => null()
2570 9 : self%gqk => null()
2571 :
2572 9 : call self%krank_kpts%free()
2573 9 : call self%krank_qpts%free()
2574 :
2575 9 : call self%cryst%free()
2576 :
2577 9 : end subroutine polstate_free
2578 : !!***
2579 :
2580 : !----------------------------------------------------------------------
2581 :
2582 : !!****f* m_varpeq/polstate_setup
2583 : !! NAME
2584 : !! polstate_setup
2585 : !!
2586 : !! FUNCTION
2587 : !! Setup optimization process at a given polaronic state.
2588 : !! This routine specfifies an initial electronic vector A_nk, either by
2589 : !! initializaing it by a pre-determined algorithm or loading from a
2590 : !! *VPQ.nc netcdf file.
2591 : !!
2592 : !! INPUTS
2593 : !! ip=Index of the polaronic state.
2594 : !! a_src(self%gqk%nb_k, self%gqk%glob_nk) [optional]=Global A_nk coefficients at
2595 : !! this state, which have to be provided if load_src=.true.
2596 : !! load_src [optional]=.true. if A_nk is initialized from an external source,
2597 : !! e.g. loaded from file. Default: .false.
2598 : !! atloc [optional]=Index of an atom where charge is localized.
2599 : !! chrgat(natom) [optional]=Charge of the atoms.
2600 : !! scell<supercell_type> [optional]=Supercell data structure.
2601 : !! cryst<crystal_t> [optional]=Crystal data structure.
2602 : !! nstep_loc [optional]=Number of iterations for charge localization.
2603 : !!
2604 : !! OUTPUT
2605 : !!
2606 : !! SOURCE
2607 :
2608 19 : subroutine polstate_setup(self, ip, a_src, load, atloc, chrgat, scell, cryst, nstep_loc)
2609 :
2610 : !Arguments ------------------------------------
2611 : class(polstate_t), target, intent(inout) :: self
2612 : integer, intent(in) :: ip
2613 : logical, optional, intent(in) :: load
2614 : complex(dp), optional, intent(in) :: a_src(self%gqk%nb_k, self%gqk%glob_nk)
2615 : integer, optional, intent(in) :: atloc
2616 : class(supercell_type), optional, intent(in) :: scell
2617 : class(crystal_t), optional, intent(in) :: cryst
2618 : real(dp), optional, intent(in) :: chrgat(:)
2619 : integer, optional, intent(in) :: nstep_loc
2620 :
2621 : !Local variables-------------------------------
2622 : real(dp) :: a_sqnorm
2623 : class(gqk_t), pointer :: gqk
2624 : !----------------------------------------------------------------------
2625 :
2626 19 : gqk => self%gqk
2627 :
2628 19 : if (present(load) .and. load) then
2629 0 : ABI_CHECK(present(a_src), "polstate_setup: A_nk is expected but not provided")
2630 : !print *, "load"
2631 0 : call self%load_a(a_src, ip)
2632 : else
2633 19 : call self%seed_a(self%aseed, ip, atloc, chrgat, scell, cryst, nstep_loc)
2634 : endif
2635 :
2636 19 : if (self%efilter > zero) call self%filter("a", ip)
2637 :
2638 : ! Orthogonalize current states to the previous ones
2639 19 : call self%ort_to_states(self%my_a(:,:,ip), 1, ip-1, ip, tr_flag=self%translate)
2640 :
2641 : ! Normalize A_nk at current polaronic state
2642 19 : a_sqnorm = self%get_sqnorm("a", ip)
2643 19 : if (a_sqnorm > tol12) then
2644 6749 : self%my_a(:,:,ip) = sqrt(self%nkbz/a_sqnorm) * self%my_a(:,:,ip)
2645 : endif
2646 :
2647 19 : end subroutine polstate_setup
2648 : !!***
2649 :
2650 : !----------------------------------------------------------------------
2651 :
2652 : !!****f* m_varpeq/polstate_update_a
2653 : !! NAME
2654 : !! polstate_update_a
2655 : !!
2656 : !! FUNCTION
2657 : !! Update the vector of electronic coefficients A_nk by line minimization in
2658 : !! the pcj direction for current polaronic state.
2659 : !!
2660 : !! INPUTS
2661 : !! ip=Index of a polaronic state.
2662 : !!
2663 : !! OUTPUT
2664 : !!
2665 : !! SOURCE
2666 :
2667 2141 : subroutine polstate_update_a(self, ip)
2668 :
2669 : !Arguments ------------------------------------
2670 : class(polstate_t), intent(inout) :: self
2671 : integer, intent(in) :: ip
2672 :
2673 : !Local variables-------------------------------
2674 : real(dp) :: theta
2675 : !----------------------------------------------------------------------
2676 :
2677 2141 : theta = self%get_lm_theta(ip)
2678 : self%my_a(:,:,ip) = &
2679 979741 : cos(theta)*self%my_a(:,:,ip) + sin(theta)*self%my_pcjgrad(:,:)
2680 :
2681 2141 : end subroutine polstate_update_a
2682 : !!***
2683 :
2684 : !----------------------------------------------------------------------
2685 :
2686 : !!****f* m_varpeq/polstate_update_pc
2687 : !! NAME
2688 : !! polstate_update_pc
2689 : !!
2690 : !! FUNCTION
2691 : !! Update the preconditioner for the present configuration. Changing this
2692 : !! procedure may significantly improve (or worsen) the optimization process.
2693 : !!
2694 : !! INPUTS
2695 : !! ip=Index of a polaronic state.
2696 : !!
2697 : !! OUTPUT
2698 : !!
2699 : !! SOURCE
2700 :
2701 2141 : subroutine polstate_update_pc(self, ip)
2702 :
2703 : !Arguments ------------------------------------
2704 : class(polstate_t), intent(inout) :: self
2705 : integer, intent(in) :: ip
2706 :
2707 : !Local variables-------------------------------
2708 : class(gqk_t), pointer :: gqk
2709 : integer :: my_ik, ik_ibz, ib
2710 : real(dp) :: eps
2711 : !----------------------------------------------------------------------
2712 :
2713 2141 : gqk => self%gqk
2714 :
2715 2141 : eps = self%enterms(4, ip)
2716 247787 : do my_ik=1,gqk%my_nk
2717 245646 : ik_ibz = gqk%my_k2ibz(1, my_ik)
2718 979741 : do ib=1,gqk%nb_k
2719 : self%my_pc(ib, my_ik) = &
2720 977600 : one/abs(self%eig(ib, ik_ibz) - two*abs(self%e_frohl) + abs(eps))
2721 : !self%my_pc(ib, my_ik) = one
2722 : enddo
2723 : enddo
2724 :
2725 2141 : end subroutine polstate_update_pc
2726 : !!***
2727 :
2728 : !----------------------------------------------------------------------
2729 :
2730 : !!****f* m_varpeq/polstate_ort_to_states
2731 : !! NAME
2732 : !! polstate_ort_to_states
2733 : !!
2734 : !! FUNCTION
2735 : !! Orthognongalize a vector wrt polaronic states using the Gram-Schmidt process.
2736 : !! The orthogonalization is performed for a range of states, specified by arguments.
2737 : !!
2738 : !! INPUTS
2739 : !! my_v(:,:)=Vetor to be orthogonalized
2740 : !! istart=Index of starting polaronic state
2741 : !! iend=Index of final polaronic state
2742 : !! this_ip=Index of current polaronic state
2743 : !! tr_flag=.true. if orthogonalization must include all states invariant by
2744 : !! translations inside a supercell
2745 : !!
2746 : !! OUTPUT
2747 : !!
2748 : !! SOURCE
2749 :
2750 6257 : subroutine polstate_ort_to_states(self, my_v, istart, iend, this_ip, tr_flag)
2751 :
2752 : !Arguments ------------------------------------
2753 : class(polstate_t), intent(inout) :: self
2754 : logical, intent(in) :: tr_flag
2755 : integer, intent(in) :: istart, iend, this_ip
2756 : complex(dp), intent(inout) :: my_v(self%gqk%nb_k, self%gqk%my_nk)
2757 :
2758 : !Local variables-------------------------------
2759 : class(gqk_t), pointer :: gqk
2760 : integer :: ip, my_ik, vx, vy, vz
2761 : complex(dp) :: phase, proj
2762 : integer :: tr_vec(3), ngkpt_tr(3)
2763 : real(dp) :: kpt(3)
2764 12514 : complex(dp) :: a_tr(self%gqk%nb_k, self%gqk%my_nk)
2765 : !----------------------------------------------------------------------
2766 :
2767 6257 : gqk => self%gqk
2768 :
2769 : !! TODO: optimize
2770 : !ngkpt_tr(:) = 1
2771 : !if (tr_flag) ngkpt_tr(:) = self%ngkpt(:)
2772 :
2773 28472 : do ip=istart,iend
2774 :
2775 22215 : if ((tr_flag) .and. ip /= this_ip) then
2776 732 : ngkpt_tr(:) = self%ngkpt(:)
2777 : else
2778 88128 : ngkpt_tr(:) = 1
2779 : endif
2780 :
2781 51053 : do vx=1,ngkpt_tr(1)
2782 22581 : tr_vec(1) = vx - 1
2783 68475 : do vy=1,ngkpt_tr(2)
2784 23679 : tr_vec(2) = vy - 1
2785 73233 : do vz=1,ngkpt_tr(3)
2786 26973 : tr_vec(3) = vz - 1
2787 :
2788 2870422 : do my_ik=1,gqk%my_nk
2789 11373796 : kpt(:) = self%my_kpts(:, my_ik)
2790 11373796 : phase = exp(j_dpc*sum(kpt(:)*tr_vec(:))*two_pi)
2791 11390801 : a_tr(:, my_ik) = phase*self%my_a(:, my_ik, ip)
2792 : enddo
2793 :
2794 26973 : proj = get_proj_(a_tr)
2795 11414480 : my_v(:,:) = my_v(:,:) - proj * a_tr(:,:)
2796 :
2797 : enddo
2798 : enddo
2799 : enddo
2800 :
2801 : enddo
2802 :
2803 : !----------------------------------------------------------------------
2804 : contains
2805 :
2806 26973 : complex(dp) function get_proj_(my_u) result(proj)
2807 :
2808 : complex(dp), intent(in) :: my_u(gqk%nb_k, gqk%my_nk)
2809 : integer :: ierr
2810 : real(dp) :: u_sqnorm
2811 : !----------------------------------------------------------------------
2812 11390801 : u_sqnorm = sum(abs(my_u(:,:))**2)
2813 26973 : call xmpi_sum(u_sqnorm, gqk%kpt_comm%value, ierr)
2814 :
2815 26973 : proj = zero
2816 26973 : if (u_sqnorm > tol12) then
2817 11390801 : proj = sum(conjg(my_u(:,:))*my_v(:,:))
2818 26973 : call xmpi_sum(proj, gqk%kpt_comm%value, ierr)
2819 26973 : proj = proj/u_sqnorm
2820 : endif
2821 26973 : end function get_proj_
2822 :
2823 : end subroutine polstate_ort_to_states
2824 : !!***
2825 :
2826 : !----------------------------------------------------------------------
2827 :
2828 : !!****f* m_varpeq/polstate_get_lm_theta
2829 : !! NAME
2830 : !! polstate_get_lm_theta
2831 : !!
2832 : !! FUNCTION
2833 : !! Calculate line minimization parameter theta that minimizes polaron binding
2834 : !! energy at the next electronic configuration:
2835 : !! A^n = A^(n-1)*cos(theta) + D^(n-1)*sin(theta), where
2836 : !! A^(n-1) and D^(n-1) are current electronic configuration and gradient.
2837 : !!
2838 : !! INPUTS
2839 : !! ip=Index of a polaronic state.
2840 : !!
2841 : !! OUTPUT
2842 : !! theta=Line minimization parameter.
2843 : !!
2844 : !! SOURCE
2845 :
2846 2141 : real(dp) function polstate_get_lm_theta(self, ip) result(theta)
2847 :
2848 : !Arguments ------------------------------------
2849 : class(polstate_t), intent(inout) :: self
2850 : integer, intent(in) :: ip
2851 :
2852 : !Local variables-------------------------------
2853 : !scalars
2854 : class(gqk_t), pointer :: gqk
2855 : real(dp) :: sqnorm, e1, e2, eps
2856 :
2857 : !----------------------------------------------------------------------
2858 :
2859 2141 : gqk => self%gqk
2860 :
2861 : ! Orthogonalize pcj direction to the current state and normalize
2862 2141 : call self%ort_to_states(self%my_pcjgrad, ip, ip, ip, tr_flag=.false.)
2863 :
2864 2141 : sqnorm = self%get_sqnorm('pcjgrad', ip)
2865 979741 : self%my_pcjgrad(:,:) = sqrt(self%nkbz/sqnorm)*self%my_pcjgrad(:,:)
2866 :
2867 : ! Calculation of theta requires globally available pcj direction
2868 2141 : call self%gather("pcjgrad", ip)
2869 :
2870 :
2871 : ! Line-minimization theta
2872 : ! E_pol(theta) = E_pol[A,B] - e1/2 + e1/2*cos(2*theta) + e2/2*sin(theta)
2873 2141 : eps = self%enterms(4, ip)
2874 : e1 = eps - self%calc_hpol(self%my_pcjgrad, self%pcjgrad_glob, &
2875 2141 : self%my_pcjgrad, self%my_b(:,:,ip))
2876 : e2 = &
2877 : self%calc_hpol(self%my_a(:,:,ip), self%a_glob(:,:,ip), self%my_pcjgrad, self%my_b(:,:,ip)) + &
2878 2141 : self%calc_hpol(self%my_pcjgrad, self%pcjgrad_glob, self%my_a(:,:,ip), self%my_b(:,:,ip))
2879 :
2880 2141 : theta = half*atan2(-e2, -e1)
2881 : !if (theta < zero) then
2882 : ! theta = theta + pi
2883 : !end if
2884 : !theta = half*theta
2885 :
2886 2141 : end function polstate_get_lm_theta
2887 : !!***
2888 :
2889 : !----------------------------------------------------------------------
2890 :
2891 : !!****f* m_varpeq/polstate_calc_pcjgrad
2892 : !! NAME
2893 : !! polstate_calc_pcjgrad
2894 : !!
2895 : !! FUNCTION
2896 : !! Calculate preconditioned conjugate gradient direction, orthogonal to all
2897 : !! already optimized polaronic states. The procedure is similar to the one
2898 : !! described in [Payne et al, Rev. Mod. Phys, 64, 4, 1045-1097 (1992)].
2899 : !!
2900 : !! INPUTS
2901 : !! ip=Index of a polaronic state.
2902 : !! ii=Iteration number.
2903 : !! nstep_ort=Iteration number, after which the orthogonality constraint
2904 : !! on all PREVIOUS states is lifted.
2905 : !!
2906 : !! OUTPUT
2907 : !!
2908 : !! SOURCE
2909 :
2910 2141 : subroutine polstate_calc_pcjgrad(self, ip, ii, nstep_ort)
2911 :
2912 : !Arguments ------------------------------------
2913 : class(polstate_t), intent(inout) :: self
2914 : integer, intent(in) :: ip, ii, nstep_ort
2915 :
2916 : !Local variables-------------------------------
2917 : class(gqk_t), pointer :: gqk
2918 : integer :: ierr
2919 : !real(dp) :: eps
2920 : complex(dp) :: beta, beta_num, beta_den
2921 : !----------------------------------------------------------------------
2922 :
2923 2141 : gqk => self%gqk
2924 :
2925 2141 : if (self%efilter > zero) call self%filter("grad", ip)
2926 :
2927 : ! Orthogonalize current gradient to all previous bands
2928 2141 : if (ii <= nstep_ort) then
2929 1956 : call self%ort_to_states(self%my_grad, 1, ip-1, ip, tr_flag=self%translate)
2930 : endif
2931 :
2932 : ! Precondtion vector
2933 979741 : self%my_pcgrad(:,:) = self%my_pc(:,:)*self%my_grad(:,:)
2934 : ! Orthogonalize to all bands
2935 2141 : if (ii <= nstep_ort) then
2936 1956 : call self%ort_to_states(self%my_pcgrad, 1, ip, ip, tr_flag=self%translate)
2937 : else
2938 185 : call self%ort_to_states(self%my_pcgrad, ip, ip, ip, tr_flag=.false.)
2939 : endif
2940 :
2941 : ! Conjugate gradient direction
2942 2141 : if (self%has_prev_grad(ip)) then
2943 : ! Polak-Ribiere coefficient
2944 : beta_num = &
2945 972992 : sum(conjg(self%my_pcgrad(:,:))*(self%my_grad(:,:) - self%my_prev_grad(:,:)))
2946 972992 : beta_den = sum(conjg(self%my_prev_pcgrad(:,:))*self%my_prev_grad(:,:))
2947 2122 : call xmpi_sum(beta_num, gqk%kpt_comm%value, ierr)
2948 2122 : call xmpi_sum(beta_den, gqk%kpt_comm%value, ierr)
2949 2122 : beta = beta_num / beta_den
2950 : !if (abs(aimag(beta)) < tol12) beta = real(beta, dp)
2951 :
2952 :
2953 972992 : self%my_pcjgrad(:,:) = self%my_pcgrad(:,:) + beta*self%my_prev_pcjgrad(:,:)
2954 : else
2955 6749 : self%my_pcjgrad(:,:) = self%my_pcgrad(:,:)
2956 : endif
2957 :
2958 : ! Save previous gradients
2959 979741 : self%my_prev_grad(:,:) = self%my_grad(:,:)
2960 979741 : self%my_prev_pcgrad(:,:) = self%my_pcgrad(:,:)
2961 979741 : self%my_prev_pcjgrad(:,:) = self%my_pcjgrad(:,:)
2962 2141 : self%has_prev_grad(ip) = .true.
2963 :
2964 2141 : end subroutine polstate_calc_pcjgrad
2965 : !!***
2966 :
2967 : !----------------------------------------------------------------------
2968 :
2969 : !!****f* m_varpeq/polstate_calc_grad
2970 : !! NAME
2971 : !! polstate_calc_grad
2972 : !!
2973 : !! FUNCTION
2974 : !! Calculate the electronic gradient D_nk at the current configuration.
2975 : !!
2976 : !! INPUTS
2977 : !! ip=Index of a polaronic state.
2978 : !!
2979 : !! OUTPUT
2980 : !!
2981 : !! SOURCE
2982 :
2983 2285 : subroutine polstate_calc_grad(self, ip)
2984 :
2985 : !Arguments ------------------------------------
2986 : class(polstate_t), intent(inout) :: self
2987 : integer, intent(in) :: ip
2988 :
2989 : !Local variables-------------------------------
2990 : !scalars
2991 : class(gqk_t), pointer :: gqk
2992 : logical :: q_gamma
2993 : integer :: ierr
2994 : integer :: my_iq, my_pert
2995 : integer :: my_ik, ik_ibz, ik_forw, ik_back, ib, jb
2996 : real(dp) :: eps, fact
2997 : complex(dp) :: a_forw, a_back
2998 : complex(dp) :: g_forw, g_back
2999 : complex(dp) :: b, g0
3000 : !arrays
3001 : real(dp) :: kpt(3), qpt(3), kpq(3), kmq(3)
3002 4570 : complex(dp) :: ak(self%gqk%nb_k), akq(self%gqk%nb_kq), akmq(self%gqk%nb_kq), bq(self%gqk%my_npert)
3003 2285 : complex(dp), allocatable :: gq_gathered(:,:,:,:)
3004 : !----------------------------------------------------------------------
3005 :
3006 2285 : gqk => self%gqk
3007 :
3008 : !ABI_MALLOC(gq_gathered, (gqk%my_npert, gqk%nb_k, gqk%nb_k, gqk%glob_nk))
3009 :
3010 : ! Scattering-dependent part
3011 1049507 : self%my_grad(:, :) = zero
3012 :
3013 2285 : fact = two * self%psign / real(self%nqbz, dp)
3014 :
3015 280113 : do my_iq=1,gqk%my_nq
3016 1111312 : qpt(:) = self%my_qpts(:, my_iq)
3017 1944796 : bq(:) = self%my_b(:, my_iq, ip)
3018 :
3019 : ! Check if q=\Gamma
3020 346465 : q_gamma = .false.
3021 346465 : if (all(abs(qpt) < tol6)) q_gamma = .true.
3022 :
3023 : ! For this q, gather all the k-distributed matrix elements
3024 277828 : call gqk%gather("q", my_iq, gq_gathered)
3025 :
3026 33678830 : do my_ik=1,gqk%my_nk
3027 130522294 : ak(:) = self%my_a(:,my_ik,ip)
3028 133604008 : kpt(:) = self%my_kpts(:, my_ik)
3029 :
3030 : ! Forward scattering
3031 : ! Find k+q-->k' index in krank_kpts
3032 133604008 : kpq(:) = kpt(:) + qpt(:)
3033 33401002 : ik_forw = self%krank_kpts%get_index(kpq)
3034 :
3035 : ! If erange filter was used in gstore, some transitions are not valid
3036 33401002 : if (ik_forw /= -1) then
3037 127899998 : akq(:) = self%a_glob(:, ik_forw, ip)
3038 :
3039 127899998 : do ib=1,gqk%nb_k
3040 :
3041 414871012 : do jb=1,gqk%nb_kq
3042 286971014 : a_forw = akq(jb)
3043 :
3044 2104607242 : do my_pert=1,gqk%my_npert
3045 1721826084 : b = bq(my_pert)
3046 :
3047 1721826084 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
3048 : ! Add long-range correction to matrix elements at Gamma
3049 1721826084 : g0 = self%my_g0(my_pert)
3050 1721826084 : if (q_gamma .and. (ib == jb)) then
3051 4704756 : g_forw = g_forw + g0
3052 : endif
3053 :
3054 : self%my_grad(ib, my_ik) = &
3055 2008797098 : self%my_grad(ib, my_ik) + a_forw*b*conjg(g_forw)
3056 :
3057 : enddo
3058 : enddo
3059 : enddo
3060 : endif ! Forward scattering
3061 :
3062 : ! Backward scattering
3063 : ! Find k-q-->k' index in krank_kpts
3064 133604008 : kmq(:) = kpt(:) - qpt(:)
3065 33401002 : ik_back = self%krank_kpts%get_index(kmq)
3066 :
3067 : ! If erange filter was used in gstore, some transitions are not valid
3068 33678830 : if (ik_back /= -1) then
3069 127899998 : akmq(:) = self%a_glob(:, ik_back, ip)
3070 :
3071 127899998 : do ib=1,gqk%nb_k
3072 :
3073 414871012 : do jb=1,gqk%nb_kq
3074 286971014 : a_back = akmq(jb)
3075 :
3076 2104607242 : do my_pert=1,gqk%my_npert
3077 1721826084 : b = bq(my_pert)
3078 :
3079 1721826084 : g_back = gq_gathered(my_pert, ib, jb, ik_back)
3080 : ! Add long-range correction to matrix elements at Gamma
3081 1721826084 : g0 = self%my_g0(my_pert)
3082 1721826084 : if (q_gamma .and. (ib == jb)) then
3083 4704756 : g_back = g_back + g0
3084 : endif
3085 :
3086 : self%my_grad(ib, my_ik) = &
3087 2008797098 : self%my_grad(ib, my_ik) + a_back*conjg(b)*g_back
3088 :
3089 : enddo
3090 : enddo
3091 : enddo
3092 : endif ! Backward scattering
3093 :
3094 : enddo
3095 :
3096 280113 : ABI_FREE(gq_gathered)
3097 : enddo
3098 2285 : call xmpi_sum(self%my_grad, gqk%qpt_pert_comm%value, ierr)
3099 : self%my_grad(:, :) = &
3100 1049507 : -two/(real(self%nkbz, dp)*real(self%nqbz, dp)) * self%my_grad(:, :)
3101 :
3102 : ! Scattering-independent part
3103 2285 : eps = self%enterms(4, ip)
3104 265381 : do my_ik=1,gqk%my_nk
3105 263096 : ik_ibz = gqk%my_k2ibz(1, my_ik)
3106 1047222 : ak(:) = self%my_a(:, my_ik, ip)
3107 :
3108 1049507 : do ib=1,gqk%nb_k
3109 : self%my_grad(ib, my_ik) = self%my_grad(ib, my_ik) + &
3110 1047222 : two/self%nkbz * (self%eig(ib, ik_ibz) - eps) * self%my_a(ib, my_ik, ip)
3111 :
3112 : enddo
3113 : enddo
3114 :
3115 : ! Here we're actually changing the gradient to the steepest descent direction
3116 : ! A bit messy but i dont want to make another array called "forces" or something
3117 : !self%my_grad(:,:) = -self%my_grad(:,:)
3118 :
3119 : !ABI_FREE(gq_gathered)
3120 :
3121 2285 : if (self%efilter > zero) call self%filter("grad", ip)
3122 :
3123 2285 : end subroutine polstate_calc_grad
3124 : !!***
3125 :
3126 : !----------------------------------------------------------------------
3127 :
3128 : !!****f* m_varpeq/polstate_localize
3129 : !! NAME
3130 : !! polstate_localize
3131 : !!
3132 : !! FUNCTION
3133 : !! Calculate vibrational coefficients, polaron energy terms and polaron energy
3134 : !! level at a current state.
3135 : !!
3136 : !! INPUTS
3137 : !! ip=Index of a polaronic state.
3138 : !! alpha=Mixing factor.
3139 : !! fix_displ [optional]=if .true., DO NOT update the phonon vector B_\qnu.
3140 : !! Relevant for hopping calculations. Defaults to. False.
3141 : !!
3142 : !! OUTPUT
3143 : !!
3144 : !! SOURCE
3145 :
3146 2285 : subroutine polstate_localize(self, ip, alpha, fix_displ)
3147 :
3148 : !Arguments ------------------------------------
3149 : class(polstate_t), intent(inout) :: self
3150 : integer, intent(in) :: ip
3151 : real(dp), intent(in) :: alpha
3152 : logical, optional, intent(in) :: fix_displ
3153 : !----------------------------------------------------------------------
3154 :
3155 : ! Calculation of B_qnu requires globally available A_nk
3156 2285 : call self%gather("a", ip)
3157 :
3158 : ! We skip B calculation ONLY if the fix_displ is present and .true.
3159 2285 : if (.not. (present(fix_displ) .and. fix_displ)) then
3160 316 : call self%calc_b_from_a(ip)
3161 : end if
3162 :
3163 : ! Mixing the previous & current vectors of vibrational coefficients
3164 2285 : if (self%has_prev_grad(ip)) then
3165 1931417 : self%my_b(:,:, ip) = (one - alpha)*self%my_b(:,:, ip) + alpha*self%my_prev_b(:,:)
3166 : end if
3167 1947081 : self%my_prev_b(:,:) = self%my_b(:,:, ip)
3168 :
3169 : ! Electronic term
3170 2285 : self%enterms(1, ip) = self%get_enel(ip)
3171 : ! Vibrational term
3172 2285 : self%enterms(2, ip) = self%get_enph(ip)
3173 : ! Electron-phonon term
3174 2285 : self%enterms(3, ip) = self%get_enelph(ip)
3175 : ! Polaron energy level
3176 2285 : self%enterms(4, ip) = self%enterms(1, ip) + self%enterms(3, ip)
3177 :
3178 2285 : end subroutine polstate_localize
3179 : !!***
3180 :
3181 : !----------------------------------------------------------------------
3182 :
3183 : !!****f* m_varpeq/polstate_get_enelph
3184 : !! NAME
3185 : !! polstate_get_enelph
3186 : !!
3187 : !! FUNCTION
3188 : !! Returns the electron-phonon term of the polaron bidning energy at a specified
3189 : !! polaronic state.
3190 : !!
3191 : !! INPUTS
3192 : !! ip=Index of a polaronic state.
3193 : !!
3194 : !! OUTPUT
3195 : !! enelph=Electron-phonon term of the polaron binding energy.
3196 : !!
3197 : !! SOURCE
3198 :
3199 2285 : real(dp) function polstate_get_enelph(self, ip) result(enelph)
3200 :
3201 : !Arguments ------------------------------------
3202 : class(polstate_t), intent(inout) :: self
3203 : integer, intent(in) :: ip
3204 :
3205 : !Local variables-------------------------------
3206 : !scalars
3207 : class(gqk_t), pointer :: gqk
3208 : logical :: q_gamma
3209 : integer :: ierr, my_iq, my_pert, my_ik, ik_forw, ib, jb
3210 : complex(dp) :: a_from, a_forw, g_forw, g0, b
3211 : !arrays
3212 : real(dp) :: kpt(3), qpt(3), kpq(3)
3213 4570 : complex(dp) :: ak(self%gqk%nb_k), akq(self%gqk%nb_kq), bq(self%gqk%my_npert)
3214 : !----------------------------------------------------------------------
3215 :
3216 2285 : gqk => self%gqk
3217 :
3218 2285 : enelph = zero
3219 265381 : do my_ik=1,gqk%my_nk
3220 1052384 : kpt(:) = self%my_kpts(:, my_ik)
3221 1047222 : ak(:) = self%my_a(:, my_ik, ip)
3222 :
3223 33666383 : do my_iq=1,gqk%my_nq
3224 133604008 : qpt(:) = self%my_qpts(:, my_iq)
3225 :
3226 : ! Find k+q-->k' index in krank_kpts
3227 133604008 : kpq(:) = kpt(:) + qpt(:)
3228 33401002 : ik_forw = self%krank_kpts%get_index(kpq)
3229 : ! If erange filter was used in gstore, some transitions are not valid
3230 33401002 : if (ik_forw == -1) cycle
3231 :
3232 : ! Check if q=\Gamma
3233 40097230 : q_gamma = .false.
3234 40097230 : if (all(abs(qpt) < tol6)) q_gamma = .true.
3235 :
3236 127899998 : akq(:) = self%a_glob(:, ik_forw, ip)
3237 224628978 : bq(:) = self%my_b(:, my_iq, ip)
3238 :
3239 128163094 : do ib=1,gqk%nb_k
3240 95810144 : a_from = ak(ib)
3241 :
3242 416182160 : do jb=1,gqk%nb_kq
3243 286971014 : a_forw = akq(jb)
3244 :
3245 2104607242 : do my_pert=1,gqk%my_npert
3246 1721826084 : b = bq(my_pert)
3247 :
3248 1721826084 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
3249 : ! Add long-range correction to matrix elements at Gamma
3250 1721826084 : g0 = self%my_g0(my_pert)
3251 1721826084 : if (q_gamma .and. (ib == jb)) then
3252 4704756 : g_forw = g_forw + g0
3253 : endif
3254 :
3255 2008797098 : enelph = enelph + real(a_from*conjg(b)*g_forw*conjg(a_forw), dp)
3256 : enddo
3257 : enddo
3258 : enddo
3259 :
3260 : enddo
3261 : enddo
3262 2285 : call xmpi_sum(enelph, gqk%comm%value, ierr)
3263 2285 : enelph = -two*enelph/(one*self%nkbz*self%nqbz)
3264 :
3265 2285 : end function polstate_get_enelph
3266 : !!***
3267 :
3268 : !----------------------------------------------------------------------
3269 :
3270 : !!****f* m_varpeq/polstate_get_enph
3271 : !! NAME
3272 : !! polstate_get_enph
3273 : !!
3274 : !! FUNCTION
3275 : !! Returns the vibrational term of the polaron bidning energy at a specified
3276 : !! polaronic state.
3277 : !!
3278 : !! INPUTS
3279 : !! ip=Index of a polaronic state.
3280 : !!
3281 : !! OUTPUT
3282 : !! enph=Vibrational term of the polaron binding energy.
3283 : !!
3284 : !! SOURCE
3285 :
3286 4570 : real(dp) function polstate_get_enph(self, ip) result(enph)
3287 :
3288 : !Arguments ------------------------------------
3289 : class(polstate_t), intent(inout) :: self
3290 : integer, intent(in) :: ip
3291 :
3292 : !Local variables-------------------------------
3293 : class(gqk_t), pointer :: gqk
3294 : integer :: ierr, my_iq, my_pert
3295 : !----------------------------------------------------------------------
3296 :
3297 2285 : gqk => self%gqk
3298 :
3299 2285 : enph = zero
3300 280113 : do my_iq=1,gqk%my_nq
3301 1947081 : do my_pert=1,gqk%my_npert
3302 : enph = enph + &
3303 1944796 : gqk%my_wnuq(my_pert, my_iq)*abs(self%my_b(my_pert, my_iq, ip))**2
3304 : enddo
3305 : enddo
3306 2285 : call xmpi_sum(enph, gqk%qpt_pert_comm%value, ierr)
3307 2285 : enph = enph/real(self%nqbz, dp)
3308 :
3309 2285 : end function polstate_get_enph
3310 : !!***
3311 :
3312 : !----------------------------------------------------------------------
3313 :
3314 : !!****f* m_varpeq/polstate_get_enel
3315 : !! NAME
3316 : !! polstate_get_enel
3317 : !!
3318 : !! FUNCTION
3319 : !! Returns the electronic term of the polaron bidning energy at a specified
3320 : !! polaronic state.
3321 : !!
3322 : !! INPUTS
3323 : !! ip=Index of a polaronic state.
3324 : !!
3325 : !! OUTPUT
3326 : !! enel=Electronic term of the polaron binding energy.
3327 : !!
3328 : !! SOURCE
3329 :
3330 4570 : real(dp) function polstate_get_enel(self, ip) result(enel)
3331 :
3332 : !Arguments ------------------------------------
3333 : class(polstate_t), intent(inout) :: self
3334 : integer, intent(in) :: ip
3335 :
3336 : !Local variables-------------------------------
3337 : class(gqk_t), pointer :: gqk
3338 : integer :: ierr, my_ik, ik_ibz, ib
3339 :
3340 : !----------------------------------------------------------------------
3341 :
3342 2285 : gqk => self%gqk
3343 :
3344 2285 : enel = zero
3345 265381 : do my_ik=1,gqk%my_nk
3346 263096 : ik_ibz = gqk%my_k2ibz(1, my_ik)
3347 1049507 : do ib=1,gqk%nb_k
3348 1047222 : enel = enel + self%eig(ib, ik_ibz)*abs(self%my_a(ib, my_ik, ip))**2
3349 : enddo
3350 : enddo
3351 2285 : call xmpi_sum(enel, gqk%kpt_comm%value, ierr)
3352 2285 : enel = enel/self%nkbz
3353 :
3354 2285 : end function polstate_get_enel
3355 : !!***
3356 :
3357 : !----------------------------------------------------------------------
3358 :
3359 : !!****f* m_varpeq/polstate_calc_b_from_a
3360 : !! NAME
3361 : !! polstate_calc_b_from_a
3362 : !!
3363 : !! FUNCTION
3364 : !! Calculate vibrational coefficients B_q\nu at the current configuration,
3365 : !! so the energy gradient wrt B_q\nu is 0.
3366 : !!
3367 : !! INPUTS
3368 : !! ip=Index of a polaronic state.
3369 : !!
3370 : !! OUTPUT
3371 : !!
3372 : !! SOURCE
3373 :
3374 316 : subroutine polstate_calc_b_from_a(self, ip)
3375 :
3376 : !Arguments ------------------------------------
3377 : class(polstate_t), intent(inout) :: self
3378 : integer, intent(in) :: ip
3379 :
3380 : !Local variables-------------------------------
3381 : !scalars
3382 : class(gqk_t), pointer :: gqk
3383 : logical :: q_gamma
3384 : integer :: ierr, my_iq, my_pert, my_ik, ik_forw, ib, jb, ik_ibz, pert_glob
3385 : real(dp) :: wqnu
3386 : complex(dp) :: a_from, a_forw, g_forw, g0, b_tmp
3387 : !arrays
3388 : real(dp) :: qpt(3), kpq(3)
3389 632 : complex(dp) :: ak(self%gqk%nb_k), akq(self%gqk%nb_kq)
3390 : !----------------------------------------------------------------------
3391 :
3392 316 : gqk => self%gqk
3393 :
3394 32019 : do my_iq=1,gqk%my_nq
3395 126812 : qpt(:) = self%my_qpts(:, my_iq)
3396 : ! Check if q=\Gamma
3397 39301 : q_gamma = .false.
3398 39301 : if (all(abs(qpt) < tol6)) q_gamma = .true.
3399 :
3400 222237 : do my_pert=1,gqk%my_npert
3401 190218 : pert_glob = gqk%my_pert_start + my_pert - 1
3402 :
3403 190218 : wqnu = gqk%my_wnuq(my_pert, my_iq)
3404 :
3405 : ! Skip acoustic modes at Gamma
3406 190218 : if (wqnu < EPHTK_WTOL) then
3407 948 : self%my_b(my_pert, my_iq, ip) = zero
3408 948 : cycle
3409 : endif
3410 189270 : g0 = self%my_g0(my_pert)
3411 :
3412 : ! For this q and perturbation, calculate linear B_q\nu sum
3413 189270 : b_tmp = zero
3414 15950619 : do my_ik=1,self%gqk%my_nk
3415 : ! Find k+q-->k' index in krank_kpts
3416 63045396 : kpq(:) = qpt(:) + self%my_kpts(:, my_ik)
3417 15761349 : ik_forw = self%krank_kpts%get_index(kpq)
3418 : ! If erange filter was used in gstore, some transitions are not valid
3419 15761349 : if (ik_forw == -1) cycle
3420 :
3421 7894461 : ik_ibz = gqk%my_k2ibz(1, my_ik)
3422 :
3423 28836822 : ak(:) = self%my_a(:, my_ik, ip)
3424 28836822 : akq(:) = self%a_glob(:, ik_forw, ip)
3425 :
3426 29026092 : do ib=1,gqk%nb_k
3427 20942361 : a_from = ak(ib)
3428 :
3429 96789771 : do jb=1,gqk%nb_kq
3430 60086061 : a_forw = akq(jb)
3431 :
3432 60086061 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
3433 : ! Add long-range correction to matrix elements at Gamma
3434 60086061 : if (q_gamma .and. (ib == jb)) then
3435 137253 : g_forw = g_forw + g0
3436 : endif
3437 :
3438 81028422 : b_tmp = b_tmp + a_from*g_forw*conjg(a_forw)
3439 :
3440 : enddo
3441 : enddo
3442 : enddo
3443 189270 : call xmpi_sum(b_tmp, gqk%kpt_comm%value, ierr)
3444 :
3445 189270 : b_tmp = b_tmp / (real(self%nkbz, dp) * wqnu)
3446 :
3447 220973 : self%my_b(my_pert, my_iq, ip) = b_tmp
3448 : enddo
3449 : enddo
3450 :
3451 316 : end subroutine polstate_calc_b_from_a
3452 : !!***
3453 :
3454 : !----------------------------------------------------------------------
3455 :
3456 : !!****f* m_varpeq/polstate_calc_b_from_displ
3457 : !! NAME
3458 : !! polstate_calc_b_from_displ
3459 : !!
3460 : !! FUNCTION
3461 : !! Calculate vibrational coefficients B_q\nu from a set of atomic
3462 : !! displacements.
3463 : !!
3464 : !! INPUTS
3465 : !! ip=Index of a polaronic state.
3466 : !! scell<supercell_type>=Supercell data structure.
3467 : !! cryst<crystal_type>=Crystal data structure.
3468 : !!
3469 : !! OUTPUT
3470 : !!
3471 : !! SOURCE
3472 :
3473 0 : subroutine polstate_calc_b_from_displ(self, ip, scell, cryst)
3474 :
3475 : !Arguments ------------------------------------
3476 : class(polstate_t), intent(inout) :: self
3477 : integer, intent(in) :: ip
3478 : class(supercell_type), intent(in) :: scell
3479 : class(crystal_t), intent(in) :: cryst
3480 :
3481 : !Local variables-------------------------------
3482 : !scalars
3483 : class(gqk_t), pointer :: gqk
3484 : integer :: iatom_sc, iatom_uc, typat, idir
3485 : integer :: my_iq, my_pert
3486 : real(dp) :: wqnu, phdispl_re, phdispl_im, amass, dtau
3487 : complex(dp) :: bqnu_tmp, my_phdispl, cphase
3488 : !arrays
3489 : integer :: cell_vec(3)
3490 : real(dp) :: qpt(3)
3491 :
3492 : !----------------------------------------------------------------------
3493 :
3494 0 : gqk => self%gqk
3495 :
3496 : ! loop over q-vectors and perturbations
3497 0 : do my_iq=1,gqk%my_nq
3498 0 : qpt(:) = self%my_qpts(:, my_iq)
3499 :
3500 0 : do my_pert=1,gqk%my_npert
3501 0 : wqnu = gqk%my_wnuq(my_pert, my_iq)
3502 :
3503 0 : bqnu_tmp = zero
3504 : ! loop over ALL atoms in the supercell & cartesian directions
3505 0 : do iatom_sc=1,scell%natom
3506 0 : iatom_uc = scell%atom_indexing(iatom_sc)
3507 0 : typat = scell%typat(iatom_sc)
3508 0 : cell_vec(:) = scell%uc_indexing(:, iatom_sc)
3509 :
3510 0 : cphase = exp(j_dpc * two_pi * dot_product(qpt, cell_vec))
3511 0 : amass = cryst%amu(typat) * amu_emass
3512 :
3513 0 : do idir=1,3
3514 0 : phdispl_re = gqk%my_displ_cart(1, idir, iatom_uc, my_pert, my_iq)
3515 0 : phdispl_im = gqk%my_displ_cart(2, idir, iatom_uc, my_pert, my_iq)
3516 0 : my_phdispl = phdispl_re + j_dpc * phdispl_im
3517 :
3518 0 : dtau = self%displ(idir, iatom_sc, ip)
3519 :
3520 0 : bqnu_tmp = bqnu_tmp + amass * my_phdispl * cphase * dtau
3521 :
3522 : enddo
3523 : enddo
3524 0 : bqnu_tmp = -bqnu_tmp * sqrt(wqnu) * sqrthalf
3525 :
3526 0 : self%my_b(my_pert, my_iq, ip) = bqnu_tmp
3527 :
3528 : enddo
3529 : enddo
3530 :
3531 0 : end subroutine polstate_calc_b_from_displ
3532 : !!***
3533 :
3534 : !----------------------------------------------------------------------
3535 :
3536 : !!****f* m_varpeq/polstate_calc_displ_from_b
3537 : !! NAME
3538 : !! polstate_calc_displ_from_b
3539 : !!
3540 : !! FUNCTION
3541 : !! Calculate atomic displacements from a set of vibrational
3542 : !! coefficients B_q\nu.
3543 : !!
3544 : !! INPUTS
3545 : !! ip=Index of a polaronic state.
3546 : !! scell<supercell_type>=Supercell data structure.
3547 : !!
3548 : !! OUTPUT
3549 : !!
3550 : !! SOURCE
3551 :
3552 0 : subroutine polstate_calc_displ_from_b(self, ip, scell)
3553 :
3554 : !Arguments ------------------------------------
3555 : class(polstate_t), intent(inout) :: self
3556 : integer, intent(in) :: ip
3557 : class(supercell_type), intent(in) :: scell
3558 :
3559 : !Local variables-------------------------------
3560 : !scalars
3561 : class(gqk_t), pointer :: gqk
3562 : integer :: ierr
3563 : integer :: iatom_sc, iatom_uc, idir
3564 : integer :: my_iq, my_pert
3565 : real(dp) :: wqnu, phdispl_re, phdispl_im
3566 : complex(dp) :: bqnu, my_phdispl, cphase
3567 : !arrays
3568 : integer :: cell_vec(3)
3569 : real(dp) :: qpt(3)
3570 :
3571 : !----------------------------------------------------------------------
3572 :
3573 0 : gqk => self%gqk
3574 :
3575 0 : self%displ(:,:,ip) = zero
3576 :
3577 : ! loop over q-vectors and perturbations
3578 0 : do my_iq=1,gqk%my_nq
3579 0 : qpt(:) = self%my_qpts(:, my_iq)
3580 :
3581 0 : do my_pert=1,gqk%my_npert
3582 0 : wqnu = gqk%my_wnuq(my_pert, my_iq)
3583 :
3584 : ! Skip acoustic modes at Gamma
3585 0 : if (abs(wqnu) < tol12) cycle
3586 :
3587 0 : bqnu = self%my_b(my_pert, my_iq, ip)
3588 :
3589 : ! loop over ALL atoms in the supercell & cartesian directions
3590 0 : do iatom_sc=1,scell%natom
3591 0 : iatom_uc = scell%atom_indexing(iatom_sc)
3592 0 : cell_vec(:) = scell%uc_indexing(:, iatom_sc)
3593 :
3594 0 : cphase = exp(j_dpc * two_pi * dot_product(qpt, cell_vec))
3595 :
3596 0 : do idir=1,3
3597 0 : phdispl_re = gqk%my_displ_cart(1, idir, iatom_uc, my_pert, my_iq)
3598 0 : phdispl_im = gqk%my_displ_cart(2, idir, iatom_uc, my_pert, my_iq)
3599 0 : my_phdispl = phdispl_re + j_dpc * phdispl_im
3600 :
3601 : self%displ(idir, iatom_sc, ip) = self%displ(idir, iatom_sc, ip) + &
3602 0 : conjg(bqnu) * my_phdispl * cphase / sqrt(wqnu)
3603 :
3604 : enddo
3605 : enddo
3606 : enddo
3607 : enddo
3608 :
3609 0 : call xmpi_sum(self%displ(:,:,ip), gqk%qpt_pert_comm%value, ierr)
3610 0 : self%displ(:,:,ip) = -sqrt2 / self%nqbz * self%displ(:,:,ip)
3611 :
3612 0 : end subroutine polstate_calc_displ_from_b
3613 : !!***
3614 :
3615 : !----------------------------------------------------------------------
3616 :
3617 : !!****f* m_varpeq/calc_clb_displ
3618 : !! NAME
3619 : !! calc_clb_displ
3620 : !!
3621 : !! FUNCTION
3622 : !! Calculate Coulomb displacements introduced by a unit charge localized
3623 : !! on a supercell atom
3624 : !! Note that internally it is assumed that the localized charge is negative,
3625 : !! regargdless of the polaron type. However, for hole polaron, the displacements
3626 : !! are filpped at post-processing and hence physically correct.
3627 : !!
3628 : !! INPUTS
3629 : !! atloc=Index of an atom where charge is localized.
3630 : !! chrgat(natom)=Charge of the atoms.
3631 : !! scell<supercell_type>=Supercell data structure.
3632 : !! cryst<crystal_t>=Crystal data structure.
3633 : !!
3634 : !! OUTPUT
3635 : !!
3636 : !! SOURCE
3637 :
3638 0 : subroutine polstate_calc_clb_displ(self, atloc, chrgat, scell, cryst)
3639 :
3640 : !Arguments ------------------------------------
3641 : class(polstate_t), intent(inout) :: self
3642 : integer, intent(in) :: atloc
3643 : class(supercell_type), intent(in) :: scell
3644 : class(crystal_t), intent(in) :: cryst
3645 : real(dp), intent(in) :: chrgat(scell%natom_primcell)
3646 :
3647 : !Local variables-------------------------------
3648 : !scalars
3649 : integer :: ip, ii, ix, iy, iz, iatom_sc, iatom_uc
3650 : !arrays
3651 : integer :: sc_ind(3)
3652 : integer :: box_dim(3), center_sc(3), vec(3)
3653 : real(dp) :: rdist(3), center_atom(3)
3654 0 : integer, allocatable :: scell_flag(:,:,:)
3655 0 : real(dp), allocatable :: scell_displ(:,:,:,:,:)
3656 :
3657 : !----------------------------------------------------------------------
3658 :
3659 : ! supercell center
3660 0 : center_sc(:) = (self%ngkpt(:) + 1) / 2
3661 : ! select atom in the unit cell
3662 0 : center_atom(:) = cryst%xcart(:, atloc)
3663 :
3664 : ! displacement box dimensions
3665 0 : do ii=1,3
3666 0 : box_dim(ii) = min(3, self%ngkpt(ii))
3667 : enddo
3668 :
3669 0 : ABI_MALLOC(scell_displ, (3, cryst%natom, self%ngkpt(1), self%ngkpt(2), self%ngkpt(3)))
3670 0 : ABI_MALLOC(scell_flag, (self%ngkpt(1), self%ngkpt(2), self%ngkpt(3)))
3671 0 : scell_flag(:,:,:) = zero
3672 0 : scell_displ(:,:,:,:,:) = zero
3673 :
3674 0 : do iz=1,box_dim(3)
3675 0 : do iy=1,box_dim(2)
3676 0 : do ix=1,box_dim(1)
3677 : ! get all displacement vectors wrt centeral supercell
3678 : vec(:) = [shift_(ix-1, box_dim(1)), &
3679 0 : shift_(iy-1, box_dim(2)), shift_(iz-1, box_dim(3))]
3680 :
3681 : ! get indices of related supercells
3682 0 : do ii=1,3
3683 : ! add 1 becasue of fortran indexing
3684 0 : sc_ind(ii) = mod(center_sc(ii) + vec(ii) - 1, self%ngkpt(ii))
3685 0 : sc_ind(ii) = sc_ind(ii) + 1
3686 : enddo
3687 0 : scell_flag(sc_ind(1), sc_ind(2), sc_ind(3)) = one
3688 :
3689 : ! calculate displacements in these supercells
3690 0 : do iatom_uc=1,cryst%natom
3691 : rdist(:) = matmul(cryst%rprimd(:,:), vec(:)) + &
3692 0 : cryst%xcart(:,iatom_uc) - center_atom(:)
3693 0 : if (norm2(rdist) > tol8) then
3694 : scell_displ(:,iatom_uc,sc_ind(1),sc_ind(2),sc_ind(3)) = &
3695 0 : -chrgat(iatom_uc) * rdist(:) / norm2(rdist)**3
3696 : endif
3697 :
3698 : enddo
3699 :
3700 : enddo
3701 : enddo
3702 : enddo
3703 :
3704 0 : self%displ(:,:,:) = zero
3705 : ! insert displacements
3706 0 : ip = 1
3707 0 : do iatom_sc=1,scell%natom
3708 0 : iatom_uc = scell%atom_indexing(iatom_sc)
3709 : ! add 1 becasue of fortran indexing
3710 0 : vec(:) = scell%uc_indexing(:, iatom_sc) + 1
3711 0 : if (scell_flag(vec(1), vec(2), vec(3)) /= 0) then
3712 0 : self%displ(:,iatom_sc,ip) = scell_displ(:,iatom_uc, vec(1), vec(2), vec(3))
3713 : endif
3714 : enddo
3715 0 : do ip=2,self%np
3716 0 : self%displ(:,iatom_sc,ip) = self%displ(:,iatom_sc,ip-1)
3717 : enddo
3718 :
3719 0 : ABI_FREE(scell_flag)
3720 0 : ABI_FREE(scell_displ)
3721 :
3722 : !----------------------------------------------------------------------
3723 :
3724 : contains
3725 0 : integer function shift_(x, n) result(x_shift)
3726 : integer, intent(in) :: x, n
3727 0 : x_shift = x
3728 0 : if (x > n/2) x_shift = -n + x
3729 : end function shift_
3730 :
3731 : end subroutine polstate_calc_clb_displ
3732 : !!***
3733 :
3734 : !----------------------------------------------------------------------
3735 :
3736 : !!****f* m_varpeq/polstate_calc_phgrad
3737 : !! NAME
3738 : !! polstate_calc_phgrad
3739 : !!
3740 : !! FUNCTION
3741 : !! Calculate the gradient of the polaronic energy wrt phonon
3742 : !! coefficients B_q\nu.
3743 : !!
3744 : !! INPUTS
3745 : !! ip=Index of a polaronic state.
3746 : !!
3747 : !! OUTPUT
3748 : !!
3749 : !! SOURCE
3750 :
3751 135 : subroutine polstate_calc_phgrad(self, ip)
3752 :
3753 : !Arguments ------------------------------------
3754 : class(polstate_t), intent(inout) :: self
3755 : integer, intent(in) :: ip
3756 :
3757 : !Local variables-------------------------------
3758 : !scalars
3759 : class(gqk_t), pointer :: gqk
3760 : logical :: q_gamma
3761 : integer :: ierr
3762 : integer :: my_iq, my_pert
3763 : integer:: my_ik, ik_forw, ib, jb
3764 : real(dp) :: wqnu
3765 : complex(dp) :: a_from, a_forw, g_forw, phgrad_tmp
3766 : complex(dp) :: g0, b
3767 : !arrays
3768 : real(dp) :: qpt(3), kpq(3)
3769 270 : complex(dp) :: ak(self%gqk%nb_k), akq(self%gqk%nb_kq)
3770 :
3771 : !----------------------------------------------------------------------
3772 :
3773 135 : gqk => self%gqk
3774 :
3775 17010 : do my_iq=1,gqk%my_nq
3776 67500 : qpt(:) = self%my_qpts(:, my_iq)
3777 : ! Check if q=\Gamma
3778 21060 : q_gamma = .false.
3779 21060 : if (all(abs(qpt) < tol6)) q_gamma = .true.
3780 :
3781 118260 : do my_pert=1,gqk%my_npert
3782 101250 : b = self%my_b(my_pert, my_iq, ip)
3783 101250 : wqnu = gqk%my_wnuq(my_pert, my_iq)
3784 101250 : g0 = self%my_g0(my_pert)
3785 :
3786 : ! For this q and perturbation, calculate B_q\nu sum
3787 101250 : phgrad_tmp = zero
3788 :
3789 12757500 : do my_ik=1,self%gqk%my_nk
3790 : ! Find k+q-->k' index in krank_kpts
3791 50625000 : kpq(:) = qpt(:) + self%my_kpts(:, my_ik)
3792 12656250 : ik_forw = self%krank_kpts%get_index(kpq)
3793 :
3794 : ! If erange filter was used in gstore, some transitions are not valid
3795 12656250 : if (ik_forw == -1) cycle
3796 :
3797 50625000 : ak(:) = self%my_a(:, my_ik, ip)
3798 50625000 : akq(:) = self%a_glob(:, ik_forw, ip)
3799 :
3800 50726250 : do ib=1,gqk%nb_k
3801 37968750 : a_from = ak(ib)
3802 :
3803 164531250 : do jb=1,gqk%nb_kq
3804 113906250 : a_forw = akq(jb)
3805 :
3806 113906250 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
3807 : ! Add long-range correction to matrix elements at Gamma
3808 113906250 : if (q_gamma .and. (ib == jb)) then
3809 303750 : g_forw = g_forw + g0
3810 : endif
3811 :
3812 151875000 : phgrad_tmp = phgrad_tmp + a_from*g_forw*conjg(a_forw)
3813 : enddo
3814 : enddo
3815 : enddo
3816 101250 : call xmpi_sum(phgrad_tmp, gqk%kpt_comm%value, ierr)
3817 :
3818 : self%my_phgrad(my_pert, my_iq, ip) = &
3819 118125 : two/self%nqbz * (b*wqnu - phgrad_tmp / self%nkbz)
3820 :
3821 : enddo
3822 : enddo
3823 :
3824 135 : end subroutine polstate_calc_phgrad
3825 : !!***
3826 :
3827 : !----------------------------------------------------------------------
3828 :
3829 : !!****f* m_varpeq/polstate_calc_hpol
3830 : !! NAME
3831 : !! polstate_calc_hpol
3832 : !!
3833 : !! FUNCTION
3834 : !! Construct and calculate the polaron Hamiltonian, <A_bra|H_pol(B)|A_ket>
3835 : !!
3836 : !! INPUTS
3837 : !! my_a_bra(gqk%nb_kq,gqk%my_nk)=Electronic coefficients for the bra vector
3838 : !! a_bra_lglob(gqk%nb_kq,gqk%glob_nk)=Electronic coefficients for the bra vectors
3839 : !! (global array)
3840 : !! my_a_ket(gqk%nb_k,gqk%my_nk)=Electronic coefficients for the ket vector
3841 : !! my_b(gqk%my_npert,gqk%my_nq)=Phonon coeffcients used to construct H_pol(B)
3842 : !! Note that these arrays are MPI-distributed as in the polstate datatype
3843 : !!
3844 : !! OUTPUT
3845 : !!
3846 : !! SOURCE
3847 :
3848 6450 : real(dp) function polstate_calc_hpol(self, my_a_bra, a_bra_glob, my_a_ket, my_b) &
3849 : result(hpol)
3850 :
3851 : !Arguments ------------------------------------
3852 : class(polstate_t), intent(inout) :: self
3853 : complex(dp), intent(in) :: my_a_bra(self%gqk%nb_kq, self%gqk%my_nk)
3854 : complex(dp), intent(in) :: a_bra_glob(self%gqk%nb_kq, self%gqk%glob_nk)
3855 : complex(dp), intent(in) :: my_a_ket(self%gqk%nb_k, self%gqk%my_nk)
3856 : complex(dp), intent(in) :: my_b(self%gqk%my_npert, self%gqk%my_nq)
3857 :
3858 : !Local variables-------------------------------
3859 : !scalars
3860 : class(gqk_t), pointer :: gqk
3861 : logical :: q_gamma
3862 : integer :: ierr, my_ik, ik_ibz, ib
3863 : integer :: my_iq, my_pert, ik_forw, jb
3864 : complex(dp) :: a_from, a_forw, g_forw, g0, b
3865 : real(dp) :: h_el, h_elph
3866 : !arrays
3867 : real(dp) :: kpt(3), qpt(3), kpq(3)
3868 12900 : complex(dp) :: a_ket(self%gqk%nb_k), a_bra(self%gqk%nb_kq), bq(self%gqk%my_npert)
3869 :
3870 : !----------------------------------------------------------------------
3871 :
3872 6450 : gqk => self%gqk
3873 :
3874 : ! <A_1|H_p(B)|A_2> = <A_1|H_el(B)|A_2> + <A_1|H_el-ph(B)|A_2>
3875 :
3876 : ! 1st term
3877 6450 : h_el = zero
3878 746763 : do my_ik=1,gqk%my_nk
3879 740313 : ik_ibz = gqk%my_k2ibz(1, my_ik)
3880 2952750 : do ib=1,gqk%nb_k
3881 : h_el = h_el + &
3882 2946300 : conjg(my_a_bra(ib, my_ik))*self%eig(ib, ik_ibz)*my_a_ket(ib, my_ik)
3883 : enddo
3884 : enddo
3885 6450 : call xmpi_sum(h_el, gqk%kpt_comm%value, ierr)
3886 6450 : h_el = h_el/self%nkbz
3887 :
3888 : ! 2nd term
3889 6450 : h_elph = zero
3890 746763 : do my_ik=1,gqk%my_nk
3891 2961252 : kpt(:) = self%my_kpts(:, my_ik)
3892 2946300 : a_ket(:) = my_a_ket(:, my_ik)
3893 :
3894 94757046 : do my_iq=1,gqk%my_nq
3895 376041132 : qpt(:) = self%my_qpts(:, my_iq)
3896 :
3897 : ! Find k+q-->k' index in krank_kpts
3898 376041132 : kpq(:) = kpt(:) + qpt(:)
3899 94010283 : ik_forw = self%krank_kpts%get_index(kpq)
3900 : ! If erange filter was used in gstore, some transitions are not valid
3901 94010283 : if (ik_forw == -1) cycle
3902 :
3903 : ! Check if q=\Gamma
3904 112728480 : q_gamma = .false.
3905 112728480 : if (all(abs(qpt) < tol6)) q_gamma = .true.
3906 :
3907 359519172 : a_bra(:) = a_bra_glob(:, ik_forw)
3908 631487325 : bq(:) = my_b(:, my_iq)
3909 :
3910 360259485 : do ib=1,gqk%nb_k
3911 269306697 : a_from = a_ket(ib)
3912 :
3913 1169906343 : do jb=1,gqk%nb_kq
3914 806589363 : a_forw = a_bra(jb)
3915 :
3916 5915432238 : do my_pert=1,gqk%my_npert
3917 4839536178 : b = bq(my_pert)
3918 :
3919 4839536178 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
3920 : ! Add long-range correction to matrix elements at Gamma
3921 4839536178 : g0 = self%my_g0(my_pert)
3922 4839536178 : if (q_gamma .and. (ib == jb)) then
3923 13235922 : g_forw = g_forw + g0
3924 : endif
3925 :
3926 5646125541 : h_elph = h_elph + real(a_from*conjg(b)*g_forw*conjg(a_forw), dp)
3927 : enddo
3928 : enddo
3929 : enddo
3930 :
3931 : enddo
3932 : enddo
3933 6450 : call xmpi_sum(h_elph, gqk%comm%value, ierr)
3934 6450 : h_elph = -two*h_elph/(one*self%nkbz*self%nqbz)
3935 :
3936 6450 : hpol = h_el + h_elph
3937 :
3938 6450 : end function polstate_calc_hpol
3939 : !!***
3940 :
3941 : !----------------------------------------------------------------------
3942 :
3943 : !!****f* m_varpeq/polstate_load_a
3944 : !! NAME
3945 : !! polstate_load_a
3946 : !!
3947 : !! FUNCTION
3948 : !! Load the initial vector of electronic coefficients A_nk from source.
3949 : !!
3950 : !! INPUTS
3951 : !! a_src(self%gqk%nb_k, self%gqk%glob_nk)=Global A_nk array to be loaded.
3952 : !!
3953 : !! OUTPUT
3954 : !!
3955 : !! SOURCE
3956 :
3957 0 : subroutine polstate_load_a(self, a_src, ip)
3958 :
3959 : !Arguments ------------------------------------
3960 : !scalars
3961 : class(polstate_t), intent(inout) :: self
3962 : integer, intent(in) :: ip
3963 : !arrays
3964 : complex(dp), intent(in) :: a_src(self%gqk%nb_k, self%gqk%glob_nk)
3965 :
3966 : !Local variables-------------------------------
3967 : class(gqk_t), pointer :: gqk
3968 : integer :: my_ik, ik_glob, ib
3969 : !----------------------------------------------------------------------
3970 :
3971 0 : gqk => self%gqk
3972 :
3973 0 : do my_ik=1,gqk%my_nk
3974 0 : ik_glob = gqk%my_kstart + my_ik - 1
3975 0 : do ib=1,gqk%nb_k
3976 0 : self%my_a(ib, my_ik, ip) = a_src(ib, ik_glob)
3977 : enddo
3978 : enddo
3979 :
3980 0 : end subroutine polstate_load_a
3981 : !!***
3982 :
3983 : !----------------------------------------------------------------------
3984 :
3985 : !!****f* m_varpeq/polstate_load_b
3986 : !! NAME
3987 : !! polstate_load_b
3988 : !!
3989 : !! FUNCTION
3990 : !! Load the initial vector of phonon coefficients B_q\nu from source.
3991 : !!
3992 : !! INPUTS
3993 : !! b_src(self%gqk%natom3, self%gqk%glob_nq)=Global B_q\nu array to be loaded.
3994 : !! ip=Index of a polaronic state.
3995 : !! trvec(3) [optional]=Translational vector.
3996 : !!
3997 : !! OUTPUT
3998 : !!
3999 : !! SOURCE
4000 :
4001 2 : subroutine polstate_load_b(self, b_src, ip, trvec)
4002 :
4003 : !Arguments ------------------------------------
4004 : !scalars
4005 : class(polstate_t), intent(inout) :: self
4006 : integer, intent(in) :: ip
4007 : !arrays
4008 : complex(dp), intent(in) :: b_src(self%gqk%natom3, self%gqk%glob_nq)
4009 : integer, optional, intent(in) :: trvec(3)
4010 :
4011 : !Local variables-------------------------------
4012 : !scalars
4013 : class(gqk_t), pointer :: gqk
4014 : integer :: my_iq, my_pert, pert_glob, iq_glob
4015 : complex(dp) :: cphase
4016 : !arrays
4017 : integer trvec_(3)
4018 : real(dp) :: qpt(3)
4019 :
4020 : !----------------------------------------------------------------------
4021 :
4022 2 : gqk => self%gqk
4023 :
4024 2 : trvec_ (:) = 0
4025 2 : if (present(trvec)) trvec_(:) = trvec(:)
4026 :
4027 252 : do my_iq=1,gqk%my_nq
4028 250 : iq_glob = gqk%my_qstart + my_iq - 1
4029 1000 : qpt(:) = self%my_qpts(:, my_iq)
4030 :
4031 1000 : cphase = exp(+j_dpc * two_pi * dot_product(qpt, trvec_))
4032 :
4033 1752 : do my_pert=1,gqk%my_npert
4034 1500 : pert_glob = gqk%my_pert_start + my_pert - 1
4035 1750 : self%my_b(my_pert, my_iq, ip) = b_src(pert_glob, iq_glob) * cphase
4036 : enddo
4037 : enddo
4038 :
4039 2 : end subroutine polstate_load_b
4040 : !!***
4041 :
4042 : !----------------------------------------------------------------------
4043 :
4044 : !!****f* m_varpeq/polstate_linterp_b
4045 : !! NAME
4046 : !! polstate_linterp_b
4047 : !!
4048 : !! FUNCTION
4049 : !! Interpolate the phonon coefficients on uniform mesh between initial
4050 : !! and final state.
4051 : !!
4052 : !! INPUTS
4053 : !!
4054 : !! OUTPUT
4055 : !!
4056 : !! SOURCE
4057 :
4058 1 : subroutine polstate_linterp_b(self)
4059 :
4060 : !Arguments ------------------------------------
4061 : class(polstate_t), intent(inout) :: self
4062 :
4063 : !Local variables-------------------------------
4064 : !scalars
4065 : class(gqk_t), pointer :: gqk
4066 : integer :: ip, from_ip, to_ip
4067 : !arrays
4068 2 : complex(dp) :: my_b_step(self%gqk%my_npert, self%gqk%my_nq)
4069 :
4070 : !----------------------------------------------------------------------
4071 :
4072 1 : gqk => self%gqk
4073 :
4074 1 : from_ip = 1; to_ip = self%np
4075 :
4076 876 : my_b_step(:,:) = (self%my_b(:,:,to_ip) - self%my_b(:,:,from_ip)) / (self%np - one)
4077 :
4078 8 : do ip=2,self%np-1
4079 6133 : self%my_b(:,:,ip) = self%my_b(:,:,ip-1) + my_b_step(:,:)
4080 : enddo
4081 :
4082 1 : end subroutine polstate_linterp_b
4083 : !!***
4084 :
4085 : !----------------------------------------------------------------------
4086 :
4087 : !!****f* m_varpeq/polstate_redistr_b
4088 : !! NAME
4089 : !! polstate_redistr_b
4090 : !!
4091 : !! FUNCTION
4092 : !! Redistribute the phonon coefficients B_\qnu with spline interpolation.
4093 : !! Assumes parametrization [0,1] -> B_\qnu.
4094 : !!
4095 : !! INPUTS
4096 : !!
4097 : !! OUTPUT
4098 : !!
4099 : !! SOURCE
4100 :
4101 15 : subroutine polstate_redistr_b(self, mesh_in, mesh_out, nimag)
4102 :
4103 : !Arguments ------------------------------------
4104 : class(polstate_t), intent(inout) :: self
4105 : integer, intent(in) :: nimag
4106 : real(dp), intent(in) :: mesh_in(nimag)
4107 : real(dp), intent(in) :: mesh_out(nimag)
4108 :
4109 : !Local variables-------------------------------
4110 : !scalars
4111 : class(gqk_t), pointer :: gqk
4112 : integer :: from_ip, to_ip
4113 : integer :: my_iq, my_pert
4114 : complex(dp) :: phgrad_qnu_from, phgrad_qnu_to
4115 : !arrays
4116 30 : complex(dp) :: my_b_mesh_in(nimag)
4117 30 : complex(dp) :: my_2der_mesh_out(nimag)
4118 30 : complex(dp) :: my_b_mesh_out(nimag)
4119 :
4120 : !----------------------------------------------------------------------
4121 :
4122 15 : gqk => self%gqk
4123 :
4124 15 : from_ip = 1; to_ip = self%np
4125 :
4126 : !call self%calc_phgrad(from_ip)
4127 : !call self%calc_phgrad(to_ip)
4128 :
4129 1890 : do my_iq=1,gqk%my_nq
4130 13140 : do my_pert=1,gqk%my_npert
4131 :
4132 : ! derivatives of the binding energy at endpoints
4133 11250 : phgrad_qnu_from = self%my_phgrad(my_pert, my_iq, from_ip)
4134 11250 : phgrad_qnu_to = self%my_phgrad(my_pert, my_iq, to_ip)
4135 :
4136 112500 : my_b_mesh_in(:) = self%my_b(my_pert, my_iq, :)
4137 :
4138 : ! spline interpolation
4139 : call spline_complex(mesh_in, my_b_mesh_in, nimag, phgrad_qnu_from, &
4140 11250 : phgrad_qnu_to, my_2der_mesh_out)
4141 : call splint_complex(nimag, mesh_in, my_b_mesh_in, my_2der_mesh_out, nimag, &
4142 11250 : mesh_out, my_b_mesh_out)
4143 :
4144 114375 : self%my_b(my_pert, my_iq, :) = my_b_mesh_out(:)
4145 :
4146 : enddo
4147 : enddo
4148 :
4149 15 : end subroutine polstate_redistr_b
4150 : !!***
4151 :
4152 : !----------------------------------------------------------------------
4153 :
4154 : !!****f* m_varpeq/polstate_calc_eff_phforce
4155 : !! NAME
4156 : !! polstate_calc_eff_phforce
4157 : !!
4158 : !! FUNCTION
4159 : !! Calculate the effective phonon force for the hopping optimization.
4160 : !!
4161 : !! INPUTS
4162 : !! mode=Select the method for effective force calculation:
4163 : !! "fw_euler" ---> Forward Euler method;
4164 : !! "rk" ---> 4th order Runge-Kutta method.
4165 : !! ip=Index of a polaronic state.
4166 : !!
4167 : !! OUTPUT
4168 : !!
4169 : !! SOURCE
4170 :
4171 135 : subroutine polstate_calc_eff_phforce(self, mode, ip)
4172 :
4173 : !Arguments ------------------------------------
4174 : class(polstate_t), intent(inout) :: self
4175 : character(len=*), intent(in) :: mode
4176 : integer, intent(in) :: ip
4177 :
4178 : !Local variables-------------------------------
4179 : !scalars
4180 : class(gqk_t), pointer :: gqk
4181 : integer :: ii
4182 : !arrays
4183 : real(dp) :: rk_shifts(3) = [half, half, one]
4184 270 : complex(dp) :: my_b_copy(self%gqk%my_npert, self%gqk%my_nq)
4185 270 : complex(dp) :: rk_coeff(self%gqk%my_npert, self%gqk%my_nq, 4)
4186 :
4187 : !----------------------------------------------------------------------
4188 :
4189 135 : gqk => self%gqk
4190 :
4191 135 : call self%calc_phgrad(ip)
4192 :
4193 135 : select case(mode)
4194 : case ("fw_euler")
4195 118260 : self%my_eff_phforce(:,:,ip) = -self%my_phgrad(:,:,ip)
4196 :
4197 : case ("rk")
4198 : ! Make a copy of current B_\qnu
4199 0 : my_b_copy(:,:) = self%my_b(:,:,ip)
4200 :
4201 : ! Calculate Runge-Kutta coefficients
4202 0 : rk_coeff(:,:,1) = self%my_phgrad(:,:,ip)
4203 :
4204 0 : do ii=1,3
4205 0 : self%my_b(:,:,ip) = my_b_copy(:,:) + rk_shifts(ii) * rk_coeff(:,:,ii)
4206 0 : call self%calc_phgrad(ip)
4207 0 : rk_coeff(:,:,ii+1) = self%my_phgrad(:,:,ip)
4208 : enddo
4209 :
4210 : self%my_eff_phforce(:,:,ip) = -sixth*(rk_coeff(:,:,1) + two*rk_coeff(:,:,2) + &
4211 0 : two*rk_coeff(:,:,3) + rk_coeff(:,:,4))
4212 :
4213 : ! Restore the gradient value
4214 0 : self%my_phgrad(:,:,ip) = rk_coeff(:,:,1)
4215 :
4216 : case default
4217 135 : ABI_ERROR(sjoin("polstate_calc_eff_phforce, unsupported mode: ", mode))
4218 : end select
4219 :
4220 135 : end subroutine polstate_calc_eff_phforce
4221 : !!***
4222 :
4223 : !----------------------------------------------------------------------
4224 :
4225 : !!****f* m_varpeq/polstate_calc_hop_timestep
4226 : !! NAME
4227 : !! polstate_calc_hop_timestep
4228 : !!
4229 : !! FUNCTION
4230 : !! Calculate the optimal timestep for hopping optimization at this state.
4231 : !!
4232 : !! INPUTS
4233 : !! ip=Index of a polaronic state.
4234 : !!
4235 : !! OUTPUT
4236 : !!
4237 : !! SOURCE
4238 :
4239 7 : subroutine polstate_calc_hop_timestep(self, ip)
4240 :
4241 : !Arguments ------------------------------------
4242 : class(polstate_t), intent(inout) :: self
4243 : integer, intent(in) :: ip
4244 :
4245 : !Local variables-------------------------------
4246 : !scalars
4247 : class(gqk_t), pointer :: gqk
4248 : logical :: q_gamma
4249 : integer :: ierr
4250 : integer :: my_iq, my_pert
4251 : integer :: my_ik, ik_forw, ib, jb
4252 : real(dp) :: wqnu
4253 : real(dp) :: ts_denom, ts_num_ph, ts_num_elph
4254 : complex(dp) :: a_from, a_forw
4255 : complex(dp) :: g_forw, g0, phforce, b
4256 : !arrays
4257 : real(dp) :: kpt(3), qpt(3), kpq(3)
4258 14 : complex(dp) :: ak(self%gqk%nb_k), akq(self%gqk%nb_kq)
4259 14 : complex(dp) :: phforce_q(self%gqk%my_npert)
4260 :
4261 : !----------------------------------------------------------------------
4262 :
4263 7 : gqk => self%gqk
4264 :
4265 7 : ts_denom = zero
4266 7 : ts_num_ph = zero ! Phonon-dependent component of the numerator
4267 :
4268 : ! Phonon-dependent components
4269 882 : do my_iq=1,gqk%my_nq
4270 6132 : do my_pert=1,gqk%my_npert
4271 5250 : wqnu = gqk%my_wnuq(my_pert, my_iq)
4272 5250 : phforce = self%my_eff_phforce(my_pert, my_iq, ip)
4273 5250 : b = self%my_b(my_pert, my_iq, ip)
4274 :
4275 5250 : ts_denom = ts_denom + wqnu*abs(phforce)**2
4276 6125 : ts_num_ph = ts_num_ph + wqnu*real(b*phforce, dp)
4277 : enddo
4278 : enddo
4279 7 : call xmpi_sum(ts_denom, gqk%qpt_pert_comm%value, ierr)
4280 7 : call xmpi_sum(ts_num_ph, gqk%qpt_pert_comm%value, ierr)
4281 :
4282 7 : ts_num_elph = zero
4283 : ! Scattering-dependent component of the numerator
4284 882 : do my_iq=1,gqk%my_nq
4285 3500 : qpt(:) = self%my_qpts(:, my_iq)
4286 6125 : phforce_q(:) = self%my_eff_phforce(:, my_iq, ip)
4287 :
4288 : ! Check if q=\Gamma
4289 1092 : q_gamma = .false.
4290 1092 : if (all(abs(qpt) < tol6)) q_gamma = .true.
4291 :
4292 110257 : do my_ik=1,gqk%my_nk
4293 437500 : kpt(:) = self%my_kpts(:, my_ik)
4294 :
4295 : ! Find k+q-->k' index in krank_kpts
4296 437500 : kpq(:) = kpt(:) + qpt(:)
4297 109375 : ik_forw = self%krank_kpts%get_index(kpq)
4298 : ! If erange filter was used in gstore, some transitions are not valid
4299 109375 : if (ik_forw == -1) cycle
4300 :
4301 437500 : ak(:) = self%my_a(:, my_ik, ip)
4302 437500 : akq(:) = self%a_glob(:, ik_forw, ip)
4303 :
4304 438375 : do ib=1,gqk%nb_k
4305 328125 : a_from = ak(ib)
4306 :
4307 1421875 : do jb=1,gqk%nb_kq
4308 984375 : a_forw = akq(jb)
4309 :
4310 7218750 : do my_pert=1,gqk%my_npert
4311 5906250 : phforce = phforce_q(my_pert)
4312 5906250 : g_forw = gqk%my_g(my_pert, jb, my_iq, ib, my_ik)
4313 :
4314 : ! Add long-range correction to matrix elements at Gamma
4315 5906250 : g0 = self%my_g0(my_pert)
4316 5906250 : if (q_gamma .and. (ib == jb)) then
4317 15750 : g_forw = g_forw + g0
4318 : endif
4319 :
4320 : ts_num_elph = &
4321 6890625 : ts_num_elph + real(a_from*g_forw*conjg(phforce)*conjg(a_forw), dp)
4322 :
4323 : enddo
4324 : enddo
4325 : enddo
4326 : enddo
4327 : enddo
4328 7 : call xmpi_sum(ts_num_elph, gqk%comm%value, ierr)
4329 7 : ts_num_elph = ts_num_elph / self%nkbz
4330 :
4331 7 : self%hop_ts(ip) = (ts_num_elph - ts_num_ph) / ts_denom
4332 :
4333 7 : end subroutine polstate_calc_hop_timestep
4334 : !!***
4335 :
4336 : !----------------------------------------------------------------------
4337 :
4338 : !!****f* m_varpeq/polstate_seed_a
4339 : !! NAME
4340 : !! polstate_seed_a
4341 : !!
4342 : !! FUNCTION
4343 : !! Specify initial vector of electronic coefficients A_nk at runtime.
4344 : !!
4345 : !! INPUTS
4346 : !! mode=Select the initialization type. Possible options:
4347 : !! "gau_energy" ---> Gaussian shape based on the energy of electronic states;
4348 : !! "gau_length" ---> Gaussian shape based on the polaron localization length;
4349 : !! "even" ---> equal contribution from each electronic state;
4350 : !! "random" ---> random initalization.
4351 : !! ip=Index of a polaronic state.
4352 : !! atloc [optional]=Index of an atom where charge is localized.
4353 : !! chrgat(natom) [optional]=Charge of the atoms.
4354 : !! scell<supercell_type> [optional]=Supercell data structure.
4355 : !! cryst<crystal_t> [optional]=Crystal data structure.
4356 : !! nstep_loc [optional]=Number of iterations for charge localization.
4357 : !!
4358 : !! OUTPUT
4359 : !!
4360 : !! SOURCE
4361 :
4362 19 : subroutine polstate_seed_a(self, mode, ip, atloc, chrgat, scell, cryst, nstep_loc)
4363 :
4364 : !Arguments ------------------------------------
4365 : class(polstate_t), intent(inout) :: self
4366 : character(len=*), intent(in) :: mode
4367 : integer, intent(in) :: ip
4368 : integer, optional, intent(in) :: atloc
4369 : class(supercell_type), optional, intent(in) :: scell
4370 : class(crystal_t), optional, intent(in) :: cryst
4371 : real(dp), optional, intent(in) :: chrgat(:)
4372 : integer, optional, intent(in) :: nstep_loc
4373 :
4374 : !Local variables-------------------------------
4375 : class(gqk_t), pointer :: gqk
4376 : logical :: flag
4377 : integer :: ii, ierr
4378 : real(dp) :: a_sqnorm, grad_sqnorm
4379 : !----------------------------------------------------------------------
4380 :
4381 19 : gqk => self%gqk
4382 :
4383 38 : select case(mode)
4384 : case ("gau_energy")
4385 19 : call gau_energy_()
4386 : case ("gau_length")
4387 0 : call gau_length_()
4388 : case ("random")
4389 0 : call random_()
4390 : case ("localize")
4391 : flag = (present(atloc) .and. present(chrgat) .and. present(scell) &
4392 0 : .and. present(cryst) .and. present(nstep_loc))
4393 0 : ABI_CHECK(flag, "polstate_seed_a: atloc, chrgat and scell are needed to localize")
4394 :
4395 : ! enforce Coulomb displacements
4396 0 : if (ip == 1) call self%calc_clb_displ(atloc, chrgat, scell, cryst)
4397 :
4398 : ! calcualte associated deformation field
4399 0 : call self%calc_b_from_displ(ip, scell, cryst)
4400 :
4401 : ! start with random charge distribution
4402 0 : self%my_a(:,:,ip) = one + j_dpc*one
4403 :
4404 0 : if (self%efilter > zero) call self%filter("a", ip)
4405 :
4406 : ! then we do miniminzation with fixed displacements to localize charge distribution
4407 0 : a_sqnorm = self%get_sqnorm("a", ip)
4408 0 : if (a_sqnorm > tol12) then
4409 0 : self%my_a(:,:,ip) = sqrt(self%nkbz/a_sqnorm) * self%my_a(:,:,ip)
4410 : endif
4411 :
4412 0 : do ii=1,nstep_loc
4413 0 : call self%localize(ip, alpha=zero, fix_displ=.true.)
4414 0 : call self%calc_grad(ip)
4415 0 : grad_sqnorm = self%get_sqnorm("grad", ip)
4416 0 : if (sqrt(grad_sqnorm) < tol6) then
4417 : exit
4418 : endif
4419 0 : call self%update_pc(ip)
4420 0 : call self%calc_pcjgrad(ip,ii,nstep_ort=0)
4421 0 : call self%update_a(ip)
4422 : enddo
4423 :
4424 : case ("even")
4425 0 : self%my_a(:,:,ip) = one + j_dpc*one
4426 : case default
4427 19 : ABI_ERROR(sjoin("polstate_seed_a, unsupported mode: ", mode))
4428 : end select
4429 :
4430 : !----------------------------------------------------------------------
4431 :
4432 : contains
4433 19 : subroutine gau_energy_()
4434 : real(dp) :: mu, sigma
4435 : integer :: my_ik, ik_ibz, ib
4436 : real(dp) :: eig
4437 19 : mu = self%gpr_energy(1); sigma = self%gpr_energy(2)
4438 19 : ABI_CHECK(sigma /= 0, "gpr_energy: standard deviation must be non-zero")
4439 1746 : do my_ik=1,gqk%my_nk
4440 1727 : ik_ibz = gqk%my_k2ibz(1, my_ik)
4441 6749 : do ib=1,gqk%nb_k
4442 5003 : eig = self%eig(ib, ik_ibz)
4443 6730 : self%my_a(ib, my_ik, ip) = exp(-half*(eig - mu)**2/sigma**2)
4444 : enddo
4445 : enddo
4446 19 : end subroutine gau_energy_
4447 :
4448 0 : subroutine gau_length_()
4449 : integer :: my_ik, ib
4450 : real(dp) :: kpt(3)
4451 0 : do my_ik=1,gqk%my_nk
4452 0 : kpt(:) = self%my_kpts(:, my_ik)
4453 0 : do ib=1,gqk%nb_k
4454 0 : self%my_a(ib, my_ik, ip) = exp(-sum((kpt(:)*self%gpr_length(:))**2))
4455 : enddo
4456 : enddo
4457 0 : end subroutine gau_length_
4458 :
4459 0 : subroutine random_()
4460 0 : real(dp), allocatable :: re_rand(:,:), im_rand(:,:)
4461 0 : ABI_MALLOC(re_rand, (gqk%nb_k, gqk%my_nk))
4462 0 : ABI_MALLOC(im_rand, (gqk%nb_k, gqk%my_nk))
4463 0 : call random_number(re_rand)
4464 0 : call random_number(im_rand)
4465 0 : self%my_a(:,:,ip) = re_rand(:,:) + j_dpc*im_rand(:,:)
4466 0 : call xmpi_sum(self%my_a(:,:,ip), gqk%qpt_pert_comm%value, ierr)
4467 0 : ABI_FREE(re_rand)
4468 0 : ABI_FREE(im_rand)
4469 0 : end subroutine random_
4470 :
4471 : end subroutine polstate_seed_a
4472 : !!***
4473 :
4474 : !----------------------------------------------------------------------
4475 :
4476 : !!****f* m_varpeq/polstate_get_sqnorm
4477 : !! NAME
4478 : !! polstate_get_sqnorm
4479 : !!
4480 : !! FUNCTION
4481 : !! Helper function that calculates squared L^2-norm of an MPI-distributed
4482 : !! array at current state.
4483 : !!
4484 : !! INPUTS
4485 : !! mode=Select which array to gather. Possible options:
4486 : !! "a" ---> A_nk coefficients for this state;
4487 : !! "b" ---> B_q\nu coefficients for this state;
4488 : !! "grad" ---> current gradient;
4489 : !! "pcjgrad" ---> current preconditioned conjugate gradient direction.
4490 : !! ip=Index of a polaronic state.
4491 : !!
4492 : !! OUTPUT
4493 : !! Squared L^2-norm of the specified array.
4494 : !!
4495 : !! SOURCE
4496 :
4497 4580 : real(dp) function polstate_get_sqnorm(self, mode, ip) result(sqnorm)
4498 :
4499 : !Arguments ------------------------------------
4500 : class(polstate_t), target, intent(inout) :: self
4501 : character(len=*), intent(in) :: mode
4502 : integer, intent(in) :: ip
4503 :
4504 : !Local variables-------------------------------
4505 : class(gqk_t), pointer :: gqk
4506 : !----------------------------------------------------------------------
4507 :
4508 4580 : gqk => self%gqk
4509 4599 : select case(mode)
4510 : case ("a")
4511 19 : sqnorm = get_sqnorm_(self%my_a(:,:,ip), gqk%kpt_comm%value)
4512 : case ("b")
4513 0 : sqnorm = get_sqnorm_(self%my_b(:,:,ip), gqk%qpt_pert_comm%value)
4514 : case ("grad")
4515 2285 : sqnorm = get_sqnorm_(self%my_grad, gqk%kpt_comm%value)
4516 : case ("pcjgrad")
4517 2141 : sqnorm = get_sqnorm_(self%my_pcjgrad, gqk%kpt_comm%value)
4518 : case ("phgrad")
4519 135 : sqnorm = get_sqnorm_(self%my_phgrad(:,:,ip), gqk%qpt_pert_comm%value)
4520 : case ("eff_phforce")
4521 0 : sqnorm = get_sqnorm_(self%my_eff_phforce(:,:,ip), gqk%qpt_pert_comm%value)
4522 : case default
4523 4580 : ABI_ERROR(sjoin("polstate_get_sqnorm, unsupported mode: ", mode))
4524 : end select
4525 :
4526 : !----------------------------------------------------------------------
4527 :
4528 : contains
4529 4580 : real(dp) function get_sqnorm_(my_arr, comm) result(sqnorm)
4530 : integer, intent(in) :: comm
4531 : complex(dp), intent(in) :: my_arr(:, :)
4532 : integer :: ierr
4533 2154257 : sqnorm = sum(abs(my_arr(:,:))**2)
4534 4580 : call xmpi_sum(sqnorm, comm, ierr)
4535 4580 : end function get_sqnorm_
4536 :
4537 : end function polstate_get_sqnorm
4538 : !!***
4539 :
4540 : !----------------------------------------------------------------------
4541 :
4542 : !!****f* m_varpeq/polstate_gather
4543 : !! NAME
4544 : !! polstate_gather
4545 : !!
4546 : !! FUNCTION
4547 : !! Helper function that gathers a MPI-distributed array from the datatype
4548 : !! into a global one.
4549 : !!
4550 : !! INPUTS
4551 : !! mode=Select which array to gather. Possible options:
4552 : !! "a" ---> my_a(:,:,ip) -- array of A_nk coefficients for this state;
4553 : !! "pcjgrad" ---> my_pcjgrad(:,:) -- gradient.
4554 : !! ip=Index of a polaronic state.
4555 : !!
4556 : !! OUTPUT
4557 : !!
4558 : !! SOURCE
4559 :
4560 4426 : subroutine polstate_gather(self, mode, ip)
4561 :
4562 : !Arguments ------------------------------------
4563 : class(polstate_t), target, intent(inout) :: self
4564 : character(len=*),intent(in) :: mode
4565 : integer, intent(in) :: ip
4566 :
4567 : !Local variables-------------------------------
4568 : class(gqk_t), pointer :: gqk
4569 : !----------------------------------------------------------------------
4570 :
4571 4426 : gqk => self%gqk
4572 6711 : select case(mode)
4573 : case ("a")
4574 : call gather_(self%my_a(:,:,ip), gqk%my_nk, gqk%my_kstart, self%a_glob(:,:,ip), &
4575 2285 : gqk%kpt_comm%value)
4576 : case ("pcjgrad")
4577 : call gather_(self%my_pcjgrad, gqk%my_nk, gqk%my_kstart, self%pcjgrad_glob, &
4578 2141 : gqk%kpt_comm%value)
4579 : case default
4580 4426 : ABI_ERROR(sjoin("polstate_gather, unsupported mode: ", mode))
4581 : end select
4582 :
4583 : !----------------------------------------------------------------------
4584 :
4585 : contains
4586 4426 : subroutine gather_(my_arr, my_nk, my_kstart, glob_arr, comm)
4587 :
4588 : integer, intent(in) :: comm, my_nk, my_kstart
4589 : complex(dp), intent(in) :: my_arr(:, :)
4590 : complex(dp), intent(out) :: glob_arr(:, :)
4591 :
4592 : integer :: ierr, my_ik, ik_glob
4593 : !----------------------------------------------------------------------
4594 :
4595 2029248 : glob_arr(:, :) = zero
4596 513168 : do my_ik=1,my_nk
4597 508742 : ik_glob = my_ik + my_kstart - 1
4598 2029248 : glob_arr(:, ik_glob) = my_arr(:, my_ik)
4599 : enddo
4600 4426 : call xmpi_sum(glob_arr, comm, ierr)
4601 :
4602 4426 : end subroutine gather_
4603 :
4604 : end subroutine polstate_gather
4605 : !!***
4606 :
4607 : !----------------------------------------------------------------------
4608 :
4609 : !!****f* m_varpeq/polstate_filter
4610 : !! NAME
4611 : !! polstate_filter
4612 : !!
4613 : !! FUNCTION
4614 : !! Helper function that filters a MPI-distributed array from the datatype
4615 : !!
4616 : !! INPUTS
4617 : !! mode=Select which array to filter. Possible options:
4618 : !! "a" ---> my_a(:,:,ip) -- array of A_nk coefficients for this state;
4619 : !! "pcjgrad" ---> my_pcjgrad(:,:) -- gradient.
4620 : !! ip=Index of a polaronic state.
4621 : !!
4622 : !! OUTPUT
4623 : !!
4624 : !! SOURCE
4625 :
4626 0 : subroutine polstate_filter(self, mode, ip)
4627 :
4628 : !Arguments ------------------------------------
4629 : class(polstate_t), target, intent(inout) :: self
4630 : character(len=*),intent(in) :: mode
4631 : integer, intent(in) :: ip
4632 :
4633 : !Local variables-------------------------------
4634 : class(gqk_t), pointer :: gqk
4635 : !----------------------------------------------------------------------
4636 :
4637 0 : gqk => self%gqk
4638 0 : select case(mode)
4639 : case ("a")
4640 0 : call filter_(self%my_a(:,:,ip))
4641 : case ("grad")
4642 0 : call filter_(self%my_grad(:,:))
4643 : case ("pcjgrad")
4644 0 : call filter_(self%my_pcjgrad(:,:))
4645 : case default
4646 0 : ABI_ERROR(sjoin("polstate_filter, unsupported mode: ", mode))
4647 : end select
4648 :
4649 : !----------------------------------------------------------------------
4650 :
4651 : contains
4652 0 : subroutine filter_(my_arr)
4653 :
4654 : complex(dp), intent(inout) :: my_arr(:, :)
4655 :
4656 : integer :: my_ik, ik_ibz, ib
4657 : !----------------------------------------------------------------------
4658 :
4659 0 : do my_ik=1,gqk%my_nk
4660 0 : ik_ibz = gqk%my_k2ibz(1, my_ik)
4661 0 : do ib=1,gqk%nb_k
4662 0 : if (self%eig(ib, ik_ibz) > self%efilter .or. self%eig(ib, ik_ibz) < -tol6) then
4663 0 : my_arr(ib, my_ik) = zero
4664 : endif
4665 : enddo
4666 : enddo
4667 :
4668 0 : end subroutine filter_
4669 :
4670 : end subroutine polstate_filter
4671 : !!***
4672 :
4673 : !----------------------------------------------------------------------
4674 :
4675 : !!****f* m_varpeq/polstate_get_krank_glob
4676 : !! NAME
4677 : !! polstate_get_krank_glob
4678 : !!
4679 : !! FUNCTION
4680 : !! Helper function that computes krank objects for the full BZ q/k-meshes.
4681 : !!
4682 : !! INPUTS
4683 : !! mode=Select reciprocal space for which krank is computed. Possible options:
4684 : !! "k" ---> k-space;
4685 : !! "q" ---> q-space.
4686 : !! kptrlatt(3,3)=Lattice specifiyng the grid.
4687 : !!
4688 : !! OUTPUT
4689 : !! krank<krank_t>=Object proiding mapping between reciprocal space points and
4690 : !! their indices.
4691 : !!
4692 : !! SOURCE
4693 :
4694 18 : type(krank_t) function polstate_get_krank_glob(self, mode, kptrlatt) result(krank_kpts)
4695 :
4696 : !Arguments ------------------------------------
4697 : !scalars
4698 : class(polstate_t), target, intent(in) :: self
4699 : character(len=*), intent(in) :: mode
4700 : !arrays
4701 : integer, intent(in) :: kptrlatt(3, 3)
4702 :
4703 : !Local variables-------------------------------
4704 : !scalars
4705 : class(gqk_t), pointer :: gqk
4706 : !----------------------------------------------------------------------
4707 :
4708 18 : gqk => self%gqk
4709 :
4710 45 : select case(mode)
4711 : case ("k")
4712 : krank_kpts = get_krank_glob_(self%my_kpts, gqk%my_nk, gqk%my_kstart, &
4713 9 : gqk%glob_nk, gqk%kpt_comm%value)
4714 : case ("q")
4715 : krank_kpts = get_krank_glob_(self%my_qpts, gqk%my_nq, gqk%my_qstart, &
4716 9 : gqk%glob_nq, gqk%qpt_comm%value)
4717 : case default
4718 18 : ABI_ERROR(sjoin("polstate_get_krank_glob, unsupported mode: ", mode))
4719 : end select
4720 :
4721 : !----------------------------------------------------------------------
4722 :
4723 : contains
4724 18 : type(krank_t) function get_krank_glob_(my_kpts, my_nk, my_kstart, glob_nk, &
4725 : comm) result(krank_kpts)
4726 :
4727 : integer, intent(in) :: comm, my_nk, my_kstart, glob_nk
4728 : real(dp), intent(in) :: my_kpts(3, my_nk)
4729 :
4730 18 : type(krank_t) :: krank_tmp
4731 : integer :: ierr, my_ik, ik_glob
4732 36 : real(dp) :: kpts(3, glob_nk)
4733 :
4734 : !----------------------------------------------------------------------
4735 :
4736 7434 : kpts(:, :) = zero
4737 1872 : do my_ik=1,my_nk
4738 1854 : ik_glob = my_ik + my_kstart - 1
4739 7434 : kpts(:, ik_glob) = my_kpts(:, my_ik)
4740 : enddo
4741 18 : call xmpi_sum(kpts, comm, ierr)
4742 :
4743 18 : call krank_tmp%from_kptrlatt(glob_nk, kpts, kptrlatt, compute_invrank=.true.)
4744 18 : krank_kpts = krank_tmp%copy()
4745 18 : call krank_tmp%free()
4746 :
4747 36 : end function get_krank_glob_
4748 :
4749 : end function polstate_get_krank_glob
4750 : !!***
4751 :
4752 : !!****f* m_varpeq/varpeq_plot
4753 : !! NAME
4754 : !! varpeq_plot
4755 : !!
4756 : !! FUNCTION
4757 : !! Compute polaron wavefunctions and atomic displacements in the supercell and write results to XSF files
4758 : !!
4759 : !! INPUTS
4760 : !! wfk0_path=String with the path to the GS unperturbed WFK file.
4761 : !! ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes.
4762 : !! dtset<dataset_type>=All input variables for this dataset.
4763 : !! ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...)
4764 : !! comm=MPI communicator.
4765 : !!
4766 : !! OUTPUT
4767 : !!
4768 : !! SOURCE
4769 :
4770 2 : subroutine varpeq_plot(wfk0_path, ngfft, dtset, dtfil, cryst, ebands, pawtab, psps, comm)
4771 :
4772 : !Arguments ------------------------------------
4773 : character(len=*),intent(in) :: wfk0_path
4774 : integer,intent(in) :: ngfft(18) !,ngfftf(18)
4775 : type(dataset_type), intent(in) :: dtset
4776 : type(datafiles_type), intent(in) :: dtfil
4777 : type(crystal_t),intent(in) :: cryst
4778 : type(ebands_t),intent(in) :: ebands
4779 : type(pseudopotential_type),intent(in) :: psps
4780 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
4781 : integer,intent(in) :: comm
4782 :
4783 : !Local variables-------------------------------
4784 : !scalars
4785 : integer,parameter :: master = 0, ndat1 = 1
4786 : integer :: my_rank, ib, irsp, ir, uc_idx, mpw, band, bstart, ir1, ir2, ir3, wp1, wp2, wp3, ds_nfft, idir
4787 : integer :: natom, natom3, nsppol, nspinor, nspden, nkibz, mband, spin, ik, sc_nfft, cnt, nproc, num_writes, ii, uc_nfft
4788 : integer :: ik_ibz, isym_k, trev_k, npw_k, istwf_k, npw_kq_ibz, istwf_k_ibz, nkpg_k, ierr, nk, spinor, spad, nkbz, ncid
4789 : integer :: nqibz, iq, iq_ibz, isym_q, trev_q, qptopt, uc_iat, sc_iat, nu, nqbz, ip, ds_iscale
4790 : logical :: isirr_k, isirr_q, have_scell_q, use_displaced_scell
4791 : real(dp) :: cpu_all, wall_all, gflops_all, spread, kdotr
4792 : real(dp) :: psign
4793 : character(len=500) :: msg
4794 : character(len=fnlen) :: path
4795 104 : type(varpeq_t) :: vpq
4796 2 : type(wfd_t) :: wfd
4797 2 : type(supercell_type), target :: scell_q, scell_k
4798 2 : type(krank_t) :: krank_ibz, qrank_ibz
4799 : complex(dp) :: a_nk, bstar_qnu, cphase, cphase_tr, c3tmp(3)
4800 : complex(gwp) :: c123, c23, c3
4801 : !arrays
4802 : integer :: sc_ngfft(18), ds_ngfft(18), mapl_k(6), kptrlatt_(3,3), qptrlatt_(3,3)
4803 : integer :: units(2), work_ngfft(18), gmax(3), g0_k(3), mapl_qq(6), g0_q(3), ngqpt(3)
4804 4 : integer,allocatable :: nband(:,:), wfd_istwfk(:), kg_k(:,:), gbound_k(:,:)
4805 : real(dp),parameter :: origin0(3) = zero
4806 : real(dp) :: kk(3), kk_ibz(3), kk_sc(3), qq(3), qq_ibz(3), center_cart(3)
4807 4 : real(dp),allocatable :: kpg_k(:,:), ug_k(:,:), work(:,:,:,:), pol_rhor(:), qibz(:,:)
4808 2 : real(dp),allocatable :: pheigvec_qibz(:,:,:,:)
4809 2 : real(dp),allocatable :: displ_cart_qbz(:,:,:,:), pheigvec_qbz(:,:,:,:) !displ_red_qbz(:,:,:,:), displ_cart_qibz(:,:,:,:),
4810 2 : real(dp),allocatable :: phfreqs_ibz(:,:), pheigvec_cart_ibz(:,:,:,:,:) !, pheigvec_cart_qbz(:,:,:,:)
4811 2 : real(dp),allocatable :: sc_displ_cart_re(:,:,:,:), sc_displ_cart_im(:,:,:,:)
4812 2 : real(dp), contiguous, pointer :: xcart_ptr(:,:)
4813 2 : logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:)
4814 2 : complex(gwp),allocatable :: ur_k(:,:), ds_ur_k(:,:), pol_wfr(:,:,:), sc_ceikr_1d(:,:)
4815 : !----------------------------------------------------------------------
4816 :
4817 6 : units = [std_out, ab_out]
4818 2 : my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm)
4819 :
4820 : ! Read A_nk and B_qnu and other useful tables from file
4821 2 : call vpq%ncread(dtfil%filvpqin, comm, keep_open=.false.)
4822 : !call wrtout(std_out, " Reading done")
4823 :
4824 2 : psign = -1
4825 2 : if (vpq%pkind == "hole") psign = 1
4826 :
4827 : ! Copy important dimensions
4828 2 : natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden
4829 2 : nkibz = ebands%nkpt; mband = ebands%mband
4830 :
4831 2 : if (dtfil%filgstorein == ABI_NOFILE) then
4832 0 : call wrtout(units, "gstore_filepath is not specified in input. Cannot compute polaron-induced displacements!")
4833 0 : have_scell_q = .false.
4834 :
4835 : else
4836 : ! Start by reading ph displacements and frequencies in the IBZ from the gstore file.
4837 : ! First compute displaced supercell then polaron wf so that we can use both when writing the XSF file.
4838 2 : call wrtout(units, sjoin(" Computing polaron-induced displacements. Reading phonons from: ", dtfil%filgstorein))
4839 2 : call cwtime(cpu_all, wall_all, gflops_all, "start")
4840 2 : have_scell_q = .true.
4841 :
4842 2 : NCF_CHECK(nctk_open_read(ncid, dtfil%filgstorein, comm))
4843 2 : NCF_CHECK(nctk_get_dim(ncid, "gstore_nqibz", nqibz))
4844 : !NCF_CHECK(nctk_get_dim(ncid, "gstore_nqbz", nqbz))
4845 :
4846 : ! TODO: Wrap phstore API?
4847 : ! Encapsulate this part as we're gonna re-use it to deal with hopping
4848 : !call gstore_read_ph_qibz(dtfil%filgstorein, ph, comm)
4849 : !call ph%free()
4850 :
4851 6 : ABI_MALLOC(qibz, (3, nqibz))
4852 8 : ABI_MALLOC(phfreqs_ibz, (natom3, nqibz))
4853 12 : ABI_MALLOC(pheigvec_cart_ibz, (2, 3, cryst%natom, cryst%natom * 3, nqibz))
4854 2 : if (nproc > 1) then
4855 0 : NCF_CHECK(nctk_set_collective(ncid, vid("gstore_qibz")))
4856 0 : NCF_CHECK(nctk_set_collective(ncid, vid("phfreqs_ibz")))
4857 0 : NCF_CHECK(nctk_set_collective(ncid, vid("pheigvec_cart_ibz")))
4858 : end if
4859 2 : NCF_CHECK(nf90_get_var(ncid, vid("gstore_qibz"), qibz))
4860 2 : NCF_CHECK(nf90_get_var(ncid, vid("phfreqs_ibz"), phfreqs_ibz))
4861 2 : NCF_CHECK(nf90_get_var(ncid, vid("pheigvec_cart_ibz"), pheigvec_cart_ibz))
4862 2 : NCF_CHECK(nf90_get_var(ncid, vid("gstore_ngqpt"), ngqpt))
4863 2 : NCF_CHECK(nf90_get_var(ncid, vid("gstore_qptopt"), qptopt))
4864 2 : NCF_CHECK(nf90_close(ncid))
4865 :
4866 8 : ABI_MALLOC(pheigvec_qibz, (2, 3, cryst%natom, natom3))
4867 8 : ABI_MALLOC(displ_cart_qbz, (2, 3, cryst%natom, cryst%natom * 3))
4868 6 : ABI_MALLOC(pheigvec_qbz, (2, 3, cryst%natom, 3*cryst%natom))
4869 :
4870 8 : nqbz = product(ngqpt)
4871 2 : call kptrlatt_from_ngkpt(ngqpt, qptrlatt_)
4872 2 : call qrank_ibz%from_kptrlatt(nqibz, qibz, qptrlatt_, compute_invrank=.false.)
4873 :
4874 2 : call scell_q%init(cryst%natom, qptrlatt_, cryst%rprimd, cryst%typat, cryst%xcart, cryst%znucl, xyz_order="xyz")
4875 :
4876 18240 : ABI_CALLOC(sc_displ_cart_re, (3, scell_q%natom, vpq%nstates, nsppol))
4877 18238 : ABI_CALLOC(sc_displ_cart_im, (3, scell_q%natom, vpq%nstates, nsppol))
4878 :
4879 2 : cnt = 0
4880 4 : do spin=1,nsppol
4881 628 : do iq=1,vpq%nq_spin(spin)
4882 624 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism inside comm.
4883 2496 : qq = vpq%qpts_spin(:, iq, spin)
4884 : ! Note symrec here
4885 624 : if (kpts_map("symrec", qptopt, cryst, qrank_ibz, 1, qq, mapl_qq) /= 0) then
4886 0 : ABI_ERROR("Cannot map qBZ to IBZ!")
4887 : end if
4888 624 : iq_ibz = mapl_qq(1); isym_q = mapl_qq(2)
4889 624 : trev_q = mapl_qq(6); g0_q = mapl_qq(3:5)
4890 : ! Don't test if umklapp == 0 because we use the periodic gauge:
4891 : !
4892 : ! phfreq(q+G) = phfreq(q) and eigvec(q) = eigvec(q+G)
4893 : !
4894 624 : isirr_q = (isym_q == 1 .and. trev_q == 0)
4895 2496 : qq_ibz = qibz(:, iq_ibz)
4896 : !if (all(abs(qq_ibz) < tol6)) cycle
4897 79872 : pheigvec_qibz = pheigvec_cart_ibz(:,:,:,:,iq_ibz)
4898 :
4899 624 : if (isirr_q) then
4900 : ! Compute phonon displacements in Cartesian coordinates
4901 34 : call phdispl_from_eigvec(cryst%natom, cryst%ntypat, cryst%typat, cryst%amu, pheigvec_qibz, displ_cart_qbz)
4902 :
4903 : else
4904 : ! Rotate phonon eigenvectors from q_ibz to q_bz.
4905 : ! This part is needed to enforce the gauge in the ph eigenvectors, including e(-q) = e(q)^*
4906 590 : call pheigvec_rotate(cryst, qq_ibz, isym_q, trev_q, pheigvec_qibz, pheigvec_qbz, displ_cart_qbz)
4907 : end if
4908 :
4909 : ! Phase due to the primitive translation (default: 0)
4910 2496 : cphase_tr = exp(-j_dpc * two_pi * dot_product(qq, dtset%vpq_trvec))
4911 :
4912 2625302 : do sc_iat=1, scell_q%natom
4913 2624676 : uc_iat = scell_q%atom_indexing(sc_iat)
4914 : ! Compute phase e^{iq.R}
4915 10498704 : cphase = exp(+j_dpc * two_pi * dot_product(qq, scell_q%uc_indexing(:, sc_iat)))
4916 5252892 : do ip=1,vpq%nstates
4917 : ! Summing over ph modes.
4918 21017820 : do nu=1,natom3
4919 : ! skip acoustic modes at Gamma
4920 15765552 : if (abs(phfreqs_ibz(nu, iq_ibz)) < tol12) cycle
4921 :
4922 15751884 : bstar_qnu = conjg(vpq%b_spin(nu, iq, ip, spin)) * cphase_tr / sqrt(phfreqs_ibz(nu, iq_ibz))
4923 63007536 : c3tmp = (displ_cart_qbz(1,:,uc_iat,nu) + j_dpc * displ_cart_qbz(2,:,uc_iat,nu)) * bstar_qnu * cphase
4924 63007536 : sc_displ_cart_re(:,sc_iat,ip,spin) = sc_displ_cart_re(:,sc_iat,ip,spin) + real(c3tmp)
4925 65635128 : sc_displ_cart_im(:,sc_iat,ip,spin) = sc_displ_cart_im(:,sc_iat,ip,spin) + aimag(c3tmp) ! This to check if the imag part is zero.
4926 : end do
4927 : end do ! ip
4928 : end do ! sc_iat
4929 :
4930 : end do ! iq
4931 : end do ! spin
4932 :
4933 2 : ABI_FREE(qibz)
4934 2 : ABI_FREE(phfreqs_ibz)
4935 2 : ABI_FREE(pheigvec_cart_ibz)
4936 2 : ABI_FREE(displ_cart_qbz)
4937 2 : ABI_FREE(pheigvec_qbz)
4938 2 : ABI_FREE(pheigvec_qibz)
4939 2 : call qrank_ibz%free()
4940 :
4941 2 : call xmpi_sum_master(sc_displ_cart_re, master, comm, ierr)
4942 2 : call xmpi_sum_master(sc_displ_cart_im, master, comm, ierr)
4943 18232 : sc_displ_cart_re = psign * sqrt2 * sc_displ_cart_re / nqbz
4944 18232 : sc_displ_cart_im = psign * sqrt2 * sc_displ_cart_im / nqbz
4945 :
4946 : ! Write polaron-induced displacements in XSF format.
4947 2 : if (my_rank == master) then
4948 : ! Handle spin-polarized case by writing two XSF files.
4949 4 : do spin=1,nsppol
4950 : ! Handle multiple polaronic states for each spin.
4951 8 : do ip=1,vpq%nstates
4952 : write(msg, "(2(a,i0),a,es16.6)") &
4953 18228 : " For spin: ", spin, ": pstate: ", ip, ": maxval(abs(sc_displ_cart_re)): ", maxval(abs(sc_displ_cart_re(:,:,ip,spin)))
4954 4 : call wrtout(units, msg)
4955 : write(msg, "(2(a,i0),a,es16.6)") &
4956 18228 : " For spin: ", spin, ": pstate: ", ip, ": maxval(abs(sc_displ_cart_im)): ", maxval(abs(sc_displ_cart_im(:,:,ip,spin)))
4957 4 : call wrtout(units, msg)
4958 :
4959 : ! Here we displace the atoms in the supercell for this spin (only master has the correct values)
4960 18232 : scell_q%xcart = scell_q%xcart_ref + sc_displ_cart_re(:,:,ip,spin)
4961 :
4962 4 : path = strcat(dtfil%filnam_ds(4), "_pstate_", itoa(ip), "_POLARON_DISPL_VECTORS.xsf")
4963 4 : if (nsppol == 2) path = strcat(dtfil%filnam_ds(4), strcat("_spin_", itoa(spin)), "_pstate_", itoa(ip), "_POLARON_DISPL_VECTORS.xsf")
4964 4 : call wrtout(units, sjoin(" Writing displacement vectors to:", path))
4965 6 : call scell_q%write_xsf(path)
4966 : end do
4967 : end do
4968 : end if
4969 :
4970 2 : ABI_FREE(sc_displ_cart_im)
4971 : call cwtime_report(" Computation of polaron-induced displacements completed:", cpu_all, wall_all, gflops_all, &
4972 4 : pre_str=ch10, end_str=ch10)
4973 : end if
4974 :
4975 2 : call wrtout(std_out, " varpeq_plot: computing polaron wavefunction in real space.", pre_newlines=1)
4976 :
4977 2 : call krank_ibz%from_kptrlatt(ebands%nkpt, ebands%kptns, ebands%kptrlatt, compute_invrank=.false.)
4978 :
4979 : ! Initialize the wave function descriptor.
4980 8 : ABI_MALLOC(nband, (nkibz, nsppol))
4981 10 : ABI_MALLOC(bks_mask, (mband, nkibz, nsppol))
4982 8 : ABI_MALLOC(keep_ur, (mband, nkibz, nsppol))
4983 1328 : nband = mband; bks_mask = .false.; keep_ur = .false.
4984 :
4985 : ! Here we use brange_spin to select the number of bands that should be read and stored in memory.
4986 : ! For the time being, spin and k-points are not MPI-distributed inside comm.
4987 4 : do spin=1,nsppol
4988 2 : bstart = vpq%brange_spin(1, spin)
4989 120 : do ik=1, vpq%nk_spin(spin)
4990 464 : kk = vpq%kpts_spin(:, ik, spin)
4991 : ! Note symrel option
4992 116 : if (kpts_map("symrel", ebands%kptopt, cryst, krank_ibz, 1, kk, mapl_k) /= 0) then
4993 0 : write(msg, '(4a)' )"k-mesh is not closed!",ch10, "k-point could not be generated from a symmetrical one.",trim(ltoa(kk))
4994 0 : ABI_ERROR(msg)
4995 : end if
4996 116 : ik_ibz = mapl_k(1)
4997 288 : do ib=1,vpq%nb_spin(spin)
4998 286 : band = bstart + ib - 1; bks_mask(band, ik_ibz, spin) = .true.
4999 : end do
5000 : end do
5001 : end do
5002 :
5003 : ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts)
5004 : ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1.
5005 6 : ABI_MALLOC(wfd_istwfk, (nkibz))
5006 90 : wfd_istwfk = 1
5007 :
5008 : call wfd%init(cryst, pawtab, psps, keep_ur, mband, nband, nkibz, nsppol, bks_mask,&
5009 : dtset%nspden, dtset%nspinor, dtset%ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft, &
5010 2 : dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm)
5011 :
5012 4 : call wfd%print([std_out], header="Wavefunctions for varpeq_plot")
5013 :
5014 2 : if (dtset%boxcutmin >= two) then
5015 2 : call wrtout(std_out, " To decrease the size of the FFT mesh and the size of the XSF file, reduce boxcutmin from 2 to e.g. 1.1")
5016 : end if
5017 :
5018 2 : ABI_FREE(nband)
5019 2 : ABI_FREE(keep_ur)
5020 2 : ABI_FREE(wfd_istwfk)
5021 2 : ABI_FREE(bks_mask)
5022 :
5023 : ! Read wavefunctions.
5024 2 : call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path))
5025 :
5026 2 : call cwtime(cpu_all, wall_all, gflops_all, "start")
5027 2 : call wrtout(std_out, " Computing mpw and gmax needed to allocate workspace array.")
5028 4 : do spin=1,nsppol
5029 2 : nk = vpq%nk_spin(spin)
5030 : call ephtk_get_mpw_gmax(nk, vpq%kpts_spin(:, 1:nk, spin), dtset%ecut, cryst%gmet, mpw, gmax, comm, &
5031 4 : init_with_zero=spin==1)
5032 : end do
5033 :
5034 : ! Init work_ngfft
5035 8 : gmax = gmax + 4 ! FIXME: this is to account for umklapp, should also consider Gamma-only and istwfk
5036 8 : gmax = 2*gmax + 1
5037 2 : call ngfft_seq(work_ngfft, gmax)
5038 : !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3)
5039 10 : ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6)))
5040 :
5041 2 : call wrtout(std_out, sjoin(" Building supercell from ngkpt:", ltoa(vpq%ngkpt)))
5042 : ! note xyz_order="xyz"
5043 2 : call kptrlatt_from_ngkpt(vpq%ngkpt, kptrlatt_)
5044 2 : call scell_k%init(cryst%natom, kptrlatt_, cryst%rprimd, cryst%typat, cryst%xcart, cryst%znucl, xyz_order="xyz")
5045 :
5046 : ! There are three meshes.
5047 : ! ngfft: mesh in the unit cell used to compute u(r) via FFT
5048 : ! ds_ngfft: coarse mesh used to downsample u(r). Note ngfft(1:3)/ds_ngfft(1:3) = ds_iscale
5049 : ! sc_ngfft: mesh in the real space supercell computed as ds_ngfft * vqp%ngkpt
5050 : ! the polaron wavefunction is defined on this mesh.
5051 :
5052 2 : ds_iscale = dtset%vpq_mesh_fact
5053 2 : ds_ngfft = ngfft
5054 8 : ds_ngfft(1:3) = ngfft(1:3) / ds_iscale
5055 14 : ABI_CHECK(all(ds_ngfft(1:6) > 0), "ds_iscale too large and ds_ngfft == 0")
5056 2 : ds_ngfft(4:6) = ds_ngfft(4:6) ! No augmentation
5057 8 : ds_nfft = product(ds_ngfft(1:3)) ! Total number of points in the supercell
5058 :
5059 8 : sc_ngfft(1:3) = vpq%ngkpt(1:3) * ds_ngfft(1:3)
5060 8 : sc_ngfft(4:6) = sc_ngfft(1:3) ! No augmentation
5061 8 : sc_nfft = product(sc_ngfft(1:3)) ! Total number of points in the supercell
5062 :
5063 2 : call wrtout(std_out, " Computing polaron wavefunction in the real-space supercell...")
5064 2 : call wrtout(std_out, sjoin(" Using vpq_mesh_fact:", itoa(ds_iscale)))
5065 : call wrtout(std_out, sjoin(" Memory required by pol_wfr:", &
5066 2 : ftoa(two*sc_nfft*nspinor*vpq%nstates*nsppol*storage_size(one)/eight*b2Mb), " (Mb) <<< MEM"))
5067 :
5068 8 : nkbz = product(vpq%ngkpt)
5069 2 : uc_nfft = wfd%nfft
5070 6 : ABI_MALLOC(kg_k, (3, mpw))
5071 8 : ABI_MALLOC(gbound_k, (2*wfd%mgfft+8, 2))
5072 6 : ABI_MALLOC(ur_k, (uc_nfft*nspinor, ndat1))
5073 2 : if (ds_iscale /= 1) then
5074 6 : ABI_MALLOC(ds_ur_k, (ds_nfft*nspinor, ndat1))
5075 : end if
5076 :
5077 747943 : ABI_CALLOC(pol_wfr, (sc_nfft*nspinor, vpq%nstates, nsppol)) ! Init output with zeros.
5078 12 : ABI_MALLOC(sc_ceikr_1d, (maxval(sc_ngfft(1:3)), 3))
5079 :
5080 2 : cnt = 0
5081 4 : do spin=1,nsppol
5082 2 : bstart = vpq%brange_spin(1, spin)
5083 2 : nk = vpq%nk_spin(spin)
5084 120 : do ik=1, nk
5085 116 : cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism inside comm.
5086 :
5087 464 : kk = vpq%kpts_spin(:, ik, spin)
5088 116 : if (kpts_map("symrel", ebands%kptopt, cryst, krank_ibz, 1, kk, mapl_k) /= 0) then
5089 0 : write(msg, '(4a)' )"k-mesh is not closed!",ch10, "k-point could not be generated from a symmetrical one.",trim(ltoa(kk))
5090 0 : ABI_ERROR(msg)
5091 : end if
5092 :
5093 464 : ik_ibz = mapl_k(1); isym_k = mapl_k(2); trev_k = mapl_k(6); g0_k = mapl_k(3:5)
5094 452 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
5095 464 : kk_ibz = ebands%kptns(:, ik_ibz)
5096 116 : istwf_k_ibz = wfd%istwfk(ik_ibz); npw_kq_ibz = wfd%npwarr(ik_ibz)
5097 :
5098 : ! Get npw_k, kg_k for this k-point.
5099 : call wfd%get_gvec_gbound(cryst%gmet, dtset%ecut, kk, ik_ibz, isirr_k, dtset%nloalg, & ! in
5100 116 : istwf_k, npw_k, kg_k, nkpg_k, kpg_k, gbound_k) ! out
5101 :
5102 348 : ABI_MALLOC(ug_k, (2, npw_k*nspinor))
5103 464 : kk_sc = kk * vpq%ngkpt
5104 :
5105 : ! Precompute 1d phases e^{ik_j R_j} on the supercell.
5106 464 : do idir=1,3
5107 24935 : do ii=0,sc_ngfft(idir) - 1
5108 24471 : kdotr = two_pi * (kk_sc(idir) * (ii / dble(sc_ngfft(idir))))
5109 24819 : sc_ceikr_1d(ii+1, idir) = dcmplx(cos(kdotr), sin(kdotr))
5110 : end do
5111 : end do
5112 :
5113 : ! Phase due to the primitive translation (default: 0)
5114 464 : cphase_tr = exp(-j_dpc * two_pi * dot_product(kk, dtset%vpq_trvec))
5115 :
5116 : ! Sum over bands
5117 286 : do ib=1,vpq%nb_spin(spin)
5118 170 : band = bstart + ib - 1
5119 : ! Get periodic part on the real space FFT mesh in ur_k
5120 : call wfd%rotate_cg(band, ndat1, spin, kk_ibz, npw_k, kg_k, istwf_k, &
5121 170 : cryst, mapl_k, gbound_k, work_ngfft, work, ug_k, urs_kbz=ur_k)
5122 : !print *, "int_omega dr |u(r)}^2:", sum(abs(ur_k(:,1)) ** 2) / uc_nfft
5123 :
5124 : ! Use linear interpolate to downsample from ngfft to ds_ngfft.
5125 170 : if (ds_iscale /= 1) then
5126 170 : call interpolate_ur(ngfft, nspinor*ndat1, ur_k, ds_ngfft, ds_ur_k)
5127 : !print *, "int_omega dr |ds_u(r)}^2:", sum(abs(ds_ur_k(:,1)) ** 2) / product(ds_ngfft(1:3))
5128 : end if
5129 :
5130 : ! Accumulate for each polaron state.
5131 618 : do ip=1,vpq%nstates
5132 332 : a_nk = vpq%a_spin(ib, ik, ip, spin) * cphase_tr
5133 834 : do spinor=1,nspinor
5134 332 : spad = (spinor - 1) * sc_nfft
5135 :
5136 664 : if (ds_iscale == 1) then
5137 : ! Loop over the points in the R supercell without downsampling.
5138 0 : ir = 0
5139 0 : do ir3=0,sc_ngfft(3)-1
5140 : ! The FFT index of the point wrapped into the unit cell.
5141 0 : wp3 = modulo(ir3, ngfft(3))
5142 0 : c3 = sc_ceikr_1d(ir3+1, 3)
5143 0 : do ir2=0,sc_ngfft(2)-1
5144 0 : wp2 = modulo(ir2, ngfft(2))
5145 0 : c23 = sc_ceikr_1d(ir2+1, 2) * c3
5146 0 : do ir1=0,sc_ngfft(1)-1
5147 0 : wp1 = modulo(ir1, ngfft(1))
5148 0 : c123 = sc_ceikr_1d(ir1+1, 1) * c23
5149 0 : uc_idx = 1 + wp1 + wp2*ngfft(1) + wp3*ngfft(1)*ngfft(2)
5150 0 : ir = ir + 1; irsp = ir + spad
5151 0 : pol_wfr(irsp, ip, spin) = pol_wfr(irsp, ip, spin) + a_nk * ur_k(uc_idx,1) * c123
5152 : end do
5153 : end do
5154 : end do
5155 : else
5156 : ! Loop over the points in the supercell with downsampling
5157 332 : ir = 0
5158 18209 : do ir3=0,sc_ngfft(3)-1
5159 : ! The FFT index of the point wrapped into the unit cell.
5160 17877 : wp3 = modulo(ir3, ds_ngfft(3))
5161 17877 : c3 = sc_ceikr_1d(ir3+1, 3)
5162 1051760 : do ir2=0,sc_ngfft(2)-1
5163 1033551 : wp2 = modulo(ir2, ds_ngfft(2))
5164 1033551 : c23 = sc_ceikr_1d(ir2+1, 2) * c3
5165 65429931 : do ir1=0,sc_ngfft(1)-1
5166 64378503 : wp1 = modulo(ir1, ds_ngfft(1))
5167 64378503 : c123 = sc_ceikr_1d(ir1+1, 1) * c23
5168 64378503 : uc_idx = 1 + wp1 + wp2*ds_ngfft(1) + wp3*ds_ngfft(1)*ds_ngfft(2)
5169 64378503 : ir = ir + 1; irsp = ir + spad
5170 : ! Note the use of the downsampled ds_ur_k here.
5171 65412054 : pol_wfr(irsp, ip, spin) = pol_wfr(irsp, ip, spin) + a_nk * ds_ur_k(uc_idx,1) * c123
5172 : end do
5173 : end do
5174 : end do
5175 : end if
5176 :
5177 : end do
5178 : end do ! ip
5179 : end do ! ib
5180 :
5181 116 : ABI_FREE(ug_k)
5182 118 : ABI_FREE(kpg_k)
5183 : end do ! ik
5184 : end do ! spin
5185 :
5186 2 : ABI_FREE(kg_k)
5187 2 : ABI_FREE(work)
5188 2 : ABI_FREE(ur_k)
5189 2 : ABI_SFREE(ds_ur_k)
5190 2 : ABI_FREE(gbound_k)
5191 2 : ABI_FREE(sc_ceikr_1d)
5192 2 : call krank_ibz%free()
5193 :
5194 : ! Collect pol_wfr on the master rank who's gonna write the polaron density in XSF format.
5195 2 : call xmpi_sum_master(pol_wfr, master, comm, ierr)
5196 :
5197 2 : if (my_rank == master) then
5198 747935 : pol_wfr = pol_wfr / (nkbz * sqrt(cryst%ucvol))
5199 : ! FIXME: Here we're gonna have another big allocation
5200 6 : ABI_MALLOC(pol_rhor, (sc_nfft))
5201 :
5202 : ! Here decide if we are gonna write the polaron wavefunctions with_diplaced atoms or not.
5203 26 : use_displaced_scell = all(kptrlatt_ == qptrlatt_) .and. have_scell_q
5204 2 : num_writes = 1; if (use_displaced_scell) num_writes = 2
5205 :
5206 2 : if (nspinor == 1) then
5207 : ! Handle spin-polarized case by writing two XSF files.
5208 4 : do spin=1,nsppol
5209 : ! Handle multiple polaronic states for each spin.
5210 8 : do ip=1,vpq%nstates
5211 : write(msg, "(2(a,i0),a,es16.6)")&
5212 506 : " For spin: ", spin, ": pstate: ", ip, ": 1/N_k \sum_nk |A_nk|^2 = ", sum(abs(vpq%a_spin(:,:,ip,spin))**2) / nkbz
5213 4 : call wrtout(units, msg)
5214 747935 : pol_rhor = abs(pol_wfr(:, ip, spin)) ** 2
5215 4 : write(msg, "(2(a,i0),a,es16.6)")" Polaron density for spin: ", spin, ": pstate: ", ip, &
5216 747947 : " integrates to: ", sum(pol_rhor) * cryst%ucvol / product(ds_ngfft(1:3))
5217 4 : call wrtout(units, msg)
5218 747935 : write(msg, "(a,es16.6)")" maxval(abs(aimag(pol_wfr))): ", maxval(abs(aimag(pol_wfr(:, ip, spin))))
5219 4 : call wrtout(std_out, msg)
5220 4 : call center_and_spread(cryst, vpq%ngkpt, sc_ngfft, pol_rhor, center_cart, spread, units)
5221 :
5222 14 : do ii=1,num_writes
5223 8 : if (ii == 1) then
5224 4 : xcart_ptr => scell_k%xcart
5225 4 : path = strcat(dtfil%filnam_ds(4), "_pstate_", itoa(ip), "_POLARON.xsf")
5226 4 : if (nsppol == 2) path = strcat(dtfil%filnam_ds(4), strcat("_spin_", itoa(spin)), "_pstate_", itoa(ip), "_POLARON.xsf")
5227 4 : call wrtout(units, strcat("- Writing the polaron wavefunction with undisplaced atoms to: ", path))
5228 : else
5229 4 : path = strcat(dtfil%filnam_ds(4), "_pstate_", itoa(ip), "_POLARON_DISPL.xsf")
5230 4 : if (nsppol == 2) path = strcat(dtfil%filnam_ds(4), strcat("_spin_", itoa(spin)), "_pstate_", itoa(ip), "_POLARON_DISPL.xsf")
5231 4 : call wrtout(units, strcat("- Writing the polaron wavefunction with displaced atoms to: ", path))
5232 :
5233 : ! Here we displace the atoms in the supercell for this spin (only master has the correct values)
5234 18232 : scell_q%xcart = scell_q%xcart_ref + sc_displ_cart_re(:,:,ip,spin)
5235 4 : xcart_ptr => scell_q%xcart
5236 : end if
5237 :
5238 : call write_xsf(path, sc_ngfft(1), sc_ngfft(2), sc_ngfft(3), pol_rhor, scell_k%rprimd, origin0, &
5239 12 : scell_k%natom, scell_k%ntypat, scell_k%typat, xcart_ptr, scell_k%znucl, 0)
5240 : end do ! ii
5241 : end do ! ip
5242 : end do ! spin
5243 :
5244 : else
5245 : ! Spinor wavefunctions.
5246 0 : do ip=1,vpq%nstates
5247 : write(msg, "(2(a,i0),a,es16.6)")&
5248 0 : " For spin: ", spin, ": pstate: ", ip, ": 1/N_k \sum_nk |A_nk|^2 = ", sum(abs(vpq%a_spin(:,:,ip,spin))**2) / nkbz
5249 0 : call wrtout(units, msg)
5250 :
5251 0 : pol_rhor(:) = abs(pol_wfr(1:sc_nfft, ip, 1)) ** 2
5252 0 : pol_rhor(:) = abs(pol_wfr(sc_nfft+1:, ip, 1)) ** 2 + pol_rhor(:)
5253 0 : write(msg, "(2(a,i0),a,es16.6)")" Polaron density for spin: ", spin, ": pstate: ", ip, &
5254 0 : " integrates to: ", sum(pol_rhor) * cryst%ucvol / product(ds_ngfft(1:3))
5255 0 : call wrtout(units, msg)
5256 0 : write(msg, "(a,es16.6)")" maxval(abs(aimag(pol_wfr))): ", maxval(abs(aimag(pol_wfr(:, ip, 1))))
5257 0 : call wrtout(units, msg)
5258 0 : call center_and_spread(cryst, vpq%ngkpt, sc_ngfft, pol_rhor, center_cart, spread, units)
5259 :
5260 0 : spin = 1
5261 0 : do ii=1,num_writes
5262 0 : if (ii == 1) then
5263 0 : xcart_ptr => scell_k%xcart
5264 0 : path = strcat(dtfil%filnam_ds(4), "_pstate_", itoa(ip), "_POLARON.xsf")
5265 0 : call wrtout(units, strcat("- Writing the polaron wavefunction with undisplaced atoms to: ", path))
5266 : else
5267 0 : path = strcat(dtfil%filnam_ds(4), "_pstate_", itoa(ip), "_POLARON_DISPL.xsf")
5268 0 : call wrtout(units, strcat("- Writing the polaron wavefunction with displaced atoms to: ", path))
5269 :
5270 : ! Here we displace the atoms in the supercell for this spin (only master has the correct values)
5271 0 : scell_q%xcart = scell_q%xcart_ref + sc_displ_cart_re(:,:,ip,spin)
5272 0 : xcart_ptr => scell_q%xcart
5273 : end if
5274 :
5275 : call write_xsf(path, sc_ngfft(1), sc_ngfft(2), sc_ngfft(3), pol_rhor, scell_k%rprimd, origin0, &
5276 0 : scell_k%natom, scell_k%ntypat, scell_k%typat, xcart_ptr, scell_k%znucl, 0)
5277 : end do ! ii
5278 : end do ! ip
5279 : end if
5280 2 : ABI_FREE(pol_rhor)
5281 : end if ! master
5282 :
5283 2 : call cwtime_report(" Computation of polaron wavefunction completed", cpu_all, wall_all, gflops_all, pre_str=ch10, end_str=ch10)
5284 :
5285 2 : ABI_SFREE(sc_displ_cart_re)
5286 2 : ABI_FREE(pol_wfr)
5287 :
5288 8 : call wfd%free(); call scell_q%free(); call scell_k%free(); call vpq%free()
5289 :
5290 : contains
5291 10 : integer function vid(var_name)
5292 : character(len=*),intent(in) :: var_name
5293 10 : vid = nctk_idname(ncid, var_name)
5294 : end function vid
5295 :
5296 : end subroutine varpeq_plot
5297 : !!***
5298 :
5299 4 : subroutine center_and_spread(prim_cryst, ncells, sc_ngfft, rhor, center_cart, spread, units)
5300 :
5301 : !Arguments ------------------------------------
5302 : type(crystal_t),intent(in) :: prim_cryst
5303 : integer,intent(in) :: ncells(3), sc_ngfft(18)
5304 : real(dp),intent(in) :: rhor(sc_ngfft(1), sc_ngfft(2), sc_ngfft(3))
5305 : !real(dp),intent(in) :: rhor(sc_ngfft(4), sc_ngfft(5), sc_ngfft(6)) ! FIXME
5306 : real(dp),intent(out) :: center_cart(3), spread
5307 : integer,intent(in) :: units(:)
5308 :
5309 : !Local variables-------------------------------
5310 : integer :: i1, i2, i3, nfft, num_cells
5311 : real(dp) :: rr(3), rcart(3), rmr0(3), sc_rprimd(3,3), center_red(3), r2_mean, fwhm
5312 : character(len=500) :: msg
5313 : !----------------------------------------------------------------------
5314 :
5315 16 : nfft = product(sc_ngfft(1:3))
5316 16 : sc_rprimd(:,1) = prim_cryst%rprimd(:,1) * ncells(1)
5317 16 : sc_rprimd(:,2) = prim_cryst%rprimd(:,2) * ncells(2)
5318 16 : sc_rprimd(:,3) = prim_cryst%rprimd(:,3) * ncells(3)
5319 16 : num_cells = product(ncells)
5320 :
5321 : ! Compute center_cart = \int r rhor(r) dr
5322 4 : center_cart = zero !; r2_mean = zero
5323 217 : do i3=1,sc_ngfft(3)
5324 213 : rr(3) = (i3 - one) / sc_ngfft(3)
5325 12376 : do i2=1,sc_ngfft(2)
5326 12159 : rr(2) = (i2 - one) / sc_ngfft(2)
5327 760299 : do i1=1,sc_ngfft(1)
5328 747927 : rr(1) = (i1 - one) / sc_ngfft(1)
5329 : ! Go to cartesian coordinates.
5330 9723051 : rcart = matmul(sc_rprimd, rr)
5331 2991708 : center_cart = center_cart + rhor(i1, i2, i3) * rcart
5332 760086 : r2_mean = r2_mean + rhor(i1, i2, i3) * dot_product(rcart, rcart)
5333 : end do
5334 : end do
5335 : end do
5336 :
5337 : !center_cart = matmul(sc_rprimd, center_cart)
5338 16 : center_cart = center_cart * num_cells * prim_cryst%ucvol / (one*nfft)
5339 : !r2_mean = r2_mean * num_cells * prim_cryst%ucvol / (one*nfft)
5340 : !spread = sqrt(r2_mean - dot_product(center_cart, center_cart))
5341 : !write(msg, "(a,2(es16.6,a))")" Polaron spread: ", spread, " (Bohr)", spread * Bohr_Ang, " (Ang)"
5342 : !call wrtout(units, msg)
5343 :
5344 : ! Compute = \int (r - center_cart)^2 rhor dr = <r^2> - <r>^2
5345 4 : spread = zero
5346 217 : do i3=1,sc_ngfft(3)
5347 213 : rr(3) = (i3 - one) / sc_ngfft(3)
5348 12376 : do i2=1,sc_ngfft(2)
5349 12159 : rr(2) = (i2 - one) / sc_ngfft(2)
5350 760299 : do i1=1,sc_ngfft(1)
5351 747927 : rr(1) = (i1 - one) / sc_ngfft(1)
5352 : ! Go to cartesian coordinates.
5353 11966832 : rmr0 = matmul(sc_rprimd, rr) - center_cart
5354 3003867 : spread = spread + rhor(i1, i2, i3) * dot_product(rmr0, rmr0)
5355 : end do
5356 : end do
5357 : end do
5358 4 : spread = sqrt(spread * num_cells * prim_cryst%ucvol / (one*nfft))
5359 :
5360 4 : call xcart2xred(1, prim_cryst%rprimd, center_cart, center_red)
5361 4 : call wrtout(units, sjoin(" Polaron center in Cartesian coordinates: ", ltoa(center_cart), " (Bohr)"))
5362 4 : call wrtout(units, sjoin(" Fractional coordinates in terms of the primitive cell:", ltoa(center_red)))
5363 4 : write(msg, "(a,2(es16.6,a))")" Polaron spread: ", spread, " (Bohr)", spread * Bohr_Ang, " (Ang)"
5364 4 : call wrtout(units, msg)
5365 : ! full width at half-maximum
5366 4 : fwhm = (two * sqrt(two * log(two))) * spread
5367 4 : write(msg, "(a,2(es16.6,a))")" Full width at half-maximum (FWHM) ", fwhm, " (Bohr) ", fwhm * Bohr_Ang, " (Ang)"
5368 4 : call wrtout(units, msg)
5369 :
5370 4 : end subroutine center_and_spread
5371 :
5372 110592 : end module m_varpeq
5373 : !!***
|
