1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Outtakes from Wannier90 code 8!> \par History 9!> 06.2016 created [JGH] 10!> \author JGH 11! ************************************************************************************************** 12!-*- mode: F90; mode: font-lock; column-number-mode: true -*-! 13! ! 14! WANNIER90 ! 15! ! 16! The Maximally-Localised Generalised ! 17! Wannier Functions Code ! 18! ! 19! Wannier90 v2.0 authors: ! 20! Arash A. Mostofi (Imperial College London) ! 21! Jonathan R. Yates (University of Oxford) ! 22! Giovanni Pizzi (EPFL, Switzerland) ! 23! Ivo Souza (Universidad del Pais Vasco) ! 24! ! 25! Contributors: ! 26! Young-Su Lee (KIST, S. Korea) ! 27! Matthew Shelley (Imperial College London) ! 28! Nicolas Poilvert (Penn State University) ! 29! Raffaello Bianco (Paris 6 and CNRS) ! 30! Gabriele Sclauzero (ETH Zurich) ! 31! ! 32! Please cite ! 33! ! 34! [ref] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, ! 35! D. Vanderbilt and N. Marzari, "Wannier90: A Tool ! 36! for Obtaining Maximally Localised Wannier ! 37! Functions", Computer Physics Communications, ! 38! 178, 685 (2008) ! 39! ! 40! in any publications arising from the use of this code. ! 41! ! 42! ! 43! Wannier90 is based on Wannier77, written by N. Marzari, ! 44! I. Souza and D. Vanderbilt. For the method please cite ! 45! ! 46! [ref] N. Marzari and D. Vanderbilt, ! 47! Phys. Rev. B 56 12847 (1997) ! 48! ! 49! [ref] I. Souza, N. Marzari and D. Vanderbilt, ! 50! Phys. Rev. B 65 035109 (2001) ! 51! ! 52! ! 53! Copyright (C) 2007-13 Jonathan Yates, Arash Mostofi, ! 54! Giovanni Pizzi, Young-Su Lee, ! 55! Nicola Marzari, Ivo Souza, David Vanderbilt ! 56! ! 57! This file is distributed under the terms of the GNU ! 58! General Public License. See the file `LICENSE' in ! 59! the root directory of the present distribution, or ! 60! http://www.gnu.org/copyleft/gpl.txt . ! 61! ! 62!------------------------------------------------------------! 63 64MODULE wannier90 65 USE kinds, ONLY: dp 66 USE physcon, ONLY: bohr 67#include "./base/base_uses.f90" 68 69 IMPLICIT NONE 70 PRIVATE 71 72 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'wannier90' 73 74 !Input 75 INTEGER :: iprint 76 CHARACTER(len=20) :: length_unit 77 INTEGER :: stdout 78 79 !parameters dervied from input 80 INTEGER :: num_kpts 81 REAL(kind=dp) :: real_lattice(3, 3) 82 REAL(kind=dp) :: recip_lattice(3, 3) 83 REAL(kind=dp) :: cell_volume 84 REAL(kind=dp), ALLOCATABLE ::kpt_cart(:, :) !kpoints in cartesians 85 REAL(kind=dp) :: lenconfac 86 INTEGER :: num_exclude_bands 87 88 REAL(kind=dp) :: kmesh_tol 89 INTEGER :: num_shells 90 INTEGER :: mp_grid(3) 91 INTEGER :: search_shells 92 INTEGER, ALLOCATABLE :: shell_list(:) 93 REAL(kind=dp), ALLOCATABLE :: kpt_latt(:, :) !kpoints in lattice vecs 94 95 ! kmesh parameters (set in kmesh) 96 INTEGER :: nnh ! the number of b-directions (bka) 97 INTEGER :: nntot ! total number of neighbours for each k-point 98 INTEGER, ALLOCATABLE :: nnlist(:, :) ! list of neighbours for each k-point 99 INTEGER, ALLOCATABLE :: neigh(:, :) 100 INTEGER, ALLOCATABLE :: nncell(:, :, :) ! gives BZ of each neighbour of each k-point 101 REAL(kind=dp) :: wbtot 102 REAL(kind=dp), ALLOCATABLE :: wb(:) ! weights associated with neighbours of each k-point 103 REAL(kind=dp), ALLOCATABLE :: bk(:, :, :) ! the b-vectors that go from each k-point to its neighbours 104 REAL(kind=dp), ALLOCATABLE :: bka(:, :) ! the b-directions from 1st k-point to its neighbours 105 106 ! The maximum number of shells we need to satisfy B1 condition in kmesh 107 INTEGER, PARAMETER :: max_shells = 6 108 INTEGER, PARAMETER :: num_nnmax = 12 109 110 INTEGER, PARAMETER :: nsupcell = 5 111 INTEGER :: lmn(3, (2*nsupcell + 1)**3) 112 113 REAL(kind=dp), PARAMETER :: eps5 = 1.0e-5_dp 114 REAL(kind=dp), PARAMETER :: eps6 = 1.0e-6_dp 115 REAL(kind=dp), PARAMETER :: eps8 = 1.0e-8_dp 116 117 PUBLIC :: w90_write_header, wannier_setup 118 119! ************************************************************************************************** 120 121CONTAINS 122 123! ************************************************************************************************** 124!> \brief ... 125!> \param mp_grid_loc ... 126!> \param num_kpts_loc ... 127!> \param real_lattice_loc ... 128!> \param recip_lattice_loc ... 129!> \param kpt_latt_loc ... 130!> \param nntot_loc ... 131!> \param nnlist_loc ... 132!> \param nncell_loc ... 133!> \param iounit ... 134! ************************************************************************************************** 135 SUBROUTINE wannier_setup(mp_grid_loc, num_kpts_loc, & 136 real_lattice_loc, recip_lattice_loc, kpt_latt_loc, & 137 nntot_loc, nnlist_loc, nncell_loc, iounit) 138 139 INTEGER, DIMENSION(3), INTENT(in) :: mp_grid_loc 140 INTEGER, INTENT(in) :: num_kpts_loc 141 REAL(kind=dp), DIMENSION(3, 3), INTENT(in) :: real_lattice_loc, recip_lattice_loc 142 REAL(kind=dp), DIMENSION(3, num_kpts_loc), & 143 INTENT(in) :: kpt_latt_loc 144 INTEGER, INTENT(out) :: nntot_loc 145 INTEGER, DIMENSION(num_kpts_loc, num_nnmax), & 146 INTENT(out) :: nnlist_loc 147 INTEGER, DIMENSION(3, num_kpts_loc, num_nnmax), & 148 INTENT(out) :: nncell_loc 149 INTEGER, INTENT(in) :: iounit 150 151 INTEGER :: nkp 152 153 ! interface uses atomic units 154 length_unit = 'bohr' 155 lenconfac = 1.0_dp/bohr 156 stdout = iounit 157 158 CALL w90_write_header(stdout) 159 160 WRITE (stdout, '(a/)') ' Setting up k-point neighbours...' 161 162 ! copy local data into module variables 163 mp_grid = mp_grid_loc 164 num_kpts = num_kpts_loc 165 real_lattice = real_lattice_loc 166 recip_lattice = recip_lattice_loc 167 ALLOCATE (kpt_latt(3, num_kpts)) 168 ALLOCATE (kpt_cart(3, num_kpts)) 169 kpt_latt(1:3, 1:num_kpts) = kpt_latt_loc(1:3, 1:num_kpts) 170 DO nkp = 1, num_kpts 171 kpt_cart(:, nkp) = MATMUL(kpt_latt(:, nkp), recip_lattice(:, :)) 172 END DO 173 174 num_shells = 0 175 ALLOCATE (shell_list(max_shells)) 176 177 cell_volume = real_lattice(1, 1)*(real_lattice(2, 2)*real_lattice(3, 3) - real_lattice(3, 2)*real_lattice(2, 3)) + & 178 real_lattice(1, 2)*(real_lattice(2, 3)*real_lattice(3, 1) - real_lattice(3, 3)*real_lattice(2, 1)) + & 179 real_lattice(1, 3)*(real_lattice(2, 1)*real_lattice(3, 2) - real_lattice(3, 1)*real_lattice(2, 2)) 180 iprint = 1 181 search_shells = 12 182 kmesh_tol = 0.000001_dp 183 num_exclude_bands = 0 184 185 CALL w90_param_write(stdout) 186 187 CALL w90_kmesh_get() 188 189 nntot_loc = nntot 190 nnlist_loc = 0 191 nnlist_loc(:, 1:nntot) = nnlist(:, 1:nntot) 192 nncell_loc = 0 193 nncell_loc(:, :, 1:nntot) = nncell(:, :, 1:nntot) 194 195 DEALLOCATE (bk, bka, wb) 196 DEALLOCATE (nncell, neigh, nnlist) 197 DEALLOCATE (kpt_latt, kpt_cart, shell_list) 198 199 WRITE (stdout, '(/a/)') ' Finished setting up k-point neighbours.' 200 201 END SUBROUTINE wannier_setup 202! ************************************************************************************************** 203!> \brief ... 204!> \param stdout ... 205! ************************************************************************************************** 206 SUBROUTINE w90_write_header(stdout) 207 INTEGER, INTENT(IN) :: stdout 208 209 WRITE (stdout, *) 210 WRITE (stdout, *) ' +---------------------------------------------------+' 211 WRITE (stdout, *) ' | |' 212 WRITE (stdout, *) ' | WANNIER90 |' 213 WRITE (stdout, *) ' | |' 214 WRITE (stdout, *) ' +---------------------------------------------------+' 215 WRITE (stdout, *) ' | |' 216 WRITE (stdout, *) ' | Welcome to the Maximally-Localized |' 217 WRITE (stdout, *) ' | Generalized Wannier Functions code |' 218 WRITE (stdout, *) ' | http://www.wannier.org |' 219 WRITE (stdout, *) ' | |' 220 WRITE (stdout, *) ' | Wannier90 v2.0 Authors: |' 221 WRITE (stdout, *) ' | Arash A. Mostofi (Imperial College London) |' 222 WRITE (stdout, *) ' | Giovanni Pizzi (EPFL) |' 223 WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |' 224 WRITE (stdout, *) ' | Jonathan R. Yates (University of Oxford) |' 225 WRITE (stdout, *) ' | |' 226 WRITE (stdout, *) ' | Wannier90 Contributors: |' 227 WRITE (stdout, *) ' | Young-Su Lee (KIST, S. Korea) |' 228 WRITE (stdout, *) ' | Matthew Shelley (Imperial College London) |' 229 WRITE (stdout, *) ' | Nicolas Poilvert (Penn State University) |' 230 WRITE (stdout, *) ' | Raffaello Bianco (Paris 6 and CNRS) |' 231 WRITE (stdout, *) ' | Gabriele Sclauzero (ETH Zurich) |' 232 WRITE (stdout, *) ' | |' 233 WRITE (stdout, *) ' | Wannier77 Authors: |' 234 WRITE (stdout, *) ' | Nicola Marzari (EPFL) |' 235 WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |' 236 WRITE (stdout, *) ' | David Vanderbilt (Rutgers University) |' 237 WRITE (stdout, *) ' | |' 238 WRITE (stdout, *) ' | Please cite |' 239 WRITE (stdout, *) ' | |' 240 WRITE (stdout, *) ' | [ref] "Wannier90: A Tool for Obtaining Maximally |' 241 WRITE (stdout, *) ' | Localised Wannier Functions" |' 242 WRITE (stdout, *) ' | A. A. Mostofi, J. R. Yates, Y.-S. Lee, |' 243 WRITE (stdout, *) ' | I. Souza, D. Vanderbilt and N. Marzari |' 244 WRITE (stdout, *) ' | Comput. Phys. Commun. 178, 685 (2008) |' 245 WRITE (stdout, *) ' | |' 246 WRITE (stdout, *) ' | in any publications arising from the use of |' 247 WRITE (stdout, *) ' | this code. For the method please cite |' 248 WRITE (stdout, *) ' | |' 249 WRITE (stdout, *) ' | [ref] "Maximally Localized Generalised Wannier |' 250 WRITE (stdout, *) ' | Functions for Composite Energy Bands" |' 251 WRITE (stdout, *) ' | N. Marzari and D. Vanderbilt |' 252 WRITE (stdout, *) ' | Phys. Rev. B 56 12847 (1997) |' 253 WRITE (stdout, *) ' | |' 254 WRITE (stdout, *) ' | [ref] "Maximally Localized Wannier Functions |' 255 WRITE (stdout, *) ' | for Entangled Energy Bands" |' 256 WRITE (stdout, *) ' | I. Souza, N. Marzari and D. Vanderbilt |' 257 WRITE (stdout, *) ' | Phys. Rev. B 65 035109 (2001) |' 258 WRITE (stdout, *) ' | |' 259 WRITE (stdout, *) ' | |' 260 WRITE (stdout, *) ' | Copyright (c) 1996-2015 |' 261 WRITE (stdout, *) ' | Arash A. Mostofi, Jonathan R. Yates, |' 262 WRITE (stdout, *) ' | Young-Su Lee, Giovanni Pizzi, Ivo Souza, |' 263 WRITE (stdout, *) ' | David Vanderbilt and Nicola Marzari |' 264 WRITE (stdout, *) ' | |' 265 WRITE (stdout, *) ' | Release: 2.0.1 2nd April 2015 |' 266 WRITE (stdout, *) ' | |' 267 WRITE (stdout, *) ' | This program is free software; you can |' 268 WRITE (stdout, *) ' | redistribute it and/or modify it under the terms |' 269 WRITE (stdout, *) ' | of the GNU General Public License as published by |' 270 WRITE (stdout, *) ' | the Free Software Foundation; either version 2 of |' 271 WRITE (stdout, *) ' | the License, or (at your option) any later version|' 272 WRITE (stdout, *) ' | |' 273 WRITE (stdout, *) ' | This program is distributed in the hope that it |' 274 WRITE (stdout, *) ' | will be useful, but WITHOUT ANY WARRANTY; without |' 275 WRITE (stdout, *) ' | even the implied warranty of MERCHANTABILITY or |' 276 WRITE (stdout, *) ' | FITNESS FOR A PARTICULAR PURPOSE. See the GNU |' 277 WRITE (stdout, *) ' | General Public License for more details. |' 278 WRITE (stdout, *) ' | |' 279 WRITE (stdout, *) ' | You should have received a copy of the GNU General|' 280 WRITE (stdout, *) ' | Public License along with this program; if not, |' 281 WRITE (stdout, *) ' | write to the Free Software Foundation, Inc., |' 282 WRITE (stdout, *) ' | 675 Mass Ave, Cambridge, MA 02139, USA. |' 283 WRITE (stdout, *) ' | |' 284 WRITE (stdout, *) ' +---------------------------------------------------+' 285 WRITE (stdout, *) '' 286 287 END SUBROUTINE w90_write_header 288 289! ************************************************************************************************** 290!> \brief ... 291!> \param stdout ... 292! ************************************************************************************************** 293 SUBROUTINE w90_param_write(stdout) 294 INTEGER, INTENT(IN) :: stdout 295 296 INTEGER :: i, nkp 297 298 ! System 299 WRITE (stdout, '(36x,a6)') '------' 300 WRITE (stdout, '(36x,a6)') 'SYSTEM' 301 WRITE (stdout, '(36x,a6)') '------' 302 WRITE (stdout, *) 303 WRITE (stdout, '(28x,a22)') 'Lattice Vectors (Bohr)' 304 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_1', (real_lattice(1, I)*lenconfac, i=1, 3) 305 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_2', (real_lattice(2, I)*lenconfac, i=1, 3) 306 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_3', (real_lattice(3, I)*lenconfac, i=1, 3) 307 WRITE (stdout, *) 308 WRITE (stdout, '(19x,a17,3x,f11.5)', advance='no') & 309 'Unit Cell Volume:', cell_volume*lenconfac**3 310 WRITE (stdout, '(2x,a8)') '(Bohr^3)' 311 WRITE (stdout, *) 312 WRITE (stdout, '(22x,a34)') 'Reciprocal-Space Vectors (Bohr^-1)' 313 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_1', (recip_lattice(1, I)/lenconfac, i=1, 3) 314 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_2', (recip_lattice(2, I)/lenconfac, i=1, 3) 315 WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_3', (recip_lattice(3, I)/lenconfac, i=1, 3) 316 WRITE (stdout, *) ' ' 317 WRITE (stdout, *) ' ' 318 ! K-points 319 WRITE (stdout, '(32x,a)') '------------' 320 WRITE (stdout, '(32x,a)') 'K-POINT GRID' 321 WRITE (stdout, '(32x,a)') '------------' 322 WRITE (stdout, *) ' ' 323 WRITE (stdout, '(13x,a,i3,1x,a1,i3,1x,a1,i3,6x,a,i5)') 'Grid size =', mp_grid(1), 'x', mp_grid(2), 'x', mp_grid(3), & 324 'Total points =', num_kpts 325 WRITE (stdout, *) ' ' 326 WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*' 327 WRITE (stdout, '(1x,a)') '| k-point Fractional Coordinate Cartesian Coordinate (Bohr^-1) |' 328 WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+' 329 DO nkp = 1, num_kpts 330 WRITE (stdout, '(1x,a1,i6,1x,3F10.5,3x,a1,1x,3F10.5,4x,a1)') '|', & 331 nkp, kpt_latt(:, nkp), '|', kpt_cart(:, nkp)/lenconfac, '|' 332 END DO 333 WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*' 334 WRITE (stdout, *) ' ' 335 336 END SUBROUTINE w90_param_write 337 338! ************************************************************************************************** 339!> \brief ... 340! ************************************************************************************************** 341 SUBROUTINE w90_kmesh_get() 342 343 ! Variables that are private 344 345 REAL(kind=dp), PARAMETER :: eta = 99999999.0_dp 346 347 INTEGER :: counter, i, ifound, j, l, loop, loop_b, loop_s, m, multi(search_shells), n, na, & 348 nap, ndnn, ndnntot, ndnnx, nkp, nkp2, nlist, nn, nnsh, nnshell(num_kpts, search_shells), & 349 nnx, shell 350 LOGICAL :: isneg, ispos 351 REAL(kind=dp) :: bb1, bbn, bk_local(3, num_nnmax, num_kpts), bweight(max_shells), ddelta, & 352 dist, dnn(search_shells), dnn0, dnn1, vkpp(3), vkpp2(3), wb_local(num_nnmax) 353 REAL(kind=dp), ALLOCATABLE :: bvec_tmp(:, :) 354 355 WRITE (stdout, '(/1x,a)') & 356 '*---------------------------------- K-MESH ----------------------------------*' 357 358 ! Sort the cell neighbours so we loop in order of distance from the home shell 359 CALL w90_kmesh_supercell_sort 360 361 ! find the distance between k-point 1 and its nearest-neighbour shells 362 ! if we have only one k-point, the n-neighbours are its periodic images 363 364 dnn0 = 0.0_dp 365 dnn1 = eta 366 ndnntot = 0 367 DO nlist = 1, search_shells 368 DO nkp = 1, num_kpts 369 DO loop = 1, (2*nsupcell + 1)**3 370 l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop) 371 ! 372 vkpp = kpt_cart(:, nkp) + MATMUL(lmn(:, loop), recip_lattice) 373 dist = SQRT((kpt_cart(1, 1) - vkpp(1))**2 & 374 + (kpt_cart(2, 1) - vkpp(2))**2 + (kpt_cart(3, 1) - vkpp(3))**2) 375 ! 376 IF ((dist .GT. kmesh_tol) .AND. (dist .GT. dnn0 + kmesh_tol)) THEN 377 IF (dist .LT. dnn1 - kmesh_tol) THEN 378 dnn1 = dist ! found a closer shell 379 counter = 0 380 END IF 381 IF (dist .GT. (dnn1 - kmesh_tol) .AND. dist .LT. (dnn1 + kmesh_tol)) THEN 382 counter = counter + 1 ! count the multiplicity of the shell 383 END IF 384 END IF 385 ENDDO 386 ENDDO 387 IF (dnn1 .LT. eta - kmesh_tol) ndnntot = ndnntot + 1 388 dnn(nlist) = dnn1 389 multi(nlist) = counter 390 dnn0 = dnn1 391 dnn1 = eta 392 ENDDO 393 394 WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+' 395 WRITE (stdout, '(1x,a)') '| Distance to Nearest-Neighbour Shells |' 396 WRITE (stdout, '(1x,a)') '| ------------------------------------ |' 397 WRITE (stdout, '(1x,a)') '| Shell Distance (Bohr^-1) Multiplicity |' 398 WRITE (stdout, '(1x,a)') '| ----- ------------------ ------------ |' 399 DO ndnn = 1, ndnntot 400 WRITE (stdout, '(1x,a,11x,i3,17x,f10.6,19x,i4,12x,a)') '|', ndnn, dnn(ndnn)/lenconfac, multi(ndnn), '|' 401 ENDDO 402 WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+' 403 404 ! Get the shell weights to satisfy the B1 condition 405 CALL kmesh_shell_automatic(multi, dnn, bweight) 406 407 WRITE (stdout, '(1x,a)', advance='no') '| The following shells are used: ' 408 DO ndnn = 1, num_shells 409 IF (ndnn .EQ. num_shells) THEN 410 WRITE (stdout, '(i3,1x)', advance='no') shell_list(ndnn) 411 ELSE 412 WRITE (stdout, '(i3,",")', advance='no') shell_list(ndnn) 413 ENDIF 414 ENDDO 415 DO l = 1, 11 - num_shells 416 WRITE (stdout, '(4x)', advance='no') 417 ENDDO 418 WRITE (stdout, '("|")') 419 420 nntot = 0 421 DO loop_s = 1, num_shells 422 nntot = nntot + multi(shell_list(loop_s)) 423 END DO 424 IF (nntot > num_nnmax) THEN 425 WRITE (stdout, '(a,i2,a)') ' **WARNING: kmesh has found >', num_nnmax, ' nearest neighbours**' 426 WRITE (stdout, '(a)') ' ' 427 WRITE (stdout, '(a)') ' This is probably caused by an error in your unit cell specification' 428 WRITE (stdout, '(a)') ' ' 429 430 ALLOCATE (bvec_tmp(3, MAXVAL(multi))) 431 bvec_tmp = 0.0_dp 432 counter = 0 433 DO shell = 1, search_shells 434 CALL kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvec_tmp(:, 1:multi(shell))) 435 DO loop = 1, multi(shell) 436 counter = counter + 1 437 WRITE (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector ', counter, ': (', & 438 bvec_tmp(:, loop)/lenconfac, ')', dnn(shell)/lenconfac, ' |' 439 END DO 440 END DO 441 WRITE (stdout, '(a)') ' ' 442 DEALLOCATE (bvec_tmp) 443 CPABORT('kmesh_get: something wrong, found too many nearest neighbours') 444 END IF 445 446 ALLOCATE (nnlist(num_kpts, nntot)) 447 ALLOCATE (neigh(num_kpts, nntot/2)) 448 ALLOCATE (nncell(3, num_kpts, nntot)) 449 450 ALLOCATE (wb(nntot)) 451 ALLOCATE (bka(3, nntot/2)) 452 ALLOCATE (bk(3, nntot, num_kpts)) 453 454 nnx = 0 455 DO loop_s = 1, num_shells 456 DO loop_b = 1, multi(shell_list(loop_s)) 457 nnx = nnx + 1 458 wb_local(nnx) = bweight(loop_s) 459 END DO 460 END DO 461 462 ! Now build up the list of nearest-neighbour shells for each k-point. 463 ! nnlist(nkp,1...nnx) points to the nnx neighbours (ordered along increa 464 ! shells) of the k-point nkp. nncell(i,nkp,nnth) tells us in which BZ is 465 ! nnth nearest-neighbour of the k-point nkp. Construct the nnx b-vectors 466 ! go from k-point nkp to each neighbour bk(1:3,nkp,1...nnx). 467 ! Comment: Now we have bk(3,nntot,num_kps) 09/04/2006 468 469 WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+' 470 WRITE (stdout, '(1x,a)') '| Shell # Nearest-Neighbours |' 471 WRITE (stdout, '(1x,a)') '| ----- -------------------- |' 472 ! 473 ! Standard routine 474 ! 475 nnshell = 0 476 DO nkp = 1, num_kpts 477 nnx = 0 478 ok: DO ndnnx = 1, num_shells 479 ndnn = shell_list(ndnnx) 480 DO loop = 1, (2*nsupcell + 1)**3 481 l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop) 482 vkpp2 = MATMUL(lmn(:, loop), recip_lattice) 483 DO nkp2 = 1, num_kpts 484 vkpp = vkpp2 + kpt_cart(:, nkp2) 485 dist = SQRT((kpt_cart(1, nkp) - vkpp(1))**2 & 486 + (kpt_cart(2, nkp) - vkpp(2))**2 + (kpt_cart(3, nkp) - vkpp(3))**2) 487 IF ((dist .GE. dnn(ndnn)*(1 - kmesh_tol)) .AND. (dist .LE. dnn(ndnn)*(1 + kmesh_tol))) THEN 488 nnx = nnx + 1 489 nnshell(nkp, ndnn) = nnshell(nkp, ndnn) + 1 490 nnlist(nkp, nnx) = nkp2 491 nncell(1, nkp, nnx) = l 492 nncell(2, nkp, nnx) = m 493 nncell(3, nkp, nnx) = n 494 bk_local(:, nnx, nkp) = vkpp(:) - kpt_cart(:, nkp) 495 ENDIF 496 !if we have the right number of neighbours we can exit 497 IF (nnshell(nkp, ndnn) == multi(ndnn)) CYCLE ok 498 ENDDO 499 ENDDO 500 ! check to see if too few neighbours here 501 END DO ok 502 503 END DO 504 505 DO ndnnx = 1, num_shells 506 ndnn = shell_list(ndnnx) 507 WRITE (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|' 508 END DO 509 WRITE (stdout, '(1x,"+",76("-"),"+")') 510 511 DO nkp = 1, num_kpts 512 nnx = 0 513 DO ndnnx = 1, num_shells 514 ndnn = shell_list(ndnnx) 515 DO nnsh = 1, nnshell(nkp, ndnn) 516 bb1 = 0.0_dp 517 bbn = 0.0_dp 518 nnx = nnx + 1 519 DO i = 1, 3 520 bb1 = bb1 + bk_local(i, nnx, 1)*bk_local(i, nnx, 1) 521 bbn = bbn + bk_local(i, nnx, nkp)*bk_local(i, nnx, nkp) 522 ENDDO 523 IF (ABS(SQRT(bb1) - SQRT(bbn)) .GT. kmesh_tol) THEN 524 WRITE (stdout, '(1x,2f10.6)') bb1, bbn 525 CPABORT('Non-symmetric k-point neighbours!') 526 ENDIF 527 ENDDO 528 ENDDO 529 ENDDO 530 531 ! now check that the completeness relation is satisfied for every kpoint 532 ! We know it is true for kpt=1; but we check the rest to be safe. 533 ! Eq. B1 in Appendix B PRB 56 12847 (1997) 534 535 DO nkp = 1, num_kpts 536 DO i = 1, 3 537 DO j = 1, 3 538 ddelta = 0.0_dp 539 nnx = 0 540 DO ndnnx = 1, num_shells 541 ndnn = shell_list(ndnnx) 542 DO nnsh = 1, nnshell(1, ndnn) 543 nnx = nnx + 1 544 ddelta = ddelta + wb_local(nnx)*bk_local(i, nnx, nkp)*bk_local(j, nnx, nkp) 545 ENDDO 546 ENDDO 547 IF ((i .EQ. j) .AND. (ABS(ddelta - 1.0_dp) .GT. kmesh_tol)) THEN 548 WRITE (stdout, '(1x,2i3,f12.8)') i, j, ddelta 549 CPABORT('Eq. (B1) not satisfied in kmesh_get (1)') 550 ENDIF 551 IF ((i .NE. j) .AND. (ABS(ddelta) .GT. kmesh_tol)) THEN 552 WRITE (stdout, '(1x,2i3,f12.8)') i, j, ddelta 553 CPABORT('Eq. (B1) not satisfied in kmesh_get (2)') 554 ENDIF 555 ENDDO 556 ENDDO 557 ENDDO 558 559 WRITE (stdout, '(1x,a)') '| Completeness relation is fully satisfied [Eq. (B1), PRB 56, 12847 (1997)] |' 560 WRITE (stdout, '(1x,"+",76("-"),"+")') 561 562 ! 563 wbtot = 0.0_dp 564 nnx = 0 565 DO ndnnx = 1, num_shells 566 ndnn = shell_list(ndnnx) 567 DO nnsh = 1, nnshell(1, ndnn) 568 nnx = nnx + 1 569 wbtot = wbtot + wb_local(nnx) 570 ENDDO 571 ENDDO 572 573 nnh = nntot/2 574 ! make list of bka vectors from neighbours of first k-point 575 ! delete any inverse vectors as you collect them 576 na = 0 577 DO nn = 1, nntot 578 ifound = 0 579 IF (na .NE. 0) THEN 580 DO nap = 1, na 581 CALL utility_compar(bka(1, nap), bk_local(1, nn, 1), ispos, isneg) 582 IF (isneg) ifound = 1 583 ENDDO 584 ENDIF 585 IF (ifound .EQ. 0) THEN 586 ! found new vector to add to set 587 na = na + 1 588 bka(1, na) = bk_local(1, nn, 1) 589 bka(2, na) = bk_local(2, nn, 1) 590 bka(3, na) = bk_local(3, nn, 1) 591 ENDIF 592 ENDDO 593 IF (na .NE. nnh) CPABORT('Did not find right number of bk directions') 594 595 WRITE (stdout, '(1x,a)') '| b_k Vectors (Bohr^-1) and Weights (Bohr^2) |' 596 WRITE (stdout, '(1x,a)') '| ------------------------------------------ |' 597 WRITE (stdout, '(1x,a)') '| No. b_k(x) b_k(y) b_k(z) w_b |' 598 WRITE (stdout, '(1x,a)') '| --- -------------------------------- -------- |' 599 DO i = 1, nntot 600 WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') & 601 i, (bk_local(j, i, 1)/lenconfac, j=1, 3), wb_local(i)*lenconfac**2 602 ENDDO 603 WRITE (stdout, '(1x,"+",76("-"),"+")') 604 WRITE (stdout, '(1x,a)') '| b_k Directions (Bohr^-1) |' 605 WRITE (stdout, '(1x,a)') '| ------------------------ |' 606 WRITE (stdout, '(1x,a)') '| No. x y z |' 607 WRITE (stdout, '(1x,a)') '| --- -------------------------------- |' 608 DO i = 1, nnh 609 WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,21x,"|")') i, (bka(j, i)/lenconfac, j=1, 3) 610 ENDDO 611 WRITE (stdout, '(1x,"+",76("-"),"+")') 612 WRITE (stdout, *) ' ' 613 614 ! find index array 615 DO nkp = 1, num_kpts 616 DO na = 1, nnh 617 ! first, zero the index array so we can check it gets filled 618 neigh(nkp, na) = 0 619 ! now search through list of neighbours of this k-point 620 DO nn = 1, nntot 621 CALL utility_compar(bka(1, na), bk_local(1, nn, nkp), ispos, isneg) 622 IF (ispos) neigh(nkp, na) = nn 623 ENDDO 624 ! check found 625 IF (neigh(nkp, na) .EQ. 0) THEN 626 WRITE (stdout, *) ' nkp,na=', nkp, na 627 CPABORT('kmesh_get: failed to find neighbours for this kpoint') 628 ENDIF 629 ENDDO 630 ENDDO 631 632 !fill in the global arrays from the local ones 633 DO loop = 1, nntot 634 wb(loop) = wb_local(loop) 635 END DO 636 637 DO loop_s = 1, num_kpts 638 DO loop = 1, nntot 639 bk(:, loop, loop_s) = bk_local(:, loop, loop_s) 640 END DO 641 END DO 642 643 END SUBROUTINE w90_kmesh_get 644 645! ************************************************************************************************** 646!> \brief ... 647! ************************************************************************************************** 648 SUBROUTINE w90_kmesh_supercell_sort 649 !==================================================================! 650 ! ! 651 ! We look for kpoint neighbours in a large supercell of reciprocal ! 652 ! unit cells. Done sequentially this is very slow. ! 653 ! Here we order the cells by the distance from the origin ! 654 ! Doing the search in this order gives a dramatic speed up ! 655 ! ! 656 !==================================================================! 657 INTEGER :: counter, indx(1), l, & 658 lmn_cp(3, (2*nsupcell + 1)**3), loop, & 659 m, n 660 REAL(kind=dp) :: dist((2*nsupcell + 1)**3), & 661 dist_cp((2*nsupcell + 1)**3), pos(3) 662 663 counter = 1 664 lmn(:, counter) = 0 665 dist(counter) = 0.0_dp 666 DO l = -nsupcell, nsupcell 667 DO m = -nsupcell, nsupcell 668 DO n = -nsupcell, nsupcell 669 IF (l == 0 .AND. m == 0 .AND. n == 0) CYCLE 670 counter = counter + 1 671 lmn(1, counter) = l; lmn(2, counter) = m; lmn(3, counter) = n 672 pos = MATMUL(lmn(:, counter), recip_lattice) 673 dist(counter) = SQRT(DOT_PRODUCT(pos, pos)) 674 END DO 675 END DO 676 END DO 677 678 DO loop = (2*nsupcell + 1)**3, 1, -1 679 indx = internal_maxloc(dist) 680 dist_cp(loop) = dist(indx(1)) 681 lmn_cp(:, loop) = lmn(:, indx(1)) 682 dist(indx(1)) = -1.0_dp 683 END DO 684 685 lmn = lmn_cp 686 dist = dist_cp 687 688 END SUBROUTINE w90_kmesh_supercell_sort 689 690! ************************************************************************************************** 691!> \brief ... 692!> \param dist ... 693!> \return ... 694! ************************************************************************************************** 695 FUNCTION internal_maxloc(dist) 696 !=========================================================================! 697 ! ! 698 ! A predictable maxloc. ! 699 ! ! 700 !=========================================================================! 701 702 REAL(kind=dp), INTENT(in) :: dist((2*nsupcell + 1)**3) 703 INTEGER :: internal_maxloc 704 705 INTEGER :: counter, guess(1), & 706 list((2*nsupcell + 1)**3), loop 707 708 list = 0 709 counter = 1 710 711 guess = MAXLOC(dist) 712 list(1) = guess(1) 713 ! look for any degenerate values 714 DO loop = 1, (2*nsupcell + 1)**3 715 IF (loop == guess(1)) CYCLE 716 IF (ABS(dist(loop) - dist(guess(1))) < eps8) THEN 717 counter = counter + 1 718 list(counter) = loop 719 ENDIF 720 END DO 721 ! and always return the lowest index 722 internal_maxloc = MINVAL(list(1:counter)) 723 724 END FUNCTION internal_maxloc 725 726! ************************************************************************************************** 727!> \brief ... 728!> \param multi ... 729!> \param dnn ... 730!> \param bweight ... 731! ************************************************************************************************** 732 SUBROUTINE kmesh_shell_automatic(multi, dnn, bweight) 733 !==========================================================================! 734 ! ! 735 ! Find the correct set of shells to satisfy B1 ! 736 ! The stratagy is: ! 737 ! Take the bvectors from the next shell ! 738 ! Reject them if they are parallel to exisiting b vectors ! 739 ! Test to see if we satisfy B1, if not add another shell and repeat ! 740 ! ! 741 !==========================================================================! 742 743 INTEGER, INTENT(in) :: multi(search_shells) 744 REAL(kind=dp), INTENT(in) :: dnn(search_shells) 745 REAL(kind=dp), INTENT(out) :: bweight(max_shells) 746 747 INTEGER, PARAMETER :: lwork = max_shells*10 748 REAL(kind=dp), PARAMETER :: TARGET(6) = (/1.0_dp, 1.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/) 749 750 INTEGER :: cur_shell, info, loop_b, loop_bn, & 751 loop_i, loop_j, loop_s, shell 752 LOGICAL :: b1sat, lpar 753 REAL(kind=dp) :: delta, work(lwork) 754 REAL(kind=dp), ALLOCATABLE :: bvector(:, :, :) 755 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: singv 756 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, smat, umat, vmat 757 758 ALLOCATE (bvector(3, MAXVAL(multi), max_shells)) 759 bvector = 0.0_dp; bweight = 0.0_dp 760 761 WRITE (stdout, '(1x,a)') '| The b-vectors are chosen automatically |' 762 763 b1sat = .FALSE. 764 DO shell = 1, search_shells 765 cur_shell = num_shells + 1 766 767 ! get the b vectors for the new shell 768 CALL kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvector(:, 1:multi(shell), cur_shell)) 769 770 ! We check that the new shell is not parrallel to an existing shell (cosine=1) 771 lpar = .FALSE. 772 IF (num_shells > 0) THEN 773 DO loop_bn = 1, multi(shell) 774 DO loop_s = 1, num_shells 775 DO loop_b = 1, multi(shell_list(loop_s)) 776 delta = DOT_PRODUCT(bvector(:, loop_bn, cur_shell), bvector(:, loop_b, loop_s))/ & 777 SQRT(DOT_PRODUCT(bvector(:, loop_bn, cur_shell), bvector(:, loop_bn, cur_shell))* & 778 DOT_PRODUCT(bvector(:, loop_b, loop_s), bvector(:, loop_b, loop_s))) 779 IF (ABS(ABS(delta) - 1.0_dp) < eps6) lpar = .TRUE. 780 END DO 781 END DO 782 END DO 783 END IF 784 785 IF (lpar) THEN 786 IF (iprint >= 3) THEN 787 WRITE (stdout, '(1x,a)') '| This shell is linearly dependent on existing shells: Trying next shell |' 788 END IF 789 CYCLE 790 END IF 791 792 num_shells = num_shells + 1 793 shell_list(num_shells) = shell 794 795 ALLOCATE (amat(max_shells, num_shells)) 796 ALLOCATE (umat(max_shells, max_shells)) 797 ALLOCATE (vmat(num_shells, num_shells)) 798 ALLOCATE (smat(num_shells, max_shells)) 799 ALLOCATE (singv(num_shells)) 800 amat = 0.0_dp; umat = 0.0_dp; vmat = 0.0_dp; smat = 0.0_dp; singv = 0.0_dp 801 802 amat = 0.0_dp 803 DO loop_s = 1, num_shells 804 DO loop_b = 1, multi(shell_list(loop_s)) 805 amat(1, loop_s) = amat(1, loop_s) + bvector(1, loop_b, loop_s)*bvector(1, loop_b, loop_s) 806 amat(2, loop_s) = amat(2, loop_s) + bvector(2, loop_b, loop_s)*bvector(2, loop_b, loop_s) 807 amat(3, loop_s) = amat(3, loop_s) + bvector(3, loop_b, loop_s)*bvector(3, loop_b, loop_s) 808 amat(4, loop_s) = amat(4, loop_s) + bvector(1, loop_b, loop_s)*bvector(2, loop_b, loop_s) 809 amat(5, loop_s) = amat(5, loop_s) + bvector(2, loop_b, loop_s)*bvector(3, loop_b, loop_s) 810 amat(6, loop_s) = amat(6, loop_s) + bvector(3, loop_b, loop_s)*bvector(1, loop_b, loop_s) 811 END DO 812 END DO 813 814 info = 0 815 CALL dgesvd('A', 'A', max_shells, num_shells, amat, max_shells, singv, umat, & 816 max_shells, vmat, num_shells, work, lwork, info) 817 IF (info < 0) THEN 818 WRITE (stdout, '(1x,a,1x,I1,1x,a)') 'kmesh_shell_automatic: Argument', ABS(info), 'of dgesvd is incorrect' 819 CPABORT('kmesh_shell_automatic: Problem with Singular Value Decomposition') 820 ELSE IF (info > 0) THEN 821 CPABORT('kmesh_shell_automatic: Singular Value Decomposition did not converge') 822 END IF 823 824 IF (ANY(ABS(singv) < eps5)) THEN 825 IF (num_shells == 1) THEN 826 CALL cp_abort(__LOCATION__, "kmesh_shell_automatic: "// & 827 "Singular Value Decomposition has found a very small singular value.") 828 ELSE 829 WRITE (stdout, '(1x,a)') '| SVD found small singular value, Rejecting this shell and trying the next |' 830 b1sat = .FALSE. 831 num_shells = num_shells - 1 832 DEALLOCATE (amat, umat, vmat, smat, singv) 833 CYCLE 834 END IF 835 END IF 836 837 smat = 0.0_dp 838 DO loop_s = 1, num_shells 839 smat(loop_s, loop_s) = 1/singv(loop_s) 840 END DO 841 842 bweight(1:num_shells) = MATMUL(TRANSPOSE(vmat), MATMUL(smat, MATMUL(TRANSPOSE(umat), TARGET))) 843 IF (iprint >= 2) THEN 844 DO loop_s = 1, num_shells 845 WRITE (stdout, '(1x,a,I2,a,f12.7,5x,a8,36x,a)') '| Shell: ', loop_s, & 846 ' w_b ', bweight(loop_s)*lenconfac**2, '('//TRIM(length_unit)//'^2)', '|' 847 END DO 848 END IF 849 850 !check b1 851 b1sat = .TRUE. 852 DO loop_i = 1, 3 853 DO loop_j = loop_i, 3 854 delta = 0.0_dp 855 DO loop_s = 1, num_shells 856 DO loop_b = 1, multi(shell_list(loop_s)) 857 delta = delta + bweight(loop_s)*bvector(loop_i, loop_b, loop_s)*bvector(loop_j, loop_b, loop_s) 858 END DO 859 END DO 860 IF (loop_i == loop_j) THEN 861 IF (ABS(delta - 1.0_dp) > kmesh_tol) b1sat = .FALSE. 862 END IF 863 IF (loop_i /= loop_j) THEN 864 IF (ABS(delta) > kmesh_tol) b1sat = .FALSE. 865 END IF 866 END DO 867 END DO 868 869 IF (.NOT. b1sat) THEN 870 IF (shell < search_shells .AND. iprint >= 3) THEN 871 WRITE (stdout, '(1x,a,24x,a1)') '| B1 condition is not satisfied: Adding another shell', '|' 872 ELSEIF (shell == search_shells) THEN 873 WRITE (stdout, *) ' ' 874 WRITE (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells' 875 WRITE (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid' 876 WRITE (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=12)' 877 WRITE (stdout, *) ' ' 878 CPABORT('kmesh_get_automatic') 879 END IF 880 END IF 881 882 DEALLOCATE (amat, umat, vmat, smat, singv) 883 884 IF (b1sat) EXIT 885 886 END DO 887 888 IF (.NOT. b1sat) THEN 889 WRITE (stdout, *) ' ' 890 WRITE (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells' 891 WRITE (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid' 892 WRITE (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=12)' 893 WRITE (stdout, *) ' ' 894 CPABORT('kmesh_get_automatic') 895 END IF 896 897 END SUBROUTINE kmesh_shell_automatic 898 899! ************************************************************************************************** 900!> \brief ... 901!> \param multi ... 902!> \param kpt ... 903!> \param shell_dist ... 904!> \param bvector ... 905! ************************************************************************************************** 906 SUBROUTINE kmesh_get_bvectors(multi, kpt, shell_dist, bvector) 907 !==================================================================! 908 ! ! 909 ! Returns the bvectors for a given shell and kpoint ! 910 ! ! 911 !=================================================================== 912 913 INTEGER, INTENT(in) :: multi, kpt 914 REAL(kind=dp), INTENT(in) :: shell_dist 915 REAL(kind=dp), INTENT(out) :: bvector(3, multi) 916 917 INTEGER :: loop, nkp2, num_bvec 918 REAL(kind=dp) :: dist, vkpp(3), vkpp2(3) 919 920 bvector = 0.0_dp 921 922 num_bvec = 0 923 ok: DO loop = 1, (2*nsupcell + 1)**3 924 vkpp2 = MATMUL(lmn(:, loop), recip_lattice) 925 DO nkp2 = 1, num_kpts 926 vkpp = vkpp2 + kpt_cart(:, nkp2) 927 dist = SQRT((kpt_cart(1, kpt) - vkpp(1))**2 & 928 + (kpt_cart(2, kpt) - vkpp(2))**2 + (kpt_cart(3, kpt) - vkpp(3))**2) 929 IF ((dist .GE. shell_dist*(1.0_dp - kmesh_tol)) .AND. dist .LE. shell_dist*(1.0_dp + kmesh_tol)) THEN 930 num_bvec = num_bvec + 1 931 bvector(:, num_bvec) = vkpp(:) - kpt_cart(:, kpt) 932 ENDIF 933 !if we have the right number of neighbours we can exit 934 IF (num_bvec == multi) CYCLE ok 935 ENDDO 936 ENDDO ok 937 938 IF (num_bvec < multi) CPABORT('kmesh_get_bvector: Not enough bvectors found') 939 940 END SUBROUTINE kmesh_get_bvectors 941 942! ************************************************************************************************** 943!> \brief Checks whether a ~= b (ispos) or a ~= -b (isneg) up to a precision of eps8 944!> \param a 3-vector 945!> \param b 3-vector 946!> \param ispos true if |a-b|^2 < eps8, otherwise false 947!> \param isneg true if |a+b|^2 < eps8, otherwise false 948! ************************************************************************************************** 949 PURE SUBROUTINE utility_compar(a, b, ispos, isneg) 950 REAL(kind=dp), DIMENSION(3), INTENT(in) :: a, b 951 LOGICAL, INTENT(out) :: ispos, isneg 952 953 ispos = SUM((a - b)**2) .LT. eps8 954 isneg = SUM((a + b)**2) .LT. eps8 955 END SUBROUTINE utility_compar 956 957! ************************************************************************************************** 958 959END MODULE wannier90 960