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 LR_Modules routines to implement the calculation of the phonon frequencies in the noncollinear magnetic case. Note: the calculation of the Born effective charges does not work in this case yet, it will be fixed by a further development

This commit is contained in:
andreaurru247 2020-01-29 11:44:02 +01:00 committed by andreaurru247
parent 3ebf9c1831
commit d7bedb5acf
10 changed files with 170 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
LR_Modules/Makefile
View File

@ -36,6 +36,7 @@ setup_offset_beta.o \
setup_nscf.o \
set_dbecsum_nc.o \
set_int3_nc.o \
set_kplusq_nc.o \
smallgq.o \
lr_sm1_psi.o \
lr_addusddens.o \

9
LR_Modules/addusdbec_nc.f90
View File

@ -6,7 +6,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------
subroutine addusdbec_nc (ik, wgt, dpsi, dbecsum_nc)
subroutine addusdbec_nc (ik, wgt, dpsi, dbecsum_nc, becp1)
!----------------------------------------------------------------------
!
! This routine adds to the dbecsum the term which correspond to this
@ -16,15 +16,14 @@ subroutine addusdbec_nc (ik, wgt, dpsi, dbecsum_nc)
USE kinds, ONLY : DP
USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp, ntyp => nsp
USE becmod, ONLY : calbec
USE becmod, ONLY : calbec, bec_type
USE wvfct, ONLY : npwx, nbnd
USE uspp, ONLY : nkb, vkb, okvan
USE noncollin_module, ONLY : noncolin, npol
USE uspp_param, ONLY : upf, nh, nhm
USE mp_bands, ONLY : intra_bgrp_comm
USE klist, ONLY : ngk
USE lrus, ONLY : becp1
USE qpoint, ONLY : ikks, ikqs
USE qpoint, ONLY : nksq, ikks, ikqs
USE control_lr, ONLY : nbnd_occ
!
IMPLICIT NONE
@ -34,6 +33,8 @@ subroutine addusdbec_nc (ik, wgt, dpsi, dbecsum_nc)
COMPLEX(DP) :: dbecsum_nc (nhm,nhm,nat,nspin), dpsi(npwx*npol,nbnd)
! inp/out: the sum kv of bec *
! input : contains delta psi
TYPE(bec_type) :: becp1(nksq)
!
INTEGER :: ik
! input: the k point
REAL(DP) :: wgt

33
LR_Modules/incdrhoscf_nc.f90
View File

