(git:374b731)
Loading...
Searching...
No Matches
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
66MODULE 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
123CONTAINS
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
961END 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