diff --git a/Doc/release-notes b/Doc/release-notes index f1e0f9dd9..63fae4185 100644 --- a/Doc/release-notes +++ b/Doc/release-notes @@ -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 diff --git a/EPW/ZG/README b/EPW/ZG/README index cdeb22bd0..353c511ae 100644 --- a/EPW/ZG/README +++ b/EPW/ZG/README @@ -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.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.out". + ii) To generate the ZG-displacement run "/path_to_your_espresso/bin/ZG.x 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.out". The output files of interest are: diff --git a/EPW/ZG/doc/README b/EPW/ZG/doc/README deleted file mode 100644 index 573d988ad..000000000 --- a/EPW/ZG/doc/README +++ /dev/null @@ -1,2 +0,0 @@ -For the exercises in the tutorial please download material from: -https://docs.epw-code.org/doc/Virtual2021.html diff --git a/EPW/ZG/doc/input_variables.pdf b/EPW/ZG/doc/input_variables.pdf deleted file mode 100644 index e98cd66e8..000000000 Binary files a/EPW/ZG/doc/input_variables.pdf and /dev/null differ diff --git a/EPW/ZG/doc/tutorial.pdf b/EPW/ZG/doc/tutorial.pdf deleted file mode 100644 index 992d94f27..000000000 Binary files a/EPW/ZG/doc/tutorial.pdf and /dev/null differ diff --git a/EPW/ZG/example.tar.gz b/EPW/ZG/example.tar.gz deleted file mode 100644 index 34d7a692b..000000000 Binary files a/EPW/ZG/example.tar.gz and /dev/null differ diff --git a/EPW/ZG/src/ZG.f90 b/EPW/ZG/src/ZG.f90 index 3fe1454e5..48516a2fc 100644 --- a/EPW/ZG/src/ZG.f90 +++ b/EPW/ZG/src/ZG.f90 @@ -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 diff --git a/EPW/ZG/src/disca.f90 b/EPW/ZG/src/disca.f90 index 51bf6ed54..bb2a0b846 100644 --- a/EPW/ZG/src/disca.f90 +++ b/EPW/ZG/src/disca.f90 @@ -27,9 +27,9 @@ PROGRAM diff_sca !----------------------------------------------------------------------- !----------------------------------------------------------------------- !----------------------------------------------------------------------- - !! authors: Marios Zacharias, Pantelis C. Kelires, Feliciano Giustino + !! authors: Marios Zacharias and Feliciano Giustino !! acknowledgement: Hyungjun Lee for help packaging this release - !! version: v0.1 + !! version: v0.2 !! license: GNU ! ! This program generates diffuse scattering maps for a list of @@ -145,10 +145,10 @@ PROGRAM diff_sca CHARACTER(LEN = 256) :: flfrc, filename CHARACTER(LEN = 10) :: asr LOGICAL :: has_zstar, q_in_cryst_coord, eigen_similarity, loto_disable - COMPLEX(DP), ALLOCATABLE :: dyn(:,:,:,:), dyn_blk(:,:,:,:), frc_ifc(:,:,:,:) - COMPLEX(DP), ALLOCATABLE :: z(:,:) - REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), freq(:,:), wq(:) - REAL(DP), ALLOCATABLE:: q_super(:,:), q_strft(:,:), q_strft_part(:,:) + COMPLEX(DP), ALLOCATABLE :: dyn(:, :, :, :), dyn_blk(:, :, :, :) + COMPLEX(DP), ALLOCATABLE :: z(:,:), frc_ifc(:, :, :, :) + REAL(DP), ALLOCATABLE:: tau(:, :), q(:, :), w2(:, :), freq(:,:), wq(:) + REAL(DP), ALLOCATABLE:: q_super(:, :), q_strft(:, :), q_strft_part(:, :) INTEGER, ALLOCATABLE:: ityp(:), itau_blk(:) REAL(DP) :: omega,alat, &! cell parameters and volume at_blk(3, 3), bg_blk(3, 3), &! original cell @@ -157,7 +157,7 @@ PROGRAM diff_sca amass(ntypx), &! atomic masses amass_blk(ntypx), &! original atomic masses atws(3, 3), &! lattice vector for WS initialization - rws(0:3, nrwsx) ! nearest neighbor list, rws(0,*) = norm^2 + rws(0 : 3, nrwsx) ! nearest neighbor list, rws(0,*) = norm^2 ! INTEGER :: nat, nat_blk, ntyp, ntyp_blk, & l1, l2, l3, &! supercell dimensions @@ -170,17 +170,17 @@ PROGRAM diff_sca ! REAL(DP) :: qhat(3), qh, E, qq REAL(DP) :: delta - REAL(DP), ALLOCATABLE :: xqaux(:,:) + REAL(DP), ALLOCATABLE :: xqaux(:, :) INTEGER, ALLOCATABLE :: nqb(:) INTEGER :: n, i, j, it, nq, nq_super, nq_strft, nqx, na, nb, nqtot LOGICAL, EXTERNAL :: has_xml - INTEGER, ALLOCATABLE :: num_rap_mode(:,:) + INTEGER, ALLOCATABLE :: num_rap_mode(:, :) LOGICAL, ALLOCATABLE :: high_sym(:) LOGICAL :: q_in_band_form ! .... variables for band plotting based on similarity of eigenvalues - COMPLEX(DP), ALLOCATABLE :: tmp_z(:,:) - REAL(DP), ALLOCATABLE :: ABS_similarity(:,:), tmp_w2(:) - COMPLEX(DP), ALLOCATABLE :: f_of_q(:,:,:,:) + COMPLEX(DP), ALLOCATABLE :: tmp_z(:, :) + REAL(DP), ALLOCATABLE :: ABS_similarity(:, :), tmp_w2(:) + COMPLEX(DP), ALLOCATABLE :: f_of_q(:, :, :, :) INTEGER :: location(1), isig CHARACTER(LEN=6) :: int_to_char LOGICAL, ALLOCATABLE :: mask(:) @@ -193,15 +193,18 @@ PROGRAM diff_sca CHARACTER(len=80) :: k_points = 'tpiba' ! mz_b COMPLEX(DP), ALLOCATABLE :: z_nq_zg(:, :, :) ! nomdes,nmodes,nq - REAL(DP), ALLOCATABLE :: q_nq_zg(:, :), qlist_str_f(:, :), qlist_str_f_cart(:, :) ! 3,nq - REAL(DP) :: atmsf_zg_a(ntypx, 5), atmsf_zg_b(ntypx, 5) - LOGICAL :: disca, q_external - LOGICAL :: zero_one_phonon, atom_resolved, full_phonon + REAL(DP), ALLOCATABLE :: q_nq_zg(:, :), qlist_strf(:, :), qlist_strf_cart(:, :) ! 3,nq + REAL(DP) :: atmsf_a(ntypx, 5), atmsf_b(ntypx, 5) + LOGICAL :: disca, q_external, print_raw_data + LOGICAL :: zero_one_phonon, atom_resolved, mode_resolved, full_phonon INTEGER :: dim1, dim2, dim3, nks1, nks2, nks3, nksf1, nksf2, nksf3 INTEGER :: lower_bnd, upper_bnd INTEGER :: T, plane_dir - REAL(DP) :: plane_val + REAL(DP) :: plane_val, eps2 CHARACTER(LEN=3) :: atm_zg(ntypx) + ! inputs for pp_disca, rotate + INTEGER :: nrots, kres1, kres2, col1, col2, Np + REAL(DP) :: kmin, kmax ! mz_e ! NAMELIST /input/ flfrc, amass, asr, at, & @@ -209,16 +212,18 @@ PROGRAM diff_sca & q_in_band_form, q_in_cryst_coord, & & eigen_similarity, na_ifc, fd, point_label_type, nosym, & & loto_2d, loto_disable, & -! mz_b we add the inputs for diffuse scattering - & disca, dim1, dim2, dim3, atm_zg, & - & T, q_external, zero_one_phonon, & + ! mz_b we add the inputs for diffuse scattering + & disca, dim1, dim2, dim3, atm_zg, mode_resolved, & + & T, q_external, zero_one_phonon, print_raw_data, & & nks1, nks2, nks3, plane_val, plane_dir, qstart, qfinal, atom_resolved, & - & full_phonon, atmsf_zg_a, atmsf_zg_b, nksf1, nksf2, nksf3 -! disca --> if true compute the diffuse_scattering -! mz_e + & full_phonon, atmsf_a, atmsf_b, nksf1, nksf2, nksf3, eps2 + ! disca --> if true compute the diffuse_scattering + !we add a new list for the inputs og rotate and pp_disca + NAMELIST /pp_disca/ nrots, kres1, kres2, kmin, kmax, col1, col2, Np + ! mz_e ! CALL mp_startup() - CALL environment_start('MATDYN') + CALL environment_start('DISCA') ! IF (ionode) CALL input_from_file ( ) ! @@ -245,9 +250,9 @@ PROGRAM diff_sca na_ifc = .FALSE. fd = .FALSE. point_label_type = 'SC' - nosym = .false. - loto_2d = .false. - loto_disable = .false. + nosym = .FALSE. + loto_2d = .FALSE. + loto_disable = .FALSE. ! mz_b disca = .TRUE. T = 0 @@ -255,12 +260,14 @@ PROGRAM diff_sca zero_one_phonon = .TRUE. full_phonon = .FALSE. atom_resolved = .FALSE. + mode_resolved = .FALSE. + print_raw_data = .FALSE. dim1 = 0 dim2 = 0 dim3 = 0 atm_zg = "Element" - atmsf_zg_a = 0.d0 - atmsf_zg_b = 0.d0 + atmsf_a = 0.d0 + atmsf_b = 0.d0 plane_val = 0.d0 plane_dir = 3 nks1 = 0 @@ -271,12 +278,25 @@ PROGRAM diff_sca nksf3 = 6 qstart = 1 qfinal = 2 + eps2 = 1.0d-15 + ! + nrots = 1 + kres1 = 250 + kres2 = 250 + kmin = -5 + kmax = 10 + col1 = 1 + col2 = 2 + Np = 100 ! mz_e ! ! IF (ionode) READ (5,input,IOSTAT=ios) CALL mp_bcast(ios, ionode_id, world_comm) CALL errore('disca', 'reading input namelist', ABS(ios)) + IF (ionode) READ (5,pp_disca,IOSTAT=ios) + CALL mp_bcast(ios, ionode_id, world_comm) + CALL errore('pp_disca', 'reading pp_disca namelist', ABS(ios)) CALL mp_bcast(nk1,ionode_id, world_comm) CALL mp_bcast(nk2,ionode_id, world_comm) CALL mp_bcast(nk3,ionode_id, world_comm) @@ -304,9 +324,11 @@ PROGRAM diff_sca CALL mp_bcast(q_external, ionode_id, world_comm) CALL mp_bcast(zero_one_phonon, ionode_id, world_comm) CALL mp_bcast(atom_resolved, ionode_id, world_comm) + CALL mp_bcast(mode_resolved, ionode_id, world_comm) + CALL mp_bcast(print_raw_data, ionode_id, world_comm) CALL mp_bcast(full_phonon, 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(dim1, ionode_id, world_comm) CALL mp_bcast(dim2, ionode_id, world_comm) CALL mp_bcast(dim3, ionode_id, world_comm) @@ -320,6 +342,16 @@ PROGRAM diff_sca CALL mp_bcast(plane_dir, ionode_id, world_comm) CALL mp_bcast(qstart, ionode_id, world_comm) CALL mp_bcast(qfinal, ionode_id, world_comm) + CALL mp_bcast(eps2, 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('disca', & 'loto_2d and loto_disable cannot be both true', 1) @@ -444,30 +476,28 @@ PROGRAM diff_sca ! CALL qpoint_gen2(dim1, dim2, dim3, nq_super, q_super) ! - !DO n = 1, nq_super - ! IF (ionode) WRITE(*,*) "q_AB", q_super(:,n) - !ENDDO - ! CALL mp_bcast(q_super, ionode_id, world_comm) ! CALL cryst_to_cart(nq_super, q_super, bg, +1) ! convert them to Cartesian ! ! we first read nq for structure factor q-points - S - ALLOCATE(qlist_str_f(dim1 * dim2 * dim3 * (nksf1 - nks1) * (nksf2 - nks2) * (nksf3 - nks3), 3)) - ALLOCATE(qlist_str_f_cart(dim1 * dim2 * dim3 * (nksf1 - nks1) * (nksf2 - nks2) * (nksf3 - nks3), 3)) + ALLOCATE(qlist_strf(dim1 * dim2 * dim3 * (nksf1 - nks1) * (nksf2 - nks2) * (nksf3 - nks3), 3)) + ALLOCATE(qlist_strf_cart(dim1 * dim2 * dim3 * (nksf1 - nks1) * (nksf2 - nks2) * (nksf3 - nks3), 3)) IF (ionode) CALL qpoint_gen_map_1(dim1, dim2, dim3, nks1, nks2, nks3, nksf1, nksf2, nksf3, & - plane_val, plane_dir, nq_strft, qlist_str_f, qlist_str_f_cart) + plane_val, plane_dir, nq_strft, qlist_strf, qlist_strf_cart) ! nq_strft is the # of q-points CALL mp_bcast(nq_strft, ionode_id, world_comm) + ! IF (ionode) WRITE (*,*) IF (ionode) WRITE (*,*) "Generating", nq_strft, "q-points for the full map ..." + ! IF (qfinal .GT. nq_strft) CALL errore('disca', 'qfinal larger than full # of q-points', qfinal) ! ALLOCATE (q_strft(3, nq_strft), q_strft_part(3, qfinal - qstart + 1 )) ! IF (ionode) CALL qpoint_gen_map_2(dim1, dim2, dim3, nks1, nks2, nks3, nksf1, nksf2, nksf3, & - plane_val, plane_dir, nq_strft, qlist_str_f, qlist_str_f_cart, q_strft) ! q_strft array with q-points - DEALLOCATE(qlist_str_f, qlist_str_f_cart) + plane_val, plane_dir, nq_strft, qlist_strf, qlist_strf_cart, q_strft) ! q_strft array with q-points + DEALLOCATE(qlist_strf, qlist_strf_cart) ! CALL mp_bcast(q_strft, ionode_id, world_comm) IF (ionode) then @@ -486,7 +516,6 @@ PROGRAM diff_sca ! convert them to Cartesian CALL cryst_to_cart(nq_strft, q_strft_part, bg, +1) ! - ! now generate normal q-points in sets AB ! Now ALLOCATE and read full q-list containing ! supercell q-points and structure factor q-points nq = nq_super + nq_strft @@ -498,6 +527,8 @@ PROGRAM diff_sca DO n = nq_super + 1, nq q(:, n) = q_strft_part(:, n - nq_super) ENDDO + ! + DEALLOCATE(q_strft_part, q_super) ELSE ! mz_ends IF (ionode) READ (5,*) nq @@ -505,10 +536,10 @@ PROGRAM diff_sca ALLOCATE ( q(3,nq) ) IF (.NOT.q_in_band_form) THEN DO n = 1,nq - ! mz_edits + ! mz_edits IF (ionode) READ (5,*) (q(i,n),i= 1,3) -! IF (ionode) READ (5,'(3F10.6)') q(:,n) - ! mz_done + ! IF (ionode) READ (5,'(3F10.6)') q(:,n) + ! mz_done ENDDO CALL mp_bcast(q, ionode_id, world_comm) ! @@ -599,18 +630,10 @@ PROGRAM diff_sca ! mz_b w2 = 0.d0 IF (disca) THEN - ALLOCATE ( z_nq_zg(3*nat,3*nat,nq),q_nq_zg(3,nq)) + ALLOCATE ( z_nq_zg(3 * nat, 3 * nat, nq),q_nq_zg(3, nq)) z_nq_zg(:,:,:)= (0.d0, 0.d0) q_nq_zg(:,:) = 0.d0 ENDIF -! WRITE atomic structure factor from input - !IF (ionode) THEN - ! DO it = 1, nat - ! na = ityp(it) - ! WRITE(*,*) atmsf_zg_a(na,:) - ! WRITE(*,*) atmsf_zg_b(na,:) - ! ENDDO - !ENDIF ! mz_e IF (xmlifc) CALL set_sym(nat, tau, ityp, nspin_mag, m_loc ) @@ -623,7 +646,6 @@ PROGRAM diff_sca CALL fkbounds( nq, lower_bnd, upper_bnd ) ! DO n= lower_bnd, upper_bnd ! 1, nq - !DO n = 1, nq dyn(:,:,:,:) = (0.d0, 0.d0) lo_to_split=.FALSE. @@ -643,7 +665,7 @@ PROGRAM diff_sca 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, & + loto_2d, & epsil, zeu, frc, nr1,nr2,nr3, has_zstar, rws, nrws, na_ifc,f_of_q,fd) IF (.not.loto_2d) THEN @@ -692,22 +714,20 @@ PROGRAM diff_sca 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 ! ENDIF ! - - !mzzz_comments out --> if(iout_dyn.ne.0) CALL WRITE_dyn_on_file(q(1,n),dyn,nat, iout_dyn) - - CALL dyndiag(nat,ntyp,amass,ityp,dyn,w2(1,n),z) - ! mz_b fill a 3D matrix with all eigenvectors - ! mz_e - !!!!! mzzz IF (ionode.and.iout_eig.ne.0) & - !!!!mz & CALL WRITE_eigenvectors(nat,ntyp,amass,ityp,q(1,n),w2(1,n),z,iout_eig) ! + ! Fill a 3D matrix with all eigenvectors + ! + IF (disca) THEN + z_nq_zg(:,:,n) = z(:,:) + q_nq_zg(:,n) = q(:,n) + ENDIF ! Cannot use the small group of \Gamma to analize the symmetry ! of the mode if there is an electric field. ! @@ -739,12 +759,7 @@ PROGRAM diff_sca tmp_z(:,:) = z(:,:) ENDIF ! - IF (disca) THEN - z_nq_zg(:,:,n) = z(:,:) - q_nq_zg(:,n) = q(:,n) - ENDIF ! - !!!!!!!!!mz IF (ionode.and.iout.ne.0) CALL WRITEmodes(nat,q(1,n),w2(1,n),z,iout) ! ENDDO !nq ! @@ -755,30 +770,18 @@ PROGRAM diff_sca DEALLOCATE (tmp_w2, ABS_similarity, mask) IF (eigen_similarity) DEALLOCATE(tmp_z) ! - ! - !ALLOCATE (freq(3*nat, nq)) - !DO n= 1,nq - ! ! freq(i,n) = frequencies in cm^(-1), with negative sign if omega^2 is negative - ! DO i= 1,3*nat - ! freq(i,n)= SQRT(ABS(w2(i,n))) * RY_TO_CMM1 - ! IF (w2(i,n) < 0.0d0) freq(i,n) = -freq(i,n) - ! ENDDO - !ENDDO - ! - ! - ! !mz_b - !IF (ionode) WRITE(*,*) "step1" - IF (disca) CALL diffuse_scattering(nq,nq_super,nq_strft, nat, ntyp, amass, ityp, q_nq_zg, w2, z_nq_zg, & - q_external, zero_one_phonon, full_phonon, atmsf_zg_a, atmsf_zg_b, & - dim1, dim2, dim3, tau, alat, atm_zg, atom_resolved, & - ntypx, at, q_in_cryst_coord, T) + IF (disca) CALL diffuse_scattering(nq, nq_super, nq_strft, nat, ntyp, amass, ityp, q_nq_zg, & + w2, z_nq_zg, q_external, zero_one_phonon, full_phonon, atmsf_a, atmsf_b, & + dim1, dim2, dim3, tau, alat, atm_zg, atom_resolved, mode_resolved, & + ntypx, at, q_in_cryst_coord, T, eps2, print_raw_data, & + nrots, kres1, kres2, kmin, kmax, col1, col2, Np) !mz_e ! ! DEALLOCATE (z, w2, dyn, dyn_blk) ! mz_b - IF (disca) DEALLOCATE (z_nq_zg,q_nq_zg) + IF (disca) DEALLOCATE (z_nq_zg, q_nq_zg) ! mz_e ! ! @@ -787,7 +790,7 @@ PROGRAM diff_sca DEALLOCATE(high_sym) ! - CALL environment_end('MATDYN') + CALL environment_end('DISCA') ! CALL mp_global_end() ! @@ -926,120 +929,10 @@ 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,fd) - !----------------------------------------------------------------------- - ! calculates the dynamical matrix at q from the (short-range part of the) - ! force constants - ! - USE kinds, ONLY : DP - USE constants, ONLY : tpi - USE io_global, ONLY : stdout - ! - IMPLICIT NONE - INTEGER nr1, nr2, nr3, nat, n1, n2, n3, nr1_, nr2_, nr3_, & - ipol, jpol, na, nb, m1, m2, m3, nint, i,j, nrws, nax - COMPLEX(DP) dyn(3,3,nat,nat), f_of_q(3,3,nat,nat) - REAL(DP) frc(nr1,nr2,nr3,3,3,nat,nat), tau(3,nat), q(3), arg, & - at(3,3), bg(3,3), r(3), weight, r_ws(3), & - total_weight, rws(0:3,nrws), alat - REAL(DP), EXTERNAL :: wsweight - REAL(DP),SAVE,ALLOCATABLE :: wscache(:,:,:,:,:) - REAL(DP), ALLOCATABLE :: ttt(:,:,:,:,:), tttx(:,:) - LOGICAL,SAVE :: first=.true. - LOGICAL :: fd - ! - nr1_=2*nr1 - nr2_=2*nr2 - nr3_=2*nr3 - FIRST_TIME : IF (first) THEN - first=.false. - ALLOCATE( wscache(-nr3_:nr3_, -nr2_:nr2_, -nr1_:nr1_, nat,nat) ) - DO na= 1, nat - DO nb= 1, nat - total_weight=0.0d0 - ! - DO n1=-nr1_,nr1_ - DO n2=-nr2_,nr2_ - DO n3=-nr3_,nr3_ - 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 - ENDDO - ENDDO - ENDDO - ENDDO - ENDIF FIRST_TIME - ! - ALLOCATE(ttt(3,nat,nr1,nr2,nr3)) - ALLOCATE(tttx(3,nat*nr1*nr2*nr3)) - ttt(:,:,:,:,:)=0.d0 - - DO na= 1, nat - DO nb= 1, nat - total_weight=0.0d0 - DO n1=-nr1_,nr1_ - DO n2=-nr2_,nr2_ - DO n3=-nr3_,nr3_ - ! - ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY SAFE RANGE! - ! - DO i= 1, 3 - r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3) - ENDDO - - weight = wscache(n3,n2,n1,nb,na) - IF (weight .GT. 0.0d0) THEN - ! - ! FIND THE VECTOR CORRESPONDING TO R IN THE ORIGINAL CELL - ! - m1 = MOD(n1+1,nr1) - IF(m1.LE.0) m1=m1+nr1 - m2 = MOD(n2+1,nr2) - IF(m2.LE.0) m2=m2+nr2 - m3 = MOD(n3+1,nr3) - IF(m3.LE.0) m3=m3+nr3 - ! WRITE(*,'(6i4)') n1,n2,n3,m1,m2,m3 - ! - ! FOURIER TRANSFORM - ! - DO i= 1,3 - ttt(i,na,m1,m2,m3)=tau(i,na)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) - ttt(i,nb,m1,m2,m3)=tau(i,nb)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) - ENDDO - - arg = tpi*(q(1)*r(1) + q(2)*r(2) + q(3)*r(3)) - DO ipol= 1, 3 - DO jpol= 1, 3 - dyn(ipol,jpol,na,nb) = & - dyn(ipol,jpol,na,nb) + & - (frc(m1,m2,m3,ipol,jpol,na,nb)+f_of_q(ipol,jpol,na,nb)) & - *CMPLX(COS(arg),-SIN(arg),kind=DP)*weight - ENDDO - ENDDO - ENDIF - total_weight=total_weight + weight - ENDDO - ENDDO - ENDDO - IF (ABS(total_weight-nr1*nr2*nr3).GT.1.0d-8) THEN - WRITE(stdout,*) total_weight - CALL errore ('frc_blk','wrong total_weight',1) - ENDIF - ENDDO - ENDDO - ! - RETURN -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,na_ifc,f_of_q,fd) + & epsil,zeu,frc,nr1,nr2,nr3,has_zstar,rws,nrws,na_ifc,f_of_q,fd) !----------------------------------------------------------------------- ! compute the dynamical matrix (the analytic part only) ! @@ -1080,7 +973,7 @@ 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,fd) - IF (has_zstar .and. .not.na_ifc) & + 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) ! LOTO 2D added celldm(1)=alat to passed arguments @@ -1967,57 +1860,10 @@ SUBROUTINE read_tau & END SUBROUTINE read_tau ! !----------------------------------------------------------------------- -SUBROUTINE gen_qpoints (ibrav, at_, bg_, nat, tau, ityp, nk1, nk2, nk3, & - nqx, nq, q, nosym) - !----------------------------------------------------------------------- - ! - USE kinds, ONLY : DP - USE cell_base, ONLY : at, bg - USE symm_base, ONLY : set_sym_bl, find_sym, s, irt, nsym, & - nrot, t_rev, time_reversal, sname - USE ktetra, ONLY : tetra_init - ! - IMPLICIT NONE - ! input - INTEGER :: ibrav, nat, nk1, nk2, nk3, ityp(*) - REAL(DP) :: at_(3,3), bg_(3,3), tau(3,nat) - LOGICAL :: nosym - ! output - INTEGER :: nqx, nq - REAL(DP) :: q(3,nqx) - ! local - REAL(DP) :: xqq(3), wk(nqx), mdum(3,nat) - LOGICAL :: magnetic_sym=.FALSE., skip_equivalence=.FALSE. - ! - time_reversal = .true. - if (nosym) time_reversal = .false. - t_rev(:) = 0 - xqq (:) =0.d0 - at = at_ - bg = bg_ - CALL set_sym_bl ( ) - ! - if (nosym) then - nrot = 1 - nsym = 1 - ENDIF - CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev, bg, nqx, & - 0,0,0, nk1,nk2,nk3, nq, q, wk) - ! - CALL find_sym ( nat, tau, ityp, .not.time_reversal, mdum ) - ! - CALL irreducible_BZ (nrot, s, nsym, time_reversal, magnetic_sym, & - at, bg, nqx, nq, q, wk, t_rev) - ! - CALL tetra_init (nsym, s, time_reversal, t_rev, at, bg, nqx, 0, 0, 0, & - nk1, nk2, nk3, nq, q) - ! - RETURN -END SUBROUTINE gen_qpoints ! !--------------------------------------------------------------------- subroutine setgam (q, gam, nat, at,bg,tau,itau_blk,nsc,alat, & - & gam_blk, nat_blk, at_blk,bg_blk,tau_blk,omega_blk, & + & gam_blk, nat_blk, at_blk,bg_blk,tau_blk,omega_blk, & & frcg, nr1,nr2,nr3, rws,nrws, fd) !----------------------------------------------------------------------- ! compute the dynamical matrix (the analytic part only) @@ -2033,14 +1879,14 @@ 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) - LOGICAL :: fd + 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 ! REAL(DP) :: arg - complex(DP) :: cfac(nat) + COMPLEX(DP) :: cfac(nat) INTEGER :: i,j,k, na,nb, na_blk, nb_blk, iq REAL(DP) :: qp(3), qbid(3,nsc) ! automatic array ! @@ -2090,130 +1936,7 @@ subroutine setgam (q, gam, nat, at,bg,tau,itau_blk,nsc,alat, & end subroutine setgam ! !-------------------------------------------------------------------- -function dos_gam (nbndx, nq, jbnd, gamma, et, ef) - !-------------------------------------------------------------------- - ! calculates weights with the tetrahedron method (Bloechl version) - ! this subroutine is based on tweights.f90 belonging to PW - ! it calculates a2F on the surface of given frequency <=> histogram - ! Band index means the frequency mode here - ! and "et" means the frequency(mode,q-point) - ! - USE kinds, ONLY: DP - USE parameters - USE ktetra, ONLY : ntetra, tetra - implicit none - ! - INTEGER :: nq, nbndx, jbnd - REAL(DP) :: et(nbndx,nq), gamma(nbndx,nq), func - - REAL(DP) :: ef - REAL(DP) :: e1, e2, e3, e4, c1, c2, c3, c4, etetra(4) - INTEGER :: ik, ibnd, nt, nk, ns, i, ik1, ik2, ik3, ik4, itetra(4) - - REAL(DP) :: f12,f13,f14,f23,f24,f34, f21,f31,f41,f42,f32,f43 - REAL(DP) :: P1,P2,P3,P4, G, o13, Y1,Y2,Y3,Y4, eps,vol, Tint - REAL(DP) :: dos_gam - - Tint = 0.0d0 - o13 = 1.0_dp/3.0_dp - eps = 1.0d-14 - vol = 1.0d0/ntetra - P1 = 0.0_dp - P2 = 0.0_dp - P3 = 0.0_dp - P4 = 0.0_dp - DO nt = 1, ntetra - ibnd = jbnd - ! - ! etetra are the energies at the vertexes of the nt-th tetrahedron - ! - DO i = 1, 4 - etetra(i) = et(ibnd, tetra(i,nt)) - ENDDO - itetra(1) = 0 - CALL hpsort (4,etetra,itetra) - ! - ! ...sort in ascending order: e1 < e2 < e3 < e4 - ! - e1 = etetra (1) - e2 = etetra (2) - e3 = etetra (3) - e4 = etetra (4) - ! - ! kp1-kp4 are the irreducible k-points corresponding to e1-e4 - ! - ik1 = tetra(itetra(1),nt) - ik2 = tetra(itetra(2),nt) - ik3 = tetra(itetra(3),nt) - ik4 = tetra(itetra(4),nt) - Y1 = gamma(ibnd,ik1)/et(ibnd,ik1) - Y2 = gamma(ibnd,ik2)/et(ibnd,ik2) - Y3 = gamma(ibnd,ik3)/et(ibnd,ik3) - Y4 = gamma(ibnd,ik4)/et(ibnd,ik4) - - IF ( e3 < ef .and. ef < e4) THEN - - f14 = (ef-e4)/(e1-e4) - f24 = (ef-e4)/(e2-e4) - f34 = (ef-e4)/(e3-e4) - - G = 3.0_dp * f14 * f24 * f34 / (e4-ef) - P1 = f14 * o13 - P2 = f24 * o13 - P3 = f34 * o13 - P4 = (3.0_dp - f14 - f24 - f34 ) * o13 - - ELSE IF ( e2 < ef .and. ef < e3 ) THEN - - f13 = (ef-e3)/(e1-e3) - f31 = 1.0_dp - f13 - f14 = (ef-e4)/(e1-e4) - f41 = 1.0_dp-f14 - f23 = (ef-e3)/(e2-e3) - f32 = 1.0_dp - f23 - f24 = (ef-e4)/(e2-e4) - f42 = 1.0_dp - f24 - - G = 3.0_dp * (f23*f31 + f32*f24) - P1 = f14 * o13 + f13*f31*f23 / G - P2 = f23 * o13 + f24*f24*f32 / G - P3 = f32 * o13 + f31*f31*f23 / G - P4 = f41 * o13 + f42*f24*f32 / G - G = G / (e4-e1) - - ELSE IF ( e1 < ef .and. ef < e2 ) THEN - - f12 = (ef-e2)/(e1-e2) - f21 = 1.0_dp - f12 - f13 = (ef-e3)/(e1-e3) - f31 = 1.0_dp - f13 - f14 = (ef-e4)/(e1-e4) - f41 = 1.0_dp - f14 - - G = 3.0_dp * f21 * f31 * f41 / (ef-e1) - P1 = o13 * (f12 + f13 + f14) - P2 = o13 * f21 - P3 = o13 * f31 - P4 = o13 * f41 - - ELSE - - G = 0.0_dp - - ENDIF - - Tint = Tint + G * (Y1*P1 + Y2*P2 + Y3*P3 + Y4*P4) * vol - - ENDDO ! ntetra - - - dos_gam = Tint !2 because DOS_ee is per 1 spin - - RETURN -end function dos_gam ! -! -!----------------------------------------------------------------------- subroutine readfg ( ifn, nr1, nr2, nr3, nat, frcg ) !----------------------------------------------------------------------- ! @@ -2300,7 +2023,7 @@ SUBROUTINE find_representations_mode_q ( nat, ntyp, xq, w2, u, tau, ityp, & ! search the symmetries only if there are no G such that Sq -> q+G ! search_sym=.TRUE. - IF ( ANY ( ft(:,1:nsymq) > 1.0d-8 ) ) THEN + IF ( ANY ( ABS(ft(:,1:nsymq)) > 1.0d-8 ) ) THEN DO isym= 1,nsymq search_sym=( search_sym.and.(ABS(gi(1,isym))<1.d-8).and. & (ABS(gi(2,isym))<1.d-8).and. & @@ -2324,64 +2047,67 @@ SUBROUTINE find_representations_mode_q ( nat, ntyp, xq, w2, u, tau, ityp, & END SUBROUTINE find_representations_mode_q !mz adds this routine -SUBROUTINE diffuse_scattering(nq,nq_super,nq_strft,nat,ntyp,amass,ityp,q,w2,z_nq_zg, & - q_external, zero_one_phonon, full_phonon, atmsf_zg_a, atmsf_zg_b, & - dim1, dim2, dim3, tau, alat, atm, atom_resolved, & - ntypx,at,q_in_cryst_coord,T) -! we start here with the WRITE_eigenvectors.f90 routine - use kinds, only: dp - use constants, only: amu_ry, ry_to_thz, ry_to_cmm1, H_PLANCK_SI, & - K_BOLTZMANN_SI, AMU_SI, pi , BOHR_RADIUS_ANGS +SUBROUTINE diffuse_scattering(nq, nq_super, nq_strft, nat, ntyp, amass, ityp, q, w2, & + z_nq_zg, q_external, zero_one_phonon, full_phonon, atmsf_a, atmsf_b, & + dim1, dim2, dim3, tau, alat, atm, atom_resolved, mode_resolved, & + ntypx, at, q_in_cryst_coord, T, eps2, print_raw_data, & + nrots, kres1, kres2, kmin, kmax, col1, col2, Np) + ! we start here with the WRITE_eigenvectors.f90 routine + USE kinds, ONLY : dp + USE constants, ONLY : amu_ry, ry_to_thz, ry_to_cmm1, H_PLANCK_SI, & + K_BOLTZMANN_SI, AMU_SI, pi, BOHR_RADIUS_ANGS USE cell_base, ONLY : bg USE io_global, ONLY : ionode, ionode_id USE mp_global, ONLY : inter_pool_comm USE mp, ONLY : mp_bcast, mp_barrier, mp_sum USE mp_world, ONLY : world_comm - implicit none + ! + IMPLICIT NONE ! input - LOGICAL, INTENT(in) :: q_in_cryst_coord, q_external - LOGICAL, INTENT(in) :: zero_one_phonon, full_phonon, atom_resolved + LOGICAL, INTENT(in) :: q_in_cryst_coord, q_external, print_raw_data + LOGICAL, INTENT(in) :: zero_one_phonon, full_phonon + LOGICAL, INTENT(in) :: atom_resolved, mode_resolved INTEGER, INTENT(in) :: dim1, dim2, dim3, T INTEGER, INTENT(in) :: nq, nat, ntyp, ntypx, nq_super, nq_strft - REAL(DP), INTENT(in) :: alat + INTEGER, INTENT(in) :: nrots, kres1, kres2, col1, col2, Np + REAL(DP), INTENT(in) :: kmin, kmax + REAL(DP), INTENT(in) :: alat, eps2 REAL(DP), INTENT(in) :: at(3, 3) - REAL(DP), intent(in) :: atmsf_zg_a(ntypx, 5), atmsf_zg_b(ntypx, 5) + REAL(DP), intent(in) :: atmsf_a(ntypx, 5), atmsf_b(ntypx, 5) REAL(DP), INTENT(in) :: q(3, nq), w2(3 * nat, nq), amass(ntyp), tau(3, nat) ! nq is the number of qpoints in sets A and B CHARACTER(LEN=3), INTENT(in) :: atm(ntypx) INTEGER ityp(nat) - complex(DP), INTENT(in) :: z_nq_zg(3*nat,3*nat,nq) + COMPLEX(DP), INTENT(in) :: z_nq_zg(3 * nat, 3 * nat, nq) ! ! local ! - INTEGER :: nat3, na, nta, ipol, i, j, k, qp, qp2, ii, p, p1, ntap,kp - INTEGER :: nq_tot + INTEGER :: nat3, na, nta, ipol, i, j, ii + INTEGER :: nq_tot, qp, qp2, p, p1, ntap, k, kp INTEGER :: ctr, ctr2, lower_bnd, upper_bnd, qlistA(nq) ! nq_tot total number of q-points (including sets A,B,C) REAL(DP) :: freq(3 * nat, nq), ph_w(3 * nat, nq),l_q(3 * nat, nq) REAL(DP) :: q_A(3), q_B(3), atm_ffct(nat, nq), p_q(3 * nat, nq), e_nq(3 * nat, nq) - REAL(DP) :: hbar, ang, dotp, dotp2, dotp3, dotp4, ph_w_mean - REAL(DP), PARAMETER :: eps = 1.0d-4 !, + REAL(DP) :: hbar, ang, dotp, dotp2, dotp3, dotp4, ph_w_mean, units + REAL(DP), PARAMETER :: eps = 1.0d-4 ! l_q is the amplitude \sigma at temperature T, e_nq --> to calculate total vibrational energy ! p_q is the momentum on the nuclei \hbar\2\l_\bq\nu \SQRT(n_{q\nu,T}+1/2) ! - complex(DP) :: imagi - complex(DP) :: z_zg(3 * nat,3 * nat, nq) - CHARACTER(len=256) :: filename, pointer_mz + COMPLEX(DP) :: imagi + COMPLEX(DP) :: z_zg(3 * nat,3 * nat, nq) + CHARACTER(len=256) :: filename, pt_mz, pt_mz2 ! ALLOCATE TABLES - !COMPLEX(DP), ALLOCATABLE :: DW_T_fact(:,:), strctr_fact(:,:,:,:), strctr_fact_q1(:,:,:) - !COMPLEX(DP), ALLOCATABLE :: DW_T_fact_q0(:,:), strctr_fact_q0(:,:) REAL(DP), ALLOCATABLE :: DW_T_fact_kkp(:, :), DW_T_fact_q0_kkp(:, :), DW_T_fact_q0(:) REAL(DP), ALLOCATABLE :: sigma_DW(:, :, :), Q_ppkkaa(:, :, :) - REAL(DP), ALLOCATABLE :: strctr_fact_kkp(:, :, :), strctr_fact(:), strctr_fact_j(:, :)!, strctr_fact_full(:) - REAL(DP), ALLOCATABLE :: strctr_fact_q0(:), strctr_fact_q0_kkp(:, :, :) - Complex(DP), ALLOCATABLE :: strctr_fact_full(:), strctr_fact_full_kk(:, :, :) + REAL(DP), ALLOCATABLE :: strf_kkp(:, :, :), strf(:), strf_j(:, :) + REAL(DP), ALLOCATABLE :: strf_q0(:), strf_q0_kkp(:, :, :), strf_rot(:,:) + COMPLEX(DP), ALLOCATABLE :: strf_full(:), strf_full_kk(:, :, :) ! for displacements - ! matrices to account for the coupling terms between different phonon branches ! REAL(DP), ALLOCATABLE :: Rlist(:, :) ! ! hbar = 0.5 * H_PLANCK_SI / pi ! reduce Plnack constant + units = DBLE(2.d0 * pi / alat / BOHR_RADIUS_ANGS) ang = 1.0E-10 ! angstrom units imagi = (0.0d0, 1.0d0) !imaginary unit ! Inititialize eigenvectors matrix @@ -2391,7 +2117,7 @@ SUBROUTINE diffuse_scattering(nq,nq_super,nq_strft,nat,ntyp,amass,ityp,q,w2,z_nq nat3 = 3 * nat ! ! - freq = 0.0 + freq = 0.d0 CALL fkbounds( nq, lower_bnd, upper_bnd ) DO qp = lower_bnd, upper_bnd ! 1, nq ! convert eigenvectors to mass-unsCALLed @@ -2412,6 +2138,15 @@ SUBROUTINE diffuse_scattering(nq,nq_super,nq_strft,nat,ntyp,amass,ityp,q,w2,z_nq CALL mp_sum(freq, inter_pool_comm) CALL mp_sum(z_zg, inter_pool_comm) CALL mp_barrier(inter_pool_comm) + ! + IF (ionode) THEN + WRITE(*,*) nq, nat3 + DO i = 1, nat3 + DO qp = 1, nq + IF (w2(i, qp) < 0.0 .AND. ionode) WRITE(*,*) "WARNING: Negative frequencies", w2(i, qp) + ENDDO + ENDDO + ENDIF ! Einstein model set polarization vectors to unity ! ! z_zg = (1.d0, 1.d0)/DBLE(SQRT(real(2.d0*nat3))) ! For Einstein model @@ -2447,7 +2182,7 @@ SUBROUTINE diffuse_scattering(nq,nq_super,nq_strft,nat,ntyp,amass,ityp,q,w2,z_nq CALL fkbounds( nq, lower_bnd, upper_bnd ) DO qp = lower_bnd, upper_bnd ! 1, nq DO i = 1, nat3 - IF (w2(i, qp) .lt. 0.0d0) THEN + IF (w2(i, qp) .lt. 0.0d0 + eps2) THEN l_q(i, qp) = 0.d0 p_q(i, qp) = 0.d0 ELSE @@ -2457,354 +2192,398 @@ SUBROUTINE diffuse_scattering(nq,nq_super,nq_strft,nat,ntyp,amass,ityp,q,w2,z_nq ! l_qn SQRT(2n_qn +1) = sigma_qn ENDDO ! First set all displ. amplitudes to zero for accoustic modes for every G - if ( (( ABS(MODULO(ABS(q(1, qp)),1.0d0)) .LT. eps) .OR. ( abs(MODULO(abs(q(1, qp)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q(2, qp)),1.0d0)) .LT. eps) .OR. ( abs(MODULO(abs(q(2, qp)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q(3, qp)),1.0d0)) .LT. eps) .OR. ( abs(MODULO(abs(q(3, qp)),1.0d0)-1.0d0) .LT. eps)) ) THEN - l_q(1, qp) = 0.0d0 - l_q(2, qp) = 0.0d0 - l_q(3, qp) = 0.0d0 - p_q(1, qp) = 0.0d0 - p_q(2, qp) = 0.0d0 - p_q(3, qp) = 0.0d0 + IF ( (( ABS(MODULO(ABS(q(1, qp)),1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q(1, qp)),1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q(2, qp)),1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q(2, qp)),1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q(3, qp)),1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q(3, qp)),1.0d0)-1.0d0) .LT. eps)) ) THEN + l_q(1, qp) = 0.0d0 + l_q(2, qp) = 0.0d0 + l_q(3, qp) = 0.0d0 + p_q(1, qp) = 0.0d0 + p_q(2, qp) = 0.0d0 + p_q(3, qp) = 0.0d0 ENDIF ENDDO ! qp -CALL mp_sum(l_q, inter_pool_comm) -CALL mp_sum(p_q, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) -!IF (ionode) WRITE(*,*) "step2" -! -! To determine which points belong in set A -qlistA = 0 -CALL fkbounds( nq, lower_bnd, upper_bnd ) - DO qp = lower_bnd, upper_bnd ! 1, nq - q_A(1) = q(1, qp) + q(1, qp) - q_A(2) = q(2, qp) + q(2, qp) - q_A(3) = q(3, qp) + q(3, qp) + CALL mp_sum(l_q, inter_pool_comm) + CALL mp_sum(p_q, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) ! - ! q_A = q(:, qp) + q(:, qp) - IF ( (( ABS(MODULO(ABS(q_A(1)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(1)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(2)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(2)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(3)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(3)),1.0d0)-1.0d0) .LT. eps)) ) THEN - qlistA(qp) = 1 - !WRITE(*,*) "qAAA", q(:, qp) - ENDIF - ENDDO -CALL mp_sum(qlistA, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) -! -! -! Exponent of DW factors -ALLOCATE(sigma_DW(nat,3,3)) -! sigma_DW thermal displacement tensor -CALL fkbounds( nq - nq_strft, lower_bnd, upper_bnd ) -!WRITE(*,*) "mz_para", lower_bnd, upper_bnd, nq - nq_strft + ! To determine which points belong in set A + qlistA = 0 + CALL fkbounds( nq, lower_bnd, upper_bnd ) + DO qp = lower_bnd, upper_bnd ! 1, nq + q_A(1) = q(1, qp) + q(1, qp) + q_A(2) = q(2, qp) + q(2, qp) + q_A(3) = q(3, qp) + q(3, qp) + ! + ! q_A = q(:, qp) + q(:, qp) + IF ( (( ABS(MODULO(ABS(q_A(1)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(1)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(2)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(2)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(3)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(3)), 1.0d0) - 1.0d0) .LT. eps)) ) THEN + qlistA(qp) = 1 + !WRITE(*,*) "qAAA", q(:, qp) + ENDIF + ENDDO + CALL mp_sum(qlistA, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) + ! + ! + ! Exponent of DW factors + ALLOCATE(sigma_DW(nat, 3, 3)) + ! sigma_DW: thermal displacement tensor + CALL fkbounds( nq - nq_strft, lower_bnd, upper_bnd ) + ! sigma_DW = 0.0d0 DO qp = lower_bnd, upper_bnd ! for q-points in supercell ! DO qp = 1, nq_super DO k= 1, nat nta = ityp(k) - DO i= 1,3 ! i is for cart directions - DO p = 1,3 ! p is also for cartesian - DO j = 1,nat3 ! modes \nu + DO i= 1, 3 ! i is for cart directions + DO p = 1, 3 ! p is also for cartesian + DO j = 1, nat3 ! modes \nu IF (qlistA(qp) .EQ. 1) THEN - sigma_DW(k,i,p) = sigma_DW(k,i,p) + & - 1.0/DBLE(nq_tot)/DBLE(2.0d0*amass(nta))*DBLE(z_zg((k-1)*3+i,j,qp) & - *CONJG(z_zg((k-1)*3+p,j,qp)))*l_q(j,qp)**2 + sigma_DW(k, i, p) = sigma_DW(k, i, p) + 1.0 / DBLE(nq_tot) / DBLE(2.0d0 * amass(nta)) * & + DBLE(z_zg((k - 1) * 3 + i, j, qp) * & + CONJG(z_zg((k - 1) * 3 + p, j, qp))) * l_q(j, qp)**2 ELSE - sigma_DW(k,i,p) = sigma_DW(k,i,p) + & - 1.0/DBLE(nq_tot)/DBLE(amass(nta))*DBLE(z_zg((k-1)*3+i,j,qp)*CONJG(z_zg((k-1)*3+p,j,qp)))*l_q(j,qp)**2 - ! an extra factor of 2 is need because we have only points in set B + sigma_DW(k, i, p) = sigma_DW(k, i, p) + 1.0/DBLE(nq_tot) / DBLE(amass(nta)) * & + DBLE(z_zg((k - 1) * 3 + i, j, qp) * & + CONJG(z_zg((k - 1) * 3 + p, j, qp))) * l_q(j, qp)**2 + ! an extra factor of 2 is need because we have only points in set B ENDIF ENDDO ! j ENDDO ! p ENDDO ! i ENDDO ! k ENDDO ! qp -! -CALL mp_sum(sigma_DW, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) -! -IF (ionode) THEN - WRITE(*,*) "DW exponent" - DO k= 1, nat - WRITE(*,*) "atom:", k - DO i= 1,3 ! i is for cart directions - WRITE(*,'(3F14.8)') (sigma_DW(k,i,p), p = 1, 3) - ENDDO ! i - ENDDO ! k -ENDIF -!!!!!!!!!!!!! -IF (ionode) THEN - atm_ffct = 0.d0 - DO k = 1, nat - nta = ityp(k) - DO qp = nq_super + 1, nq !1, qpts_strf - DO ii = 1, 5 - q_A = q(:, qp) !* ( 2.d0 * pi / alat / BOHR_RADIUS_ANGS) - CALL cryst_to_cart(1, q_A, bg, +1) - q_A = q_A * ( 2.d0 * pi / alat / BOHR_RADIUS_ANGS) - atm_ffct(k, qp) = atm_ffct(k, qp) + (atmsf_zg_a(nta, ii) * EXP(-atmsf_zg_b(nta, ii) * & - (NORM2(q_A) / 4.d0 / pi)**2)) - ENDDO - ENDDO - ENDDO -ENDIF !(ionode) -CALL mp_bcast(atm_ffct, ionode_id, world_comm) -! -!! atm_ffct = 1 - ! convert back to crystal -!!!!!!!!!!!!! -IF (zero_one_phonon) THEN - ALLOCATE(DW_T_fact_q0_kkp(nq, nat),DW_T_fact_q0(nq)) - ALLOCATE(strctr_fact_q0(nq)) - IF (atom_resolved) ALLOCATE(strctr_fact_q0_kkp(nq, nat, nat)) - !CALL cryst_to_cart(nq,q,at,1) - DW_T_fact_q0 = 0.0d0!(0.0d0, 0.0d0) - DW_T_fact_q0_kkp = 0.0d0!(0.0d0, 0.0d0) - strctr_fact_q0 = 0.0d0!(0.0d0, 0.0d0) - ! - CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) - DO qp = lower_bnd + nq_super, upper_bnd + nq_super - q_A(1) = q(1, qp) !+ q(1,qp) - q_A(2) = q(2, qp) !+ q(2,qp) - q_A(3) = q(3, qp) !+ q(3,qp) - - ! q_A = q(:,qp) + q(:,qp) - if ( (( ABS(MODULO(ABS(q_A(1)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(1)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(2)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(2)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(3)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(3)),1.0d0)-1.0d0) .LT. eps)) ) THEN - !WRITE(*,*) "q_mz", q(:,qp), l_q(:,qp) - ! - CALL cryst_to_cart(1, q(:,qp), bg, +1) - DO k= 1 , nat ! k represents the atom - nta = ityp(k) - dotp = 0.0d0 - DO i= 1,3 ! i is for cart directions - DO p = 1,3 ! p is also for cartesian - dotp = dotp + q(i, qp) * q(p, qp) * sigma_DW(k, i, p) * (2 * pi / alat / BOHR_RADIUS_ANGS)**2 - ! this will be complex otherwise - ENDDO ! p - ENDDO ! i - DW_T_fact_q0_kkp(qp, k) = DW_T_fact_q0_kkp(qp, k) + dotp ! DW factor for each q - !WRITE(*,*) "mzole", dotp, DW_T_fact_q0_kkp(qp,k) - ENDDO ! k - CALL cryst_to_cart(1, q(:,qp), at, -1) - ENDIF - ENDDO ! qp -CALL mp_sum(DW_T_fact_q0_kkp, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) ! - ! evaluate structure factor for zero-phonon contribution -DO qp = lower_bnd + nq_super, upper_bnd + nq_super - q_A(1) = q(1,qp) !+ q(1,qp) - q_A(2) = q(2,qp) !+ q(2,qp) - q_A(3) = q(3,qp) !+ q(3,qp) - !! q_A = q(:,qp) + q(:,qp) - if ( (( ABS(MODULO(ABS(q_A(1)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(1)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(2)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(2)),1.0d0)-1.0d0) .LT. eps)) .AND. & - (( ABS(MODULO(ABS(q_A(3)),1.0d0)) .LT. eps) .OR. & - ( ABS(MODULO(ABS(q_A(3)),1.0d0)-1.0d0) .LT. eps)) ) THEN - CALL cryst_to_cart(1,q(:,qp),bg,+1) - DO k = 1 , nat ! k represents the atom - DO kp = 1 , nat ! k represents the atom - dotp = 0.0d0 - DO ii = 1, 3 - dotp = dotp + q(ii,qp) * (tau(ii,k) - tau(ii,kp))! - ENDDO - !WRITE(*,*) "mz_ole2", q(:,qp), tau(:,k), cos(2*pi*dotp) - IF (atom_resolved) THEN - strctr_fact_q0_kkp(qp,k,kp) = EXP(-DW_T_fact_q0_kkp(qp, k))*EXP(-DW_T_fact_q0_kkp(qp, kp))& - !*nq_tot*nq_tot*EXP(imagi*2*pi*dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) - * nq_tot * nq_tot * COS(2*pi*dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) - ENDIF - strctr_fact_q0(qp) = strctr_fact_q0(qp) + EXP(-DW_T_fact_q0_kkp(qp,k))*EXP(-DW_T_fact_q0_kkp(qp,kp))& - * nq_tot * nq_tot * COS(2*pi*dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) !+ strctr_fact_q0_kkp(qp,k,kp) - ENDDO ! kp - ENDDO ! k - CALL cryst_to_cart(1, q(:,qp), at, -1) - ENDIF -ENDDO ! qp -IF (atom_resolved) CALL mp_sum(strctr_fact_q0_kkp, inter_pool_comm) -CALL mp_sum(strctr_fact_q0, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) -ENDIF ! if zero_one_phonon -! -! convert to Cartesian -CALL cryst_to_cart(nq, q, bg, +1) -! now for structure factor one and full phonon contribution -ALLOCATE(DW_T_fact_kkp(nq, nat)) -DW_T_fact_kkp = 0.0d0!(0.0d0, 0.0d0) -CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) -DO qp = lower_bnd + nq_super, upper_bnd + nq_super - DO k = 1, nat ! k represents the atom - dotp = 0.0d0 - !CALL cryst_to_cart(1,q(:,qp),bg,+1) ! unne - DO i= 1,3 ! i is for cart directions - DO p = 1,3 ! p is also for cartesian - dotp = dotp + q(i, qp) * q(p, qp) * sigma_DW(k, i, p) * (2.0d0 * pi / alat / BOHR_RADIUS_ANGS)**2 - ctr2 = ctr2 + 1 - ! this will be complex otherwise - ENDDO ! p - ENDDO ! i - DW_T_fact_kkp(qp, k) = DW_T_fact_kkp(qp, k) + dotp - !CALL cryst_to_cart(1,q(:,qp),at,-1) ! unne - ENDDO ! k - !WRITE(*,*) "check_point", DW_T_fact_kkp(qp,k) , EXP(-DW_T_fact_kkp(qp,k)) -ENDDO ! qp -CALL mp_sum(DW_T_fact_kkp, inter_pool_comm) -CALL mp_barrier(inter_pool_comm) -DEALLOCATE(sigma_DW) -! -! compute one_phonon contribution -IF (zero_one_phonon) THEN -! Here q_pts should be in cartessian - ALLOCATE(strctr_fact_j(nq,nat3)) - IF (atom_resolved) ALLOCATE(strctr_fact_kkp(nq, nat, nat)) - !WRITE(*,*) "mz6" - IF (atom_resolved) strctr_fact_kkp = 0.0d0 - strctr_fact_j = 0.0d0 - - CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) - DO qp = lower_bnd + nq_super, upper_bnd + nq_super - DO k = 1 , nat ! k represents the atom - nta = ityp(k) - DO kp = 1 , nat ! k represents the atom - ntap = ityp(kp) - dotp2 = 0.0d0 - DO ii= 1,3 - dotp2 = dotp2 + q(ii,qp)*(tau(ii,k)-tau(ii,kp))! - ENDDO - DO j = 1, nat3 ! modes \nu - dotp = 0.0d0 - DO i = 1, 3 ! i is for cart directions - DO p = 1, 3 ! p is also for cartesian - dotp = dotp + q(i, qp) * q(p, qp) * (2.0d0 * pi / alat / BOHR_RADIUS_ANGS)**2 * & ! - DBLE(z_zg((k-1)*3+i,j,qp)*CONJG(z_zg((kp-1)*3+p,j,qp)) * EXP(-imagi * 2.0d0 * pi * dotp2) ) - - ENDDO ! p - ENDDO ! i - ! I do not really need to store this matrix - IF (atom_resolved) THEN - strctr_fact_kkp(qp, k, kp) = strctr_fact_kkp(qp, k, kp) & - + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) & - * nq_tot / DBLE(SQRT(amass(nta) * amass(ntap))) * dotp * l_q(j, qp)**2 & - * atm_ffct(k, qp) * atm_ffct(kp, qp) - ENDIF - ! IF (qlistA(qp) .EQ. 1) THEN - ! strctr_fact_j(qp, j) = strctr_fact_j(qp, j) + 2.0d0 * strctr_fact_kkp(qp, k, kp) - ! ELSE - strctr_fact_j(qp, j) = strctr_fact_j(qp, j) + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) & - * nq_tot / DBLE(SQRT(amass(nta) * amass(ntap))) * dotp * l_q(j, qp)**2 & - * atm_ffct(k, qp) * atm_ffct(kp, qp) - ENDDO ! j - ENDDO ! kp - ENDDO ! k - ENDDO ! qp - IF (atom_resolved) CALL mp_sum(strctr_fact_kkp, inter_pool_comm) - CALL mp_sum(strctr_fact_j, inter_pool_comm) + CALL mp_sum(sigma_DW, inter_pool_comm) CALL mp_barrier(inter_pool_comm) -ENDIF ! zero_one_phonon -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -IF (zero_one_phonon) THEN - ALLOCATE(strctr_fact(nq)) - strctr_fact = 0.0d0 - DO qp = nq_super + 1, nq - DO j = 1,nat3 - strctr_fact(qp) = strctr_fact(qp) + strctr_fact_j(qp,j) - ! - strctr_fact_j(qp,j) = strctr_fact_j(qp,j) + strctr_fact_q0(qp) - ENDDO - strctr_fact(qp) = strctr_fact(qp) + strctr_fact_q0(qp) - IF (atom_resolved) strctr_fact_kkp(qp, :, :) = strctr_fact_kkp(qp, :, :) + strctr_fact_q0_kkp(qp, :, :) - ENDDO - - IF (ionode) THEN - filename = 'Bragg_scattering.dat' - OPEN (unit = 90, file = filename, status = 'unknown', form = 'formatted') - filename = 'structure_factor_q_nu_one-phonon.dat' - OPEN (unit = 95, file = filename, status = 'unknown', form = 'formatted') - filename = 'structure_factor_one-phonon.dat' - OPEN (unit = 97, file = filename, status = 'unknown', form = 'formatted') - ! - DO qp = nq_super + 1, nq - WRITE (90, '(40f26.6)') q(:,qp)*(2*pi/alat/BOHR_RADIUS_ANGS), strctr_fact_q0(qp) - WRITE (95, '(40f26.6)') q(:,qp)*(2*pi/alat/BOHR_RADIUS_ANGS), ( strctr_fact_j(qp, j), j = 1, nat3) - WRITE (97, '(40f26.6)') q(:,qp)*(2*pi/alat/BOHR_RADIUS_ANGS), strctr_fact(qp) - ! - ENDDO ! - CLOSE(90) - CLOSE(95) - CLOSE(97) - CLOSE(98) - CLOSE(99) + IF (ionode) THEN + WRITE(*,*) "Mean-squared displacement (Ang**2)" + DO k= 1, nat + WRITE(*,*) "atom:", k + DO i= 1,3 ! i is for cart directions + WRITE(*,'(3F14.8)') (2.0d0 * sigma_DW(k, i, p), p = 1, 3) + ENDDO ! i + ENDDO ! k + ENDIF + !!!!!!!!!!!!! + IF (ionode) THEN + atm_ffct = 0.d0 + DO k = 1, nat + nta = ityp(k) + DO qp = nq_super + 1, nq !1, qpts_strf + DO ii = 1, 5 + q_A = q(:, qp) + CALL cryst_to_cart(1, q_A, bg, +1) + q_A = q_A * units + atm_ffct(k, qp) = atm_ffct(k, qp) + (atmsf_a(nta, ii) * & + EXP(-atmsf_b(nta, ii) * (NORM2(q_A) / 4.d0 / pi)**2)) + ENDDO + ENDDO + ENDDO + ENDIF !(ionode) + CALL mp_bcast(atm_ffct, ionode_id, world_comm) + ! + IF (zero_one_phonon) THEN + ALLOCATE(DW_T_fact_q0_kkp(nq, nat),DW_T_fact_q0(nq)) + ALLOCATE(strf_q0(nq)) + IF (atom_resolved) ALLOCATE(strf_q0_kkp(nq, nat, nat)) + DW_T_fact_q0_kkp = 0.0d0 + DW_T_fact_q0 = 0.0d0 + strf_q0 = 0.0d0 + ! + CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) + DO qp = lower_bnd + nq_super, upper_bnd + nq_super + q_A(1) = q(1, qp) + q_A(2) = q(2, qp) + q_A(3) = q(3, qp) + ! + IF ( (( ABS(MODULO(ABS(q_A(1)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(1)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(2)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(2)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(3)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(3)), 1.0d0) - 1.0d0) .LT. eps)) ) THEN + ! + CALL cryst_to_cart(1, q(:, qp), bg, +1) + DO k= 1 , nat ! k represents the atom + nta = ityp(k) + dotp = 0.0d0 + DO i= 1, 3 ! i is for cart directions + DO p = 1, 3 ! p is also for cartesian + dotp = dotp + q(i, qp) * q(p, qp) * sigma_DW(k, i, p) * units**2 + ! this will be COMPLEX otherwise + ENDDO ! p + ENDDO ! i + DW_T_fact_q0_kkp(qp, k) = DW_T_fact_q0_kkp(qp, k) + dotp ! DW factor for each q + ENDDO ! k + CALL cryst_to_cart(1, q(:,qp), at, -1) + ENDIF + ENDDO ! qp + CALL mp_sum(DW_T_fact_q0_kkp, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) + ! + ! evaluate structure factor for zero-phonon contribution + DO qp = lower_bnd + nq_super, upper_bnd + nq_super + q_A(1) = q(1, qp) + q_A(2) = q(2, qp) + q_A(3) = q(3, qp) + ! + IF ( (( ABS(MODULO(ABS(q_A(1)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(1)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(2)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(2)), 1.0d0) - 1.0d0) .LT. eps)) .AND. & + (( ABS(MODULO(ABS(q_A(3)), 1.0d0)) .LT. eps) .OR. & + ( ABS(MODULO(ABS(q_A(3)), 1.0d0) - 1.0d0) .LT. eps)) ) THEN + CALL cryst_to_cart(1, q(:, qp), bg, +1) + DO k = 1, nat ! k for the atom + DO kp = 1, nat ! kp for the atom + dotp = 0.0d0 + DO ii = 1, 3 + dotp = dotp + q(ii, qp) * (tau(ii, k) - tau(ii, kp))! + ENDDO + IF (atom_resolved) THEN + strf_q0_kkp(qp, k, kp) = EXP(-DW_T_fact_q0_kkp(qp, k))*EXP(-DW_T_fact_q0_kkp(qp, kp))& + * nq_tot * nq_tot * COS(2.d0 * pi * dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) + !*nq_tot*nq_tot*EXP(imagi*2*pi*dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) + ENDIF + strf_q0(qp) = strf_q0(qp) + EXP(-DW_T_fact_q0_kkp(qp,k))*EXP(-DW_T_fact_q0_kkp(qp,kp))& + * nq_tot * nq_tot * COS(2.d0 * pi * dotp) * atm_ffct(k, qp) * atm_ffct(kp, qp) + ENDDO ! kp + ENDDO ! k + CALL cryst_to_cart(1, q(:,qp), at, -1) + ENDIF + ENDDO ! qp + ! + IF (atom_resolved) CALL mp_sum(strf_q0_kkp, inter_pool_comm) + ! + CALL mp_sum(strf_q0, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) + ENDIF ! if zero_one_phonon + ! + ! convert to Cartesian + CALL cryst_to_cart(nq, q, bg, +1) + ! now for structure factor one and full phonon contribution + ALLOCATE(DW_T_fact_kkp(nq, nat)) + DW_T_fact_kkp = 0.0d0!(0.0d0, 0.0d0) + CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) + DO qp = lower_bnd + nq_super, upper_bnd + nq_super + DO k = 1, nat ! k represents the atom + dotp = 0.0d0 + DO i= 1, 3 ! i is for Cart. directions + DO p = 1, 3 ! p is also for Cart. directions + dotp = dotp + q(i, qp) * q(p, qp) * sigma_DW(k, i, p) * units**2 + ctr2 = ctr2 + 1 + ENDDO ! p + ENDDO ! i + DW_T_fact_kkp(qp, k) = DW_T_fact_kkp(qp, k) + dotp + ENDDO ! k + ENDDO ! qp + ! + CALL mp_sum(DW_T_fact_kkp, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) + ! + DEALLOCATE(sigma_DW) + ! + ! compute one_phonon contribution + IF (zero_one_phonon) THEN + ! Here q_pts should be in cartessian + ALLOCATE(strf_j(nq,nat3)) + IF (atom_resolved) ALLOCATE(strf_kkp(nq, nat, nat)) + !WRITE(*,*) "mz6" + IF (atom_resolved) strf_kkp = 0.0d0 + strf_j = 0.0d0 + + CALL fkbounds( nq - nq_super, lower_bnd, upper_bnd ) + DO qp = lower_bnd + nq_super, upper_bnd + nq_super + DO k = 1, nat ! k represents the atom + nta = ityp(k) + DO kp = 1, nat ! k represents the atom + ntap = ityp(kp) + dotp2 = 0.0d0 + DO ii = 1, 3 + dotp2 = dotp2 + q(ii, qp) * (tau(ii, k) - tau(ii, kp))! + ENDDO + DO j = 1, nat3 ! modes \nu + dotp = 0.0d0 + DO i = 1, 3 ! i is for cart directions + DO p = 1, 3 ! p is also for cartesian + dotp = dotp + q(i, qp) * q(p, qp) * units**2 * & ! + DBLE(z_zg((k - 1) * 3 + i, j, qp) * CONJG(z_zg((kp - 1) * 3 + p, j, qp)) * & + EXP(-imagi * 2.0d0 * pi * dotp2)) + + ENDDO ! p + ENDDO ! i + IF (atom_resolved) THEN + strf_kkp(qp, k, kp) = strf_kkp(qp, k, kp) + & + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) * & + nq_tot / DBLE(SQRT(amass(nta) * amass(ntap))) * dotp * & + l_q(j, qp)**2 * atm_ffct(k, qp) * atm_ffct(kp, qp) + ENDIF + ! + strf_j(qp, j) = strf_j(qp, j) + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) * & + nq_tot / DBLE(SQRT(amass(nta) * amass(ntap))) * dotp * & + l_q(j, qp)**2 * atm_ffct(k, qp) * atm_ffct(kp, qp) + ENDDO ! j + ENDDO ! kp + ENDDO ! k + ENDDO ! qp + IF (atom_resolved) CALL mp_sum(strf_kkp, inter_pool_comm) + CALL mp_sum(strf_j, inter_pool_comm) + CALL mp_barrier(inter_pool_comm) + ENDIF ! zero_one_phonon + ! + ! + IF (zero_one_phonon) THEN + ALLOCATE(strf(nq)) + strf = 0.0d0 + DO qp = nq_super + 1, nq + DO j = 1,nat3 + strf(qp) = strf(qp) + strf_j(qp, j) + ! + strf_j(qp, j) = strf_j(qp, j) + strf_q0(qp) + ENDDO + strf(qp) = strf(qp) + strf_q0(qp) + IF (atom_resolved) strf_kkp(qp, :, :) = strf_kkp(qp, :, :) + strf_q0_kkp(qp, :, :) + ENDDO + ! APPLY BROADENING and print outputs + ctr = 0 + DO p = nq_super + 1, nq + IF ((q(col1, p) .GT. 0.d0 - eps) .AND. (q(col2, p) .GT. 0.d0 - eps)) THEN + ctr = ctr + 1 + ENDIF + ENDDO + ALLOCATE(strf_rot(ctr * nrots,4)) + ! + strf_rot = 0.d0 + IF (ionode) CALL rotate(strf, q, nq, nq_super, 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_one-ph_broad.dat') + ! for mode_resolved + IF (mode_resolved) THEN + DO j = 1, nat3 + strf_rot = 0.d0 + WRITE( pt_mz,'(i2.2)') j + filename = 'strf_mode_'//TRIM(pt_mz)//'_broad.dat' + IF (ionode) CALL rotate(strf_j(:, j), q, nq, nq_super, 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, filename) + ENDDO + ENDIF ! mode_resolved ! IF (atom_resolved) THEN DO k = 1, nat - WRITE( pointer_mz,'(i1.1)') k - filename = 'structure_factor_q_zero_one_k_resolved_' // TRIM( pointer_mz ) //'.dat' !'.fp' - OPEN (unit = 66, file = filename, status = 'unknown', form = 'formatted') - DO qp = nq_super + 1, nq - WRITE (66, '(40f26.6)') q(:, qp) * (2*pi/alat/BOHR_RADIUS_ANGS), & - ( DBLE((strctr_fact_kkp(qp, k, kp) + strctr_fact_kkp(qp, kp, k))/2.0d0), kp = 1, k) - ! NOTE: We should take the real part if and only if we sum - ! [strctr_fact_full_kk(qp, 2, 1) + strctr_fact_full_kk(qp, 1, 2)]/2. - ! strctr_fact_full_kk(qp, 2, 1) is not necessarilly real - !, aimag(strctr_fact_full(qp)) - ENDDO ! qp - CLOSE(66) + DO kp = 1, k + strf_rot = 0.d0 + WRITE( pt_mz,'(i1.1)') k + WRITE( pt_mz2,'(i1.1)') kp + filename = 'strf_one_k_' // TRIM(pt_mz) // '_' & + // TRIM(pt_mz2) // '_broad.dat' + ! + IF (ionode) CALL rotate(strf_kkp(:, k, kp), q, nq, nq_super, 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, filename) + !'zero_one_k_'//TRIM()//'_'//TRIM(pt_mz2)//'_broad.dat') + ENDDO ! kp ENDDO - ENDIF ! k-resolved + ENDIF + ! + DEALLOCATE(strf_rot) + ! + IF (print_raw_data) THEN + IF (ionode) THEN + filename = 'Bragg_scattering.dat' + OPEN (unit = 90, file = filename, status = 'unknown', form = 'formatted') + filename = 'strf_q_nu_one-phonon.dat' + OPEN (unit = 95, file = filename, status = 'unknown', form = 'formatted') + filename = 'strf_one-phonon.dat' + OPEN (unit = 97, file = filename, status = 'unknown', form = 'formatted') + ! + DO qp = nq_super + 1, nq + WRITE (90, '(40f26.6)') q(:,qp) * units, strf_q0(qp) + WRITE (95, '(40f26.6)') q(:,qp) * units, ( strf_j(qp, j), j = 1, nat3) + WRITE (97, '(40f26.6)') q(:,qp) * units, strf(qp) + ! + ENDDO + ! + CLOSE(90) + CLOSE(95) + CLOSE(97) + ! + IF (atom_resolved) THEN + DO k = 1, nat + WRITE( pt_mz,'(i1.1)') k + filename = 'strf_one_k_resolved_' // TRIM( pt_mz ) //'.dat' !'.fp' + OPEN (unit = 66, file = filename, status = 'unknown', form = 'formatted') + DO qp = nq_super + 1, nq + WRITE (66, '(40f26.6)') q(:, qp) * units, & + ( DBLE((strf_kkp(qp, k, kp) + strf_kkp(qp, kp, k)) / 2.0d0), kp = 1, k) + ! NOTE: We should take the real part if and only if we sum + ! [strf_full_kk(qp, 2, 1) + strf_full_kk(qp, 1, 2)]/2. + ! strf_full_kk(qp, 2, 1) is not necessarilly real + ENDDO ! qp + CLOSE(66) + ENDDO + ENDIF ! k-resolved + ! + ENDIF + ENDIF ! print_raw_data + ENDIF ! zero_one_phonon + ! + ! + IF (zero_one_phonon) THEN + DEALLOCATE(DW_T_fact_q0_kkp, DW_T_fact_q0) + DEALLOCATE(strf_j, strf,strf_q0) + IF (atom_resolved) DEALLOCATE(strf_kkp, strf_q0_kkp) + ENDIF + ! Compute Full strf + ! + IF (full_phonon) THEN + ! Generate lattice vectors in crystal coordinates + ALLOCATE(Rlist(3, nq_tot)) ! - ENDIF -ENDIF ! zero_one_phonon -! -! -IF (zero_one_phonon) THEN - DEALLOCATE(DW_T_fact_q0_kkp, DW_T_fact_q0) - DEALLOCATE(strctr_fact_j, strctr_fact,strctr_fact_q0) - IF (atom_resolved) DEALLOCATE(strctr_fact_kkp, strctr_fact_q0_kkp) -ENDIF -! Compute Full structure_factor -! -IF (full_phonon) THEN - ! Generate lattice vectors in crystal coordinates - ALLOCATE(Rlist(3, nq_tot)) - ! - ctr2 = 1 - DO i = 0, dim3 - 1 - DO j = 0, dim2 - 1 - DO k= 0, dim1 - 1 - Rlist(1, ctr2) = k - Rlist(2, ctr2) = j - Rlist(3, ctr2) = i !(/ k,j,i /) - !WRITE(*,*) Rlist(ctr2,:) - ctr2= ctr2 + 1 + ctr2 = 1 + DO i = 0, dim3 - 1 + DO j = 0, dim2 - 1 + DO k= 0, dim1 - 1 + Rlist(1, ctr2) = k + Rlist(2, ctr2) = j + Rlist(3, ctr2) = i + ! + ctr2= ctr2 + 1 + ENDDO + ENDDO ENDDO - ENDDO - ENDDO - ! - ! Convert to Cartesian - CALL cryst_to_cart(nq_tot, Rlist, at, +1) - ! - ! Convert to crystal - ALLOCATE(strctr_fact_full(nq)) - IF (atom_resolved) ALLOCATE(strctr_fact_full_kk(nq, nat, nat)) - ! - CALL fkbounds( nq_tot, lower_bnd, upper_bnd ) - ! - ALLOCATE(Q_ppkkaa(upper_bnd - lower_bnd + 1, nat3, nat3)) - ! allocate array in parallel / minimize memory burden - Q_ppkkaa = 0.0d0 - ! DO qp2 = lower_bnd, upper_bnd ! for q-points in supercell + ! + ! Convert to Cartesian + CALL cryst_to_cart(nq_tot, Rlist, at, +1) + ! + ! Convert to crystal + ALLOCATE(strf_full(nq)) + IF (atom_resolved) ALLOCATE(strf_full_kk(nq, nat, nat)) + ! + CALL fkbounds( nq_tot, lower_bnd, upper_bnd ) + ! + ALLOCATE(Q_ppkkaa(upper_bnd - lower_bnd + 1, nat3, nat3)) + ! allocate array in parallel / minimize memory burden + Q_ppkkaa = 0.0d0 + ! DO qp2 = lower_bnd, upper_bnd ! for q-points in supercell ctr = 1 DO p1 = lower_bnd, upper_bnd !1, nq_tot DO qp2 = 1, nq_super !lower_bnd, upper_bnd ! for q-points in supercell @@ -2822,13 +2601,17 @@ IF (full_phonon) THEN DO kp = 1, nat ! k represents the atom ntap = ityp(kp) IF (qlistA(qp2) .EQ. 1) THEN - Q_ppkkaa(ctr, (k-1)*3+i, (kp-1)*3+p) = Q_ppkkaa(ctr, (k-1)*3+i, (kp-1)*3+p) & - + 1.0d0 / DBLE(SQRT(amass(nta)*amass(ntap))) / nq_tot * l_q(j, qp2)**2.0d0 & - * DBLE(z_zg((k-1)*3+i,j,qp2)*CONJG(z_zg((kp-1)*3+p,j,qp2))*EXP(imagi*2.0d0*pi*dotp4)) + Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) = & + Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) + & + 1.0d0 / DBLE(SQRT(amass(nta) * amass(ntap))) / nq_tot * l_q(j, qp2)**2.0d0 * & + DBLE(z_zg((k - 1) * 3 + i, j, qp2) * CONJG(z_zg((kp - 1) * 3 + p, j, qp2)) * & + EXP(imagi * 2.0d0 * pi * dotp4)) ELSE - Q_ppkkaa(ctr, (k-1)*3+i, (kp-1)*3+p) = Q_ppkkaa(ctr, (k-1)*3+i, (kp-1)*3+p) & - + 2.0d0 / DBLE(SQRT(amass(nta)*amass(ntap))) / nq_tot * l_q(j, qp2)**2.0d0 & - * DBLE(z_zg((k-1)*3+i,j,qp2)*CONJG(z_zg((kp-1)*3+p,j,qp2))*EXP(imagi*2.0d0*pi*dotp4)) + Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) = & + Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) + & + 2.0d0 / DBLE(SQRT(amass(nta) * amass(ntap))) / nq_tot * l_q(j, qp2)**2.0d0 * & + DBLE(z_zg((k - 1) * 3 + i, j, qp2) * CONJG(z_zg((kp - 1) * 3 + p, j, qp2)) * & + EXP(imagi * 2.0d0 * pi * dotp4)) ENDIF ENDDO ! kp ENDDO ! k @@ -2837,91 +2620,126 @@ IF (full_phonon) THEN ENDDO ! j ENDDO ! qp2 ctr = ctr + 1 - ENDDO ! p1 - ! - strctr_fact_full = 0.0d0 - ! - IF (atom_resolved) strctr_fact_full_kk = 0.0d0 - ! - ctr = 1 - DO p1 = lower_bnd, upper_bnd !1, nq_tot - DO qp = nq_super + 1, nq !lower_bnd + nq_super, upper_bnd + nq_super - dotp3 = 0.0d0 - DO ii = 1, 3 - dotp3 = dotp3 + q(ii, qp) * Rlist(ii, p1)! to get EXP (iS.(Rp-Rp')) - ENDDO - DO k = 1, nat ! k represents the atom - DO kp = 1, nat ! kp represents the atom - dotp2 = 0.0d0 - DO ii = 1, 3 - dotp2 = dotp2 + q(ii, qp) * (tau(ii, k) - tau(ii, kp))! - ENDDO - dotp = 0.0d0 - DO i = 1, 3 ! i is for cart directions - DO p = 1, 3 ! p is also for cartesian - dotp = dotp + q(i, qp) * q(p, qp) * (2.0d0 * pi / alat / BOHR_RADIUS_ANGS)**2.0d0 & - * Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) - ! dotp is P_p,kk' + ENDDO ! p1 + ! + strf_full = 0.0d0 + ! + IF (atom_resolved) strf_full_kk = 0.0d0 + ! + ctr = 1 + DO p1 = lower_bnd, upper_bnd !1, nq_tot + DO qp = nq_super + 1, nq !lower_bnd + nq_super, upper_bnd + nq_super + dotp3 = 0.0d0 + DO ii = 1, 3 + dotp3 = dotp3 + q(ii, qp) * Rlist(ii, p1)! to get EXP (iS.(Rp-Rp')) + ENDDO + DO k = 1, nat ! k represents the atom + DO kp = 1, nat ! kp represents the atom + dotp2 = 0.0d0 + DO ii = 1, 3 + dotp2 = dotp2 + q(ii, qp) * (tau(ii, k) - tau(ii, kp))! + ENDDO + dotp = 0.0d0 + DO i = 1, 3 ! i is for cart directions + DO p = 1, 3 ! p is also for cartesian + dotp = dotp + q(i, qp) * q(p, qp) * units**2.0d0 & ! dotp is P_p,kk' + * Q_ppkkaa(ctr, (k - 1) * 3 + i, (kp - 1) * 3 + p) + ENDDO ENDDO - ENDDO - ! - strctr_fact_full(qp) = strctr_fact_full(qp) + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp))& - * EXP(dotp) * EXP(imagi * 2.0d0 * pi * dotp2) * EXP(imagi * 2.0d0 * pi * dotp3) & - * atm_ffct(k, qp) * atm_ffct(kp, qp) * nq_tot - ! strctr_fact_full(qp) = strctr_fact_full(qp) + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp))& - ! * EXP(dotp) * EXP(imagi * 2.0d0 * pi * dotp2) * EXP(imagi * 2.0d0 * pi * dotp3) & - ! * atm_ffct(k, qp) * atm_ffct(kp, qp) - IF (atom_resolved) THEN - strctr_fact_full_kk(qp, k, kp) = strctr_fact_full_kk(qp, k, kp) & - + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp))& - * EXP(dotp) * EXP(imagi * 2.0d0 * pi * dotp2) * EXP(imagi * 2.0d0 * pi * dotp3) & - * atm_ffct(k, qp) * atm_ffct(kp, qp) * nq_tot - ENDIF - ! WRITE(*,*) "step", p2 - ENDDO ! kp - ENDDO ! k - ENDDO ! qp - ctr = ctr + 1 - ENDDO ! p1 - CALL mp_sum( strctr_fact_full, inter_pool_comm ) - CALL mp_barrier(inter_pool_comm) - IF (atom_resolved) CALL mp_sum( strctr_fact_full_kk, inter_pool_comm ) - CALL mp_barrier(inter_pool_comm) - !!!!!!!!!!!!!!!!!!!!!!!!!111 - IF (ionode) THEN - filename = 'structure_factor_all-phonon.dat' - OPEN (unit = 66, file = filename, status = 'unknown', form = 'formatted') - DO qp = nq_super + 1, nq - WRITE (66, '(40f26.6)') q(:,qp)*(2*pi/alat/BOHR_RADIUS_ANGS), DBLE(strctr_fact_full(qp)) !, aimag(strctr_fact_full(qp)) + ! + strf_full(qp) = strf_full(qp) + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) & + * EXP(dotp) * EXP(imagi * 2.0d0 * pi * dotp2) * EXP(imagi * 2.0d0 * pi * dotp3) & + * atm_ffct(k, qp) * atm_ffct(kp, qp) * nq_tot + ! here we can replace EXP(dotp) with 1 + dotp + IF (atom_resolved) THEN + strf_full_kk(qp, k, kp) = strf_full_kk(qp, k, kp) & + + EXP(-DW_T_fact_kkp(qp, k)) * EXP(-DW_T_fact_kkp(qp, kp)) & + * EXP(dotp) * EXP(imagi * 2.0d0 * pi * dotp2) * EXP(imagi * 2.0d0 * pi * dotp3) & + * atm_ffct(k, qp) * atm_ffct(kp, qp) * nq_tot + ENDIF + ENDDO ! kp + ENDDO ! k + ENDDO ! qp + ctr = ctr + 1 + ENDDO ! p1 + CALL mp_sum( strf_full, inter_pool_comm ) + CALL mp_barrier(inter_pool_comm) + IF (atom_resolved) CALL mp_sum( strf_full_kk, inter_pool_comm ) + CALL mp_barrier(inter_pool_comm) + ! To allocate a new matrix without negative values for q coordinates + ctr2 = 0 + DO p = nq_super + 1, nq + IF ((q(col1, p) .GT. 0.d0 - eps) .AND. (q(col2, p) .GT. 0.d0 - eps)) THEN + ctr2 = ctr2 + 1 + ENDIF ENDDO - CLOSE(66) + ! To rotate, apply broadening, and print the map + ALLOCATE(strf_rot(ctr2 * nrots, 4)) + ! + strf_rot = 0.d0 + ! + ! strf_full = strf_full - strf_q0 + IF (ionode) CALL rotate(DBLE(strf_full), q, nq, nq_super, nrots, ctr2, & + strf_rot, col1, col2) + CALL mp_bcast(strf_rot, ionode_id, world_comm) + CALL disca_broadening(strf_rot, ctr2 * nrots, kres1, kres2, alat, & + kmin, kmax, col1, col2, Np, 'strf_all-ph_broad.dat') + ! IF (atom_resolved) THEN DO k = 1, nat - WRITE( pointer_mz,'(i1.1)') k - filename = 'structure_factor_q_full_k_resolved_' // TRIM( pointer_mz ) // '.dat' !'.fp' + DO kp = 1, k + strf_rot = 0.d0 + WRITE( pt_mz,'(i1.1)') k + WRITE( pt_mz2,'(i1.1)') kp + filename = 'strf_all_k_'//TRIM(pt_mz)//'_'//TRIM(pt_mz2)//'_broad.dat' + ! + IF (ionode) CALL rotate(DBLE(strf_full_kk(:, k, kp)), & + q, nq, nq_super, nrots, ctr2, strf_rot, col1, col2) + CALL mp_bcast(strf_rot, ionode_id, world_comm) + CALL disca_broadening(strf_rot, ctr2 * nrots, kres1, kres2, alat, & + kmin, kmax, col1, col2, Np, filename) + !'zero_one_k_'//TRIM()//'_'//TRIM(pt_mz2)//'_broad.dat') + ENDDO ! kp + ENDDO + ENDIF + ! + DEALLOCATE(strf_rot) + ! + ! + IF (print_raw_data) THEN + IF (ionode) THEN + filename = 'strf_all-phonon.dat' OPEN (unit = 66, file = filename, status = 'unknown', form = 'formatted') DO qp = nq_super + 1, nq - WRITE (66, '(40f26.6)') q(:, qp) * (2.0d0 * pi / alat / BOHR_RADIUS_ANGS), & - ( DBLE((strctr_fact_full_kk(qp, k, kp) + strctr_fact_full_kk(qp, kp, k)) / 2.0d0), kp = 1, k) - ! NOTE: We should take the real part if and only if we sum - ! [strctr_fact_full_kk(qp, 2, 1) + strctr_fact_full_kk(qp, 1, 2)]/2. - ! strctr_fact_full_kk(qp, 2, 1) is not necessarilly real - !, aimag(strctr_fact_full(qp)) - ENDDO ! qp + WRITE (66, '(40f26.6)') q(:,qp) * units, DBLE(strf_full(qp)) !, aimag(strf_full(qp)) + ENDDO CLOSE(66) - ENDDO - ENDIF ! atom_resolved - ENDIF - ! - DEALLOCATE(Rlist, strctr_fact_full, Q_ppkkaa) - IF (atom_resolved) DEALLOCATE(strctr_fact_full_kk) -! -ENDIF ! full_phonon -! -DEALLOCATE(DW_T_fact_kkp) -RETURN ! if we DO not want to go further + IF (atom_resolved) THEN + DO k = 1, nat + WRITE( pt_mz,'(i1.1)') k + filename = 'strf_q_full_k_resolved_' // TRIM( pt_mz ) // '.dat' !'.fp' + OPEN (unit = 66, file = filename, status = 'unknown', form = 'formatted') + DO qp = nq_super + 1, nq + WRITE (66, '(40f26.6)') q(:, qp) * units, & + ( DBLE((strf_full_kk(qp, k, kp) + strf_full_kk(qp, kp, k)) / 2.0d0), kp = 1, k) + ! We take the real part iff we sum [strf_full_kk(qp, 2, 1) + strf_full_kk(qp, 1, 2)]/2. + ! strf_full_kk(qp, 2, 1) is not necessarilly real + ENDDO ! qp + CLOSE(66) + ENDDO + ENDIF ! atom_resolved + ENDIF + ENDIF ! print_raw_data + ! + DEALLOCATE(Rlist, strf_full, Q_ppkkaa) + IF (atom_resolved) DEALLOCATE(strf_full_kk) ! + ENDIF ! full_phonon ! + DEALLOCATE(DW_T_fact_kkp) + RETURN ! if we DO not want to go further + ! + ! END SUBROUTINE diffuse_scattering ! !SUBROUTINE qpoint_gen1_serial(dim1, dim2, dim3, ctrAB) @@ -2948,9 +2766,9 @@ END SUBROUTINE diffuse_scattering ! ! this is nothing but consecutive ordering ! n = (k - 1) + (j - 1) * dim3 + (i - 1) * dim2 * dim3 + 1 ! ! q_all are the components of the complete grid in crystal axis -! q_all(1, n) = dble(i - 1) / dim1 ! + dble(k1)/2/dim1 -! q_all(2, n) = dble(j - 1) / dim2 ! + dble(k2)/2/dim2 -! q_all(3, n) = dble(k - 1) / dim3 ! + dble(k3)/2/dim3 ! k1 , k2 , k3 is for the shift +! q_all(1, n) = DBLE(i - 1) / dim1 ! + DBLE(k1)/2/dim1 +! q_all(2, n) = DBLE(j - 1) / dim2 ! + DBLE(k2)/2/dim2 +! q_all(3, n) = DBLE(k - 1) / dim3 ! + DBLE(k3)/2/dim3 ! k1 , k2 , k3 is for the shift ! ENDDO ! ENDDO ! ENDDO @@ -3005,9 +2823,9 @@ END SUBROUTINE diffuse_scattering ! ! this is nothing but consecutive ordering ! n = (k - 1) + (j - 1) * dim3 + (i - 1) * dim2 * dim3 + 1 ! ! q_all are the components of the complete grid in crystal axis -! q_all(1, n) = dble(i - 1) / dim1 ! + dble(k1)/2/dim1 -! q_all(2, n) = dble(j - 1) / dim2 ! + dble(k2)/2/dim2 -! q_all(3, n) = dble(k - 1) / dim3 ! + dble(k3)/2/dim3 ! k1 , k2 , k3 is for the shift +! q_all(1, n) = DBLE(i - 1) / dim1 ! + DBLE(k1)/2/dim1 +! q_all(2, n) = DBLE(j - 1) / dim2 ! + DBLE(k2)/2/dim2 +! q_all(3, n) = DBLE(k - 1) / dim3 ! + DBLE(k3)/2/dim3 ! k1 , k2 , k3 is for the shift ! ENDDO ! ENDDO ! ENDDO @@ -3101,7 +2919,7 @@ END SUBROUTINE diffuse_scattering ! SUBROUTINE qpoint_gen_map_1(nk1, nk2, nk3, k1, k2, k3, kf1, kf2, kf3, & - plane_val, plane_dir, nq_strft, qlist_str_f, qlist_str_f_cart) + plane_val, plane_dir, nq_strft, qlist_strf, qlist_strf_cart) ! USE kinds, ONLY: DP USE cell_base, ONLY : bg @@ -3110,8 +2928,8 @@ IMPLICIT NONE INTEGER, intent(in) :: nk1, nk2, nk3, k1, k2, k3, kf1, kf2, kf3, plane_dir REAL(DP), intent(in) :: plane_val INTEGER, intent(out) :: nq_strft - REAL(DP), intent(out) :: qlist_str_f(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) - REAL(DP), intent(out) :: qlist_str_f_cart(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) + REAL(DP), intent(out) :: qlist_strf(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) + REAL(DP), intent(out) :: qlist_strf_cart(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) INTEGER :: i, j, k, n, ctr, p, nk_tot REAL(DP), ALLOCATABLE :: xkg(:, :), xkg_or(:, :) REAL(DP) :: q_B(3), q_A(3), eps, B(3, 3) @@ -3119,9 +2937,6 @@ IMPLICIT NONE ! Reciprocal axis, as printed form espresso but transpose ! B = (bg) -!WRITE(*,*) B(1,:) -!WRITE(*,*) B(2,:) -!WRITE(*,*) B(3,:) ! ! eps = 1.0E-05 @@ -3134,11 +2949,11 @@ B = (bg) DO j = 1, nk2 DO k = 1, nk3 ! this is nothing but consecutive ordering - n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 + n = (k - 1) + (j - 1) * nk3 + (i - 1) * nk2 * nk3 + 1 ! xkg are the components of the complete grid in crystal axis - xkg(1,n) = dble(i-1)/nk1 ! + dble(k1)/2/nk1 - xkg(2,n) = dble(j-1)/nk2 ! + dble(k2)/2/nk2 - xkg(3,n) = dble(k-1)/nk3 ! + dble(k3)/2/nk3 ! k1 , k2 , k3 is for the shift + xkg(1, n) = DBLE(i - 1) / nk1 + xkg(2, n) = DBLE(j - 1) / nk2 + xkg(3, n) = DBLE(k - 1) / nk3 ENDDO ENDDO ENDDO @@ -3149,12 +2964,12 @@ B = (bg) DO p = k1, kf1 - 1 DO j = k2, kf2 - 1 DO k = k3, kf3 - 1 - xkg(1,:) = xkg_or(1,:) + p - xkg(2,:) = xkg_or(2,:) + j - xkg(3,:) = xkg_or(3,:) + k + xkg(1, :) = xkg_or(1, :) + p + xkg(2, :) = xkg_or(2, :) + j + xkg(3, :) = xkg_or(3, :) + k DO i = 1, nk_tot - qlist_str_f(ctr,:) = xkg(:,i) - qlist_str_f_cart(ctr,:) = B(:, 1) * xkg(1, i) + B(:, 2) * xkg(2, i) + B(:, 3) * xkg(3, i) + qlist_strf(ctr, :) = xkg(:, i) + qlist_strf_cart(ctr, :) = B(:, 1) * xkg(1, i) + B(:, 2) * xkg(2, i) + B(:, 3) * xkg(3, i) ctr= ctr+1 END DO END DO @@ -3164,10 +2979,10 @@ B = (bg) ! To allocate matrix size nq_strft = 0 DO i = 1, nk_tot * (kf1 - k1) * (kf2 - k2) * (kf3 - k3) - IF ( (qlist_str_f_cart(i, plane_dir) .LT. plane_val + eps) .AND. & - (qlist_str_f_cart(i, plane_dir) .GT. plane_val - eps) ) THEN - nq_strft = nq_strft + 1 - ENDIF + IF ( (qlist_strf_cart(i, plane_dir) .LT. plane_val + eps) .AND. & + (qlist_strf_cart(i, plane_dir) .GT. plane_val - eps) ) THEN + nq_strft = nq_strft + 1 + ENDIF ENDDO ! DEALLOCATE(xkg, xkg_or) @@ -3176,39 +2991,37 @@ END SUBROUTINE qpoint_gen_map_1 ! ! SUBROUTINE qpoint_gen_map_2(nk1, nk2, nk3, k1, k2, k3, kf1, kf2, kf3, & - plane_val, plane_dir, nq_strft, qlist_str_f, qlist_str_f_cart, qlist_str_f_cryst_z) + plane_val, plane_dir, nq_strft, qlist_strf, qlist_strf_cart, qlist_strf_cryst_z) ! USE kinds, ONLY: DP USE cell_base, ONLY : at IMPLICIT NONE INTEGER, intent(in) :: nk1, nk2, nk3, k1, k2, k3, nq_strft, plane_dir, kf1, kf2, kf3 - REAL(DP), intent(in) :: plane_val, qlist_str_f(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) - REAL(DP), intent(in) :: qlist_str_f_cart(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) - REAL(DP), intent(out) :: qlist_str_f_cryst_z(3, nq_strft) + REAL(DP), intent(in) :: plane_val, qlist_strf(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) + REAL(DP), intent(in) :: qlist_strf_cart(nk1 * nk2 * nk3 * (kf1 - k1) * (kf2 - k2) * (kf3 - k3), 3) + REAL(DP), intent(out) :: qlist_strf_cryst_z(3, nq_strft) INTEGER :: i, ctr, nk_tot REAL(DP) :: eps, invB(3, 3) ! invB = TRANSPOSE(at) -!WRITE(*,*) invB(1,:) -!WRITE(*,*) invB(2,:) -!WRITE(*,*) invB(3,:) ! Crystal axis, as printed form espresso but NOT transpose ! ! eps = 1.0E-05 ! - nk_tot= nk1 * nk2 * nk3 ! nk1, nk2 ,nk3 define the resolution of the q-grid + nk_tot= nk1 * nk2 * nk3 + ! nk1, nk2 ,nk3 define the resolution of the q-grid ! ctr = 0 DO i = 1, nk_tot * (kf1 - k1) * (kf2 - k2) * (kf3 - k3) - IF ( (qlist_str_f_cart(i, plane_dir) .LT. plane_val + eps) .AND. & - (qlist_str_f_cart(i, plane_dir) .GT. plane_val - eps) ) THEN + IF ( (qlist_strf_cart(i, plane_dir) .LT. plane_val + eps) .AND. & + (qlist_strf_cart(i, plane_dir) .GT. plane_val - eps) ) THEN ctr = ctr + 1 -! qlist_str_f_cart_z(:, ctr) = qlist_str_f_cart(i,:) - qlist_str_f_cryst_z(:, ctr) = invB(:, 1) * qlist_str_f_cart(i, 1) + & - invB(:, 2) * qlist_str_f_cart(i, 2) + & - invB(:, 3) * qlist_str_f_cart(i, 3) +! qlist_strf_cart_z(:, ctr) = qlist_strf_cart(i,:) + qlist_strf_cryst_z(:, ctr) = invB(:, 1) * qlist_strf_cart(i, 1) + & + invB(:, 2) * qlist_strf_cart(i, 2) + & + invB(:, 3) * qlist_strf_cart(i, 3) ENDIF ENDDO ! @@ -3226,7 +3039,6 @@ SUBROUTINE qpoint_gen1(dim1, dim2, dim3, ctrAB) ! input 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 @@ -3247,9 +3059,9 @@ SUBROUTINE qpoint_gen1(dim1, dim2, dim3, ctrAB) ! this is nothing but consecutive ordering n = (k - 1) + (j - 1) * dim3 + (i - 1) * dim2 * dim3 + 1 ! q_all are the components of the complete grid in crystal axis - q_all(1, n) = dble(i - 1) / dim1 ! + dble(k1)/2/dim1 - q_all(2, n) = dble(j - 1) / dim2 ! + dble(k2)/2/dim2 - q_all(3, n) = dble(k - 1) / dim3 ! + dble(k3)/2/dim3 ! k1 , k2 , k3 is for the shift + q_all(1, n) = DBLE(i - 1) / dim1 ! + DBLE(k1)/2/dim1 + q_all(2, n) = DBLE(j - 1) / dim2 ! + DBLE(k2)/2/dim2 + q_all(3, n) = DBLE(k - 1) / dim3 ! + DBLE(k3)/2/dim3 ! k1 , k2 , k3 is for the shift ENDDO ENDDO ENDDO @@ -3315,9 +3127,9 @@ SUBROUTINE qpoint_gen2(dim1, dim2, dim3, ctrAB, q_AB) ! this is nothing but consecutive ordering n = (k - 1) + (j - 1) * dim3 + (i - 1) * dim2 * dim3 + 1 ! q_all are the components of the complete grid in crystal axis - q_all(1, n) = dble(i - 1) / dim1 ! + dble(k1)/2/dim1 - q_all(2, n) = dble(j - 1) / dim2 ! + dble(k2)/2/dim2 - q_all(3, n) = dble(k - 1) / dim3 ! + dble(k3)/2/dim3 ! k1 , k2 , k3 is for the shift + q_all(1, n) = DBLE(i - 1) / dim1 + q_all(2, n) = DBLE(j - 1) / dim2 + q_all(3, n) = DBLE(k - 1) / dim3 ENDDO ENDDO ENDDO @@ -3334,7 +3146,6 @@ SUBROUTINE qpoint_gen2(dim1, dim2, dim3, ctrAB, q_AB) ((ABS(q_A(2)) .LT. eps) .OR. (ABS(ABS(q_A(2)) - 1) .LT. eps)) .AND. & ((ABS(q_A(3)) .LT. eps) .OR. (ABS(ABS(q_A(3)) - 1) .LT. eps))) THEN q_AB_TMP(:, i) = q_all(:, i) - ! WRITE(*,*) "A", q_AB(:, ctr) ELSE DO j = i + 1, nqs q_B = q_all(:, i) + q_all(:, j) @@ -3342,7 +3153,6 @@ SUBROUTINE qpoint_gen2(dim1, dim2, dim3, ctrAB, q_AB) ((ABS(q_B(2)) .LT. eps) .OR. (ABS(ABS(q_B(2)) - 1) .LT. eps)) .AND. & ((ABS(q_B(3)) .LT. eps) .OR. (ABS(ABS(q_B(3)) - 1) .LT. eps))) THEN q_AB_TMP(:, i) = q_all(:, i) - ! WRITE(*,*) q_AB(:, ctr) END IF END DO END IF @@ -3364,3 +3174,257 @@ SUBROUTINE qpoint_gen2(dim1, dim2, dim3, ctrAB, q_AB) ! END SUBROUTINE qpoint_gen2 +! +!----------------------------------------------------------------------- +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 + ! + USE kinds, ONLY : DP + USE constants, ONLY : tpi + USE io_global, ONLY : stdout + ! + IMPLICIT NONE + INTEGER nr1, nr2, nr3, nat, n1, n2, n3, nr1_, nr2_, nr3_, & + ipol, jpol, na, nb, m1, m2, m3, nint, i,j, nrws, nax + COMPLEX(DP) dyn(3,3,nat,nat), f_of_q(3,3,nat,nat) + REAL(DP) frc(nr1,nr2,nr3,3,3,nat,nat), tau(3,nat), q(3), arg, & + at(3,3), bg(3,3), r(3), weight, r_ws(3), & + total_weight, rws(0:3,nrws), alat + REAL(DP), EXTERNAL :: wsweight + REAL(DP),SAVE,ALLOCATABLE :: wscache(:,:,:,:,:) + REAL(DP), ALLOCATABLE :: ttt(:,:,:,:,:), tttx(:,:) + LOGICAL,SAVE :: first=.true. + LOGICAL :: fd + ! + nr1_=2*nr1 + nr2_=2*nr2 + nr3_=2*nr3 + FIRST_TIME : IF (first) THEN + first=.false. + ALLOCATE( wscache(-nr3_:nr3_, -nr2_:nr2_, -nr1_:nr1_, nat,nat) ) + DO na=1, nat + DO nb=1, nat + total_weight=0.0d0 + ! + ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY VERY SAFE RANGE! + ! + DO n1=-nr1_,nr1_ + DO n2=-nr2_,nr2_ + DO n3=-nr3_,nr3_ + 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) + END DO + wscache(n3,n2,n1,nb,na) = wsweight(r_ws,rws,nrws) + total_weight=total_weight + wscache(n3,n2,n1,nb,na) + ENDDO + ENDDO + ENDDO + IF (ABS(total_weight-nr1*nr2*nr3).GT.1.0d-8) THEN + WRITE(stdout,*) na,nb,total_weight + CALL errore ('frc_blk','wrong total_weight',1) + END IF + ENDDO + ENDDO + ENDIF FIRST_TIME + ! + ALLOCATE(ttt(3,nat,nr1,nr2,nr3)) + ALLOCATE(tttx(3,nat*nr1*nr2*nr3)) + ttt(:,:,:,:,:)=0.d0 + + DO na=1, nat + DO nb=1, nat + DO n1=-nr1_,nr1_ + DO n2=-nr2_,nr2_ + DO n3=-nr3_,nr3_ + ! + ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY SAFE RANGE! + ! + DO i=1, 3 + r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3) + END DO + + weight = wscache(n3,n2,n1,nb,na) + IF (weight .GT. 0.0d0) THEN + ! + ! FIND THE VECTOR CORRESPONDING TO R IN THE ORIGINAL CELL + ! + m1 = MOD(n1+1,nr1) + IF(m1.LE.0) m1=m1+nr1 + m2 = MOD(n2+1,nr2) + IF(m2.LE.0) m2=m2+nr2 + m3 = MOD(n3+1,nr3) + IF(m3.LE.0) m3=m3+nr3 + ! write(*,'(6i4)') n1,n2,n3,m1,m2,m3 + ! + ! FOURIER TRANSFORM + ! + DO i=1,3 + ttt(i,na,m1,m2,m3)=tau(i,na)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) + ttt(i,nb,m1,m2,m3)=tau(i,nb)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) + END DO + + arg = tpi*(q(1)*r(1) + q(2)*r(2) + q(3)*r(3)) + DO ipol=1, 3 + DO jpol=1, 3 + dyn(ipol,jpol,na,nb) = dyn(ipol,jpol,na,nb) + & + (frc(m1,m2,m3,ipol,jpol,na,nb)+f_of_q(ipol,jpol,na,nb)) & + *CMPLX(COS(arg),-SIN(arg),kind=DP)*weight + END DO + END DO + + END IF + END DO + END DO + END DO + END DO + END DO + ! + RETURN +END SUBROUTINE frc_blk +! +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, BOHR_RADIUS_ANGS +! +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, units !, pi + REAL(DP), ALLOCATABLE :: kgridx(:), kgridy(:), strf_out(:, :) +! +! +! +units = DBLE(2.d0*pi/alat/BOHR_RADIUS_ANGS) +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) * units - kgridx(ik))**2.d0 / sf_smearingx**2.d0 / 2.d0))*& + (EXP(-(strf_rot(ii, col2) * units - 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 diff --git a/EPW/ZG/src/pp_disca.f90 b/EPW/ZG/src/pp_disca.f90 index f881e615c..3f02cdc27 100644 --- a/EPW/ZG/src/pp_disca.f90 +++ b/EPW/ZG/src/pp_disca.f90 @@ -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 diff --git a/EPW/ZG/tutorial.tar.gz b/EPW/ZG/tutorial.tar.gz new file mode 100644 index 000000000..0703be865 Binary files /dev/null and b/EPW/ZG/tutorial.tar.gz differ diff --git a/EPW/examples/diamond/epw/epw4.in b/EPW/examples/diamond/epw/epw4.in new file mode 100644 index 000000000..ac3867735 --- /dev/null +++ b/EPW/examples/diamond/epw/epw4.in @@ -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 + / diff --git a/EPW/src/bcast_epw_input.f90 b/EPW/src/bcast_epw_input.f90 index 847190887..2ef5c8e32 100644 --- a/EPW/src/bcast_epw_input.f90 +++ b/EPW/src/bcast_epw_input.f90 @@ -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 ! diff --git a/EPW/src/cum_mod.f90 b/EPW/src/cum_mod.f90 index 4001279ce..a27e97dbf 100644 --- a/EPW/src/cum_mod.f90 +++ b/EPW/src/cum_mod.f90 @@ -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, 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(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 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 !----------------------------------------------------------------------- diff --git a/EPW/src/elph2.f90 b/EPW/src/elph2.f90 index 25b9a20a6..04194f05f 100644 --- a/EPW/src/elph2.f90 +++ b/EPW/src/elph2.f90 @@ -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 diff --git a/EPW/src/ephwann_shuffle.f90 b/EPW/src/ephwann_shuffle.f90 index b8f005982..bfbaf7917 100644 --- a/EPW/src/ephwann_shuffle.f90 +++ b/EPW/src/ephwann_shuffle.f90 @@ -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 diff --git a/EPW/src/epw.f90 b/EPW/src/epw.f90 index 944104aea..aa7aab878 100644 --- a/EPW/src/epw.f90 +++ b/EPW/src/epw.f90 @@ -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.) ! diff --git a/EPW/src/epw_readin.f90 b/EPW/src/epw_readin.f90 index 88e4c9f3f..8fc298bc8 100644 --- a/EPW/src/epw_readin.f90 +++ b/EPW/src/epw_readin.f90 @@ -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 diff --git a/EPW/src/epwcom.f90 b/EPW/src/epwcom.f90 index 3ec6afc2c..6f109917b 100644 --- a/EPW/src/epwcom.f90 +++ b/EPW/src/epwcom.f90 @@ -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 diff --git a/EPW/src/io_eliashberg.f90 b/EPW/src/io_eliashberg.f90 index 3197146ec..475cc685c 100644 --- a/EPW/src/io_eliashberg.f90 +++ b/EPW/src/io_eliashberg.f90 @@ -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 ! diff --git a/EPW/src/io_transport.f90 b/EPW/src/io_transport.f90 index cf24700c1..0c12708c7 100644 --- a/EPW/src/io_transport.f90 +++ b/EPW/src/io_transport.f90 @@ -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 ! diff --git a/EPW/src/low_lvl.f90 b/EPW/src/low_lvl.f90 index b50057b32..8da53ad49 100644 --- a/EPW/src/low_lvl.f90 +++ b/EPW/src/low_lvl.f90 @@ -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) diff --git a/EPW/src/printing.f90 b/EPW/src/printing.f90 index c813a0dd0..415f938dd 100644 --- a/EPW/src/printing.f90 +++ b/EPW/src/printing.f90 @@ -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) ! diff --git a/EPW/src/supercond.f90 b/EPW/src/supercond.f90 index d439df223..3f16579a3 100644 --- a/EPW/src/supercond.f90 +++ b/EPW/src/supercond.f90 @@ -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 ! diff --git a/EPW/src/supercond_aniso.f90 b/EPW/src/supercond_aniso.f90 index aeda80eee..63dafd499 100644 --- a/EPW/src/supercond_aniso.f90 +++ b/EPW/src/supercond_aniso.f90 @@ -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 diff --git a/EPW/src/wan2bloch.f90 b/EPW/src/wan2bloch.f90 index baf8b319b..d1e68bd22 100644 --- a/EPW/src/wan2bloch.f90 +++ b/EPW/src/wan2bloch.f90 @@ -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 diff --git a/test-suite/epw_mob_ibte_sym/benchmark.out.git.inp=epw9.in.args=5 b/test-suite/epw_mob_ibte_sym/benchmark.out.git.inp=epw9.in.args=5 new file mode 100644 index 000000000..29cd813fe --- /dev/null +++ b/test-suite/epw_mob_ibte_sym/benchmark.out.git.inp=epw9.in.args=5 @@ -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. + =============================================================================== + diff --git a/test-suite/epw_mob_ibte_sym/epw9.in b/test-suite/epw_mob_ibte_sym/epw9.in new file mode 100644 index 000000000..319743a53 --- /dev/null +++ b/test-suite/epw_mob_ibte_sym/epw9.in @@ -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 + / diff --git a/test-suite/jobconfig b/test-suite/jobconfig index bd0496c02..28e3c2b97 100644 --- a/test-suite/jobconfig +++ b/test-suite/jobconfig @@ -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