scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Developments on the PH package routines to implement the calculation of phonon frequencies in the noncollinear magnetic case

This commit is contained in:
andreaurru247 2020-01-29 11:37:15 +01:00 committed by andreaurru247
parent 0f2315a927
commit 3ebf9c1831
30 changed files with 1139 additions and 307 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

127
PHonon/PH/adddvscf_ph_mag.f90 Normal file
View File

@ -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

25
PHonon/PH/allocate_phq.f90
View File

@ -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.

56
PHonon/PH/apply_trev.f90 Normal file
View File

@ -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

179
PHonon/PH/c_bands_ph.f90 Normal file
View File

@ -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

23
PHonon/PH/deallocate_phq.f90
View File

@ -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)

8
PHonon/PH/do_phonon.f90
View File

@ -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

35
PHonon/PH/drhodv.f90
View File

@ -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,9 +97,11 @@ 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 isolv=1, nsolv
do mu = 1, nper
nrec = (mu - 1) * nksq + ik
if (nksq > 1 .or. nper > 1) call get_buffer(dpsi, lrdwf, iudwf, nrec)
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)
@ -117,7 +126,22 @@ subroutine drhodv (nu_i0, nper, drhoscf)
CALL becscal(fact,dalpq(ipol,ipert),nkb,nbnd)
ENDDO
ENDDO
call drhodvnl (ik, ikk, nper, nu_i0, dynwrk, dbecq, dalpq)
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
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
!

11
PHonon/PH/drhodvnl.f90
View File

@ -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 <dpsi|dv-e ds|psi> 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)

7
PHonon/PH/dvqpsi_us.f90
View File

@ -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

9
PHonon/PH/dvqpsi_us_only.f90
View File

@ -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
!

4
PHonon/PH/dynmatrix.f90
View File

@ -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
!

7
PHonon/PH/elphon.f90
View File

@ -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
!

5
PHonon/PH/ep_matrix_element_wannier.f90
View File

@ -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

2
PHonon/PH/incdrhous_nc.f90
View File

@ -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)

29
PHonon/PH/initialize_ph.f90
View File

@ -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,6 +26,18 @@ SUBROUTINE initialize_ph()
!
IF ( lgamma ) THEN
!
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))
@ -30,9 +45,22 @@ SUBROUTINE initialize_ph()
ikks(ik) = ik
ikqs(ik) = ik
ENDDO
END IF
!
ELSE
!
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))
@ -40,6 +68,7 @@ SUBROUTINE initialize_ph()
ikks(ik) = 2 * ik - 1
ikqs(ik) = 2 * ik
ENDDO
END IF
!
END IF
!

112
PHonon/PH/non_scf_ph.f90 Normal file
View File

@ -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

1
PHonon/PH/openfilq.f90
View File

@ -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

28
PHonon/PH/phcom.f90
View File

@ -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

32
PHonon/PH/phq_init.f90
View File

@ -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
!
@ -190,9 +192,13 @@ SUBROUTINE phq_init()
!
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)
! 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
@ -222,6 +228,24 @@ SUBROUTINE phq_init()
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
@ -242,7 +266,7 @@ SUBROUTINE phq_init()
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 davcio (evq, lrwfc, iunwfcwann, ikqg, - 1)
call calculate_and_apply_phase(ik, ikqg, igqg, &
npwq_refolded, g_kpq, xk_gamma, evq, .false.)
ENDIF

16
PHonon/PH/phq_setup.f90
View File

@ -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

4
PHonon/PH/raman_mat.f90
View File

@ -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)

13
PHonon/PH/run_nscf.f90
View File

@ -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
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

26
PHonon/PH/set_int12_nc.f90
View File

@ -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

4
PHonon/PH/set_irr.f90
View File

@ -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
!

11
PHonon/PH/set_irr_sym.f90
View File

@ -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) &
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

196
PHonon/PH/solve_linter.f90
View File

@ -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,31 +274,44 @@ 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)
!
! Start the loop on the two linear systems, one at B and one at
! -B
!
DO isolv=1, nsolv
IF (isolv==2) THEN
ikmk = ikmks(ik)
ikmkmq = ikmkmqs(ik)
rsign=-1.0_DP
ELSE
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
nrec = (ipert - 1) * nksq + ik + (isolv-1) * npe * nksq
!
! and now adds the contribution of the self consistent term
!
@ -298,6 +325,16 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! 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))
@ -329,19 +366,31 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! selfconsist term which comes from the dependence of D on
! V_{eff} on the bare change of the potential
!
call adddvscf (ipert, ik)
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.
!
call dvqpsi_us (ik, u (1, mode),.false. )
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
@ -349,6 +398,17 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
!
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
@ -356,7 +416,7 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! Ortogonalize dvpsi to valence states: ps = <evq|dvpsi>
! Apply -P_c^+.
!
CALL orthogonalize(dvpsi, evq, ikk, ikq, dpsi, npwq, .false.)
CALL orthogonalize(dvpsi, evq, ikmk, ikmkmq, dpsi, npwq, .false.)
!
if (where_rec=='solve_lint'.or.iter > 1) then
!
@ -378,17 +438,25 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
!
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,ikk), dvpsi, dpsi, &
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, &
@ -403,15 +471,18 @@ SUBROUTINE solve_linter (irr, imode0, npe, drhoscf)
! 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), dpsi)
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 (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')

12
PHonon/PH/sym_dmag.f90
View File

@ -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)

31
PHonon/PH/symdvscf.f90
View File

@ -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
@ -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) * &
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

17
PHonon/PH/symdynph_gq.f90
View File

@ -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
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

4
PHonon/PH/zstar_eu.f90
View File

@ -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
!