1!-*- mode: F90 -*-! 2!------------------------------------------------------------! 3! This file is distributed as part of the Wannier90 code and ! 4! under the terms of the GNU General Public License. See the ! 5! file `LICENSE' in the root directory of the Wannier90 ! 6! distribution, or http://www.gnu.org/copyleft/gpl.txt ! 7! ! 8! The webpage of the Wannier90 code is www.wannier.org ! 9! ! 10! The Wannier90 code is hosted on GitHub: ! 11! ! 12! https://github.com/wannier-developers/wannier90 ! 13!------------------------------------------------------------! 14 15module w90_spin 16 !! Module to compute spin 17 18 use w90_constants, only: dp 19 20 implicit none 21 22 private 23 24 public :: spin_get_moment, spin_get_nk, spin_get_S 25 26contains 27 28 !===========================================================! 29 ! PUBLIC PROCEDURES ! 30 !===========================================================! 31 32 subroutine spin_get_moment 33 !============================================================! 34 ! ! 35 !! Computes the spin magnetic moment by Wannier interpolation 36 ! ! 37 !============================================================! 38 39 use w90_constants, only: dp, pi, cmplx_i 40 use w90_comms, only: on_root, my_node_id, num_nodes, comms_reduce 41 use w90_io, only: io_error, stdout 42 use w90_postw90_common, only: num_int_kpts_on_node, int_kpts, weight 43 use w90_parameters, only: spin_kmesh, wanint_kpoint_file, & 44 nfermi, fermi_energy_list 45 use w90_get_oper, only: get_HH_R, get_SS_R 46 47 integer :: loop_x, loop_y, loop_z, loop_tot 48 real(kind=dp) :: kweight, kpt(3), spn_k(3), spn_all(3), & 49 spn_mom(3), magnitude, theta, phi, conv 50 51 if (nfermi > 1) call io_error('Routine spin_get_moment requires nfermi=1') 52 53 call get_HH_R 54 call get_SS_R 55 56 if (on_root) then 57 write (stdout, '(/,/,1x,a)') '------------' 58 write (stdout, '(1x,a)') 'Calculating:' 59 write (stdout, '(1x,a)') '------------' 60 write (stdout, '(/,3x,a)') '* Spin magnetic moment' 61 end if 62 63 spn_all = 0.0_dp 64 if (wanint_kpoint_file) then 65 66 if (on_root) then 67 write (stdout, '(/,1x,a)') 'Sampling the irreducible BZ only' 68 write (stdout, '(5x,a)') & 69 'WARNING: - IBZ implementation is currently limited to simple cases:' 70 write (stdout, '(5x,a)') & 71 ' Check results against a full BZ calculation!' 72 end if 73 ! 74 ! Loop over k-points on the irreducible wedge of the Brillouin zone, 75 ! read from file 'kpoint.dat' 76 ! 77 do loop_tot = 1, num_int_kpts_on_node(my_node_id) 78 kpt(:) = int_kpts(:, loop_tot) 79 kweight = weight(loop_tot) 80 call spin_get_moment_k(kpt, fermi_energy_list(1), spn_k) 81 spn_all = spn_all + spn_k*kweight 82 end do 83 84 else 85 86 if (on_root) & 87 write (stdout, '(/,1x,a)') 'Sampling the full BZ (not using symmetry)' 88 kweight = 1.0_dp/real(PRODUCT(spin_kmesh), kind=dp) 89 do loop_tot = my_node_id, PRODUCT(spin_kmesh) - 1, num_nodes 90 loop_x = loop_tot/(spin_kmesh(2)*spin_kmesh(3)) 91 loop_y = (loop_tot - loop_x*(spin_kmesh(2)*spin_kmesh(3)))/spin_kmesh(3) 92 loop_z = loop_tot - loop_x*(spin_kmesh(2)*spin_kmesh(3)) & 93 - loop_y*spin_kmesh(3) 94 kpt(1) = (real(loop_x, dp)/real(spin_kmesh(1), dp)) 95 kpt(2) = (real(loop_y, dp)/real(spin_kmesh(2), dp)) 96 kpt(3) = (real(loop_z, dp)/real(spin_kmesh(3), dp)) 97 call spin_get_moment_k(kpt, fermi_energy_list(1), spn_k) 98 spn_all = spn_all + spn_k*kweight 99 end do 100 101 end if 102 103 ! Collect contributions from all nodes 104 ! 105 call comms_reduce(spn_all(1), 3, 'SUM') 106 107 ! No factor of g=2 because the spin variable spans [-1,1], not 108 ! [-1/2,1/2] (i.e., it is really the Pauli matrix sigma, not S) 109 ! 110 spn_mom(1:3) = -spn_all(1:3) 111 112 if (on_root) then 113 write (stdout, '(/,1x,a)') 'Spin magnetic moment (Bohr magn./cell)' 114 write (stdout, '(1x,a,/)') '====================' 115 write (stdout, '(1x,a18,f11.6)') 'x component:', spn_mom(1) 116 write (stdout, '(1x,a18,f11.6)') 'y component:', spn_mom(2) 117 write (stdout, '(1x,a18,f11.6)') 'z component:', spn_mom(3) 118 119 ! Polar and azimuthal angles of the magnetization (defined as in pwscf) 120 ! 121 conv = 180.0_dp/pi 122 magnitude = sqrt(spn_mom(1)**2 + spn_mom(2)**2 + spn_mom(3)**2) 123 theta = acos(spn_mom(3)/magnitude)*conv 124 phi = atan(spn_mom(2)/spn_mom(1))*conv 125 write (stdout, '(/,1x,a18,f11.6)') 'Polar theta (deg):', theta 126 write (stdout, '(1x,a18,f11.6)') 'Azim. phi (deg):', phi 127 end if 128 129 end subroutine spin_get_moment 130 131! ========================================================================= 132 133 subroutine spin_get_nk(kpt, spn_nk) 134 !=============================================================! 135 ! ! 136 !! Computes <psi_{mk}^(H)|S.n|psi_{mk}^(H)> (m=1,...,num_wann) 137 !! where S.n = n_x.S_x + n_y.S_y + n_z.Z_z 138 !! 139 !! S_i are the Pauli matrices and n=(n_x,n_y,n_z) is the unit 140 !! vector along the chosen spin quantization axis 141 ! ! 142 !============================================================ ! 143 144 use w90_constants, only: dp, pi, cmplx_0, cmplx_i 145 use w90_io, only: io_error 146 use w90_utility, only: utility_diagonalize, utility_rotate_diag 147 use w90_parameters, only: num_wann, spin_axis_polar, & 148 spin_axis_azimuth 149 use w90_postw90_common, only: pw90common_fourier_R_to_k 150 use w90_get_oper, only: HH_R, SS_R 151 152 ! Arguments 153 ! 154 real(kind=dp), intent(in) :: kpt(3) 155 real(kind=dp), intent(out) :: spn_nk(num_wann) 156 157 ! Physics 158 ! 159 complex(kind=dp), allocatable :: HH(:, :) 160 complex(kind=dp), allocatable :: UU(:, :) 161 complex(kind=dp), allocatable :: SS(:, :, :), SS_n(:, :) 162 163 ! Misc/Dummy 164 ! 165 integer :: is 166 real(kind=dp) :: eig(num_wann), alpha(3), conv 167 168 allocate (HH(num_wann, num_wann)) 169 allocate (UU(num_wann, num_wann)) 170 allocate (SS(num_wann, num_wann, 3)) 171 allocate (SS_n(num_wann, num_wann)) 172 173 call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0) 174 call utility_diagonalize(HH, num_wann, eig, UU) 175 176 do is = 1, 3 177 call pw90common_fourier_R_to_k(kpt, SS_R(:, :, :, is), SS(:, :, is), 0) 178 enddo 179 180 ! Unit vector along the magnetization direction 181 ! 182 conv = 180.0_dp/pi 183 alpha(1) = sin(spin_axis_polar/conv)*cos(spin_axis_azimuth/conv) 184 alpha(2) = sin(spin_axis_polar/conv)*sin(spin_axis_azimuth/conv) 185 alpha(3) = cos(spin_axis_polar/conv) 186 187 ! Vector of spin matrices projected along the quantization axis 188 ! 189 SS_n(:, :) = alpha(1)*SS(:, :, 1) + alpha(2)*SS(:, :, 2) + alpha(3)*SS(:, :, 3) 190 191 spn_nk(:) = real(utility_rotate_diag(SS_n, UU, num_wann), dp) 192 193 end subroutine spin_get_nk 194 195 !===========================================================! 196 ! PRIVATE PROCEDURES ! 197 !===========================================================! 198 199 subroutine spin_get_moment_k(kpt, ef, spn_k) 200 !! Computes the spin magnetic moment by Wannier interpolation 201 !! at the specified k-point 202 use w90_constants, only: dp, cmplx_0, cmplx_i 203 use w90_io, only: io_error 204 use w90_utility, only: utility_diagonalize, utility_rotate_diag 205 use w90_parameters, only: num_wann 206 use w90_postw90_common, only: pw90common_fourier_R_to_k, pw90common_get_occ 207 use w90_get_oper, only: HH_R, SS_R 208 ! Arguments 209 ! 210 real(kind=dp), intent(in) :: kpt(3) 211 real(kind=dp), intent(in) :: ef 212 real(kind=dp), intent(out) :: spn_k(3) 213 214 ! Physics 215 ! 216 complex(kind=dp), allocatable :: HH(:, :) 217 complex(kind=dp), allocatable :: SS(:, :, :) 218 complex(kind=dp), allocatable :: UU(:, :) 219 real(kind=dp) :: spn_nk(num_wann, 3) 220 221 ! Misc/Dummy 222 ! 223 integer :: i, is 224 real(kind=dp) :: eig(num_wann), occ(num_wann) 225 226 allocate (HH(num_wann, num_wann)) 227 allocate (UU(num_wann, num_wann)) 228 allocate (SS(num_wann, num_wann, 3)) 229 230 call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0) 231 call utility_diagonalize(HH, num_wann, eig, UU) 232 call pw90common_get_occ(eig, occ, ef) 233 234 spn_k(1:3) = 0.0_dp 235 do is = 1, 3 236 call pw90common_fourier_R_to_k(kpt, SS_R(:, :, :, is), SS(:, :, is), 0) 237 spn_nk(:, is) = aimag(cmplx_i*utility_rotate_diag(SS(:, :, is), UU, num_wann)) 238 do i = 1, num_wann 239 spn_k(is) = spn_k(is) + occ(i)*spn_nk(i, is) 240 end do 241 enddo 242 243 end subroutine spin_get_moment_k 244 245 subroutine spin_get_S(kpt, S) 246 !===========================================================! 247 ! ! 248 ! Computes <psi_{nk}^(H)|S|psi_{nk}^(H)> (n=1,...,num_wann) ! 249 ! where S = (S_x,S_y,S_z) is the vector of Pauli matrices ! 250 ! ! 251 !========================================================== ! 252 253 use w90_constants, only: dp, pi, cmplx_0, cmplx_i 254 use w90_io, only: io_error 255 use w90_utility, only: utility_diagonalize, utility_rotate_diag 256 use w90_parameters, only: num_wann 257 use w90_postw90_common, only: pw90common_fourier_R_to_k 258 use w90_get_oper, only: HH_R, SS_R 259 260 ! Arguments 261 ! 262 real(kind=dp), intent(in) :: kpt(3) 263 real(kind=dp), intent(out) :: S(num_wann, 3) 264 265 ! Physics 266 ! 267 complex(kind=dp), allocatable :: HH(:, :) 268 complex(kind=dp), allocatable :: UU(:, :) 269 complex(kind=dp), allocatable :: SS(:, :, :) 270 real(kind=dp) :: eig(num_wann) 271 272 ! Misc/Dummy 273 ! 274 integer :: i 275 276 allocate (HH(num_wann, num_wann)) 277 allocate (UU(num_wann, num_wann)) 278 allocate (SS(num_wann, num_wann, 3)) 279 280 call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0) 281 call utility_diagonalize(HH, num_wann, eig, UU) 282 283 do i = 1, 3 284 call pw90common_fourier_R_to_k(kpt, SS_R(:, :, :, i), SS(:, :, i), 0) 285 S(:, i) = real(utility_rotate_diag(SS(:, :, i), UU, num_wann), dp) 286 enddo 287 288 end subroutine spin_get_S 289 290end module w90_spin 291