(git:1f285aa)
wannier90.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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  SUBROUTINE wannier_setup(mp_grid_loc, num_kpts_loc, &
138  real_lattice_loc, recip_lattice_loc, kpt_latt_loc, &
139  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  length_unit = 'bohr'
157  lenconfac = 1.0_dp/bohr
158  stdout = iounit
159 
160  CALL w90_write_header(stdout)
161 
162  WRITE (stdout, '(a/)') ' Setting up k-point neighbours...'
163 
164  ! copy local data into module variables
165  mp_grid = mp_grid_loc
166  num_kpts = num_kpts_loc
167  real_lattice = real_lattice_loc
168  recip_lattice = recip_lattice_loc
169  ALLOCATE (kpt_latt(3, num_kpts))
170  ALLOCATE (kpt_cart(3, num_kpts))
171  kpt_latt(1:3, 1:num_kpts) = kpt_latt_loc(1:3, 1:num_kpts)
172  DO nkp = 1, num_kpts
173  kpt_cart(:, nkp) = matmul(kpt_latt(:, nkp), recip_lattice(:, :))
174  END DO
175 
176  num_shells = 0
177  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  real_lattice(1, 3)*(real_lattice(2, 1)*real_lattice(3, 2) - real_lattice(3, 1)*real_lattice(2, 2))
182  iprint = 1
183  search_shells = 12
184  kmesh_tol = 0.000001_dp
185  num_exclude_bands = 0
186 
187  CALL w90_param_write(stdout)
188 
189  CALL w90_kmesh_get()
190 
191  nntot_loc = nntot
192  nnlist_loc = 0
193  nnlist_loc(:, 1:nntot) = nnlist(:, 1:nntot)
194  nncell_loc = 0
195  nncell_loc(:, :, 1:nntot) = nncell(:, :, 1:nntot)
196 
197  DEALLOCATE (bk, bka, wb)
198  DEALLOCATE (nncell, neigh, nnlist)
199  DEALLOCATE (kpt_latt, kpt_cart, shell_list)
200 
201  WRITE (stdout, '(/a/)') ' Finished setting up k-point neighbours.'
202 
203  END SUBROUTINE wannier_setup
204 ! **************************************************************************************************
205 !> \brief ...
206 !> \param stdout ...
207 ! **************************************************************************************************
208  SUBROUTINE w90_write_header(stdout)
209  INTEGER, INTENT(IN) :: stdout
210 
211  WRITE (stdout, *)
212  WRITE (stdout, *) ' +---------------------------------------------------+'
213  WRITE (stdout, *) ' | |'
214  WRITE (stdout, *) ' | WANNIER90 |'
215  WRITE (stdout, *) ' | |'
216  WRITE (stdout, *) ' +---------------------------------------------------+'
217  WRITE (stdout, *) ' | |'
218  WRITE (stdout, *) ' | Welcome to the Maximally-Localized |'
219  WRITE (stdout, *) ' | Generalized Wannier Functions code |'
220  WRITE (stdout, *) ' | http://www.wannier.org |'
221  WRITE (stdout, *) ' | |'
222  WRITE (stdout, *) ' | Wannier90 v2.0 Authors: |'
223  WRITE (stdout, *) ' | Arash A. Mostofi (Imperial College London) |'
224  WRITE (stdout, *) ' | Giovanni Pizzi (EPFL) |'
225  WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |'
226  WRITE (stdout, *) ' | Jonathan R. Yates (University of Oxford) |'
227  WRITE (stdout, *) ' | |'
228  WRITE (stdout, *) ' | Wannier90 Contributors: |'
229  WRITE (stdout, *) ' | Young-Su Lee (KIST, S. Korea) |'
230  WRITE (stdout, *) ' | Matthew Shelley (Imperial College London) |'
231  WRITE (stdout, *) ' | Nicolas Poilvert (Penn State University) |'
232  WRITE (stdout, *) ' | Raffaello Bianco (Paris 6 and CNRS) |'
233  WRITE (stdout, *) ' | Gabriele Sclauzero (ETH Zurich) |'
234  WRITE (stdout, *) ' | |'
235  WRITE (stdout, *) ' | Wannier77 Authors: |'
236  WRITE (stdout, *) ' | Nicola Marzari (EPFL) |'
237  WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |'
238  WRITE (stdout, *) ' | David Vanderbilt (Rutgers University) |'
239  WRITE (stdout, *) ' | |'
240  WRITE (stdout, *) ' | Please cite |'
241  WRITE (stdout, *) ' | |'
242  WRITE (stdout, *) ' | [ref] "Wannier90: A Tool for Obtaining Maximally |'
243  WRITE (stdout, *) ' | Localised Wannier Functions" |'
244  WRITE (stdout, *) ' | A. A. Mostofi, J. R. Yates, Y.-S. Lee, |'
245  WRITE (stdout, *) ' | I. Souza, D. Vanderbilt and N. Marzari |'
246  WRITE (stdout, *) ' | Comput. Phys. Commun. 178, 685 (2008) |'
247  WRITE (stdout, *) ' | |'
248  WRITE (stdout, *) ' | in any publications arising from the use of |'
249  WRITE (stdout, *) ' | this code. For the method please cite |'
250  WRITE (stdout, *) ' | |'
251  WRITE (stdout, *) ' | [ref] "Maximally Localized Generalised Wannier |'
252  WRITE (stdout, *) ' | Functions for Composite Energy Bands" |'
253  WRITE (stdout, *) ' | N. Marzari and D. Vanderbilt |'
254  WRITE (stdout, *) ' | Phys. Rev. B 56 12847 (1997) |'
255  WRITE (stdout, *) ' | |'
256  WRITE (stdout, *) ' | [ref] "Maximally Localized Wannier Functions |'
257  WRITE (stdout, *) ' | for Entangled Energy Bands" |'
258  WRITE (stdout, *) ' | I. Souza, N. Marzari and D. Vanderbilt |'
259  WRITE (stdout, *) ' | Phys. Rev. B 65 035109 (2001) |'
260  WRITE (stdout, *) ' | |'
261  WRITE (stdout, *) ' | |'
262  WRITE (stdout, *) ' | Copyright (c) 1996-2015 |'
263  WRITE (stdout, *) ' | Arash A. Mostofi, Jonathan R. Yates, |'
264  WRITE (stdout, *) ' | Young-Su Lee, Giovanni Pizzi, Ivo Souza, |'
265  WRITE (stdout, *) ' | David Vanderbilt and Nicola Marzari |'
266  WRITE (stdout, *) ' | |'
267  WRITE (stdout, *) ' | Release: 2.0.1 2nd April 2015 |'
268  WRITE (stdout, *) ' | |'
269  WRITE (stdout, *) ' | This program is free software; you can |'
270  WRITE (stdout, *) ' | redistribute it and/or modify it under the terms |'
271  WRITE (stdout, *) ' | of the GNU General Public License as published by |'
272  WRITE (stdout, *) ' | the Free Software Foundation; either version 2 of |'
273  WRITE (stdout, *) ' | the License, or (at your option) any later version|'
274  WRITE (stdout, *) ' | |'
275  WRITE (stdout, *) ' | This program is distributed in the hope that it |'
276  WRITE (stdout, *) ' | will be useful, but WITHOUT ANY WARRANTY; without |'
277  WRITE (stdout, *) ' | even the implied warranty of MERCHANTABILITY or |'
278  WRITE (stdout, *) ' | FITNESS FOR A PARTICULAR PURPOSE. See the GNU |'
279  WRITE (stdout, *) ' | General Public License for more details. |'
280  WRITE (stdout, *) ' | |'
281  WRITE (stdout, *) ' | You should have received a copy of the GNU General|'
282  WRITE (stdout, *) ' | Public License along with this program; if not, |'
283  WRITE (stdout, *) ' | write to the Free Software Foundation, Inc., |'
284  WRITE (stdout, *) ' | 675 Mass Ave, Cambridge, MA 02139, USA. |'
285  WRITE (stdout, *) ' | |'
286  WRITE (stdout, *) ' +---------------------------------------------------+'
287  WRITE (stdout, *) ''
288 
289  END SUBROUTINE w90_write_header
290 
291 ! **************************************************************************************************
292 !> \brief ...
293 !> \param stdout ...
294 ! **************************************************************************************************
295  SUBROUTINE w90_param_write(stdout)
296  INTEGER, INTENT(IN) :: stdout
297 
298  INTEGER :: i, nkp
299 
300  ! System
301  WRITE (stdout, '(36x,a6)') '------'
302  WRITE (stdout, '(36x,a6)') 'SYSTEM'
303  WRITE (stdout, '(36x,a6)') '------'
304  WRITE (stdout, *)
305  WRITE (stdout, '(28x,a22)') 'Lattice Vectors (Bohr)'
306  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_1', (real_lattice(1, i)*lenconfac, i=1, 3)
307  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_2', (real_lattice(2, i)*lenconfac, i=1, 3)
308  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_3', (real_lattice(3, i)*lenconfac, i=1, 3)
309  WRITE (stdout, *)
310  WRITE (stdout, '(19x,a17,3x,f11.5)', advance='no') &
311  'Unit Cell Volume:', cell_volume*lenconfac**3
312  WRITE (stdout, '(2x,a8)') '(Bohr^3)'
313  WRITE (stdout, *)
314  WRITE (stdout, '(22x,a34)') 'Reciprocal-Space Vectors (Bohr^-1)'
315  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_1', (recip_lattice(1, i)/lenconfac, i=1, 3)
316  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_2', (recip_lattice(2, i)/lenconfac, i=1, 3)
317  WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_3', (recip_lattice(3, i)/lenconfac, i=1, 3)
318  WRITE (stdout, *) ' '
319  WRITE (stdout, *) ' '
320  ! K-points
321  WRITE (stdout, '(32x,a)') '------------'
322  WRITE (stdout, '(32x,a)') 'K-POINT GRID'
323  WRITE (stdout, '(32x,a)') '------------'
324  WRITE (stdout, *) ' '
325  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  'Total points =', num_kpts
327  WRITE (stdout, *) ' '
328  WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*'
329  WRITE (stdout, '(1x,a)') '| k-point Fractional Coordinate Cartesian Coordinate (Bohr^-1) |'
330  WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
331  DO nkp = 1, num_kpts
332  WRITE (stdout, '(1x,a1,i6,1x,3F10.5,3x,a1,1x,3F10.5,4x,a1)') '|', &
333  nkp, kpt_latt(:, nkp), '|', kpt_cart(:, nkp)/lenconfac, '|'
334  END DO
335  WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*'
336  WRITE (stdout, *) ' '
337 
338  END SUBROUTINE w90_param_write
339 
340 ! **************************************************************************************************
341 !> \brief ...
342 ! **************************************************************************************************
343  SUBROUTINE w90_kmesh_get()
344 
345  ! Variables that are private
346 
347  REAL(kind=dp), PARAMETER :: eta = 99999999.0_dp
348 
349  INTEGER :: counter, i, ifound, j, l, loop, loop_b, loop_s, m, multi(search_shells), n, na, &
350  nap, ndnn, ndnntot, ndnnx, nkp, nkp2, nlist, nn, nnsh, nnshell(num_kpts, search_shells), &
351  nnx, shell
352  LOGICAL :: isneg, ispos
353  REAL(kind=dp) :: bb1, bbn, bk_local(3, num_nnmax, num_kpts), bweight(max_shells), ddelta, &
354  dist, dnn(search_shells), dnn0, dnn1, vkpp(3), vkpp2(3), wb_local(num_nnmax)
355  REAL(kind=dp), ALLOCATABLE :: bvec_tmp(:, :)
356 
357  WRITE (stdout, '(/1x,a)') &
358  '*---------------------------------- K-MESH ----------------------------------*'
359 
360  ! Sort the cell neighbours so we loop in order of distance from the home shell
361  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  dnn0 = 0.0_dp
367  dnn1 = eta
368  ndnntot = 0
369  DO nlist = 1, search_shells
370  DO nkp = 1, num_kpts
371  DO loop = 1, (2*nsupcell + 1)**3
372  l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
373  !
374  vkpp = kpt_cart(:, nkp) + matmul(lmn(:, loop), recip_lattice)
375  dist = sqrt((kpt_cart(1, 1) - vkpp(1))**2 &
376  + (kpt_cart(2, 1) - vkpp(2))**2 + (kpt_cart(3, 1) - vkpp(3))**2)
377  !
378  IF ((dist .GT. kmesh_tol) .AND. (dist .GT. dnn0 + kmesh_tol)) THEN
379  IF (dist .LT. dnn1 - kmesh_tol) THEN
380  dnn1 = dist ! found a closer shell
381  counter = 0
382  END IF
383  IF (dist .GT. (dnn1 - kmesh_tol) .AND. dist .LT. (dnn1 + kmesh_tol)) THEN
384  counter = counter + 1 ! count the multiplicity of the shell
385  END IF
386  END IF
387  END DO
388  END DO
389  IF (dnn1 .LT. eta - kmesh_tol) ndnntot = ndnntot + 1
390  dnn(nlist) = dnn1
391  multi(nlist) = counter
392  dnn0 = dnn1
393  dnn1 = eta
394  END DO
395 
396  WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
397  WRITE (stdout, '(1x,a)') '| Distance to Nearest-Neighbour Shells |'
398  WRITE (stdout, '(1x,a)') '| ------------------------------------ |'
399  WRITE (stdout, '(1x,a)') '| Shell Distance (Bohr^-1) Multiplicity |'
400  WRITE (stdout, '(1x,a)') '| ----- ------------------ ------------ |'
401  DO ndnn = 1, ndnntot
402  WRITE (stdout, '(1x,a,11x,i3,17x,f10.6,19x,i4,12x,a)') '|', ndnn, dnn(ndnn)/lenconfac, multi(ndnn), '|'
403  END DO
404  WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
405 
406  ! Get the shell weights to satisfy the B1 condition
407  CALL kmesh_shell_automatic(multi, dnn, bweight)
408 
409  WRITE (stdout, '(1x,a)', advance='no') '| The following shells are used: '
410  DO ndnn = 1, num_shells
411  IF (ndnn .EQ. num_shells) THEN
412  WRITE (stdout, '(i3,1x)', advance='no') shell_list(ndnn)
413  ELSE
414  WRITE (stdout, '(i3,",")', advance='no') shell_list(ndnn)
415  END IF
416  END DO
417  DO l = 1, 11 - num_shells
418  WRITE (stdout, '(4x)', advance='no')
419  END DO
420  WRITE (stdout, '("|")')
421 
422  nntot = 0
423  DO loop_s = 1, num_shells
424  nntot = nntot + multi(shell_list(loop_s))
425  END DO
426  IF (nntot > num_nnmax) THEN
427  WRITE (stdout, '(a,i2,a)') ' **WARNING: kmesh has found >', num_nnmax, ' nearest neighbours**'
428  WRITE (stdout, '(a)') ' '
429  WRITE (stdout, '(a)') ' This is probably caused by an error in your unit cell specification'
430  WRITE (stdout, '(a)') ' '
431 
432  ALLOCATE (bvec_tmp(3, maxval(multi)))
433  bvec_tmp = 0.0_dp
434  counter = 0
435  DO shell = 1, search_shells
436  CALL kmesh_get_bvectors(multi(shell), 1, dnn(shell), bvec_tmp(:, 1:multi(shell)))
437  DO loop = 1, multi(shell)
438  counter = counter + 1
439  WRITE (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector ', counter, ': (', &
440  bvec_tmp(:, loop)/lenconfac, ')', dnn(shell)/lenconfac, ' |'
441  END DO
442  END DO
443  WRITE (stdout, '(a)') ' '
444  DEALLOCATE (bvec_tmp)
445  cpabort('kmesh_get: something wrong, found too many nearest neighbours')
446  END IF
447 
448  ALLOCATE (nnlist(num_kpts, nntot))
449  ALLOCATE (neigh(num_kpts, nntot/2))
450  ALLOCATE (nncell(3, num_kpts, nntot))
451 
452  ALLOCATE (wb(nntot))
453  ALLOCATE (bka(3, nntot/2))
454  ALLOCATE (bk(3, nntot, num_kpts))
455 
456  nnx = 0
457  DO loop_s = 1, num_shells
458  DO loop_b = 1, multi(shell_list(loop_s))
459  nnx = nnx + 1
460  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  WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
472  WRITE (stdout, '(1x,a)') '| Shell # Nearest-Neighbours |'
473  WRITE (stdout, '(1x,a)') '| ----- -------------------- |'
474  !
475  ! Standard routine
476  !
477  nnshell = 0
478  DO nkp = 1, num_kpts
479  nnx = 0
480  ok: DO ndnnx = 1, num_shells
481  ndnn = shell_list(ndnnx)
482  DO loop = 1, (2*nsupcell + 1)**3
483  l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
484  vkpp2 = matmul(lmn(:, loop), recip_lattice)
485  DO nkp2 = 1, num_kpts
486  vkpp = vkpp2 + kpt_cart(:, nkp2)
487  dist = sqrt((kpt_cart(1, nkp) - vkpp(1))**2 &
488  + (kpt_cart(2, nkp) - vkpp(2))**2 + (kpt_cart(3, nkp) - vkpp(3))**2)
489  IF ((dist .GE. dnn(ndnn)*(1 - kmesh_tol)) .AND. (dist .LE. dnn(ndnn)*(1 + kmesh_tol))) THEN
490  nnx = nnx + 1
491  nnshell(nkp, ndnn) = nnshell(nkp, ndnn) + 1
492  nnlist(nkp, nnx) = nkp2
493  nncell(1, nkp, nnx) = l
494  nncell(2, nkp, nnx) = m
495  nncell(3, nkp, nnx) = n
496  bk_local(:, nnx, nkp) = vkpp(:) - kpt_cart(:, nkp)
497  END IF
498  !if we have the right number of neighbours we can exit
499  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  DO ndnnx = 1, num_shells
508  ndnn = shell_list(ndnnx)
509  WRITE (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|'
510  END DO
511  WRITE (stdout, '(1x,"+",76("-"),"+")')
512 
513  DO nkp = 1, num_kpts
514  nnx = 0
515  DO ndnnx = 1, num_shells
516  ndnn = shell_list(ndnnx)
517  DO nnsh = 1, nnshell(nkp, ndnn)
518  bb1 = 0.0_dp
519  bbn = 0.0_dp
520  nnx = nnx + 1
521  DO i = 1, 3
522  bb1 = bb1 + bk_local(i, nnx, 1)*bk_local(i, nnx, 1)
523  bbn = bbn + bk_local(i, nnx, nkp)*bk_local(i, nnx, nkp)
524  END DO
525  IF (abs(sqrt(bb1) - sqrt(bbn)) .GT. kmesh_tol) THEN
526  WRITE (stdout, '(1x,2f10.6)') bb1, bbn
527  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  DO nkp = 1, num_kpts
538  DO i = 1, 3
539  DO j = 1, 3
540  ddelta = 0.0_dp
541  nnx = 0
542  DO ndnnx = 1, num_shells
543  ndnn = shell_list(ndnnx)
544  DO nnsh = 1, nnshell(1, ndnn)
545  nnx = nnx + 1
546  ddelta = ddelta + wb_local(nnx)*bk_local(i, nnx, nkp)*bk_local(j, nnx, nkp)
547  END DO
548  END DO
549  IF ((i .EQ. j) .AND. (abs(ddelta - 1.0_dp) .GT. kmesh_tol)) THEN
550  WRITE (stdout, '(1x,2i3,f12.8)') i, j, ddelta
551  cpabort('Eq. (B1) not satisfied in kmesh_get (1)')
552  END IF
553  IF ((i .NE. j) .AND. (abs(ddelta) .GT. kmesh_tol)) THEN
554  WRITE (stdout, '(1x,2i3,f12.8)') i, j, ddelta
555  cpabort('Eq. (B1) not satisfied in kmesh_get (2)')
556  END IF
557  END DO
558  END DO
559  END DO
560 
561  WRITE (stdout, '(1x,a)') '| Completeness relation is fully satisfied [Eq. (B1), PRB 56, 12847 (1997)] |'
562  WRITE (stdout, '(1x,"+",76("-"),"+")')
563 
564  !
565  wbtot = 0.0_dp
566  nnx = 0
567  DO ndnnx = 1, num_shells
568  ndnn = shell_list(ndnnx)
569  DO nnsh = 1, nnshell(1, ndnn)
570  nnx = nnx + 1
571  wbtot = wbtot + wb_local(nnx)
572  END DO
573  END DO
574 
575  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  na = 0
579  DO nn = 1, nntot
580  ifound = 0
581  IF (na .NE. 0) THEN
582  DO nap = 1, na
583  CALL utility_compar(bka(1, nap), bk_local(1, nn, 1), ispos, isneg)
584  IF (isneg) ifound = 1
585  END DO
586  END IF
587  IF (ifound .EQ. 0) THEN
588  ! found new vector to add to set
589  na = na + 1
590  bka(1, na) = bk_local(1, nn, 1)
591  bka(2, na) = bk_local(2, nn, 1)
592  bka(3, na) = bk_local(3, nn, 1)
593  END IF
594  END DO
595  IF (na .NE. nnh) cpabort('Did not find right number of bk directions')
596 
597  WRITE (stdout, '(1x,a)') '| b_k Vectors (Bohr^-1) and Weights (Bohr^2) |'
598  WRITE (stdout, '(1x,a)') '| ------------------------------------------ |'
599  WRITE (stdout, '(1x,a)') '| No. b_k(x) b_k(y) b_k(z) w_b |'
600  WRITE (stdout, '(1x,a)') '| --- -------------------------------- -------- |'
601  DO i = 1, nntot
602  WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
603  i, (bk_local(j, i, 1)/lenconfac, j=1, 3), wb_local(i)*lenconfac**2
604  END DO
605  WRITE (stdout, '(1x,"+",76("-"),"+")')
606  WRITE (stdout, '(1x,a)') '| b_k Directions (Bohr^-1) |'
607  WRITE (stdout, '(1x,a)') '| ------------------------ |'
608  WRITE (stdout, '(1x,a)') '| No. x y z |'
609  WRITE (stdout, '(1x,a)') '| --- -------------------------------- |'
610  DO i = 1, nnh
611  WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,21x,"|")') i, (bka(j, i)/lenconfac, j=1, 3)
612  END DO
613  WRITE (stdout, '(1x,"+",76("-"),"+")')
614  WRITE (stdout, *) ' '
615 
616  ! find index array
617  DO nkp = 1, num_kpts
618  DO na = 1, nnh
619  ! first, zero the index array so we can check it gets filled
620  neigh(nkp, na) = 0
621  ! now search through list of neighbours of this k-point
622  DO nn = 1, nntot
623  CALL utility_compar(bka(1, na), bk_local(1, nn, nkp), ispos, isneg)
624  IF (ispos) neigh(nkp, na) = nn
625  END DO
626  ! check found
627  IF (neigh(nkp, na) .EQ. 0) THEN
628  WRITE (stdout, *) ' nkp,na=', nkp, na
629  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  DO loop = 1, nntot
636  wb(loop) = wb_local(loop)
637  END DO
638 
639  DO loop_s = 1, num_kpts
640  DO loop = 1, nntot
641  bk(:, loop, loop_s) = bk_local(:, loop, loop_s)
642  END DO
643  END DO
644 
645  END SUBROUTINE w90_kmesh_get
646 
647 ! **************************************************************************************************
648 !> \brief ...
649 ! **************************************************************************************************
650  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  counter = 1
666  lmn(:, counter) = 0
667  dist(counter) = 0.0_dp
668  DO l = -nsupcell, nsupcell
669  DO m = -nsupcell, nsupcell
670  DO n = -nsupcell, nsupcell
671  IF (l == 0 .AND. m == 0 .AND. n == 0) cycle
672  counter = counter + 1
673  lmn(1, counter) = l; lmn(2, counter) = m; lmn(3, counter) = n
674  pos = matmul(lmn(:, counter), recip_lattice)
675  dist(counter) = sqrt(dot_product(pos, pos))
676  END DO
677  END DO
678  END DO
679 
680  DO loop = (2*nsupcell + 1)**3, 1, -1
681  indx = internal_maxloc(dist)
682  dist_cp(loop) = dist(indx(1))
683  lmn_cp(:, loop) = lmn(:, indx(1))
684  dist(indx(1)) = -1.0_dp
685  END DO
686 
687  lmn = lmn_cp
688  dist = dist_cp
689 
690  END SUBROUTINE w90_kmesh_supercell_sort
691 
692 ! **************************************************************************************************
693 !> \brief ...
694 !> \param dist ...
695 !> \return ...
696 ! **************************************************************************************************
697  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  list = 0
711  counter = 1
712 
713  guess = maxloc(dist)
714  list(1) = guess(1)
715  ! look for any degenerate values
716  DO loop = 1, (2*nsupcell + 1)**3
717  IF (loop == guess(1)) cycle
718  IF (abs(dist(loop) - dist(guess(1))) < eps8) THEN
719  counter = counter + 1
720  list(counter) = loop
721  END IF
722  END DO
723  ! and always return the lowest index
724  internal_maxloc = minval(list(1:counter))
725 
726  END FUNCTION internal_maxloc
727 
728 ! **************************************************************************************************
729 !> \brief ...
730 !> \param multi ...
731 !> \param dnn ...
732 !> \param bweight ...
733 ! **************************************************************************************************
734  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  REAL(kind=dp), ALLOCATABLE :: bvector(:, :, :)
757  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: singv
758  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, smat, umat, vmat
759 
760  ALLOCATE (bvector(3, maxval(multi), max_shells))
761  bvector = 0.0_dp; bweight = 0.0_dp
762 
763  WRITE (stdout, '(1x,a)') '| The b-vectors are chosen automatically |'
764 
765  b1sat = .false.
766  DO shell = 1, search_shells
767  cur_shell = num_shells + 1
768 
769  ! get the b vectors for the new shell
770  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  lpar = .false.
774  IF (num_shells > 0) THEN
775  DO loop_bn = 1, multi(shell)
776  DO loop_s = 1, num_shells
777  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  dot_product(bvector(:, loop_b, loop_s), bvector(:, loop_b, loop_s)))
781  IF (abs(abs(delta) - 1.0_dp) < eps6) lpar = .true.
782  END DO
783  END DO
784  END DO
785  END IF
786 
787  IF (lpar) THEN
788  IF (iprint >= 3) THEN
789  WRITE (stdout, '(1x,a)') '| This shell is linearly dependent on existing shells: Trying next shell |'
790  END IF
791  cycle
792  END IF
793 
794  num_shells = num_shells + 1
795  shell_list(num_shells) = shell
796 
797  ALLOCATE (amat(max_shells, num_shells))
798  ALLOCATE (umat(max_shells, max_shells))
799  ALLOCATE (vmat(num_shells, num_shells))
800  ALLOCATE (smat(num_shells, max_shells))
801  ALLOCATE (singv(num_shells))
802  amat = 0.0_dp; umat = 0.0_dp; vmat = 0.0_dp; smat = 0.0_dp; singv = 0.0_dp
803 
804  amat = 0.0_dp
805  DO loop_s = 1, num_shells
806  DO loop_b = 1, multi(shell_list(loop_s))
807  amat(1, loop_s) = amat(1, loop_s) + bvector(1, loop_b, loop_s)*bvector(1, loop_b, loop_s)
808  amat(2, loop_s) = amat(2, loop_s) + bvector(2, loop_b, loop_s)*bvector(2, loop_b, loop_s)
809  amat(3, loop_s) = amat(3, loop_s) + bvector(3, loop_b, loop_s)*bvector(3, loop_b, loop_s)
810  amat(4, loop_s) = amat(4, loop_s) + bvector(1, loop_b, loop_s)*bvector(2, loop_b, loop_s)
811  amat(5, loop_s) = amat(5, loop_s) + bvector(2, loop_b, loop_s)*bvector(3, loop_b, loop_s)
812  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  info = 0
817  CALL dgesvd('A', 'A', max_shells, num_shells, amat, max_shells, singv, umat, &
818  max_shells, vmat, num_shells, work, lwork, info)
819  IF (info < 0) THEN
820  WRITE (stdout, '(1x,a,1x,I1,1x,a)') 'kmesh_shell_automatic: Argument', abs(info), 'of dgesvd is incorrect'
821  cpabort('kmesh_shell_automatic: Problem with Singular Value Decomposition')
822  ELSE IF (info > 0) THEN
823  cpabort('kmesh_shell_automatic: Singular Value Decomposition did not converge')
824  END IF
825 
826  IF (any(abs(singv) < eps5)) THEN
827  IF (num_shells == 1) THEN
828  CALL cp_abort(__location__, "kmesh_shell_automatic: "// &
829  "Singular Value Decomposition has found a very small singular value.")
830  ELSE
831  WRITE (stdout, '(1x,a)') '| SVD found small singular value, Rejecting this shell and trying the next |'
832  b1sat = .false.
833  num_shells = num_shells - 1
834  DEALLOCATE (amat, umat, vmat, smat, singv)
835  cycle
836  END IF
837  END IF
838 
839  smat = 0.0_dp
840  DO loop_s = 1, num_shells
841  smat(loop_s, loop_s) = 1/singv(loop_s)
842  END DO
843 
844  bweight(1:num_shells) = matmul(transpose(vmat), matmul(smat, matmul(transpose(umat), TARGET)))
845  IF (iprint >= 2) THEN
846  DO loop_s = 1, num_shells
847  WRITE (stdout, '(1x,a,I2,a,f12.7,5x,a8,36x,a)') '| Shell: ', loop_s, &
848  ' 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  DO loop_i = 1, 3
855  DO loop_j = loop_i, 3
856  delta = 0.0_dp
857  DO loop_s = 1, num_shells
858  DO loop_b = 1, multi(shell_list(loop_s))
859  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  IF (loop_i == loop_j) THEN
863  IF (abs(delta - 1.0_dp) > kmesh_tol) b1sat = .false.
864  END IF
865  IF (loop_i /= loop_j) THEN
866  IF (abs(delta) > kmesh_tol) b1sat = .false.
867  END IF
868  END DO
869  END DO
870 
871  IF (.NOT. b1sat) THEN
872  IF (shell < search_shells .AND. iprint >= 3) THEN
873  WRITE (stdout, '(1x,a,24x,a1)') '| B1 condition is not satisfied: Adding another shell', '|'
874  ELSEIF (shell == search_shells) THEN
875  WRITE (stdout, *) ' '
876  WRITE (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells'
877  WRITE (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid'
878  WRITE (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=12)'
879  WRITE (stdout, *) ' '
880  cpabort('kmesh_get_automatic')
881  END IF
882  END IF
883 
884  DEALLOCATE (amat, umat, vmat, smat, singv)
885 
886  IF (b1sat) EXIT
887 
888  END DO
889 
890  IF (.NOT. b1sat) THEN
891  WRITE (stdout, *) ' '
892  WRITE (stdout, '(1x,a,i3,a)') 'Unable to satisfy B1 with any of the first ', search_shells, ' shells'
893  WRITE (stdout, '(1x,a)') 'Your cell might be very long, or you may have an irregular MP grid'
894  WRITE (stdout, '(1x,a)') 'Try increasing the parameter search_shells in the win file (default=12)'
895  WRITE (stdout, *) ' '
896  cpabort('kmesh_get_automatic')
897  END IF
898 
899  END SUBROUTINE kmesh_shell_automatic
900 
901 ! **************************************************************************************************
902 !> \brief ...
903 !> \param multi ...
904 !> \param kpt ...
905 !> \param shell_dist ...
906 !> \param bvector ...
907 ! **************************************************************************************************
908  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  bvector = 0.0_dp
923 
924  num_bvec = 0
925  ok: DO loop = 1, (2*nsupcell + 1)**3
926  vkpp2 = matmul(lmn(:, loop), recip_lattice)
927  DO nkp2 = 1, num_kpts
928  vkpp = vkpp2 + kpt_cart(:, nkp2)
929  dist = sqrt((kpt_cart(1, kpt) - vkpp(1))**2 &
930  + (kpt_cart(2, kpt) - vkpp(2))**2 + (kpt_cart(3, kpt) - vkpp(3))**2)
931  IF ((dist .GE. shell_dist*(1.0_dp - kmesh_tol)) .AND. dist .LE. shell_dist*(1.0_dp + kmesh_tol)) THEN
932  num_bvec = num_bvec + 1
933  bvector(:, num_bvec) = vkpp(:) - kpt_cart(:, kpt)
934  END IF
935  !if we have the right number of neighbours we can exit
936  IF (num_bvec == multi) cycle ok
937  END DO
938  END DO ok
939 
940  IF (num_bvec < multi) cpabort('kmesh_get_bvector: Not enough bvectors found')
941 
942  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  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  ispos = sum((a - b)**2) .LT. eps8
956  isneg = sum((a + b)**2) .LT. eps8
957  END SUBROUTINE utility_compar
958 
959 ! **************************************************************************************************
960 
961 END MODULE wannier90
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
wannier90
Outtakes from Wannier90 code.
Definition: wannier90.F:66
wannier90::wannier_setup
subroutine, public wannier_setup(mp_grid_loc, num_kpts_loc, real_lattice_loc, recip_lattice_loc, kpt_latt_loc, nntot_loc, nnlist_loc, nncell_loc, iounit)
...
Definition: wannier90.F:140
wannier90::w90_write_header
subroutine, public w90_write_header(stdout)
...
Definition: wannier90.F:209