@ -6,7 +6,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi, rsign_in)
!-----------------------------------------------------------------------
!
! This routine computes the change of the charge density due to the
@ -27,6 +27,8 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
USE klist, ONLY : ngk,igk_k
USE qpoint, ONLY : ikks, ikqs
USE control_lr, ONLY : nbnd_occ
USE qpoint_aux, ONLY : becpt
USE lrus, ONLY : becp1
USE mp_bands, ONLY : me_bgrp, inter_bgrp_comm, ntask_groups
USE mp, ONLY : mp_sum
USE fft_helper_subroutines
@ -37,6 +39,7 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
INTEGER, INTENT(IN) :: ik
! input: the k point
REAL(DP), INTENT(IN) :: weight
REAL(DP), INTENT(IN), OPTIONAL:: rsign_in
! input: the weight of the k point
COMPLEX(DP), INTENT(IN) :: dpsi(npwx*npol,nbnd)
! input: the perturbed wfcs at the given k point
@ -45,6 +48,8 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
!
! here the local variable
!
REAL(DP) :: rsign
! the sign in front of the response of the magnetization density
REAL(DP) :: wgt
! the effective weight of the k point
!
@ -59,6 +64,14 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
INTEGER :: ntgrp, right_inc
! counters
!
!
IF (PRESENT(rsign_in)) THEN
rsign = rsign_in
ELSE
rsign = 1.
END IF
!
!
CALL start_clock ('incdrhoscf')
!
ALLOCATE (dpsic(dffts%nnr, npol))
@ -130,11 +143,11 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
enddo
IF (domag) THEN
do ir = 1, dffts%nr1x * dffts%nr2x * dffts%my_nr3p
tg_drho(ir,2)= tg_drho(ir,2) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
tg_drho(ir,2)= tg_drho(ir,2) + rsign *wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
+ CONJG(tg_psi(ir,2))*tg_dpsi(ir,1) )
tg_drho(ir,3)= tg_drho(ir,3) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
tg_drho(ir,3)= tg_drho(ir,3) + rsign *wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
- CONJG(tg_psi(ir,2))*tg_dpsi(ir,1) ) * (0.d0,-1.d0)
tg_drho(ir,4)= tg_drho(ir,4) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,1) &
tg_drho(ir,4)= tg_drho(ir,4) + rsign *wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,1) &
- CONJG(tg_psi(ir,2))*tg_dpsi(ir,2) )
enddo
ENDIF
@ -176,11 +189,11 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
enddo
IF (domag) THEN
do ir = 1, dffts%nnr
drhoscf(ir,2)=drhoscf (ir,2) + wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
drhoscf(ir,2)=drhoscf (ir,2) + rsign *wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
+ CONJG(psi(ir,2))*dpsic(ir,1) )
drhoscf(ir,3)=drhoscf (ir,3) + wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
drhoscf(ir,3)=drhoscf (ir,3) + rsign *wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
- CONJG(psi(ir,2))*dpsic(ir,1) ) * (0.d0,-1.d0)
drhoscf(ir,4)=drhoscf (ir,4) + wgt * (CONJG(psi(ir,1))*dpsic(ir,1) &
drhoscf(ir,4)=drhoscf (ir,4) + rsign *wgt * (CONJG(psi(ir,1))*dpsic(ir,1) &
- CONJG(psi(ir,2))*dpsic(ir,2) )
enddo
END IF
@ -192,7 +205,11 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
! Ultrasoft contribution
! Calculate dbecsum_nc = <evc|vkb><vkb|dpsi>
!
CALL addusdbec_nc (ik, weight, dpsi, dbecsum)
IF (rsign==1) THEN
CALL addusdbec_nc (ik, weight, dpsi, dbecsum, becp1)
ELSE
CALL addusdbec_nc (ik, weight, dpsi, dbecsum, becpt)
ENDIF
!
DEALLOCATE(psi)
DEALLOCATE(dpsic)

16
LR_Modules/lrcom.f90
View File

@ -35,6 +35,22 @@ MODULE qpoint
!
END MODULE qpoint
!
!
!
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 control_lr
!
USE kinds, ONLY : DP

81
LR_Modules/set_kplusq_nc.f90 Normal file
View File

@ -0,0 +1,81 @@
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine set_kplusq_nc (xk, wk, xq, nks, npk)
!-----------------------------------------------------------------------
! This routine sets the k and k+q points (with zero weight) used in
! the preparatory run for a linear response calculation.
!
! on input: xk and wk contain k-points and corresponding weights
!
! on output: the number of points is doubled and xk and wk in the
! odd positions are the original ones; those in the
! even positions are the corresponding k+q values.
! the gamma point is treated in a special way. No change is done
! to the k-points
!
USE kinds, only : DP
implicit none
!
! First the dummy variables
!
integer :: npk, nks
! input-output: maximum allowed number of k
! input-output: starting and ending number of
real(DP) :: xk (3, npk), wk (npk), eps, xq (3)
! input-output: coordinates of k points
! input-output: weights of k points
! the smallest xq
! input: coordinates of a q-point
!
! And then the local variables
!
logical :: lgamma
! true if xq is the gamma point
integer :: ik, j
! counter on k
! counter
!
eps = 1.d-12
!
! shift the k points in the odd positions and fill the even ones with k+
!
lgamma = abs (xq (1) ) .lt.eps.and.abs (xq (2) ) .lt.eps.and.abs ( &
xq (3) ) .lt.eps
if (.not.lgamma) then
if (4 * nks.gt.npk) call errore ('set_kplusq', 'too many k points', &
& nks)
do ik = nks, 1, - 1
xk (:, 4 * ik - 3) = xk (:, ik)
xk (:, 4 * ik - 2) = xk (:, ik) + xq (:)
xk (:, 4 * ik - 1) = -xk (:, ik)
xk (:, 4 * ik ) = -xk (:, ik) - xq (:)
wk (4 * ik - 3) = wk (ik)
wk (4 * ik - 2) = 0.0_DP
wk (4 * ik - 1) = 0.0_DP
wk (4 * ik ) = 0.d0
enddo
nks = 4 * nks
else
if (2 * nks.gt.npk) call errore ('set_kplusq', 'too many k points', &
& nks)
do ik = nks, 1, - 1
xk (:, 2 * ik - 1) = xk (:, ik)
xk (:, 2 * ik ) = -xk (:, ik)
wk (2 * ik - 1) = wk (ik)
wk (2 * ik ) = 0.0_DP
enddo
nks = 2 * nks
ENDIF
return
end subroutine set_kplusq_nc

2
LR_Modules/set_small_group_of_q.f90
View File

@ -83,6 +83,7 @@ SUBROUTINE smallg_q (xq, modenum, at, bg, nrot, s, sym, minus_q)
! input-output variables
!
USE kinds, ONLY : DP
USE symm_base, ONLY : t_rev
implicit none
@ -146,6 +147,7 @@ SUBROUTINE smallg_q (xq, modenum, at, bg, nrot, s, sym, minus_q)
raq(ipol) = raq(ipol) + DBLE( s(ipol,jpol,irot) ) * aq( jpol)
enddo
enddo
IF (t_rev(irot)==1) raq=-raq
sym (irot) = eqvect (raq, aq, zero, accep)
!
! if "iswitch.le.-3" (modenum.ne.0) S must be such that Sq=q exactly !

18
LR_Modules/setup_nscf.f90
View File

@ -24,7 +24,6 @@ SUBROUTINE setup_nscf ( newgrid, xq, elph_mat )
!
USE kinds, ONLY : DP
USE parameters, ONLY : npk
USE io_global, ONLY : stdout
USE constants, ONLY : pi, degspin
USE cell_base, ONLY : at, bg
USE ions_base, ONLY : nat, tau, ityp, zv
@ -54,6 +53,7 @@ SUBROUTINE setup_nscf ( newgrid, xq, elph_mat )
LOGICAL, INTENT (IN) :: elph_mat ! used to be passed through a module.
!
REAL (DP), ALLOCATABLE :: rtau (:,:,:)
INTEGER :: t_rev_eff(48), ik
LOGICAL :: magnetic_sym, sym(48)
LOGICAL :: skip_equivalence
LOGICAL, EXTERNAL :: check_para_diag
@ -109,7 +109,9 @@ SUBROUTINE setup_nscf ( newgrid, xq, elph_mat )
! (and possibly in other cases as well) the k-points should not be reduced
!
skip_equivalence = elph_mat
CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev, &
t_rev_eff=0
! Ancora da capire se t_rev_eff sia veramente necessario
CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev_eff, &
bg, nk1*nk2*nk3, k1,k2,k3, nk1,nk2,nk3, nkstot, xk, wk)
endif
@ -123,14 +125,24 @@ SUBROUTINE setup_nscf ( newgrid, xq, elph_mat )
!
! ... add k+q to the list of k
!
CALL set_kplusq( xk, wk, xq, nkstot, npk )
IF (noncolin.AND.domag) THEN
CALL set_kplusq_nc( xk, wk, xq, nkstot, npk)
ELSE
CALL set_kplusq( xk, wk, xq, nkstot, npk)
ENDIF
!
! ... set the granularity for k-point distribution
!
IF ( lgamma ) THEN
!
kunit = 1
IF (noncolin.AND.domag) kunit = 2
!
ELSE
!
kunit = 2
IF (noncolin.AND.domag) kunit = 4
!
ENDIF
!
! ... Map each k point in the irr.-BZ into tetrahedra

