From 3ebf9c1831a4a56452ff215aef8df13559125348 Mon Sep 17 00:00:00 2001 From: andreaurru247 Date: Wed, 29 Jan 2020 11:37:15 +0100 Subject: [PATCH] Developments on the PH package routines to implement the calculation of phonon frequencies in the noncollinear magnetic case --- PHonon/PH/adddvscf_ph_mag.f90 | 127 +++++++ PHonon/PH/allocate_phq.f90 | 25 +- PHonon/PH/apply_trev.f90 | 56 +++ PHonon/PH/c_bands_ph.f90 | 179 +++++++++ PHonon/PH/deallocate_phq.f90 | 23 +- PHonon/PH/do_phonon.f90 | 8 + PHonon/PH/drhodv.f90 | 73 ++-- PHonon/PH/drhodvnl.f90 | 11 +- PHonon/PH/dvqpsi_us.f90 | 7 +- PHonon/PH/dvqpsi_us_only.f90 | 9 +- PHonon/PH/dynmatrix.f90 | 4 +- PHonon/PH/elphon.f90 | 7 +- PHonon/PH/ep_matrix_element_wannier.f90 | 5 +- PHonon/PH/incdrhous_nc.f90 | 2 +- PHonon/PH/initialize_ph.f90 | 57 ++- PHonon/PH/non_scf_ph.f90 | 112 ++++++ PHonon/PH/openfilq.f90 | 13 +- PHonon/PH/phcom.f90 | 28 +- PHonon/PH/phq_init.f90 | 90 +++-- PHonon/PH/phq_setup.f90 | 16 +- PHonon/PH/raman_mat.f90 | 4 +- PHonon/PH/run_nscf.f90 | 15 +- PHonon/PH/set_int12_nc.f90 | 26 +- PHonon/PH/set_irr.f90 | 4 +- PHonon/PH/set_irr_sym.f90 | 13 +- PHonon/PH/solve_linter.f90 | 458 +++++++++++++++--------- PHonon/PH/sym_dmag.f90 | 12 +- PHonon/PH/symdvscf.f90 | 35 +- PHonon/PH/symdynph_gq.f90 | 23 +- PHonon/PH/zstar_eu.f90 | 4 +- 30 files changed, 1139 insertions(+), 307 deletions(-) create mode 100644 PHonon/PH/adddvscf_ph_mag.f90 create mode 100644 PHonon/PH/apply_trev.f90 create mode 100644 PHonon/PH/c_bands_ph.f90 create mode 100644 PHonon/PH/non_scf_ph.f90 diff --git a/PHonon/PH/adddvscf_ph_mag.f90 b/PHonon/PH/adddvscf_ph_mag.f90 new file mode 100644 index 000000000..e42b98184 --- /dev/null +++ b/PHonon/PH/adddvscf_ph_mag.f90 @@ -0,0 +1,127 @@ +! +! 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 adddvscf_ph_mag (ipert, ik, becp1) + !---------------------------------------------------------------------- + ! + ! This routine computes the contribution of the selfconsistent + ! change of the potential to the known part of the linear + ! system and adds it to dvpsi. + ! It implements the second term in Eq. B30 of PRB 64, 235118 (2001). + ! + USE kinds, ONLY : DP + USE uspp_param, ONLY : upf, nh + USE uspp, ONLY : vkb, okvan +! modules from pwcom + USE lsda_mod, ONLY : lsda, current_spin, isk + USE ions_base, ONLY : ntyp => nsp, nat, ityp + USE wvfct, ONLY : nbnd, npwx + USE klist, ONLY : ngk + USE noncollin_module, ONLY : noncolin, npol + USE becmod, ONLY : bec_type +! modules from phcom + USE lrus, ONLY : int3, int3_nc + USE qpoint, ONLY : nksq, ikks, ikqs + USE eqv, ONLY : dvpsi + + implicit none + ! + ! The dummy variables + ! + integer :: ik, ipert + ! input: the k point + ! input: the perturbation + ! + TYPE(bec_type) :: becp1(nksq) + ! And the local variables + ! + integer :: na, nt, ibnd, ih, jh, ijkb0, ikb, jkb, is, js, ijs + ! counter on atoms + ! counter on atomic types + ! counter on bands + ! counter on beta functions + ! counter on beta functions + ! auxiliary variable for indexing + ! counter on the k points + ! counter on vkb + ! counter on vkb + integer :: ikk, ikq, npwq + ! index of the point k + ! index of the point k+q + ! number of the plane-waves at point k+q + complex(DP) :: sum, sum_nc(npol) + ! auxiliary variable + + if (.not.okvan) return + ! + call start_clock ('adddvscf') + ! + ikk = ikks(ik) + ikq = ikqs(ik) + npwq = ngk(ikq) + ! + if (lsda) current_spin = isk(ikk) + ! + ijkb0 = 0 + do nt = 1, ntyp + if (upf(nt)%tvanp ) then + do na = 1, nat + if (ityp (na) .eq.nt) then + ! + ! we multiply the integral for the becp term and the beta_n + ! + do ibnd = 1, nbnd + do ih = 1, nh (nt) + ikb = ijkb0 + ih + IF (noncolin) THEN + sum_nc = (0.d0, 0.d0) + ELSE + sum = (0.d0, 0.d0) + END IF + do jh = 1, nh (nt) + jkb = ijkb0 + jh + IF (noncolin) THEN + ijs=0 + do is=1,npol + do js=1,npol + ijs=ijs+1 + sum_nc(is)=sum_nc(is)+ & + int3_nc(ih,jh,na,ijs,ipert)* & + becp1(ik)%nc(jkb, js, ibnd) + enddo + enddo + ELSE + sum = sum + int3 (ih, jh, na, current_spin, ipert)*& + becp1(ik)%k(jkb, ibnd) + END IF + enddo + IF (noncolin) THEN + call zaxpy(npwq,sum_nc(1),vkb(1,ikb),1,dvpsi(1,ibnd),1) + call zaxpy(npwq,sum_nc(2),vkb(1,ikb),1, & + dvpsi(1+npwx,ibnd),1) + ELSE + call zaxpy(npwq,sum,vkb(1,ikb),1,dvpsi(1,ibnd),1) + END IF + enddo + enddo + ijkb0 = ijkb0 + nh (nt) + endif + enddo + else + do na = 1, nat + if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt) + enddo + endif + enddo + ! + call stop_clock ('adddvscf') + ! + return + ! +end subroutine adddvscf_ph_mag diff --git a/PHonon/PH/allocate_phq.f90 b/PHonon/PH/allocate_phq.f90 index 6ba8da1f6..a3537a215 100644 --- a/PHonon/PH/allocate_phq.f90 +++ b/PHonon/PH/allocate_phq.f90 @@ -22,7 +22,8 @@ subroutine allocate_phq USE noncollin_module, ONLY : noncolin, npol, nspin_mag USE fft_base, ONLY : dfftp USE wavefunctions, ONLY : evc - USE spin_orb, ONLY : lspinorb + USE spin_orb, ONLY : lspinorb, domag + USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save USE becmod, ONLY : bec_type, becp, allocate_bec_type USE uspp, ONLY : okvan, nkb, vkb USE paw_variables, ONLY : okpaw @@ -47,6 +48,7 @@ subroutine allocate_phq dwfcatomk, sdwfcatomk, wfcatomkpq, dwfcatomkpq, & swfcatomk, swfcatomkpq, sdwfcatomkpq, dvkb, vkbkpq, & dvkbkpq + USE qpoint_aux, ONLY : becpt, alphapt IMPLICIT NONE INTEGER :: ik, ipol, ldim @@ -93,6 +95,25 @@ subroutine allocate_phq zstareu0=(0.0_DP,0.0_DP) zstarue0=(0.0_DP,0.0_DP) zstarue0_rec=(0.0_DP,0.0_DP) + IF (noncolin.AND.domag) THEN + ALLOCATE (becpt(nksq)) + ALLOCATE (alphapt(3,nksq)) + DO ik=1,nksq + CALL allocate_bec_type ( nkb, nbnd, becpt(ik) ) + DO ipol=1,3 + CALL allocate_bec_type ( nkb, nbnd, alphapt(ipol,ik) ) + ENDDO + ENDDO + IF (okvan) THEN + ALLOCATE(int1_nc_save( nhm, nhm, 3, nat, nspin, 2)) + ALLOCATE (deeq_nc_save( nhm, nhm, nat, nspin, 2)) + ENDIF + ! AAA Nota: da definire this_pcxpsi_is_on_file_tpw + !ALLOCATE (this_pcxpsi_is_on_file_tpw(nksq,3,2)) + ! ELSE + !ALLOCATE (this_pcxpsi_is_on_file_tpw(nksq,3,1)) + ENDIF + if (okvan) then allocate (int1 ( nhm, nhm, 3, nat, nspin_mag)) allocate (int2 ( nhm , nhm , 3 , nat , nat)) @@ -113,6 +134,8 @@ subroutine allocate_phq allocate(dpqq_so( nhm, nhm, nspin, 3, ntyp)) END IF END IF + + allocate (alphasum ( nhm * (nhm + 1)/2 , 3 , nat , nspin_mag)) allocate (this_dvkb3_is_on_file(nksq)) this_dvkb3_is_on_file(:)=.false. diff --git a/PHonon/PH/apply_trev.f90 b/PHonon/PH/apply_trev.f90 new file mode 100644 index 000000000..fdc099f9a --- /dev/null +++ b/PHonon/PH/apply_trev.f90 @@ -0,0 +1,56 @@ +! +! Copyright (C) 2018 Andrea Dal Corso +! 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 +! evc at the k point ikk_evc and puts the output in evc with the order +! of G vectors of 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/PHonon/PH/c_bands_ph.f90 b/PHonon/PH/c_bands_ph.f90 new file mode 100644 index 000000000..5f81af40a --- /dev/null +++ b/PHonon/PH/c_bands_ph.f90 @@ -0,0 +1,179 @@ +! +! 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, U_projection, wfcU + 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 + USE spin_orb, ONLY : domag + USE save_ph, ONLY : tmp_dir_save + USE io_files, ONLY : tmp_dir, prefix + ! + IMPLICIT NONE + ! + REAL(DP) :: avg_iter, ethr_ + ! average number of H*psi products + INTEGER :: ik_, ik, nkdum, ios, iuawfc, lrawfc + ! 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 + IF (tmp_dir /= tmp_dir_save) THEN + iuawfc = 20 + lrawfc = nbnd * npwx * npol + CALL open_buffer (iuawfc, 'wfc', lrawfc, io_level, exst_mem, exst, & + tmp_dir_save) + IF (.NOT.exst.AND..NOT.exst_mem) THEN + CALL errore ('c_bands_ph', 'file '//trim(prefix)//'.wfc not found', 1) + END IF + ENDIF + ! + ! ... 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 + 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. (U_projection .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 + IF (tmp_dir /= tmp_dir_save) CALL close_buffer(iuawfc,'keep') + ! + 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/PHonon/PH/deallocate_phq.f90 b/PHonon/PH/deallocate_phq.f90 index bff4a9a27..6a72246b1 100644 --- a/PHonon/PH/deallocate_phq.f90 +++ b/PHonon/PH/deallocate_phq.f90 @@ -42,6 +42,9 @@ subroutine deallocate_phq wfcatomk, swfcatomk, dwfcatomk, sdwfcatomk, & wfcatomkpq, dwfcatomkpq, swfcatomkpq, sdwfcatomkpq, & dvkb, vkbkpq, dvkbkpq + USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt, alphapt + USE becmod, ONLY : deallocate_bec_type + USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save IMPLICIT NONE INTEGER :: ik, ipol @@ -61,6 +64,8 @@ subroutine deallocate_phq ! if(allocated(ikks)) deallocate (ikks) if(allocated(ikqs)) deallocate (ikqs) + IF(ALLOCATED(ikmks)) DEALLOCATE(ikmks) + IF(ALLOCATED(ikmkmqs)) DEALLOCATE(ikmkmqs) if(allocated(eigqts)) deallocate (eigqts) if(allocated(rtau)) deallocate (rtau) if(associated(u)) deallocate (u) @@ -98,11 +103,13 @@ subroutine deallocate_phq if(allocated(int2_so)) deallocate(int2_so) if(allocated(int5_so)) deallocate(int5_so) if(allocated(dpqq_so)) deallocate(dpqq_so) - + IF (ALLOCATED(int1_nc_save)) DEALLOCATE (int1_nc_save) + IF (ALLOCATED(deeq_nc_save)) DEALLOCATE (deeq_nc_save) if(allocated(alphasum)) deallocate (alphasum) if(allocated(this_dvkb3_is_on_file)) deallocate (this_dvkb3_is_on_file) if(allocated(this_pcxpsi_is_on_file)) deallocate (this_pcxpsi_is_on_file) + !IF (ALLOCATED(this_pcxpsi_is_on_file_tpw)) DEALLOCATE(this_pcxpsi_is_on_file_tpw) !!! AAA Da definire chi รจ this_pcxpsi_is_on_file_tpw if(allocated(alphap)) then do ik=1,nksq do ipol=1,3 @@ -117,6 +124,20 @@ subroutine deallocate_phq end do deallocate(becp1) end if + IF (ALLOCATED(alphapt)) THEN + DO ik=1,nksq + DO ipol=1,3 + CALL deallocate_bec_type ( alphapt(ipol,ik) ) + ENDDO + ENDDO + DEALLOCATE (alphapt) + ENDIF + IF (ALLOCATED(becpt)) THEN + DO ik=1, nksq + CALL deallocate_bec_type ( becpt(ik) ) + ENDDO + DEALLOCATE(becpt) + ENDIF call deallocate_bec_type ( becp ) if(allocated(el_ph_mat)) deallocate (el_ph_mat) diff --git a/PHonon/PH/do_phonon.f90 b/PHonon/PH/do_phonon.f90 index f1c86a103..8b05e94c2 100644 --- a/PHonon/PH/do_phonon.f90 +++ b/PHonon/PH/do_phonon.f90 @@ -41,6 +41,7 @@ SUBROUTINE do_phonon(auxdyn) ! USE elph_tetra_mod, ONLY : elph_tetra, elph_tetra_lambda, elph_tetra_gamma USE elph_scdft_mod, ONLY : elph_scdft + USE io_global, ONLY : stdout IMPLICIT NONE ! @@ -58,6 +59,13 @@ SUBROUTINE do_phonon(auxdyn) ! ! If necessary the bands are recalculated ! + ! Note (A. Urru): This has still to be cleaned (setup_pw + ! should be correctly set by prepare_q: here we force it + ! to be .true. in order for the code to work properly in + ! the case SO-MAG). + ! + setup_pw=.true. + WRITE(stdout,*) 'setup_pw', setup_pw IF (setup_pw) CALL run_nscf(do_band, iq) ! ! If only_wfc=.TRUE. the code computes only the wavefunctions diff --git a/PHonon/PH/drhodv.f90 b/PHonon/PH/drhodv.f90 index 61709f6eb..89b266ead 100644 --- a/PHonon/PH/drhodv.f90 +++ b/PHonon/PH/drhodv.f90 @@ -25,7 +25,7 @@ subroutine drhodv (nu_i0, nper, drhoscf) USE cell_base, ONLY : tpiba USE lsda_mod, ONLY : current_spin, lsda, isk, nspin USE wvfct, ONLY : npwx, nbnd - USE uspp, ONLY : nkb, vkb + USE uspp, ONLY : nkb, vkb, deeq_nc, okvan USE becmod, ONLY : calbec, bec_type, becscal, allocate_bec_type, & deallocate_bec_type USE fft_base, ONLY : dfftp @@ -40,6 +40,11 @@ subroutine drhodv (nu_i0, nper, drhoscf) USE eqv, ONLY : dpsi USE qpoint, ONLY : nksq, ikks, ikqs USE control_lr, ONLY : lgamma + USE spin_orb, ONLY : domag + USE lrus, ONLY : becp1 + USE phus, ONLY : alphap, int1_nc + USE qpoint_aux, ONLY : becpt, alphapt + USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save USE mp_pools, ONLY : inter_pool_comm USE mp, ONLY : mp_sum @@ -54,7 +59,7 @@ subroutine drhodv (nu_i0, nper, drhoscf) ! the change of density due to perturbations integer :: mu, ik, ikq, ig, nu_i, nu_j, na_jcart, ibnd, nrec, & - ipol, ikk, ipert, npw, npwq + ipol, ikk, ipert, npw, npwq, isolv, nsolv ! counters ! ikk: record position for wfc at k @@ -83,6 +88,8 @@ subroutine drhodv (nu_i0, nper, drhoscf) ! ! We need a sum over all k points ... ! + nsolv=1 + IF (noncolin.AND.domag) nsolv=2 do ik = 1, nksq ikk = ikks(ik) ikq = ikqs(ik) @@ -90,34 +97,51 @@ subroutine drhodv (nu_i0, nper, drhoscf) npwq= ngk(ikq) if (lsda) current_spin = isk (ikk) call init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) - do mu = 1, nper - nrec = (mu - 1) * nksq + ik - if (nksq > 1 .or. nper > 1) call get_buffer(dpsi, lrdwf, iudwf, nrec) - call calbec (npwq, vkb, dpsi, dbecq(mu) ) - do ipol = 1, 3 - aux=(0.d0,0.d0) - do ibnd = 1, nbnd - do ig = 1, npwq - aux (ig, ibnd) = dpsi (ig, ibnd) * & - (xk (ipol, ikq) + g (ipol, igk_k(ig,ikq) ) ) - enddo - if (noncolin) then + DO isolv=1, nsolv + do mu = 1, nper + nrec = (mu - 1) * nksq + ik + (isolv-1) * nksq * nper + if (nksq > 1 .or. nper > 1 .OR. nsolv==2) & + call get_buffer(dpsi, lrdwf, iudwf, nrec) + call calbec (npwq, vkb, dpsi, dbecq(mu) ) + do ipol = 1, 3 + aux=(0.d0,0.d0) + do ibnd = 1, nbnd do ig = 1, npwq - aux (ig+npwx, ibnd) = dpsi (ig+npwx, ibnd) * & - (xk (ipol, ikq) + g (ipol, igk_k(ig,ikq) ) ) + aux (ig, ibnd) = dpsi (ig, ibnd) * & + (xk (ipol, ikq) + g (ipol, igk_k(ig,ikq) ) ) enddo - endif + if (noncolin) then + do ig = 1, npwq + aux (ig+npwx, ibnd) = dpsi (ig+npwx, ibnd) * & + (xk (ipol, ikq) + g (ipol, igk_k(ig,ikq) ) ) + enddo + endif + enddo + call calbec (npwq, vkb, aux, dalpq(ipol,mu) ) enddo - call calbec (npwq, vkb, aux, dalpq(ipol,mu) ) enddo - enddo - fact = CMPLX(0.d0, tpiba,kind=DP) - DO ipert=1,nper - DO ipol=1,3 - CALL becscal(fact,dalpq(ipol,ipert),nkb,nbnd) + fact = CMPLX(0.d0, tpiba,kind=DP) + DO ipert=1,nper + DO ipol=1,3 + CALL becscal(fact,dalpq(ipol,ipert),nkb,nbnd) + ENDDO ENDDO + IF (isolv==1) THEN + call drhodvnl (ik, ikk, nper, nu_i0, dynwrk, becp1, alphap, & + dbecq, dalpq) + ELSE + IF (okvan) THEN + deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2) + int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,2) + ENDIF + call drhodvnl (ik, ikk, nper, nu_i0, dynwrk, becpt, alphapt, & + dbecq, dalpq) + IF (okvan) THEN + deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) + int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,1) + ENDIF + ENDIF ENDDO - call drhodvnl (ik, ikk, nper, nu_i0, dynwrk, dbecq, dalpq) enddo ! ! put in the basis of the modes @@ -136,6 +160,7 @@ subroutine drhodv (nu_i0, nper, drhoscf) ! collect contributions from all pools (sum over k-points) ! call mp_sum ( wdyn, inter_pool_comm ) + IF (nsolv==2) wdyn=wdyn/2.0_DP ! ! add the contribution of the local part of the perturbation ! diff --git a/PHonon/PH/drhodvnl.f90 b/PHonon/PH/drhodvnl.f90 index ad2bd75e2..b9bec0657 100644 --- a/PHonon/PH/drhodvnl.f90 +++ b/PHonon/PH/drhodvnl.f90 @@ -7,7 +7,8 @@ ! ! !----------------------------------------------------------------------- -subroutine drhodvnl (ik, ikk, nper, nu_i0, wdyn, dbecq, dalpq) +subroutine drhodvnl (ik, ikk, nper, nu_i0, wdyn, becp1, alphap, & + dbecq, dalpq) !----------------------------------------------------------------------- ! ! This subroutine computes the electronic term 2 of @@ -27,10 +28,8 @@ subroutine drhodvnl (ik, ikk, nper, nu_i0, wdyn, dbecq, dalpq) USE klist, ONLY : wk USE lsda_mod, ONLY : current_spin, nspin USE spin_orb, ONLY : lspinorb - USE phus, ONLY : int1, int1_nc, int2, int2_so, alphap - - USE lrus, ONLY : becp1 - + USE phus, ONLY : int1, int1_nc, int2, int2_so + USE qpoint, ONLY : nksq USE mp_bands, ONLY: intra_bgrp_comm USE mp, ONLY: mp_sum @@ -40,7 +39,7 @@ subroutine drhodvnl (ik, ikk, nper, nu_i0, wdyn, dbecq, dalpq) ! input: the number of perturbations ! input: the initial mode - TYPE(bec_type) :: dbecq(nper), dalpq(3,nper) + TYPE(bec_type) :: dbecq(nper), dalpq(3,nper), becp1(nksq), alphap(3,nksq) ! input: the becp with psi_{k+q} ! input: the alphap with psi_{k} complex(DP) :: wdyn (3 * nat, 3 * nat) diff --git a/PHonon/PH/dvqpsi_us.f90 b/PHonon/PH/dvqpsi_us.f90 index bf9f44b7f..72b8a1e5c 100644 --- a/PHonon/PH/dvqpsi_us.f90 +++ b/PHonon/PH/dvqpsi_us.f90 @@ -7,7 +7,7 @@ ! ! !---------------------------------------------------------------------- -subroutine dvqpsi_us (ik, uact, addnlcc) +subroutine dvqpsi_us (ik, uact, addnlcc, becp1, alphap) !---------------------------------------------------------------------- ! ! This routine calculates dV_bare/dtau * psi for one perturbation @@ -43,6 +43,8 @@ subroutine dvqpsi_us (ik, uact, addnlcc) USE Coul_cut_2D, ONLY: do_cutoff_2D USE Coul_cut_2D_ph, ONLY : cutoff_localq + USE qpoint, ONLY : nksq + USE becmod, ONLY : bec_type ! IMPLICIT NONE ! @@ -83,6 +85,7 @@ subroutine dvqpsi_us (ik, uact, addnlcc) !! INTEGER :: ip !! + TYPE(bec_type) :: becp1(nksq), alphap(3,nksq) ! complex(DP) :: gtau, gu, fact, u1, u2, u3, gu0 complex(DP) , allocatable :: aux (:,:) @@ -243,7 +246,7 @@ subroutine dvqpsi_us (ik, uact, addnlcc) ! First a term similar to the KB case. ! Then a term due to the change of the D coefficients. ! - call dvqpsi_us_only (ik, uact) + call dvqpsi_us_only (ik, uact, becp1, alphap) call stop_clock ('dvqpsi_us') return diff --git a/PHonon/PH/dvqpsi_us_only.f90 b/PHonon/PH/dvqpsi_us_only.f90 index 10081215e..04a61ff98 100644 --- a/PHonon/PH/dvqpsi_us_only.f90 +++ b/PHonon/PH/dvqpsi_us_only.f90 @@ -7,7 +7,7 @@ ! ! !---------------------------------------------------------------------- -subroutine dvqpsi_us_only (ik, uact) +subroutine dvqpsi_us_only (ik, uact, becp1, alphap) !---------------------------------------------------------------------- ! ! This routine calculates dV_bare/dtau * psi for one perturbation @@ -29,10 +29,10 @@ subroutine dvqpsi_us_only (ik, uact) USE noncollin_module, ONLY : noncolin, npol USE uspp, ONLY: okvan, nkb, vkb USE uspp_param, ONLY: nh, nhm - USE phus, ONLY : int1, int1_nc, int2, int2_so, alphap + USE phus, ONLY : int1, int1_nc, int2, int2_so - USE lrus, ONLY : becp1 - USE qpoint, ONLY : ikks, ikqs + USE qpoint, ONLY : nksq, ikks, ikqs + USE becmod, ONLY : bec_type USE eqv, ONLY : dvpsi USE control_lr, ONLY : lgamma @@ -45,6 +45,7 @@ subroutine dvqpsi_us_only (ik, uact) ! input: the k point complex(DP) :: uact (3 * nat) ! input: the pattern of displacements + TYPE(bec_type) :: becp1(nksq), alphap(3,nksq) ! ! And the local variables ! diff --git a/PHonon/PH/dynmatrix.f90 b/PHonon/PH/dynmatrix.f90 index e93bff303..d618c342f 100644 --- a/PHonon/PH/dynmatrix.f90 +++ b/PHonon/PH/dynmatrix.f90 @@ -21,7 +21,7 @@ subroutine dynmatrix_new(iq_) USE io_global, ONLY : stdout USE control_flags, ONLY : modenum USE cell_base, ONLY : at, bg, celldm, ibrav, omega - USE symm_base, ONLY : s, sr, irt, nsym, invs + USE symm_base, ONLY : s, sr, irt, nsym, invs, t_rev USE dynmat, ONLY : dyn, w2 USE noncollin_module, ONLY : nspin_mag USE modes, ONLY : u, nmodes, npert, nirr, num_rap_mode @@ -145,7 +145,7 @@ subroutine dynmatrix_new(iq_) ! ! Generates the star of q ! - call star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, .TRUE. ) + call star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, .TRUE., t_rev ) ! ! write on file information on the system ! diff --git a/PHonon/PH/elphon.f90 b/PHonon/PH/elphon.f90 index eed355d6c..4607b09bf 100644 --- a/PHonon/PH/elphon.f90 +++ b/PHonon/PH/elphon.f90 @@ -36,8 +36,7 @@ SUBROUTINE elphon() 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 lrus, ONLY : int3, int3_nc, int3_paw, becp1 USE qpoint, ONLY : xq ! IMPLICIT NONE @@ -310,6 +309,8 @@ SUBROUTINE elphel (irr, npe, imode0, dvscfins) 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 ! @@ -410,7 +411,7 @@ SUBROUTINE elphel (irr, npe, imode0, dvscfins) ELSE mode = imode0 + ipert ! FIXME: .false. or .true. ??? - CALL dvqpsi_us (ik, u (1, mode), .FALSE. ) + CALL dvqpsi_us (ik, u (1, mode), .FALSE., becp1, alphap) ! ! DFPT+U: calculate the bare derivative of the Hubbard potential in el-ph ! diff --git a/PHonon/PH/ep_matrix_element_wannier.f90 b/PHonon/PH/ep_matrix_element_wannier.f90 index b8da7428a..d388dd419 100644 --- a/PHonon/PH/ep_matrix_element_wannier.f90 +++ b/PHonon/PH/ep_matrix_element_wannier.f90 @@ -31,7 +31,6 @@ SUBROUTINE ep_matrix_element_wannier() USE paw_variables, ONLY : okpaw USE uspp_param, ONLY : nhm USE lsda_mod, ONLY : nspin - USE lrus, ONLY : int3, int3_nc, int3_paw USE qpoint, ONLY : xq, nksq, ikks ! @@ -351,6 +350,8 @@ SUBROUTINE elphel_refolded (npe, imode0, dvscfins) USE eqv, ONLY : dvpsi!, evq USE qpoint, ONLY : nksq, ikks, ikqs USE control_lr, ONLY : lgamma + USE lrus, ONLY : becp1 + USE phus, ONLY : alphap IMPLICIT NONE ! @@ -438,7 +439,7 @@ SUBROUTINE elphel_refolded (npe, imode0, dvscfins) ELSE mode = imode0 + ipert ! FIXME : .false. or .true. ??? - CALL dvqpsi_us (ik, u (1, mode), .FALSE. ) + CALL dvqpsi_us (ik, u (1, mode), .FALSE., becp1, alphap) ENDIF ! ! calculate dvscf_q*psi_k diff --git a/PHonon/PH/incdrhous_nc.f90 b/PHonon/PH/incdrhous_nc.f90 index 785212177..138178d59 100644 --- a/PHonon/PH/incdrhous_nc.f90 +++ b/PHonon/PH/incdrhous_nc.f90 @@ -168,7 +168,7 @@ subroutine incdrhous_nc (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, & enddo enddo - call addusdbec_nc (ik, weight, dpsi, dbecsum) + call addusdbec_nc (ik, weight, dpsi, dbecsum, becp1) deallocate (ps1) deallocate (dpsir) diff --git a/PHonon/PH/initialize_ph.f90 b/PHonon/PH/initialize_ph.f90 index 9eee2d7b6..4b6d8ae45 100644 --- a/PHonon/PH/initialize_ph.f90 +++ b/PHonon/PH/initialize_ph.f90 @@ -14,7 +14,10 @@ SUBROUTINE initialize_ph() USE klist, ONLY : nks, nkstot ! USE qpoint, ONLY : nksq, nksqtot, ikks, ikqs + USE qpoint_aux, ONLY : ikmks, ikmkmqs USE control_lr, ONLY : lgamma + USE noncollin_module, ONLY : noncolin + USE spin_orb, ONLY : domag ! IMPLICIT NONE INTEGER :: ik @@ -23,23 +26,49 @@ SUBROUTINE initialize_ph() ! IF ( lgamma ) THEN ! - nksq = nks - nksqtot = nkstot - ALLOCATE(ikks(nksq), ikqs(nksq)) - DO ik=1,nksq - ikks(ik) = ik - ikqs(ik) = ik - ENDDO + IF (noncolin.AND.domag) THEN + nksq = nks/2 + nksqtot = nkstot/2 + ALLOCATE(ikks(nksq), ikqs(nksq)) + ALLOCATE(ikmks(nksq), ikmkmqs(nksq)) + DO ik=1,nksq + ikks(ik) = 2*ik-1 + ikqs(ik) = 2*ik-1 + ikmks(ik) = 2*ik + ikmkmqs(ik) = 2*ik + ENDDO + ELSE + nksq = nks + nksqtot = nkstot + ALLOCATE(ikks(nksq), ikqs(nksq)) + DO ik=1,nksq + ikks(ik) = ik + ikqs(ik) = ik + ENDDO + END IF ! ELSE ! - nksq = nks / 2 - nksqtot = nkstot / 2 - ALLOCATE(ikks(nksq), ikqs(nksq)) - DO ik=1,nksq - ikks(ik) = 2 * ik - 1 - ikqs(ik) = 2 * ik - ENDDO + IF (noncolin.AND.domag) THEN + nksq = nks / 4 + nksqtot = nkstot / 4 + ALLOCATE(ikks(nksq), ikqs(nksq)) + ALLOCATE(ikmks(nksq), ikmkmqs(nksq)) + DO ik=1,nksq + ikks(ik) = 4 * ik - 3 + ikqs(ik) = 4 * ik - 2 + ikmks(ik) = 4 * ik - 1 + ikmkmqs(ik) = 4 * ik + ENDDO + ELSE + nksq = nks / 2 + nksqtot = nkstot / 2 + ALLOCATE(ikks(nksq), ikqs(nksq)) + DO ik=1,nksq + ikks(ik) = 2 * ik - 1 + ikqs(ik) = 2 * ik + ENDDO + END IF ! END IF ! diff --git a/PHonon/PH/non_scf_ph.f90 b/PHonon/PH/non_scf_ph.f90 new file mode 100644 index 000000000..640950a21 --- /dev/null +++ b/PHonon/PH/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 + diff --git a/PHonon/PH/openfilq.f90 b/PHonon/PH/openfilq.f90 index f3848e6ce..4a2abfd35 100644 --- a/PHonon/PH/openfilq.f90 +++ b/PHonon/PH/openfilq.f90 @@ -12,16 +12,16 @@ SUBROUTINE openfilq() ! ... This subroutine opens all the files necessary for the phononq ! ... calculation. ! - USE kinds, ONLY : DP - USE control_flags, ONLY : io_level, modenum - USE units_ph, ONLY : iudwf, iubar, iucom, iudvkb3, & + USE kinds, ONLY : DP + USE control_flags, ONLY : io_level, modenum + USE units_ph, ONLY : iudwf, iubar, iucom, iudvkb3, & iudrhous, iuebar, iudrho, iudyn, iudvscf, & lrdwf, lrbar, lrcom, lrdvkb3, & lrdrhous, lrebar, lrdrho, lint3paw, iuint3paw, & iundnsscf - USE units_lr, ONLY : iuwfc, lrwfc - USE io_files, ONLY : tmp_dir, diropn, seqopn, nwordwfcU - USE control_ph, ONLY : epsil, zue, ext_recover, trans, & + USE units_lr, ONLY : iuwfc, lrwfc + USE io_files, ONLY : tmp_dir, diropn, seqopn, nwordwfcU + USE control_ph, ONLY : epsil, zue, ext_recover, trans, & tmp_dir_phq, start_irr, last_irr, xmldyn, & all_done, newgrid USE save_ph, ONLY : tmp_dir_save @@ -82,6 +82,7 @@ SUBROUTINE openfilq() ELSE ! this is the standard treatment IF (lgamma.AND.modenum==0.AND..NOT.newgrid ) tmp_dir=tmp_dir_save + IF ((noncolin.AND.domag).OR.lsda) tmp_dir=tmp_dir_phq ENDIF !!!!!!!!!!!!!!!!!!!!!!!! END OF ACFDT TEST !!!!!!!!!!!!!!!! iuwfc = 20 diff --git a/PHonon/PH/phcom.f90 b/PHonon/PH/phcom.f90 index ea6b4a20e..5c41a75c0 100644 --- a/PHonon/PH/phcom.f90 +++ b/PHonon/PH/phcom.f90 @@ -230,8 +230,9 @@ MODULE control_ph only_init=.FALSE., &! if .TRUE. computes only initial stuff with_ext_images=.FALSE., & ! if .TRUE. use an external driver ! to decide what each image does. - always_run=.FALSE., & ! if .TRUE. the code do not stop after + !always_run=.FALSE., & ! if .TRUE. the code do not stop after ! doing partial representations + always_run=.TRUE., & ! only for testing purposes recover, &! if .TRUE. the run restarts low_directory_check=.FALSE., & ! if .TRUE. search on the phsave ! directory only the representations requested @@ -426,6 +427,29 @@ MODULE ldaU_ph ! END MODULE ldaU_ph +MODULE nc_mag_aux + USE kinds, ONLY : DP + SAVE + + COMPLEX (DP), ALLOCATABLE :: & + deeq_nc_save(:,:,:,:,:), & + int1_nc_save(:,:,:,:,:,:), & + int3_save(:, :, :, :, :, :) +END MODULE nc_mag_aux + +!MODULE qpoint_aux +! USE kinds, ONLY : DP +! USE becmod, ONLY : bec_type +! SAVE + +! INTEGER, ALLOCATABLE :: ikmks(:) ! index of -k for magnetic calculations + +! INTEGER, ALLOCATABLE :: ikmkmqs(:) ! index of -k-q for magnetic calculations + +! TYPE(bec_type), ALLOCATABLE :: becpt(:), alphapt(:,:) + +!END MODULE qpoint_aux + MODULE phcom USE modes USE dynmat @@ -442,4 +466,6 @@ MODULE phcom USE disp USE grid_irr_iq USE ldaU_ph + USE nc_mag_aux +! USE qpoint_aux END MODULE phcom diff --git a/PHonon/PH/phq_init.f90 b/PHonon/PH/phq_init.f90 index a0bcfa928..6dd535558 100644 --- a/PHonon/PH/phq_init.f90 +++ b/PHonon/PH/phq_init.f90 @@ -44,7 +44,7 @@ SUBROUTINE phq_init() USE io_global, ONLY : stdout USE atom, ONLY : msh, rgrid USE vlocal, ONLY : strf - USE spin_orb, ONLY : lspinorb + USE spin_orb, ONLY : lspinorb, domag USE wvfct, ONLY : npwx, nbnd USE gvecw, ONLY : gcutw USE wavefunctions, ONLY : evc @@ -66,6 +66,7 @@ SUBROUTINE phq_init() USE Coul_cut_2D_ph, ONLY : cutoff_lr_Vlocq , cutoff_fact_qg USE lrus, ONLY : becp1, dpqq, dpqq_so USE qpoint, ONLY : xq, nksq, eigqts, ikks, ikqs + USE qpoint_aux, ONLY : becpt, alphapt, ikmks USE eqv, ONLY : vlocq, evq USE control_lr, ONLY : nbnd_occ, lgamma USE ldaU, ONLY : lda_plus_u @@ -87,7 +88,7 @@ SUBROUTINE phq_init() INTEGER :: npw, npwq REAL(DP) :: arg ! the argument of the phase - COMPLEX(DP), ALLOCATABLE :: aux1(:,:) + COMPLEX(DP), ALLOCATABLE :: aux1(:,:), tevc(:,:) ! used to compute alphap ! ! @@ -157,6 +158,7 @@ SUBROUTINE phq_init() endif ! ALLOCATE( aux1( npwx*npol, nbnd ) ) + IF (noncolin.AND.domag) ALLOCATE(tevc(npwx*npol,nbnd)) ! DO ik = 1, nksq ! @@ -188,19 +190,23 @@ SUBROUTINE phq_init() ! ! ... read the wavefunctions at k ! - if(elph_mat) then + if(elph_mat) then call read_wfc_rspace_and_fwfft( evc, ik, lrwfcr, iunwfcwann, npw, igk_k(1,ikk) ) -! CALL davcio (evc, lrwfc, iunwfcwann, ik, - 1) - else - CALL get_buffer( evc, lrwfc, iuwfc, ikk ) - endif + ! CALL davcio (evc, lrwfc, iunwfcwann, ik, - 1) + else + CALL get_buffer( evc, lrwfc, iuwfc, ikk ) + IF (noncolin.AND.domag) THEN + CALL get_buffer( tevc, lrwfc, iuwfc, ikmks(ik) ) + CALL calbec (npw, vkb, tevc, becpt(ik) ) + ENDIF + endif ! ! ... e) we compute the becp terms which are used in the rest of ! ... the code ! - + CALL calbec (npw, vkb, evc, becp1(ik) ) - + ! ! ... e') we compute the derivative of the becp term with respect to an ! atomic displacement @@ -210,43 +216,61 @@ SUBROUTINE phq_init() DO ibnd = 1, nbnd DO ig = 1, npw aux1(ig,ibnd) = evc(ig,ibnd) * tpiba * ( 0.D0, 1.D0 ) * & - ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) + ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) END DO IF (noncolin) THEN DO ig = 1, npw aux1(ig+npwx,ibnd)=evc(ig+npwx,ibnd)*tpiba*(0.D0,1.D0)*& - ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) + ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) END DO END IF END DO CALL calbec (npw, vkb, aux1, alphap(ipol,ik) ) END DO ! + IF (noncolin.AND.domag) THEN + DO ipol = 1, 3 + aux1=(0.d0,0.d0) + DO ibnd = 1, nbnd + DO ig = 1, npw + aux1(ig,ibnd) = tevc(ig,ibnd) * tpiba * ( 0.D0, 1.D0 ) * & + ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) + END DO + IF (noncolin) THEN + DO ig = 1, npw + aux1(ig+npwx,ibnd)=tevc(ig+npwx,ibnd)*tpiba*(0.D0,1.D0)*& + ( xk(ipol,ikk) + g(ipol,igk_k(ig,ikk)) ) + END DO + END IF + END DO + CALL calbec (npw, vkb, aux1, alphapt(ipol,ik) ) + END DO + ENDIF !!!!!!!!!!!!!!!!!!!!!!!! ACFDT TEST !!!!!!!!!!!!!!!! - IF (acfdt_is_active) THEN - ! ACFDT -test always read calculated wcf from non_scf calculation - IF(acfdt_num_der) then - CALL get_buffer( evq, lrwfc, iuwfc, ikq ) + IF (acfdt_is_active) THEN + ! ACFDT -test always read calculated wcf from non_scf calculation + IF(acfdt_num_der) then + CALL get_buffer( evq, lrwfc, iuwfc, ikq ) + ELSE + IF ( .NOT. lgamma ) & + CALL get_buffer( evq, lrwfc, iuwfc, ikq ) + ENDIF ELSE - IF ( .NOT. lgamma ) & - CALL get_buffer( evq, lrwfc, iuwfc, ikq ) + ! this is the standard treatment + IF ( .NOT. lgamma .and..not. elph_mat )then + CALL get_buffer( evq, lrwfc, iuwfc, ikq ) + ELSEIF(.NOT. lgamma .and. elph_mat) then + ! + ! I read the wavefunction in real space and fwfft it + ! + ikqg = kpq(ik) + call read_wfc_rspace_and_fwfft( evq, ikqg, lrwfcr, iunwfcwann, npwq, & + igk_k(1,ikq) ) + ! CALL davcio (evq, lrwfc, iunwfcwann, ikqg, - 1) + call calculate_and_apply_phase(ik, ikqg, igqg, & + npwq_refolded, g_kpq, xk_gamma, evq, .false.) + ENDIF ENDIF - ELSE - ! this is the standard treatment - IF ( .NOT. lgamma .and..not. elph_mat )then - CALL get_buffer( evq, lrwfc, iuwfc, ikq ) - ELSEIF(.NOT. lgamma .and. elph_mat) then - ! - ! I read the wavefunction in real space and fwfft it - ! - ikqg = kpq(ik) - call read_wfc_rspace_and_fwfft( evq, ikqg, lrwfcr, iunwfcwann, npwq, & - igk_k(1,ikq) ) -! CALL davcio (evq, lrwfc, iunwfcwann, ikqg, - 1) - call calculate_and_apply_phase(ik, ikqg, igqg, & - npwq_refolded, g_kpq, xk_gamma, evq, .false.) - ENDIF - ENDIF !!!!!!!!!!!!!!!!!!!!!!!! END OF ACFDT TEST !!!!!!!!!!!!!!!! ! END DO diff --git a/PHonon/PH/phq_setup.f90 b/PHonon/PH/phq_setup.f90 index abbb5ff4e..b3bf338b4 100644 --- a/PHonon/PH/phq_setup.f90 +++ b/PHonon/PH/phq_setup.f90 @@ -57,13 +57,14 @@ subroutine phq_setup USE klist, ONLY : xk, nks, nkstot USE lsda_mod, ONLY : nspin, starting_magnetization USE scf, ONLY : v, vrs, vltot, kedtau, rho + USE dfunct, ONLY : newd USE fft_base, ONLY : dfftp USE gvect, ONLY : ngm USE gvecs, ONLY : doublegrid USE symm_base, ONLY : nrot, nsym, s, irt, t_rev, time_reversal, & sr, invs, inverse_s, d1, d2, d3 USE uspp_param, ONLY : upf - USE uspp, ONLY : nlcc_any + USE uspp, ONLY : nlcc_any, deeq_nc, okvan USE spin_orb, ONLY : domag USE noncollin_module, ONLY : noncolin, m_loc, angle1, angle2, ux USE nlcc_ph, ONLY : drc @@ -95,6 +96,7 @@ subroutine phq_setup USE mp, ONLY : mp_max, mp_min USE lr_symm_base, ONLY : gi, gimq, irotmq, minus_q, invsymq, nsymq, rtau USE qpoint, ONLY : xq, xk_col + USE nc_mag_aux, ONLY : deeq_nc_save USE control_lr, ONLY : lgamma USE ldaU, ONLY : lda_plus_u, Hubbard_U, Hubbard_J0 USE ldaU_ph, ONLY : effU @@ -167,6 +169,18 @@ subroutine phq_setup END DO ux=0.0_DP if (dft_is_gradient()) call compute_ux(m_loc,ux,nat) + IF (okvan) THEN +! +! Change the sign of the magnetic field in the screened US coefficients +! and save also the coefficients computed with -B_xc. +! + deeq_nc_save(:,:,:,:,1)=deeq_nc(:,:,:,:) + v%of_r(:,2:4)=-v%of_r(:,2:4) + CALL newd() + v%of_r(:,2:4)=-v%of_r(:,2:4) + deeq_nc_save(:,:,:,:,2)=deeq_nc(:,:,:,:) + deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) + ENDIF ENDIF ! ! 3) Computes the derivative of the XC potential diff --git a/PHonon/PH/raman_mat.f90 b/PHonon/PH/raman_mat.f90 index fb7b92dbb..a07b40f95 100644 --- a/PHonon/PH/raman_mat.f90 +++ b/PHonon/PH/raman_mat.f90 @@ -155,7 +155,7 @@ subroutine raman_mat enddo do imod = 1, 3 * nat - call dvqpsi_us (ik, uact (1, imod),.false. ) + call dvqpsi_us (ik, uact (1, imod),.false., becp1, alphap) do ipa = 1, 6 tmp = 0.d0 do ibnd = 1, nbnd_occ (ik) @@ -198,7 +198,7 @@ subroutine raman_mat call calbec (npw, vkb, aux1, alphap (ipb,ik) ) enddo - call dvqpsi_us (ik, uact (1, imod),.false. ) + call dvqpsi_us (ik, uact (1, imod),.false., becp1, alphap ) do ipb = 1, ipa tmp = 0.d0 do ibnd = 1, nbnd_occ (ik) diff --git a/PHonon/PH/run_nscf.f90 b/PHonon/PH/run_nscf.f90 index c8da93d77..1d430d903 100644 --- a/PHonon/PH/run_nscf.f90 +++ b/PHonon/PH/run_nscf.f90 @@ -31,7 +31,7 @@ SUBROUTINE run_nscf(do_band, iq) ext_restart, bands_computed, newgrid, qplot, & only_wfc USE io_global, ONLY : stdout - USE save_ph, ONLY : tmp_dir_save + !USE save_ph, ONLY : tmp_dir_save ! USE grid_irr_iq, ONLY : done_bands USE acfdtest, ONLY : acfdt_is_active, acfdt_num_der, ir_point, delta_vrs @@ -41,6 +41,8 @@ SUBROUTINE run_nscf(do_band, iq) USE lr_symm_base, ONLY : minus_q, nsymq, invsymq USE control_lr, ONLY : ethr_nscf USE qpoint, ONLY : xq + USE noncollin_module,ONLY : noncolin + USE spin_orb, ONLY : domag USE klist, ONLY : qnorm, nelec USE el_phon, ONLY : elph_mat ! @@ -63,10 +65,11 @@ SUBROUTINE run_nscf(do_band, iq) ! FIXME: kunit is set here: in this case we do not go through setup_nscf ! FIXME: and read_file calls divide_et_impera that needs kunit ! FIXME: qnorm (also set in setup_nscf) is needed by allocate_nlpot - IF ( lgamma_iq(iq) ) THEN - kunit = 1 - ELSE - kunit = 2 + kunit = 2 + IF ( lgamma_iq(iq) ) kunit = 1 + IF (noncolin.AND.domag) THEN + kunit = 4 + IF (lgamma_iq(iq)) kunit=2 ENDIF qnorm = SQRT(xq(1)**2+xq(2)**2+xq(3)**2) ! @@ -109,7 +112,7 @@ SUBROUTINE run_nscf(do_band, iq) ENDIF !!!!!!!!!!!!!!!!!!!!!!!!END OF ACFDT TEST !!!!!!!!!!!!!!!! ! - IF (do_band) CALL non_scf ( ) + IF (do_band) CALL non_scf_ph ( ) IF ( check_stop_now() ) THEN diff --git a/PHonon/PH/set_int12_nc.f90 b/PHonon/PH/set_int12_nc.f90 index 5b6ca9808..e9775254b 100644 --- a/PHonon/PH/set_int12_nc.f90 +++ b/PHonon/PH/set_int12_nc.f90 @@ -13,13 +13,35 @@ SUBROUTINE set_int12_nc(iflag) ! by the Pauli matrices the integrals. ! USE ions_base, ONLY : nat, ntyp => nsp, ityp -USE spin_orb, ONLY : lspinorb +USE spin_orb, ONLY : lspinorb, domag +USE noncollin_module, ONLY : noncolin USE uspp_param, only: upf USE phus, ONLY : int1, int2, int1_nc, int2_so +USE nc_mag_aux, ONLY : int1_nc_save IMPLICIT NONE INTEGER :: iflag INTEGER :: np, na +IF (noncolin.AND.domag) THEN + int1_nc=(0.d0,0.d0) + int1 ( :, :, :, :, 2:4)=-int1 ( :, :, :, :, 2:4) + DO np = 1, ntyp + IF ( upf(np)%tvanp ) THEN + DO na = 1, nat + IF (ityp(na)==np) THEN + IF (upf(np)%has_so) THEN + CALL transform_int1_so(int1,na,iflag) + ELSE + CALL transform_int1_nc(int1,na,iflag) + END IF + END IF + END DO + END IF + END DO + int1 ( :, :, :, :, 2:4)=-int1 ( :, :, :, :, 2:4) + int1_nc_save(:,:,:,:,:,2) = int1_nc +END IF + int1_nc=(0.d0,0.d0) IF (lspinorb) int2_so=(0.d0,0.d0) DO np = 1, ntyp @@ -37,5 +59,7 @@ DO np = 1, ntyp END DO END IF END DO +IF (noncolin.AND.domag) int1_nc_save(:,:,:,:,:,1) = int1_nc + END SUBROUTINE set_int12_nc diff --git a/PHonon/PH/set_irr.f90 b/PHonon/PH/set_irr.f90 index 7ac7d156e..5991ac781 100644 --- a/PHonon/PH/set_irr.f90 +++ b/PHonon/PH/set_irr.f90 @@ -139,7 +139,7 @@ subroutine set_irr_new (xq, u, npert, nirr, eigen) end do end if - IF (search_sym) THEN + IF (search_sym.AND.nspin_mag/=4) THEN CALL find_mode_sym_new (u, eigen, tau, nat, nsymq, s, sr, irt, xq, & rtau, amass, ntyp, ityp, 0, .FALSE., .TRUE., num_rap_mode, ierr) @@ -246,7 +246,7 @@ subroutine set_irr_new (xq, u, npert, nirr, eigen) endif enddo - IF (search_sym) THEN + IF (search_sym.AND.nspin_mag/=4) THEN ! ! Here we set the name of the representation for each mode ! diff --git a/PHonon/PH/set_irr_sym.f90 b/PHonon/PH/set_irr_sym.f90 index 38d32e986..fb6624af7 100644 --- a/PHonon/PH/set_irr_sym.f90 +++ b/PHonon/PH/set_irr_sym.f90 @@ -18,7 +18,7 @@ subroutine set_irr_sym_new ( t, tmq, npertx ) USE constants, ONLY: tpi USE ions_base, ONLY : nat USE cell_base, ONLY : at, bg - USE symm_base, ONLY : s, irt + USE symm_base, ONLY : s, irt, t_rev USE modes, ONLY : u, nirr, npert USE control_flags, ONLY : modenum USE mp, ONLY : mp_bcast @@ -99,7 +99,7 @@ subroutine set_irr_sym_new ( t, tmq, npertx ) arg = arg + xq (ipol) * rtau (ipol, irot, na) enddo arg = arg * tpi - if (isymq.eq.nsymtot.and.minus_q) then + if ((isymq.eq.nsymtot.and.minus_q).OR.(t_rev(isymq)==1)) then fase = CMPLX (cos (arg), sin (arg), KIND=dp ) else fase = CMPLX (cos (arg),-sin (arg), KIND=dp ) @@ -130,8 +130,13 @@ subroutine set_irr_sym_new ( t, tmq, npertx ) tmq (jpert, ipert, irr) = tmq (jpert, ipert, irr) + CONJG(u ( & jmode, imode) * wrk_ru (ipol, na) ) else - t (jpert, ipert, irot, irr) = t (jpert, ipert, irot, irr) & - + CONJG(u (jmode, imode) ) * wrk_ru (ipol, na) + IF (t_rev(irot)==1) THEN + t (jpert,ipert,irot,irr)=t(jpert,ipert,irot,irr) & + + CONJG(u (jmode, imode) * wrk_ru (ipol, na)) + ELSE + t (jpert,ipert,irot,irr)=t(jpert,ipert,irot,irr) & + + CONJG(u (jmode, imode) ) * wrk_ru (ipol, na) + ENDIF endif enddo enddo diff --git a/PHonon/PH/solve_linter.f90 b/PHonon/PH/solve_linter.f90 index 620387bba..8cb9cb420 100644 --- a/PHonon/PH/solve_linter.f90 +++ b/PHonon/PH/solve_linter.f90 @@ -32,19 +32,17 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) USE io_files, ONLY : prefix, diropn USE check_stop, ONLY : check_stop_now USE wavefunctions, ONLY : evc - USE constants, ONLY : degspin - USE cell_base, ONLY : at, tpiba2 - USE klist, ONLY : ltetra, lgauss, degauss, ngauss, & - xk, wk, ngk, igk_k + USE cell_base, ONLY : at + USE klist, ONLY : ltetra, lgauss, xk, wk, ngk, igk_k USE gvect, ONLY : g USE gvecs, ONLY : doublegrid USE fft_base, ONLY : dfftp, dffts USE lsda_mod, ONLY : lsda, nspin, current_spin, isk USE spin_orb, ONLY : domag - USE wvfct, ONLY : nbnd, npwx, g2kin, et - USE scf, ONLY : rho - USE uspp, ONLY : okvan, vkb - USE uspp_param, ONLY : upf, nhm, nh + USE wvfct, ONLY : nbnd, npwx, et + USE scf, ONLY : rho, vrs + USE uspp, ONLY : okvan, vkb, deeq_nc + USE uspp_param, ONLY : upf, nhm USE noncollin_module, ONLY : noncolin, npol, nspin_mag USE paw_variables, ONLY : okpaw USE paw_onecenter, ONLY : paw_dpotential @@ -60,7 +58,7 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) iudvscf, iuint3paw, lint3paw USE units_lr, ONLY : iuwfc, lrwfc USE output, ONLY : fildrho, fildvscf - USE phus, ONLY : becsumort + USE phus, ONLY : becsumort, alphap, int1_nc USE modes, ONLY : npertx, npert, u, t, tmq USE recover_mod, ONLY : read_rec, write_rec ! used to write fildrho: @@ -68,18 +66,20 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) USE save_ph, ONLY : tmp_dir_save ! used oly to write the restart file USE mp_pools, ONLY : inter_pool_comm - USE mp_bands, ONLY : intra_bgrp_comm, ntask_groups, me_bgrp + USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp USE mp, ONLY : mp_sum USE efermi_shift, ONLY : ef_shift, ef_shift_paw, def - USE lrus, ONLY : int3_paw + USE lrus, ONLY : int3_paw, becp1, int3_nc USE lr_symm_base, ONLY : irotmq, minus_q, nsymq, rtau USE eqv, ONLY : dvpsi, dpsi, evq USE qpoint, ONLY : xq, nksq, ikks, ikqs - USE control_lr, ONLY : alpha_pv, nbnd_occ, lgamma - USE dv_of_drho_lr + USE qpoint_aux, ONLY : ikmks, ikmkmqs, becpt, alphapt + USE control_lr, ONLY : nbnd_occ, lgamma + USE dv_of_drho_lr, ONLY : dv_of_drho USE fft_helper_subroutines USE fft_interfaces, ONLY : fft_interpolate USE ldaU, ONLY : lda_plus_u + USE nc_mag_aux, ONLY : int1_nc_save, deeq_nc_save, int3_save implicit none @@ -93,11 +93,12 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) real(DP) , allocatable :: h_diag (:,:) ! h_diag: diagonal part of the Hamiltonian - real(DP) :: thresh, anorm, averlt, dr2 + real(DP) :: thresh, anorm, averlt, dr2, rsign ! thresh: convergence threshold ! anorm : the norm of the error ! averlt: average number of iterations ! dr2 : self-consistency error + ! rsign : sign or the term in the magnetization real(DP) :: dos_ef, weight, aux_avg (2) ! Misc variables for metals ! dos_ef: density of states at Ef @@ -110,8 +111,8 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) ! change of rho / scf potential (output) ! change of scf potential (output) complex(DP), allocatable :: ldos (:,:), ldoss (:,:), mixin(:), mixout(:), & - dbecsum (:,:,:,:), dbecsum_nc(:,:,:,:,:), aux1 (:,:), tg_dv(:,:), & - tg_psic(:,:), aux2(:,:), drhoc(:) + dbecsum (:,:,:,:), dbecsum_nc(:,:,:,:,:,:), aux1 (:,:), tg_dv(:,:), & + tg_psic(:,:), aux2(:,:), drhoc(:), dbecsum_aux (:,:,:,:) ! Misc work space ! ldos : local density of states af Ef ! ldoss: as above, without augmentation charges @@ -141,7 +142,11 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) incr, & ! used for tg v_siz, & ! size of the potential ipol, & ! counter on polarization - mode ! mode index + mode, & ! mode index + isolv, & ! counter on linear systems + nsolv, & ! number of linear systems + ikmk, & ! index of mk + ikmkmq ! index of mk-mq integer :: npw, npwq integer :: iq_dummy @@ -156,6 +161,9 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) ! ! This routine is task group aware ! + nsolv=1 + IF (noncolin.AND.domag) nsolv=2 + allocate (dvscfin ( dfftp%nnr , nspin_mag , npe)) if (doublegrid) then allocate (dvscfins (dffts%nnr , nspin_mag , npe)) @@ -169,12 +177,18 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) allocate (mixin(dfftp%nnr*nspin_mag*npe+(nhm*(nhm+1)*nat*nspin_mag*npe)/2) ) allocate (mixout(dfftp%nnr*nspin_mag*npe+(nhm*(nhm+1)*nat*nspin_mag*npe)/2) ) mixin=(0.0_DP,0.0_DP) + ELSE + ALLOCATE(mixin(1)) ENDIF - IF (noncolin) allocate (dbecsum_nc (nhm,nhm, nat , nspin , npe)) + IF (noncolin) allocate (dbecsum_nc (nhm,nhm, nat , nspin , npe, nsolv)) allocate (aux1 ( dffts%nnr, npol)) allocate (h_diag ( npwx*npol, nbnd)) allocate (aux2(npwx*npol, nbnd)) allocate (drhoc(dfftp%nnr)) + IF (noncolin.AND.domag.AND.okvan) THEN + ALLOCATE (int3_save( nhm, nhm, nat, nspin_mag, npe, 2)) + ALLOCATE (dbecsum_aux ( (nhm * (nhm + 1))/2 , nat , nspin_mag , npe)) + ENDIF incr=1 IF ( dffts%has_task_groups ) THEN ! @@ -260,158 +274,215 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) ! if (lsda) current_spin = isk (ikk) ! - ! read unperturbed wavefunctions 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 - ! ! compute beta functions and kinetic energy for k-point ikq ! needed by h_psi, called by ch_psi_all, called by cgsolve_all ! CALL init_us_2 (npwq, igk_k(1,ikq), xk (1, ikq), vkb) CALL g2_kin (ikq) ! - ! compute preconditioning matrix h_diag used by cgsolve_all + ! Start the loop on the two linear systems, one at B and one at + ! -B ! - CALL h_prec (ik, evq, h_diag) - ! - do ipert = 1, npe - mode = imode0 + ipert - nrec = (ipert - 1) * nksq + ik - ! - ! and now adds the contribution of the self consistent term - ! - if (where_rec =='solve_lint'.or.iter>1) then - ! - ! After the first iteration dvbare_q*psi_kpoint is read from file - ! - call get_buffer (dvpsi, lrbar, iubar, nrec) - ! - ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint - ! dvscf_q from previous iteration (mix_potential) - ! - call start_clock ('vpsifft') - 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 = 1, nbnd_occ (ikk), incr - IF( dffts%has_task_groups ) THEN - call cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd_occ (ikk) ) - call apply_dpot(v_siz, tg_psic, tg_dv, 1) - call cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd_occ (ikk)) - 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 stop_clock ('vpsifft') - ! - ! In the case of US pseudopotentials there is an additional - ! selfconsist term which comes from the dependence of D on - ! V_{eff} on the bare change of the potential - ! - call adddvscf (ipert, ik) - ! - ! DFPT+U: add to dvpsi the scf part of the response - ! Hubbard potential dV_hub - ! - if (lda_plus_u) call adddvhubscf (ipert, ik) - ! - else - ! - ! At the first iteration dvbare_q*psi_kpoint is calculated - ! and written to file. - ! - call dvqpsi_us (ik, u (1, mode),.false. ) - ! - ! DFPT+U: At the first ph iteration the bare perturbed - ! Hubbard potential dvbare_hub_q * psi_kpoint - ! is calculated and added to dvpsi. - ! - if (lda_plus_u) call dvqhub_barepsi_us (ik, u(1,mode)) - ! - call save_buffer (dvpsi, lrbar, iubar, nrec) - ! - endif - ! - ! Ortogonalize dvpsi to valence states: ps = - ! Apply -P_c^+. - ! - CALL orthogonalize(dvpsi, evq, ikk, ikq, dpsi, npwq, .false.) - ! - if (where_rec=='solve_lint'.or.iter > 1) then - ! - ! starting value for delta_psi is read from iudwf - ! - call get_buffer( dpsi, lrdwf, iudwf, nrec) - ! - ! threshold for iterative solution of the linear system - ! - thresh = min (1.d-1 * sqrt (dr2), 1.d-2) - else - ! - ! At the first iteration dpsi and dvscfin are set to zero - ! - dpsi(:,:) = (0.d0, 0.d0) - dvscfin (:, :, ipert) = (0.d0, 0.d0) - ! - ! starting threshold for iterative solution of the linear system - ! - thresh = 1.0d-2 - endif - - ! - ! iterative solution of the linear system (H-eS)*dpsi=dvpsi, - ! dvpsi=-P_c^+ (dvbare+dvscf)*psi , dvscf fixed. - ! - conv_root = .true. - - call cgsolve_all (ch_psi_all, cg_psi, et(1,ikk), dvpsi, dpsi, & - h_diag, npwx, npwq, thresh, ik, lter, conv_root, & - anorm, nbnd_occ(ikk), npol ) - - ltaver = ltaver + lter - lintercall = lintercall + 1 - if (.not.conv_root) WRITE( stdout, '(5x,"kpoint",i4," ibnd",i4, & - & " solve_linter: root not converged ",es10.3)') & - & ik , ibnd, anorm - ! - ! writes delta_psi on iunit iudwf, k=kpoint, - ! - ! if (nksq.gt.1 .or. npert(irr).gt.1) - call save_buffer (dpsi, lrdwf, iudwf, nrec) - ! - ! calculates dvscf, sum over k => dvscf_q_ipert - ! - weight = wk (ikk) - IF (noncolin) THEN - call incdrhoscf_nc(drhoscf(1,1,ipert),weight,ik, & - dbecsum_nc(1,1,1,1,ipert), dpsi) + DO isolv=1, nsolv + IF (isolv==2) THEN + ikmk = ikmks(ik) + ikmkmq = ikmkmqs(ik) + rsign=-1.0_DP ELSE - call incdrhoscf (drhoscf(1,current_spin,ipert), weight, ik, & - dbecsum(1,1,current_spin,ipert), dpsi) - END IF - ! on perturbations - enddo + ikmk=ikk + ikmkmq=ikq + rsign=1.0_DP + ENDIF + ! + ! read unperturbed wavefunctions psi(k) and psi(k+q) + ! + if (nksq.gt.1.OR.nsolv==2) then + if (lgamma) then + call get_buffer (evc, lrwfc, iuwfc, ikmk) + else + call get_buffer (evc, lrwfc, iuwfc, ikmk) + call get_buffer (evq, lrwfc, iuwfc, ikmkmq) + end if + endif + ! + ! compute preconditioning matrix h_diag used by cgsolve_all + ! + CALL h_prec (ik, evq, h_diag) + ! + do ipert = 1, npe + mode = imode0 + ipert + nrec = (ipert - 1) * nksq + ik + (isolv-1) * npe * nksq + ! + ! and now adds the contribution of the self consistent term + ! + if (where_rec =='solve_lint'.or.iter>1) then + ! + ! After the first iteration dvbare_q*psi_kpoint is read from file + ! + call get_buffer (dvpsi, lrbar, iubar, nrec) + ! + ! calculates dvscf_q*psi_k in G_space, for all bands, k=kpoint + ! dvscf_q from previous iteration (mix_potential) + ! + call start_clock ('vpsifft') + ! + ! change the sign of the magnetic field if required + ! + IF (isolv==2) THEN + dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert) + IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,2) + ENDIF + ! + ! Set the potential for task groups + ! + 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 = 1, nbnd_occ (ikk), incr + IF( dffts%has_task_groups ) THEN + call cft_wave_tg (ik, evc, tg_psic, 1, v_siz, ibnd, nbnd_occ (ikk) ) + call apply_dpot(v_siz, tg_psic, tg_dv, 1) + call cft_wave_tg (ik, aux2, tg_psic, -1, v_siz, ibnd, nbnd_occ (ikk)) + 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 stop_clock ('vpsifft') + ! + ! In the case of US pseudopotentials there is an additional + ! selfconsist term which comes from the dependence of D on + ! V_{eff} on the bare change of the potential + ! + IF (isolv==1) THEN + call adddvscf_ph_mag (ipert, ik, becp1) + ! + ! DFPT+U: add to dvpsi the scf part of the response + ! Hubbard potential dV_hub + ! + if (lda_plus_u) call adddvhubscf (ipert, ik) + ELSE + call adddvscf_ph_mag (ipert, ik, becpt) + END IF + ! + ! reset the original magnetic field if it was changed + ! + IF (isolv==2) THEN + dvscfins(:,2:4,ipert)=-dvscfins(:,2:4,ipert) + IF (okvan) int3_nc(:,:,:,:,ipert)=int3_save(:,:,:,:,ipert,1) + ENDIF + ! + else + ! + ! At the first iteration dvbare_q*psi_kpoint is calculated + ! and written to file. + ! + IF (isolv==1) THEN + call dvqpsi_us (ik, u (1, mode),.false., becp1, alphap ) + ! + ! DFPT+U: At the first ph iteration the bare perturbed + ! Hubbard potential dvbare_hub_q * psi_kpoint + ! is calculated and added to dvpsi. + ! + if (lda_plus_u) call dvqhub_barepsi_us (ik, u(1,mode)) + ! + ELSE + IF (okvan) THEN + deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2) + int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,2) + ENDIF + call dvqpsi_us (ik, u (1, mode),.false., becpt, alphapt) + IF (okvan) THEN + deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) + int1_nc(:,:,:,:,:)=int1_nc_save(:,:,:,:,:,1) + ENDIF + ENDIF + call save_buffer (dvpsi, lrbar, iubar, nrec) + ! + endif + ! + ! Ortogonalize dvpsi to valence states: ps = + ! Apply -P_c^+. + ! + CALL orthogonalize(dvpsi, evq, ikmk, ikmkmq, dpsi, npwq, .false.) + ! + if (where_rec=='solve_lint'.or.iter > 1) then + ! + ! starting value for delta_psi is read from iudwf + ! + call get_buffer( dpsi, lrdwf, iudwf, nrec) + ! + ! threshold for iterative solution of the linear system + ! + thresh = min (1.d-1 * sqrt (dr2), 1.d-2) + else + ! + ! At the first iteration dpsi and dvscfin are set to zero + ! + dpsi(:,:) = (0.d0, 0.d0) + dvscfin (:, :, ipert) = (0.d0, 0.d0) + ! + ! starting threshold for iterative solution of the linear system + ! + thresh = 1.0d-2 + endif + ! + ! iterative solution of the linear system (H-eS)*dpsi=dvpsi, + ! dvpsi=-P_c^+ (dvbare+dvscf)*psi , dvscf fixed. + ! + IF (isolv==2) THEN + vrs(:,2:4)=-vrs(:,2:4) + IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,2) + ENDIF + conv_root = .true. + + call cgsolve_all (ch_psi_all, cg_psi, et(1,ikmk), dvpsi, dpsi, & + h_diag, npwx, npwq, thresh, ik, lter, conv_root, & + anorm, nbnd_occ(ikk), npol ) + + IF (isolv==2) THEN + vrs(:,2:4)=-vrs(:,2:4) + IF (okvan) deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1) + ENDIF + + ltaver = ltaver + lter + lintercall = lintercall + 1 + if (.not.conv_root) WRITE( stdout, '(5x,"kpoint",i4," ibnd",i4, & + & " solve_linter: root not converged ",es10.3)') & + & ik , ibnd, anorm + ! + ! writes delta_psi on iunit iudwf, k=kpoint, + ! + ! if (nksq.gt.1 .or. npert(irr).gt.1) + call save_buffer (dpsi, lrdwf, iudwf, nrec) + ! + ! calculates dvscf, sum over k => dvscf_q_ipert + ! + weight = wk (ikk) + IF (nsolv==2) weight=weight/2.0_DP + IF (noncolin) THEN + call incdrhoscf_nc(drhoscf(1,1,ipert),weight,ik, & + dbecsum_nc(1,1,1,1,ipert,isolv), dpsi, rsign) + ELSE + call incdrhoscf (drhoscf(1,current_spin,ipert), weight, ik, & + dbecsum(1,1,current_spin,ipert), dpsi) + END IF + ! on perturbations + enddo + ! on isolv + END DO ! on k-points enddo ! @@ -436,7 +507,15 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) ! ! In the noncolinear, spin-orbit case rotate dbecsum ! - IF (noncolin.and.okvan) CALL set_dbecsum_nc(dbecsum_nc, dbecsum, npe) + IF (noncolin.and.okvan) THEN + CALL set_dbecsum_nc(dbecsum_nc, dbecsum, npe) + IF (nsolv==2) THEN + dbecsum_aux=(0.0_DP,0.0_DP) + CALL set_dbecsum_nc(dbecsum_nc(1,1,1,1,1,2), dbecsum_aux, npe) + dbecsum(:,:,1,:)=dbecsum(:,:,1,:)+dbecsum_aux(:,:,1,:) + dbecsum(:,:,2:4,:)=dbecsum(:,:,2:4,:)-dbecsum_aux(:,:,2:4,:) + ENDIF + ENDIF ! ! Now we compute for all perturbations the total charge and potential ! @@ -539,12 +618,49 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) ! calculate here the change of the D1-~D1 coefficients due to the phonon ! perturbation ! - IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) - ! - ! with the new change of the potential we compute the integrals - ! of the change of potential and Q - ! - call newdq (dvscfin, npe) + IF (okvan) THEN + IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) + ! + ! with the new change of the potential we compute the integrals + ! of the change of potential and Q + ! + call newdq (dvscfin, npe) + ! + ! In the noncollinear magnetic case computes the int3 coefficients with + ! the opposite sign of the magnetic field. They are saved in int3_save, + ! that must have been allocated by the calling routine + ! + IF (noncolin.AND.domag) THEN + int3_save(:,:,:,:,:,1)=int3_nc(:,:,:,:,:) + IF (okpaw) rho%bec(:,:,2:4)=-rho%bec(:,:,2:4) + DO ipert=1,npe + dvscfin(:,2:4,ipert)=-dvscfin(:,2:4,ipert) + IF (okpaw) dbecsum(:,:,2:4,ipert)=-dbecsum(:,:,2:4,ipert) + ENDDO + ! + ! if needed recompute the paw coeffients with the opposite sign of + ! the magnetic field + ! + IF (okpaw) CALL PAW_dpotential(dbecsum,rho%bec,int3_paw,npe) + ! + ! here compute the int3 integrals + ! + CALL newdq (dvscfin, npe) + int3_save(:,:,:,:,:,2)=int3_nc(:,:,:,:,:) + ! + ! restore the correct sign of the magnetic field. + ! + DO ipert=1,npe + dvscfin(:,2:4,ipert)=-dvscfin(:,2:4,ipert) + IF (okpaw) dbecsum(:,:,2:4,ipert)=-dbecsum(:,:,2:4,ipert) + ENDDO + IF (okpaw) rho%bec(:,:,2:4)=-rho%bec(:,:,2:4) + ! + ! put into int3_nc the coefficient with +B + ! + int3_nc(:,:,:,:,:)=int3_save(:,:,:,:,:,1) + ENDIF + END IF #if defined(__MPI) aux_avg (1) = DBLE (ltaver) aux_avg (2) = DBLE (lintercall) @@ -621,6 +737,10 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf) DEALLOCATE( tg_dv ) DEALLOCATE( tg_psic ) ENDIF + IF (noncolin.AND.domag.AND.okvan) THEN + DEALLOCATE (int3_save) + DEALLOCATE (dbecsum_aux) + ENDIF call stop_clock ('solve_linter') diff --git a/PHonon/PH/sym_dmag.f90 b/PHonon/PH/sym_dmag.f90 index cfba17daa..8dc8a03a3 100644 --- a/PHonon/PH/sym_dmag.f90 +++ b/PHonon/PH/sym_dmag.f90 @@ -141,10 +141,13 @@ subroutine sym_dmag (nper, irr, dmagtosym) g2 (isym) = 0.d0 g3 (isym) = 0.d0 do ipol = 1, 3 - g1 (isym) = g1 (isym) + gi (ipol, isym) * in1 * at (ipol, 1) - g2 (isym) = g2 (isym) + gi (ipol, isym) * in2 * at (ipol, 2) - g3 (isym) = g3 (isym) + gi (ipol, isym) * in3 * at (ipol, 3) + g1 (isym) = g1 (isym) + gi (ipol, isym) * at (ipol, 1) + g2 (isym) = g2 (isym) + gi (ipol, isym) * at (ipol, 2) + g3 (isym) = g3 (isym) + gi (ipol, isym) * at (ipol, 3) enddo + g1 (isym) = NINT(g1(isym))*in1 + g2 (isym) = NINT(g2(isym))*in2 + g3 (isym) = NINT(g3(isym))*in3 term (1, isym) = CMPLX(cos (g1 (isym) ), sin (g1 (isym) ) ,kind=DP) term (2, isym) = CMPLX(cos (g2 (isym) ), sin (g2 (isym) ) ,kind=DP) term (3, isym) = CMPLX(cos (g3 (isym) ), sin (g3 (isym) ) ,kind=DP) @@ -189,6 +192,9 @@ subroutine sym_dmag (nper, irr, dmagtosym) at(kpol,2)*magrot(2) + & at(kpol,3)*magrot(3) enddo + if (t_rev(isym) == 1) then + mag(:) = conjg(mag(:)) + end if dmagsym(i,j,k,1,ipert)=dmagsym(i,j,k,1,ipert)+mag(1) dmagsym(i,j,k,2,ipert)=dmagsym(i,j,k,2,ipert)+mag(2) dmagsym(i,j,k,3,ipert)=dmagsym(i,j,k,3,ipert)+mag(3) diff --git a/PHonon/PH/symdvscf.f90 b/PHonon/PH/symdvscf.f90 index 6b2f4a7b6..2a15b107d 100644 --- a/PHonon/PH/symdvscf.f90 +++ b/PHonon/PH/symdvscf.f90 @@ -16,7 +16,7 @@ subroutine symdvscf (nper, irr, dvtosym) USE constants, ONLY: tpi USE fft_base, ONLY: dfftp USE cell_base, ONLY : at - USE symm_base, ONLY : s, ft + USE symm_base, ONLY : s, ft, t_rev USE noncollin_module, ONLY : nspin_lsda, nspin_mag USE modes, ONLY : t, tmq @@ -27,7 +27,7 @@ subroutine symdvscf (nper, irr, dvtosym) integer :: nper, irr ! the number of perturbations ! the representation under conside - integer :: ftau(3,48) + integer :: ft(3,48) complex(DP) :: dvtosym (dfftp%nr1x, dfftp%nr2x, dfftp%nr3x, nspin_mag, nper) ! the potential to be symmetrized @@ -37,7 +37,7 @@ subroutine symdvscf (nper, irr, dvtosym) ! counters real(DP) :: gf(3), n(3) ! temp variables - complex(DP), allocatable :: dvsym (:,:,:,:) + complex(DP), allocatable :: dvsym (:,:,:,:), add_dvsym(:) ! the symmetrized potential complex(DP) :: aux2, term (3, 48), phase (48) ! auxiliary space @@ -48,15 +48,16 @@ subroutine symdvscf (nper, irr, dvtosym) call start_clock ('symdvscf') allocate (dvsym( dfftp%nr1x , dfftp%nr2x , dfftp%nr3x , nper)) + allocate (add_dvsym(nper)) ! ! if necessary we symmetrize with respect to S(irotmq)*q = -q + Gi ! n(1) = tpi / DBLE (dfftp%nr1) n(2) = tpi / DBLE (dfftp%nr2) n(3) = tpi / DBLE (dfftp%nr3) - ftau(1,1:nsymq) = NINT ( ft(1,1:nsymq)*dfftp%nr1 ) - ftau(2,1:nsymq) = NINT ( ft(2,1:nsymq)*dfftp%nr2 ) - ftau(3,1:nsymq) = NINT ( ft(3,1:nsymq)*dfftp%nr3 ) + ft(1,1:nsymq) = NINT ( ft(1,1:nsymq)*dfftp%nr1 ) + ft(2,1:nsymq) = NINT ( ft(2,1:nsymq)*dfftp%nr2 ) + ft(3,1:nsymq) = NINT ( ft(3,1:nsymq)*dfftp%nr3 ) if (minus_q) then gf(:) = gimq (1) * at (1, :) * n(:) + & gimq (2) * at (2, :) * n(:) + & @@ -67,7 +68,7 @@ subroutine symdvscf (nper, irr, dvtosym) do k = 1, dfftp%nr3 do j = 1, dfftp%nr2 do i = 1, dfftp%nr1 - CALL ruotaijk (s(1,1,irotmq), ftau(1,irotmq), i, j, k, & + CALL ruotaijk (s(1,1,irotmq), ft(1,irotmq), i, j, k, & dfftp%nr1, dfftp%nr2, dfftp%nr3, ri, rj, rk) do ipert = 1, nper @@ -99,7 +100,7 @@ subroutine symdvscf (nper, irr, dvtosym) gi (3,isym) * at (3, :) * n(:) term (:, isym) = CMPLX(cos (gf (:) ), sin (gf (:) ) ,kind=DP) enddo - + do is = 1, nspin_lsda dvsym(:,:,:,:) = (0.d0, 0.d0) do isym = 1, nsymq @@ -110,16 +111,23 @@ subroutine symdvscf (nper, irr, dvtosym) do i = 1, dfftp%nr1 do isym = 1, nsymq irot = isym - CALL ruotaijk (s(1,1,irot), ftau(1,irot), i, j, k, & + CALL ruotaijk (s(1,1,irot), ft(1,irot), i, j, k, & dfftp%nr1, dfftp%nr2, dfftp%nr3, ri, rj, rk) - + add_dvsym(:) = (0.d0, 0.d0) do ipert = 1, nper do jpert = 1, nper - dvsym (i, j, k, ipert) = dvsym (i, j, k, ipert) + & - t (jpert, ipert, irot, irr) * & - dvtosym (ri, rj, rk, is, jpert) * phase (isym) + add_dvsym(ipert) = add_dvsym(ipert) + t (jpert, ipert, irot, irr) * & + dvtosym (ri, rj, rk, is, jpert) * phase (isym) + !dvsym (i, j, k, ipert) = dvsym (i, j, k, ipert) + & + ! t (jpert, ipert, irot, irr) * & + ! dvtosym (ri, rj, rk, is, jpert) * phase (isym) enddo enddo + if (t_rev(isym)==0) then + dvsym (i, j, k, :) = dvsym (i, j, k, :) + add_dvsym(:) + else + dvsym (i, j, k, :) = dvsym (i, j, k, :) + conjg(add_dvsym(:)) + end if enddo do isym = 1, nsymq phase (isym) = phase (isym) * term (1, isym) @@ -140,6 +148,7 @@ subroutine symdvscf (nper, irr, dvtosym) enddo deallocate (dvsym) + deallocate (add_dvsym) call stop_clock ('symdvscf') return diff --git a/PHonon/PH/symdynph_gq.f90 b/PHonon/PH/symdynph_gq.f90 index 181f61d47..e2449df16 100644 --- a/PHonon/PH/symdynph_gq.f90 +++ b/PHonon/PH/symdynph_gq.f90 @@ -18,6 +18,7 @@ subroutine symdynph_gq_new (xq, phi, s, invs, rtau, irt, nsymq, & ! USE kinds, only : DP USE constants, ONLY: tpi + USE symm_base, ONLY : t_rev implicit none ! ! The dummy variables @@ -128,9 +129,15 @@ subroutine symdynph_gq_new (xq, phi, s, invs, rtau, irt, nsymq, & do jpol = 1, 3 do kpol = 1, 3 do lpol = 1, 3 - work (ipol, jpol) = work (ipol, jpol) + & - s (ipol, kpol, irot) * s (jpol, lpol, irot) & - * phi (kpol, lpol, sna, snb) * faseq (isymq) + IF (t_rev(isymq)==1) THEN + work (ipol, jpol) = work (ipol, jpol) + & + s (ipol, kpol, irot) * s (jpol, lpol, irot) & + * CONJG(phi (kpol, lpol, sna, snb) * faseq (isymq)) + ELSE + work (ipol, jpol) = work (ipol, jpol) + & + s (ipol, kpol, irot) * s (jpol, lpol, irot) & + * phi (kpol, lpol, sna, snb) * faseq (isymq) + ENDIF enddo enddo enddo @@ -145,9 +152,15 @@ subroutine symdynph_gq_new (xq, phi, s, invs, rtau, irt, nsymq, & phi (ipol, jpol, sna, snb) = (0.d0, 0.d0) do kpol = 1, 3 do lpol = 1, 3 - phi (ipol, jpol, sna, snb) = phi (ipol, jpol, sna, snb) & - + s (ipol, kpol, invs (irot) ) * s (jpol, lpol, invs (irot) ) & + IF (t_rev(isymq)==1) THEN + phi(ipol,jpol,sna,snb)=phi(ipol,jpol,sna,snb) & + + s(ipol,kpol,invs(irot))*s(jpol,lpol,invs(irot))& + * CONJG(work (kpol, lpol)*faseq (isymq)) + ELSE + phi(ipol,jpol,sna,snb)=phi(ipol,jpol,sna,snb) & + + s(ipol,kpol,invs(irot))*s(jpol,lpol,invs(irot))& * work (kpol, lpol) * CONJG(faseq (isymq) ) + ENDIF enddo enddo enddo diff --git a/PHonon/PH/zstar_eu.f90 b/PHonon/PH/zstar_eu.f90 index f664387d6..fe340e226 100644 --- a/PHonon/PH/zstar_eu.f90 +++ b/PHonon/PH/zstar_eu.f90 @@ -37,6 +37,8 @@ subroutine zstar_eu USE mp_bands, ONLY : intra_bgrp_comm USE mp, ONLY : mp_sum USE ldaU, ONLY : lda_plus_u + USE lrus, ONLY : becp1 + USE phus, ONLY : alphap implicit none @@ -66,7 +68,7 @@ subroutine zstar_eu ! ! recalculate DeltaV*psi(ion) for mode nu ! - call dvqpsi_us (ik, u (1, mode), .not.okvan) + call dvqpsi_us (ik, u (1, mode), .not.okvan, becp1, alphap) ! ! DFPT+U: add the bare variation of the Hubbard potential !