scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'EPW-5.4.1' into 'develop' 

EPW v.5.4.1

See merge request QEF/q-e!1644
This commit is contained in:
giannozz 2021-11-26 07:17:10 +00:00
parent 1d69093b87 0326659b6f
commit bcb15072ce
28 changed files with 2603 additions and 1378 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Doc/release-notes
View File

@ -4,6 +4,12 @@ New in development version:
* projwfc.x can be used to compute the PDOS in a local basis (I. Timrov)
* Charge self-consistent DFT+DMFT calculations with the TRIQS software package
via Wannier90 (S. Beck, CCQ Flatiron Institute)
* EPW:
(1) Addition of spectral decomposition capability in the transport module
(2) Support for the frmsf file format of FermiSurfer
(3) The cumulant module updated
For the full list of new features, bug fixes, and changes leading to backward incompatibility issues,
please visit the Releases page of the EPW documentation site [https://docs.epw-code.org/doc/Releases.html].
Fixed in development version:
* Possible out-of-bound error (gfortran only) could crash DFT+U

38
EPW/ZG/README
View File

@ -1,4 +1,4 @@
Marios Zacharias [1] & Feliciano Giustino [2,3], July 2021
Marios Zacharias [1] & Feliciano Giustino [2,3], November 2021
[1] Department of Mechanical and Materials Science Engineering, Cyprus University of Technology,
P.O. Box 50329, 3603 Limassol, Cyprus
@ -12,6 +12,11 @@ If you use ZG.x or bands_unfold.x please do cite the following work:
[1] Marios Zacharias and Feliciano Giustino, Phys. Rev. Research 2, 013357, (2020).
[2] Marios Zacharias and Feliciano Giustino, Phys. Rev. B 94, 075125, (2016).
If you use disca.x please do cite the following work:
[1] Marios Zacharias, Feliciano Giustino et al., Phys. Rev. B 104, 205109, (2021)
[2] Marios Zacharias, Feliciano Giustino et al., Phys. Rev. Lett. 127, 207401, (2021)
---------------------------------------------------------------------------------------------------
Acknowledgement: We thank Hyungjun Lee, Oden Institute for Computational Engineering and Sciences,
@ -27,9 +32,15 @@ pp_spctrlfn.x ---> for obtaining the electron spectral function after bands_u
epsilon_Gaus.x ---> for calculating optical properties as in epsilon.x but Gaussian broadening
disca.x ---> for calculating one-phonon and all-phonon inelastic scattering intensities
pp_disca.x ---> for applying broadening and setting a resolution of scattering patterns
src/local folder ---> fortran routines for post-processing. Compile them by "./compile_gfortran.sh"
src/local folder ---> for post-processing. Compile them by "./compile_gfortran.sh"
---------------------------------------------------------------------------------------------------
For full instructions on how to run the exercises in the tarball "tutorial.tar.gz"
please refer to the EPW documentation site, or send an email to Marios Zacharias:
zachariasmarios@gmail.com
---------------------------------------------------------------------------------------------------
Instructions for the construction of the Zacharias-Giustino "ZG" displacement following
Eq. (2) of Phys. Rev. Research 2, 013357, 2020. The approach for generating the ZG-displacement
@ -55,8 +66,8 @@ STEPS for generating the "ZG-displacement" for the calculation of temperature-de
for a relatively "dense" q-grid to obtain a converged phonon dispersion.
3. Run "/path_to_your_espresso/bin/q2r.x <q2r.in> q2r.out".
The input "q2r.in" has the standard format of Quantum Espresso for
calculating the interatomic force constant matrix elements to be used to construct the
The input "q2r.in" has the standard format of Quantum Espresso for calculating the
interatomic force constant matrix elements to be used to construct the
"ZG-configuration.dat" file.
4. Calculate the phonon dispersion with "matdyn.x" and make sure that your phonon dispersion is
@ -66,12 +77,10 @@ STEPS for generating the "ZG-displacement" for the calculation of temperature-de
included.
5. Now one needs to decide on the size of the supercell configuration to be used for calculating
temperature dependent properties. For help, please see the example folder by "tar -xvf example.tar.gz".
In file "example/silicon/ZG_structure/inputs/ZG_444.in" we show the example for
constructing a 4x4x4 ZG-configuration. One could potentially generate any supercell size
by simply changing "dim1","dim2","dim3", and the list of q-points (optional, see below).
"ZG.in" has the standard format as a "matdyn.in" file for Quantum Espresso.
Here we use the following input parameters:
temperature dependent properties. For help, please see exercise1 in tutorial tarball.
One could potentially generate any supercell size by simply changing "dim1","dim2",
"dim3", and the list of q-points (optional, see below). "ZG.in" has the standard
format as a "matdyn.in" file for Quantum Espresso. Here we use the following input parameters:
---------------------------------------------------------------------------------------
i) "ZG_conf" : Logical flag that enables the creation of the ZG-displacement.
(default .true.)
@ -102,7 +111,7 @@ STEPS for generating the "ZG-displacement" for the calculation of temperature-de
sizes are considered for which the error is minimized by the first set of signs,
(ii) only single phonon displacements are of interest (see below)
(default .true.)
"threshold" : Real number indicating the error at which the algorithm stops while it's
"error_thresh" : Real number indicating the error at which the algorithm stops while it's
looking for possible combinations of signs. Once this limit is reached,
the ZG-displacement is constructed. The threshold is usually chosen
to be less than 5% of the diagonal terms, i.e. those terms that contribute
@ -130,12 +139,11 @@ STEPS for generating the "ZG-displacement" for the calculation of temperature-de
Paste the output of "qlist_AB.txt" to "ZG.in" after namelist &input. Set the flag
q_external = .true. for the code to read the list.
ii) To generate the ZG-displacement run "/path_to_your_espresso/bin/ZG.x <ZG_444.in> ZG_444.out".
ii) To generate the ZG-displacement run "/path_to_your_espresso/bin/ZG.x <ZG.in> ZG.out".
This generates three output files: the "equil_pos.dat", "ZG-configuration.dat" and "ZG-velocities.dat".
The first file has the equilibrium coordinates of the nuclei and the second has the optimum set of nuclear coordinates
that define the ZG-displacement for a particular temperature and supercell size.
The third one has the ZG-velocities or momenta of the nuclei generated in the same spirit with the ZG displacement.
The outputs for a supercell size of 4x4x4 are in the folder "example/silicon/ZG_structure/output/".
iii) The calculation of the ZG-displacement should usually takes a few seconds to few minutes with one processor.
6. VERY IMPORTANT NOTE: It is perfectly reasonable to find different ZG-displacements / ZG-configurations,
@ -159,7 +167,7 @@ STEPS for generating the "ZG-displacement" for the calculation of temperature-de
----------------------------------------------------------------------------------------------------------
Example using the ZG-displacement
Example using the ZG-displacement and JDOS
Here we show a simple example on how to use the ZG-displacement for calculating the temperature dependence
of the JDOS and eventually extract the zero-point renormalization and temperature dependence of the band gap.
@ -206,7 +214,7 @@ Other capabilities of ZG.x:
1. ZG.x can also be used to generate displaced configurations along single phonon modes.
Please see "single_phonon_displ" flag.
Find an example in:
"example/silicon/Displacement_along_single_phonon_modes/inputs/".
"example/silicon/single_ph_displ/444/inputs/".
To run the example use:
"/path_to_your_espresso/bin/ZG.x <ZG_444.in> ZG_444.out".
The output files of interest are:

2
EPW/ZG/doc/README
View File

@ -1,2 +0,0 @@
For the exercises in the tutorial please download material from:
https://docs.epw-code.org/doc/Virtual2021.html

BIN
EPW/ZG/doc/input_variables.pdf
View File

Binary file not shown.

BIN
EPW/ZG/doc/tutorial.pdf
View File

Binary file not shown.

BIN
EPW/ZG/example.tar.gz
View File

Binary file not shown.

412
EPW/ZG/src/ZG.f90
View File

