1 ! 2 ! Copyright (C) 2010-2019 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino 3 ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino 4 ! 5 ! This file is distributed under the terms of the GNU General Public 6 ! License. See the file `LICENSE' in the root directory of the 7 ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . 8 ! 9 !---------------------------------------------------------------------- 10 MODULE wigner 11 !---------------------------------------------------------------------- 12 !! 13 !! This module contains all the routines linked creation of Wigner-Seitz cell 14 !! 15 IMPLICIT NONE 16 ! 17 CONTAINS 18 !----------------------------------------------------------------- 19 SUBROUTINE wigner_seitz_wrap(nkc1, nkc2, nkc3, nqc1, nqc2, nqc3, & 20 irvec_k, irvec_q, irvec_g, & 21 ndegen_k, ndegen_q, ndegen_g, & 22 wslen_k, wslen_q, wslen_g, & 23 w_centers, dims, tau, dims2) 24 !----------------------------------------------------------------- 25 !! 26 !! June 2018 - SP - CV 27 !! 28 !! This routine wrap the call to three Wigner-Seitz routines: 29 !! wigner_seitzk : Contruct a grid of points that fall inside of (and eventually 30 !! on the surface of) the Wigner-Seitz supercell centered on the 31 !! origin of the Bravais lattice. Use for electronic properties. 32 !! wigner_seitzq : Creates a set of WS vectors for each pair of atoms tau(nb)-tau(na) 33 !! On exiting, ndegen_q contains the degeneracies of each pairs 34 !! of atoms while irvec_q contains the minimal communal sets of WS vectors. 35 !! Used for phonon properties 36 !! wigner_seitzg : Creates a set of WS vector for each atoms tau(na). 37 !! On exiting, ndegen_g contains the degeneracies of each atoms 38 !! while irvec_g contains the minimal communal sets of WS vectors. 39 !! Used for electron-phonon properties. 40 !! 41 !! Note 1: ndegen_k is always > 0 while ndegen_q and ndegen_g might contains 0 weigths. 42 !! Note 2: No sorting of vectors is needed anymore 43 !! Note 3: The dimension 20*nkc1*nkc2*nkc3 should be safe enough. 44 !! Note 4: The Wigner-Seitz construction in EPW was done by constructing a cell 45 !! centred unit cell. This is fine for electronic properties (this is what is done in wannier90). 46 !! However for phonon or electron-phonon properties, one can have issues when the cell 47 !! is tilted for example. 48 !! The proper way is to construct a set of WS vectors centred on pairs of atoms (phonons) 49 !! or atoms (el-ph). 50 !! In the matdyn code, a FT grid is constructed with weights centred on pairs of atoms 51 !! and zeros everywhere else. 52 !! EPW now reproduced exactly the results of matdyn for the interpolated phonons at a 53 !! lower computation cost. Indeed we minimize the number of zeros by keeping the union 54 !! of values between all the cells. 55 !! In both cases this is very fast anyway but is important for el-ph properties. 56 !! 57 !----------------------------------------------------------------- 58 USE kinds, ONLY : DP 59 ! 60 IMPLICIT NONE 61 ! 62 INTEGER, INTENT(in) :: nkc1 63 !! size of the uniform k mesh 64 INTEGER, INTENT(in) :: nkc2 65 !! size of the uniform k mesh 66 INTEGER, INTENT(in) :: nkc3 67 !! size of the uniform k mesh 68 INTEGER, INTENT(in) :: nqc1 69 !! size of the uniform q mesh 70 INTEGER, INTENT(in) :: nqc2 71 !! size of the uniform q mesh 72 INTEGER, INTENT(in) :: nqc3 73 !! size of the uniform q mesh 74 INTEGER, INTENT(in) :: dims 75 !! Number of bands in the Wannier space 76 INTEGER, INTENT(in) :: dims2 77 !! Number of atoms 78 INTEGER, ALLOCATABLE, INTENT(out) :: irvec_k(:, :) 79 !! INTEGER components of the ir-th Wigner-Seitz grid point in the basis 80 !! of the lattice vectors for electrons 81 INTEGER, ALLOCATABLE, INTENT(out) :: irvec_q(:, :) 82 !! INTEGER components of the ir-th Wigner-Seitz grid point for phonons 83 INTEGER, ALLOCATABLE, INTENT(out) :: irvec_g(:, :) 84 !! INTEGER components of the ir-th Wigner-Seitz grid point for electron-phonon 85 INTEGER, ALLOCATABLE, INTENT(out) :: ndegen_k(:, :, :) 86 !! Wigner-Seitz number of degenerescence (weights) for the electrons grid 87 INTEGER, ALLOCATABLE, INTENT(out) :: ndegen_q(:, :, :) 88 !! Wigner-Seitz weights for the phonon grid that depend on 89 !! atomic positions $R + \tau(nb) - \tau(na)$ 90 INTEGER, ALLOCATABLE, INTENT(out) :: ndegen_g(:, :, :, :) 91 !! Wigner-Seitz weights for the electron-phonon grid that depend on 92 !! atomic positions $R - \tau(na)$ 93 REAL(KIND = DP), ALLOCATABLE, INTENT(out) :: wslen_k(:) 94 !! real-space length for electrons, in units of alat 95 REAL(KIND = DP), ALLOCATABLE, INTENT(out) :: wslen_q(:) 96 !! real-space length for phonons, in units of alat 97 REAL(KIND = DP), ALLOCATABLE, INTENT(out) :: wslen_g(:) 98 !! real-space length for electron-phonons, in units of alat 99 REAL(KIND = DP), INTENT(in) :: w_centers(3, dims) 100 !! Wannier centers used for the creation of electron and el-ph WS cells 101 REAL(KIND = DP), INTENT(in) :: tau(3, dims2) 102 !! Atomic position in the unit cell. 103 ! 104 ! Work Variables 105 INTEGER :: ir 106 !! Index for WS vectors 107 INTEGER :: nrr_k 108 !! maximum number of WS vectors for the electrons 109 INTEGER :: nrr_q 110 !! maximum number of WS vectors for the phonons 111 INTEGER :: nrr_g 112 !! maximum number of WS vectors for the electron-phonon 113 INTEGER :: ierr 114 !! Error status 115 INTEGER :: irvec_kk (3, 20 * nkc1 * nkc2 * nkc3) 116 !! local INTEGER components of the ir-th Wigner-Seitz grid point 117 !! in the basis of the lattice vectors for electrons 118 INTEGER :: irvec_qq(3, 20 * nqc1 * nqc2 * nqc3) 119 !! local INTEGER components of the ir-th Wigner-Seitz grid point for phonons 120 INTEGER :: irvec_gg(3, 20 * nqc1 * nqc2 * nqc3) 121 !! local INTEGER components of the ir-th Wigner-Seitz grid point for electron-phonons 122 !! We use nkc1 instead of nqc1 because the k-grid is always larger or equal to q-grid. 123 INTEGER :: ndegen_kk(20 * nkc1 * nkc2 * nkc3, dims, dims) 124 !! local Wigner-Seitz number of degenerescence (weights) for the electrons grid 125 INTEGER :: ndegen_qq(20 * nqc1 * nqc2 * nqc3, dims2, dims2) 126 !! local Wigner-Seitz number of degenerescence (weights) for the phonons grid 127 INTEGER :: ndegen_gg(20 * nqc1 * nqc2 * nqc3, dims2, dims, dims) 128 !! local Wigner-Seitz number of degenerescence (weights) for the electron-phonons grid 129 REAL(KIND = DP) :: wslen_kk(20 * nkc1 * nkc2 * nkc3) 130 !! local real-space length for electrons, in units of alat 131 REAL(KIND = DP) :: wslen_qq(20 * nqc1 * nqc2 * nqc3) 132 !! local real-space length for phonons, in units of alat 133 REAL(KIND = DP) :: wslen_gg(20 * nqc1 * nqc2 * nqc3) 134 !! local real-space length for electron-phonon, in units of alat 135 ! 136 ! Check the bounds 137 IF (nqc1 > nkc1 .OR. nqc2 > nkc2 .OR. nqc3 > nkc3 ) call errore & 138 ('wigner_seitz_wrap', ' the phonon grid should be smaller than electron grid', 1) 139 ! 140 ! Now generated the un-sorted points for the electrons, phonons and electron-phonon 141 ! 142 ! If dims > 1, it includes the position of Wannier-Centers 143 CALL wigner_seitzkq(nkc1, nkc2, nkc3, irvec_kk, ndegen_kk, wslen_kk, nrr_k, w_centers, dims) 144 CALL wigner_seitzkq(nqc1, nqc2, nqc3, irvec_qq, ndegen_qq, wslen_qq, nrr_q, tau, dims2) 145 CALL wigner_seitzg(nqc1, nqc2, nqc3, irvec_gg, ndegen_gg, wslen_gg, nrr_g, w_centers, tau, dims, dims2) 146 ! 147 ALLOCATE(irvec_k(3, nrr_k), STAT = ierr) 148 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_k', 1) 149 ALLOCATE(irvec_q(3, nrr_q), STAT = ierr) 150 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_q', 1) 151 ALLOCATE(irvec_g(3, nrr_g), STAT = ierr) 152 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating irvec_g', 1) 153 ALLOCATE(ndegen_k(nrr_k, dims, dims), STAT = ierr) 154 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_k', 1) 155 ALLOCATE(ndegen_q(nrr_q, dims2, dims2), STAT = ierr) 156 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_q', 1) 157 ALLOCATE(ndegen_g(nrr_g, dims2, dims, dims), STAT = ierr) 158 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating ndegen_g', 1) 159 ALLOCATE(wslen_k(nrr_k), STAT = ierr) 160 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_k', 1) 161 ALLOCATE(wslen_q(nrr_q), STAT = ierr) 162 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_q', 1) 163 ALLOCATE(wslen_g(nrr_g), STAT = ierr) 164 IF (ierr /= 0) CALL errore('wigner_seitz_wrap', 'Error allocating wslen_g', 1) 165 ! 166 ! Create vectors with correct size. 167 DO ir = 1, nrr_k 168 ndegen_k(ir, :, :) = ndegen_kk(ir, :, :) 169 irvec_k(:, ir) = irvec_kk(:, ir) 170 wslen_k(ir) = wslen_kk(ir) 171 ENDDO 172 DO ir = 1, nrr_q 173 ndegen_q(ir, :, :) = ndegen_qq(ir, :, :) 174 irvec_q(:, ir) = irvec_qq(:, ir) 175 wslen_q(ir) = wslen_qq(ir) 176 ENDDO 177 DO ir = 1, nrr_g 178 ndegen_g(ir, :, :, :) = ndegen_gg(ir, :, :, :) 179 irvec_g(:, ir) = irvec_gg(:, ir) 180 wslen_g(ir) = wslen_gg(ir) 181 ENDDO 182 ! 183 !----------------------------------------------------------------------------- 184 END SUBROUTINE wigner_seitz_wrap 185 !----------------------------------------------------------------------------- 186 ! 187 !----------------------------------------------------------------------------- 188 SUBROUTINE wigner_seitzkq(nc1, nc2, nc3, irvec, ndegen, wslen, nrr, shift, dims) 189 !----------------------------------------------------------------------------- 190 !! 191 !! Calculates a grid of points that fall inside of (and eventually 192 !! on the surface of) the Wigner-Seitz supercell centered on the 193 !! origin of the Bravais lattice with primitive translations 194 !! nkc1*a_1+nkc2*a_2+nkc3*a_3 195 !! 196 !----------------------------------------------------------------- 197 USE kinds, ONLY : DP 198 USE cell_base, ONLY : at, bg 199 USE constants_epw, ONLY : eps6 200 USE low_lvl, ONLY : hpsort_eps_epw 201 ! 202 IMPLICIT NONE 203 ! 204 INTEGER, INTENT(in) :: nc1 205 !! size of the uniform k mesh 206 INTEGER, INTENT(in) :: nc2 207 !! size of the uniform k mesh 208 INTEGER, INTENT(in) :: nc3 209 !! size of the uniform k mesh 210 INTEGER, INTENT(in) :: dims 211 !! dim of shift(3,dims) 212 INTEGER, INTENT(out) :: irvec(3, 20 * nc1 * nc2 * nc3) 213 !! INTEGER components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors 214 INTEGER, INTENT(out) :: ndegen(20 * nc1 * nc2 * nc3, dims, dims) 215 !! Number of degeneracies 216 INTEGER, INTENT(out) :: nrr 217 !! number of Wigner-Seitz grid points 218 REAL(KIND = DP), INTENT(out) :: wslen(20 * nc1 * nc2 * nc3) 219 !! real-space length, in units of alat 220 REAL(KIND = DP), INTENT(in) :: shift(3, dims) 221 !! Wannier centers (electron WS cell) or atomic positions (lattice WS cell) 222 ! 223 ! Local variables 224 LOGICAL :: found 225 !! True if the vector has been found 226 INTEGER :: n1, n2, n3 227 !! Index for the larger supercell 228 INTEGER :: i1, i2, i3 229 !! Index to compute |r-R| distance 230 INTEGER :: i, ir, irtot, iw, iw2 231 !! Iterative index 232 INTEGER :: ipol, jpol 233 !! Cartesian direction 234 INTEGER, ALLOCATABLE :: ind(:) 235 !! Index of sorting 236 INTEGER :: nind 237 !! The metric tensor 238 INTEGER :: ierr 239 !! Error status 240 INTEGER :: nrr_tmp(dims, dims) 241 !! Temporary WS matrix 242 INTEGER :: irvec_tmp(3, 20 * nc1 * nc2 * nc3, dims, dims) 243 !! Temp 244 INTEGER :: ndegen_tmp(20 * nc1 * nc2 * nc3, dims, dims) 245 !! Number of degenerate points 246 REAL(KIND = DP) :: adot(3, 3) 247 !! Dot product between lattice vector 248 REAL(KIND = DP) :: dist(125) 249 !! Contains the distance squared |r-R|^2 250 REAL(KIND = DP) :: mindist 251 !! Minimum distance 252 REAL(KIND = DP) :: tot, tot2 253 !! Sum of all the weigths 254 REAL(KIND = DP) :: ndiff(3) 255 !! Distances. Must be now real because of fractional atoms 256 ! 257 nind = 20 * nc1 * nc2 * nc3 258 IF (nind < 125) THEN 259 nind = 125 260 ENDIF 261 ALLOCATE(ind(nind), STAT = ierr) 262 IF (ierr /= 0) CALL errore('wigner_seitzkq', 'Error allocating ind', 1) 263 ! 264 DO ipol = 1, 3 265 DO jpol = 1, 3 266 adot(ipol, jpol) = DOT_PRODUCT(at(:, ipol), at(:, jpol)) 267 ENDDO 268 ENDDO 269 ! 270 CALL cryst_to_cart(dims, shift(:, :), bg, -1) 271 ! 272 ndegen_tmp(:, :, :) = 0 273 irvec_tmp(:, :, :, :) = 0 274 DO iw = 1, dims 275 DO iw2 = 1, dims 276 nrr_tmp(iw, iw2) = 0 277 DO n1 = -2 * nc1, 2 * nc1 278 DO n2 = -2 * nc2, 2 * nc2 279 DO n3 = -2 * nc3, 2 * nc3 280 ! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63 281 i = 0 282 dist(:) = 0.d0 283 DO i1 = -2, 2 284 DO i2 = -2, 2 285 DO i3 = -2, 2 286 i = i + 1 287 ! Calculate distance squared |r-R|^2 288 ndiff(1) = n1 - i1 * nc1 + shift(1, iw2) - shift(1, iw) 289 ndiff(2) = n2 - i2 * nc2 + shift(2, iw2) - shift(2, iw) 290 ndiff(3) = n3 - i3 * nc3 + shift(3, iw2) - shift(3, iw) 291 !ndiff(1) = n1 - i1*nc1 + (shift(1,iw2) + shift(1,iw)) / 2.d0 292 !ndiff(2) = n2 - i2*nc2 + (shift(2,iw2) + shift(2,iw)) / 2.d0 293 !ndiff(3) = n3 - i3*nc3 + (shift(3,iw2) + shift(3,iw)) / 2.d0 294 DO ipol = 1, 3 295 DO jpol = 1, 3 296 dist(i) = dist(i) + DBLE(ndiff(ipol)) * adot(ipol, jpol) * DBLE(ndiff(jpol)) 297 ENDDO 298 ENDDO 299 ! 300 ENDDO ! i3 301 ENDDO ! i2 302 ENDDO ! i1 303 !IF(iw==iw2) print*,'iw dist ',iw,dist(63) 304 ! 305 ! Sort the 125 vectors R by increasing value of |r-R|^2 306 ind(1) = 0 ! required for hpsort_eps (see the subroutine) 307 CALL hpsort_eps_epw(125, dist, ind, eps6) 308 ! 309 ! Find all the vectors R with the (same) smallest |r-R|^2; 310 ! if R=0 is one of them, then the current point r belongs to 311 ! Wignez-Seitz cell => set found to true 312 ! 313 found = .FALSE. 314 i = 1 315 mindist = dist(1) 316 DO WHILE (ABS(dist(i) - mindist) < eps6 .AND. i < 125) 317 IF (ind(i) == 63) found = .TRUE. 318 i = i + 1 319 ENDDO 320 ! 321 IF (found) THEN 322 nrr_tmp(iw, iw2) = nrr_tmp(iw, iw2) + 1 323 ndegen_tmp(nrr_tmp(iw, iw2), iw, iw2) = i - 1 324 irvec_tmp(:, nrr_tmp(iw, iw2), iw, iw2) = (/n1, n2, n3/) 325 ENDIF 326 ENDDO ! n3 327 ENDDO ! n2 328 ENDDO ! n1 329 ENDDO ! iw2 330 ENDDO ! iw 331 ! 332 nrr = nrr_tmp(1, 1) 333 irvec(:, :) = irvec_tmp(:, :, 1, 1) 334 DO iw = 1, dims 335 DO iw2 = 1, dims 336 DO ir = 1, nrr_tmp(iw, iw2) 337 found = .FALSE. 338 DO irtot = 1, nrr 339 IF (ALL(irvec_tmp(:, ir, iw, iw2) == irvec(:, irtot))) THEN 340 found = .TRUE. 341 ENDIF 342 ENDDO !nrr 343 IF (.NOT. found) THEN 344 nrr = nrr + 1 345 irvec(:, nrr) = irvec_tmp(:, ir, iw, iw2) 346 ENDIF 347 ENDDO ! ir 348 ENDDO ! nb 349 ENDDO ! na 350 ! 351 ! Creates a pair of atoms depended degeneracy array but with a number of WS 352 ! vectors per pair that is equal to the global set. Populate with zero weights 353 ! the one that are not part of that pair set. 354 ndegen(:, :, :) = 0 355 DO iw = 1, dims 356 DO iw2 = 1, dims 357 DO ir = 1, nrr_tmp(iw, iw2) 358 DO irtot = 1, nrr 359 IF (ALL(irvec(:, irtot) == irvec_tmp(:, ir, iw, iw2))) THEN 360 ndegen(irtot, iw, iw2) = ndegen_tmp(ir, iw, iw2) 361 ENDIF 362 ENDDO 363 ENDDO 364 ENDDO 365 ENDDO 366 ! 367 DO iw = 1, dims 368 DO iw2 = 1, dims 369 tot = 0.d0 370 tot2 = 0.d0 371 DO i = 1, nrr 372 IF (ndegen(i, iw, iw2) > 0) THEN 373 tot2 = tot2 + 1.d0 / DBLE(ndegen(i, iw, iw2)) 374 ENDIF 375 ENDDO 376 DO i = 1, nrr_tmp(iw, iw2) 377 tot = tot + 1.d0 / DBLE(ndegen_tmp(i, iw, iw2)) 378 ENDDO 379 ! 380 IF (ABS(tot - DBLE(nc1*nc2*nc3)) > eps6) CALL errore & 381 ('wigner_seitzkq',' weights do not add up to nc1*nc2*nc3', 1) 382 IF (ABS(tot - tot2) > eps6) CALL errore & 383 ('wigner_seitzkq', ' weigths of pair of atoms is not equal to global weights', 1) 384 ENDDO 385 ENDDO 386 ! 387 IF (nrr > 20 * nc1 * nc2 * nc3) CALL errore & 388 ('wigner_seitzkq', 'too many WS points, try to increase the bound 20*nc1*nc2*nc3', 1) 389 ! 390 wslen(:) = 0.d0 391 DO i = 1, nrr 392 DO ipol = 1, 3 393 DO jpol = 1, 3 394 wslen(i) = wslen(i) + DBLE(irvec(ipol, i)) * adot(ipol, jpol) * DBLE(irvec(jpol, i)) 395 ENDDO 396 ENDDO 397 wslen(i) = DSQRT(wslen(i)) 398 ENDDO 399 ! 400 CALL cryst_to_cart(dims, shift(:, :), at, 1) 401 ! 402 DEALLOCATE(ind, STAT = ierr) 403 IF (ierr /= 0) CALL errore('wigner_seitzkq', 'Error deallocating ind', 1) 404 ! 405 !----------------------------------------------------------------------------- 406 END SUBROUTINE wigner_seitzkq 407 !----------------------------------------------------------------------------- 408 ! 409 !----------------------------------------------------------------- 410 SUBROUTINE wigner_seitzg(nc1, nc2, nc3, irvec, ndegen, wslen, nrr, w_centers, tau, dims, dims2) 411 !----------------------------------------------------------------- 412 !! 413 !! Calculates a grid of points that fall inside of (and eventually 414 !! on the surface of) the Wigner-Seitz supercell centered on 415 !! each atoms. 416 !! Follows Eq. 66 of PRB 55, 10355 (1997). 417 !! We are part of the WS if $R_b - \tau_{\kappa}$ is inside the 418 !! supercell. 419 !! 420 !----------------------------------------------------------------- 421 USE kinds, ONLY : DP 422 USE cell_base, ONLY : at, bg 423 USE constants_epw, ONLY : eps6 424 USE low_lvl, ONLY : hpsort_eps_epw 425 ! 426 IMPLICIT NONE 427 ! 428 INTEGER, INTENT(in) :: nc1 429 !! size of the uniform k mesh 430 INTEGER, INTENT(in) :: nc2 431 !! size of the uniform k mesh 432 INTEGER, INTENT(in) :: nc3 433 !! size of the uniform k mesh 434 INTEGER, INTENT(in) :: dims 435 !! Dims is either nbndsub or 1 depending on use_ws 436 INTEGER, INTENT(in) :: dims2 437 !! Number of atoms 438 INTEGER, INTENT(out) :: irvec(3, 20 * nc1 * nc2 * nc3) 439 !! INTEGER components of the ir-th Wigner-Seitz grid point in the basis of the lattice vectors 440 INTEGER, INTENT(out) :: ndegen(20 * nc1 * nc2 * nc3, dims2, dims, dims) 441 !! Number of degeneracies 442 INTEGER, INTENT(out) :: nrr 443 !! number of Wigner-Seitz grid points 444 REAL(KIND = DP), INTENT(in) :: w_centers(3, dims) 445 !! Wannier centers 446 REAL(KIND = DP), INTENT(in) :: tau(3, dims2) 447 !! Atomic positions 448 REAL(KIND = DP), INTENT(out) :: wslen(20 * nc1 * nc2 * nc3) 449 !! real-space length, in units of alat 450 ! 451 ! Local variables 452 LOGICAL :: found 453 !! True if the vector has been foun 454 INTEGER :: n1, n2, n3 455 !! Index for the larger supercell 456 INTEGER :: i1, i2, i3 457 !! Index to compute |r-R| distance 458 INTEGER :: na, iw, iw2 459 !! Atom index 460 INTEGER :: i 461 !! Iterative index 462 INTEGER :: ir 463 !! Iterative index on the pair of atoms 464 INTEGER :: irtot 465 !! Iterative index on the combined pair of atoms 466 INTEGER :: ipol, jpol 467 !! Cartesian direction 468 INTEGER, ALLOCATABLE :: ind(:) 469 !! Index of sorting 470 INTEGER :: nind 471 !! The metric tensor 472 INTEGER :: ierr 473 !! Error status 474 INTEGER :: nrr_tmp(dims2, dims, dims) 475 !! Temporary array that contains the max number of WS vectors 476 !! for a pair of atoms. 477 INTEGER :: irvec_tmp(3, 20 * nc1 * nc2 * nc3, dims2, dims, dims) 478 !! Temporary WS vectors for each atoms 479 INTEGER :: ndegen_tmp(20 * nc1 * nc2 * nc3, dims2, dims, dims) 480 !! Temporary WS vectors weigths for each atoms 481 REAL(KIND = DP) :: adot(3, 3) 482 !! Dot product between lattice vector 483 REAL(KIND = DP) :: dist(125) 484 !! Contains the distance squared |r-R|^2 485 REAL(KIND = DP) :: mindist 486 !! Minimum distance 487 REAL(KIND = DP) :: tot, tot2 488 !! Sum of all the weigths 489 REAL(KIND = DP) :: ndiff(3) 490 !! Distances. Must be now real because of fractional atoms 491 ! 492 nind = 20 * nc1 * nc2 * nc3 493 IF (nind < 125) THEN 494 nind = 125 495 ENDIF 496 ALLOCATE(ind(nind), STAT = ierr) 497 IF (ierr /= 0) CALL errore('wigner_seitzg', 'Error allocating ind', 1) 498 ! 499 DO ipol = 1, 3 500 DO jpol = 1, 3 501 adot(ipol, jpol) = DOT_PRODUCT(at(:, ipol), at(:, jpol)) 502 ENDDO 503 ENDDO 504 ! 505 CALL cryst_to_cart(dims2, tau(:, :), bg, -1) 506 CALL cryst_to_cart(dims, w_centers(:, :), bg, -1) 507 ! 508 ! Loop over grid points r on a unit cell that is 8 times larger than a 509 ! primitive supercell. In the end nrr contains the total number of grids 510 ! points that have been found in the Wigner-Seitz cell 511 ! 512 nrr_tmp(:, :, :) = 0 513 DO iw = 1, dims 514 DO iw2 = 1, dims 515 DO na = 1, dims2 516 DO n1 = -2 * nc1, 2 * nc1 517 DO n2 = -2 * nc2, 2 * nc2 518 DO n3 = -2 * nc3, 2 * nc3 519 ! 520 ! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63 521 ! 522 i = 0 523 dist(:) = 0.d0 524 DO i1 = -2, 2 525 DO i2 = -2, 2 526 DO i3 = -2, 2 527 i = i + 1 528 ! 529 ! Calculate distance squared |r-R|^2 530 ! 531 !ndiff(1) = n1 - i1*nc1 + ( tau(1,na) + w_centers(1,iw2) + w_centers(1,iw) ) / 3.d0 532 !ndiff(2) = n2 - i2*nc2 + ( tau(2,na) + w_centers(2,iw2) + w_centers(2,iw) ) / 3.d0 533 !ndiff(3) = n3 - i3*nc3 + ( tau(3,na) + w_centers(3,iw2) + w_centers(3,iw) ) / 3.d0 534 ndiff(1) = n1 - i1 * nc1 + tau(1, na) - (w_centers(1, iw2) + w_centers(1, iw)) / 2.d0 535 ndiff(2) = n2 - i2 * nc2 + tau(2, na) - (w_centers(2, iw2) + w_centers(2, iw)) / 2.d0 536 ndiff(3) = n3 - i3 * nc3 + tau(3, na) - (w_centers(3, iw2) + w_centers(3, iw)) / 2.d0 537 DO ipol = 1, 3 538 DO jpol = 1, 3 539 dist(i) = dist(i) + DBLE(ndiff(ipol)) * adot(ipol, jpol) * DBLE(ndiff(jpol)) 540 ENDDO 541 ENDDO 542 ! 543 ENDDO 544 ENDDO 545 ENDDO 546 ! 547 ! Sort the 125 vectors R by increasing value of |r-R|^2 548 ind(1) = 0 ! required for hpsort_eps (see the subroutine) 549 CALL hpsort_eps_epw(125, dist, ind, eps6) 550 ! 551 ! Find all the vectors R with the (same) smallest |r-R|^2; 552 ! if R=0 is one of them, then the current point r belongs to 553 ! Wignez-Seitz cell => set found to true 554 ! 555 found = .FALSE. 556 i = 1 557 mindist = dist(1) 558 DO WHILE (ABS(dist(i) - mindist) < eps6 .AND. i < 125) 559 IF (ind(i) == 63) found = .TRUE. 560 i = i + 1 561 ENDDO 562 ! 563 IF (found) THEN 564 nrr_tmp(na, iw, iw2) = nrr_tmp(na, iw, iw2) + 1 565 ndegen_tmp(nrr_tmp(na, iw, iw2), na, iw, iw2) = i - 1 566 irvec_tmp (:, nrr_tmp(na, iw, iw2), na, iw, iw2) = (/n1, n2, n3/) 567 ENDIF 568 ENDDO ! n3 569 ENDDO ! n2 570 ENDDO ! n3 571 ENDDO ! na 572 ENDDO ! iw2 573 ENDDO ! iw 574 ! 575 ! Now creates a global set of WS vectors from all the atoms pair. 576 ! Also remove the duplicated ones. 577 nrr = nrr_tmp(1, 1, 1) 578 irvec(:, :) = irvec_tmp(:, :, 1, 1, 1) 579 DO iw = 1, dims 580 DO iw2 = 1, dims 581 DO na = 1, dims2 582 DO ir = 1, nrr_tmp(na, iw, iw2) 583 found = .FALSE. 584 DO irtot = 1, nrr 585 IF (ALL(irvec_tmp(:, ir, na, iw, iw2) == irvec(:, irtot))) THEN 586 found = .TRUE. 587 ENDIF 588 ENDDO !nrr 589 IF (.NOT. found) THEN 590 nrr = nrr + 1 591 irvec(:, nrr) = irvec_tmp(:, ir, na, iw, iw2) 592 ENDIF 593 ENDDO ! ir 594 ENDDO ! na 595 ENDDO ! iw2 596 ENDDO ! iw 597 ! 598 ! Creates a pair of atoms-dependent degeneracy array but with a number of WS 599 ! vectors per pair that is equal to the global set. Populate with zero weights 600 ! the one that are not part of that pair set. 601 ndegen(:, :, :, :) = 0 602 DO iw = 1, dims 603 DO iw2 = 1, dims 604 DO na = 1, dims2 605 DO ir = 1, nrr_tmp(na, iw, iw2) 606 DO irtot = 1, nrr 607 IF (ALL(irvec(:, irtot) == irvec_tmp(:, ir, na, iw, iw2))) THEN 608 ndegen(irtot, na, iw, iw2) = ndegen_tmp(ir, na, iw, iw2) 609 ENDIF 610 ENDDO 611 ENDDO 612 ENDDO 613 ENDDO 614 ENDDO 615 ! 616 DO iw = 1, dims 617 DO iw2 = 1, dims 618 DO na = 1, dims2 619 tot = 0.d0 620 tot2 = 0.d0 621 DO i = 1, nrr 622 IF (ndegen(i, na, iw, iw2) > 0) THEN 623 tot2 = tot2 + 1.d0 / DBLE(ndegen(i, na, iw, iw2)) 624 ENDIF 625 ENDDO 626 DO i = 1, nrr_tmp(na, iw, iw2) 627 tot = tot + 1.d0 / DBLE(ndegen_tmp(i, na, iw, iw2)) 628 ENDDO 629 ! 630 IF (ABS(tot - DBLE(nc1 * nc2 * nc3)) > eps6) CALL errore & 631 ('wigner_seitzg', ' weights do not add up to nqc1*nqc2*nqc3', 1) 632 IF (ABS(tot - tot2) > eps6) CALL errore & 633 ('wigner_seitzg', ' weigths of pair of atoms is not equal to global weights', 1) 634 ENDDO 635 ENDDO 636 ENDDO 637 ! 638 IF (nrr > 20 * nc1 * nc2 * nc3) CALL errore('wigner_seitzg', 'too many WS points, try to increase the bound 20*nc1*nc2*nc3', 1) 639 ! 640 ! Now we compute the WS distances 641 wslen(:) = 0.d0 642 DO i = 1, nrr 643 DO ipol = 1, 3 644 DO jpol = 1, 3 645 wslen(i) = wslen(i) + DBLE(irvec(ipol, i)) * adot(ipol, jpol) * DBLE(irvec(jpol, i)) 646 ENDDO 647 ENDDO 648 wslen(i) = DSQRT(wslen(i)) 649 ENDDO 650 ! 651 CALL cryst_to_cart(dims2, tau(:, :), at, 1) 652 CALL cryst_to_cart(dims, w_centers(:, :), at, 1) 653 ! 654 DEALLOCATE(ind, STAT = ierr) 655 IF (ierr /= 0) CALL errore('wigner_seitzg', 'Error deallocating ind', 1) 656 ! 657 !------------------------------------------------------------------------------------------ 658 END SUBROUTINE wigner_seitzg 659 !------------------------------------------------------------------------------------------ 660 ! 661 ! ----------------------------------------------------------------------------------------- 662 SUBROUTINE backtows(v, ws, rws, nrwsx, nrws) 663 ! ----------------------------------------------------------------------------------------- 664 !! 03/2019 - F. Macheda and S. Ponce 665 !! This subroutines takes a vector "v" like a k-point and bring it back to the 666 !! first Wigner-Seitz cell. 667 ! ----------------------------------------------------------------------------------------- 668 USE kinds, ONLY : DP 669 USE cell_base, ONLY : bg 670 USE constants_epw, ONLY : zero, eps8 671 USE low_lvl, ONLY : find_minimum 672 ! 673 IMPLICIT NONE 674 ! 675 REAL(KIND = DP), INTENT(in) :: v(3) 676 !! Input vector 677 INTEGER, INTENT(in) :: nrws 678 !! Number of WS vectors 679 INTEGER, INTENT(in) :: nrwsx 680 !! Maximum number of WS vector 681 REAL(KIND = DP), INTENT(in) :: rws(4, nrwsx) 682 !! List of WS vectors 683 REAL(KIND = DP), INTENT(out) :: ws(3) 684 !! Vector back into the first WS cell. 685 ! 686 ! Local variables 687 INTEGER :: L, M, N 688 !! l,m,n directions 689 INTEGER :: i 690 !! Index on number of WS vector 691 INTEGER :: minpos 692 !! Minimum position 693 INTEGER, PARAMETER :: nn = 3 694 !! number of neighbours 695 REAL(KIND = DP) :: dist 696 !! Distance 697 REAL(KIND = DP) :: distances(nrws) 698 !! Dsitance array 699 ! 700 ws(:) = zero 701 distances(:) = zero 702 ! 703 alpha : DO L = -nn, nn 704 DO M = -nn, nn 705 DO N = -nn, nn 706 ws(:) = v(:) + L * bg(:, 1) + M * bg(:, 2) + N * bg(:, 3) 707 DO i = 1, nrws 708 distances(i) = DSQRT((ws(1) - rws(2, i))**2 + (ws(2) - rws(3, i))**2 + (ws(3) - rws(4, i))**2) 709 ENDDO 710 minpos = find_minimum(distances, nrws) 711 dist = DSQRT(ws(1)**2 + ws(2)**2 + ws(3)**2) 712 IF (dist < distances(minpos) + eps8) EXIT alpha 713 ENDDO ! N 714 ENDDO ! M 715 ENDDO alpha ! L 716 ! 717 !----------------------------------------------------------------------------------------- 718 END SUBROUTINE backtoWS 719 !----------------------------------------------------------------------------------------- 720 ! 721 !------------------------------------------------------------------------------------------- 722 END MODULE wigner 723 !------------------------------------------------------------------------------------------- 724