! !----------------------------------------------------------------------- subroutine elphon_shuffle_wrap !----------------------------------------------------------------------- ! ! Electron-phonon calculation with Wannier functions: load all phonon q's ! ! Explanations and comments to come... ! !----------------------------------------------------------------------- ! #include "f_defs.h" ! #ifdef __PARA use para USE io_global, ONLY : ionode_id USE mp_global, ONLY : nproc, my_pool_id, nproc_pool, & intra_image_comm, inter_pool_comm, me_pool, root_pool, & mpime, intra_pool_comm USE mp, ONLY : mp_barrier, mp_bcast #endif use io_files, only : prefix, tmp_dir USE wavefunctions_module, ONLY: evc USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau use control_flags, ONLY : iswitch, iverbosity, modenum, nosym, noinv use io_global, only : stdout USE kinds, only : DP use pwcom, only : nr1, nr2, nr3, et, xk, irt, s, nsym, & ibrav, ftau, symm_type, igk, at, bg, sname, nks, nbnd, nkstot, ngm use phcom, only : npert, u, gi, nirr, t, irotmq, max_irr_dim, tmq, dpsi, dvpsi, & phonselfen, elecselfen, nbndsub, gimq, igkq, evq, nq1, ephwanout, nq3, nq2, & elinterp, nmodes, trans, phinterp, iunepmatf, minus_q, rtau, nsymq, irgq, & lgamma, iuetf, invs, xq, epbread, epbwrite, epwread, specfun use el_phon, only : epmatq, dynq, et_all, xk_all implicit none real(kind=DP), allocatable :: xqc_irr(:,:), wqlist_irr(:), xqc(:,:), wqlist(:) ! the qpoints in the irr wedge ! the corresponding weigths ! the qpoints in the uniform mesh ! the corresponding weigths integer :: nqc_irr, nqc, iuepb ! number of qpoints in the irreducible wedge ! number of qpoints on the uniform grid ! unit for .epb files ! symmetry-related variables ! integer :: gmapsym(ngm,48) ! correspondence G -> S(G) complex(kind=DP) :: eigv (ngm, 48) ! e^{ iGv} for 1...nsym (v the fractional translation) complex(kind=DP) :: cz1( nmodes, nmodes), cz2(nmodes, nmodes) ! the eigenvectors for the first q in the star ! the rotated eigenvectors, for the current q in the star ! integer :: nq, isq (48), imq, na, nb, nt, nta, ntb, jmode0, irr, jrr ! degeneracy of the star of q ! index of q in the star of a given sym.op. ! index of -q in the star of q (0 if not present) integer :: sym_sgq(48) ! the symmetries giving the q point iq in the star real(kind=DP) :: sxq (3, 48) ! list of vectors in the star of q integer :: ibnd, jbnd, ik, i, j, imode0, iq, iq_irr, isym, isgq, & imode, jmode, ipool, iq_first, ig, jsym, ism1, isg, itmp, nsq, ipol real(kind=DP) tmp, xq0(3), aq(3), saq(3), raq(3), sr(3,3), ft1, ft2, ft3 logical :: sym(48), found, eqvect_strict, nog, symmo character (len=256) :: tempfile character (len=3) :: filelab ! ! unit for the .epf and .etf file[s], temporary storage for interpolation ! and for the .epb files (e-ph matrix in bloch representation) ! iunepmatf = 76 iuetf = 78 iuepb = 85 ! ! if we read all quantities in the Wannier representation from file, ! we skip all this and we jump directly to the interpolation routine ! if ( .not. epwread ) then ! if ( elinterp .and. (.not.phinterp ) ) call errore & ('elphon_shuffle_wrap','elinterp requires phinterp' ,1) if ( trans ) call errore & ('elphon_shuffle_wrap','el-phon interp requires trans = F' ,1) ! ! print out infos for electron-phonon Wannier interpolation (including wfs) ! if ( ephwanout ) then call wanout endif ! ! read qpoint list from stdin ! #ifdef __PARA if (mpime.eq.ionode_id) & #endif read(5,*) nqc_irr call mp_bcast (nqc_irr, ionode_id, inter_pool_comm) call mp_bcast (nqc_irr, root_pool, intra_pool_comm) allocate ( xqc_irr(3,nqc_irr), wqlist_irr(nqc_irr) ) allocate ( xqc(3,nq1*nq2*nq3), wqlist(nq1*nq2*nq3) ) ! #ifdef __PARA if (mpime.eq.ionode_id) then #endif do iq = 1, nqc_irr read (5,*) xqc_irr (:,iq), wqlist_irr (iq) enddo if ( abs(sum(wqlist_irr)-2.d0) .gt. 1d-5 ) & call errore ('elphon_shuffle_wrap','q point weights do not add up to 2',1) #ifdef __PARA endif call mp_bcast (xqc_irr, ionode_id, inter_pool_comm) call mp_bcast (xqc_irr, root_pool, intra_pool_comm) #endif ! ! gather electronic eigenvalues for subsequent shuffle ! allocate ( xk_all( 3, nkstot) , et_all( nbnd, nkstot) ) call poolgather ( 3, nkstot, nks, xk(:,1:nks), xk_all) call poolgather ( nbnd, nkstot, nks, et(1:nbnd,1:nks), et_all) ! ! if we start with a gamma point calculation, ../PW/set_kplusq.f90 ! is not active and the gmap has not been produced... ! if (lgamma) call errore & ('elphon_shuffle_wrap','tshuffle2 requires q!=0 starting nscf calculation',1) ! ! allocate dynamical matrix and ep matrix for all q's ! allocate ( dynq (nmodes, nmodes, nq1*nq2*nq3), & epmatq (nbnd, nbnd, nks, nmodes, nq1*nq2*nq3) ) ! ! NB: the symmetries are q-dependent and we should reset them ! for every q in the star. However, since we use the rotation of the ! wavefunctions (see notes), the dvscf_q_nu must be exactly the same ! as for the originating qpoints, and so the patterns etc. ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ftau comes from the modules. Problem: when I call star_q to reset ! symmetries, in star_q.f90 ftau is a local variable. Therefore, in the ! call to sgama.f90 within star_q.f90 (which calls sgama_at.f90 where ! teh fractional translations are defined) the ftau are reset but this ! subroutine does not know that. And here I am getting the ftau corresponding ! to the weird q I use in the nscf calculation. This means that we have to ! get the ftau's from star_q. !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! if we do not have epmatq already on file, we calculate the sandwiches here ! if ( .not. epbread) then ! nqc = 0 iq_first = 1 do iq_irr = 1, nqc_irr ! write(6,'(//5x,a)') repeat('=',67) write(6,'(5x,"irreducible q point # ",i4)') iq_irr write(6,'(5x,a/)') repeat('=',67) ! xq = xqc_irr(:,iq_irr) ! ! store the first q in the star xq0 = xq ! ! if the nscf run is with q!=0, we will stay with lgamma = .false. throughout ! (it does not change anywhere) ! ! reset symmetries of the crystal (not the small group of q) ! and determine the star of q ! call star_q2 (xq, at, bg, ibrav, symm_type, nat, tau, ityp, nr1, nr2, nr3, & nsym, s, invs, irt, rtau, nq, sxq, isq, imq, noinv, modenum, ftau) ! ! determine the G vector map S(G) -> G ! NB: after star_q we have the global symms, not those of the small group of q: ! do not move this call. ! if (iq_irr.eq.1) then call gmap_sym ( nsym, s, ftau, gmapsym, eigv, invs) ! I checked that gmapsym(gmapsym(ig,isym),invs(isym)) = ig endif ! call readmat_shuffle2 ( iq_irr, nqc_irr, nq, iq_first, sxq, imq) ! now dynq is the cartesian dyn mat (NOT divided by the masses) ! ! re-set the variables needed for the pattern representation ! and the symmetries of the small group of irr-q ! (from phq_setup.f90) ! do isym = 1, nsym sym (isym) = .true. enddo call sgam_ph (at, bg, nsym, s, irt, tau, rtau, nat, sym) minus_q = (iswitch .gt. -3) ! call set_irr & (nat, at, bg, xq, s, invs, nsym, rtau, irt, irgq, nsymq, minus_q, & irotmq, t, tmq, max_irr_dim, u, npert, nirr, gi, gimq, iverbosity) ! ! if smallgq.f90 determined nsymq.eq.1, then the original phonon ! calculation has been done with the patterns set to the identity ! matrix, and we need to use the same convention here ! if ( nsymq .eq. 1 ) call set_irr_nosym & (nat, at, bg, xq, s, invs, nsym, rtau, irt, irgq, nsymq, minus_q, & irotmq, t, tmq, max_irr_dim, u, npert, nirr, gi, gimq, iverbosity) ! ! re-initialize the q-dependent quantities for the local ! and for the nonlocal part of the pseudopotentials ! call phq_init ! ! loop over the q points of the star ! do iq = 1, nq ! ! retrieve the q in the star xq = sxq(:,iq) ! ! and popolate the uniform grid nqc = nqc + 1 xqc(:,nqc) = xq ! if (iq.eq.1) write(6,*) write( stdout, 5) nqc, xq ! ! prepare the gmap for the refolding ! call createkmap ( xq ) ! if (iverbosity.eq.1) then ! ! description of symmetries ! WRITE( stdout, '(36x,"s",24x,"frac. trans.")') do isym = 1, nsym WRITE( stdout, '(/6x,"isym = ",i2,5x,a45/)') isym, sname(isym) call s_axis_to_cart (s(1,1,isym), sr, at, bg) if (ftau(1,isym).ne.0.or.ftau(2,isym).ne.0.or.ftau(3,isym).ne.0) then ft1 = at(1,1)*ftau(1,isym)/nr1 + at(1,2)*ftau(2,isym)/nr2 + & at(1,3)*ftau(3,isym)/nr3 ft2 = at(2,1)*ftau(1,isym)/nr1 + at(2,2)*ftau(2,isym)/nr2 + & at(2,3)*ftau(3,isym)/nr3 ft3 = at(3,1)*ftau(1,isym)/nr1 + at(3,2)*ftau(2,isym)/nr2 + & at(3,3)*ftau(3,isym)/nr3 WRITE( stdout, '(1x,"cryst.",3x,"s(",i2,") = (",3(i6,5x), & & " ) f =( ",f10.7," )")') & isym, (s(1,ipol,isym),ipol=1,3), dble(ftau(1,isym))/dble(nr1) WRITE( stdout, '(17x," (",3(i6,5x), " ) ( ",f10.7," )")') & (s(2,ipol,isym),ipol=1,3), dble(ftau(2,isym))/dble(nr2) WRITE( stdout, '(17x," (",3(i6,5x), " ) ( ",f10.7," )"/)') & (s(3,ipol,isym),ipol=1,3), dble(ftau(3,isym))/dble(nr3) WRITE( stdout, '(1x,"cart. ",3x,"s(",i2,") = (",3f11.7, & & " ) f =( ",f10.7," )")') & isym, (sr(1,ipol),ipol=1,3), ft1 WRITE( stdout, '(17x," (",3f11.7, " ) ( ",f10.7," )")') & (sr(2,ipol),ipol=1,3), ft2 WRITE( stdout, '(17x," (",3f11.7, " ) ( ",f10.7," )"/)') & (sr(3,ipol),ipol=1,3), ft3 else WRITE( stdout, '(1x,"cryst.",3x,"s(",i2,") = (",3(i6,5x), " )")') & isym, (s (1, ipol, isym) , ipol = 1,3) WRITE( stdout, '(17x," (",3(i6,5x)," )")') (s(2,ipol,isym), ipol=1,3) WRITE( stdout, '(17x," (",3(i6,5x)," )"/)') (s(3,ipol,isym), ipol=1,3) WRITE( stdout, '(1x,"cart. ",3x,"s(",i2,") = (",3f11.7," )")') & isym, (sr (1, ipol) , ipol = 1, 3) WRITE( stdout, '(17x," (",3f11.7," )")') (sr (2, ipol) , ipol = 1, 3) WRITE( stdout, '(17x," (",3f11.7," )"/)') (sr (3, ipol) , ipol = 1, 3) endif enddo ! endif ! ! isq(isym)=iq means: when we apply symmetry isym to the originating q ! of the star, we get the iq-th member of the star. There are as many ! matches as the degeneracy of the star. ! ! pick up the first q in the small group of q* ! [it is important to select the first q in the small group, since ! it really corresponds to Sq. If we choose another element in the small group ! the actual q-point may be Sq+G and we screw up the q-vector below to generate ! k+q from k and for the KB projectors ! nsq = 0 ! nsq is the degeneracy of the small group for this iq in the star ! do jsym = 1, nsym if ( isq(jsym) .eq. iq ) then nsq = nsq + 1 sym_sgq(nsq) = jsym endif enddo if ( nsq*nq .ne. nsym ) call errore ('elphon_shuffle_wrap', 'wrong degeneracy', iq) ! if (iverbosity.eq.1) then ! write(6,*) 'iq, i, isym, nog, symmo' do i = 1, nsq ! isym = sym_sgq(i) ism1 = invs (isym) ! ! check for G such that Sq = q* + G ! aq = xq0 saq = xq call cryst_to_cart (1, aq, at, -1) do j = 1, 3 raq (j) = s (j, 1, ism1) * aq (1) & + s (j, 2, ism1) * aq (2) & + s (j, 3, ism1) * aq (3) enddo call cryst_to_cart (1, saq, at, -1) nog = eqvect_strict( raq, saq) ! ! check whether the symmetry belongs to a symmorphic group ! symmo = (ftau(1,isym).eq.0 .and. ftau(2,isym).eq.0 .and. ftau(3,isym).eq.0) ! write(6,'(3i5,2a)') iq, i, isym, nog, symmo ! enddo ! endif ! isym = sym_sgq (1) ! primo q del piccolo gruppo, assicura che G=0 ! ! calculate the sandwiches ! ! a more accurate way of doing is to symmetrize the matrix element w.r.t. ! the small group of the given q in the star. I'm not doint this here. ! (but I checked that even without symm the result of full zone and irr zone ! are equal to 5+ digits) ! For any volunteers, please write to giustino@civet.berkeley.edu ! call elphon_shuffle ( iq_irr, nqc_irr, nqc, gmapsym, eigv, isym, invs, xq0, .false. ) ! ! bring epmatq in the mode representation of iq_first, ! and then in the cartesian representation of iq ! call rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, & rtau, xq, isq, cz1, cz2, .false. ) ! call rotate_epmat ( cz1, cz2, xq, nqc ) ! if (imq .eq. 0) then ! ! -q is not in the star: we introduce it here using time-reversal symmetry ! [Eq. (4.4) Maradudin&Vsoko RMP] ! call errore('elphon_shuffle_wrap','-q is not in the star of q: never tested',0) ! ! retrieve the -q in the star xq = - xq nqc = nqc + 1 xqc(:,nqc) = xq ! write( stdout, 5) nqc, xq ! ! prepare the gmap for the refolding ! call createkmap ( xq ) ! call elphon_shuffle ( iq_irr, nqc_irr, nqc, gmapsym, eigv, isym, invs, xq0, .true. ) ! ! bring epmatq in the mode representation of iq_first, ! and then in the cartesian representation of iq ! taking into account time reversal ! call rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, & rtau, xq, isq, cz1, cz2, .true. ) ! call rotate_epmat ( cz1, cz2, xq, nqc ) ! endif ! call flush(stdout) ! opteron only ! enddo iq_first = iq_first + nq ! enddo ! endif ! if ( epbread .or. epbwrite ) then ! ! read/write the e-ph matrix elements and other info in the Bloch representation ! in .epb files (one for each pool) ! call set_ndnmbr ( 0, mypool, 1, npool, filelab) tempfile = trim(tmp_dir) // trim(prefix) // '.epb' // filelab open (iuepb, file = tempfile, form = 'unformatted') ! if (epbwrite) then write(6,'(/5x,"Writing epmatq on .epb files"/)') write (iuepb) nqc, xqc, et, dynq, epmatq endif if (epbread) then write(6,'(/5x,"Reading epmatq from .epb files"/)') read (iuepb) nqc, xqc, et, dynq, epmatq endif ! close (iuepb) endif #ifdef __PARA call mp_barrier() #endif ! if (nqc.ne.nq1*nq2*nq3) & call errore ('elphon_shuffle_wrap','nqc .ne. nq1*nq2*nq3',nqc) wqlist = float(1)/float(nqc) ! ! now dynq is the cartesian dyn mat ( NOT divided by the masses) ! and epmatq is the epmat in cartesian representation (rotation in elphon_shuffle) ! ! modify default (in openfilq.f90 we did not know nbnd yet) ! if (nbndsub.eq.0) nbndsub = nbnd ! ! ! endif .not. epwread endif ! nqc = nq1*nq2*nq3 ! just in case we skipped the calculation of mat. el. ! ! free up some memory ! if ( associated (evq) ) nullify (evq) if ( allocated (evc) ) deallocate (evc) if ( associated (igkq) ) nullify (igkq) if ( allocated (igk) ) deallocate (igk) if ( allocated (dvpsi)) deallocate (dvpsi) if ( allocated (dpsi) ) deallocate (dpsi) ! ! the electron-phonon wannier interpolation ! call ephwann_shuffle ( nqc, xqc ) ! ! the calculation of the self energy ! if (elecselfen) call selfen_elec ! if (phonselfen) call selfen_phon ! if (specfun) call spectral_func ! 5 format (8x,"q(",i5," ) = (",3f12.7," )") ! return end subroutine elphon_shuffle_wrap ! !--------------------------------------------------------------------------- subroutine irotate ( x, s, sx) !--------------------------------------------------------------------------- ! ! a simple symmetry operation in crystal coordinates ( s is integer!) ! USE kinds, only : DP implicit none real(kind=DP) x(3), sx(3) integer :: s(3,3),i ! do i = 1, 3 sx (i) = float(s (i, 1)) * x (1) & + float(s (i, 2)) * x (2) & + float(s (i, 3)) * x (3) enddo ! return end !--------------------------------------------------------------------------- logical function eqvect_strict (x, y) !----------------------------------------------------------------------- ! ! This function test if two tridimensional vectors are equal ! USE kinds implicit none real(kind=DP) :: x (3), y (3) ! input: input vector ! input: second input vector real(kind=DP) :: accep ! acceptance parameter parameter (accep = 1.0d-5) ! eqvect_strict = abs( x(1)-y(1) ).lt.accep .and. & abs( x(2)-y(2) ).lt.accep .and. & abs( x(3)-y(3) ).lt.accep return end function eqvect_strict !