diff --git a/test-suite/not_epw_comp/do_phonon.f90 b/test-suite/not_epw_comp/do_phonon.f90 deleted file mode 100644 index b48496c00..000000000 --- a/test-suite/not_epw_comp/do_phonon.f90 +++ /dev/null @@ -1,142 +0,0 @@ -! -! Copyright (C) 2001-2013 Quantum ESPRESSO group -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! -!----------------------------------------------------------------------- -SUBROUTINE do_phonon(auxdyn) - !----------------------------------------------------------------------- - ! - ! ... This is the main driver of the phonon code. - ! ... It assumes that the preparatory stuff has been already done. - ! ... When the code calls this routine it has already read input - ! ... decided which irreducible representations have to be calculated - ! ... and it has set the variables that decide which work this routine - ! ... will do. The parallel stuff has been already setup by the calling - ! ... codes. This routine makes the two loops over - ! ... the q points and the irreps and does only the calculations - ! ... that have been decided by the driver routine. - ! ... At a generic q, if necessary it recalculates the band structure - ! ... calling pwscf again. - ! ... Then it can calculate the response to an atomic displacement, - ! ... the dynamical matrix at that q, and the electron-phonon - ! ... interaction at that q. At q=0 it can calculate the linear response - ! ... to an electric field perturbation and hence the dielectric - ! ... constant, the Born effective charges and the polarizability - ! ... at imaginary frequencies. - ! ... At q=0, from the second order response to an electric field, - ! ... it can calculate also the electro-optic and the raman tensors. - ! - - USE disp, ONLY : nqs - USE control_ph, ONLY : epsil, trans, qplot, only_init, & - only_wfc, rec_code, where_rec - USE el_phon, ONLY : elph, elph_mat, elph_simple, elph_epa - ! - ! YAMBO > - USE YAMBO, ONLY : elph_yambo - ! YAMBO < - ! - USE elph_tetra_mod, ONLY : elph_tetra, elph_tetra_lambda, elph_tetra_gamma - USE elph_scdft_mod, ONLY : elph_scdft - - IMPLICIT NONE - ! - CHARACTER (LEN=256), INTENT(IN) :: auxdyn - INTEGER :: iq - LOGICAL :: do_band, do_iq, setup_pw - ! - DO iq = 1, nqs - ! - CALL prepare_q(auxdyn, do_band, do_iq, setup_pw, iq) - ! - ! If this q is not done in this run, cycle - ! - IF (.NOT.do_iq) CYCLE - ! - ! If necessary the bands are recalculated - ! - IF (setup_pw) CALL run_nscf(do_band, iq) - ! - ! If only_wfc=.TRUE. the code computes only the wavefunctions - ! - IF (only_wfc) THEN - where_rec='only_wfc' - rec_code=-1000 - GOTO 100 - ENDIF - ! - ! Initialize the quantities which do not depend on - ! the linear response of the system - ! - CALL initialize_ph() - ! - ! electric field perturbation - ! - IF (epsil) CALL phescf() - ! - ! IF only_init is .true. the code computes only the - ! initialization parts. - ! - IF (only_init) THEN - where_rec='only_init' - rec_code=-1000 - GOTO 100 - ENDIF - ! - ! phonon perturbation - ! - IF ( trans ) THEN - ! - CALL phqscf() - CALL dynmatrix_new(iq) - ! - END IF - ! - CALL rotate_dvscf_star(iq) - ! - ! electron-phonon interaction - ! - IF ( elph ) THEN - ! - IF ( .NOT. trans ) THEN - ! - CALL dvanqq() - IF ( elph_mat ) THEN - CALL ep_matrix_element_wannier() - ELSE - CALL elphon() - END IF - ! - END IF - ! - !IF ( elph_mat ) THEN - ! CALL elphsum_wannier(iq) - !ELSEIF( elph_simple ) THEN - ! CALL elphsum_simple() - !ELSEIF( elph_epa ) THEN - ! CALL elphfil_epa(iq) - !ELSEIF( elph_yambo ) THEN - ! CALL elph_yambo_eval_and_IO() - !ELSEIF(elph_tetra == 1) THEN - ! CALL elph_tetra_lambda() - !ELSEIF(elph_tetra == 2) THEN - ! CALL elph_tetra_gamma() - !ELSEIF(elph_tetra == 3) THEN - ! CALL elph_scdft() - !ELSE - !CALL elphsum() - CALL elphsum2() - !END IF - ! - END IF - ! - ! ... cleanup of the variables for the next q point - ! -100 CALL clean_pw_ph(iq) - ! - END DO - -END SUBROUTINE do_phonon diff --git a/test-suite/not_epw_comp/elphon.f90 b/test-suite/not_epw_comp/elphon.f90 deleted file mode 100644 index 72edb06a0..000000000 --- a/test-suite/not_epw_comp/elphon.f90 +++ /dev/null @@ -1,1930 +0,0 @@ -! -! Copyright (C) 2001-2018 Quantum ESPRESSO group -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! -! -!----------------------------------------------------------------------- -SUBROUTINE elphon() - !----------------------------------------------------------------------- - ! - ! Electron-phonon calculation from data saved in fildvscf - ! - USE kinds, ONLY : DP - USE constants, ONLY : amu_ry, RY_TO_THZ, RY_TO_CMM1 - USE cell_base, ONLY : celldm, omega, ibrav, at, bg - USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, amass - USE gvecs, ONLY: doublegrid - USE fft_base, ONLY : dfftp, dffts - USE fft_interfaces, ONLY : fft_interpolate - USE noncollin_module, ONLY : nspin_mag, noncolin, m_loc - USE lsda_mod, ONLY : nspin - USE uspp, ONLY: okvan - USE paw_variables, ONLY : okpaw - USE el_phon, ONLY : done_elph - USE dynmat, ONLY : dyn, w2 - USE modes, ONLY : npert, nirr, u - USE uspp_param, ONLY : nhm - USE control_ph, ONLY : trans, xmldyn - USE output, ONLY : fildyn,fildvscf - USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, & - read_dyn_mat, read_dyn_mat_tail - USE units_ph, ONLY : iudyn, lrdrho, iudvscf, iuint3paw, lint3paw - USE dfile_star, ONLY : dvscf_star - USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, root_bgrp - USE mp, ONLY : mp_bcast - USE io_global, ONLY: stdout - - USE lrus, ONLY : int3, int3_nc, int3_paw - USE qpoint, ONLY : xq - ! - IMPLICIT NONE - ! - INTEGER :: irr, imode0, ipert, is, npe - ! counter on the representations - ! counter on the modes - ! the change of Vscf due to perturbations - INTEGER :: i,j - COMPLEX(DP), POINTER :: dvscfin(:,:,:), dvscfins (:,:,:) - COMPLEX(DP), allocatable :: phip (:, :, :, :) - - INTEGER :: ntyp_, nat_, ibrav_, nspin_mag_, mu, nu, na, nb, nta, ntb, nqs_ - REAL(DP) :: celldm_(6), w1 - CHARACTER(LEN=3) :: atm(ntyp) - - CALL start_clock ('elphon') - - if(dvscf_star%basis.eq.'cartesian') then - write(stdout,*) 'Setting patterns to identity' - u=(0.d0,0.d0) - do irr=1,3*nat - u(irr,irr)=(1.d0,0.d0) - enddo - endif - ! - ! read Delta Vscf and calculate electron-phonon coefficients - ! - imode0 = 0 - WRITE (6, '(5x,a)') "Reading dVscf from file "//trim(fildvscf) - DO irr = 1, nirr - npe=npert(irr) - ALLOCATE (dvscfin (dfftp%nnr, nspin_mag , npe) ) - IF (okvan) THEN - ALLOCATE (int3 ( nhm, nhm, nat, nspin_mag, npe)) - IF (okpaw) ALLOCATE (int3_paw (nhm, nhm, nat, nspin_mag, npe)) - IF (noncolin) ALLOCATE(int3_nc( nhm, nhm, nat, nspin, npe)) - ENDIF - DO ipert = 1, npe - CALL davcio_drho ( dvscfin(1,1,ipert), lrdrho, iudvscf, & - imode0 + ipert, -1 ) - IF (okpaw .AND. me_bgrp==0) & - CALL davcio( int3_paw(:,:,:,:,ipert), lint3paw, & - iuint3paw, imode0 + ipert, - 1 ) - END DO - IF (okpaw) CALL mp_bcast(int3_paw, root_bgrp, intra_bgrp_comm) - IF (doublegrid) THEN - ALLOCATE (dvscfins (dffts%nnr, nspin_mag , npert(irr)) ) - DO is = 1, nspin_mag - DO ipert = 1, npe - CALL fft_interpolate (dfftp, dvscfin(:,is,ipert), dffts, dvscfins(:,is,ipert)) - ENDDO - ENDDO - ELSE - dvscfins => dvscfin - ENDIF - CALL newdq (dvscfin, npert(irr)) - CALL elphel (irr, npert (irr), imode0, dvscfins) - ! - imode0 = imode0 + npe - IF (doublegrid) DEALLOCATE (dvscfins) - DEALLOCATE (dvscfin) - IF (okvan) THEN - DEALLOCATE (int3) - IF (okpaw) DEALLOCATE (int3_paw) - IF (noncolin) DEALLOCATE(int3_nc) - ENDIF - ENDDO - ! - ! now read the eigenvalues and eigenvectors of the dynamical matrix - ! calculated in a previous run - ! - IF (.NOT.trans) THEN - IF (.NOT. xmldyn) THEN - WRITE (6, '(5x,a)') "Reading dynamics matrix from file "//trim(fildyn) - CALL readmat (iudyn, ibrav, celldm, nat, ntyp, & - ityp, omega, amass, tau, xq, w2, dyn) - ELSE - allocate( phip(3,3,nat,nat) ) - CALL read_dyn_mat_param(fildyn, ntyp_, nat_) - IF ( ntyp_ /= ntyp .OR. nat_ /= nat ) & - CALL errore('elphon','uncorrect nat or ntyp',1) - - CALL read_dyn_mat_header(ntyp, nat, ibrav_, nspin_mag_, & - celldm_, at, bg, omega, atm, amass, tau, ityp, & - m_loc, nqs_) - - IF (ibrav_.NE.ibrav .OR. ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5 & - .OR. (nspin_mag_ /= nspin_mag ) ) CALL errore ('elphon', & - 'inconsistent data', 1) - - CALL read_dyn_mat(nat,1,xq,phip) - ! - ! Diagonalize the dynamical matrix - ! - - - DO i=1,3 - do na=1,nat - nta = ityp (na) - mu=3*(na-1)+i - do j=1,3 - do nb=1,nat - nu=3*(nb-1)+j - ntb = ityp (nb) - dyn (mu, nu) = phip (i, j, na, nb) / & - sqrt( amass(nta)*amass(ntb))/amu_ry - enddo - enddo - enddo - enddo - - ! - CALL cdiagh (3 * nat, dyn, 3 * nat, w2, dyn) - ! - ! divide by sqrt(mass) to get displacements - ! - DO nu = 1, 3 * nat - DO mu = 1, 3 * nat - na = (mu - 1) / 3 + 1 - dyn (mu, nu) = dyn (mu, nu) / SQRT ( amu_ry * amass (ityp (na) ) ) - ENDDO - ENDDO - - CALL read_dyn_mat_tail(nat) - - deallocate( phip ) - ENDIF - ! - ! Write phonon frequency to stdout - ! - WRITE( stdout, 8000) (xq (i), i = 1, 3) - ! - DO nu = 1, 3 * nat - w1 = SQRT( ABS( w2(nu) ) ) - if (w2(nu) < 0.d0) w1 = - w1 - WRITE( stdout, 8010) nu, w1 * RY_TO_THZ, w1 * RY_TO_CMM1 - ENDDO - ! - WRITE( stdout, '(1x,74("*"))') - ! - ENDIF ! .NOT. trans - ! - CALL stop_clock ('elphon') - ! -8000 FORMAT(/,5x,'Diagonalizing the dynamical matrix', & - & //,5x,'q = ( ',3f14.9,' ) ',//,1x,74('*')) -8010 FORMAT (5x,'freq (',i5,') =',f15.6,' [THz] =',f15.6,' [cm-1]') - ! - RETURN -END SUBROUTINE elphon -! -!----------------------------------------------------------------------- -SUBROUTINE readmat (iudyn, ibrav, celldm, nat, ntyp, ityp, omega, & - amass, tau, q, w2, dyn) - !----------------------------------------------------------------------- - ! - USE kinds, ONLY : DP - USE constants, ONLY : amu_ry - IMPLICIT NONE - ! Input - INTEGER :: iudyn, ibrav, nat, ntyp, ityp (nat) - REAL(DP) :: celldm (6), amass (ntyp), tau (3, nat), q (3), & - omega - ! output - REAL(DP) :: w2 (3 * nat) - COMPLEX(DP) :: dyn (3 * nat, 3 * nat) - ! local (control variables) - INTEGER :: ntyp_, nat_, ibrav_, ityp_ - REAL(DP) :: celldm_ (6), amass_, tau_ (3), q_ (3) - ! local - REAL(DP) :: dynr (2, 3, nat, 3, nat) - CHARACTER(len=80) :: line - CHARACTER(len=3) :: atm - INTEGER :: nt, na, nb, naa, nbb, nu, mu, i, j - ! - ! - REWIND (iudyn) - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (iudyn, * ) ntyp_, nat_, ibrav_, celldm_ - IF ( ntyp.NE.ntyp_ .OR. nat.NE.nat_ .OR.ibrav_.NE.ibrav .OR. & - ABS ( celldm_ (1) - celldm (1) ) > 1.0d-5) & - CALL errore ('readmat', 'inconsistent data', 1) - IF ( ibrav_ == 0 ) THEN - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - END IF - DO nt = 1, ntyp - READ (iudyn, * ) i, atm, amass_ - IF ( nt.NE.i .OR. ABS (amass_ - amu_ry*amass (nt) ) > 1.0d-5) & - CALL errore ( 'readmat', 'inconsistent data', 1 + nt) - ENDDO - DO na = 1, nat - READ (iudyn, * ) i, ityp_, tau_ - IF (na.NE.i.OR.ityp_.NE.ityp (na) ) CALL errore ('readmat', & - 'inconsistent data', 10 + na) - ENDDO - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (iudyn, '(a)') line - READ (line (11:80), * ) (q_ (i), i = 1, 3) - READ (iudyn, '(a)') line - DO na = 1, nat - DO nb = 1, nat - READ (iudyn, * ) naa, nbb - IF (na.NE.naa.OR.nb.NE.nbb) CALL errore ('readmat', 'error reading & - &file', nb) - READ (iudyn, * ) ( (dynr (1, i, na, j, nb), dynr (2, i, na, j, nb) & - , j = 1, 3), i = 1, 3) - ENDDO - ENDDO - ! - ! divide the dynamical matrix by the (input) masses (in amu) - ! - DO nb = 1, nat - DO j = 1, 3 - DO na = 1, nat - DO i = 1, 3 - dynr (1, i, na, j, nb) = dynr (1, i, na, j, nb) / SQRT (amass ( & - ityp (na) ) * amass (ityp (nb) ) ) / amu_ry - dynr (2, i, na, j, nb) = dynr (2, i, na, j, nb) / SQRT (amass ( & - ityp (na) ) * amass (ityp (nb) ) ) / amu_ry - ENDDO - ENDDO - ENDDO - ENDDO - ! - ! solve the eigenvalue problem. - ! NOTA BENE: eigenvectors are overwritten on dyn - ! - CALL cdiagh (3 * nat, dynr, 3 * nat, w2, dyn) - ! - ! divide by sqrt(mass) to get displacements - ! - DO nu = 1, 3 * nat - DO mu = 1, 3 * nat - na = (mu - 1) / 3 + 1 - dyn (mu, nu) = dyn (mu, nu) / SQRT ( amu_ry * amass (ityp (na) ) ) - ENDDO - ENDDO - ! - ! - RETURN -END SUBROUTINE readmat -! -!----------------------------------------------------------------------- -SUBROUTINE elphel (irr, npe, imode0, dvscfins) - !----------------------------------------------------------------------- - ! - ! Calculation of the electron-phonon matrix elements el_ph_mat - ! <\psi(k+q)|dV_{SCF}/du^q_{i a}|\psi(k)> - ! Original routine written by Francesco Mauri - ! Modified by A. Floris and I. Timrov to include Hubbard U (01.10.2018) - ! - USE kinds, ONLY : DP - USE fft_base, ONLY : dffts - USE ions_base, ONLY : nat, ityp - USE control_flags, ONLY : iverbosity - USE wavefunctions, ONLY : evc - USE buffers, ONLY : get_buffer - USE klist, ONLY : xk, ngk, igk_k - USE lsda_mod, ONLY : lsda, current_spin, isk, nspin - USE noncollin_module, ONLY : noncolin, npol, nspin_mag - USE wvfct, ONLY : nbnd, npwx - USE buffers, ONLY : get_buffer - USE uspp, ONLY : vkb - USE el_phon, ONLY : el_ph_mat, el_ph_mat_rec, el_ph_mat_rec_col, & - comp_elph, done_elph, elph_nbnd_min, elph_nbnd_max - USE modes, ONLY : u, nmodes - USE units_ph, ONLY : iubar, lrbar, iundnsscf - USE units_lr, ONLY : iuwfc, lrwfc - USE control_ph, ONLY : trans, current_iq - USE ph_restart, ONLY : ph_writefile - USE spin_orb, ONLY : domag - USE mp_bands, ONLY : intra_bgrp_comm, ntask_groups - USE mp_pools, ONLY : npool - USE mp, ONLY : mp_sum, mp_bcast - USE mp_world, ONLY : world_comm - USE elph_tetra_mod, ONLY : elph_tetra - USE eqv, ONLY : dvpsi, evq - USE qpoint, ONLY : nksq, ikks, ikqs, nksqtot - USE control_lr, ONLY : lgamma - USE fft_helper_subroutines - USE ldaU, ONLY : lda_plus_u, Hubbard_lmax - USE ldaU_ph, ONLY : dnsscf_all_modes, dnsscf - USE io_global, ONLY : ionode, ionode_id - USE lrus, ONLY : becp1 - USE phus, ONLY : alphap - - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: irr, npe, imode0 - COMPLEX(DP), INTENT(IN) :: dvscfins (dffts%nnr, nspin_mag, npe) - ! LOCAL variables - INTEGER :: npw, npwq, nrec, ik, ikk, ikq, ipert, mode, ibnd, jbnd, ir, ig, & - ipol, ios, ierr - COMPLEX(DP) , ALLOCATABLE :: aux1 (:,:), elphmat (:,:,:), tg_dv(:,:), & - tg_psic(:,:), aux2(:,:) - INTEGER :: v_siz, incr - COMPLEX(DP), EXTERNAL :: zdotc - integer :: ibnd_fst, ibnd_lst - ! - if(elph_tetra == 0) then - ibnd_fst = 1 - ibnd_lst = nbnd - else - ibnd_fst = elph_nbnd_min - ibnd_lst = elph_nbnd_max - end if - ! - IF (.NOT. comp_elph(irr) .OR. done_elph(irr)) RETURN - - ALLOCATE (aux1 (dffts%nnr, npol)) - ALLOCATE (elphmat ( nbnd , nbnd , npe)) - ALLOCATE( el_ph_mat_rec (nbnd,nbnd,nksq,npe) ) - el_ph_mat_rec=(0.0_DP,0.0_DP) - ALLOCATE (aux2(npwx*npol, nbnd)) - incr=1 - IF ( dffts%has_task_groups ) THEN - ! - v_siz = dffts%nnr_tg - ALLOCATE( tg_dv ( v_siz, nspin_mag ) ) - ALLOCATE( tg_psic( v_siz, npol ) ) - incr = fftx_ntgrp(dffts) - ! - ENDIF - ! - ! DFPT+U case - ! - IF (lda_plus_u) THEN - ! - ! Allocate and re-read dnsscf_all_modes from file - ! - ALLOCATE (dnsscf_all_modes(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, nmodes)) - dnsscf_all_modes = (0.d0, 0.d0) - IF (ionode) READ(iundnsscf,*) dnsscf_all_modes - CALL mp_bcast(dnsscf_all_modes, ionode_id, world_comm) - REWIND(iundnsscf) - ! - ! Check whether the re-read is correct - ! - IF (iverbosity==1) CALL elphel_read_dnsscf_check() - ! - ! Allocate dnsscf - ! - ALLOCATE (dnsscf(2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, npe)) - dnsscf = (0.d0, 0.d0) - ! - ENDIF - ! - ! Start the loops over the k-points - ! - DO ik = 1, nksq - ! - ! ik = counter of k-points with vector k - ! ikk= index of k-point with vector k - ! ikq= index of k-point with vector k+q - ! k and k+q are alternated if q!=0, are the same if q=0 - ! - ikk = ikks(ik) - ikq = ikqs(ik) - IF (lsda) current_spin = isk (ikk) - npw = ngk(ikk) - npwq= ngk(ikq) - ! - CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) - ! - ! read unperturbed wavefuctions psi(k) and psi(k+q) - ! - IF (nksq.GT.1) THEN - IF (lgamma) THEN - CALL get_buffer(evc, lrwfc, iuwfc, ikk) - ELSE - CALL get_buffer (evc, lrwfc, iuwfc, ikk) - CALL get_buffer (evq, lrwfc, iuwfc, ikq) - ENDIF - ENDIF - ! - DO ipert = 1, npe - nrec = (ipert - 1) * nksq + ik - ! - ! dvbare_q*psi_kpoint is read from file (if available) or recalculated - ! - IF (trans) THEN - CALL get_buffer (dvpsi, lrbar, iubar, nrec) - ELSE - mode = imode0 + ipert - ! FIXME: .false. or .true. ??? - CALL dvqpsi_us (ik, u (1, mode), .FALSE., becp1, alphap) - ! - ! DFPT+U: calculate the bare derivative of the Hubbard potential in el-ph - ! - IF (lda_plus_u) CALL dvqhub_barepsi_us (ik, u(1,mode)) - ! - ENDIF - ! - ! calculate dvscf_q*psi_k - ! - IF ( dffts%has_task_groups ) THEN - IF (noncolin) THEN - CALL tg_cgather( dffts, dvscfins(:,1,ipert), tg_dv(:,1)) - IF (domag) THEN - DO ipol=2,4 - CALL tg_cgather( dffts, dvscfins(:,ipol,ipert), tg_dv(:,ipol)) - ENDDO - ENDIF - ELSE - CALL tg_cgather( dffts, dvscfins(:,current_spin,ipert), tg_dv(:,1)) - ENDIF - ENDIF - aux2=(0.0_DP,0.0_DP) - DO ibnd = ibnd_fst, ibnd_lst, incr - IF ( dffts%has_task_groups ) THEN - CALL cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd ) - CALL apply_dpot(v_siz, tg_psic, tg_dv, 1) - CALL cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd) - ELSE - CALL cft_wave (ik, evc(1, ibnd), aux1, +1) - CALL apply_dpot(dffts%nnr, aux1, dvscfins(1,1,ipert), current_spin) - CALL cft_wave (ik, aux2(1, ibnd), aux1, -1) - ENDIF - ENDDO - dvpsi=dvpsi+aux2 - ! - CALL adddvscf (ipert, ik) - ! - ! DFPT+U: add to dvpsi the scf part of the perturbed Hubbard potential - ! - IF (lda_plus_u) THEN - dnsscf(:,:,:,:,ipert) = dnsscf_all_modes(:,:,:,:,mode) - CALL adddvhubscf (ipert, ik) - ENDIF - ! - ! calculate elphmat(j,i)= for this pertur - ! - DO ibnd = ibnd_fst, ibnd_lst - DO jbnd = ibnd_fst, ibnd_lst - elphmat (jbnd, ibnd, ipert) = zdotc (npwq, evq (1, jbnd), 1, & - dvpsi (1, ibnd), 1) - IF (noncolin) & - elphmat (jbnd, ibnd, ipert) = elphmat (jbnd, ibnd, ipert)+ & - zdotc (npwq, evq(npwx+1,jbnd),1,dvpsi(npwx+1,ibnd), 1) - ENDDO - ENDDO - ENDDO - ! - CALL mp_sum (elphmat, intra_bgrp_comm) - ! - ! save all e-ph matrix elements into el_ph_mat - ! - DO ipert = 1, npe - DO jbnd = ibnd_fst, ibnd_lst - DO ibnd = ibnd_fst, ibnd_lst - el_ph_mat (ibnd, jbnd, ik, ipert + imode0) = elphmat (ibnd, jbnd, ipert) - el_ph_mat_rec (ibnd, jbnd, ik, ipert ) = elphmat (ibnd, jbnd, ipert) - ENDDO - ENDDO - ENDDO - ENDDO - ! - done_elph(irr)=.TRUE. - if(elph_tetra == 0) then - IF (npool>1) THEN - ALLOCATE(el_ph_mat_rec_col(nbnd,nbnd,nksqtot,npe)) - CALL el_ph_collect(npe,el_ph_mat_rec,el_ph_mat_rec_col,nksqtot,nksq) - ELSE - el_ph_mat_rec_col => el_ph_mat_rec - ENDIF - CALL ph_writefile('el_phon',current_iq,irr,ierr) - IF (npool > 1) DEALLOCATE(el_ph_mat_rec_col) - end if - DEALLOCATE(el_ph_mat_rec) - ! - DEALLOCATE (elphmat) - DEALLOCATE (aux1) - DEALLOCATE (aux2) - IF ( dffts%has_task_groups ) THEN - DEALLOCATE( tg_dv ) - DEALLOCATE( tg_psic ) - ENDIF - ! - IF (lda_plus_u) THEN - DEALLOCATE (dnsscf_all_modes) - DEALLOCATE (dnsscf) - ENDIF - ! - RETURN - ! -END SUBROUTINE elphel -! -!------------------------------------------------------------------------ -SUBROUTINE elphel_read_dnsscf_check() - ! - ! DFPT+U: This subroutine checks whether dnsscf_all_modes was - ! read correctly from file. - ! - USE kinds, ONLY : DP - USE ions_base, ONLY : nat, ityp - USE modes, ONLY : u, nmodes - USE lsda_mod, ONLY : nspin - USE ldaU, ONLY : Hubbard_l, is_hubbard, Hubbard_lmax - USE ldaU_ph, ONLY : dnsscf_all_modes - USE io_global, ONLY : stdout - ! - IMPLICIT NONE - ! - COMPLEX(DP), ALLOCATABLE :: dnsscf_all_modes_cart(:,:,:,:,:) - INTEGER :: na_icart, nah, is, m1, m2, na, icart, nt, na_icar, imode - ! - ALLOCATE(dnsscf_all_modes_cart (2*Hubbard_lmax+1, 2*Hubbard_lmax+1, nspin, nat, nmodes)) - dnsscf_all_modes_cart = (0.d0, 0.d0) - ! - ! Transform dnsscf_all_modes from pattern to cartesian coordinates - ! - DO na_icart = 1, 3*nat - DO imode = 1, nmodes - DO nah = 1, nat - nt = ityp(nah) - IF (is_hubbard(nt)) THEN - DO is = 1, nspin - DO m1 = 1, 2*Hubbard_l(nt) + 1 - DO m2 = 1, 2*Hubbard_l(nt) + 1 - ! - dnsscf_all_modes_cart (m1, m2, is, nah, na_icart) = & - dnsscf_all_modes_cart (m1, m2, is, nah, na_icart) + & - dnsscf_all_modes (m1, m2, is, nah, imode) * & - CONJG(u(na_icart,imode)) - ! - ENDDO - ENDDO - ENDDO - ENDIF - ENDDO - ENDDO - ENDDO - ! - ! Write dnsscf in cartesian coordinates - ! - WRITE(stdout,*) - WRITE(stdout,*) 'DNS_SCF SYMMETRIZED IN CARTESIAN COORDINATES' - ! - DO na = 1, nat - DO icart = 1, 3 - WRITE( stdout,'(a,1x,i2,2x,a,1x,i2)') 'displaced atom L =', na, 'ipol=', icart - na_icart = 3*(na-1) + icart - DO nah = 1, nat - nt = ityp(nah) - IF (is_hubbard(nt)) THEN - DO is = 1, nspin - WRITE(stdout,'(a,1x,i2,2x,a,1x,i2)') ' Hubbard atom', nah, 'spin', is - DO m1 = 1, 2*Hubbard_l(nt) + 1 - WRITE(stdout,'(14(f15.10,1x))') dnsscf_all_modes_cart (m1,:,is,nah,na_icart) - ENDDO - ENDDO - ENDIF - ENDDO - ENDDO - ENDDO - WRITE(stdout,*) - ! - DEALLOCATE(dnsscf_all_modes_cart) - ! - RETURN - ! -END SUBROUTINE elphel_read_dnsscf_check -!------------------------------------------------------------------------ - -!------------------------------------------------------------------------ -SUBROUTINE elphsum ( ) - !----------------------------------------------------------------------- - ! - ! Sum over BZ of the electron-phonon matrix elements el_ph_mat - ! Original routine written by Francesco Mauri, modified by PG - ! New version by Malgorzata Wierzbowska - ! - USE kinds, ONLY : DP - USE constants, ONLY : pi, rytoev, ry_to_cmm1, ry_to_ghz, degspin - USE ions_base, ONLY : nat, ityp, tau - USE cell_base, ONLY : at, bg - USE lsda_mod, ONLY: isk, nspin - USE klist, ONLY: nks, nkstot, xk, wk, nelec - USE start_k, ONLY: nk1, nk2, nk3 - USE symm_base, ONLY: s, irt, nsym, invs - USE noncollin_module, ONLY: nspin_lsda, nspin_mag - USE wvfct, ONLY: nbnd, et - USE parameters, ONLY : npk - USE el_phon, ONLY : el_ph_mat, done_elph, el_ph_nsigma, el_ph_ngauss, & - el_ph_sigma - USE modes, ONLY : u, nirr - USE dynmat, ONLY : dyn, w2 - USE io_global, ONLY : stdout, ionode, ionode_id - USE mp_pools, ONLY : my_pool_id, npool, kunit - USE mp_images, ONLY : intra_image_comm - USE mp, ONLY : mp_bcast - USE control_ph, ONLY : tmp_dir_phq, xmldyn, current_iq - USE save_ph, ONLY : tmp_dir_save - USE io_files, ONLY : prefix, tmp_dir, seqopn, create_directory - ! - USE lr_symm_base, ONLY : minus_q, nsymq, rtau - USE qpoint, ONLY : xq, nksq - USE control_lr, ONLY : lgamma - ! - IMPLICIT NONE - ! epsw = 20 cm^-1, in Ry - REAL(DP), PARAMETER :: epsw = 20.d0 / ry_to_cmm1 - REAL(DP), PARAMETER :: eps = 1.0d-6 - ! - INTEGER :: iuna2Fsave = 40 - ! - REAL(DP), allocatable :: gam(:,:), lamb(:,:) - ! - ! Quantities ending with "fit" are relative to the "dense" grid - ! - REAL(DP), allocatable :: xkfit(:,:) - REAL(DP), allocatable, target :: etfit(:,:), wkfit(:) - INTEGER :: nksfit, nk1fit, nk2fit, nk3fit, nkfit, nksfit_real - INTEGER, allocatable :: eqkfit(:), eqqfit(:), sfit(:) - ! - integer :: nq, isq (48), imq - ! nq : degeneracy of the star of q - ! isq: index of q in the star of a given sym.op. - ! imq: index of -q in the star of q (0 if not present) - real(DP) :: sxq (3, 48) - ! list of vectors in the star of q - ! - ! workspace used for symmetrisation - ! - COMPLEX(DP), allocatable :: g1(:,:,:), g2(:,:,:), g0(:,:), gf(:,:,:) - COMPLEX(DP), allocatable :: point(:), noint(:), ctemp(:) - COMPLEX(DP) :: dyn22(3*nat,3*nat) - ! - INTEGER :: ik, ikk, ikq, isig, ibnd, jbnd, ipert, jpert, nu, mu, & - vu, ngauss1, iuelph, ios, i,k,j, ii, jj - INTEGER :: nkBZ, nti, ntj, ntk, nkr, itemp1, itemp2, nn, & - qx,qy,qz,iq,jq,kq - INTEGER, ALLOCATABLE :: eqBZ(:), sBZ(:) - REAL(DP) :: weight, wqa, w0g1, w0g2, degauss1, dosef, & - ef1, lambda, gamma - REAL(DP), ALLOCATABLE :: deg(:), effit(:), dosfit(:) - REAL(DP) :: etk, etq - REAL(DP), EXTERNAL :: dos_ef, efermig, w0gauss - character(len=80) :: name - LOGICAL :: exst, xmldyn_save - ! - COMPLEX(DP) :: el_ph_sum (3*nat,3*nat) - - COMPLEX(DP), POINTER :: el_ph_mat_collect(:,:,:,:) - REAL(DP), ALLOCATABLE :: xk_collect(:,:) - REAL(DP), POINTER :: wkfit_dist(:), etfit_dist(:,:) - INTEGER :: nksfit_dist, rest, kunit_save - INTEGER :: nks_real, ispin, nksqtot, irr - CHARACTER(LEN=256) :: elph_dir - CHARACTER(LEN=6) :: int_to_char - ! - ! - ! - ! If the electron phonon matrix elements have not been calculated for - ! all representations this routine exit - ! - DO irr=1,nirr - IF (.NOT.done_elph(irr)) RETURN - ENDDO - elph_dir='elph_dir/' - IF (ionode) INQUIRE(file=TRIM(elph_dir), EXIST=exst) - CALL mp_bcast(exst, ionode_id, intra_image_comm) - IF (.NOT.exst) CALL create_directory( elph_dir ) - WRITE (6, '(5x,"electron-phonon interaction ..."/)') - ngauss1 = 0 - - ALLOCATE(xk_collect(3,nkstot)) - - ALLOCATE(deg(el_ph_nsigma)) - ALLOCATE(effit(el_ph_nsigma)) - ALLOCATE(dosfit(el_ph_nsigma)) - - IF (npool==1) THEN -! -! no pool, just copy old variables on the new ones -! - nksqtot=nksq - xk_collect(:,1:nks) = xk(:,1:nks) - el_ph_mat_collect => el_ph_mat - ELSE -! -! pools, allocate new variables and collect the results. All the rest -! remain unchanged. -! - IF (lgamma) THEN - nksqtot=nkstot - ELSE - nksqtot=nkstot/2 - ENDIF - CALL poolcollect(3, nks, xk, nkstot, xk_collect) - ALLOCATE(el_ph_mat_collect(nbnd,nbnd,nksqtot,3*nat)) - ! FIXME: this routine should be replaced by a generic routine - CALL el_ph_collect(3*nat,el_ph_mat,el_ph_mat_collect,nksqtot,nksq) - ENDIF - ! - ! read eigenvalues for the dense grid - ! FIXME: this should be done from the xml file, not from a specialized file - ! parallel case: only first node reads - ! - IF ( ionode ) THEN - tmp_dir=tmp_dir_save - CALL seqopn( iuna2Fsave, 'a2Fsave', 'FORMATTED', exst ) - tmp_dir=tmp_dir_phq - READ(iuna2Fsave,*) ibnd, nksfit - END IF - ! - CALL mp_bcast (ibnd, ionode_id, intra_image_comm) - CALL mp_bcast (nksfit, ionode_id, intra_image_comm) - if ( ibnd /= nbnd ) call errore('elphsum','wrong file read',iuna2Fsave) - allocate (etfit(nbnd,nksfit), xkfit(3,nksfit), wkfit(nksfit)) - ! - IF ( ionode ) THEN - READ(iuna2Fsave,*) etfit - READ(iuna2Fsave,*) ((xkfit(i,ik), i=1,3), ik=1,nksfit) - READ(iuna2Fsave,*) wkfit - READ(iuna2Fsave,*) nk1fit, nk2fit, nk3fit - CLOSE( UNIT = iuna2Fsave, STATUS = 'KEEP' ) - END IF - ! - ! broadcast all variables read - ! - CALL mp_bcast (etfit, ionode_id, intra_image_comm) - CALL mp_bcast (xkfit, ionode_id, intra_image_comm) - CALL mp_bcast (wkfit, ionode_id, intra_image_comm) - CALL mp_bcast (nk1fit, ionode_id, intra_image_comm) - CALL mp_bcast (nk2fit, ionode_id, intra_image_comm) - CALL mp_bcast (nk3fit, ionode_id, intra_image_comm) - ! - nkfit=nk1fit*nk2fit*nk3fit - ! - ! efermig and dos_ef require scattered points and eigenvalues - ! isk is neither read nor used. phonon with two Fermi energies is - ! not yet implemented. - ! - nksfit_dist = ( nksfit / npool ) - rest = ( nksfit - nksfit_dist * npool ) - IF ( ( my_pool_id + 1 ) <= rest ) nksfit_dist = nksfit_dist + 1 - kunit_save=kunit - kunit=1 - -#if defined(__MPI) - ALLOCATE(etfit_dist(nbnd,nksfit_dist)) - ALLOCATE(wkfit_dist(nksfit_dist)) - CALL poolscatter( 1, nksfit, wkfit, nksfit_dist, wkfit_dist ) - CALL poolscatter( nbnd, nksfit, etfit, nksfit_dist, etfit_dist ) -#else - wkfit_dist => wkfit - etfit_dist => etfit -#endif - ! - do isig=1,el_ph_nsigma - ! - ! recalculate Ef = effit and DOS at Ef N(Ef) = dosfit using dense grid - ! for value "deg" of gaussian broadening - ! - deg(isig) = isig * el_ph_sigma - ! - effit(isig) = efermig & - ( etfit_dist, nbnd, nksfit_dist, nelec, wkfit_dist, & - deg(isig), ngauss1, 0, isk) - dosfit(isig) = dos_ef ( ngauss1, deg(isig), effit(isig), etfit_dist, & - wkfit_dist, nksfit_dist, nbnd) / 2.0d0 - enddo -#if defined(__MPI) - DEALLOCATE(etfit_dist) - DEALLOCATE(wkfit_dist) -#endif - kunit=kunit_save - allocate (eqkfit(nkfit), eqqfit(nkfit), sfit(nkfit)) - ! - ! map k-points in the IBZ to k-points in the complete uniform grid - ! - nksfit_real=nksfit/nspin_lsda - call lint ( nsym, s, .true., at, bg, npk, 0,0,0, & - nk1fit,nk2fit,nk3fit, nksfit_real, xkfit, 1, nkfit, eqkfit, sfit) - deallocate (sfit, xkfit, wkfit) - ! - ! find epsilon(k+q) in the dense grid - ! - call cryst_to_cart (1, xq, at, -1) - qx = nint(nk1fit*xq(1)) - if (abs(qx-nk1fit*xq(1)) > eps) & - call errore('elphsum','q is not a vector in the dense grid',1) - if (qx < 0) qx = qx + nk1fit - if (qx > nk1fit) qx = qx - nk1fit - qy = nint(nk2fit*xq(2)) - if (abs(qy-nk2fit*xq(2)) > eps) & - call errore('elphsum','q is not a vector in the dense grid',2) - if (qy < 0) qy = qy + nk2fit - if (qy > nk2fit) qy = qy - nk2fit - qz = nint(nk3fit*xq(3)) - if (abs(qz-nk3fit*xq(3)) > eps) & - call errore('elphsum','q is not a vector in the dense grid',3) - if (qz < 0) qz = qz + nk3fit - if (qz > nk3fit) qz = qz - nk3fit - call cryst_to_cart (1, xq, bg, 1) - ! - eqqfit(:) = 0 - do i=1,nk1fit - do j=1,nk2fit - do k=1,nk3fit - ik = k-1 + (j-1)*nk3fit + (i-1)*nk2fit*nk3fit + 1 - iq = i+qx - if (iq > nk1fit) iq = iq - nk1fit - jq = j+qy - if (jq > nk2fit) jq = jq - nk2fit - kq = k+qz - if (kq > nk3fit) kq = kq - nk3fit - nn = (kq-1)+(jq-1)*nk3fit+(iq-1)*nk2fit*nk3fit + 1 - eqqfit(ik) = eqkfit(nn) - enddo - enddo - enddo - ! - ! calculate the electron-phonon coefficient using the dense grid - ! - nti = nk1fit/nk1 - ntj = nk2fit/nk2 - ntk = nk3fit/nk3 - nkBZ = nk1*nk2*nk3 - allocate (eqBZ(nkBZ), sBZ(nkBZ)) - ! - nks_real=nkstot/nspin_lsda - IF ( lgamma ) THEN - call lint ( nsymq, s, minus_q, at, bg, npk, 0,0,0, & - nk1,nk2,nk3, nks_real, xk_collect, 1, nkBZ, eqBZ, sBZ) - ELSE - call lint ( nsymq, s, minus_q, at, bg, npk, 0,0,0, & - nk1,nk2,nk3, nks_real, xk_collect, 2, nkBZ, eqBZ, sBZ) - END IF - ! - allocate (gf(3*nat,3*nat,el_ph_nsigma)) - gf = (0.0d0,0.0d0) - ! - wqa = 1.0d0/nkfit - IF (nspin==1) wqa=degspin*wqa - ! - do ibnd = 1, nbnd - do jbnd = 1, nbnd - allocate (g2(nkBZ*nspin_lsda,3*nat,3*nat)) - allocate (g1(nksqtot,3*nat,3*nat)) - do ik = 1, nksqtot - do ii = 1, 3*nat - do jj = 1, 3*nat - g1(ik,ii,jj)=CONJG(el_ph_mat_collect(jbnd,ibnd,ik,ii))* & - el_ph_mat_collect(jbnd,ibnd,ik,jj) - enddo ! ipert - enddo !jpert - enddo ! ik - ! - allocate (g0(3*nat,3*nat)) - do i=1,nk1 - do j=1,nk2 - do k=1,nk3 - do ispin=1,nspin_lsda - nn = k-1 + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 - itemp1 = eqBZ(nn) - if (ispin==2) itemp1=itemp1+nksqtot/2 - g0(:,:) = g1(itemp1,:,:) - itemp2 = sBZ(nn) - call symm ( g0, u, xq, s, itemp2, rtau, irt, & - at, bg, nat) - if (ispin==2) nn=nn+nkBZ - g2(nn,:,:) = g0(:,:) - enddo - enddo ! k - enddo !j - enddo !i - deallocate (g0) - deallocate (g1) - ! - allocate ( point(nkBZ), noint(nkfit), ctemp(nkfit) ) - do jpert = 1, 3 * nat - do ipert = 1, 3 * nat - do ispin=1,nspin_lsda - ! - point(1:nkBZ) = & - g2(1+nkBZ*(ispin-1):nkBZ+nkBZ*(ispin-1),ipert,jpert) - ! - CALL clinear(nk1,nk2,nk3,nti,ntj,ntk,point,noint) - ! - do isig = 1, el_ph_nsigma - degauss1 = deg(isig) - do ik=1,nkfit - etk = etfit(ibnd,eqkfit(ik)+nksfit*(ispin-1)/2) - etq = etfit(jbnd,eqqfit(ik)+nksfit*(ispin-1)/2) - w0g1 = w0gauss( (effit(isig)-etk) & - / degauss1,ngauss1) / degauss1 - w0g2 = w0gauss( (effit(isig)-etq) & - / degauss1,ngauss1) / degauss1 - ctemp(ik) = noint(ik)* wqa * w0g1 * w0g2 - enddo - gf(ipert,jpert,isig) = gf(ipert,jpert,isig) + & - SUM (ctemp) - enddo ! isig - enddo ! ispin - enddo ! ipert - enddo !jpert - deallocate (point, noint, ctemp) - deallocate (g2) - ! - enddo ! ibnd - enddo ! jbnd - - deallocate (eqqfit, eqkfit) - deallocate (etfit) - deallocate (eqBZ, sBZ) -! - allocate (gam(3*nat,el_ph_nsigma), lamb(3*nat,el_ph_nsigma)) - lamb(:,:) = 0.0d0 - gam (:,:) = 0.0d0 - do isig= 1, el_ph_nsigma - do nu = 1,3*nat - gam(nu,isig) = 0.0d0 - do mu = 1, 3 * nat - do vu = 1, 3 * nat - gam(nu,isig) = gam(nu,isig) + DBLE(conjg(dyn(mu,nu)) * & - gf(mu,vu,isig) * dyn(vu,nu)) - enddo - enddo - gam(nu,isig) = gam(nu,isig) * pi/2.0d0 - ! - ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that appears - ! in the definition of the electron-phonon matrix element g - ! The sqrt(1/M) factor is actually hidden into the normal modes - ! - ! gamma = \pi \sum_k\sum_{i,j} \delta(e_{k,i}-Ef) \delta(e_{k+q,j}-Ef) - ! | \sum_mu z(mu,nu) |^2 - ! where z(mu,nu) is the mu component of normal mode nu (z = dyn) - ! gamma(nu) is the phonon linewidth of mode nu - ! - ! The factor N(Ef)^2 that appears in most formulations of el-ph interact - ! is absent because we sum, not average, over the Fermi surface. - ! The factor 2 is provided by the sum over spins - ! - if (sqrt(abs(w2(nu))) > epsw) then - ! lambda is the adimensional el-ph coupling for mode nu: - ! lambda(nu)= gamma(nu)/(pi N(Ef) \omega_{q,nu}^2) - lamb(nu,isig) = gam(nu,isig)/pi/w2(nu)/dosfit(isig) - else - lamb(nu,isig) = 0.0d0 - endif - gam(nu,isig) = gam(nu,isig)*ry_to_ghz - enddo !nu - enddo ! isig - ! - do isig= 1,el_ph_nsigma - WRITE (6, 9000) deg(isig), ngauss1 - WRITE (6, 9005) dosfit(isig), effit(isig) * rytoev - do nu=1,3*nat - WRITE (6, 9010) nu, lamb(nu,isig), gam(nu,isig) - enddo - enddo - ! Isaev: save files in suitable format for processing by lambda.x - name=TRIM(elph_dir)// 'elph.inp_lambda.' //TRIM(int_to_char(current_iq)) - - IF (ionode) THEN - open(unit=12, file=TRIM(name), form='formatted', status='unknown', & - iostat=ios) - - write(12, "(5x,3f14.6,2i6)") xq(1),xq(2),xq(3), el_ph_nsigma, 3*nat - write(12, "(6e14.6)") (w2(i), i=1,3*nat) - - do isig= 1, el_ph_nsigma - WRITE (12, 9000) deg(isig), ngauss1 - WRITE (12, 9005) dosfit(isig), effit(isig) * rytoev - do nu=1,3*nat - WRITE (12, 9010) nu, lamb(nu,isig), gam(nu,isig) - enddo - enddo - close (unit=12,status='keep') - ENDIF - ! Isaev end - CALL mp_bcast(ios, ionode_id, intra_image_comm) - IF (ios /= 0) CALL errore('elphsum','problem opening file'//TRIM(name),1) - deallocate (gam) - deallocate (lamb) - write(stdout,*) - ! - ! Prepare interface to q2r and matdyn - ! - call star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, .TRUE. ) - ! - do isig=1,el_ph_nsigma - name=TRIM(elph_dir)//'a2Fq2r.'// TRIM(int_to_char(50 + isig)) & - //'.'//TRIM(int_to_char(current_iq)) - if (ionode) then - iuelph = 4 - open(iuelph, file=name, STATUS = 'unknown', FORM = 'formatted', & - iostat=ios) - else - ! - ! this node doesn't write: unit 6 is redirected to /dev/null - ! - iuelph =6 - end if - call mp_bcast(ios, ionode_id, intra_image_comm) - IF (ios /= 0) call errore('elphsum','opening output file '// TRIM(name),1) - dyn22(:,:) = gf(:,:,isig) - write(iuelph,*) deg(isig), effit(isig), dosfit(isig) - IF ( imq == 0 ) THEN - write(iuelph,*) 2*nq - ELSE - write(iuelph,*) nq - ENDIF - xmldyn_save=xmldyn - xmldyn=.FALSE. - call q2qstar_ph (dyn22, at, bg, nat, nsym, s, invs, & - irt, rtau, nq, sxq, isq, imq, iuelph) - xmldyn=xmldyn_save - if (ionode) CLOSE( UNIT = iuelph, STATUS = 'KEEP' ) - enddo - deallocate (gf) - DEALLOCATE( deg ) - DEALLOCATE( effit ) - DEALLOCATE( dosfit ) - DEALLOCATE( xk_collect ) - IF (npool /= 1) DEALLOCATE(el_ph_mat_collect) - - ! -9000 FORMAT(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) -9005 FORMAT(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=', & - & f10.6,' eV') -9006 FORMAT(5x,'double delta at Ef =',f10.6) -9010 FORMAT(5x,'lambda(',i5,')=',f8.4,' gamma=',f8.2,' GHz') - ! - RETURN -END SUBROUTINE elphsum - -!----------------------------------------------------------------------- -subroutine elphsum2 - !----------------------------------------------------------------------- - ! - ! Sum over BZ of the electron-phonon matrix elements el_ph_mat - ! Original routine written by Francesco Mauri - ! Routine to symmetrize and print gkk (C. Verdi and S. Ponce) - ! -!#include "f_defs.h" - USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau - USE cell_base, ONLY : at, bg - use pwcom - USE symm_base, ONLY : s, irt, nsym, invs - USE kinds, only : DP - use klist, only : xk - use phcom - use el_phon - USE lr_symm_base, ONLY : rtau, nsymq, irotmq, minus_q - USE io_global, ONLY : stdout - USE control_lr, ONLY: lgamma - USE qpoint, ONLY : xq - ! - implicit none - integer :: n, ik, ikk, ikq, pbnd, ibnd, jbnd, ipert, jpert, nu, mu, vu - complex(kind=DP) :: el_ph_sum (3*nat,3*nat) - real(kind=DP) :: g2, gamma, w, epc(nbnd,nbnd,3*nat), w_1, w_2, epc_sym(nbnd,nbnd,3*nat) - real(kind=DP), parameter :: ryd2ev = 13.6058, ryd2mev = 13605.8, eps = 0.01/ryd2mev - ! - ! now read the eigenvalues and eigenvectors of the dynamical matrix - ! (already done in elphon subr) - !call readmat (iudyn, ibrav, celldm, nat, ntyp, ityp, omega, amass, tau, xq, - !w2, dyn) - ! - write(6,*) 'We here write only the matrix elements for k in the first proc:' - write(6,*) - ! - ! consider only initial state k=0 - ik = 2 - ! - if (lgamma) then - ikk = ik - ikq = ik - else - ikk = 2 * ik - 1 - ikq = ikk + 1 - endif - ! - WRITE(6,'(a,3f10.6)') 'xk ', (xk(n,ikk),n=1,3) - do ibnd = 1, nbnd - do jbnd = 1, nbnd - ! - do jpert = 1, 3 * nat - do ipert = 1, 3 * nat - el_ph_sum (ipert, jpert) = conjg (el_ph_mat (jbnd, ibnd, ik,ipert) ) * & - el_ph_mat (jbnd, ibnd, ik, jpert) - enddo - enddo - ! - ! from pert to cart - ! - call symdyn_munu2 (el_ph_sum, u, xq, s, invs, rtau, irt, at, & - bg, nsymq, nat, irotmq, minus_q) - ! - do nu = 1, nmodes - gamma = 0.0 - do mu = 1, 3 * nat - do vu = 1, 3 * nat - gamma = gamma + real (conjg (dyn (mu, nu) ) * el_ph_sum (mu,vu) & - * dyn (vu, nu) ) - enddo - enddo - gamma = gamma / 2.d0 - ! - ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that - ! appears - ! in the definition of the electron-phonon matrix element g - ! The sqrt(1/M) factor is actually hidden into the normal modes - ! we still need to divide by the phonon frequency in Ry - ! - if (w2(nu).gt.0.d0) then - w = sqrt(w2(nu)) - gamma = gamma / w - else - w = sqrt(-w2(nu)) - gamma = 0.d0 - endif - ! - if (gamma.lt.0.d0) gamma = 0.d0 - gamma = sqrt(gamma) - ! - ! gamma = |g| [Ry] - ! - epc(ibnd,jbnd,nu) = gamma - ! - enddo - ! - enddo - enddo - ! - ! HERE WE "SYMMETRIZE": actually we simply take the averages over - ! degenerate states, it is only a convention because g is gauge-dependent! - ! - ! first the phonons - ! - do ibnd = 1, nbnd - do jbnd = 1, nbnd - ! - do nu = 1, nmodes - ! - w_1 = sqrt(abs(w2(nu))) - g2 = 0.d0 - n = 0 - ! - do mu = 1, nmodes - ! - w_2 = sqrt(abs(w2(mu))) - ! - if ( abs(w_2-w_1).lt.eps ) then - n = n + 1 - g2 = g2 + epc(ibnd,jbnd,mu)*epc(ibnd,jbnd,mu) - endif - ! - enddo - ! - g2 = g2 / float(n) - epc_sym (ibnd, jbnd, nu) = sqrt (g2) - ! - enddo - ! - enddo - enddo - epc = epc_sym - ! - ! then the k electrons - ! - do nu = 1, nmodes - do jbnd = 1, nbnd - ! - do ibnd = 1, nbnd - ! - w_1 = et (ibnd, ikk) - g2 = 0.d0 - n = 0 - ! - do pbnd = 1, nbnd - ! - w_2 = et(pbnd, ikk) - ! - if ( abs(w_2-w_1).lt.eps ) then - n = n + 1 - g2 = g2 + epc(pbnd,jbnd,nu)*epc(pbnd,jbnd,nu) - endif - ! - enddo - ! - g2 = g2 / float(n) - epc_sym (ibnd, jbnd, nu) = sqrt (g2) - ! - enddo - ! - enddo - enddo - epc = epc_sym - ! - ! and finally the k+q electrons - ! - do nu = 1, nmodes - do ibnd = 1, nbnd - ! - do jbnd = 1, nbnd - ! - w_1 = et (jbnd, ikq) - g2 = 0.d0 - n = 0 - ! - do pbnd = 1, nbnd - ! - w_2 = et(pbnd, ikq) - ! - if ( abs(w_2-w_1).lt.eps ) then - n = n + 1 - g2 = g2 + epc(ibnd,pbnd,nu)*epc(ibnd,pbnd,nu) - endif - ! - enddo - ! - g2 = g2 / float(n) - epc_sym (ibnd, jbnd, nu) = sqrt (g2) - ! - enddo - ! - enddo - enddo - epc = epc_sym - ! - !if (my_pool_id.eq.0) then - ! - do ibnd = 1, nbnd - do jbnd = 1, nbnd - do nu = 1, nmodes - ! - if (w2(nu).gt.0.d0) then - w = sqrt( w2(nu)) - else - w = sqrt(-w2(nu)) - endif - ! - write(stdout,'(3i5,4f15.6)') ibnd, jbnd, nu, & - ryd2ev * et (ibnd, ikk), ryd2ev * et (jbnd, ikq), & - ryd2mev * w, ryd2mev * epc(ibnd,jbnd,nu) - ! - enddo - enddo - enddo - ! - !endif - ! - ! - return -end subroutine elphsum2 -! -!----------------------------------------------------------------------- -SUBROUTINE elphsum_simple - !----------------------------------------------------------------------- - ! - ! Sum over BZ of the electron-phonon matrix elements el_ph_mat - ! Original routine written by Francesco Mauri - ! Rewritten by Matteo Calandra - !----------------------------------------------------------------------- - USE kinds, ONLY : DP - USE constants, ONLY : pi, ry_to_cmm1, ry_to_ghz, rytoev - USE ions_base, ONLY : nat - USE cell_base, ONLY : at, bg - USE symm_base, ONLY : s, irt, nsym, invs - USE klist, ONLY : xk, nelec, nks, wk - USE wvfct, ONLY : nbnd, et - USE el_phon, ONLY : el_ph_mat, el_ph_nsigma, el_ph_ngauss, el_ph_sigma - USE mp_pools, ONLY : inter_pool_comm, npool - USE mp_images, ONLY : intra_image_comm - USE output, ONLY : fildyn - USE dynmat, ONLY : dyn, w2 - USE modes, ONLY : u, nirr - USE control_ph, only : current_iq, qplot - USE lsda_mod, only : isk - USE el_phon, ONLY : done_elph, gamma_disp - USE io_global, ONLY : stdout, ionode, ionode_id - USE mp, ONLY: mp_sum, mp_bcast - - USE lr_symm_base, ONLY : rtau, nsymq, irotmq, minus_q - USE qpoint, ONLY : xq, nksq, ikks, ikqs - ! - IMPLICIT NONE - REAL(DP), PARAMETER :: eps = 20_dp/ry_to_cmm1 ! eps = 20 cm^-1, in Ry - ! - INTEGER :: ik, ikk, ikq, isig, ibnd, jbnd, ipert, jpert, nu, mu, & - vu, ngauss1, iuelph, ios, irr - INTEGER :: nmodes - REAL(DP) :: weight, w0g1, w0g2, w0gauss, wgauss, degauss1, dosef, & - ef1, phase_space, lambda, gamma - REAL(DP), EXTERNAL :: dos_ef, efermig - character(len=80) :: filelph - CHARACTER(len=256) :: file_elphmat - ! - COMPLEX(DP) :: el_ph_sum (3*nat,3*nat), dyn_corr(3*nat,3*nat) - - INTEGER, EXTERNAL :: find_free_unit - CHARACTER(LEN=6) :: int_to_char - - - DO irr=1,nirr - IF (.NOT.done_elph(irr)) RETURN - ENDDO - - nmodes=3*nat - - filelph=TRIM(fildyn)//'.elph.'//TRIM(int_to_char(current_iq)) - - ! parallel case: only first node writes - IF ( ionode ) THEN - ! - iuelph = find_free_unit() - OPEN (unit = iuelph, file = TRIM(filelph), status = 'unknown', err = & - 100, iostat = ios) - REWIND (iuelph) - ELSE - iuelph = 0 - ! - END IF -100 CONTINUE - CALL mp_bcast(ios,ionode_id,intra_image_comm) - CALL errore ('elphsum_simple', 'opening file '//filelph, ABS (ios) ) - - IF (ionode) THEN - WRITE (iuelph, '(3f15.8,2i8)') xq, el_ph_nsigma, 3 * nat - WRITE (iuelph, '(6e14.6)') (w2 (nu) , nu = 1, nmodes) - ENDIF - - - ngauss1=0 - DO isig = 1, el_ph_nsigma - ! degauss1 = 0.01 * isig - degauss1 = el_ph_sigma * isig - el_ph_sum(:,:) = (0.d0, 0.d0) - phase_space = 0.d0 - ! - ! Recalculate the Fermi energy Ef=ef1 and the DOS at Ef, dosef = N(Ef) - ! for this gaussian broadening - ! - ! Note that the weights of k+q points must be set to zero for the - ! following call to yield correct results - ! - - ef1 = efermig (et, nbnd, nks, nelec, wk, degauss1, el_ph_ngauss, 0, isk) - dosef = dos_ef (el_ph_ngauss, degauss1, ef1, et, wk, nks, nbnd) - ! N(Ef) is the DOS per spin, not summed over spin - dosef = dosef / 2.d0 - ! - ! Sum over bands with gaussian weights - ! - - DO ik = 1, nksq - - ! - ! see subroutine elphel for the logic of indices - ! - ikk = ikks(ik) - ikq = ikqs(ik) - DO ibnd = 1, nbnd - w0g1 = w0gauss ( (ef1 - et (ibnd, ikk) ) / degauss1, ngauss1) & - / degauss1 - DO jbnd = 1, nbnd - w0g2 = w0gauss ( (ef1 - et (jbnd, ikq) ) / degauss1, ngauss1) & - / degauss1 - ! note that wk(ikq)=wk(ikk) - weight = wk (ikk) * w0g1 * w0g2 - DO jpert = 1, 3 * nat - DO ipert = 1, 3 * nat - el_ph_sum (ipert, jpert) = el_ph_sum (ipert, jpert) + weight * & - CONJG (el_ph_mat (jbnd, ibnd, ik, ipert) ) * & - el_ph_mat (jbnd, ibnd, ik, jpert) - ENDDO - ENDDO - phase_space = phase_space+weight - ENDDO - ENDDO - - ENDDO - - ! el_ph_sum(mu,nu)=\sum_k\sum_{i,j}[ - ! x - ! x \delta(e_{k,i}-Ef) \delta(e_{k+q,j} - ! - ! collect contributions from all pools (sum over k-points) - ! - - call mp_sum ( el_ph_sum , inter_pool_comm ) - call mp_sum ( phase_space , inter_pool_comm ) - - ! - ! symmetrize el_ph_sum(mu,nu) : it transforms as the dynamical matrix - ! - - CALL symdyn_munu_new (el_ph_sum, u, xq, s, invs, rtau, irt, at, & - bg, nsymq, nat, irotmq, minus_q) - ! - WRITE (stdout, *) - WRITE (stdout, 9000) degauss1, ngauss1 - WRITE (stdout, 9005) dosef, ef1 * rytoev - WRITE (stdout, 9006) phase_space - IF (ionode) THEN - WRITE (iuelph, 9000) degauss1, ngauss1 - WRITE (iuelph, 9005) dosef, ef1 * rytoev - ENDIF - - DO nu = 1, nmodes - gamma = 0.d0 - DO mu = 1, 3 * nat - DO vu = 1, 3 * nat - gamma = gamma + DBLE (CONJG (dyn (mu, nu) ) * el_ph_sum (mu, vu)& - * dyn (vu, nu) ) - ENDDO - ENDDO - gamma = pi * gamma / 2.d0 - ! - ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that appears - ! in the definition of the electron-phonon matrix element g - ! The sqrt(1/M) factor is actually hidden into the normal modes - ! - ! gamma = \pi \sum_k\sum_{i,j} \delta(e_{k,i}-Ef) \delta(e_{k+q,j}-Ef) - ! | \sum_mu z(mu,nu) |^2 - ! where z(mu,nu) is the mu component of normal mode nu (z = dyn) - ! gamma(nu) is the phonon linewidth of mode nu - ! - ! The factor N(Ef)^2 that appears in most formulations of el-ph interact - ! is absent because we sum, not average, over the Fermi surface. - ! The factor 2 is provided by the sum over spins - ! - IF (SQRT (ABS (w2 (nu) ) ) > eps) THEN - ! lambda is the adimensional el-ph coupling for mode nu: - ! lambda(nu)= gamma(nu)/(pi N(Ef) \omega_{q,nu}^2) - lambda = gamma / pi / w2 (nu) / dosef - ELSE - lambda = 0.d0 - ENDIF - WRITE (stdout, 9010) nu, lambda, gamma * ry_to_gHz - IF (ionode) WRITE (iuelph, 9010) nu, lambda, gamma * ry_to_gHz - IF (qplot) gamma_disp(nu,isig,current_iq) = gamma * ry_to_gHz - ENDDO - ENDDO - - -9000 FORMAT(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) -9005 FORMAT(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=', & - & f10.6,' eV') -9006 FORMAT(5x,'double delta at Ef =',f10.6) -9010 FORMAT(5x,'lambda(',i5,')=',f8.4,' gamma=',f8.2,' GHz') - ! - ! - IF (ionode) CLOSE (unit = iuelph) - RETURN - - - - ! call star_q(x_q(1,iq), at, bg, nsym , s , invs , nq, sxq, & - ! isq, imq, .FALSE. ) - - -END SUBROUTINE elphsum_simple - -!----------------------------------------------------------------------- -SUBROUTINE elphfil_epa(iq) - !----------------------------------------------------------------------- - ! - ! EPA (electron-phonon-averaged) approximation - ! Writes electron-phonon matrix elements to a file - ! Written by Georgy Samsonidze on 2015-01-28 - ! - ! Adv. Energy Mater. 2018, 1800246 - ! doi:10.1002/aenm.201800246 - ! https://doi.org/10.1002/aenm.201800246 - ! - !----------------------------------------------------------------------- - USE cell_base, ONLY : ibrav, alat, omega, tpiba, at, bg - USE disp, ONLY : nq1, nq2, nq3, nqs, x_q, wq, lgamma_iq - USE dynmat, ONLY : dyn, w2 - USE el_phon, ONLY : el_ph_mat, done_elph - USE fft_base, ONLY : dfftp, dffts, dfftb - USE gvect, ONLY : ngm_g, ecutrho - USE io_global, ONLY : ionode, ionode_id - USE ions_base, ONLY : nat, nsp, atm, ityp, tau - USE kinds, ONLY : DP - USE klist, ONLY : xk, wk, nelec, nks, nkstot, ngk - USE lsda_mod, ONLY : nspin, isk - USE modes, ONLY : nirr, nmodes, npert, npertx, u, t, tmq, & - name_rap_mode, num_rap_mode - USE lr_symm_base, ONLY : irgq, nsymq, irotmq, rtau, gi, gimq, & - minus_q, invsymq - USE mp, ONLY : mp_bcast, mp_sum - USE mp_images, ONLY : intra_image_comm - USE mp_pools, ONLY : npool, intra_pool_comm - USE qpoint, ONLY : nksq, nksqtot, ikks, ikqs, eigqts - USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3 - USE symm_base, ONLY : s, invs, ft, nrot, nsym, nsym_ns, & - nsym_na, ft, sr, sname, t_rev, irt, time_reversal, & - invsym, nofrac, allfrac, nosym, nosym_evc, no_t_rev - USE wvfct, ONLY : nbnd, et, wg - USE gvecw, ONLY : ecutwfc - USE io_files, ONLY : prefix - - IMPLICIT NONE - - INTEGER, INTENT(IN) :: iq - - INTEGER :: iuelph, ios, irr, ii, jj, kk, ll - character :: cdate*9, ctime*9, sdate*32, stime*32, & - stitle*32, myaccess*10, mystatus*7 - CHARACTER(LEN=80) :: filelph - - REAL(DP), ALLOCATABLE :: xk_collect(:,:), wk_collect(:) - REAL(DP), ALLOCATABLE :: et_collect(:,:), wg_collect(:,:) - INTEGER, ALLOCATABLE :: ngk_collect(:) - INTEGER, ALLOCATABLE :: ikks_collect(:), ikqs_collect(:) - COMPLEX(DP), ALLOCATABLE :: el_ph_mat_collect(:,:,:,:) - INTEGER :: ftau(3,48) - INTEGER, EXTERNAL :: find_free_unit, atomic_number - - filelph = TRIM(prefix) // '.epa.k' - - DO irr = 1, nirr - IF (.NOT. done_elph(irr)) RETURN - ENDDO - - IF (iq .EQ. 1) THEN - myaccess = 'sequential' - mystatus = 'replace' - ELSE - myaccess = 'append' - mystatus = 'old' - ENDIF - IF (ionode) THEN - iuelph = find_free_unit() - OPEN(unit = iuelph, file = TRIM(filelph), form = 'unformatted', & - access = myaccess, status = mystatus, iostat = ios) - ELSE - iuelph = 0 - ENDIF - CALL mp_bcast(ios, ionode_id, intra_image_comm) - CALL errore('elphfil_epa', 'opening file ' // filelph, ABS(ios)) - - IF (iq .EQ. 1) THEN - CALL date_and_tim(cdate, ctime) - WRITE(sdate, '(A2,"-",A3,"-",A4,21X)') cdate(1:2), cdate(3:5), cdate(6:9) - WRITE(stime, '(A8,24X)') ctime(1:8) - WRITE(stitle, '("EPA-Complex",21X)') - CALL cryst_to_cart(nqs, x_q, at, -1) - ! write header - IF (ionode) THEN - WRITE(iuelph) stitle, sdate, stime - WRITE(iuelph) ibrav, nat, nsp, nrot, nsym, nsym_ns, nsym_na, & - ngm_g, nspin, nbnd, nmodes, nqs - WRITE(iuelph) nq1, nq2, nq3, nk1, nk2, nk3, k1, k2, k3 - WRITE(iuelph) time_reversal, invsym, nofrac, allfrac, nosym, & - nosym_evc, no_t_rev - WRITE(iuelph) alat, omega, tpiba, nelec, ecutrho, ecutwfc - WRITE(iuelph) dfftp%nr1, dfftp%nr2, dfftp%nr3 - WRITE(iuelph) dffts%nr1, dffts%nr2, dffts%nr3 - WRITE(iuelph) dfftb%nr1, dfftb%nr2, dfftb%nr3 - WRITE(iuelph) ((at(ii, jj), ii = 1, 3), jj = 1, 3) - WRITE(iuelph) ((bg(ii, jj), ii = 1, 3), jj = 1, 3) - WRITE(iuelph) (atomic_number(atm(ii)), ii = 1, nsp) - WRITE(iuelph) (ityp(ii), ii = 1, nat) - WRITE(iuelph) ((tau(ii, jj), ii = 1, 3), jj = 1, nat) - WRITE(iuelph) ((x_q(ii, jj), ii = 1, 3), jj = 1, nqs) - WRITE(iuelph) (wq(ii), ii = 1, nqs) - WRITE(iuelph) (lgamma_iq(ii), ii = 1, nqs) - ENDIF - CALL cryst_to_cart(nqs, x_q, bg, 1) - ENDIF - - ! collect data for current q-point - ALLOCATE(xk_collect(3, nkstot)) - ALLOCATE(wk_collect(nkstot)) - ALLOCATE(et_collect(nbnd, nkstot)) - ALLOCATE(wg_collect(nbnd, nkstot)) - ALLOCATE(ngk_collect(nkstot)) - ALLOCATE(ikks_collect(nksqtot)) - ALLOCATE(ikqs_collect(nksqtot)) - ALLOCATE(el_ph_mat_collect(nbnd, nbnd, nksqtot, nmodes)) - IF (npool > 1) THEN - CALL poolcollect(3, nks, xk, nkstot, xk_collect) - CALL poolcollect(1, nks, wk, nkstot, wk_collect) - CALL poolcollect(nbnd, nks, et, nkstot, et_collect) - CALL poolcollect(nbnd, nks, wg, nkstot, wg_collect) - CALL ipoolcollect(1, nks, ngk, nkstot, ngk_collect) - CALL jpoolcollect(1, nksq, ikks, nksqtot, ikks_collect) - CALL jpoolcollect(1, nksq, ikqs, nksqtot, ikqs_collect) - CALL el_ph_collect(nmodes, el_ph_mat, el_ph_mat_collect, nksqtot, nksq) - ELSE - xk_collect(1:3, 1:nks) = xk(1:3, 1:nks) - wk_collect(1:nks) = wk(1:nks) - et_collect(1:nbnd, 1:nks) = et(1:nbnd, 1:nks) - wg_collect(1:nbnd, 1:nks) = wg(1:nbnd, 1:nks) - ngk_collect(1:nks) = ngk(1:nks) - ikks_collect(1:nksq) = ikks(1:nksq) - ikqs_collect(1:nksq) = ikqs(1:nksq) - el_ph_mat_collect(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes) = & - el_ph_mat(1:nbnd, 1:nbnd, 1:nksq, 1:nmodes) - ENDIF - CALL cryst_to_cart(nkstot, xk_collect, at, -1) - ! write data for current q-point - IF (ionode) THEN - WRITE(iuelph) nsymq, irotmq, nirr, npertx, nkstot, nksqtot - WRITE(iuelph) minus_q, invsymq - WRITE(iuelph) (irgq(ii), ii = 1, 48) - WRITE(iuelph) (npert(ii), ii = 1, nmodes) - WRITE(iuelph) (((rtau(ii, jj, kk), ii = 1, 3), jj = 1, 48), & - kk = 1, nat) - WRITE(iuelph) ((gi(ii, jj), ii = 1, 3), jj = 1, 48) - WRITE(iuelph) (gimq(ii), ii = 1, 3) - WRITE(iuelph) ((u(ii, jj), ii = 1, nmodes), jj = 1, nmodes) - WRITE(iuelph) ((((t(ii, jj, kk, ll), ii = 1, npertx), & - jj = 1, npertx), kk = 1, 48), ll = 1, nmodes) - WRITE(iuelph) (((tmq(ii, jj, kk), ii = 1, npertx), & - jj = 1, npertx), kk = 1, nmodes) - WRITE(iuelph) (name_rap_mode(ii), ii = 1, nmodes) - WRITE(iuelph) (num_rap_mode(ii), ii = 1, nmodes) - WRITE(iuelph) (((s(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48) - WRITE(iuelph) (invs(ii), ii = 1, 48) - ! FIXME: should disappear - ftau(1,1:48) = NINT(ft(1,1:48)*dfftp%nr1) - ftau(2,1:48) = NINT(ft(2,1:48)*dfftp%nr2) - ftau(3,1:48) = NINT(ft(3,1:48)*dfftp%nr3) - WRITE(iuelph) ((ftau(ii, jj), ii = 1, 3), jj = 1, 48) - ! end FIXME - WRITE(iuelph) ((ft(ii, jj), ii = 1, 3), jj = 1, 48) - WRITE(iuelph) (((sr(ii, jj, kk), ii = 1, 3), jj = 1, 3), kk = 1, 48) - WRITE(iuelph) (sname(ii), ii = 1, 48) - WRITE(iuelph) (t_rev(ii), ii = 1, 48) - WRITE(iuelph) ((irt(ii, jj), ii = 1, 48), jj = 1, nat) - WRITE(iuelph) ((xk_collect(ii, jj), ii = 1, 3), jj = 1, nkstot) - WRITE(iuelph) (wk_collect(ii), ii = 1, nkstot) - WRITE(iuelph) ((et_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot) - WRITE(iuelph) ((wg_collect(ii, jj), ii = 1, nbnd), jj = 1, nkstot) - WRITE(iuelph) (isk(ii), ii = 1, nkstot) - WRITE(iuelph) (ngk_collect(ii), ii = 1, nkstot) - WRITE(iuelph) (ikks_collect(ii), ii = 1, nksqtot) - WRITE(iuelph) (ikqs_collect(ii), ii = 1, nksqtot) - WRITE(iuelph) (eigqts(ii), ii = 1, nat) - WRITE(iuelph) (w2(ii), ii = 1, nmodes) - WRITE(iuelph) ((dyn(ii, jj), ii = 1, nmodes), jj = 1, nmodes) - WRITE(iuelph) ((((el_ph_mat_collect(ii, jj, kk, ll), ii = 1, nbnd), & - jj = 1, nbnd), kk = 1, nksqtot), ll = 1, nmodes) - CLOSE (unit = iuelph, status = 'keep') - ENDIF - CALL cryst_to_cart(nkstot, xk_collect, bg, 1) - DEALLOCATE(xk_collect) - DEALLOCATE(wk_collect) - DEALLOCATE(et_collect) - DEALLOCATE(wg_collect) - DEALLOCATE(ngk_collect) - DEALLOCATE(ikks_collect) - DEALLOCATE(ikqs_collect) - DEALLOCATE(el_ph_mat_collect) - - RETURN - -END SUBROUTINE elphfil_epa - -!---------------------------------------------------------------------------- -SUBROUTINE ipoolcollect( length, nks, f_in, nkstot, f_out ) - !---------------------------------------------------------------------------- - ! - ! ... as poolcollect, for an integer vector - ! - USE mp_pools, ONLY : my_pool_id, npool, kunit, & - inter_pool_comm, intra_pool_comm - USE mp, ONLY : mp_sum - ! - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: length, nks, nkstot - ! first dimension of arrays - ! number of k-points per pool - ! total number of k-points - INTEGER, INTENT(IN) :: f_in (length,nks) - ! pool-distributed function - INTEGER, INTENT(OUT) :: f_out(length,nkstot) - ! pool-collected function - ! - INTEGER :: nbase, rest, nks1 - ! - nks1 = kunit * ( nkstot / kunit / npool ) - ! - rest = ( nkstot - nks1 * npool ) / kunit - ! - IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + kunit - ! - IF (nks1.ne.nks) & - call errore('ipoolcollect','inconsistent number of k-points',1) - ! - ! ... calculates nbase = the position in the list of the first point that - ! ... belong to this npool - 1 - ! - nbase = nks * my_pool_id - ! - IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest * kunit - ! - ! copy the original points in the correct position of the list - ! - f_out=0 - f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks) - ! - CALL mp_sum( f_out, inter_pool_comm ) - ! - RETURN - ! -END SUBROUTINE ipoolcollect - -!---------------------------------------------------------------------------- -SUBROUTINE jpoolcollect( length, nks, f_in, nkstot, f_out ) - !---------------------------------------------------------------------------- - ! - ! ... as ipoolcollect, without kunit and with an index shift - ! - USE mp_pools, ONLY : my_pool_id, npool, kunit, & - inter_pool_comm, intra_pool_comm - USE mp, ONLY : mp_sum - ! - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: length, nks, nkstot - ! first dimension of arrays - ! number of k-points per pool - ! total number of k-points - INTEGER, INTENT(IN) :: f_in (length,nks) - ! pool-distributed function - INTEGER, INTENT(OUT) :: f_out(length,nkstot) - ! pool-collected function - ! - INTEGER :: nbase, rest, nks1 - ! - nks1 = ( nkstot / npool ) - ! - rest = ( nkstot - nks1 * npool ) - ! - IF ( ( my_pool_id + 1 ) <= rest ) nks1 = nks1 + 1 - ! - IF (nks1.ne.nks) & - call errore('jpoolcollect','inconsistent number of k-points',1) - ! - ! ... calculates nbase = the position in the list of the first point that - ! ... belong to this npool - 1 - ! - nbase = nks * my_pool_id - ! - IF ( ( my_pool_id + 1 ) > rest ) nbase = nbase + rest - ! - ! copy the original points in the correct position of the list - ! - f_out=0 - f_out(:,nbase+1:nbase+nks) = f_in(:,1:nks) + nbase * kunit - ! - CALL mp_sum( f_out, inter_pool_comm ) - ! - RETURN - ! -END SUBROUTINE jpoolcollect - -!----------------------------------------------------------------------- -FUNCTION dos_ef (ngauss, degauss, ef, et, wk, nks, nbnd) - !----------------------------------------------------------------------- - ! - USE kinds, ONLY : DP - USE mp_pools, ONLY : inter_pool_comm - USE mp, ONLY : mp_sum - IMPLICIT NONE - REAL(DP) :: dos_ef - INTEGER :: ngauss, nbnd, nks - REAL(DP) :: et (nbnd, nks), wk (nks), ef, degauss - ! - INTEGER :: ik, ibnd - REAL(DP), EXTERNAL :: w0gauss - ! - ! Compute DOS at E_F (states per Ry per unit cell) - ! - dos_ef = 0.0d0 - DO ik = 1, nks - DO ibnd = 1, nbnd - dos_ef = dos_ef + wk (ik) * w0gauss ( (et (ibnd, ik) - ef) & - / degauss, ngauss) / degauss - ENDDO - ENDDO - ! - ! Collects partial sums on k-points from all pools - ! - CALL mp_sum ( dos_ef, inter_pool_comm ) - ! - RETURN -END FUNCTION dos_ef - -!a2F -subroutine lint ( nsym, s, minus_q, at, bg, npk, k1,k2,k3, & - nk1,nk2,nk3, nks, xk, kunit, nkBZ, eqBZ, sBZ) - !----------------------------------------------------------------------- - ! - ! Find which k-points of a uniform grid are in the IBZ - ! - use kinds, only : DP - implicit none - integer, intent (IN) :: nks, nsym, s(3,3,48), npk, k1, k2, k3, & - nk1, nk2, nk3, kunit, nkBZ - logical, intent (IN) :: minus_q - real(kind=DP), intent(IN):: at(3,3), bg(3,3), xk(3,npk) - integer, INTENT(OUT) :: eqBZ(nkBZ), sBZ(nkBZ) - ! - real(kind=DP) :: xkr(3), deltap(3), deltam(3) - real(kind=DP), parameter:: eps=1.0d-5 - real(kind=DP), allocatable :: xkg(:,:), xp(:,:) - integer :: i,j,k, ns, n, nk - integer :: nkh - ! - ! Re-generate a uniform grid of k-points xkg - ! - allocate (xkg( 3,nkBZ)) - ! - if(kunit < 1 .or. kunit > 2) call errore('lint','bad kunit value',kunit) - ! - ! kunit=2: get only "true" k points, not k+q points, from the list - ! - nkh = nks/kunit - allocate (xp(3,nkh)) - if (kunit == 1) then - xp(:,1:nkh) = xk(:,1:nkh) - else - do j=1,nkh - xp(:,j) = xk(:,2*j-1) - enddo - end if - do i=1,nk1 - do j=1,nk2 - do k=1,nk3 - n = (k-1) + (j-1)*nk3 + (i-1)*nk2*nk3 + 1 - 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 - end do - end do - end do - - call cryst_to_cart (nkh,xp,at,-1) - - do nk=1,nkBZ - do n=1,nkh - do ns=1,nsym - do i=1,3 - xkr(i) = s(i,1,ns) * xp(1,n) + & - s(i,2,ns) * xp(2,n) + & - s(i,3,ns) * xp(3,n) - end do - do i=1,3 - deltap(i) = xkr(i)-xkg(i,nk) - nint (xkr(i)-xkg(i,nk) ) - deltam(i) = xkr(i)+xkg(i,nk) - nint (xkr(i)+xkg(i,nk) ) - end do - if ( sqrt ( deltap(1)**2 + & - deltap(2)**2 + & - deltap(3)**2 ) < eps .or. ( minus_q .and. & - sqrt ( deltam(1)**2 + & - deltam(2)**2 + & - deltam(3)**2 ) < eps ) ) then - eqBZ(nk) = n - sBZ(nk) = ns - go to 15 - end if - end do - end do - call errore('lint','cannot locate k point xk',nk) -15 continue - end do - - do n=1,nkh - do nk=1,nkBZ - if (eqBZ(nk) == n) go to 20 - end do - ! this failure of the algorithm may indicate that the displaced grid - ! (with k1,k2,k3.ne.0) does not have the full symmetry of the lattice - call errore('lint','cannot remap grid on k-point list',n) -20 continue - end do - - deallocate(xkg) - deallocate(xp) - - return -end subroutine lint diff --git a/test-suite/not_epw_comp/phq_readin.f90 b/test-suite/not_epw_comp/phq_readin.f90 deleted file mode 100644 index f98686fe4..000000000 --- a/test-suite/not_epw_comp/phq_readin.f90 +++ /dev/null @@ -1,871 +0,0 @@ -! -! Copyright (C) 2001-2018 Quantum ESPRESSO group -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! -! -!---------------------------------------------------------------------------- -SUBROUTINE phq_readin() - !---------------------------------------------------------------------------- - ! - ! This routine reads the control variables for the program phononq. - ! from standard input (unit 5). - ! A second routine readfile reads the variables saved on a file - ! by the self-consistent program. - ! - ! - USE kinds, ONLY : DP - USE ions_base, ONLY : nat, ntyp => nsp - USE mp, ONLY : mp_bcast - USE mp_world, ONLY : world_comm - USE ions_base, ONLY : amass, atm - USE check_stop, ONLY : max_seconds - USE start_k, ONLY : reset_grid - USE klist, ONLY : xk, nks, nkstot, lgauss, two_fermi_energies, ltetra - USE control_flags, ONLY : gamma_only, tqr, restart, io_level, & - ts_vdw, ldftd3, lxdm, isolve - USE funct, ONLY : dft_is_meta, dft_is_hybrid - USE uspp, ONLY : okvan - USE fixed_occ, ONLY : tfixed_occ - USE lsda_mod, ONLY : lsda, nspin - USE fft_base, ONLY : dffts - USE spin_orb, ONLY : domag - USE cellmd, ONLY : lmovecell - USE run_info, ONLY : title - USE control_ph, ONLY : maxter, alpha_mix, lgamma_gamma, epsil, & - zue, zeu, xmldyn, newgrid, & - trans, reduce_io, tr2_ph, niter_ph, & - nmix_ph, ldisp, recover, lnoloc, start_irr, & - last_irr, start_q, last_q, current_iq, tmp_dir_ph, & - ext_recover, ext_restart, u_from_file, ldiag, & - search_sym, lqdir, electron_phonon, tmp_dir_phq, & - rec_code_read, qplot, only_init, only_wfc, & - low_directory_check, nk1, nk2, nk3, k1, k2, k3 - USE save_ph, ONLY : tmp_dir_save, save_ph_input_variables - USE gamma_gamma, ONLY : asr - USE partial, ONLY : atomo, nat_todo, nat_todo_input - USE output, ONLY : fildyn, fildvscf, fildrho - USE disp, ONLY : nq1, nq2, nq3, x_q, wq, nqs, lgamma_iq - USE io_files, ONLY : tmp_dir, prefix, postfix, create_directory, & - check_tempdir, xmlpun_schema - USE noncollin_module, ONLY : i_cons, noncolin - USE control_flags, ONLY : iverbosity, modenum - USE io_global, ONLY : meta_ionode, meta_ionode_id, ionode, ionode_id, stdout - USE mp_images, ONLY : nimage, my_image_id, intra_image_comm, & - me_image, nproc_image - USE mp_pools, ONLY : npool - USE paw_variables, ONLY : okpaw - USE ramanm, ONLY : eth_rps, eth_ns, lraman, elop, dek - USE freq_ph, ONLY : fpol, fiu, nfs - USE cryst_ph, ONLY : magnetic_sym - USE ph_restart, ONLY : ph_readfile - USE el_phon, ONLY : elph,elph_mat,elph_simple,elph_epa,elph_nbnd_min, elph_nbnd_max, & - el_ph_sigma, el_ph_nsigma, el_ph_ngauss,auxdvscf - USE dfile_star, ONLY : drho_star, dvscf_star - - USE qpoint, ONLY : nksq, xq - USE control_lr, ONLY : lgamma, lrpa - ! YAMBO > - USE YAMBO, ONLY : elph_yambo,dvscf_yambo - ! YAMBO < - USE elph_tetra_mod,ONLY : elph_tetra, lshift_q, in_alpha2f - USE ktetra, ONLY : tetra_type - USE ldaU, ONLY : lda_plus_u, U_projection, lda_plus_u_kind - USE ldaU_ph, ONLY : read_dns_bare, d2ns_type - ! - IMPLICIT NONE - ! - CHARACTER(LEN=256), EXTERNAL :: trimcheck - ! - INTEGER :: ios, ipol, iter, na, it, ierr, ierr1 - ! integer variable for I/O control - ! counter on polarizations - ! counter on iterations - ! counter on atoms - ! counter on types - REAL(DP), ALLOCATABLE :: amass_input(:) - ! save masses read from input here - CHARACTER (LEN=256) :: outdir, filename - CHARACTER (LEN=8) :: verbosity - CHARACTER(LEN=80) :: card - CHARACTER(LEN=6) :: int_to_char - INTEGER :: i - LOGICAL :: nogg - LOGICAL :: q2d, q_in_band_form - INTEGER, EXTERNAL :: atomic_number - REAL(DP), EXTERNAL :: atom_weight - LOGICAL, EXTERNAL :: imatches - LOGICAL, EXTERNAL :: has_xml - LOGICAL :: exst, parallelfs - REAL(DP), ALLOCATABLE :: xqaux(:,:) - INTEGER, ALLOCATABLE :: wqaux(:) - INTEGER :: nqaux, iq - CHARACTER(len=80) :: diagonalization='david' - ! - NAMELIST / INPUTPH / tr2_ph, amass, alpha_mix, niter_ph, nmix_ph, & - nat_todo, verbosity, iverbosity, outdir, epsil, & - trans, zue, zeu, max_seconds, reduce_io, & - modenum, prefix, fildyn, fildvscf, fildrho, & - ldisp, nq1, nq2, nq3, & - eth_rps, eth_ns, lraman, elop, dek, recover, & - fpol, asr, lrpa, lnoloc, start_irr, last_irr, & - start_q, last_q, nogg, ldiag, search_sym, lqdir, & - nk1, nk2, nk3, k1, k2, k3, & - drho_star, dvscf_star, only_init, only_wfc, & - elph_nbnd_min, elph_nbnd_max, el_ph_ngauss, & - el_ph_nsigma, el_ph_sigma, electron_phonon, & - q_in_band_form, q2d, qplot, low_directory_check, & - lshift_q, read_dns_bare, d2ns_type, diagonalization - - ! tr2_ph : convergence threshold - ! amass : atomic masses - ! alpha_mix : the mixing parameter - ! niter_ph : maximum number of iterations - ! nmix_ph : number of previous iterations used in mixing - ! nat_todo : number of atom to be displaced - ! verbosity : verbosity control (iverbosity is obsolete) - ! outdir : directory where input, output, temporary files reside - ! epsil : if true calculate dielectric constant - ! trans : if true calculate phonon - ! electron-phonon : select the kind of electron-phonon calculation - ! elph : if true calculate electron-phonon coefficients - ! elph_mat : if true eph coefficients for wannier - ! zue : if .true. calculate effective charges ( d force / dE ) - ! zeu : if .true. calculate effective charges ( d P / du ) - ! lraman : if true calculate raman tensor - ! elop : if true calculate electro-optic tensor - ! max_seconds : maximum cputime for this run - ! reduce_io : reduce I/O to the strict minimum - ! modenum : single mode calculation - ! prefix : the prefix of files produced by pwscf - ! fildyn : output file for the dynamical matrix - ! fildvscf : output file containing deltavsc - ! fildrho : output file containing deltarho - ! fildrho_dir : directory where fildrho files will be stored (default: outdir or ESPRESSO_FILDRHO_DIR variable) - ! eth_rps : threshold for calculation of Pc R |psi> (Raman) - ! eth_ns : threshold for non-scf wavefunction calculation (Raman) - ! dek : delta_xk used for wavefunctions derivation (Raman) - ! recover : recover=.true. to restart from an interrupted run - ! asr : in the gamma_gamma case apply acoustic sum rule - ! start_q : in q list does the q points from start_q to last_q - ! last_q : - ! start_irr : does the irred. representation from start_irr to last_irr - ! last_irr : - ! nogg : if .true. lgamma_gamma tricks are not used - ! ldiag : if .true. force diagonalization of the dyn mat - ! lqdir : if .true. each q writes in its own directory - ! search_sym : if .true. analyze symmetry if possible - ! nk1,nk2,nk3, k1,k2,k3 : - ! when specified in input, the phonon run uses a different - ! k-point mesh from that used for the charge density. - ! - ! dvscf_star%open : if .true. write in dvscf_star%dir the dvscf_q - ! 'for all q' in the star of q with suffix dvscf_star%ext. - ! The dvscf_q' is written in the basis dvscf_star%basis; - ! if dvscf_star%pat is .true. also save a pattern file. - ! dvscf_star%dir, dvscf_star%ext, dvscf_star%basis : see dvscf_star%open - ! drho_star%open : like dvscf_star%open but for drho_q - ! drho_star%dir, drho_star%ext, drho_star%basis : see drho_star%open - ! - ! elph_nbnd_min, - ! elph_nbnd_max: if (elph_mat=.true.) it dumps the eph matrix element from elph_nbnd_min - ! to elph_nbnd_max - ! el_ph_ngauss, - ! el_ph_nsigma, - ! el_ph_sigma : if (elph_mat=.true.) it defines the kind and the val-ue of the - ! smearing to be used in the eph coupling calculation. - ! qplot, : if true a list of q points is given in input - ! q_in_band_form: if true the input list of q points defines paths - ! q2d, : if .true. the q list define a mesh in a square. - ! low_directory_check : if .true. only the requested representations - ! are searched on file - ! - ! read_dns_bare : If .true. the code tries to read three files in DFPT+U calculations: - ! dnsorth, dnsbare, d2nsbare - ! d2ns_type : DFPT+U - the 2nd bare derivative of occupation matrices ns - ! (d2ns_bare matrix). Experimental! This is why it is not documented in Doc. - ! d2ns_type='full': matrix calculated with no approximation. - ! d2ns_type='fmmp': assume a m <=> m' symmetry. - ! d2ns_type='diag': if okvan=.true. the matrix is calculated retaining only - ! for <\beta_J|\phi_I> products where for J==I. - ! d2ns_type='dmmp': same as 'diag', but also assuming a m <=> m'. - ! diagonalization : diagonalization method used in the nscf calc - ! - ! Note: meta_ionode is a single processor that reads the input - ! (ionode is also a single processor but per image) - ! Data read from input is subsequently broadcast to all processors - ! from ionode_id (using the default communicator world_comm) - ! - IF (meta_ionode) THEN - ! - ! ... Input from file ? - ! - CALL input_from_file ( ) - ! - ! ... Read the first line of the input file - ! - READ( 5, '(A)', IOSTAT = ios ) title - ! - ENDIF - ! - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore( 'phq_readin', 'reading title ', ABS( ios ) ) - CALL mp_bcast(title, meta_ionode_id, world_comm ) - ! - ! Rewind the input if the title is actually the beginning of inputph namelist - ! - IF( imatches("&inputph", title) ) THEN - WRITE(*, '(6x,a)') "Title line not specified: using 'default'." - title='default' - IF (meta_ionode) REWIND(5, iostat=ios) - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore('phq_readin', 'Title line missing from input.', abs(ios)) - ENDIF - ! - ! ... set default values for variables in namelist - ! - tr2_ph = 1.D-12 - eth_rps = 1.D-9 - eth_ns = 1.D-12 - amass(:) = 0.D0 - alpha_mix(:) = 0.D0 - alpha_mix(1) = 0.7D0 - niter_ph = maxter - nmix_ph = 4 - nat_todo = 0 - modenum = 0 - iverbosity = 1234567 - verbosity = 'default' - trans = .TRUE. - lrpa = .FALSE. - lnoloc = .FALSE. - epsil = .FALSE. - zeu = .TRUE. - zue = .FALSE. - fpol = .FALSE. - electron_phonon=' ' - elph_nbnd_min = 1 - elph_nbnd_max = 0 - el_ph_sigma = 0.02 - el_ph_nsigma = 10 - el_ph_ngauss = 1 - lraman = .FALSE. - elop = .FALSE. - max_seconds = 1.E+7_DP - reduce_io = .FALSE. - CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) - IF ( TRIM( outdir ) == ' ' ) outdir = './' - prefix = 'pwscf' - fildyn = 'matdyn' - fildrho = ' ' - fildvscf = ' ' - ldisp = .FALSE. - nq1 = 0 - nq2 = 0 - nq3 = 0 - dek = 1.0d-3 - nogg = .FALSE. - recover = .FALSE. - asr = .FALSE. - start_irr = 0 - last_irr = -1000 - start_q = 1 - last_q =-1000 - ldiag =.FALSE. - lqdir =.FALSE. - qplot =.FALSE. - q_in_band_form=.FALSE. - q2d = .FALSE. - only_init = .FALSE. - only_wfc = .FALSE. - low_directory_check=.FALSE. - search_sym =.TRUE. - nk1 = 0 - nk2 = 0 - nk3 = 0 - k1 = 0 - k2 = 0 - k3 = 0 - ! - drho_star%open = .FALSE. - drho_star%basis = 'modes' - drho_star%pat = .TRUE. - drho_star%ext = 'drho' - CALL get_environment_variable( 'ESPRESSO_FILDRHO_DIR', drho_star%dir) - IF ( TRIM( drho_star%dir ) == ' ' ) & - drho_star%dir = TRIM(outdir)//"/Rotated_DRHO/" - ! - dvscf_star%open = .FALSE. - dvscf_star%basis = 'modes' - dvscf_star%pat = .FALSE. - dvscf_star%ext = 'dvscf' - CALL get_environment_variable( 'ESPRESSO_FILDVSCF_DIR', dvscf_star%dir) - IF ( TRIM( dvscf_star%dir ) == ' ' ) & - dvscf_star%dir = TRIM(outdir)//"/Rotated_DVSCF/" - ! - lshift_q = .false. - read_dns_bare =.false. - d2ns_type = 'full' - ! - ! ... reading the namelist inputph - ! - IF (meta_ionode) THEN - READ( 5, INPUTPH, ERR=30, IOSTAT = ios ) - ! - ! ... iverbosity/verbosity hack - ! - IF ( iverbosity == 1234567 ) THEN - SELECT CASE( trim( verbosity ) ) - CASE( 'debug', 'high', 'medium' ) - iverbosity = 1 - CASE( 'low', 'default', 'minimal' ) - iverbosity = 0 - CASE DEFAULT - iverbosity = 0 - END SELECT - ELSE - ios = 1234567 - END IF - END IF -30 CONTINUE - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - IF ( ios == 1234567 ) THEN - CALL infomsg( 'phq_readin' , & - 'iverbosity is obsolete, use "verbosity" instead' ) - ELSE IF ( ABS(ios) /= 0 ) THEN - CALL errore( 'phq_readin', 'reading inputph namelist', ABS( ios ) ) - END IF - ! diagonalization option - SELECT CASE(TRIM(diagonalization)) - CASE ('david','davidson') - isolve = 0 - CASE ('cg') - isolve = 1 - CASE DEFAULT - CALL errore('phq_readin','diagonalization '//trim(diagonalization)//' not implemented',1) - END SELECT - - ! - ! ... broadcast all input variables - ! - tmp_dir = trimcheck (outdir) - CALL bcast_ph_input ( ) - CALL mp_bcast(nogg, meta_ionode_id, world_comm ) - CALL mp_bcast(q2d, meta_ionode_id, world_comm ) - CALL mp_bcast(q_in_band_form, meta_ionode_id, world_comm ) - ! - drho_star%dir=trimcheck(drho_star%dir) - dvscf_star%dir=trimcheck(dvscf_star%dir) - ! filename for the star must always be automatically generated: - IF(drho_star%ext(1:5)/='auto:') drho_star%ext = 'auto:'//drho_star%ext - IF(dvscf_star%ext(1:5)/='auto:') dvscf_star%ext = 'auto:'//dvscf_star%ext - ! - ! ... Check all namelist variables - ! - IF (tr2_ph <= 0.D0) CALL errore (' phq_readin', ' Wrong tr2_ph ', 1) - IF (eth_rps<= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_rps', 1) - IF (eth_ns <= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_ns ', 1) - - DO iter = 1, maxter - IF (alpha_mix (iter) .LT.0.D0.OR.alpha_mix (iter) .GT.1.D0) CALL & - errore ('phq_readin', ' Wrong alpha_mix ', iter) - ENDDO - IF (niter_ph.LT.1.OR.niter_ph.GT.maxter) CALL errore ('phq_readin', & - ' Wrong niter_ph ', 1) - IF (nmix_ph.LT.1.OR.nmix_ph.GT.5) CALL errore ('phq_readin', ' Wrong & - &nmix_ph ', 1) - IF (iverbosity.NE.0.AND.iverbosity.NE.1) CALL errore ('phq_readin', & - &' Wrong iverbosity ', 1) - IF (fildyn.EQ.' ') CALL errore ('phq_readin', ' Wrong fildyn ', 1) - IF (max_seconds.LT.0.1D0) CALL errore ('phq_readin', ' Wrong max_seconds', 1) - - IF (only_init.AND.only_wfc) CALL errore('phq_readin', & - 'only_init or only_wfc can be .true.', 1) - - IF (modenum < 0) CALL errore ('phq_readin', ' Wrong modenum ', 1) - IF (dek <= 0.d0) CALL errore ( 'phq_readin', ' Wrong dek ', 1) - ! - ! - elph_tetra = 0 - SELECT CASE( trim( electron_phonon ) ) - CASE( 'simple' ) - elph=.true. - elph_mat=.false. - elph_simple=.true. - elph_epa=.false. - CASE( 'epa' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - elph_epa=.true. - CASE( 'Wannier' ) - elph=.true. - elph_mat=.true. - elph_simple=.false. - elph_epa=.false. - auxdvscf=trim(fildvscf) - CASE( 'interpolated' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - elph_epa=.false. - ! YAMBO > - CASE( 'yambo' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - elph_epa=.false. - elph_yambo=.true. - nogg=.true. - auxdvscf=trim(fildvscf) - CASE( 'dvscf' ) - elph=.false. - elph_mat=.false. - elph_simple=.false. - elph_epa=.false. - elph_yambo=.false. - dvscf_yambo=.true. - nogg=.true. - auxdvscf=trim(fildvscf) - ! YAMBO < - CASE( 'lambda_tetra' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - trans = .false. - elph_tetra = 1 - CASE( 'gamma_tetra' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - trans = .false. - elph_tetra = 2 - CASE( 'scdft_input' ) - elph=.true. - elph_mat=.false. - elph_simple=.false. - trans = .false. - elph_tetra = 3 - CASE DEFAULT - elph=.false. - elph_mat=.false. - elph_simple=.false. - elph_epa=.false. - END SELECT - ! YAMBO > - IF (.not.elph_yambo) then - ! YAMBO < - IF (elph.AND.qplot) & - CALL errore('phq_readin', 'qplot and elph not implemented',1) - ENDIF - - ! YAMBO > - IF (.not.elph_yambo.and..not.dvscf_yambo) then - ! YAMBO < - IF (qplot.AND..NOT.ldisp) CALL errore('phq_readin','qplot requires ldisp=.true.',1) - ! - ENDIF - - IF (ldisp.AND.only_init.AND.(.NOT.lqdir)) & - CALL errore('phq_readin', & - 'only_init=.TRUE. requires lqdir=.TRUE. or data are lost',1) - - epsil = epsil .OR. lraman .OR. elop - - IF (modenum /= 0) search_sym=.FALSE. - - IF (elph_mat) THEN - trans=.false. - ELSEIF (.NOT. elph) THEN - trans = trans .OR. ldisp - ENDIF - ! - ! Set default value for fildrho and fildvscf if they are required - IF ( (lraman.OR.elop.OR.drho_star%open) .AND. fildrho == ' ') fildrho = 'drho' - IF ( (elph_mat.OR.dvscf_star%open) .AND. fildvscf == ' ') fildvscf = 'dvscf' - ! - ! We can calculate dielectric, raman or elop tensors and no Born effective - ! charges dF/dE, but we cannot calculate Born effective charges dF/dE - ! without epsil. - ! - IF (zeu) zeu = epsil - ! - ! reads the q point (just if ldisp = .false.) - ! - IF (meta_ionode) THEN - ios = 0 - IF (qplot) THEN - READ (5, *, iostat = ios) nqaux - ELSE - IF (.NOT. ldisp) READ (5, *, iostat = ios) (xq (ipol), ipol = 1, 3) - ENDIF - END IF - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore ('phq_readin', 'reading xq', ABS (ios) ) - IF (qplot) THEN - CALL mp_bcast(nqaux, meta_ionode_id, world_comm ) - ALLOCATE(xqaux(3,nqaux)) - ALLOCATE(wqaux(nqaux)) - IF (meta_ionode) THEN - DO iq=1, nqaux - READ (5, *, iostat = ios) (xqaux (ipol,iq), ipol = 1, 3), wqaux(iq) - ENDDO - ENDIF - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore ('phq_readin', 'reading xq', ABS (ios) ) - CALL mp_bcast(xqaux, meta_ionode_id, world_comm ) - CALL mp_bcast(wqaux, meta_ionode_id, world_comm ) - ELSE - CALL mp_bcast(xq, meta_ionode_id, world_comm ) - ENDIF - - IF (.NOT.ldisp) THEN - lgamma = xq (1) .EQ.0.D0.AND.xq (2) .EQ.0.D0.AND.xq (3) .EQ.0.D0 - IF ( (epsil.OR.zue.or.lraman.or.elop) .AND..NOT.lgamma) & - CALL errore ('phq_readin', 'gamma is needed for elec.field', 1) - ENDIF - IF (zue.AND..NOT.trans) CALL errore ('phq_readin', 'trans must be & - &.t. for Zue calc.', 1) - - IF (trans.AND.(lrpa.OR.lnoloc)) CALL errore('phq_readin', & - 'only dielectric constant with lrpa or lnoloc',1) - IF (lrpa.or.lnoloc) THEN - zeu=.FALSE. - lraman=.FALSE. - elop = .FALSE. - ENDIF - ! - ! reads the frequencies ( just if fpol = .true. ) - ! - IF ( fpol ) THEN - IF ( .NOT. epsil) CALL errore ('phq_readin', & - 'fpol=.TRUE. needs epsil=.TRUE.', 1 ) - nfs=0 - IF (meta_ionode) THEN - READ (5, *, iostat = ios) card - IF ( TRIM(card)=='FREQUENCIES'.OR. & - TRIM(card)=='frequencies'.OR. & - TRIM(card)=='Frequencies') THEN - READ (5, *, iostat = ios) nfs - ENDIF - ENDIF - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore ('phq_readin', 'reading number of FREQUENCIES', ABS(ios) ) - CALL mp_bcast(nfs, meta_ionode_id, world_comm ) - if (nfs < 1) call errore('phq_readin','Too few frequencies',1) - ALLOCATE(fiu(nfs)) - IF (meta_ionode) THEN - IF ( TRIM(card) == 'FREQUENCIES' .OR. & - TRIM(card) == 'frequencies' .OR. & - TRIM(card) == 'Frequencies' ) THEN - DO i = 1, nfs - READ (5, *, iostat = ios) fiu(i) - END DO - END IF - END IF - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore ('phq_readin', 'reading FREQUENCIES card', ABS(ios) ) - CALL mp_bcast(fiu, meta_ionode_id, world_comm ) - ELSE - nfs=1 - ALLOCATE(fiu(1)) - fiu=0.0_DP - END IF - ! - ! - ! Here we finished the reading of the input file. - ! Now allocate space for pwscf variables, read and check them. - ! - ! amass will be read again from file: - ! save its content in auxiliary variables - ! - ALLOCATE ( amass_input( SIZE(amass) ) ) - amass_input(:)= amass(:) - ! - tmp_dir_save=tmp_dir - tmp_dir_ph= trimcheck( TRIM (tmp_dir) // '_ph' // int_to_char(my_image_id) ) - CALL check_tempdir ( tmp_dir_ph, exst, parallelfs ) - tmp_dir_phq=tmp_dir_ph - - ext_restart=.FALSE. - ext_recover=.FALSE. - rec_code_read=-1000 - IF (recover) THEN -! -! With a recover run we read here the mesh of q points, the current iq, -! and the current frequency -! - CALL ph_readfile('init', 0, 0, ierr) - CALL ph_readfile('status_ph', 0, 0, ierr1) -! -! If some error occured here, we cannot recover the run -! - IF (ierr /= 0 .OR. ierr1 /= 0 ) THEN - write(stdout,'(5x,"Run is not recoverable starting from scratch")') - recover=.FALSE. - goto 1001 - ENDIF -! -! We check if the bands and the information on the pw run are in the directory -! written by the phonon code for the current q point. If the file exists -! we read from there, otherwise use the information in tmp_dir. -! - IF (lqdir) THEN - tmp_dir_phq= trimcheck ( TRIM(tmp_dir_ph) // TRIM(prefix) // & - & '.q_' // int_to_char(current_iq) ) - CALL check_restart_recover(ext_recover, ext_restart) - IF (.NOT.ext_recover.AND..NOT.ext_restart) tmp_dir_phq=tmp_dir_ph - ENDIF - ! - filename=TRIM(tmp_dir_phq)//TRIM(prefix)//postfix//xmlpun_schema - IF (ionode) inquire (file =TRIM(filename), exist = exst) - ! - CALL mp_bcast( exst, ionode_id, intra_image_comm ) - ! - ! If this file exist we use it to recover the pw.x informations - ! - IF (exst) tmp_dir=tmp_dir_phq - u_from_file=.true. - ENDIF -1001 CONTINUE - - IF (qplot.AND..NOT.recover) THEN - IF (q2d) THEN - nqs=wqaux(2)*wqaux(3) - ALLOCATE(x_q(3,nqs)) - ALLOCATE(wq(nqs)) - CALL generate_k_in_plane(nqaux, xqaux, wqaux, x_q, wq, nqs) - ELSEIF (q_in_band_form) THEN - nqs=SUM(wqaux(1:nqaux-1))+1 - DO i=1,nqaux-1 - IF (wqaux(i)==0) nqs=nqs+1 - ENDDO - ALLOCATE(x_q(3,nqs)) - ALLOCATE(wq(nqs)) - CALL generate_k_along_lines(nqaux, xqaux, wqaux, x_q, wq, nqs) - ELSE - nqs=nqaux - ALLOCATE(x_q(3,nqs)) - ALLOCATE(wq(nqs)) - wq(:)=wqaux(:) - x_q(:,1:nqs)=xqaux(:,1:nqs) - ENDIF - DEALLOCATE(xqaux) - DEALLOCATE(wqaux) - ALLOCATE(lgamma_iq(nqs)) - DO iq=1, nqs - lgamma_iq(iq)= ( ABS(x_q(1,iq)) .LT. 1.0e-10_dp ) .AND. & - ( ABS(x_q(2,iq)) .LT. 1.0e-10_dp ) .AND. & - ( ABS(x_q(3,iq)) .LT. 1.0e-10_dp ) - ENDDO - WRITE(stdout, '(//5x,"Dynamical matrices for q-points given in input")') - WRITE(stdout, '(5x,"(",i4,"q-points):")') nqs - WRITE(stdout, '(5x," N xq(1) xq(2) xq(3) " )') - DO iq = 1, nqs - WRITE(stdout, '(5x,i3, 3f14.9)') iq, x_q(1,iq), x_q(2,iq), x_q(3,iq) - END DO - ENDIF - ! - ! DFPT+U: the occupation matrix ns is read via read_file - ! - CALL read_file ( ) - - magnetic_sym=noncolin .AND. domag - ! - ! init_start_grid returns .true. if a new k-point grid is set from values - ! read from input (this happens if nk1*nk2*nk3 > 0; otherwise reset_grid - ! returns .false., leaves the current values, read in read_file, unchanged) - ! - newgrid = reset_grid (nk1, nk2, nk3, k1, k2, k3) - ! - tmp_dir=tmp_dir_save - ! - IF (modenum > 3*nat) CALL errore ('phq_readin', ' Wrong modenum ', 2) - - IF (gamma_only) CALL errore('phq_readin',& - 'cannot start from pw.x data file using Gamma-point tricks',1) - - IF (lda_plus_u) THEN - ! - WRITE(stdout,'(/5x,a)') "Phonon calculation with DFPT+U; please cite" - WRITE(stdout,'(5x,a)') "A. Floris et al., Phys. Rev. B 84, 161102(R) (2011)" - WRITE(stdout,'(5x,a)') "in publications or presentations arising from this work." - ! - IF (U_projection.NE."atomic") CALL errore("phq_readin", & - " The phonon code for this U_projection_type is not implemented",1) - IF (lda_plus_u_kind.NE.0) CALL errore("phq_readin", & - " The phonon code for this lda_plus_u_kind is not implemented",1) - IF (elph) CALL errore("phq_readin", & - " Electron-phonon with Hubbard U is not supported",1) - IF (lraman) CALL errore("phq_readin", & - " The phonon code with Raman and Hubbard U is not implemented",1) - ! - ENDIF - - IF (ts_vdw) CALL errore('phq_readin',& - 'The phonon code with TS-VdW is not yet available',1) - - IF (ldftd3) CALL errore('phq_readin',& - 'The phonon code with Grimme''s DFT-D3 is not yet available',1) - - IF ( dft_is_meta() ) CALL errore('phq_readin',& - 'The phonon code with meta-GGA functionals is not yet available',1) - - IF ( dft_is_hybrid() ) CALL errore('phq_readin',& - 'The phonon code with hybrid functionals is not yet available',1) - - IF (okpaw.and.(lraman.or.elop)) CALL errore('phq_readin',& - 'The phonon code with paw and raman or elop is not yet available',1) - - IF (okpaw.and.noncolin.and.domag) CALL errore('phq_readin',& - 'The phonon code with paw and domag is not available yet',1) - - IF (okvan.and.(lraman.or.elop)) CALL errore('phq_readin',& - 'The phonon code with US-PP and raman or elop not yet available',1) - - IF (noncolin.and.(lraman.or.elop)) CALL errore('phq_readin', & - 'lraman, elop, and noncolin not programed',1) - - IF (lmovecell) CALL errore('phq_readin', & - 'The phonon code is not working after vc-relax',1) - - IF (reduce_io) io_level=1 - - if(elph_mat.and.fildvscf.eq.' ') call errore('phq_readin',& - 'el-ph with wannier requires fildvscf',1) - - IF(elph_mat.and.npool.ne.1) call errore('phq_readin',& - 'el-ph with wannier : pools not implemented',1) - - IF(elph.and.nimage>1) call errore('phq_readin',& - 'el-ph with images not implemented',1) - - IF( fildvscf /= ' ' .and. nimage > 1 ) call errore('phq_readin',& - 'saving dvscf to file images not implemented',1) - - IF (elph.OR.fildvscf /= ' ') lqdir=.TRUE. - - IF(dvscf_star%open.and.nimage>1) CALL errore('phq_readin',& - 'dvscf_star with image parallelization is not yet available',1) - IF(drho_star%open.and.nimage>1) CALL errore('phq_readin',& - 'drho_star with image parallelization is not yet available',1) - - IF (lda_plus_u .AND. read_dns_bare .AND. ldisp) lqdir=.TRUE. - - IF (.NOT.ldisp) lqdir=.FALSE. - - IF (i_cons /= 0) & - CALL errore('phq_readin',& - 'The phonon code with constrained magnetization is not yet available',1) - - IF (two_fermi_energies .AND. (ltetra .OR. lgauss)) & - CALL errore('phq_readin',& - 'The phonon code with two fermi energies is not available for metals',1) - - IF (tqr) CALL errore('phq_readin',& - 'The phonon code with Q in real space not available',1) - - IF (start_irr < 0 ) CALL errore('phq_readin', 'wrong start_irr',1) - ! - IF (start_q <= 0 ) CALL errore('phq_readin', 'wrong start_q',1) - ! - ! the dynamical matrix is written in xml format if fildyn ends in - ! .xml or in the noncollinear case. - ! - xmldyn=has_xml(fildyn) - !IF (noncolin) xmldyn=.TRUE. - ! - ! If a band structure calculation needs to be done do not open a file - ! for k point - ! - restart = recover - ! - ! set masses to values read from input, if available; - ! leave values read from file otherwise - ! - DO it = 1, ntyp - IF (amass_input(it) < 0.0_DP) amass_input(it)= & - atom_weight(atomic_number(TRIM(atm(it)))) - IF (amass_input(it) > 0.D0) amass(it) = amass_input(it) - IF (amass(it) <= 0.D0) CALL errore ('phq_readin', 'Wrong masses', it) - ENDDO - DEALLOCATE (amass_input) - lgamma_gamma=.FALSE. - IF (.NOT.ldisp) THEN - IF (nkstot==1.OR.(nkstot==2.AND.nspin==2)) THEN - lgamma_gamma=(lgamma.AND.(ABS(xk(1,1))<1.D-12) & - .AND.(ABS(xk(2,1))<1.D-12) & - .AND.(ABS(xk(3,1))<1.D-12) ) - ENDIF - IF (nogg) lgamma_gamma=.FALSE. - IF ((nat_todo /= 0) .and. lgamma_gamma) CALL errore( & - 'phq_readin', 'gamma_gamma tricks with nat_todo & - & not available. Use nogg=.true.', 1) - ! - IF (nimage > 1 .AND. lgamma_gamma) CALL errore( & - 'phq_readin','gamma_gamma tricks with images not implemented',1) - IF (lgamma) THEN - nksq = nks - ELSE - nksq = nks / 2 - ENDIF - ENDIF - IF (lgamma_gamma.AND.ldiag) CALL errore('phq_readin','incompatible flags',1) - IF (lgamma_gamma.AND.elph ) CALL errore('phq_readin','lgamma_gamma and electron-phonon are incompatible',1) - ! - IF (tfixed_occ) & - CALL errore('phq_readin','phonon with arbitrary occupations not tested',1) - ! - !YAMBO > - ! DBSP - !IF (elph.AND..NOT.(lgauss .or. ltetra).and..NOT.elph_yambo) CALL errore ('phq_readin', 'Electron-& - ! &phonon only for metals', 1) - !YAMBO < - IF (elph.AND.fildvscf.EQ.' ') CALL errore ('phq_readin', 'El-ph needs & - &a DeltaVscf file', 1) - ! There might be other variables in the input file which describe - ! partial computation of the dynamical matrix. Read them here - ! - CALL allocate_part ( nat ) - ! - IF ( nat_todo < 0 .OR. nat_todo > nat ) & - CALL errore ('phq_readin', 'nat_todo is wrong', 1) - IF (nat_todo.NE.0) THEN - IF (meta_ionode) READ (5, *, iostat = ios) (atomo (na), na = 1, nat_todo) - CALL mp_bcast(ios, meta_ionode_id, world_comm ) - CALL errore ('phq_readin', 'reading atoms', ABS (ios) ) - CALL mp_bcast(atomo, meta_ionode_id, world_comm ) - ENDIF - nat_todo_input=nat_todo - - IF (epsil.AND.(lgauss .OR. ltetra)) & - CALL errore ('phq_readin', 'no elec. field with metals', 1) - IF (modenum > 0) THEN - IF ( ldisp ) & - CALL errore('phq_readin','Dispersion calculation and & - & single mode calculation not possibile !',1) - nat_todo = 0 - ENDIF - ! - ! Tetrahedron method for DFPT and El-Ph - ! - IF(ltetra .AND. tetra_type == 0) & - & CALL errore ('phq_readin', 'DFPT with the Blochl correction is not implemented', 1) - ! - IF(.NOT. ltetra .AND. elph_tetra /= 0) & - & CALL errore ('phq_readin', '"electron_phonon" and "occupation" are inconsistent.', 1) - ! - IF (modenum > 0 .OR. lraman ) lgamma_gamma=.FALSE. - IF (.NOT.lgamma_gamma) asr=.FALSE. - ! - IF ((ldisp.AND..NOT.qplot) .AND. & - (nq1 .LE. 0 .OR. nq2 .LE. 0 .OR. nq3 .LE. 0)) & - CALL errore('phq_readin','nq1, nq2, and nq3 must be greater than 0',1) - - CALL save_ph_input_variables() - ! - RETURN - ! -END SUBROUTINE phq_readin diff --git a/test-suite/not_epw_comp/supp.f90 b/test-suite/not_epw_comp/supp.f90 index f4e0f4095..df164cf0d 100644 --- a/test-suite/not_epw_comp/supp.f90 +++ b/test-suite/not_epw_comp/supp.f90 @@ -97,37 +97,37 @@ SUBROUTINE elphsum2( ) et2 = et ELSE IF (found .AND. (me_bgrp == root_bgrp) .AND. (.NOT. ionode)) THEN - CALL MPI_SEND(ik1, 1, MPI_INTEGER, root_image, me_image, intra_image_comm, ierr) - CALL MPI_SEND(ikk, 1, MPI_INTEGER, root_image, me_image + nproc_image, & + CALL MPI_SEND(ik1, 1, MPI_INTEGER, root_image, 0, intra_image_comm, ierr) + CALL MPI_SEND(ikk, 1, MPI_INTEGER, root_image, nproc_image, & intra_image_comm, ierr) - CALL MPI_SEND(ikq, 1, MPI_INTEGER, root_image, me_image + 2 * nproc_image, & + CALL MPI_SEND(ikq, 1, MPI_INTEGER, root_image, 2 * nproc_image, & intra_image_comm, ierr) - CALL MPI_SEND(nksq, 1, MPI_INTEGER, root_image, me_image + 3 * nproc_image, & + CALL MPI_SEND(nksq, 1, MPI_INTEGER, root_image, 3 * nproc_image, & intra_image_comm, ierr) CALL MPI_SEND(el_ph_mat, nbnd * nbnd * nksq * 3 * nat, MPI_DOUBLE_COMPLEX, & - root_image, me_image + 4 * nproc_image, intra_image_comm, ierr) - CALL MPI_SEND(SIZE(et, 2), 1, MPI_INTEGER, root_image, me_image + 5 * nproc_image, & + root_image, 4 * nproc_image, intra_image_comm, ierr) + CALL MPI_SEND(SIZE(et, 2), 1, MPI_INTEGER, root_image, 5 * nproc_image, & intra_image_comm, ierr) CALL MPI_SEND(et, nbnd * SIZE(et, 2), MPI_DOUBLE_PRECISION, root_image, & - me_image + 6 * nproc_image, intra_image_comm, ierr) + 6 * nproc_image, intra_image_comm, ierr) ELSEIF (ionode) THEN IF (.NOT. found) THEN - CALL MPI_RECV(ik2, 1, MPI_INTEGER, MPI_ANY_SOURCE, me_image, intra_image_comm, & + CALL MPI_RECV(ik2, 1, MPI_INTEGER, MPI_ANY_SOURCE, 0, intra_image_comm, & istatus, ierr ) - CALL MPI_RECV(ikk2, 1, MPI_INTEGER, MPI_ANY_SOURCE, me_image + nproc_image, & + CALL MPI_RECV(ikk2, 1, MPI_INTEGER, MPI_ANY_SOURCE, nproc_image, & intra_image_comm, istatus, ierr ) - CALL MPI_RECV(ikq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, me_image + 2 * nproc_image, & + CALL MPI_RECV(ikq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, 2 * nproc_image, & intra_image_comm, istatus, ierr ) - CALL MPI_RECV(nksq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, me_image + 3 * nproc_image, & + CALL MPI_RECV(nksq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, 3 * nproc_image, & intra_image_comm, istatus, ierr ) ALLOCATE(el_ph_mat2(nbnd, nbnd, nksq2, 3 * nat)) CALL MPI_RECV(el_ph_mat2, nbnd * nbnd * nksq2 * 3 * nat, MPI_DOUBLE_COMPLEX, & - MPI_ANY_SOURCE, me_image + 4 * nproc_image, intra_image_comm, istatus, ierr ) - CALL MPI_RECV(nkq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, me_image + 5 * nproc_image, & + MPI_ANY_SOURCE, 4 * nproc_image, intra_image_comm, istatus, ierr ) + CALL MPI_RECV(nkq2, 1, MPI_INTEGER, MPI_ANY_SOURCE, 5 * nproc_image, & intra_image_comm, istatus, ierr ) ALLOCATE(et2(nbnd, nkq2)) CALL MPI_RECV(et2, nbnd * nkq2, MPI_DOUBLE_PRECISION, MPI_ANY_SOURCE, & - me_image + 6 * nproc_image, intra_image_comm, istatus, ierr ) + 6 * nproc_image, intra_image_comm, istatus, ierr ) ELSE ik2 = ik1 ikk2 = ikk diff --git a/test-suite/not_epw_comp/symdyn_munu.f90 b/test-suite/not_epw_comp/symdyn_munu.f90 deleted file mode 100644 index b722f3f0c..000000000 --- a/test-suite/not_epw_comp/symdyn_munu.f90 +++ /dev/null @@ -1,193 +0,0 @@ -! -! Copyright (C) 2001-2012 PWSCF group -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! -!----------------------------------------------------------------------- -subroutine symdyn_munu_new (dyn, u, xq, s, invs, rtau, irt, at, & - bg, nsymq, nat, irotmq, minus_q) - !----------------------------------------------------------------------- - ! - ! This routine symmetrize the dynamical matrix written in the basis - ! of the modes - ! - ! - USE kinds, only : DP - implicit none - integer :: nat, s (3, 3, 48), irt (48, nat), invs (48), & - nsymq, irotmq - ! input: the number of atoms - ! input: the symmetry matrices - ! input: the rotated of each atom - ! input: the small group of q - ! input: the inverse of each matrix - ! input: the order of the small gro - ! input: the symmetry q -> -q+G - - real(DP) :: xq (3), rtau (3, 48, nat), at (3, 3), bg (3, 3) - ! input: the coordinates of q - ! input: the R associated at each r - ! input: direct lattice vectors - ! input: reciprocal lattice vectors - - logical :: minus_q - ! input: if true symmetry sends q-> - - complex(DP) :: dyn (3 * nat, 3 * nat), u (3 * nat, 3 * nat) - ! inp/out: matrix to symmetrize - ! input: the patterns - - integer :: i, j, icart, jcart, na, nb, mu, nu - ! counter on modes - ! counter on modes - ! counter on cartesian coordinates - ! counter on cartesian coordinates - ! counter on atoms - ! counter on atoms - ! counter on modes - ! counter on modes - - complex(DP) :: work, phi (3, 3, nat, nat) - ! auxiliary variable - ! the dynamical matrix - ! - ! First we transform in the cartesian coordinates - ! - CALL dyn_pattern_to_cart(nat, u, dyn, phi) - ! - ! Then we transform to the crystal axis - ! - do na = 1, nat - do nb = 1, nat - call trntnsc (phi (1, 1, na, nb), at, bg, - 1) - enddo - enddo - ! - ! And we symmetrize in this basis - ! - call symdynph_gq_new (xq, phi, s, invs, rtau, irt, nsymq, nat, & - irotmq, minus_q) - ! - ! Back to cartesian coordinates - ! - do na = 1, nat - do nb = 1, nat - call trntnsc (phi (1, 1, na, nb), at, bg, + 1) - enddo - enddo - ! - ! rewrite the dynamical matrix on the array dyn with dimension 3nat x 3nat - ! - CALL compact_dyn(nat, dyn, phi) - - return -end subroutine symdyn_munu_new -! -!----------------------------------------------------------------------- -subroutine symdyn_munu2 (dyn, u, xq, s, invs, rtau, irt, at, & - bg, nsymq, nat, irotmq, minus_q) - !----------------------------------------------------------------------- - ! - ! This routine symmetrize the dynamical matrix written in the basis - ! of the modes - ! - ! -!#include "f_defs.h" - USE kinds, only : DP - implicit none - integer :: nat, s (3, 3, 48), irt (48, nat), invs (48), & - nsymq, irotmq - ! input: the number of atoms - ! input: the symmetry matrices - ! input: the rotated of each atom - ! input: the small group of q - ! input: the inverse of each matrix - ! input: the order of the small gro - ! input: the symmetry q -> -q+G - - real(kind=DP) :: xq (3), rtau (3, 48, nat), at (3, 3), bg (3, 3) - ! input: the coordinates of q - ! input: the R associated at each r - ! input: direct lattice vectors - ! input: reciprocal lattice vectors - - logical :: minus_q - ! input: if true symmetry sends q-> - - complex(kind=DP) :: dyn (3 * nat, 3 * nat), u (3 * nat, 3 * nat) - ! inp/out: matrix to symmetrize - ! input: the patterns - - integer :: i, j, icart, jcart, na, nb, mu, nu - ! counter on modes - ! counter on modes - ! counter on cartesian coordinates - ! counter on cartesian coordinates - ! counter on atoms - ! counter on atoms - ! counter on modes - ! counter on modes - - complex(kind=DP) :: work, wrk (3, 3), phi (3, 3, nat, nat) - ! auxiliary variable - ! auxiliary variable - ! the dynamical matrix - ! - ! First we transform in the cartesian coordinates - ! - do i = 1, 3 * nat - na = (i - 1) / 3 + 1 - icart = i - 3 * (na - 1) - do j = 1, 3 * nat - nb = (j - 1) / 3 + 1 - jcart = j - 3 * (nb - 1) - work = (0.d0, 0.d0) - do mu = 1, 3 * nat - do nu = 1, 3 * nat - work = work + u (i, mu) * dyn (mu, nu) * conjg (u (j, nu) ) - enddo - enddo - phi (icart, jcart, na, nb) = work - enddo - enddo - ! -!@ ! Then we transform to the crystal axis -!@ ! -!@ do na = 1, nat -!@ do nb = 1, nat -!@ call trntnsc (phi (1, 1, na, nb), at, bg, - 1) -!@ enddo -!@ enddo -!@ ! -!@ ! And we symmetrize in this basis -!@ ! -!@ call symdynph_gq (xq, phi, s, invs, rtau, irt, irgq, nsymq, nat, & -!@ irotmq, minus_q) -!@ ! -!@ ! Back to cartesian coordinates -!@ ! -!@ do na = 1, nat -!@ do nb = 1, nat -!@ call trntnsc (phi (1, 1, na, nb), at, bg, + 1) -!@ enddo -!@ enddo - ! - ! rewrite the dynamical matrix on the array dyn with dimension 3nat x 3 - ! - do i = 1, 3 * nat - na = (i - 1) / 3 + 1 - icart = i - 3 * (na - 1) - do j = 1, 3 * nat - nb = (j - 1) / 3 + 1 - jcart = j - 3 * (nb - 1) - dyn (i, j) = phi (icart, jcart, na, nb) - enddo - - enddo - return -end subroutine symdyn_munu2 - - -