LCOV - code coverage report
Current view: top level - src - wannier90.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 477 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 9 0

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

Generated by: LCOV version 2.0-1