@ -80,8 +80,12 @@ PROGRAM ZG
! are given. See below. (default: .FALSE.).
! q_in_cryst_coord IF .TRUE. input q points are in crystalline
! coordinates (default: .FALSE.)
! fd (logical) if .t. the ifc come from the finite displacement calculation
! na_ifc (logical) add non analitic contributions to the interatomic force
! constants if finite displacement method is used (as in Wang et al.
! Phys. Rev. B 85, 224303 (2012)) [to be used in conjunction with fd.x]
! loto_2d set to .true. to activate two-dimensional treatment of LO-TO
! siplitting.
! siplitting for q= 0. (default: .false.)
!
! IF (q_in_band_form) THEN
! nq ! number of q points
@ -128,7 +132,7 @@ PROGRAM ZG
! (i) large supercell sizes are considered for which the error is minimized by
! the first set of signs, and (ii) only single phonon displacements are of interest (see below)
! (default .true.)
! "threshold" : Real number indicating the error at which the algorithm stops while it's
! "error_thresh" : Real number indicating the error at which the algorithm stops while it's
! looking for possible combinations of signs. Once this limit is reached
! the ZG-displacement is constructed. The threshold is usually chosen
! to be less than 5% of the diagonal terms, i.e. those terms that contribute
@ -187,7 +191,7 @@ PROGRAM ZG
INTEGER :: nr1, nr2, nr3, nsc, ibrav
CHARACTER(LEN=256) :: flfrc, filename
CHARACTER(LEN= 10) :: asr
LOGICAL :: has_zstar, q_in_cryst_coord
LOGICAL :: has_zstar, q_in_cryst_coord, loto_disable
COMPLEX(DP), ALLOCATABLE :: dyn(:, :, :, :), dyn_blk(:, :, :, :), frc_ifc(:, :, :, :)
COMPLEX(DP), ALLOCATABLE :: z(:, :)
REAL(DP), ALLOCATABLE:: tau(:, :), q(:, :), w2(:, :), wq(:)
@ -208,9 +212,9 @@ PROGRAM ZG
INTEGER :: nspin_mag, nqs, ios
!
LOGICAL :: xmlifc, lo_to_split, loto_2d
LOGICAL :: xmlifc, lo_to_split, loto_2d, na_ifc, fd, nosym
!
REAL(DP) :: qhat(3), qh, E
REAL(DP) :: qhat(3), qh, E, qq
REAL(DP) :: delta
REAL(DP), ALLOCATABLE :: xqaux(:, :)
INTEGER, ALLOCATABLE :: nqb(:)
@ -237,19 +241,22 @@ PROGRAM ZG
LOGICAL :: ZG_strf, compute_error, single_ph_displ
INTEGER :: dim1, dim2, dim3, niters, qpts_strf
REAL(DP) :: error_thresh, T
REAL(DP) :: atmsf_zg_a(ntypx,5), atmsf_zg_b(ntypx,5)
REAL(DP) :: atmsf_a(ntypx,5), atmsf_b(ntypx,5)
REAL(DP), ALLOCATABLE :: q_nq_zg(:, :) ! 3, nq
COMPLEX(DP), ALLOCATABLE :: z_nq_zg(:, :, :) ! nomdes, nmodes, nq
!
INTEGER :: nrots, kres1, kres2, col1, col2, Np
REAL(DP) :: kmin, kmax
!
NAMELIST /input/ flfrc, amass, asr, at, ntyp, loto_2d, &
NAMELIST /input/ flfrc, amass, asr, at, ntyp, loto_2d, loto_disable, &
& q_in_band_form, q_in_cryst_coord, point_label_type, &
& na_ifc, fd, &
! we add the inputs for generating the ZG-configuration
& ZG_conf, dim1, dim2, dim3, niters, error_thresh, q_external, &
& compute_error, synch, atm_zg, T, incl_qA, single_ph_displ, &
& atmsf_zg_a, atmsf_zg_b, ZG_strf, qpts_strf
! ZG_conf --> IF TRUE compute the ZG_configuration
! Last line of inputs are for the structure factor calculation
& compute_error, synch, atm_zg, T, incl_qA, single_ph_displ, ZG_strf
NAMELIST /strf_ZG/ atmsf_a, atmsf_b, qpts_strf, &
nrots, kres1, kres2, kmin, kmax, col1, col2, Np
! Last line of inputs are for the ZG structure factor calculation
!
CALL mp_startup()
CALL environment_start('ZG')
@ -272,7 +279,10 @@ PROGRAM ZG
q_in_band_form = .FALSE.
q_in_cryst_coord = .FALSE.
point_label_type = 'SC'
na_ifc = .FALSE.
fd = .FALSE.
loto_2d = .FALSE.
loto_disable = .FALSE.
!
ZG_conf = .TRUE.
compute_error = .TRUE.
@ -288,25 +298,40 @@ PROGRAM ZG
niters = 15000
atm_zg = "Element"
ZG_strf = .FALSE.
!
nrots = 1
kres1 = 250
kres2 = 250
kmin = -5
kmax = 10
col1 = 1
col2 = 2
Np = 100
qpts_strf = 0
atmsf_zg_a = 0.d0
atmsf_zg_b = 0.d0
atmsf_a = 0.d0
atmsf_b = 0.d0
!
!
!
IF (ionode) READ (5, input,IOSTAT=ios)
IF (ionode) READ (5, input, IOSTAT = ios)
CALL mp_bcast(ios, ionode_id, world_comm)
CALL errore('ZG', 'reading input namelist', ABS(ios))
IF ((ionode) .AND. (ZG_strf)) READ (5, strf_ZG , IOSTAT = ios)
CALL mp_bcast(ios, ionode_id, world_comm)
CALL errore('strf_ZG', 'reading strf_ZG namelist', ABS(ios))
CALL mp_bcast(asr, ionode_id, world_comm)
CALL mp_bcast(flfrc, ionode_id, world_comm)
CALL mp_bcast(amass, ionode_id, world_comm)
CALL mp_bcast(amass_blk, ionode_id, world_comm)
CALL mp_bcast(at, ionode_id, world_comm)
CALL mp_bcast(ntyp, ionode_id, world_comm)
CALL mp_bcast(na_ifc,ionode_id, world_comm)
CALL mp_bcast(fd,ionode_id, world_comm)
CALL mp_bcast(q_in_band_form, ionode_id, world_comm)
CALL mp_bcast(q_in_cryst_coord, ionode_id, world_comm)
CALL mp_bcast(point_label_type, ionode_id, world_comm)
CALL mp_bcast(loto_2d, ionode_id, world_comm)
CALL mp_bcast(loto_disable,ionode_id, world_comm)
!
CALL mp_bcast(ZG_conf, ionode_id, world_comm)
CALL mp_bcast(compute_error, ionode_id, world_comm)
@ -322,10 +347,21 @@ PROGRAM ZG
CALL mp_bcast(niters, ionode_id, world_comm)
CALL mp_bcast(atm_zg, ionode_id, world_comm)
CALL mp_bcast(ZG_strf, ionode_id, world_comm)
!
CALL mp_bcast(qpts_strf, ionode_id, world_comm)
CALL mp_bcast(atmsf_zg_a, ionode_id, world_comm)
CALL mp_bcast(atmsf_zg_b, ionode_id, world_comm)
CALL mp_bcast(atmsf_a, ionode_id, world_comm)
CALL mp_bcast(atmsf_b, ionode_id, world_comm)
CALL mp_bcast(nrots, ionode_id, world_comm)
CALL mp_bcast(kres1, ionode_id, world_comm)
CALL mp_bcast(kres2, ionode_id, world_comm)
CALL mp_bcast(kmin, ionode_id, world_comm)
CALL mp_bcast(kmax, ionode_id, world_comm)
CALL mp_bcast(col1, ionode_id, world_comm)
CALL mp_bcast(col2, ionode_id, world_comm)
CALL mp_bcast(Np, ionode_id, world_comm)
!
IF (loto_2d .AND. loto_disable) CALL errore('ZG', &
'loto_2d and loto_disable cannot be both true', 1)
! To check that user specifies supercell dimensions
IF (ZG_conf) THEN
IF ((dim1 < 1) .OR. (dim2 < 1) .OR. (dim3 < 1)) CALL errore('ZG', 'reading supercell size', dim1)
@ -540,6 +576,8 @@ PROGRAM ZG
ALLOCATE ( dyn(3, 3, nat, nat), dyn_blk(3, 3, nat_blk, nat_blk) )
ALLOCATE ( z(3 * nat, 3 * nat), w2(3 * nat, nq), f_of_q(3, 3, nat, nat) )
!
! Have to initialize w2
w2 = 0.d0
IF (ZG_conf) THEN
ALLOCATE ( z_nq_zg(3 * nat, 3 * nat, nq), q_nq_zg(3, nq))
z_nq_zg(:, :, :) = (0.d0, 0.d0)
@ -554,8 +592,6 @@ PROGRAM ZG
num_rap_mode=- 1
high_sym=.TRUE.
!
! Have to initialize w2
w2 = 0.d0
CALL fkbounds( nq, lower_bnd, upper_bnd )
!
DO n = lower_bnd, upper_bnd ! 1, nq
@ -563,11 +599,23 @@ PROGRAM ZG
lo_to_split = .FALSE.
f_of_q(:, :, :, :) = (0.d0, 0.d0)
IF(na_ifc) THEN
qq=SQRT(q(1,n)**2+q(2,n)**2+q(3,n)**2)
if(ABS(qq) < 1d-8) qq= 1.0
qhat(1)=q(1,n)/qq
qhat(2)=q(2,n)/qq
qhat(3)=q(3,n)/qq
CALL nonanal_ifc (nat,nat_blk,itau_blk,epsil,qhat,zeu,omega,dyn, &
nr1, nr2, nr3,f_of_q)
ENDIF
CALL setupmat (q(1, n), dyn, nat, at, bg, tau, itau_blk, nsc, alat, &
dyn_blk, nat_blk, at_blk, bg_blk, tau_blk, omega_blk, &
loto_2d, &
epsil, zeu, frc, nr1,nr2,nr3, has_zstar, rws, nrws, f_of_q)
epsil, zeu, frc, nr1,nr2,nr3, has_zstar, rws, nrws, na_ifc, f_of_q, fd)
IF (.not.loto_2d) THEN
qhat(1) = q(1, n) * at(1, 1) + q(2, n) * at(2, 1) + q(3, n) * at(3, 1)
@ -608,11 +656,14 @@ PROGRAM ZG
CALL infomsg &
('ZG','Z* not found in file '//TRIM(flfrc)// &
', TO-LO splitting at q = 0 will be absent!')
ELSEIF (loto_disable) THEN
CALL infomsg('ZG', &
'loto_disable is true. Disable LO-TO splitting at q=0.')
ELSE
lo_to_split=.TRUE.
ENDIF
!
CALL nonanal (nat, nat_blk, itau_blk, epsil, qhat, zeu, omega, dyn)
IF (lo_to_split) CALL nonanal (nat, nat_blk, itau_blk, epsil, qhat, zeu, omega, dyn)
!
ENDIF
!
@ -660,9 +711,9 @@ PROGRAM ZG
synch, tau, alat, atm_zg, ntypx, at, &
q_in_cryst_coord, q_external, T, incl_qA, &
compute_error, single_ph_displ, &
ZG_strf, qpts_strf, atmsf_zg_a, atmsf_zg_b)
ZG_strf, qpts_strf, atmsf_a, atmsf_b, &
nrots, kres1, kres2, kmin, kmax, col1, col2, Np)
!
!
DEALLOCATE (z, w2, dyn, dyn_blk)
!
IF (ZG_conf) DEALLOCATE (z_nq_zg, q_nq_zg)
@ -813,7 +864,8 @@ SUBROUTINE readfc ( flfrc, nr1, nr2, nr3, epsil, nat, &
END SUBROUTINE readfc
!
!-----------------------------------------------------------------------
SUBROUTINE frc_blk(dyn,q,tau, nat, nr1, nr2, nr3,frc, at, bg, rws, nrws,f_of_q)
SUBROUTINE frc_blk(dyn,q,tau, nat, nr1, nr2, nr3, frc, &
at, bg, rws, nrws, f_of_q, fd)
!-----------------------------------------------------------------------
! calculates the dynamical matrix at q from the (short-range part of the)
! force constants
@ -833,6 +885,7 @@ SUBROUTINE frc_blk(dyn,q,tau, nat, nr1, nr2, nr3,frc, at, bg, rws, nrws,f_of_q)
REAL(DP),SAVE,ALLOCATABLE :: wscache(:, :, :, :, :)
REAL(DP), ALLOCATABLE :: ttt(:, :, :, :, :), tttx(:, :)
LOGICAL,SAVE :: first=.TRUE.
LOGICAL :: fd
!
nr1_=2*nr1
nr2_=2*nr2
@ -850,6 +903,7 @@ SUBROUTINE frc_blk(dyn,q,tau, nat, nr1, nr2, nr3,frc, at, bg, rws, nrws,f_of_q)
DO i = 1, 3
r(i) = n1*at(i, 1) +n2*at(i, 2) +n3*at(i, 3)
r_ws(i) = r(i) + tau(i, na) -tau(i, nb)
if (fd) r_ws(i) = r(i) + tau(i,nb)-tau(i,na)
ENDDO
wscache(n3, n2, n1, nb, na) = wsweight(r_ws, rws, nrws)
ENDDO
@ -924,7 +978,7 @@ END SUBROUTINE frc_blk
SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, &
& dyn_blk,nat_blk,at_blk,bg_blk,tau_blk,omega_blk, &
& loto_2d, &
& epsil,zeu,frc,nr1,nr2,nr3,has_zstar,rws,nrws,f_of_q)
& epsil,zeu,frc,nr1,nr2,nr3,has_zstar,rws,nrws,na_ifc,f_of_q,fd)
!-----------------------------------------------------------------------
! compute the dynamical matrix (the analytic part only)
!
@ -944,7 +998,7 @@ SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, &
REAL(DP) :: tau_blk(3, nat_blk), at_blk(3, 3), bg_blk(3, 3), omega_blk
COMPLEX(DP) dyn_blk(3, 3, nat_blk, nat_blk), f_of_q(3, 3, nat, nat)
COMPLEX(DP) :: dyn(3, 3, nat, nat)
LOGICAL :: has_zstar
LOGICAL :: has_zstar, na_ifc, fd
!
! local variables
!
@ -965,8 +1019,8 @@ SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, &
!
dyn_blk(:, :, :, :) = (0.d0,0.d0)
CALL frc_blk (dyn_blk, qp,tau_blk, nat_blk, &
& nr1, nr2, nr3,frc, at_blk, bg_blk, rws, nrws,f_of_q)
IF (has_zstar) &
& nr1, nr2, nr3,frc, at_blk, bg_blk, rws, nrws,f_of_q,fd)
IF (has_zstar .and. .not. na_ifc) &
CALL rgd_blk(nr1, nr2, nr3, nat_blk, dyn_blk, qp,tau_blk, &
epsil, zeu, bg_blk, omega_blk, celldm(1), loto_2d, +1.d0)
!
@ -1013,10 +1067,10 @@ SUBROUTINE set_asr (asr, nr1, nr2, nr3, frc, zeu, nat, ibrav, tau)
USE io_global, ONLY : ionode, stdout
!
IMPLICIT NONE
CHARACTER (LEN= 10), intent(in) :: asr
INTEGER, intent(in) :: nr1, nr2, nr3, nat, ibrav
REAL(DP), intent(in) :: tau(3, nat)
REAL(DP), intent(inout) :: frc(nr1, nr2, nr3, 3, 3, nat, nat), zeu(3, 3, nat)
CHARACTER (LEN= 10), INTENT(in) :: asr
INTEGER, INTENT(in) :: nr1, nr2, nr3, nat, ibrav
REAL(DP), INTENT(in) :: tau(3, nat)
REAL(DP), INTENT(inout) :: frc(nr1, nr2, nr3, 3, 3, nat, nat), zeu(3, 3, nat)
!
INTEGER :: axis, n, i, j, na, nb, n1, n2, n3, m, p, k,l,q, r, i1, j1, na1
REAL(DP) :: zeu_new(3, 3, nat)
@ -1858,7 +1912,7 @@ END SUBROUTINE read_tau
!-----------------------------------------------------------------------
subroutine setgam (q, gam, nat, at, bg,tau, itau_blk, nsc, alat, &
& gam_blk, nat_blk, at_blk, bg_blk,tau_blk, omega_blk, &
& frcg, nr1, nr2, nr3, rws, nrws)
& frcg, nr1, nr2, nr3, rws, nrws, fd)
!-----------------------------------------------------------------------
! compute the dynamical matrix (the analytic part only)
!
@ -1873,8 +1927,9 @@ subroutine setgam (q, gam, nat, at, bg,tau, itau_blk, nsc, alat, &
REAL(DP) :: q(3), tau(3, nat), at(3, 3), bg(3, 3), alat, rws(0:3, nrws)
REAL(DP) :: tau_blk(3, nat_blk), at_blk(3, 3), bg_blk(3, 3), omega_blk, &
frcg(nr1, nr2, nr3, 3, 3, nat_blk, nat_blk)
COMPLEX(DP) :: gam_blk(3, 3, nat_blk, nat_blk),f_of_q(3, 3, nat, nat)
COMPLEX(DP) :: gam(3, 3, nat, nat)
COMPLEX(DP) :: gam_blk(3, 3, nat_blk, nat_blk),f_of_q(3, 3, nat, nat)
COMPLEX(DP) :: gam(3, 3, nat, nat)
LOGICAL :: fd
!
! local variables
!
@ -1895,7 +1950,7 @@ subroutine setgam (q, gam, nat, at, bg,tau, itau_blk, nsc, alat, &
!
gam_blk(:, :, :, :) = (0.d0,0.d0)
CALL frc_blk (gam_blk, qp,tau_blk, nat_blk, &
nr1, nr2, nr3,frcg, at_blk, bg_blk, rws, nrws,f_of_q)
nr1, nr2, nr3,frcg, at_blk, bg_blk, rws, nrws,f_of_q, fd)
!
DO na= 1, nat
na_blk = itau_blk(na)
@ -2062,8 +2117,8 @@ subroutine readfg ( ifn, nr1, nr2, nr3, nat, frcg )
USE mp_world, ONLY : world_comm
implicit none
! I/O variable
integer, intent(in) :: nr1, nr2, nr3, nat
REAL(DP), intent(out) :: frcg(nr1, nr2, nr3, 3, 3, nat, nat)
integer, INTENT(in) :: nr1, nr2, nr3, nat
REAL(DP), INTENT(out) :: frcg(nr1, nr2, nr3, 3, 3, nat, nat)
! local variables
integer i, j, na, nb, m1,m2,m3, ifn
integer ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid
@ -2170,9 +2225,9 @@ SUBROUTINE qpoint_gen1(dim1, dim2, dim3, ctrAB)
IMPLICIT NONE
! input
INTEGER, intent(in) :: dim1, dim2, dim3
INTEGER, intent(out) :: ctrAB
!! REAL(DP), intent(out) :: q_AB(:,:)
INTEGER, INTENT(in) :: dim1, dim2, dim3
INTEGER, INTENT(out) :: ctrAB
!! REAL(DP), INTENT(out) :: q_AB(:,:)
! local
INTEGER :: i, j, k, n, nqs
INTEGER :: lower_bnd, upper_bnd
@ -2240,8 +2295,8 @@ SUBROUTINE qpoint_gen2(dim1, dim2, dim3, ctrAB, q_AB)
IMPLICIT NONE
! input
INTEGER, intent(in) :: dim1, dim2, dim3, ctrAB
REAL(DP), intent(out) :: q_AB(3, ctrAB)
INTEGER, INTENT(in) :: dim1, dim2, dim3, ctrAB
REAL(DP), INTENT(out) :: q_AB(3, ctrAB)
! local
INTEGER :: i, j, k, n, ctr, nqs
INTEGER :: lower_bnd, upper_bnd
@ -2319,11 +2374,12 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
dim1, dim2, dim3, niters, error_thresh, synch, tau, alat, atm, &
ntypx, at, q_in_cryst_coord, q_external, T, incl_qA, &
compute_error, single_ph_displ, &
ZG_strf, qpts_strf, atmsf_zg_a, atmsf_zg_b)
ZG_strf, qpts_strf, atmsf_a, atmsf_b, &
nrots, kres1, kres2, kmin, kmax, col1, col2, Np)
!
use kinds, only: dp
use constants, only: amu_ry, ry_to_thz, ry_to_cmm1, H_PLANCK_SI, &
K_BOLTZMANN_SI, AMU_SI, pi
USE kinds, ONLY : dp
USE constants, ONLY : amu_ry, ry_to_thz, ry_to_cmm1, H_PLANCK_SI, &
K_BOLTZMANN_SI, AMU_SI, pi
USE cell_base, ONLY : bg
USE io_global, ONLY : ionode, ionode_id, stdout
USE mp_world, ONLY : world_comm
@ -2331,17 +2387,19 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
IMPLICIT NONE
! input
CHARACTER(LEN=3), intent(in) :: atm(ntypx)
LOGICAL, intent(in) :: synch, q_in_cryst_coord, q_external, ZG_strf
LOGICAL, intent(in) :: incl_qA, compute_error, single_ph_displ
INTEGER, intent(in) :: dim1, dim2, dim3, niters, qpts_strf
INTEGER, intent(in) :: nq, nat, ntyp, ios, ntypx
INTEGER, intent(in) :: ityp(nat)
CHARACTER(LEN=3), INTENT(in) :: atm(ntypx)
LOGICAL, INTENT(in) :: synch, q_in_cryst_coord, q_external, ZG_strf
LOGICAL, INTENT(in) :: incl_qA, compute_error, single_ph_displ
INTEGER, INTENT(in) :: dim1, dim2, dim3, niters, qpts_strf
INTEGER, INTENT(in) :: nq, nat, ntyp, ios, ntypx
! nq is the number of qpoints in sets A and B
REAL(DP), intent(in) :: error_thresh, alat, T
REAL(DP), intent(in) :: at(3, 3), atmsf_zg_a(ntypx,5), atmsf_zg_b(ntypx,5)
REAL(DP), intent(in) :: q(3, nq), w2(3 * nat, nq), amass(ntyp), tau(3, nat)
COMPLEX(DP), intent(in) :: z_nq_zg(3 * nat, 3 * nat, nq)
INTEGER, INTENT(in) :: ityp(nat)
INTEGER, INTENT(in) :: nrots, kres1, kres2, col1, col2, Np
REAL(DP), INTENT(in) :: kmin, kmax
REAL(DP), INTENT(in) :: error_thresh, alat, T
REAL(DP), INTENT(in) :: at(3, 3), atmsf_a(ntypx,5), atmsf_b(ntypx,5)
REAL(DP), INTENT(in) :: q(3, nq), w2(3 * nat, nq), amass(ntyp), tau(3, nat)
COMPLEX(DP), INTENT(in) :: z_nq_zg(3 * nat, 3 * nat, nq)
!
! local
CHARACTER(len=256) :: filename
@ -2375,7 +2433,8 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
! for momenta/velocities
REAL(DP), ALLOCATABLE :: Cpx_matA(:, :), Cpx_matB(:, :), Cpx_matAB(:, :)
! matrices to account for the coupling terms between different phonon branches !
REAL(DP), ALLOCATABLE :: sum_error_D(:, :), sum_diag_D(:, :), sum_error_B(:), sum_diag_B(:), sum_error_B2(:), sum_diag_B2(:)
REAL(DP), ALLOCATABLE :: sum_error_D(:, :), sum_diag_D(:, :), sum_error_B(:)
REAL(DP), ALLOCATABLE :: sum_diag_B(:), sum_error_B2(:), sum_diag_B2(:)
REAL(DP), ALLOCATABLE :: D_tau(:, :, :), P_tau(:, :, :), ratio_zg(:)! displacements and velocities
REAL(DP), ALLOCATABLE :: R_mat(:, :), E_vect(:, :), D_vect(:, :), F_vect(:, :)
! D_tau : atomic displacements
@ -2748,13 +2807,14 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
! i.e. [1 1 1 1 1 1] gives same result to [- 1 -1 -1 -1 -1 -1].
!
ALLOCATE(M_mat(2 * pn, nat3), Mx_mat(pn, nat3), Mx_mat_or(pn, nat3), V_mat(2))
M_mat = 1 ! initialize M_mat
V_mat = (/ 1, -1/) ! initialize V_mat whose entries will generate the sign matrices
DO i = 1, nat3
ctr = 1
DO p = 1, 2**(i - 1)
DO qp = 1, 2
DO k = 1, 2**(nat3 - i)
IF (ctr > 2 * pn) EXIT ! I DO this in case there many branches in the system and
IF (ctr > 2 * pn) EXIT ! in case there many branches in the system and
! in that case we DO not need to ALLOCATE more signs
M_mat(ctr, i) = V_mat(qp)
ctr = ctr + 1
@ -2819,7 +2879,7 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
DO kk = 1, niters
! Allocate original matrices ! half the entries of M_mat
! We also make the inherent choice that each column of Mx_mat_or has the same number of positive and negative signs
Mx_mat_or = 0
Mx_mat_or = 1
DO i = 1, 2 * pn / 4, 2
Mx_mat_or(i, :) = M_mat(i, :)
ENDDO
@ -3140,8 +3200,8 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
IF (ionode) THEN
DO p = 1, nq_tot
DO k = 1, nat ! k represents the atom
WRITE(80,'(A6, 3F13.8)') atm(ityp(k)), D_tau(p, k, :)
WRITE(*,'(A10, A6, 3F13.8)') "ZG_conf:", atm(ityp(k)), D_tau(p, k, :)
WRITE(80,'(A6, 3F16.8)') atm(ityp(k)), D_tau(p, k, :)
WRITE(*,'(A10, A6, 3F16.8)') "ZG_conf:", atm(ityp(k)), D_tau(p, k, :)
WRITE(81,'(A6, 3F15.8)') atm(ityp(k)), P_tau(p, k, :) * 1.0E-12 ! multiply to obtain picoseconds
ENDDO
ENDDO
@ -3251,7 +3311,8 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
WRITE(*,*) "Computing ZG structure factor ..."
IF ( ZG_strf .AND. ( SUM(ABS(ratio_zg)) / nat3 < error_thresh) ) &
call ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
nat, alat, ityp, ntypx, atmsf_zg_a, atmsf_zg_b)
nat, alat, ityp, ntypx, atmsf_a, atmsf_b, &
nrots, kres1, kres2, kmin, kmax, col1, col2, Np)
!
IF (ionode) WRITE(*,*)
IF (SUM(ABS(ratio_zg)) / nat3 > error_thresh .AND. ionode ) &
@ -3272,7 +3333,8 @@ SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, &
END SUBROUTINE ZG_configuration
SUBROUTINE ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
nat, alat, ityp, ntypx, atmsf_zg_a, atmsf_zg_b)
nat, alat, ityp, ntypx, atmsf_a, atmsf_b, &
nrots, kres1, kres2, kmin, kmax, col1, col2, Np)
!
USE kinds, ONLY : DP
USE cell_base, ONLY : bg
@ -3284,34 +3346,36 @@ SUBROUTINE ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
!
IMPLICIT NONE
!
INTEGER, intent(in) :: nq_tot, nat, qpts_strf, ntypx
REAL(DP), intent(in) :: D_tau(nq_tot, nat, 3), equil_p(nq_tot, nat, 3), alat
REAL(DP), intent(in) :: atmsf_zg_a(ntypx, 5), atmsf_zg_b(ntypx, 5)
INTEGER, INTENT(in) :: nq_tot, nat, qpts_strf, ntypx
INTEGER, INTENT(in) :: nrots, kres1, kres2, col1, col2, Np
REAL(DP), INTENT(in) :: kmin, kmax
REAL(DP), INTENT(in) :: D_tau(nq_tot, nat, 3), equil_p(nq_tot, nat, 3), alat
REAL(DP), INTENT(in) :: atmsf_a(ntypx, 5), atmsf_b(ntypx, 5)
INTEGER :: i, k, p, j, kk, pp, nta, ii
INTEGER :: lower_bnd, upper_bnd
INTEGER :: lower_bnd, upper_bnd, ctr
INTEGER ityp(nat)
REAL(DP) :: qpts_strf_list(3, qpts_strf), atomic_form_factor(nat, qpts_strf)
REAL(DP) :: dotp, dotpp
REAL(DP) :: q_strf(3, qpts_strf), atomic_form_factor(nat, qpts_strf)
REAL(DP) :: dotp, dotpp, eps
REAL(DP), ALLOCATABLE :: strf_rot(:, :)
COMPLEX(DP) :: imagi, strf_map(qpts_strf)
!
!
eps = 1d-5
imagi = (0.0d0, 1.0d0) !imaginary unit
!
!
qpts_strf_list = 0.d0
q_strf = 0.d0
IF (ionode) THEN
OPEN (unit = 99, file = './qpts_strf.dat', status = 'unknown', form = 'formatted')
!
DO k = 1, qpts_strf
READ(99,*) qpts_strf_list(:, k)
! WRITE(*,*) "aa", qpts_strf_list(k,:)
CALL cryst_to_cart(1, qpts_strf_list(:, k), bg, +1)
qpts_strf_list(:, k) = qpts_strf_list(:, k) * ( 2.d0 * pi / alat / BOHR_RADIUS_ANGS ) ! / 0.138933 * 0.073520
READ(99,*) q_strf(:, k)
! WRITE(*,*) "aa", q_strf(k,:)
CALL cryst_to_cart(1, q_strf(:, k), bg, +1)
q_strf(:, k) = q_strf(:, k) * ( 2.d0 * pi / alat / BOHR_RADIUS_ANGS ) ! / 0.138933 * 0.073520
ENDDO
!
CLOSE(99)
ENDIF
CALL mp_bcast(qpts_strf_list, ionode_id, world_comm)
CALL mp_bcast(q_strf, ionode_id, world_comm)
!
!
atomic_form_factor = 0.d0
@ -3320,8 +3384,8 @@ SUBROUTINE ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
nta = ityp(k)
DO i = 1, qpts_strf
DO ii = 1, 5
atomic_form_factor(k, i) = atomic_form_factor(k, i) + (atmsf_zg_a(nta, ii) * &
EXP(-atmsf_zg_b(nta, ii) * (NORM2(qpts_strf_list(:, i)) / 4.d0 / pi)**2))
atomic_form_factor(k, i) = atomic_form_factor(k, i) + (atmsf_a(nta, ii) * &
EXP(-atmsf_b(nta, ii) * (NORM2(q_strf(:, i)) / 4.d0 / pi)**2))
ENDDO
ENDDO
ENDDO
@ -3336,8 +3400,8 @@ SUBROUTINE ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
! DO pp = 1, nq_tot !!
dotp = 0.0d0
DO j = 1, 3
dotp = dotp + qpts_strf_list(j, i) * D_tau(p, k, j)
! dotp = dotp + qpts_strf_list(j, i) * (D_tau(p, k, j) - D_tau(pp, !kk, j)) !!
dotp = dotp + q_strf(j, i) * D_tau(p, k, j)
! dotp = dotp + q_strf(j, i) * (D_tau(p, k, j) - D_tau(pp, kk, j)) !!
ENDDO
strf_map(i) = strf_map(i) + EXP(imagi * dotp) * atomic_form_factor(k, i)
!strf_map(i) = strf_map(i) + EXP(imagi * dotp) * atomic_form_factor(k, i) * atomic_form_factor(kk, i)
@ -3349,12 +3413,32 @@ SUBROUTINE ZG_structure_factor(qpts_strf, D_tau, equil_p, nq_tot, &
CALL mp_sum(strf_map, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
! APPLY BROADENING and print outputs
ctr = 0
DO k = 1, qpts_strf
IF ((q_strf(col1, k) .GT. 0.d0 - eps) .AND. &
(q_strf(col2, k) .GT. 0.d0 - eps)) THEN
ctr = ctr + 1
ENDIF
ENDDO
ALLOCATE(strf_rot(ctr * nrots, 4))
!
strf_rot = 0.d0
IF (ionode) CALL rotate(DBLE(ABS(strf_map)**2), q_strf, qpts_strf, 0, nrots, &
ctr, strf_rot, col1, col2)
!
CALL mp_bcast(strf_rot, ionode_id, world_comm)
CALL disca_broadening(strf_rot, ctr * nrots, kres1, kres2, alat, &
kmin, kmax, col1, col2, Np, 'strf_ZG_broad.dat')
!
DEALLOCATE(strf_rot)
!
IF (ionode) THEN
OPEN (unit = 98, file = './structure_factor_ZG.dat', status = 'unknown', form = 'formatted')
OPEN (unit = 98, file = './structure_factor_ZG_raw.dat', status = 'unknown', form = 'formatted')
!
DO i = 1, qpts_strf
!
WRITE(98,'(4f26.6)') qpts_strf_list(:, i), abs(strf_map(i))**2
WRITE(98,'(4f26.6)') q_strf(:, i), ABS(strf_map(i))**2
!abs(strf_map(i))**2 ! real(strf_map(i)) !
!
ENDDO
@ -3374,9 +3458,9 @@ SUBROUTINE create_supercell(at,tau, alat, dim1, dim2, dim3, nat, equil_p)
IMPLICIT NONE
!
!
INTEGER, intent(in) :: dim1, dim2, dim3, nat
REAL(DP), intent(in) :: tau(3, nat), at(3, 3), alat
REAL(DP), intent(out) :: equil_p(dim1 * dim2 * dim3, nat, 3)
INTEGER, INTENT(in) :: dim1, dim2, dim3, nat
REAL(DP), INTENT(in) :: tau(3, nat), at(3, 3), alat
REAL(DP), INTENT(out) :: equil_p(dim1 * dim2 * dim3, nat, 3)
INTEGER :: i, j, k, ctr, p
REAL(DP) :: alat_ang, crystal_pos(dim1 * dim2 * dim3, nat, 3), abc(3, nat)
alat_ang = alat * BOHR_RADIUS_ANGS !bohr_to_angst ! to convert them in angstrom !
@ -3436,13 +3520,13 @@ SUBROUTINE single_phonon(nq_tot, nat, ctrB, ctrA, nat3, ityp, ntyp, &
IMPLICIT NONE
!
!
CHARACTER(LEN=3), intent(in) :: atm(ntypx)
INTEGER, intent(in) :: nq_tot, nat, ctrB, ctrA, nat3, ntyp, ntypx
INTEGER, intent(in) :: Rlist(nq_tot, 3)
INTEGER, intent(in) :: ityp(nat)
REAL(DP), intent(in) :: qA(ctrA, 3), qB(ctrB, 3), amass(ntyp), equil_p(nq_tot, nat, 3)
REAL(DP), intent(in) :: Cx_matB(nat3, ctrB), Cx_matA(nat3, ctrA), Cpx_matA(nat3, ctrA), Cpx_matB(nat3, ctrB)
COMPLEX(DP), intent(in) :: z_nq_A(nat3, nat3, ctrA), z_nq_B(nat3, nat3, ctrB)
CHARACTER(LEN=3), INTENT(in) :: atm(ntypx)
INTEGER, INTENT(in) :: nq_tot, nat, ctrB, ctrA, nat3, ntyp, ntypx
INTEGER, INTENT(in) :: Rlist(nq_tot, 3)
INTEGER, INTENT(in) :: ityp(nat)
REAL(DP), INTENT(in) :: qA(ctrA, 3), qB(ctrB, 3), amass(ntyp), equil_p(nq_tot, nat, 3)
REAL(DP), INTENT(in) :: Cx_matB(nat3, ctrB), Cx_matA(nat3, ctrA), Cpx_matA(nat3, ctrA), Cpx_matB(nat3, ctrB)
COMPLEX(DP), INTENT(in) :: z_nq_A(nat3, nat3, ctrA), z_nq_B(nat3, nat3, ctrB)
CHARACTER(len=256) :: filename
!
INTEGER :: i, j, k, p, ii, ctr, nta, qp
@ -3600,3 +3684,145 @@ SUBROUTINE fkbounds( nktot, lower_bnd, upper_bnd )
RETURN
!
END SUBROUTINE fkbounds
!
SUBROUTINE rotate(strf, q, nq, nq_super, nrots, &
ctr, strf_rot, col1, col2)
!
USE kinds, ONLY : DP
USE constants, ONLY : tpi
USE io_global, ONLY : stdout
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: nq, nrots, col1, col2, ctr, nq_super
REAL(DP), INTENT(in) :: strf(nq), q(3, nq)
REAL(DP), INTENT(out) :: strf_rot(ctr * nrots, 4)
INTEGER :: i, p, ctr2
REAL(DP) :: str_f(ctr,4), str_fp(ctr, 4)
!
REAL(DP) :: Rmat(2,2), theta, eps
!
eps = 1d-5
!
!WRITE(*,*) "Number of rotations provided:", nrots
theta = tpi / FLOAT(nrots)
!
!
str_f = 0.d0
ctr2 = 0
DO p = nq_super + 1, nq
IF ((q(col1, p) .GT. 0.0 - eps) .AND. (q(col2, p) .GT. 0.0 - eps)) THEN
ctr2 = ctr2 + 1
str_f(ctr2, 1:3) = q(:, p)
str_f(ctr2, 4) = strf(p)
ENDIF
ENDDO
!
! To remove double contribution upon rotation
!
str_fp = 0.d0
str_fp(1, :) = str_f(1, :)
DO p = 2, ctr
IF (ATAN(str_f(p, 1) / str_f(p, 2)) .LT. ( tpi / float(nrots) - eps) ) THEN
str_fp(p, :) = str_f(p, :)
ENDIF
ENDDO
!
!
strf_rot = 0.d0
ctr2 = 1
DO i = 0, nrots - 1
Rmat(1, :) = (/ COS(i * theta), -SIN(i * theta) /)
Rmat(2, :) = (/ SIN(i * theta), COS(i * theta) /)
DO p = 1, ctr
strf_rot(ctr2, 1 : 2) = MATMUL(Rmat, str_fp(p, 1 : 2))
strf_rot(ctr2, 3) = str_fp(p, 3)
strf_rot(ctr2, 4) = str_fp(p, 4)
ctr2 = ctr2 + 1
END DO
ENDDO
END SUBROUTINE
!
SUBROUTINE disca_broadening(strf_rot, steps, kres1, kres2, alat, &
kmin, kmax, col1, col2, Np, flstrfout)
!-------------------------------------------------------------------------
!! authors: Marios Zacharias, Feliciano Giustino
!! acknowledgement: Hyungjun Lee for help packaging this release
!! version: v0.1
!! license: GNU
!
USE kinds, ONLY : dp
USE mp_global, ONLY : inter_pool_comm
USE mp_world, ONLY : world_comm
USE mp, ONLY : mp_bcast, mp_barrier, mp_sum
USE io_global, ONLY : ionode, ionode_id, stdout
USE constants, ONLY : pi
!
IMPLICIT NONE
!
CHARACTER(LEN=256), INTENT(IN) :: flstrfout
INTEGER, INTENT(IN) :: steps, kres1, kres2, Np, col1, col2
REAL(DP), INTENT(IN) :: kmin, kmax, alat
REAL(DP), INTENT(IN) :: strf_rot(steps, 4)
INTEGER :: ii, lower_bnd, upper_bnd, ik, iky
REAL(DP) :: jump, sf_smearingx, sf_smearingy, maxv !, pi
REAL(DP), ALLOCATABLE :: kgridx(:), kgridy(:), strf_out(:, :)
!
!
!
ALLOCATE( kgridx(kres1), kgridy(kres2))
ALLOCATE(strf_out(kres1,kres2))
!kmin = -10.0
!kmax = 10.0
jump = (kmax - kmin) / DBLE(kres1 - 1)
DO ik = 1, kres1
kgridx(ik) = kmin + (ik - 1) * jump
ENDDO
sf_smearingx = (kmax - kmin) / DBLE(kres1)
!
jump = (kmax - kmin) / DBLE(kres2 - 1)
DO ik = 1, kres2
kgridy(ik) = kmin + (ik- 1)*jump
ENDDO
sf_smearingy = (kmax - kmin) / DBLE(kres2)
!!
!
strf_out = 0.d0
!
CALL fkbounds( steps, lower_bnd, upper_bnd )
!
DO ii = lower_bnd, upper_bnd
DO ik = 1, kres1 !
DO iky = 1, kres2
!
strf_out(ik, iky) = strf_out(ik, iky) + &
strf_rot(ii, 4) / sf_smearingx / SQRT(2.0d0 * pi) / sf_smearingy / SQRT(2.0d0 * pi)* &
(EXP(-(strf_rot(ii, col1) - kgridx(ik))**2.d0 / sf_smearingx**2.d0 / 2.d0))*&
(EXP(-(strf_rot(ii, col2) - kgridy(iky))**2.d0 / sf_smearingy**2.d0 / 2.d0))
!
ENDDO
ENDDO
ENDDO
!
CALL mp_sum(strf_out, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
IF (ionode) THEN
OPEN(46,FILE=flstrfout)
maxv = maxval(strf_out)
WRITE(46,*) "#", maxv, maxval(strf_out)
DO ik = 1, kres1
DO iky = 1, kres2
WRITE(46,'(3F28.12)') kgridx(ik), kgridy(iky), strf_out(ik,iky) * Np**(-2.0d0) !
!Np**(-2.0d0) ! / maxv
ENDDO
WRITE(46,*)
ENDDO
CLOSE(46)
ENDIF
!
DEALLOCATE(strf_out, kgridx, kgridy)
!
!
END SUBROUTINE

1780
EPW/ZG/src/disca.f90
View File

File diff suppressed because it is too large Load Diff

58
EPW/ZG/src/pp_disca.f90
View File

@ -17,14 +17,14 @@ USE constants, ONLY : pi
IMPLICIT NONE
!
CHARACTER(LEN=256) :: flstrfin, flstrfout
INTEGER :: ksteps1, ksteps2, ik, iky, dim1, dim2
INTEGER :: kres1, kres2, ik, iky, col1, col2
INTEGER :: steps, ii, lower_bnd, upper_bnd, Np, ios
REAL(DP) :: kmin, kmax, maxv
REAL(DP) :: jump, sf_smearingx, sf_smearingy !, pi
REAL(DP), ALLOCATABLE :: kgridx(:), kgridy(:), structure_fact(:, :), structure_fact_out(:, :)
!
NAMELIST /input/ steps, ksteps1, ksteps2, kmin, kmax, &
& dim1, dim2, Np, flstrfin, flstrfout
NAMELIST /input/ steps, kres1, kres2, kmin, kmax, &
& col1, col2, Np, flstrfin, flstrfout
!
CALL mp_startup()
CALL environment_start('DISCA_BROADENING')
@ -36,12 +36,12 @@ CALL environment_start('DISCA_BROADENING')
! set namelist default
!
steps = 10000
ksteps1 = 250
ksteps2 = 250
kres1 = 250
kres2 = 250
kmin = -5
kmax = 5
dim1 = 1
dim2 = 2
col1 = 1
col2 = 2
Np = 400000
flstrfin = 'input_data.dat'
flstrfout = 'output_data.dat'
@ -50,12 +50,12 @@ CALL environment_start('DISCA_BROADENING')
CALL mp_bcast(ios, ionode_id, world_comm)
CALL errore('disca', 'reading input namelist', ABS(ios))
CALL mp_bcast(steps, ionode_id, world_comm)
CALL mp_bcast(ksteps1, ionode_id, world_comm)
CALL mp_bcast(ksteps2, ionode_id, world_comm)
CALL mp_bcast(kres1, ionode_id, world_comm)
CALL mp_bcast(kres2, ionode_id, world_comm)
CALL mp_bcast(kmin, ionode_id, world_comm)
CALL mp_bcast(kmax, ionode_id, world_comm)
CALL mp_bcast(dim1, ionode_id, world_comm)
CALL mp_bcast(dim2, ionode_id, world_comm)
CALL mp_bcast(col1, ionode_id, world_comm)
CALL mp_bcast(col2, ionode_id, world_comm)
CALL mp_bcast(Np, ionode_id, world_comm)
CALL mp_bcast(flstrfin, ionode_id, world_comm)
CALL mp_bcast(flstrfout, ionode_id, world_comm)
@ -63,31 +63,31 @@ CALL environment_start('DISCA_BROADENING')
!OPEN(46,FILE=flstrfin)
! READ(46,*) steps
! READ(46,*) ksteps1
! READ(46,*) ksteps2
! READ(46,*) kres1
! READ(46,*) kres2
! READ(46,*) kmin
! READ(46,*) kmax
! READ(46,*) dim1
! READ(46,*) dim2
! READ(46,*) col1
! READ(46,*) col2
! READ(46,*) Np
!CLOSE(46)
ALLOCATE(structure_fact(steps,4), kgridx(ksteps1), kgridy(ksteps2))
ALLOCATE(structure_fact_out(ksteps1,ksteps2))
ALLOCATE(structure_fact(steps,4), kgridx(kres1), kgridy(kres2))
ALLOCATE(structure_fact_out(kres1,kres2))
!kmin = -10.0
!kmax = 10.0
jump = (kmax - kmin) / dble(ksteps1 - 1)
DO ik = 1, ksteps1
jump = (kmax - kmin) / dble(kres1 - 1)
DO ik = 1, kres1
kgridx(ik) = kmin + (ik - 1) * jump
ENDDO
sf_smearingx = (kmax - kmin)/dble(ksteps1)
sf_smearingx = (kmax - kmin)/dble(kres1)
!
jump = (kmax - kmin) / dble(ksteps2 - 1)
DO ik = 1, ksteps2
jump = (kmax - kmin) / dble(kres2 - 1)
DO ik = 1, kres2
kgridy(ik) = kmin + (ik- 1)*jump
ENDDO
sf_smearingy = (kmax-kmin)/dble(ksteps2)
sf_smearingy = (kmax-kmin)/dble(kres2)
!!
OPEN(46, FILE = flstrfin)
DO ii = 1, steps, 1
@ -100,13 +100,13 @@ structure_fact_out = 0.d0
CALL fkbounds( steps, lower_bnd, upper_bnd )
!
DO ii = lower_bnd, upper_bnd
DO ik = 1, ksteps1 !
DO iky = 1, ksteps2
DO ik = 1, kres1 !
DO iky = 1, kres2
!
structure_fact_out(ik, iky) = structure_fact_out(ik, iky) + &
structure_fact(ii, 4) / sf_smearingx / sqrt(2.0d0 * pi) / sf_smearingy / sqrt(2.0d0 * pi)* &
(EXP(-(structure_fact(ii, dim1)- kgridx(ik))**2.d0 / sf_smearingx**2.d0 / 2.d0))*&
(EXP(-(structure_fact(ii, dim2)- kgridy(iky))**2.d0 / sf_smearingy**2.d0 / 2.d0))
(EXP(-(structure_fact(ii, col1)- kgridx(ik))**2.d0 / sf_smearingx**2.d0 / 2.d0))*&
(EXP(-(structure_fact(ii, col2)- kgridy(iky))**2.d0 / sf_smearingy**2.d0 / 2.d0))
!
ENDDO
ENDDO
@ -120,8 +120,8 @@ IF (ionode) THEN
! maxv = DBLE(6.162035539820569E+019)
maxv = maxval(structure_fact_out)
WRITE(46,*) "#", maxv, maxval(structure_fact_out)
DO ik = 1, ksteps1
DO iky = 1, ksteps2
DO ik = 1, kres1
DO iky = 1, kres2
WRITE(46,'(3F28.12)') kgridx(ik), kgridy(iky), structure_fact_out(ik,iky) * Np**(-2.0d0) !1.169035831194574E+019
!/ !maxval(abs(structure_fact_out)) !!* Np**(-2.0d0) ! / maxv
ENDDO

BIN
EPW/ZG/tutorial.tar.gz Normal file
View File

Binary file not shown.

58
EPW/examples/diamond/epw/epw4.in Normal file
View File

@ -0,0 +1,58 @@
--
&inputepw
prefix = 'diam'
amass(1) = 12.01078
outdir = './'
iverbosity = 0
ep_coupling = .false.
elph = .false.
epbwrite = .false.
epbread = .false.
epwwrite = .false.
epwread = .true.
efermi_read = .true.
fermi_energy = 13.209862
nbndsub = 1
wannierize = .false.
num_iter = 300
iprint = 2
dis_win_max = 12
dis_froz_max= 7
proj(1) = 'f=0,0,0:l=-3'
elecselfen = .true.
phonselfen = .false.
a2f = .false.
specfun_el = .true.
wmin_specfun = -4
wmax_specfun = 4
nw_specfun = 801
cumulant = .true.
bnd_cum = 1
fsthick = 1.36056981 ! eV
temps = 300 ! K
degaussw = 0.005 ! eV
dvscf_dir = '../phonons/save'
filukk = './diam.ukk'
filkf = 'meshes/path.dat'
nqf1 = 50
nqf2 = 50
nqf3 = 50
nk1 = 6
nk2 = 6
nk3 = 6
nq1 = 6
nq2 = 6
nq3 = 6
/

4
EPW/src/bcast_epw_input.f90
View File

@ -53,7 +53,7 @@
wannier_plot_supercell, wannier_plot_radius, &
fixsym, epw_no_t_rev, epw_tr, epw_nosym, epw_noinv, &
epw_crysym, bfieldx, bfieldy, bfieldz, tc_linear, &
tc_linear_solver
tc_linear_solver, mob_maxfreq, mob_nfreq
USE elph2, ONLY : elph
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
@ -200,6 +200,7 @@
CALL mp_bcast(bnd_cum , meta_ionode_id, world_comm)
CALL mp_bcast(mob_maxiter , meta_ionode_id, world_comm)
CALL mp_bcast(wannier_plot_supercell, meta_ionode_id, world_comm)
CALL mp_bcast(mob_nfreq , meta_ionode_id, world_comm)
!
! REAL*8
!
@ -245,6 +246,7 @@
CALL mp_bcast(bfieldx , meta_ionode_id, world_comm)
CALL mp_bcast(bfieldy , meta_ionode_id, world_comm)
CALL mp_bcast(bfieldz , meta_ionode_id, world_comm)
CALL mp_bcast(mob_maxfreq , meta_ionode_id, world_comm)
!
! characters
!

671
EPW/src/cum_mod.f90
View File

@ -11,7 +11,16 @@
MODULE cum_mod
!----------------------------------------------------------------------
!!
!! This module contains the various routines use for cumulant expansion
!! In this module, the hole / electron cumulant is calculated using the
!! convolution definition of the cumulant. This definition allows to
!! print out the quasi-particle and satellite contributions separately.
!!
!! The implementation of the Fourier-transformed time-ordered cumulant
!! has been retired.
!!
!! This module was initially implemented by C. Verdi and F. Caruso and
!! later updated by S. Ponce'. The latest version from July 2021 is by
!! Nikolaus Kandolf.
!!
IMPLICIT NONE
!
@ -21,21 +30,28 @@
SUBROUTINE spectral_cumulant()
!-----------------------------------------------------------------------
!!
!! Compute the electron spectral function including the electron-
!! phonon interaction using the (retarded) cumulant expansion method.
!! Use the frequency-dependent self-energy from the one-shot Migdal approximation
!! read from specfun_sup.elself.
!! Compute the second-order cumulant expansion spectral function using
!! the electron self-energy from the file specfun_sup.elself calculated
!! earlier. For occupied states e_{nk}<E_F, this routine calculates the
!! hole (or lesser) cumulant, for unoccupied states e_{nk}>E_F, the
!! electron (or greater) cumulant.
!!
!! See e.g. PRL 77, 2268 (1996) and PRB 90, 085112 (2014):
!! the cumulant spectral function can be calculated using a series of convolutions
!! in frequency space, or using FFT.
!! To converge the FFT with a small dt, we use zero padding of ImSigma outside the
!! calculated frequency range. If the convergence is not satisfactory, one can
!! set 'fact' to a larger value (6 is used as default, see the SUBROUTINE cumulant_time).
!! The implementation follows essentially the derivations given in
!! Ref.[1]. For a simpler implementation, the quasi-particle function
!! is written in terms of the full time-ordered self-energy as in
!! Ref.[2], and not in terms of the lesser and greater self-energy as
!! in Ref.[1]. Satellites are constructed from the imaginary part of
!! the lesser or greater self-energy; see also documentation on
!! https://epw-code.org
!!
!! [1] Gumhalter et al., PRB 94, 035103 (2016)
!! [2] Aryasetiawan, in Strong coulomb correlations in electronic
!! structure calculations edited by V.I. Anisimov (2000),
!! chapter 1
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP, i4b
USE constants_epw, ONLY : kelvin2eV, two, zero, ryd2ev, ryd2mev, ci
USE constants_epw, ONLY : kelvin2eV, zero, two, ryd2ev, ryd2mev, ci
USE constants, ONLY : pi
USE io_global, ONLY : stdout
USE io_var, ONLY : iospectral_sup, iospectral_cum
@ -73,14 +89,8 @@
!! Counter on temperatures
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: dw
!! Freq. increment
REAL(KIND = DP) :: e_thresh
!! Do the cumulant only for states with energy below e_thresh
REAL(KIND = DP) :: a1, a2, a3, a4
!! Auxiliary variables
REAL(KIND = DP) :: ekk
!! K-S energy
REAL(KIND = DP) :: zeta
!! Z factor
REAL(KIND = DP), ALLOCATABLE :: ww(:)
@ -91,51 +101,62 @@
!! Real self-energy for bnd_cum
REAL(KIND = DP), ALLOCATABLE :: sigmai(:, :)
!! Imaginary self-energy for bnd_cum
REAL(KIND = DP), ALLOCATABLE :: a_mig(:, :)
!! Migdal spectral function (same as in spectral_func.f90 )
REAL(KIND = DP), ALLOCATABLE :: a_cw(:)
!! Cumulant spectral function (convolutions)
REAL(KIND = DP), ALLOCATABLE :: a_ct(:)
!! Cumulant spectral function (FFT)
REAL(KIND = DP), ALLOCATABLE :: a_tmp(:)
!! Temporary spectral function
!
! e_thresh can be changed if one needs e.g. only states below the Fermi level
e_thresh = 10.d0 / ryd2ev ! referred to the Fermi level
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1)
REAL(KIND = DP), ALLOCATABLE :: a_ce(:)
!! Full cumulant spectral function
REAL(KIND = DP), ALLOCATABLE :: a_qp(:)
!! Quasi-particle part of the cumulant spectrum
REAL(KIND = DP), ALLOCATABLE :: a_s1(:)
!! First satellite of the cumulant spectrum
!
WRITE(stdout, '(/5x,a)') REPEAT('=',75)
WRITE(stdout, '(5x,a)') 'Performing the CUMULANT EXPANSION of the retarded Greens function to obtain'
WRITE(stdout, '(5x,a)') 'the electron spectral function.'
WRITE(stdout, '(5x,a)') 'Calculating the cumulant expansion spectral function from the self-energy in specfun_sup.elself'
WRITE(stdout, '(5x,a)') 'There is no spin degeneracy factor, as in spectral_func.f90'
WRITE(stdout, '(5x,a)') 'Warning: the routine is sequential but very fast.'
WRITE(stdout, '(5x,a/)') REPEAT('=',75)
!
!
! Read in self-energy from specfun_sup.elself
!
DO itemp = 1, nstemp
!
WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV
!
filespecsup = 'specfun_sup.elself.' // trim(adjustl(tp)) // 'K'
OPEN (UNIT = iospectral_sup, FILE = filespecsup, STATUS = 'old', IOSTAT = ios)
!
IF (ios /= 0) CALL errore('spectral_cumulant', 'opening file specfun_sup.elself', ABS(ios))
!
! determine number of k points, ibndmin, ibndmax
!
! Determine number of k points, ibndmin, ibndmax
!
DO im = 1, 6
!
READ(iospectral_sup, '(a)') line
!
ENDDO
!
DO im = 1, maxrecs
!
READ (iospectral_sup, *, IOSTAT = ios) i1, i2
!
IF (im == 1) ibndmin = i2
IF (ios /= 0) EXIT
IF (im == maxrecs) CALL errore('spectral_cumulant', 'increase maxrecs', 1)
!
ENDDO
!
REWIND(iospectral_sup)
!
nk = i1
ibndmax = i2
!
WRITE(stdout, '(5x,a/)') "Read self-energy from file specfun_sup.elself"
WRITE(stdout, '(5x,a,i4,a,i4,a,i4,a,f12.6/)') "Check: nk = ", nk, &
", ibndmin = ", ibndmin, ", ibndmax = ", ibndmax, " kbT (eV) = ", gtemp(itemp) * ryd2ev
!
!
! Allocate arrays for frequency, KS energy, self-energy and spectral functions
!
ALLOCATE(ww(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating ww', 1)
ALLOCATE(ek(nk), STAT = ierr)
@ -144,102 +165,142 @@
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmar', 1)
ALLOCATE(sigmai(nk, nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating sigmai', 1)
ALLOCATE(a_mig(nw_specfun, nk), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_mig', 1)
ALLOCATE(a_cw(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_cw', 1)
ALLOCATE(a_ct(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_ct', 1)
ALLOCATE(a_tmp(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_tmp', 1)
ALLOCATE(a_ce(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_ce', 1)
ALLOCATE(a_qp(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_qp', 1)
ALLOCATE(a_s1(nw_specfun), STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error allocating a_s1', 1)
!
!
! Read in KS energy, frequency grid and self-energy for designated band
!
! read and store Kohn-Sham energy, energy grid, real and im sigma for designated band
DO im = 1,6
!
READ(iospectral_sup, '(a)') line
!
ENDDO
!
!
!
! Read in KS energy and frequency in eV, self-energy in meV
!
DO ibnd = 1, ibndmax - ibndmin + 1
!
DO ik = 1, nk
!
DO iw = 1, nw_specfun
!
READ(iospectral_sup,*) i1, i2, a1, a2, a3, a4
!
IF (i2 == bnd_cum) THEN
! ek, w read in eV; Sigma read in meV
ek(ik) = a1 / ryd2ev
ww(iw) = a2 / ryd2ev
sigmar(ik, iw) = a3 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 )
sigmai(ik, iw) = a4 / ryd2mev ! / ( EXP(ww(iw)/eptemp )+1.d0 )
! spec func as in spectral_func.f90
a_mig(iw, ik) = ABS(sigmai(ik, iw)) / pi / ((ww(iw) - ek(ik) - sigmar(ik, iw))**two + (sigmai(ik, iw) )**two)
!
ek(ik) = DBLE(a1 / ryd2ev)
ww(iw) = DBLE(a2 / ryd2ev)
sigmar(ik, iw) = DBLE(a3 / ryd2mev) ! / ( EXP(ww(iw)/eptemp )+1.d0 )
sigmai(ik, iw) = DBLE(a4 / ryd2mev) ! / ( EXP(ww(iw)/eptemp )+1.d0 )
!
ENDIF
!
ENDDO
!
ENDDO
!
ENDDO
!
CLOSE(iospectral_sup)
!
! open file for cumulant spectral function
WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV
IF (bnd_cum < 10) THEN
WRITE(filespec, '(a,i1,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
ELSEIF (bnd_cum > 9 .AND. bnd_cum < 100) THEN
WRITE(filespec, '(a,i2,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
ELSE
WRITE(filespec, '(a,i3,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
!
! Calculated self-energy is retarded:
! Expecting Im Sigma_k(w) < 0 for all momenta and frequencies
!
IF( ANY(sigmai < zero ) ) THEN
!
WRITE(stdout,'(5x,a/)') 'WARNING: Some values of Im Sigma are larger than zero!'
WRITE(stdout,'(5x,a/)') ' Check validity of spectral function result!'
!
ENDIF
OPEN(UNIT = iospectral_cum, FILE = filespec)
!
WRITE(iospectral_cum, '(a)') '# k Energy [eV] A(k,w) [meV^-1] Z-factor '
WRITE(iospectral_cum, '(a)') '# with convolutions | using FFT '
WRITE(stdout, '(8x,a)') 'k Energy [eV] A(k,w) [meV^-1] Z-factor '
WRITE(stdout, '(8x,a)') ' with convolutions | using FFT '
!
! define index corresponding to omega=0 (Fermi level)
i0 = MINLOC(ABS(ww(:)), DIM = 1)
IF (ABS(ww(i0)) > dw) CALL errore('spectral_cumulant', 'w=0 needs to be included in [wmin:wmax]', 1)
a_cw = zero
a_ct = zero
a_tmp = zero
! Check if KS energy lies within the spectral function frequency window
!
DO ik = 1, nk
IF (ek(ik) < e_thresh) THEN
!
IF (ek(ik) < ww(1) .OR. ek(ik) > ww(nw_specfun-1)) THEN
!
ekk = ek(ik)
WRITE(stdout, '(5x,a,i5/)') 'WARNING: Frequency window too small at ik = ', ik
WRITE(stdout, '(5x,a,i5/)') ' Frequency window for spectral function must contain KS energy.'
WRITE(stdout, '(5x,a/)') ' Skipping this k point'
!
! Cumulant from convolutions in frequency space (ImSigma>0 in EPW)
CALL cumulant_conv(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_cw, zeta)
ENDIF
!
ENDDO
!
! Open file for cumulant spectral function output
!
WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV
!
IF (bnd_cum < 10) THEN
!
WRITE(filespec, '(a,i1,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
!
ELSEIF (bnd_cum > 9 .AND. bnd_cum < 100) THEN
!
WRITE(filespec, '(a,i2,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
!
ELSE
!
WRITE(filespec, '(a,i3,a,a,a)') 'specfun_cum', bnd_cum, '.elself.', trim(adjustl(tp)), 'K'
!
ENDIF
!
OPEN(UNIT = iospectral_cum, FILE = filespec)
!
! Prepare output to specfun_cum[bnd_cum].elself.[T]K and to standard output
!
WRITE(iospectral_cum, '(6x,a)') &
'#k ww-E_F[eV] A(k,w)[meV^{-1}] A^{QP}[meV^{-1}] A^{S1}[meV^{-1}] Z-factor'
WRITE(stdout,'(6x,a)') &
'#k ww-E_F[eV] A(k,w)[meV^{-1}] A^{QP}[meV^{-1}] A^{S1}[meV^{-1}] Z-factor'
!
!
DO ik = 1, nk ! main k loop
!
a_ce = zero
a_qp = zero
a_s1 = zero
!
! Call cumulant convolution subroutine for k points whose
! KS energy lies within the frequency range of the
! spectral function.
!
IF (ek(ik) > ww(1) .AND. ek(ik) < ww(nw_specfun-1)) THEN
CALL cumulant_conv(ek(ik), ww, sigmar(ik, :), sigmai(ik, :), a_ce, a_qp, a_s1, zeta)
ENDIF
!
! Write data to specfun_cum[bnd_cum].elself.[T]K and to standard output
!
DO iw = 1, nw_specfun
!
! Cumulant calculated in time domain + FFT (ImSigma>0 in EPW)
CALL cumulant_time(ekk, ww, sigmar(ik, :), -sigmai(ik, :), a_mig(:, ik), a_tmp)
!
DO iw = 1, nw_specfun
IF (iw == 1) THEN
!
! map the indices of the FFT frequency grid onto the original one
IF (iw >= i0) THEN
a_ct(iw) = a_tmp(iw - i0 + 1)
ELSE
a_ct(iw) = a_tmp(iw + nw_specfun - i0 + 1)
ENDIF
WRITE(iospectral_cum, '(x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,e16.7,3x,f8.4)') &
ik, ww(iw) * ryd2ev, a_ce(iw) / ryd2mev, a_qp(iw) / ryd2mev, a_s1(iw) / ryd2mev, zeta
WRITE(stdout,'(x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,e16.7,3x,f8.4)') &
ik, ww(iw) * ryd2ev, a_ce(iw) / ryd2mev, a_qp(iw) / ryd2mev, a_s1(iw) / ryd2mev, zeta
!
! write cumulant spectral function on file (in meV^-1, as in spectral_func.f90)
! 3rd column: A_cum using convolutions; 4th column: A_cum using FFT
IF (iw == 1) THEN
WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') &
ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta
WRITE(stdout,'(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,f8.4)') &
ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev, zeta
ELSE
WRITE(iospectral_cum, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') &
ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 )
WRITE(stdout, '(2x,i7,2x,f10.5,3x,e16.7,3x,e16.7)') &
ik, ww(iw) * ryd2ev, a_cw(iw) / ryd2mev, a_ct(iw) / ryd2mev !/ ( EXP(ww(iw)/eptemp )+1.d0 )
! uncomment to multiply by Fermi occupation factor
ENDIF
ELSE
!
ENDDO
WRITE(iospectral_cum, '(x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,e16.7)') &
ik, ww(iw) * ryd2ev, a_ce(iw) / ryd2mev, a_qp(iw) / ryd2mev, a_s1(iw) / ryd2mev
WRITE(stdout, '(x,i7,2x,f10.5,3x,e16.7,3x,e16.7,3x,e16.7)') &
ik, ww(iw) * ryd2ev, a_ce(iw) / ryd2mev, a_qp(iw) / ryd2mev, a_s1(iw) / ryd2mev
!
ENDIF
!
WRITE(iospectral_cum, '(a)') ' '
WRITE(stdout, '(a)') ' '
!
ENDIF ! only states below energy e_thresh
ENDDO
!
WRITE(iospectral_cum, '(a)') ' '
WRITE(stdout, '(a)') ' '
!
ENDDO ! main loop k
!
@ -255,37 +316,39 @@
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmar', 1)
DEALLOCATE(sigmai, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating sigmai', 1)
DEALLOCATE(a_mig, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_mig', 1)
DEALLOCATE(a_cw, STAT = ierr)
DEALLOCATE(a_ce, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_cw', 1)
DEALLOCATE(a_ct, STAT = ierr)
DEALLOCATE(a_qp, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_ct', 1)
DEALLOCATE(a_tmp, STAT = ierr)
DEALLOCATE(a_s1, STAT = ierr)
IF (ierr /= 0) CALL errore('spectral_cumulant', 'Error deallocating a_tmp', 1)
ENDDO !itemp
!
ENDDO ! Loop over temperatures
!
! Deallocate temperature if supercond = .FALSE.
!
! Deallocate temperature when no supercond
IF (.NOT. eliashberg) THEN
!
DEALLOCATE(gtemp, STAT = ierr)
IF (ierr /= 0) CALL errore('cum_mod', 'Error deallocating gtemp', 1)
!
ENDIF
!
!-----------------------------------------------------------------------
END SUBROUTINE spectral_cumulant
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE cumulant_conv(ek, ww, sigmar, sigmai, a_mig, a_cum, zeta)
SUBROUTINE cumulant_conv(ek, ww, sigmar, sigmai, a_ce, a_qp, a_s1, zeta)
!-----------------------------------------------------------------------
!!
!! Convolution for the cumulant expansion.
!! Calculated cumulant spectral function from convolution
!!
!
USE kinds, ONLY : DP
USE constants_epw, ONLY : two, zero, ci, ryd2ev
USE constants_epw, ONLY : zero, one, two, ci, ryd2ev, ryd2mev
USE constants, ONLY : pi
USE epwcom, ONLY : degaussw, wmin_specfun, wmax_specfun, nw_specfun
USE io_global, ONLY : stdout
!
IMPLICIT NONE
!
@ -295,12 +358,14 @@
!! Frequency variable
REAL(KIND = DP), INTENT(in) :: sigmar(nw_specfun)
!! Real self-energy
REAL(KIND = DP), INTENT(in) :: sigmai(nw_specfun)
REAL(KIND = DP), INTENT(INOUT) :: sigmai(nw_specfun)
!! Imaginary self-energy
REAL(KIND = DP), INTENT(in) :: a_mig(nw_specfun)
!! Migdal spectral function
REAL(KIND = DP), INTENT(out) :: a_cum(nw_specfun)
REAL(KIND = DP), INTENT(out) :: a_ce(nw_specfun)
!! Cumulant spectral function
REAL(KIND = DP), INTENT(out) :: a_qp(nw_specfun)
!! Quasi-particle part of spectral function
REAL(KIND = DP), INTENT(out) :: a_s1(nw_specfun)
!! First satellite of spectral function
REAL(KIND = DP), INTENT(out) :: zeta
!! Z factor
!
@ -309,255 +374,201 @@
!! Frequency counter
INTEGER :: iw2
!! Second freq. counter for the convolution
INTEGER :: indw
INTEGER :: icoul
!! Auxillary indices for convolution
INTEGER :: ind1
!! Auxillary indices for convolution
INTEGER :: isat
!! Index to number the satellite
INTEGER :: iqp
!! Energy index of renormalized quasiparticle
INTEGER :: iks
INTEGER :: iek
!! Energy index for the KS quasiparticle
INTEGER :: i0
!! Energy index of Fermi level
REAL(KIND = DP) :: dw
!! Freq. increment
REAL(KIND = DP) :: eqp
!! Renpormalized quasiparticle energy
REAL(KIND = DP) :: si_ks
REAL(KIND = DP) :: rek
!! ReSigma at KS energy
REAL(KIND = DP) :: imk
!! ImSigma at KS energy
REAL(KIND = DP) :: si_qp
!! ImSigma at renormalized energy
REAL(KIND = DP) :: diS
!! Derivative of ImSigma at KS energy
REAL(KIND = DP) :: drS
!! Derivative of ReSigma at qp energy
REAL(KIND = DP) :: a_qp(nw_specfun)
!! Quasiparticle contribution to the spectral function
REAL(KIND = DP) :: Gammak
!! +/- ImSigma at KS energy
REAL(KIND = DP) :: a_k
!! alpha_k factor for A^{QP}
REAL(KIND = DP) :: drek
!! Derivative of ReSigma at ek
REAL(KIND = DP) :: dimk
!! Derivative of ImSigma at ek
REAL(KIND = DP) :: d2imk
!! Second-order derivative of ImSigma at ek
REAL(KIND = DP) :: ww_as(nw_specfun)
!! Auxiliary frequency grid to construct satellite function
REAL(KIND = DP) :: a_s(nw_specfun)
!! Temporary quantity needed to compute the satelite
REAL(KIND = DP) :: conv(nw_specfun)
!! Temporary quantity to compute the convolution
REAL(KIND = DP) :: a_s1(nw_specfun), a_s2(nw_specfun), a_s3(nw_specfun)
!! satellites
!
! Initialization
diS = zero
drS = zero
REAL(KIND = DP) :: a_s2(nw_specfun), a_s3(nw_specfun)
!! Second and third satellites
!
! Extract indices for E_F (w=0) and the IP energy
!
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1)
i0 = MINLOC(ABS(ww(:)), DIM = 1)
! define index and energy of renormalized qp (note: qp energy needs to be the absolute max of a_mig)
iqp = MAXLOC(a_mig(:), DIM = 1)
eqp = ww(iqp)
iek = MINLOC(ABS(ww(:) - ek), DIM = 1)
!
! define index corresponding to unrenormalized (Kohn-Sham) energy
iks = MINLOC(ABS(ww(:) - ek), DIM = 1)
! Construct auxiliary frequency grid for construction of satellite
! function, in which ww[iek] = ek is contained exactly.
!
si_qp = ABS(sigmai(iqp))
si_ks = ABS(sigmai(iks))
! finite difference derivatives of Im and Re Sigma
IF (iks - 1 > 0) diS = (ABS(sigmai(iks + 1)) - ABS(sigmai(iks - 1))) / (2.d0 * dw)
IF (iqp - 1 > 0) drS = (sigmar(iqp + 1) - sigmar(iqp - 1)) / (2.d0 * dw)
zeta = EXP(drS)
IF ( ww(iek) > ek ) THEN
!
ww_as = ww - ww(iek) + ek
!
ELSE
!
ww_as = ww + ww(iek) - ek
!
ENDIF
!
!
! Initialise on-the-mass-shell self-energy: incoming self-energy
! is retarded ( Im(Sigma)>0 everywhere), need lesser and greater
! self-energy ( Im(Sigma)>0 for k>k_F, Im(Sigma)<0 for k<k_F)
!
rek = sigmar(iek)
!
IF (ek < 0) THEN ! Lesser self-energy
!
sigmai(i0:) = zero
imk = sigmai(iek) ! This is Im Sigma^<(ek)
Gammak = MAX(sigmai(iek),degaussw) ! This is \Gamma in Eq.(51) in Ref.[1]
!
ELSE ! Greater self-energy
!
sigmai(i0:) = -sigmai(i0:)
sigmai(1:i0) = zero
imk = sigmai(iek) ! This is Im Sigma^>(ek)
Gammak = MAX(-sigmai(iek),degaussw) ! This is \Gamma in Eq.(51) in Ref.[1]
!
ENDIF
!
!
! Calculate finite difference derivatives of Re and Im Sigma
!
dw = ABS(ww_as(1)-ww_as(2))
!
drek = (sigmar(iek+1) - sigmar(iek-1)) / (two * dw)
dimk = (sigmai(iek+1) - sigmai(iek-1)) / (two * dw)
d2imk = (sigmai(iek+1) - two * sigmai(iek) + sigmai(iek-1)) / &
(dw**two)
!
!
! Calculate Z_k and alpha_k factors
!
zeta = DEXP(drek)
a_k = -dimk
!
!
! Calculate A^{QP}
!
a_qp = zero
!
! calculate Aqp and As1
DO iw = 1, nw_specfun
!
a_qp(iw) = zeta * si_qp / pi / ((ww(iw) - eqp)**two + (si_qp)**two)
conv(iw) = a_qp(iw)
a_qp(iw) = zeta / pi * &
(Gammak * COS(a_k) - (ww(iw) - ek - rek) * SIN(a_k)) / &
((ww(iw) - ek - rek)**two + Gammak**two)
!
ind1 = iks + iw - i0
IF (ind1 > 0 .AND. ind1 < nw_specfun) THEN !RC
a_s(iw) = (ABS(sigmai(ind1)) - si_ks - (ww(iw)) * diS) * REAL(1.d0 / (ww(iw) - ci * degaussw)**2.d0) / pi
ENDDO
!
!
! Calculate A^{S}
!
a_s = zero
!
DO iw = 1, nw_specfun
!
IF (iw == iek) THEN
!
a_s(iw) = one / two / pi * d2imk
!
ELSE
a_s(iw) = (-si_ks - (ww(iw)) * diS) * REAL(1.d0 / (ww(iw) - ci * degaussw)**2.d0) / pi
!
IF (ek < 0) THEN ! Lesser self-energy
!
! This is Eq.(67) in [1] for k<k_F
!
IF (sigmai(iw-iek+i0) > 0.d0) THEN
a_s(iw) = (sigmai(iw) - (imk + (ww_as(iw) - ek) * dimk ) ) / &
(ww_as(iw) - ek)**two / pi
ENDIF
!
ELSE ! Greater self-energy
!
! This is Eq.(67) in [1] for k>k_F
!
IF (sigmai(iw-i0+iek) < 0.d0) THEN
a_s(iw) = (-sigmai(iw) - (-imk - (ww_as(iw) - ek) * dimk ) ) / &
(ww_as(iw) - ek)**two / pi
ENDIF
!
ENDIF
!
ENDIF
!
ENDDO
!
DO isat = 1, 3
!
a_cum = zero
!
DO iw = 1, nw_specfun
!
DO iw2 = 1, nw_specfun
!
indw = i0 + iw - iw2
IF (indw <= nw_specfun .AND. indw > 0) THEN
a_cum(iw) = a_cum(iw) + ABS(a_s(iw2) * conv(indw)) * dw
ENDIF
!
ENDDO
!
ENDDO
!
IF (isat == 1) a_s1 = a_cum
IF (isat == 2) a_s2 = a_cum / 2.d0
IF (isat == 3) a_s3 = a_cum / 6.d0
conv = a_cum
!
ENDDO ! isat
! Calculate first satellite
!
a_s1 = zero
!
DO iw = 1, nw_specfun
!
a_cum(iw) = a_qp(iw) + a_s1(iw) + a_s2(iw) + a_s3(iw)
DO iw2 = 1, nw_specfun
!
icoul = i0 + iw - iw2
!
a_s1(iw) = a_s1(iw) + a_qp(iw2) * a_s(icoul) * dw
!
ENDDO
!
ENDDO
!
! Calculate second satellite
!
a_s2 = zero
!
DO iw = 1, nw_specfun
!
DO iw2 = 1, nw_specfun
!
icoul = i0 + iw - iw2
!
a_s2(iw) = a_s2(iw) + a_s1(iw2) * a_s(icoul) * dw
!
ENDDO
!
ENDDO
!
! Calculate third satellite
!
a_s3 = zero
!
DO iw = 1, nw_specfun
!
DO iw2 = 1, nw_specfun
!
icoul = i0 + iw - iw2
!
a_s3(iw) = a_s3(iw) + a_s2(iw2) * a_s(icoul) * dw
!
ENDDO
!
ENDDO
!
a_ce = zero
!
a_ce = a_qp + a_s1 + a_s2 / two + a_s3 / 6.d0
!
!-----------------------------------------------------------------------
END SUBROUTINE cumulant_conv
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE cumulant_time(ek, ww, sigmar, sigmai, a_mig, a_cum)
!-----------------------------------------------------------------------
!!
!! Cumulant calculated in the time domain
!!
!
USE kinds, ONLY : DP
USE constants_epw, ONLY : two, zero, czero, ci, ryd2ev
USE constants, ONLY : pi
USE fft_scalar, ONLY : cfft3d
USE epwcom, ONLY : degaussw, wmin_specfun, wmax_specfun, nw_specfun
!
IMPLICIT NONE
!
REAL(KIND = DP), INTENT(in) :: ek
!! K-S energy
REAL(KIND = DP), INTENT(in) :: ww(nw_specfun)
!! Frequency variable
REAL(KIND = DP), INTENT(in) :: sigmar(nw_specfun)
!! Real self-energy
REAL(KIND = DP), INTENT(in) :: sigmai(nw_specfun)
!! Imaginary self-energy
REAL(KIND = DP), INTENT(in) :: a_mig(nw_specfun)
!! Migdal spectral function
REAL(KIND = DP), INTENT(out) :: a_cum(nw_specfun)
!! Cumulant spectral function
!
! local variables
!
INTEGER :: iw
!! Frequency counters
INTEGER :: it
!! Time counters
INTEGER :: iqp
!! Energy index of renormalized quasiparticle
INTEGER :: iks
!! Energy index for the KS quasiparticle
INTEGER :: i0
!! Energy index of Fermi level
INTEGER :: ind1, ind2
!! aux indices
INTEGER :: nw_new
!! number of points used for FFT
!! Zero padding is used for ImSigma outside the energy range computed.
INTEGER :: ierr
!! Error status
REAL(KIND = DP) :: fact
!! factor used to increase the number of points for the FFT;
!! fact=4 gives good convergence but it can be increased if needed.
REAL(KIND = DP) :: dw
!! Freq. increment
REAL(KIND = DP) :: eqp
!! Renormalized quasiparticle energy
REAL(KIND = DP) :: drS, zeta
!! Derivative of ReSigma at qp energy, Z factor
REAL(KIND = DP) :: smeart
!! small broadening
REAL(KIND = DP) :: dt
!! time increment
REAL(KIND = DP) :: tmin, tmax
!! min and max of time interval
REAL(KIND = DP) :: tt
!! time variable
COMPLEX(KIND = DP), ALLOCATABLE :: cumS(:)
!! satellite part of the cumulant function
COMPLEX(KIND = DP), ALLOCATABLE :: cum(:)
!! complete cumulant function
COMPLEX(KIND = DP) :: qpfac
!! quasiparticle factor
!
fact = 6.d0
! fact=1 corresponds to FFT with the original no. of points nw_specfun
dw = (wmax_specfun - wmin_specfun) / DBLE(nw_specfun - 1)
smeart = -degaussw ! correct sign for RC
!
tmin = 0.d0
tmax = two * pi / dw
dt = two * pi / ((wmax_specfun - wmin_specfun) * fact)
nw_new = INT(fact * (nw_specfun - 1) + 1) ! to be consistent with dt above
!
ALLOCATE(cumS(nw_new), STAT = ierr)
IF (ierr /= 0) CALL errore('cumulant_time', 'Error allocating cumS', 1)
ALLOCATE(cum(nw_new), STAT = ierr)
IF (ierr /= 0) CALL errore('cumulant_time', 'Error allocating cum', 1)
!
i0 = MINLOC(ABS(ww(:)), DIM = 1)
!
! define index and energy of renormalized qp (note: qp energy needs to be the absolute max of a_mig)
iqp = MAXLOC(a_mig(:), DIM = 1)
eqp = ww(iqp)
!
! define index corresponding to unrenormalized (Kohn-Sham) energy
iks = MINLOC(ABS(ww(:) - ek), DIM = 1)
!
drS = (sigmar(iqp + 1) - sigmar(iqp - 1)) / (two * dw)
qpfac = EXP(-ci * (ek + sigmar(iqp)) + smeart * 0.5d0)
zeta = EXP(drS)
!
cumS = czero
DO iw = 1, nw_new
!
! the w shift is needed because FFT uses positive w \in [0:Omega], Omega=wmax-wmin
IF (iw <= (nw_specfun - i0 + 1)) THEN
ind1 = iw + i0 - 1
cumS(iw) = dw * ABS(sigmai(ind1)) / pi * REAL(1.d0 / (ek - ww(ind1) - ci * smeart)**two)
ELSE IF (iw > (nw_new - i0 + 1)) THEN
ind2 = iw - nw_new + i0 - 1
cumS(iw) = dw * ABS(sigmai(ind2)) / pi * REAL(1.d0 / (ek - ww(ind2) - ci * smeart)**two)
ENDIF
!
ENDDO
!
CALL cfft3d(cumS(:), nw_new, 1, 1, nw_new, 1, 1, 1, -1)
!this is needed because cfft3d(...,-1) carries a renomalization factor 1/nw_new
cumS = cumS * nw_new
!
DO it = 1, nw_new
!
tt = tmin + DBLE(it - 1) * dt
cumS(it) = cumS(it) * EXP(ci * ek * tt)
cum(it) = zeta * (qpfac**tt) * EXP(cumS(it))
!
ENDDO
!
CALL cfft3d(cum(:), nw_new, 1, 1, nw_new, 1, 1, 1, 1)
cum = cum * dt / pi
!
! extract the spectral function a_cum on the original w FFT grid (nw_specfun points)
DO iw = 1, nw_specfun
IF (iw <= (nw_specfun - i0 + 1)) THEN
a_cum(iw) = REAL(REAL(cum(iw)))
ELSE
ind1 = iw + nw_new - nw_specfun
a_cum(iw) = REAL(REAL(cum(ind1)))
ENDIF
ENDDO
!
DEALLOCATE(cumS, STAT = ierr)
IF (ierr /= 0) CALL errore('cumulant_time', 'Error deallocating cumS', 1)
DEALLOCATE(cum, STAT = ierr)
IF (ierr /= 0) CALL errore('cumulant_time', 'Error deallocating cum', 1)
!
!-----------------------------------------------------------------------
END SUBROUTINE cumulant_time
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
END MODULE cum_mod
!-----------------------------------------------------------------------

2
EPW/src/elph2.f90
View File

@ -124,6 +124,8 @@
inv_tau_allcb(:, :, :), &! Conduction band scattering rate
inv_tau_all_mode(:, :, :, :), &! mode resolved scattering rate
inv_tau_allcb_mode(:, :, :, :), &! mode resolved conduction band scattering rate
inv_tau_all_freq(:, :, :), &! Scattering rate spectral decomposition (for one temperature)
inv_tau_allcb_freq(:, :, :), &! Scattering rate conduction band spectral decomposition (for one temperature)
zi_allvb(:, :, :), &! Z-factor in scattering rate
zi_allcb(:, :, :), &! Second Z-factor in scattering rate (for both VB and CB calculations)
ifc(:,:,:,:,:,:,:), &! Interatomic force constant in real space

26
EPW/src/ephwann_shuffle.f90
View File

@ -35,8 +35,8 @@
specfun_pl, lindabs, use_ws, epbread, fermi_plot, &
epmatkqread, selecqread, restart_step, nsmear, &
nqc1, nqc2, nqc3, nkc1, nkc2, nkc3, assume_metal, &
cumulant, eliashberg, nomega, &
omegamin, omegamax, omegastep, neta
cumulant, eliashberg, nomega, mob_maxfreq, neta, &
omegamin, omegamax, omegastep, mob_nfreq
USE control_flags, ONLY : iverbosity
USE noncollin_module, ONLY : noncolin
USE constants_epw, ONLY : ryd2ev, ryd2mev, one, two, zero, czero, eps40, &
@ -59,7 +59,8 @@
nktotf, gtemp, xkq, dos, nbndskip, nbndep, &
inv_tau_all_mode, inv_tau_allcb_mode, qrpl, Qmat, &
ef0_fca, epsilon2_abs, epsilon2_abs_lorenz, &
epsilon2_abs_all, epsilon2_abs_lorenz_all
epsilon2_abs_all, epsilon2_abs_lorenz_all, &
inv_tau_all_freq, inv_tau_allcb_freq
USE wan2bloch, ONLY : dmewan2bloch, hamwan2bloch, dynwan2bloch, &
ephwan2blochp, ephwan2bloch, vmewan2bloch, &
dynifc2blochf, vmewan2blochp
@ -136,8 +137,6 @@
!! If the file exist
LOGICAL :: first_cycle
!! Check wheter this is the first cycle after a restart.
LOGICAL :: first_time
!! Check wheter this is the first timeafter a restart.
LOGICAL :: homogeneous
!! Check if the k and q grids are homogenous and commensurate.
INTEGER :: ios
@ -855,8 +854,14 @@
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating inv_tau_all_mode', 1)
ALLOCATE(inv_tau_allcb_mode(nmodes, nbndfst, nktotf, nstemp), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating inv_tau_allcb_mode', 1)
ALLOCATE(inv_tau_all_freq(mob_nfreq, nbndfst, nktotf), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating inv_tau_all_freq', 1)
ALLOCATE(inv_tau_allcb_freq(mob_nfreq, nbndfst, nktotf), STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error allocating inv_tau_allcb_freq', 1)
inv_tau_all_mode(:, :, :, :) = zero
inv_tau_allcb_mode(:, :, :, :) = zero
inv_tau_all_freq(:, :, :) = zero
inv_tau_allcb_freq(:, :, :) = zero
ENDIF
! We save matrix elements that are smaller than machine precision (1d-16).
! The sum of all the elements must be smaller than that
@ -895,6 +900,10 @@
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating inv_tau_all_mode', 1)
DEALLOCATE(inv_tau_allcb_mode, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating inv_tau_allcb_mode', 1)
DEALLOCATE(inv_tau_all_freq, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating inv_tau_all_freq', 1)
DEALLOCATE(inv_tau_allcb_freq, STAT = ierr)
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'Error deallocating inv_tau_allcb_freq', 1)
ENDIF
!
ELSE ! (iterative_bte .AND. epmatkqread)
@ -972,7 +981,6 @@
! -----------------------------------------------------------------------
iq_restart = 1
first_cycle = .FALSE.
first_time = .TRUE.
!
! Fine mesh set of g-matrices. It is large for memory storage
ALLOCATE(epf17(nbndfst, nbndfst, nmodes, nkf), STAT = ierr)
@ -1179,9 +1187,11 @@
#endif
IF (ierr /= 0) CALL errore('ephwann_shuffle', 'error in MPI_BCAST', 1)
!
IF(ephwrite .AND. iq_restart > 1) THEN
IF(iq_restart > 1) THEN
first_cycle = .TRUE.
CALL check_restart_ephwrite
IF (ephwrite) THEN
CALL check_restart_ephwrite
ENDIF
ENDIF
!
! Now, the iq_restart point has been done, so we need to do the next

4
EPW/src/epw.f90
View File

@ -10,7 +10,7 @@
PROGRAM epw
!-----------------------------------------------------------------------
!! author: Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
!! version: v5.4
!! version: v5.4.1
!! license: GNU
!! summary: EPW main driver
!!
@ -38,7 +38,7 @@
CHARACTER(LEN = 12) :: code = 'EPW'
!! Name of the program
!
version_number = '5.4'
version_number = '5.4.1'
!
CALL init_clocks(.TRUE.)
!

10
EPW/src/epw_readin.f90
View File

@ -60,7 +60,8 @@
wannier_plot_supercell, wannier_plot_scale, reduce_unk, &
wannier_plot_radius, fermi_plot, fixsym, epw_no_t_rev, &
epw_tr, epw_nosym, epw_noinv, epw_crysym, &
bfieldx, bfieldy, bfieldz, tc_linear, tc_linear_solver
bfieldx, bfieldy, bfieldz, tc_linear, tc_linear_solver, &
mob_maxfreq, mob_nfreq
USE klist_epw, ONLY : xk_all, xk_loc, xk_cryst, isk_all, isk_loc, et_all, et_loc
USE elph2, ONLY : elph, num_wannier_plot, wanplotlist, gtemp
USE constants_epw, ONLY : ryd2mev, ryd2ev, ev2cmm1, kelvin2eV, zero, eps20, ang2m
@ -161,7 +162,7 @@
scdm_sigma, assume_metal, wannier_plot, wannier_plot_list, reduce_unk, &
wannier_plot_supercell, wannier_plot_scale, wannier_plot_radius, &
fixsym, epw_no_t_rev, epw_tr, epw_nosym, epw_noinv, epw_crysym, &
tc_linear, tc_linear_solver, &
tc_linear, tc_linear_solver, mob_maxfreq, mob_nfreq, &
!---------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
! Shell implementation for future use.
@ -578,6 +579,8 @@
bfieldx = 0.d0 ! Tesla
bfieldy = 0.d0 ! Tesla
bfieldz = 0.d0 ! Tesla
mob_maxfreq = 100 ! Maximum frequency for spectral decomposition in meV
mob_nfreq = 100 ! Number of frequency for the spectral decomposition
!
! --------------------------------------------------------------------------------
! Added for polaron calculations. Originally by Danny Sio, modified by Chao Lian.
@ -888,6 +891,9 @@
degaussq = degaussq / ryd2mev
delta_qsmear = delta_qsmear / ryd2mev
!
! Max frequency for the spectral decomposition of mobility from meV to Ry.
mob_maxfreq = mob_maxfreq / ryd2mev
!
! fermi_energy read from the input file from eV to Ryd
IF (efermi_read) THEN
fermi_energy = fermi_energy / ryd2ev

4
EPW/src/epwcom.f90
View File

@ -331,6 +331,8 @@
!! max frequency in electron spectral function due to e-p `interaction
!
! Conductivity
INTEGER :: mob_nfreq
!! Number of frequency for the spectral decomposition of mobility
REAL(KIND = DP) :: scissor
!! Value of the scissor shift in eV.
REAL(KIND = DP) :: ncarrier
@ -343,6 +345,8 @@
!! Magnetic field along the y-direction
REAL(KIND = DP) :: bfieldz
!! Magnetic field along the z-direction
REAL(KIND = DP) :: mob_maxfreq
!! Maximum frequency for the spectral decomposition of mobility. Typically that freq. is the highest phonon freq.
!
! Plasmon
REAL(KIND = DP) :: nel

26
EPW/src/io_eliashberg.f90
View File

@ -114,7 +114,7 @@
! anisotropic case
IF (temp < 10.d0) THEN
WRITE(name1, 101) TRIM(prefix), '.imag_aniso_00', temp
ELSEIF (temp >= 10.d0) THEN
ELSEIF (temp >= 10.d0 .AND. temp < 100.d0) THEN
WRITE(name1, 102) TRIM(prefix), '.imag_aniso_0', temp
ELSEIF (temp >= 100.d0) THEN
WRITE(name1, 103) TRIM(prefix), '.imag_aniso_', temp
@ -2202,6 +2202,27 @@
CLOSE(iufilgapFS)
ENDDO
!
! HP: Write in .frmsf format compatible with fermisurfer program
IF (temp < 10.d0) THEN
WRITE(name1, 113) TRIM(prefix), '.', cname, '_aniso_gap0_00', temp, '.frmsf'
ELSEIF (temp < 100.d0 .AND. temp > 9.9999d0) THEN
WRITE(name1, 114) TRIM(prefix), '.', cname, '_aniso_gap0_0', temp, '.frmsf'
ELSEIF (temp < 1000.d0 .AND. temp > 99.9999d0) THEN
WRITE(name1, 115) TRIM(prefix), '.', cname, '_aniso_gap0_', temp, '.frmsf'
ENDIF
OPEN(UNIT = iufilgapFS, FILE = name1, STATUS = 'unknown', FORM = 'formatted', IOSTAT = ios)
IF (ios /= 0) CALL errore('gap_FS', 'error opening file ' // name1, iufilgapFS)
!
WRITE(iufilgapFS, '(3i5)') nkf1, nkf2, nkf3
WRITE(iufilgapFS, '(i5)') 1
WRITE(iufilgapFS, '(i5)') nbndfs
WRITE(iufilgapFS, '(3f12.6)') (bg(i, 1), i = 1, 3)
WRITE(iufilgapFS, '(3f12.6)') (bg(i, 2), i = 1, 3)
WRITE(iufilgapFS, '(3f12.6)') (bg(i, 3), i = 1, 3)
WRITE(iufilgapFS, '(6f12.6)') ((ekfs(ibnd, ixkff(ik)) - ef0, ik = 1, nkf1 * nkf2 * nkf3), ibnd = 1, nbndfs)
WRITE(iufilgapFS, '(6f12.6)') ((agap_tmp(ibnd, ixkff(ik)) * 1000.d0, ik = 1, nkf1 * nkf2 * nkf3), ibnd = 1, nbndfs)
CLOSE(iufilgapFS)
!
ENDIF
!
! SP & RM : Write on file the superconducting gap close to the Fermi surface
@ -2257,6 +2278,9 @@
110 FORMAT(a, a1, a4, a16, f4.2)
111 FORMAT(a, a1, a4, a15, f5.2)
112 FORMAT(a, a1, a4, a14, f6.2)
113 FORMAT(a, a1, a4, a14, f4.2, a6)
114 FORMAT(a, a1, a4, a13, f5.2, a6)
115 FORMAT(a, a1, a4, a12, f6.2, a6)
!
RETURN
!

79
EPW/src/io_transport.f90
View File

@ -30,13 +30,15 @@
USE io_global, ONLY : stdout
USE modes, ONLY : nmodes
USE epwcom, ONLY : fsthick, eps_acustic, degaussw, nstemp, vme, ncarrier, &
assume_metal, etf_mem, nqf1, nqf2, nqf3, system_2d
assume_metal, etf_mem, nqf1, nqf2, nqf3, system_2d, &
mob_maxfreq, mob_nfreq
USE pwcom, ONLY : ef
USE elph2, ONLY : ibndmin, etf, nkf, dmef, vmef, wf, wqf, &
epf17, inv_tau_all, inv_tau_allcb, adapt_smearing, &
wkf, dmef, vmef, eta, gtemp, lower_bnd, dos, &
nbndfst, nktotf, vkk_all, carrier_density, &
inv_tau_all_mode, inv_tau_allcb_mode
inv_tau_all_mode, inv_tau_allcb_mode, &
inv_tau_all_freq, inv_tau_allcb_freq
USE constants_epw, ONLY : zero, one, two, pi, ryd2mev, kelvin2eV, ryd2ev, eps4, eps8, &
eps6, eps20, bohr2ang, ang2cm, hbarJ, eps160
USE constants, ONLY : electronvolt_si
@ -117,6 +119,8 @@
!! Pool index
INTEGER :: i,j
!! Cartesian index
INTEGER :: ifreq
!! Index on frequency
INTEGER :: ind(npool)
!! Nb of Matrix elements that are non-zero
INTEGER :: indcb(npool)
@ -204,6 +208,10 @@
!! Auxiliary variables
REAL(KIND = DP) :: inv_tau_allcb_MPI(nbndfst, nktotf, nstemp)
!! Auxiliary variables
REAL(KIND = DP) :: step
!! Energy step in Ry for the spectral decomposition
REAL(KIND = DP) :: mobility_freq(3, 3, mob_nfreq)
!! Spectral decomposition of mobility
REAL(KIND = DP), EXTERNAL :: efermig
!! Function that returns the Fermi energy
REAL(KIND = DP), EXTERNAL :: wgauss
@ -222,7 +230,12 @@
wqf_loc = wqf(iq)
ENDIF
!
IF (iqq == 1) THEN
IF (iverbosity == 3) THEN
! Energy steps for spectral decomposition
step = mob_maxfreq / mob_nfreq
ENDIF
!
IF (iqq == 1 .OR. first_cycle) THEN
!
WRITE(stdout, '(/5x,a)') REPEAT('=',67)
WRITE(stdout, '(5x,"Scattering rate for IBTE")')
@ -282,7 +295,7 @@
WRITE(iufilmu_q, '(a)') '# \mu(alpha,beta) = 1.0 / (sum_{\nu q} T_{\nu q}(alpha,beta))'
ENDIF
ENDIF ! iverbosity
ENDIF ! iqq == 1
ENDIF ! iqq == 1 .OR. first_cycle
IF (iverbosity == 3) mobilityq(:, :, :, :) = zero
!
! In the case of a restart do not add the first step
@ -495,6 +508,16 @@
inv_tau_all_mode(imode, ibnd, ik + lower_bnd - 1, itemp) = &
inv_tau_all_mode(imode, ibnd, ik + lower_bnd - 1, itemp) + &
two * pi * wqf_loc * g2 * ((fmkq + wgq(imode)) * w0g1 + (one - fmkq + wgq(imode)) * w0g2)
!
! Spectral decomposition histogram
ifreq = NINT(wq(imode) / step)
IF (step > eps20) THEN
ifreq = NINT(wq(imode) / step) + 1
ELSE
ifreq = 1
ENDIF
inv_tau_all_freq(ifreq, ibnd, ik + lower_bnd - 1) = inv_tau_all_freq(ifreq, ibnd, ik + lower_bnd - 1) &
+ two * pi * wqf_loc * g2 * ((fmkq + wgq(imode)) * w0g1 + (one - fmkq + wgq(imode)) * w0g2)
ENDDO !imode
ENDDO ! jbnd
DO j = 1, 3
@ -509,7 +532,7 @@
ENDDO ! i
ENDDO ! j
ENDDO ! ibnd
ENDIF ! iverbsoity
ENDIF ! iverbosity
ENDIF ! ctype
!
! In this case we are also computing the scattering rate for another Fermi level position
@ -573,6 +596,16 @@
inv_tau_allcb_mode(imode, ibnd, ik + lower_bnd - 1, itemp) = &
inv_tau_allcb_mode(imode, ibnd, ik + lower_bnd - 1, itemp) + &
two * pi * wqf_loc * g2 * ((fmkq + wgq(imode)) * w0g1 + (one - fmkq + wgq(imode)) * w0g2)
!
! Spectral decomposition histogram
ifreq = NINT(wq(imode) / step)
IF (step > eps20) THEN
ifreq = NINT(wq(imode) / step) + 1
ELSE
ifreq = 1
ENDIF
inv_tau_allcb_freq(ifreq, ibnd, ik + lower_bnd - 1) = inv_tau_allcb_freq(ifreq, ibnd, ik + lower_bnd - 1) &
+ two * pi * wqf_loc * g2 * ((fmkq + wgq(imode)) * w0g1 + (one - fmkq + wgq(imode)) * w0g2)
ENDDO !imode
ENDDO ! jbnd
DO j = 1, 3
@ -587,7 +620,7 @@
ENDDO
ENDDO
ENDDO ! ibnd
ENDIF ! iverbsoity
ENDIF ! iverbosity
ENDIF ! ctype
ENDIF ! endif fsthick
ENDDO ! end loop on k
@ -740,6 +773,8 @@
IF (iverbosity == 3) THEN
CALL mp_sum(inv_tau_all_mode, world_comm)
CALL mp_sum(inv_tau_allcb_mode, world_comm)
CALL mp_sum(inv_tau_all_freq, world_comm)
CALL mp_sum(inv_tau_allcb_freq, world_comm)
ENDIF
!
IF (my_pool_id == 0) THEN
@ -823,7 +858,37 @@
ENDDO
ENDDO
CLOSE(iufilibtev_sup)
ENDIF
!
! Compute and save spectral decomposition
OPEN(iufilibtev_sup, FILE = 'inv_tau_freq.fmt', FORM = 'formatted')
WRITE(iufilibtev_sup, '(a)') '# Electron relaxation time [Multiply the relaxation time by 20670.6944033 to get 1/ps] '
WRITE(iufilibtev_sup, '(a)') '# kpt ibnd energy [Ry] freq (meV) relaxation time [Ry]'
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO ifreq = 1, mob_nfreq
IF (inv_tau_all_freq(ifreq, ibnd, ik) > eps160) THEN
WRITE(iufilibtev_sup, '(i8,i6,3E22.12)') ik, ibnd, etf_all(ibnd, ik), ifreq * step * ryd2mev, &
inv_tau_all_freq(ifreq, ibnd, ik)
ENDIF
ENDDO ! ifreq
ENDDO ! ibnd
ENDDO ! ik
CLOSE(iufilibtev_sup)
OPEN(iufilibtev_sup, FILE = 'inv_taucb_freq.fmt', FORM = 'formatted')
WRITE(iufilibtev_sup, '(a)') '# Electron relaxation time [Multiply the relaxation time by 20670.6944033 to get 1/ps] '
WRITE(iufilibtev_sup, '(a)') '# kpt ibnd energy [Ry] freq (meV) relaxation time [Ry]'
DO ik = 1, nktotf
DO ibnd = 1, nbndfst
DO ifreq = 1, mob_nfreq
IF (inv_tau_allcb_freq(ifreq, ibnd, ik) > eps160) THEN
WRITE(iufilibtev_sup, '(i8,i6,3E22.12)') ik, ibnd, etf_all(ibnd, ik), ifreq * step * ryd2mev, &
inv_tau_allcb_freq(ifreq, ibnd, ik)
ENDIF
ENDDO ! ifreq
ENDDO ! ibnd
ENDDO ! ik
CLOSE(iufilibtev_sup)
ENDIF ! iverbosity
!
ENDIF ! master
!

13
EPW/src/low_lvl.f90
View File

@ -1248,7 +1248,7 @@
INTEGER, INTENT(in) :: v(size_v)
!! Vector to bisect
INTEGER, INTENT(in) :: n_intval
!! Number of intervals
!! Number of intervals indexes (so there are (n_intval - 1) number of intervals)
INTEGER, INTENT(out) :: val_intval(n_intval)
!! Value of the first element of each intervals
INTEGER, INTENT(out) :: pos_intval(n_intval)
@ -1270,7 +1270,7 @@
nkl = size_v / (n_intval - 1)
nkr = size_v - nkl * (n_intval - 1)
!
! the reminder goes to the first nkr intervals (0...nkr-1)
! The reminder goes to the first nkr intervals (0...nkr-1)
!
DO i = 1, n_intval
pos_intval(i) = nkl * (i - 1)
@ -1278,6 +1278,13 @@
IF (i >= nkr) pos_intval(i) = pos_intval(i) + nkr
ENDDO
!
! In case the reminder is 0
IF (nkr == 0) THEN
DO i = 1, (n_intval - 1)
pos_intval(i) = pos_intval(i) + 1
ENDDO
ENDIF
!
DO i = 1, n_intval
val_intval(i) = v(pos_intval(i))
ENDDO
@ -1303,7 +1310,7 @@
INTEGER, INTENT(in) :: v(size_v)
!! Vector to bisect
INTEGER, INTENT(in) :: n_intval
!! Number of intervals for the pre search
!! Number of intervals indexes (so there are (n_intval - 1) number of intervals
INTEGER, INTENT(in) :: val_intval(n_intval)
!! Value of the first element of each intervals
INTEGER, INTENT(in) :: pos_intval(n_intval)

20
EPW/src/printing.f90
View File

@ -1269,6 +1269,26 @@
CLOSE(iufilFS)
!
ENDDO
!
! HP: Write in .frmsf format compatible with fermisurfer program
WRITE(name1, '(a, a3, a6)') TRIM(prefix), '.fs', '.frmsf'
OPEN(iufilFS, FILE = name1, STATUS = 'unknown', FORM = 'formatted', IOSTAT = ios)
IF (ios /= 0) CALL errore('plot_fermisurface', 'error opening file ' // name1, iufilFS)
WRITE(iufilFS, '(3i5)') nkf1, nkf2, nkf3
WRITE(iufilFS, '(i5)') 1
WRITE(iufilFS, '(i5)') ibndmax - ibndmin + 1
WRITE(iufilFS, '(3f12.6)') (bg(i, 1), i = 1, 3)
WRITE(iufilFS, '(3f12.6)') (bg(i, 2), i = 1, 3)
WRITE(iufilFS, '(3f12.6)') (bg(i, 3), i = 1, 3)
IF (mp_mesh_k) THEN
WRITE(iufilFS, '(6f12.6)') (((etf_all(ibnd, 2 * bztoibz(ik) - 1) - ef) * ryd2ev, &
ik = 1, nkf1 * nkf2 * nkf3), ibnd = ibndmin, ibndmax)
ELSE
WRITE(iufilFS, '(6f12.6)') (((etf_all(ibnd, ik * 2 - 1) - ef) * ryd2ev, &
ik = 1, nktotf), ibnd = ibndmin, ibndmax)
ENDIF
CLOSE(iufilFS)
!
ENDIF
CALL mp_barrier(inter_pool_comm)
!

14
EPW/src/supercond.f90
View File

@ -609,6 +609,20 @@
WRITE(iufillambdaFS, '(6f12.6)') (lambda_k(ixkff(ik), ibnd), ik = 1, nkf1 * nkf2 * nkf3)
CLOSE(iufillambdaFS)
ENDDO
! HP: Write in .frmsf format compatible with fermisurfer program
WRITE(name1, '(a, a13)') TRIM(prefix), '.lambda.frmsf'
OPEN(UNIT = iufillambdaFS, FILE = name1, STATUS = 'unknown', FORM = 'formatted', IOSTAT = ios)
IF (ios /= 0) CALL errore('evaluate_a2f_lambda', 'error opening file ' // name1, iufillambdaFS)
!
WRITE(iufillambdaFS, '(3i5)') nkf1, nkf2, nkf3
WRITE(iufillambdaFS, '(i5)') 1
WRITE(iufillambdaFS, '(i5)') nbndfs
WRITE(iufillambdaFS, '(3f12.6)') (bg(i, 1), i = 1, 3)
WRITE(iufillambdaFS, '(3f12.6)') (bg(i, 2), i = 1, 3)
WRITE(iufillambdaFS, '(3f12.6)') (bg(i, 3), i = 1, 3)
WRITE(iufillambdaFS, '(6f12.6)') ((ekfs(ibnd, ixkff(ik)) - ef0, ik = 1, nkf1 * nkf2 * nkf3), ibnd = 1, nbndfs)
WRITE(iufillambdaFS, '(6f12.6)') ((lambda_k(ixkff(ik), ibnd), ik = 1, nkf1 * nkf2 * nkf3), ibnd = 1, nbndfs)
CLOSE(iufillambdaFS)
!
ENDIF
!

25
EPW/src/supercond_aniso.f90
View File

@ -661,19 +661,19 @@
!! bose_einstein(w') + fermi_dirac(-w + w')
REAL(KIND = DP) :: absdelta, reldelta, errdelta
!! Errors in supercond. gap
REAL(KIND = DP) :: a2f_
!! Temporary variable for Eliashberg spectral function
REAL(KIND = DP) :: weight
!! Factor in supercond. equations
REAL(KIND = DP), EXTERNAL :: w0gauss
!! The derivative of wgauss: an approximation to the delta function
!!RM - updated
REAL(KIND = DP) :: a2f_tmp(nqstep)
!! Temporary variable for Eliashberg spectral function
!
COMPLEX(KIND = DP) :: az2, ad2, esqrt, root
!! Temporary variables
COMPLEX(KIND = DP), ALLOCATABLE, SAVE :: deltaold(:)
!! supercond. gap from previous iteration
!
a2f_ = zero
inv_degaussq = one / degaussq
!
IF (iter == 1) THEN
@ -741,6 +741,8 @@
znorm(:) = czero
adelta(:, :, :) = czero
aznorm(:, :, :) = czero
!!RM - updated
a2f_tmp(:) = zero
!
CALL fkbounds(nkfs, lower_bnd, upper_bnd)
!
@ -753,15 +755,19 @@
DO jbnd = 1, nbndfs
IF (ABS(ekfs(jbnd, ixkqf(ik, iq0)) - ef0) < fsthick) THEN
!
!!RM - updated
IF (lacon_fly) THEN ! evaluate a2fij on the fly
a2f_tmp(:) = zero
DO imode = 1, nmodes
IF (wf(imode, iq0) > eps_acustic) THEN
DO iwph = 1, nqstep
weight = w0gauss((wsph(iwph) - wf(imode, iq0)) * inv_degaussq, 0) * inv_degaussq
a2f_ = weight * dosef * g2(ik, iq, ibnd, jbnd, imode)
a2f_tmp(iwph) = a2f_tmp(iwph) + weight * dosef * g2(ik, iq, ibnd, jbnd, imode)
ENDDO ! iwph
ENDIF ! wf
ENDDO ! imode
ELSE
a2f_tmp(:) = a2fij(ik, iq, ibnd, jbnd, :)
ENDIF ! lacon_fly
!
weight = wqf(iq) * w0g(jbnd, ixkqf(ik, iq0)) / dosef
@ -778,11 +784,9 @@
ELSE
esqrt = aznormp(jbnd, ixkqf(ik, iq0), i) / root
ENDIF
IF (lacon_fly) THEN
esqrt = esqrt * weight * gp(iw, iwp) * a2f_
ELSE
esqrt = esqrt * weight * gp(iw, iwp) * a2fij(ik, iq, ibnd, jbnd, iwp)
ENDIF
!!RM - updated
esqrt = esqrt * weight * gp(iw, iwp) * a2f_tmp(iwp)
!
aznorm(ibnd, ik, iw) = aznorm(ibnd, ik, iw) - ws(i) * esqrt
adelta(ibnd, ik, iw) = adelta(ibnd, ik, iw) - adeltap(jbnd, ixkqf(ik, iq0), i) * esqrt
ENDIF
@ -796,7 +800,8 @@
ELSE
esqrt = aznormp(jbnd, ixkqf(ik, iq0), i) / root
ENDIF
esqrt = esqrt * weight * gm(iw, iwp) * a2fij(ik, iq, ibnd, jbnd, iwp)
!!RM - updated
esqrt = esqrt * weight * gm(iw, iwp) * a2f_tmp(iwp)
IF (iw < iwp) THEN
aznorm(ibnd, ik, iw) = aznorm(ibnd, ik, iw) - ws(i) * esqrt
ELSE

28
EPW/src/wan2bloch.f90
View File

@ -1762,8 +1762,6 @@
!! Real space WS index
INTEGER :: iw
!! Wannier function index
INTEGER :: iw2
!! Wannier function index
INTEGER :: irn
!! Combined WS and atom index
INTEGER :: ir_start
@ -1870,8 +1868,8 @@
na = MOD(irn - 1, nat) + 1
!
DO iw = 1, dims
CALL ZAXPY(nrr_k * 3, cfac(iw, na, ir), epmatwp(iw, iw2, :, 3 * (na - 1) + 1:3 * na, ir), 1, &
eptmp(iw, iw2, :, 3 * (na - 1) + 1:3 * na), 1)
CALL ZAXPY(nrr_k * 3, cfac(iw, na, ir), epmatwp(iw, :, :, 3 * (na - 1) + 1:3 * na, ir), 1, &
eptmp(iw, :, :, 3 * (na - 1) + 1:3 * na), 1)
ENDDO
ENDDO
ELSE ! use_ws
@ -2245,7 +2243,7 @@
!! Is equal to the number of atoms if use_ws == .TRUE. or 1 otherwise.
INTEGER, INTENT(in) :: irvec_g(3, nrr_g)
!! Coordinates of WS points
INTEGER, INTENT(in) :: ndegen_g(nrr_g, nat, dims)
INTEGER, INTENT(in) :: ndegen_g(dims, nrr_g, nat)
!! Number of degeneracy of WS points
INTEGER, INTENT(in) :: nbnd
!! Number of bands
@ -2267,8 +2265,6 @@
!! Ending ir for this pool
INTEGER :: iw
!! Counter on Wannier functions
INTEGER :: iw2
!! Counter on Wannier functions
INTEGER :: iunepmatwp2
!! Return the file unit
INTEGER :: na
@ -2289,7 +2285,7 @@
!
REAL(KIND = DP) :: rdotk
!! Exponential for the FT
COMPLEX(KIND = DP) :: cfac(nat, nrr_g, dims)
COMPLEX(KIND = DP) :: cfac(nrr_g, dims)
!! Factor for the FT
COMPLEX(KIND = DP), ALLOCATABLE :: epmatw(:, :, :)
!! El-ph matrix elements
@ -2313,17 +2309,18 @@
IF (ierr /= 0) CALL errore('ephwan2blochp_mem', 'error in MPI_FILE_OPEN', 1)
#endif
!
cfac(:, :, :) = czero
cfac(:, :) = czero
!
IF (use_ws) THEN
na = (imode - 1) / 3 + 1
!
DO ir = ir_start, ir_stop
!
! note xxq is assumed to be already in cryst coord
rdotk = twopi * DOT_PRODUCT(xxq, DBLE(irvec_g(:, ir)))
na = (imode - 1) / 3 + 1
DO iw = 1, dims
IF (ndegen_g(ir, na, iw) > 0) &
cfac(na, ir, iw) = EXP(ci * rdotk) / DBLE(ndegen_g(ir, na, iw))
IF (ndegen_g(iw, ir, na) > 0) &
cfac(ir, iw) = EXP(ci * rdotk) / DBLE(ndegen_g(iw, ir, na))
ENDDO
ENDDO
ELSE
@ -2331,7 +2328,7 @@
!
! note xxq is assumed to be already in cryst coord
rdotk = twopi * DOT_PRODUCT(xxq, DBLE(irvec_g(:, ir)))
cfac(1, ir, 1) = EXP(ci * rdotk) / DBLE(ndegen_g(ir, 1, 1))
cfac(ir, 1) = EXP(ci * rdotk) / DBLE(ndegen_g(1, ir, 1))
ENDDO
ENDIF
!
@ -2376,12 +2373,11 @@
#endif
!
IF (use_ws) THEN
na = (imode - 1) / 3 + 1
DO iw = 1, dims
CALL ZAXPY(nbnd * nrr_k, cfac(na, ir, iw), epmatw(iw, :, :), 1, epmatf(iw, :, :), 1)
CALL ZAXPY(nbnd * nrr_k, cfac(ir, iw), epmatw(iw, :, :), 1, epmatf(iw, :, :), 1)
ENDDO
ELSE
CALL ZAXPY(nbnd * nbnd * nrr_k, cfac(1, ir, 1), epmatw, 1, epmatf, 1)
CALL ZAXPY(nbnd * nbnd * nrr_k, cfac(ir, 1), epmatw, 1, epmatf, 1)
ENDIF
!
ENDDO

607
test-suite/epw_mob_ibte_sym/benchmark.out.git.inp=epw9.in.args=5 Normal file
View File

@ -0,0 +1,607 @@
``:oss/
`.+s+. .+ys--yh+ `./ss+.
-sh//yy+` +yy +yy -+h+-oyy
-yh- .oyy/.-sh. .syo-.:sy- /yh
`.-.` `yh+ -oyyyo. `/syys: oys `.`
`/+ssys+-` `sh+ ` oys` .:osyo`
-yh- ./syyooyo` .sys+/oyo--yh/
`yy+ .-:-. `-/+/:` -sh-
/yh. oys
``..---hho---------` .---------..` `.-----.` -hd+---.
`./osmNMMMMMMMMMMMMMMMs. +NNMMMMMMMMNNmh+. yNMMMMMNm- oNMMMMMNmo++:`
+sy--/sdMMMhyyyyyyyNMMh- .oyNMMmyyyyyhNMMm+` -yMMMdyyo:` .oyyNMMNhs+syy`
-yy/ /MMM+.`-+/``mMMy- `mMMh:`````.dMMN:` `MMMy-`-dhhy```mMMy:``+hs
-yy+` /MMMo:-mMM+`-oo/. mMMh: `dMMN/` dMMm:`dMMMMy..MMMo-.+yo`
.sys`/MMMMNNMMMs- mMMmyooooymMMNo: oMMM/sMMMMMM++MMN//oh:
`sh+/MMMhyyMMMs- `-` mMMMMMMMMMNmy+-` -MMMhMMMsmMMmdMMd/yy+
`-/+++oyy-/MMM+.`/hh/.`mNm:` mMMd+/////:-.` NMMMMMd/:NMMMMMy:/yyo/:.`
+os+//:-..-oMMMo:--:::-/MMMo. .-mMMd+---` hMMMMN+. oMMMMMo. `-+osyso:`
syo `mNMMMMMNNNNNNNNMMMo.oNNMMMMMNNNN:` +MMMMs:` dMMMN/` ``:syo
/yh` :syyyyyyyyyyyyyyyy+.`+syyyyyyyyo:` .oyys:` .oyys:` +yh
-yh- ```````````````` ````````` `` `` oys
-+h/------------------------::::::::://////++++++++++++++++++++++///////::::/yd:
shdddddddddddddddddddddddddddddhhhhhhhhyyyyyssssssssssssssssyyyyyyyhhhhhhhddddh`
S. Ponce, E. R. Margine, C. Verdi, and F. Giustino,
Comput. Phys. Commun. 209, 116 (2016)
Program EPW v.5.4 starts on 14Oct2021 at 18: 4:42
This program is part of the open-source Quantum ESPRESSO suite
for quantum simulation of materials; please cite
"P. Giannozzi et al., J. Phys.:Condens. Matter 21 395502 (2009);
"P. Giannozzi et al., J. Phys.:Condens. Matter 29 465901 (2017);
"P. Giannozzi et al., J. Chem. Phys. 152 154105 (2020);
URL http://www.quantum-espresso.org",
in publications or presentations arising from this work. More details at
http://www.quantum-espresso.org/quote
Parallel version (MPI), running on 1 processors
MPI processes distributed on 1 nodes
22693 MiB available memory on the printing compute node when the environment starts
Waiting for input...
Reading input from standard input
WARNING: The use of etf_mem == 3 has been tested and validated for cubic and hexagonal materials.
For other materials, use with care and possibly use etf_mem == 1.
Reading supplied temperature list.
------------------------------------------------------------------------
RESTART - RESTART - RESTART - RESTART
Restart is done without reading PWSCF save file.
Be aware that some consistency checks are therefore not done.
------------------------------------------------------------------------
--
bravais-lattice index = 0
lattice parameter (a_0) = 0.0000 a.u.
unit-cell volume = 0.0000 (a.u.)^3
number of atoms/cell = 0
number of atomic types = 0
kinetic-energy cut-off = 0.0000 Ry
charge density cut-off = 0.0000 Ry
Exchange-correlation= not set
( -1 -1 -1 -1 -1 -1 -1)
celldm(1)= 0.00000 celldm(2)= 0.00000 celldm(3)= 0.00000
celldm(4)= 0.00000 celldm(5)= 0.00000 celldm(6)= 0.00000
crystal axes: (cart. coord. in units of a_0)
a(1) = ( 0.0000 0.0000 0.0000 )
a(2) = ( 0.0000 0.0000 0.0000 )
a(3) = ( 0.0000 0.0000 0.0000 )
reciprocal axes: (cart. coord. in units 2 pi/a_0)
b(1) = ( 0.0000 0.0000 0.0000 )
b(2) = ( 0.0000 0.0000 0.0000 )
b(3) = ( 0.0000 0.0000 0.0000 )
Atoms inside the unit cell:
Cartesian axes
site n. atom mass positions (a_0 units)
No symmetry!
G cutoff = 0.0000 ( 0 G-vectors) FFT grid: ( 0, 0, 0)
number of k points= 0
cart. coord. in units 2pi/a_0
EPW : 0.00s CPU 0.00s WALL
EPW : 0.00s CPU 0.00s WALL
-------------------------------------------------------------------
Using si.ukk from disk
-------------------------------------------------------------------
Symmetries of Bravais lattice: 48
Symmetries of crystal: 48
Do not need to read .epb files; read .fmt files
Band disentanglement is used: nbndsub = 8
Use zone-centred Wigner-Seitz cells
Number of WS vectors for electrons 279
Number of WS vectors for phonons 19
Number of WS vectors for electron-phonon 19
Maximum number of cores for efficient parallelization 114
Results may improve by using use_ws == .TRUE.
Reading Hamiltonian, Dynamical matrix and EP vertex in Wann rep from file
Finished reading Wann rep data from file
===================================================================
Memory usage: VmHWM = 15Mb
VmPeak = 329Mb
===================================================================
Number of k-points inside fsthick * 1.2 in the full BZ: 63
Size of k point mesh for interpolation: 10
Max number of k points per pool: 10
Fermi energy coarse grid = 6.255484 eV
===================================================================
Fermi energy is read from the input file: Ef = 6.863669 eV
===================================================================
ibndmin = 5 ebndmin = 6.964 eV
ibndmax = 6 ebndmax = 7.207 eV
Number of ep-matrix elements per pool : 120 ~= 0.94 Kb (@ 8 bytes/ DP)
File ./Fepmatkqcb1/si.epmatkqcb1_0 deleted, as requested
File ./Fsparsecb/sparsecb_0 deleted, as requested
A selecq.fmt file was found but re-created because selecqread == .FALSE.
Number selected, total 5 14
Number selected, total 10 27
Number selected, total 15 131
Number selected, total 20 157
Number selected, total 25 585
Number selected, total 30 703
Number selected, total 35 728
Number selected, total 40 792
Number selected, total 45 816
Number selected, total 50 870
Number selected, total 55 884
Number selected, total 60 898
Number selected, total 65 948
Number selected, total 70 961
Number selected, total 75 1003
Number selected, total 80 1040
Number selected, total 85 1082
Number selected, total 90 1107
Number selected, total 95 1196
Number selected, total 100 1253
Number selected, total 105 1611
We only need to compute 106 q-points
Valence band maximum = 3.799268 eV
Conduction band minimum = 6.963669 eV
Temperature 300.000 K
Mobility CB Fermi level = 6.532632 eV
===================================================================
Scattering rate for IBTE
===================================================================
Restart and restart_step inputs deactivated (restart point at every q-points).
No intermediate mobility will be shown.
Fermi Surface thickness = 0.400000 eV
This is computed with respect to the fine Fermi level 6.863669 eV
Only states between 6.463669 eV and 7.263669 eV will be included
Save matrix elements larger than threshold: 0.837244941701E-23
Progression iq (fine) = 5/ 106
Adaptative smearing = Min: 7.602499 meV
Max: 346.323656 meV
Progression iq (fine) = 10/ 106
Adaptative smearing = Min: 6.500839 meV
Max: 132.547660 meV
Progression iq (fine) = 15/ 106
Adaptative smearing = Min: 48.340388 meV
Max: 214.771190 meV
Progression iq (fine) = 20/ 106
Adaptative smearing = Min: 42.084725 meV
Max: 283.841697 meV
Progression iq (fine) = 25/ 106
Adaptative smearing = Min: 43.313615 meV
Max: 220.412235 meV
Progression iq (fine) = 30/ 106
Adaptative smearing = Min: 46.475521 meV
Max: 215.881949 meV
Progression iq (fine) = 35/ 106
Adaptative smearing = Min: 6.615569 meV
Max: 215.881949 meV
Progression iq (fine) = 40/ 106
Adaptative smearing = Min: 65.168124 meV
Max: 346.104135 meV
Progression iq (fine) = 45/ 106
Adaptative smearing = Min: 6.296537 meV
Max: 353.697477 meV
Progression iq (fine) = 50/ 106
Adaptative smearing = Min: 6.296537 meV
Max: 353.926338 meV
Progression iq (fine) = 55/ 106
Adaptative smearing = Min: 41.972377 meV
Max: 353.926338 meV
Progression iq (fine) = 60/ 106
Adaptative smearing = Min: 41.088312 meV
Max: 353.222822 meV
Progression iq (fine) = 65/ 106
Adaptative smearing = Min: 43.862625 meV
Max: 352.477995 meV
Progression iq (fine) = 70/ 106
Adaptative smearing = Min: 42.660539 meV
Max: 353.018064 meV
Progression iq (fine) = 75/ 106
Adaptative smearing = Min: 42.660539 meV
Max: 285.169721 meV
Progression iq (fine) = 80/ 106
Adaptative smearing = Min: 7.660970 meV
Max: 284.181305 meV
Progression iq (fine) = 85/ 106
Adaptative smearing = Min: 43.282283 meV
Max: 219.462039 meV
Progression iq (fine) = 90/ 106
Adaptative smearing = Min: 47.687616 meV
Max: 353.842363 meV
Progression iq (fine) = 95/ 106
Adaptative smearing = Min: 48.890472 meV
Max: 213.355380 meV
Progression iq (fine) = 100/ 106
Adaptative smearing = Min: 47.034421 meV
Max: 215.393039 meV
Progression iq (fine) = 105/ 106
Adaptative smearing = Min: 40.047972 meV
Max: 284.291242 meV
300.000 6.5326 0.999999E+13
epmatkqread automatically changed to .TRUE. as all scattering have been computed.
===================================================================
Memory usage: VmHWM = 21Mb
VmPeak = 355Mb
===================================================================
Number of elements per core 220
Symmetry mapping finished
=============================================================================================
BTE in the self-energy relaxation time approximation (SERTA)
=============================================================================================
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 0.00000E+00 0.222638E+04 0.686733E-18 0.137347E-17
0.686733E-18 0.222638E+04 0.000000E+00
0.137347E-17 0.000000E+00 0.222638E+04
=============================================================================================
Start solving iterative Boltzmann Transport Equation
=============================================================================================
Iteration number: 1
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 0.00000E+00 0.215699E+04 0.686733E-18 -0.137347E-17
0.686733E-18 0.215699E+04 0.000000E+00
-0.137347E-17 0.000000E+00 0.215699E+04
0.215699E+04 Max error
Iteration number: 2
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 -0.10340E-24 0.215915E+04 0.000000E+00 -0.137347E-17
0.000000E+00 0.215915E+04 0.000000E+00
-0.137347E-17 0.000000E+00 0.215915E+04
0.216311E+01 Max error
Iteration number: 3
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 -0.10340E-24 0.215909E+04 0.686733E-18 0.137347E-17
0.686733E-18 0.215909E+04 0.000000E+00
0.137347E-17 0.000000E+00 0.215909E+04
0.674494E-01 Max error
Iteration number: 4
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 -0.10340E-24 0.215909E+04 0.000000E+00 0.137347E-17
0.000000E+00 0.215909E+04 0.000000E+00
0.137347E-17 0.000000E+00 0.215909E+04
0.210444E-02 Max error
Iteration number: 5
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 -0.20680E-24 0.215909E+04 0.000000E+00 -0.137347E-17
0.686733E-18 0.215909E+04 0.000000E+00
-0.137347E-17 0.000000E+00 0.215909E+04
0.657414E-04 Max error
Iteration number: 6
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 0.00000E+00 0.215909E+04 0.000000E+00 0.137347E-17
0.000000E+00 0.215909E+04 0.000000E+00
0.137347E-17 0.000000E+00 0.215909E+04
0.205918E-05 Max error
Iteration number: 7
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.10000E+14 0.10340E-24 0.215909E+04 0.686733E-18 -0.137347E-17
0.686733E-18 0.215909E+04 0.000000E+00
-0.137347E-17 0.000000E+00 0.215909E+04
0.648611E-07 Max error
Printing the state-resolved mobility to mobility_nk.fmt file
=============================================================================================
BTE in the SERTA with B-field
=============================================================================================
Number of contributing elements for the master core 10560
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 0.48468E-25 0.222640E+04 0.727949E-07 0.000000E+00
-0.727949E-07 0.222640E+04 0.214606E-19
0.686740E-18 -0.107303E-18 0.222640E+04
0.673156E+02 Max error
=============================================================================================
Start solving iterative Boltzmann Transport Equation with B-field
=============================================================================================
Iteration number: 1
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.21493E-24 0.215701E+04 0.669155E-07 0.137348E-17
-0.669155E-07 0.215701E+04 -0.772582E-18
0.000000E+00 -0.214606E-18 0.215701E+04
0.693921E+02 Max error
Iteration number: 2
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.14432E-25 0.215917E+04 0.581196E-07 -0.686740E-18
-0.581196E-07 0.215917E+04 0.000000E+00
0.686740E-18 0.879885E-18 0.215917E+04
0.216313E+01 Max error
Iteration number: 3
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.12962E-24 0.215911E+04 0.593862E-07 -0.137348E-17
-0.593862E-07 0.215911E+04 -0.772582E-18
-0.343370E-17 -0.193146E-18 0.215911E+04
0.674500E-01 Max error
Iteration number: 4
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.87106E-25 0.215911E+04 0.592587E-07 0.206022E-17
-0.592587E-07 0.215911E+04 0.145932E-17
-0.686740E-18 -0.643818E-19 0.215911E+04
0.210446E-02 Max error
Iteration number: 5
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.88231E-25 0.215911E+04 0.592699E-07 0.686740E-18
-0.592699E-07 0.215911E+04 -0.128764E-17
0.137348E-17 0.815503E-18 0.215911E+04
0.657420E-04 Max error
Iteration number: 6
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.12535E-25 0.215911E+04 0.592690E-07 0.686740E-18
-0.592690E-07 0.215911E+04 -0.171685E-18
0.137348E-17 0.214606E-18 0.215911E+04
0.205920E-05 Max error
Iteration number: 7
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.20324E-24 0.215911E+04 0.592691E-07 0.343370E-17
-0.592691E-07 0.215911E+04 -0.343370E-18
0.137348E-17 0.429212E-19 0.215911E+04
0.648602E-07 Max error
Iteration number: 8
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.16983E-24 0.215911E+04 0.592691E-07 -0.686740E-18
-0.592691E-07 0.215911E+04 0.686740E-18
-0.137348E-17 -0.622358E-18 0.215911E+04
0.206728E-08 Max error
Iteration number: 9
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 -0.45012E-25 0.215911E+04 0.592691E-07 -0.686740E-18
-0.592691E-07 0.215911E+04 0.171685E-18
-0.206022E-17 0.236067E-18 0.215911E+04
0.654836E-10 Max error
Iteration number: 10
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 0.32637E-25 0.215911E+04 0.592691E-07 0.206022E-17
-0.592691E-07 0.215911E+04 0.171685E-18
-0.274696E-17 -0.472133E-18 0.215911E+04
0.136424E-11 Max error
Iteration number: 11
=============================================================================================
Temp Fermi Elec density Population SR Elec mobility
[K] [eV] [cm^-3] [e per cell] [cm^2/Vs]
=============================================================================================
300.000 6.5326 0.99999E+13 0.10190E-24 0.215911E+04 0.592691E-07 0.686740E-18
-0.592691E-07 0.215911E+04 -0.343370E-18
0.686740E-18 0.214606E-18 0.215911E+04
0.000000E+00 Max error
=============================================================================================
Summary and Hall factor
=============================================================================================
=============================================================================================
BTE in the self-energy relaxation time approximation (SERTA)
=============================================================================================
Temperature: 300.0000 K
Conductivity tensor without magnetic field | with magnetic field [Siemens/m]
0.35671E+00 0.11003E-21 0.22005E-21 | 0.35671E+00 -0.11663E-10 0.11003E-21
0.11003E-21 0.35671E+00 0.00000E+00 | 0.11663E-10 0.35671E+00 -0.17192E-22
0.22005E-21 0.00000E+00 0.35671E+00 | 0.00000E+00 0.34383E-23 0.35671E+00
Mobility tensor without magnetic field | with magnetic field [cm^2/Vs]
0.22264E+04 0.68674E-18 0.13735E-17 | 0.22264E+04 -0.72795E-07 0.68674E-18
0.68674E-18 0.22264E+04 0.00000E+00 | 0.72795E-07 0.22264E+04 -0.10730E-18
0.13735E-17 0.00000E+00 0.22264E+04 | 0.00000E+00 0.21461E-19 0.22264E+04
Hall factor
0.119268E+05 -0.389959E-06 -0.110365E-16
0.389959E-06 0.119268E+05 -0.574818E-18
-0.147154E-16 0.114964E-18 0.119268E+05
=============================================================================================
BTE
=============================================================================================
Temperature: 300.0000 K
Conductivity tensor without magnetic field | with magnetic field [Siemens/m]
0.34592E+00 0.11003E-21 -0.22005E-21 | 0.34592E+00 -0.94958E-11 0.11003E-21
0.11003E-21 0.34592E+00 0.00000E+00 | 0.94958E-11 0.34592E+00 0.34383E-22
-0.22005E-21 0.00000E+00 0.34592E+00 | 0.11003E-21 -0.55013E-22 0.34592E+00
Mobility tensor without magnetic field | with magnetic field [cm^2/Vs]
0.21591E+04 0.68674E-18 -0.13735E-17 | 0.21591E+04 -0.59269E-07 0.68674E-18
0.68674E-18 0.21591E+04 0.00000E+00 | 0.59269E-07 0.21591E+04 0.21461E-18
-0.13735E-17 0.00000E+00 0.21591E+04 | 0.68674E-18 -0.34337E-18 0.21591E+04
Hall factor
0.122985E+05 -0.337602E-06 0.195587E-16
0.337602E-06 0.122985E+05 0.122242E-17
0.195587E-16 -0.195587E-17 0.122985E+05
Unfolding on the coarse grid
elphon_wrap : 0.01s CPU 0.01s WALL ( 1 calls)
INITIALIZATION:
Electron-Phonon interpolation
ephwann : 0.89s CPU 1.28s WALL ( 1 calls)
ep-interp : 0.55s CPU 0.94s WALL ( 106 calls)
DynW2B : 0.00s CPU 0.00s WALL ( 106 calls)
HamW2B : 0.10s CPU 0.11s WALL ( 2808 calls)
ephW2Bp : 0.39s CPU 0.73s WALL ( 106 calls)
ephW2B : 0.02s CPU 0.03s WALL ( 156 calls)
print_ibte : 0.02s CPU 0.02s WALL ( 106 calls)
vmewan2bloch : 0.05s CPU 0.06s WALL ( 418 calls)
vmewan2bloch : 0.05s CPU 0.06s WALL ( 418 calls)
Total program execution
EPW : 0.90s CPU 1.29s WALL
===============================================================================
The functionality-dependent EPW.bib file was created with suggested citations.
Please consider citing the papers listed in EPW.bib.
===============================================================================

92
test-suite/epw_mob_ibte_sym/epw9.in Normal file
View File

@ -0,0 +1,92 @@
--
&inputepw
prefix = 'si'
amass(1) = 28.0855
outdir = './'
elph = .true.
epbwrite = .false.
epbread = .false.
epwwrite = .false.
epwread = .true.
etf_mem = 3
mp_mesh_k = .true.
vme = 'wannier'
use_ws = .false.
iverbosity = 3
mob_maxfreq = 70
mob_nfreq = 280
scattering = .true.
scattering_serta = .true.
int_mob = .false.
carrier = .true.
ncarrier = 1E13
iterative_bte = .true.
epmatkqread = .false.
mob_maxiter = 20
broyden_beta= 1.0
bfieldx = 0.0d0
bfieldy = 0.0d0
bfieldz = 1.0d-10
nstemp = 1
temps = 300
restart = .true.
restart_step = 5
selecqread = .false.
lpolar = .false.
lphase = .true.
nbndsub = 8
wannierize = .false.
num_iter = 1500
iprint = 2
dis_win_max = 18
dis_froz_max= 8.5
proj(1) = 'Si : sp3'
wdata(1) = 'bands_plot = .true.'
wdata(2) = 'begin kpoint_path'
wdata(3) = 'L 0.50 0.00 0.00 G 0.00 0.00 0.00'
wdata(4) = 'G 0.00 0.00 0.00 X 0.50 0.50 0.00'
wdata(5) = 'end kpoint_path'
wdata(6) = 'bands_plot_format = gnuplot'
wdata(7) = 'guiding_centres = .true.'
wdata(8) = 'dis_num_iter = 500'
wdata(9) = 'num_print_cycles = 10'
wdata(10) = 'dis_mix_ratio = 1.0'
wdata(11) = 'conv_tol = 1E-9'
wdata(12) = 'conv_window = 4'
elecselfen = .false.
phonselfen = .false.
a2f = .false.
fsthick = 0.4
degaussw = 0.0 ! eV
dvscf_dir = './save/'
efermi_read = .true
fermi_energy = 6.863669
nkf1 = 12
nkf2 = 12
nkf3 = 12
nqf1 = 12
nqf2 = 12
nqf3 = 12
nk1 = 6
nk2 = 6
nk3 = 6
nq1 = 2
nq2 = 2
nq3 = 2
/

2
test-suite/jobconfig
View File

@ -159,7 +159,7 @@ inputs_args = ('scf.in', '1'), ('ph.in', '2'), ('q2r.in', '4'), ('scf.in', '1'),
[epw_mob_ibte_sym/]
program = EPW
inputs_args = ('scf.in', '1'), ('ph.in', '2'), ('q2r.in', '4'), ('scf.in', '1'), ('nscf.in', '1'), ('epw1.in', '3'), ('epw2.in', '3'), ('epw3.in', '3'), ('epw4.in', '5'), ('epw5.in', '3'), ('epw6.in', '5'), ('epw7.in', '3'), ('epw8.in', '5')
inputs_args = ('scf.in', '1'), ('ph.in', '2'), ('q2r.in', '4'), ('scf.in', '1'), ('nscf.in', '1'), ('epw1.in', '3'), ('epw2.in', '3'), ('epw3.in', '3'), ('epw4.in', '5'), ('epw5.in', '3'), ('epw6.in', '5'), ('epw7.in', '3'), ('epw8.in', '5'), ('epw9.in', '5')
[epw_mob_polar/]
program = EPW