8
LR_Modules/smallgq.f90
View File

@ -16,6 +16,7 @@ subroutine set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq)
USE kinds, ONLY : DP
USE cell_base, ONLY : bg, at
USE control_lr, ONLY : lgamma
USE symm_base, ONLY : t_rev
IMPLICIT NONE
@ -73,10 +74,15 @@ subroutine set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq)
aq (jpol)
enddo
enddo
IF (t_rev(isym)==1) raq=-raq
if (.NOT. eqvect (raq, aq, zero, accep) ) CALL errore('set_giq',&
'problems with the input group',1)
do ipol = 1, 3
wrk (ipol) = raq (ipol) - aq (ipol)
IF (t_rev(isym)==1) THEN
wrk (ipol) = aq (ipol) - raq (ipol)
ELSE
wrk (ipol) = raq (ipol) - aq (ipol)
ENDIF
enddo
call cryst_to_cart (1, wrk, bg, 1)
gi (:, isym) = wrk (:)

16
LR_Modules/star_q.f90
View File

@ -6,7 +6,7 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity )
subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity, t_rev_in )
!-----------------------------------------------------------------------
! generate the star of q vectors that are equivalent to the input one
! NB: input s(:,:,1:nsym) must contain all crystal symmetries,
@ -19,6 +19,7 @@ subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity )
real(DP), parameter :: accep=1.e-5_dp
integer, intent(in) :: nsym, s (3, 3, 48), invs(48)
integer, intent(in), optional:: t_rev_in(48)
! nsym matrices of symmetry operations
! invs: list of inverse operation indices
real(DP), intent(in) :: xq (3), at (3, 3), bg (3, 3)
@ -36,11 +37,13 @@ subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity )
logical, intent(in) :: verbosity
! if true prints several messages.
!
integer :: nsq (48), isym, ism1, iq, i
integer :: nsq (48), isym, ism1, iq, i, t_rev(48)
! number of symmetry ops. of bravais lattice
! counters on symmetry ops.
! index of inverse of isym
! counters
! t_rev variable, says if a given simmetry operation
! needs the time reversal to be a symmetry of the crystal.
real(DP) :: saq (3, 48), aq (3), raq (3), zero (3)
! auxiliary list of q (crystal coordinates)
! input q in crystal coordinates
@ -51,6 +54,14 @@ subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity )
logical, external :: eqvect
! function used to compare two vectors
!
!
if (present(t_rev_in)) then
t_rev = t_rev_in
else
t_rev = 0
end if
!
!
zero(:) = 0.d0
!
! go to crystal coordinates
@ -73,6 +84,7 @@ subroutine star_q (xq, at, bg, nsym, s, invs, nq, sxq, isq, imq, verbosity )
+ s (i, 2, ism1) * aq (2) &
+ s (i, 3, ism1) * aq (3)
enddo
IF (t_rev(isym)==1) raq = -raq
do i = 1, 3
sxq (i, 48) = bg (i, 1) * raq (1) &
+ bg (i, 2) * raq (2) &

4
PHonon/PH/Makefile
View File

@ -14,6 +14,7 @@ add_for_charges.o \
add_zstar_ue.o \
add_zstar_ue_us.o \
addcore.o \
adddvscf_ph_mag.o \
addnlcc.o \
addnlcc_zstar_eu_us.o \
addusddens.o \
@ -23,8 +24,10 @@ adddvhubscf.o \
allocate_part.o \
allocate_pert.o \
allocate_phq.o \
apply_trev.o \
yambo.o \
bcast_ph_input.o \
c_bands_ph.o \
cch_psi_all.o \
ccg_psi.o \
check_if_partial_dyn.o \
@ -97,6 +100,7 @@ init_representations.o \
io_dyn_mat.o \
io_dyn_mat_old.o \
io_pattern.o \
non_scf_ph.o \
obsolete.o \
openfilq.o \
phcom.o \