diff --git a/LR_Modules/apply_trev.f90 b/LR_Modules/apply_trev.f90 new file mode 100644 index 000000000..722d17cf0 --- /dev/null +++ b/LR_Modules/apply_trev.f90 @@ -0,0 +1,56 @@ +! +! Copyright (C) 2018 Andrea Dal Corso and Andrea Urru +! 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 apply_trev(evc, ikk_evc, ikk_tevc) +! +!! This routine applies the time reversal operator to the wavefunctions +!! \(\text{evc}\) at the k point \(\text{ikk_evc}\) and puts the output +!! in \(\text{evc}\) with the order of G vectors of \(\text{ikk_tevc}\). +! +USE kinds, ONLY : DP +USE wvfct, ONLY : nbnd, npwx +USE klist, ONLY : ngk, igk_k +USE fft_base, ONLY: dffts +USE fft_interfaces, ONLY: invfft, fwfft +USE noncollin_module, ONLY : npol + +IMPLICIT NONE +INTEGER :: ikk_evc, ikk_tevc +COMPLEX(DP) :: evc(npwx*npol,nbnd) +COMPLEX(DP), ALLOCATABLE :: aux2(:,:) + +INTEGER :: npw, npwt, ibnd, ig + + +npw=ngk(ikk_evc) +npwt=ngk(ikk_tevc) + +ALLOCATE(aux2(dffts%nnr,2)) + +DO ibnd=1,nbnd + aux2(:,:) = (0.D0, 0.D0) + DO ig = 1, npw + aux2 (dffts%nl (igk_k (ig, ikk_evc) ), 1 ) = evc (ig, ibnd) + aux2 (dffts%nl (igk_k (ig, ikk_evc) ), 2 ) = evc (npwx+ig, ibnd) + END DO + CALL invfft ('Wave', aux2(:,1), dffts) + CALL invfft ('Wave', aux2(:,2), dffts) + aux2=CONJG(aux2) + CALL fwfft ('Wave', aux2(:,1), dffts) + CALL fwfft ('Wave', aux2(:,2), dffts) + evc(:,ibnd)=(0.0_DP,0.0_DP) + DO ig = 1, npwt + evc(ig,ibnd)=-aux2(dffts%nl(igk_k (ig, ikk_tevc)),2) + evc(ig+npwx,ibnd)=aux2(dffts%nl(igk_k (ig, ikk_tevc)),1) + END DO +END DO + +DEALLOCATE(aux2) + +RETURN +END SUBROUTINE apply_trev diff --git a/LR_Modules/c_bands_ph.f90 b/LR_Modules/c_bands_ph.f90 new file mode 100644 index 000000000..9ca2ede27 --- /dev/null +++ b/LR_Modules/c_bands_ph.f90 @@ -0,0 +1,174 @@ +! +! 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 . +! +! +!---------------------------------------------------------------------------- +! Questa e' una copia di c_bands_nscf intesa per un confronto con +! thermo_pw. +! +! +SUBROUTINE c_bands_nscf_ph( ) + !---------------------------------------------------------------------------- + !! Driver routine for Hamiltonian diagonalization routines + !! specialized to non-self-consistent calculations (no electric field). + ! + USE kinds, ONLY : DP + USE io_global, ONLY : stdout + USE io_files, ONLY : iunhub, iunwfc, nwordwfc, nwordwfcU + USE buffers, ONLY : get_buffer, save_buffer, close_buffer, open_buffer + USE basis, ONLY : starting_wfc + USE klist, ONLY : nkstot, nks, xk, ngk, igk_k + USE uspp, ONLY : vkb, nkb + USE gvect, ONLY : g + USE wvfct, ONLY : et, nbnd, npwx, current_k + USE control_lr, ONLY : lgamma + USE control_flags, ONLY : ethr, restart, isolve, io_level, iverbosity + USE ldaU, ONLY : lda_plus_u, Hubbard_projectors, wfcU, lda_plus_u_kind + USE lsda_mod, ONLY : current_spin, lsda, isk + USE wavefunctions, ONLY : evc + USE mp_pools, ONLY : npool, kunit, inter_pool_comm + USE mp, ONLY : mp_sum + USE check_stop, ONLY : check_stop_now + USE noncollin_module, ONLY : noncolin, npol, domag + ! USE save_ph, ONLY : tmp_dir_save + USE io_files, ONLY : tmp_dir, prefix + USE uspp_init, ONLY : init_us_2 + ! + IMPLICIT NONE + ! + REAL(DP) :: avg_iter, ethr_ + ! average number of H*psi products + INTEGER :: ik_, ik, nkdum, ios + ! ik_: k-point already done in a previous run + ! ik : counter on k points + LOGICAL :: exst, exst_mem + ! + REAL(DP), EXTERNAL :: get_clock + ! + CALL start_clock( 'c_bands' ) + ! + ik_ = 0 + avg_iter = 0.D0 + IF ( restart ) CALL restart_in_cbands(ik_, ethr, avg_iter, et ) + ! + ! ... If restarting, calculated wavefunctions have to be read from file + ! + DO ik = 1, ik_ + CALL get_buffer ( evc, nwordwfc, iunwfc, ik ) + END DO + ! + IF ( isolve == 0 ) THEN + WRITE( stdout, '(5X,"Davidson diagonalization with overlap")' ) + ELSE IF ( isolve == 1 ) THEN + WRITE( stdout, '(5X,"CG style diagonalization")') + ELSE IF ( isolve == 2 ) THEN + WRITE( stdout, '(5X,"PPCG style diagonalization")') + ELSE + CALL errore ( 'c_bands', 'invalid type of diagonalization', isolve) + END IF + ! + ! ... For each k point (except those already calculated if restarting) + ! ... diagonalizes the hamiltonian + ! + k_loop: DO ik = ik_+1, nks + ! + ! ... Set k-point, spin, kinetic energy, needed by Hpsi + ! + current_k = ik + ! + ! ---------- LUCA (spawoc) ----------------- + IF (lda_plus_u .AND. lda_plus_u_kind.EQ.2) CALL phase_factor(ik) + ! ------------------------------------ + ! + IF ( lsda ) current_spin = isk(ik) + call g2_kin( ik ) + ! + ! ... More stuff needed by the hamiltonian: nonlocal projectors + ! + IF ( nkb > 0 ) CALL init_us_2( ngk(ik), igk_k(1,ik), xk(1,ik), vkb ) + ! + ! ... Needed for LDA+U + ! + IF ( nks > 1 .AND. lda_plus_u .AND. (Hubbard_projectors .NE. 'pseudo') ) & + CALL get_buffer ( wfcU, nwordwfcU, iunhub, ik ) + ! + ! ... calculate starting wavefunctions + ! + IF ( iverbosity > 0 ) WRITE( stdout, 9001 ) ik + ! + IF ( TRIM(starting_wfc) == 'file' ) THEN + ! + CALL get_buffer ( evc, nwordwfc, iunwfc, ik ) + ! + ELSE + ! + CALL init_wfc ( ik ) + ! + END IF + ! + ! ... diagonalization of bands for k-point ik + ! + call diag_bands ( 1, ik, avg_iter ) + ! + ! In the noncolinear magnetic case we have k, k+q, -k -k-q and + ! to the last two wavefunctions we must apply t_rev. + ! When lgamma is true we have only k and -k + ! + IF (noncolin.AND.domag) THEN + IF (lgamma.AND. MOD(ik,2)==0) THEN + CALL apply_trev(evc, ik, ik-1) + ELSEIF (.NOT.lgamma.AND.(MOD(ik,4)==3.OR.MOD(ik,4)==0)) THEN + CALL apply_trev(evc, ik, ik-2) + ENDIF + ENDIF + ! + ! ... save wave-functions (unless disabled in input) + ! + IF ( io_level > -1 ) CALL save_buffer ( evc, nwordwfc, iunwfc, ik ) + ! + ! ... beware: with pools, if the number of k-points on different + ! ... pools differs, make sure that all processors are still in + ! ... the loop on k-points before checking for stop condition + ! + nkdum = kunit * ( nkstot / kunit / npool ) + IF (ik .le. nkdum) THEN + ! + ! ... stop requested by user: save restart information, + ! ... save wavefunctions to file + ! + IF (check_stop_now()) THEN + CALL save_in_cbands(ik, ethr, avg_iter, et ) + RETURN + END IF + ENDIF + ! + ! report about timing + ! + IF ( iverbosity > 0 ) THEN + WRITE( stdout, 9000 ) get_clock( 'PWSCF' ) + FLUSH( stdout ) + ENDIF + ! + END DO k_loop + ! + CALL mp_sum( avg_iter, inter_pool_comm ) + avg_iter = avg_iter / nkstot + ! + WRITE( stdout, '(/,5X,"ethr = ",1PE9.2,", avg # of iterations =",0PF5.1)' ) & + ethr, avg_iter + ! + CALL stop_clock( 'c_bands' ) + ! + RETURN + ! + ! formats + ! +9001 FORMAT(/' Computing kpt #: ',I5 ) +9000 FORMAT( ' total cpu time spent up to now is ',F10.1,' secs' ) + ! +END SUBROUTINE c_bands_nscf_ph + diff --git a/LR_Modules/non_scf_ph.f90 b/LR_Modules/non_scf_ph.f90 new file mode 100644 index 000000000..383d4a896 --- /dev/null +++ b/LR_Modules/non_scf_ph.f90 @@ -0,0 +1,112 @@ +! +! 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 non_scf_ph ( ) + !----------------------------------------------------------------------- + !! Diagonalization of the KS hamiltonian in the non-scf case. + ! + USE kinds, ONLY : DP + USE bp, ONLY : lelfield, lberry, lorbm + USE check_stop, ONLY : stopped_by_user + USE control_flags, ONLY : io_level, conv_elec, lbands + USE ener, ONLY : ef + USE io_global, ONLY : stdout, ionode + USE io_files, ONLY : iunwfc, nwordwfc, iunefield + USE buffers, ONLY : save_buffer + USE klist, ONLY : xk, wk, nks, nkstot + USE lsda_mod, ONLY : lsda, nspin + USE wvfct, ONLY : nbnd, et, npwx + USE wavefunctions, ONLY : evc + ! + IMPLICIT NONE + ! + ! ... local variables + ! + INTEGER :: iter, i + REAL(DP), EXTERNAL :: get_clock + ! + ! + CALL start_clock( 'electrons' ) + iter = 1 + ! + WRITE( stdout, 9002 ) + FLUSH( stdout ) + ! + IF ( lelfield) THEN + ! + CALL c_bands_efield ( iter ) + ! + ELSE + ! + CALL c_bands_nscf_ph ( ) + ! + END IF + ! + ! ... check if calculation was stopped in c_bands + ! + IF ( stopped_by_user ) THEN + conv_elec=.FALSE. + RETURN + END IF + ! + ! ... xk, wk, isk, et, wg are distributed across pools; + ! ... the first node has a complete copy of xk, wk, isk, + ! ... while eigenvalues et and weights wg must be + ! ... explicitly collected to the first node + ! ... this is done here for et, in weights () for wg + ! + CALL poolrecover( et, nbnd, nkstot, nks ) + ! + ! ... calculate weights of Kohn-Sham orbitals (only weights, not Ef, + ! ... for a "bands" calculation where Ef is read from data file) + ! ... may be needed in further calculations such as phonon + ! + IF ( lbands ) THEN + CALL weights_only ( ) + ELSE + CALL weights ( ) + END IF + ! + ! ... Note that if you want to use more k-points for the phonon + ! ... calculation then those needed for self-consistency, you can, + ! ... by performing a scf with less k-points, followed by a non-scf + ! ... one with additional k-points, whose weight on input is set to zero + ! + WRITE( stdout, 9000 ) get_clock( 'PWSCF' ) + ! + WRITE( stdout, 9102 ) + ! + ! ... write band eigenvalues (conv_elec is used in print_ks_energies) + ! + conv_elec = .true. + CALL print_ks_energies ( ) + ! + ! ... save converged wfc if they have not been written previously + ! ... FIXME: it shouldn't be necessary to do this here + ! + IF ( nks == 1 .AND. (io_level < 2) .AND. (io_level > -1) ) & + CALL save_buffer ( evc, nwordwfc, iunwfc, nks ) + ! + ! ... do a Berry phase polarization calculation if required + ! + IF ( lberry ) CALL c_phase() + ! + ! ... do an orbital magnetization (Kubo terms) calculation + ! + IF ( lorbm ) CALL orbm_kubo() + ! + CALL stop_clock( 'electrons' ) + ! +9000 FORMAT(/' total cpu time spent up to now is ',F10.1,' secs' ) +9002 FORMAT(/' Band Structure Calculation' ) +9102 FORMAT(/' End of band structure calculation' ) + ! +END SUBROUTINE non_scf_ph + +