diff --git a/PW/irrek.f90 b/PW/irrek.f90 index 2c8a0094b..1257598ba 100644 --- a/PW/irrek.f90 +++ b/PW/irrek.f90 @@ -17,6 +17,8 @@ subroutine irreducible_BZ (nrot, s, nsym, minus_q, at, bg, npk, nks, & ! of the point group of the Bravais lattice. ! USE kinds, only : DP + USE noncollin_module, ONLY : noncolin + USE spin_orb, ONLY : domag implicit none ! integer, intent(in) :: nrot, nsym, npk, s(3,3,48), t_rev(48) @@ -46,14 +48,20 @@ subroutine irreducible_BZ (nrot, s, nsym, minus_q, at, bg, npk, nks, & ! ! Find the coset in the point group of the Bravais lattice ! - sym(1:nsym) = .true. - sym(nsym+1:)= .false. - call coset (nrot, table, sym, nsym, irg) - ! - ! here we set the k-points in the irreducible wedge of the point grou - ! of the crystal - ! - call irrek (at, bg, nrot, invs, nsym, irg, minus_q, npk, nks, xk, wk, t_rev) + IF (noncolin.AND.domag) THEN + call irrek_nc(at, bg, nrot, invs, nsym, irg, minus_q, npk, nks, xk, & + wk, t_rev) + ELSE + sym(1:nsym) = .true. + sym(nsym+1:)= .false. + call coset (nrot, table, sym, nsym, irg) + ! + ! here we set the k-points in the irreducible wedge of the point grou + ! of the crystal + ! + call irrek (at, bg, nrot, invs, nsym, irg, minus_q, npk, nks, xk, & + wk, t_rev) + ENDIF ! return ! @@ -204,3 +212,152 @@ subroutine irrek (at, bg, nrot, invs, nsym, irg, minus_q, npk, & ! return end subroutine irrek +!----------------------------------------------------------------------- +subroutine irrek_nc (at, bg, nrot, invs, nsym, irg, minus_q, npk, & + nks, xk, wk, t_rev) + !----------------------------------------------------------------------- + ! + ! Given a set of special points in the Irreducible Wedge of some + ! group, finds the equivalent special points in the IW of one of + ! its subgroups. + ! + USE kinds, only : DP + implicit none + ! + integer, intent(in) :: npk, nrot, nsym, invs (3, 3, 48), irg (nrot) + ! maximum number of special points + ! order of the parent point group + ! order of the subgroup + ! inverse of the elements of the symmetry group + ! partition of the elements of the symmetry group into left cosets, + ! as given by SUBROUTINE COSET + integer, intent(inout) :: nks + ! number of special points + integer, intent(in) :: t_rev(48) + real(DP), intent(in) :: at (3, 3), bg (3, 3) + ! basis vectors of the Bravais and reciprocal lattice + real(DP), intent(inout) :: xk (3, npk), wk (npk) + ! special points and weights + logical, intent(in) :: minus_q + ! .true. if symmetries q = -q+G are acceptable + ! + ! here the local variables + ! + integer :: nks0, jk, kpol, irot, jrot, isym, ik, iks, start_k + ! nks0: used to save the initial number of k-points + ! ncos: total number of cosets + real(DP) :: xkg (3), xks (3), xkn(3), one, xk_new(3,npk), wk_new(npk), & + xk_cart(3) + ! coordinates of the k point in crystal axis + ! coordinates of the rotated k point + ! weight of each coset + ! buffer which contains the weight of k points + ! total weight of k-points + logical :: satm + ! true if equivalent point found + + nks0 = nks + nks=0 + start_k=0 + DO jk = 1, nks0 + ! + ! The k point is first computed in crystal axis + ! + ! xkg are the components of xk in the crystal base + xkg (:) = at (1, :) * xk (1, jk) + & + at (2, :) * xk (2, jk) + & + at (3, :) * xk (3, jk) + ! + ! Then it is rotated with each symmetry of the global group. + ! + DO irot = 1, nrot + xks (:) = invs (:, 1, irot) * xkg (1) + & + invs (:, 2, irot) * xkg (2) + & + invs (:, 3, irot) * xkg (3) + ! + ! Now check if there is an operation of the subgroup that + ! makes xks equivalent to some other already found k point + ! + DO jrot=1,nsym + xkn (:) = invs (:, 1, jrot) * xks (1) + & + invs (:, 2, jrot) * xks (2) + & + invs (:, 3, jrot) * xks (3) + IF (t_rev(jrot)==1) xkn =-xkn + DO ik = start_k+1, nks + satm = abs (xk_new (1, ik) - xkn (1) - & + nint (xk_new (1, ik) - xkn (1) ) ) < 1.0d-5 .and. & + abs (xk_new (2, ik) - xkn (2) - & + nint (xk_new (2, ik) - xkn (2) ) ) < 1.0d-5 .and. & + abs (xk_new (3, ik) - xkn (3) - & + nint (xk_new (3, ik) - xkn (3) ) ) < 1.0d-5 + ! + ! .... or equivalent to minus each other when minus_q=.t. + ! + if (minus_q) satm = satm .or. & + abs (xk_new (1, ik) + xkn (1) - & + nint (xk_new (1, ik) + xkn (1) ) ) < 1.0d-5 .and. & + abs (xk_new (2, ik) + xkn (2) - & + nint (xk_new (2, ik) + xkn (2) ) ) < 1.0d-5 .and. & + abs (xk_new (3, ik) + xkn (3) - & + nint (xk_new (3, ik) + xkn (3) ) ) < 1.0d-5 + IF ( satm ) THEN + wk_new(ik) = wk_new(ik) + wk(jk) + GOTO 100 + ENDIF + END DO + END DO + nks=nks+1 + IF (nks > npk) CALL errore('irrek_nc','too many k points',1) + xk_new(:,nks)=xks + wk_new(nks)=wk(jk) +100 CONTINUE + ENDDO + start_k=nks + ENDDO + ! + ! The order of the original k points is preserved + ! + iks=nks0 + DO ik = 1, nks + ! + ! for each new k point found, check if it was in the original list + ! + DO jk=1, nks0 + xkg (:) = at (1, :) * xk (1, jk) + & + at (2, :) * xk (2, jk) + & + at (3, :) * xk (3, jk) + satm = abs (xk_new (1, ik) - xkg (1) - & + nint (xk_new (1, ik) - xkg (1) ) ) < 1.0d-5 .and. & + abs (xk_new (2, ik) - xkg (2) - & + nint (xk_new (2, ik) - xkg (2) ) ) < 1.0d-5 .and. & + abs (xk_new (3, ik) - xkg (3) - & + nint (xk_new (3, ik) - xkg (3) ) ) < 1.0d-5 + IF (satm) THEN +! +! If it was, just update the weight +! + wk(jk)=wk_new(ik) + goto 200 + ENDIF + ENDDO +! +! If it was not, bring xk_new in cartesian coodinates and copy it in the +! first free place available +! + iks=iks+1 + xk_cart (:) = bg (:, 1) * xk_new (1, ik) + & + bg (:, 2) * xk_new (2, ik) + & + bg (:, 3) * xk_new (3, ik) + xk(:,iks)=xk_cart(:) + wk(iks)=wk_new(ik) +200 CONTINUE + ENDDO + IF (iks /= nks ) CALL errore('irrek_nc','Internal problem with k points',1) + ! + ! normalize weights to one + ! + one = SUM (wk(1:nks)) + IF ( one > 0.d0 ) wk(1:nks) = wk(1:nks) / one + ! + RETURN +END SUBROUTINE irrek_nc