Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 16 : SUBROUTINE wannier_setup(mp_grid_loc, num_kpts_loc, &
138 16 : real_lattice_loc, recip_lattice_loc, kpt_latt_loc, &
139 16 : 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 16 : length_unit = 'bohr'
157 16 : lenconfac = 1.0_dp/bohr
158 16 : stdout = iounit
159 :
160 16 : CALL w90_write_header(stdout)
161 :
162 16 : WRITE (stdout, '(a/)') ' Setting up k-point neighbours...'
163 :
164 : ! copy local data into module variables
165 16 : mp_grid = mp_grid_loc
166 16 : num_kpts = num_kpts_loc
167 16 : real_lattice = real_lattice_loc
168 16 : recip_lattice = recip_lattice_loc
169 48 : ALLOCATE (kpt_latt(3, num_kpts))
170 32 : ALLOCATE (kpt_cart(3, num_kpts))
171 976 : kpt_latt(1:3, 1:num_kpts) = kpt_latt_loc(1:3, 1:num_kpts)
172 256 : DO nkp = 1, num_kpts
173 6256 : kpt_cart(:, nkp) = MATMUL(kpt_latt(:, nkp), recip_lattice(:, :))
174 : END DO
175 :
176 16 : num_shells = 0
177 16 : 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 16 : real_lattice(1, 3)*(real_lattice(2, 1)*real_lattice(3, 2) - real_lattice(3, 1)*real_lattice(2, 2))
182 16 : iprint = 1
183 16 : search_shells = 12
184 16 : kmesh_tol = 0.000001_dp
185 16 : num_exclude_bands = 0
186 :
187 16 : CALL w90_param_write(stdout)
188 :
189 16 : CALL w90_kmesh_get()
190 :
191 16 : nntot_loc = nntot
192 3088 : nnlist_loc = 0
193 1552 : nnlist_loc(:, 1:nntot) = nnlist(:, 1:nntot)
194 11728 : nncell_loc = 0
195 5872 : nncell_loc(:, :, 1:nntot) = nncell(:, :, 1:nntot)
196 :
197 16 : DEALLOCATE (bk, bka, wb)
198 16 : DEALLOCATE (nncell, neigh, nnlist)
199 16 : DEALLOCATE (kpt_latt, kpt_cart, shell_list)
200 :
201 16 : WRITE (stdout, '(/a/)') ' Finished setting up k-point neighbours.'
202 :
203 16 : END SUBROUTINE wannier_setup
204 : ! **************************************************************************************************
205 : !> \brief ...
206 : !> \param stdout ...
207 : ! **************************************************************************************************
208 16 : SUBROUTINE w90_write_header(stdout)
209 : INTEGER, INTENT(IN) :: stdout
210 :
211 16 : WRITE (stdout, *)
212 16 : WRITE (stdout, *) ' +---------------------------------------------------+'
213 16 : WRITE (stdout, *) ' | |'
214 16 : WRITE (stdout, *) ' | WANNIER90 |'
215 16 : WRITE (stdout, *) ' | |'
216 16 : WRITE (stdout, *) ' +---------------------------------------------------+'
217 16 : WRITE (stdout, *) ' | |'
218 16 : WRITE (stdout, *) ' | Welcome to the Maximally-Localized |'
219 16 : WRITE (stdout, *) ' | Generalized Wannier Functions code |'
220 16 : WRITE (stdout, *) ' | http://www.wannier.org |'
221 16 : WRITE (stdout, *) ' | |'
222 16 : WRITE (stdout, *) ' | Wannier90 v2.0 Authors: |'
223 16 : WRITE (stdout, *) ' | Arash A. Mostofi (Imperial College London) |'
224 16 : WRITE (stdout, *) ' | Giovanni Pizzi (EPFL) |'
225 16 : WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |'
226 16 : WRITE (stdout, *) ' | Jonathan R. Yates (University of Oxford) |'
227 16 : WRITE (stdout, *) ' | |'
228 16 : WRITE (stdout, *) ' | Wannier90 Contributors: |'
229 16 : WRITE (stdout, *) ' | Young-Su Lee (KIST, S. Korea) |'
230 16 : WRITE (stdout, *) ' | Matthew Shelley (Imperial College London) |'
231 16 : WRITE (stdout, *) ' | Nicolas Poilvert (Penn State University) |'
232 16 : WRITE (stdout, *) ' | Raffaello Bianco (Paris 6 and CNRS) |'
233 16 : WRITE (stdout, *) ' | Gabriele Sclauzero (ETH Zurich) |'
234 16 : WRITE (stdout, *) ' | |'
235 16 : WRITE (stdout, *) ' | Wannier77 Authors: |'
236 16 : WRITE (stdout, *) ' | Nicola Marzari (EPFL) |'
237 16 : WRITE (stdout, *) ' | Ivo Souza (Universidad del Pais Vasco) |'
238 16 : WRITE (stdout, *) ' | David Vanderbilt (Rutgers University) |'
239 16 : WRITE (stdout, *) ' | |'
240 16 : WRITE (stdout, *) ' | Please cite |'
241 16 : WRITE (stdout, *) ' | |'
242 16 : WRITE (stdout, *) ' | [ref] "Wannier90: A Tool for Obtaining Maximally |'
243 16 : WRITE (stdout, *) ' | Localised Wannier Functions" |'
244 16 : WRITE (stdout, *) ' | A. A. Mostofi, J. R. Yates, Y.-S. Lee, |'
245 16 : WRITE (stdout, *) ' | I. Souza, D. Vanderbilt and N. Marzari |'
246 16 : WRITE (stdout, *) ' | Comput. Phys. Commun. 178, 685 (2008) |'
247 16 : WRITE (stdout, *) ' | |'
248 16 : WRITE (stdout, *) ' | in any publications arising from the use of |'
249 16 : WRITE (stdout, *) ' | this code. For the method please cite |'
250 16 : WRITE (stdout, *) ' | |'
251 16 : WRITE (stdout, *) ' | [ref] "Maximally Localized Generalised Wannier |'
252 16 : WRITE (stdout, *) ' | Functions for Composite Energy Bands" |'
253 16 : WRITE (stdout, *) ' | N. Marzari and D. Vanderbilt |'
254 16 : WRITE (stdout, *) ' | Phys. Rev. B 56 12847 (1997) |'
255 16 : WRITE (stdout, *) ' | |'
256 16 : WRITE (stdout, *) ' | [ref] "Maximally Localized Wannier Functions |'
257 16 : WRITE (stdout, *) ' | for Entangled Energy Bands" |'
258 16 : WRITE (stdout, *) ' | I. Souza, N. Marzari and D. Vanderbilt |'
259 16 : WRITE (stdout, *) ' | Phys. Rev. B 65 035109 (2001) |'
260 16 : WRITE (stdout, *) ' | |'
261 16 : WRITE (stdout, *) ' | |'
262 16 : WRITE (stdout, *) ' | Copyright (c) 1996-2015 |'
263 16 : WRITE (stdout, *) ' | Arash A. Mostofi, Jonathan R. Yates, |'
264 16 : WRITE (stdout, *) ' | Young-Su Lee, Giovanni Pizzi, Ivo Souza, |'
265 16 : WRITE (stdout, *) ' | David Vanderbilt and Nicola Marzari |'
266 16 : WRITE (stdout, *) ' | |'
267 16 : WRITE (stdout, *) ' | Release: 2.0.1 2nd April 2015 |'
268 16 : WRITE (stdout, *) ' | |'
269 16 : WRITE (stdout, *) ' | This program is free software; you can |'
270 16 : WRITE (stdout, *) ' | redistribute it and/or modify it under the terms |'
271 16 : WRITE (stdout, *) ' | of the GNU General Public License as published by |'
272 16 : WRITE (stdout, *) ' | the Free Software Foundation; either version 2 of |'
273 16 : WRITE (stdout, *) ' | the License, or (at your option) any later version|'
274 16 : WRITE (stdout, *) ' | |'
275 16 : WRITE (stdout, *) ' | This program is distributed in the hope that it |'
276 16 : WRITE (stdout, *) ' | will be useful, but WITHOUT ANY WARRANTY; without |'
277 16 : WRITE (stdout, *) ' | even the implied warranty of MERCHANTABILITY or |'
278 16 : WRITE (stdout, *) ' | FITNESS FOR A PARTICULAR PURPOSE. See the GNU |'
279 16 : WRITE (stdout, *) ' | General Public License for more details. |'
280 16 : WRITE (stdout, *) ' | |'
281 16 : WRITE (stdout, *) ' | You should have received a copy of the GNU General|'
282 16 : WRITE (stdout, *) ' | Public License along with this program; if not, |'
283 16 : WRITE (stdout, *) ' | write to the Free Software Foundation, Inc., |'
284 16 : WRITE (stdout, *) ' | 675 Mass Ave, Cambridge, MA 02139, USA. |'
285 16 : WRITE (stdout, *) ' | |'
286 16 : WRITE (stdout, *) ' +---------------------------------------------------+'
287 16 : WRITE (stdout, *) ''
288 :
289 16 : END SUBROUTINE w90_write_header
290 :
291 : ! **************************************************************************************************
292 : !> \brief ...
293 : !> \param stdout ...
294 : ! **************************************************************************************************
295 16 : SUBROUTINE w90_param_write(stdout)
296 : INTEGER, INTENT(IN) :: stdout
297 :
298 : INTEGER :: i, nkp
299 :
300 : ! System
301 16 : WRITE (stdout, '(36x,a6)') '------'
302 16 : WRITE (stdout, '(36x,a6)') 'SYSTEM'
303 16 : WRITE (stdout, '(36x,a6)') '------'
304 16 : WRITE (stdout, *)
305 16 : WRITE (stdout, '(28x,a22)') 'Lattice Vectors (Bohr)'
306 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_1', (real_lattice(1, I)*lenconfac, i=1, 3)
307 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_2', (real_lattice(2, I)*lenconfac, i=1, 3)
308 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'a_3', (real_lattice(3, I)*lenconfac, i=1, 3)
309 16 : WRITE (stdout, *)
310 : WRITE (stdout, '(19x,a17,3x,f11.5)', advance='no') &
311 16 : 'Unit Cell Volume:', cell_volume*lenconfac**3
312 16 : WRITE (stdout, '(2x,a8)') '(Bohr^3)'
313 16 : WRITE (stdout, *)
314 16 : WRITE (stdout, '(22x,a34)') 'Reciprocal-Space Vectors (Bohr^-1)'
315 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_1', (recip_lattice(1, I)/lenconfac, i=1, 3)
316 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_2', (recip_lattice(2, I)/lenconfac, i=1, 3)
317 64 : WRITE (stdout, '(20x,a3,2x,3F11.6)') 'b_3', (recip_lattice(3, I)/lenconfac, i=1, 3)
318 16 : WRITE (stdout, *) ' '
319 16 : WRITE (stdout, *) ' '
320 : ! K-points
321 16 : WRITE (stdout, '(32x,a)') '------------'
322 16 : WRITE (stdout, '(32x,a)') 'K-POINT GRID'
323 16 : WRITE (stdout, '(32x,a)') '------------'
324 16 : WRITE (stdout, *) ' '
325 16 : 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 32 : 'Total points =', num_kpts
327 16 : WRITE (stdout, *) ' '
328 16 : WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*'
329 16 : WRITE (stdout, '(1x,a)') '| k-point Fractional Coordinate Cartesian Coordinate (Bohr^-1) |'
330 16 : WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
331 256 : DO nkp = 1, num_kpts
332 240 : WRITE (stdout, '(1x,a1,i6,1x,3F10.5,3x,a1,1x,3F10.5,4x,a1)') '|', &
333 1216 : nkp, kpt_latt(:, nkp), '|', kpt_cart(:, nkp)/lenconfac, '|'
334 : END DO
335 16 : WRITE (stdout, '(1x,a)') '*----------------------------------------------------------------------------*'
336 16 : WRITE (stdout, *) ' '
337 :
338 16 : END SUBROUTINE w90_param_write
339 :
340 : ! **************************************************************************************************
341 : !> \brief ...
342 : ! **************************************************************************************************
343 16 : SUBROUTINE w90_kmesh_get()
344 :
345 : ! Variables that are private
346 :
347 : REAL(kind=dp), PARAMETER :: eta = 99999999.0_dp
348 :
349 32 : INTEGER :: counter, i, ifound, j, l, loop, loop_b, loop_s, m, multi(search_shells), n, na, &
350 32 : nap, ndnn, ndnntot, ndnnx, nkp, nkp2, nlist, nn, nnsh, nnshell(num_kpts, search_shells), &
351 : nnx, shell
352 : LOGICAL :: isneg, ispos
353 32 : REAL(kind=dp) :: bb1, bbn, bk_local(3, num_nnmax, num_kpts), bweight(max_shells), ddelta, &
354 32 : dist, dnn(search_shells), dnn0, dnn1, vkpp(3), vkpp2(3), wb_local(num_nnmax)
355 16 : REAL(kind=dp), ALLOCATABLE :: bvec_tmp(:, :)
356 :
357 : WRITE (stdout, '(/1x,a)') &
358 16 : '*---------------------------------- K-MESH ----------------------------------*'
359 :
360 : ! Sort the cell neighbours so we loop in order of distance from the home shell
361 16 : 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 16 : dnn0 = 0.0_dp
367 16 : dnn1 = eta
368 16 : ndnntot = 0
369 208 : DO nlist = 1, search_shells
370 3072 : DO nkp = 1, num_kpts
371 3836352 : DO loop = 1, (2*nsupcell + 1)**3
372 3833280 : l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
373 : !
374 72832320 : vkpp = kpt_cart(:, nkp) + MATMUL(lmn(:, loop), recip_lattice)
375 : dist = SQRT((kpt_cart(1, 1) - vkpp(1))**2 &
376 3833280 : + (kpt_cart(2, 1) - vkpp(2))**2 + (kpt_cart(3, 1) - vkpp(3))**2)
377 : !
378 3836160 : IF ((dist > kmesh_tol) .AND. (dist > dnn0 + kmesh_tol)) THEN
379 3818272 : IF (dist < dnn1 - kmesh_tol) THEN
380 364 : dnn1 = dist ! found a closer shell
381 364 : counter = 0
382 : END IF
383 3818272 : IF (dist > (dnn1 - kmesh_tol) .AND. dist < (dnn1 + kmesh_tol)) THEN
384 4484 : counter = counter + 1 ! count the multiplicity of the shell
385 : END IF
386 : END IF
387 : END DO
388 : END DO
389 192 : IF (dnn1 < eta - kmesh_tol) ndnntot = ndnntot + 1
390 192 : dnn(nlist) = dnn1
391 192 : multi(nlist) = counter
392 192 : dnn0 = dnn1
393 208 : dnn1 = eta
394 : END DO
395 :
396 16 : WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
397 16 : WRITE (stdout, '(1x,a)') '| Distance to Nearest-Neighbour Shells |'
398 16 : WRITE (stdout, '(1x,a)') '| ------------------------------------ |'
399 16 : WRITE (stdout, '(1x,a)') '| Shell Distance (Bohr^-1) Multiplicity |'
400 16 : WRITE (stdout, '(1x,a)') '| ----- ------------------ ------------ |'
401 208 : DO ndnn = 1, ndnntot
402 208 : WRITE (stdout, '(1x,a,11x,i3,17x,f10.6,19x,i4,12x,a)') '|', ndnn, dnn(ndnn)/lenconfac, multi(ndnn), '|'
403 : END DO
404 16 : WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
405 :
406 : ! Get the shell weights to satisfy the B1 condition
407 16 : CALL kmesh_shell_automatic(multi, dnn, bweight)
408 :
409 16 : WRITE (stdout, '(1x,a)', advance='no') '| The following shells are used: '
410 32 : DO ndnn = 1, num_shells
411 32 : IF (ndnn == num_shells) THEN
412 16 : 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 176 : DO l = 1, 11 - num_shells
418 176 : WRITE (stdout, '(4x)', advance='no')
419 : END DO
420 16 : WRITE (stdout, '("|")')
421 :
422 16 : nntot = 0
423 32 : DO loop_s = 1, num_shells
424 32 : nntot = nntot + multi(shell_list(loop_s))
425 : END DO
426 16 : 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 64 : ALLOCATE (nnlist(num_kpts, nntot))
449 64 : ALLOCATE (neigh(num_kpts, nntot/2))
450 64 : ALLOCATE (nncell(3, num_kpts, nntot))
451 :
452 48 : ALLOCATE (wb(nntot))
453 48 : ALLOCATE (bka(3, nntot/2))
454 64 : ALLOCATE (bk(3, nntot, num_kpts))
455 :
456 16 : nnx = 0
457 32 : DO loop_s = 1, num_shells
458 128 : DO loop_b = 1, multi(shell_list(loop_s))
459 96 : nnx = nnx + 1
460 112 : 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 16 : WRITE (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
472 16 : WRITE (stdout, '(1x,a)') '| Shell # Nearest-Neighbours |'
473 16 : WRITE (stdout, '(1x,a)') '| ----- -------------------- |'
474 : !
475 : ! Standard routine
476 : !
477 3088 : nnshell = 0
478 256 : DO nkp = 1, num_kpts
479 240 : nnx = 0
480 496 : ok: DO ndnnx = 1, num_shells
481 240 : ndnn = shell_list(ndnnx)
482 1524 : DO loop = 1, (2*nsupcell + 1)**3
483 1284 : l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
484 20544 : vkpp2 = MATMUL(lmn(:, loop), recip_lattice)
485 42212 : DO nkp2 = 1, num_kpts
486 164672 : vkpp = vkpp2 + kpt_cart(:, nkp2)
487 : dist = SQRT((kpt_cart(1, nkp) - vkpp(1))**2 &
488 41168 : + (kpt_cart(2, nkp) - vkpp(2))**2 + (kpt_cart(3, nkp) - vkpp(3))**2)
489 41168 : IF ((dist >= dnn(ndnn)*(1 - kmesh_tol)) .AND. (dist <= dnn(ndnn)*(1 + kmesh_tol))) THEN
490 1440 : nnx = nnx + 1
491 1440 : nnshell(nkp, ndnn) = nnshell(nkp, ndnn) + 1
492 1440 : nnlist(nkp, nnx) = nkp2
493 1440 : nncell(1, nkp, nnx) = l
494 1440 : nncell(2, nkp, nnx) = m
495 1440 : nncell(3, nkp, nnx) = n
496 5760 : bk_local(:, nnx, nkp) = vkpp(:) - kpt_cart(:, nkp)
497 : END IF
498 : !if we have the right number of neighbours we can exit
499 42212 : 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 32 : DO ndnnx = 1, num_shells
508 16 : ndnn = shell_list(ndnnx)
509 32 : WRITE (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|'
510 : END DO
511 16 : WRITE (stdout, '(1x,"+",76("-"),"+")')
512 :
513 256 : DO nkp = 1, num_kpts
514 240 : nnx = 0
515 496 : DO ndnnx = 1, num_shells
516 240 : ndnn = shell_list(ndnnx)
517 1920 : DO nnsh = 1, nnshell(nkp, ndnn)
518 1440 : bb1 = 0.0_dp
519 1440 : bbn = 0.0_dp
520 1440 : nnx = nnx + 1
521 5760 : DO i = 1, 3
522 4320 : bb1 = bb1 + bk_local(i, nnx, 1)*bk_local(i, nnx, 1)
523 5760 : bbn = bbn + bk_local(i, nnx, nkp)*bk_local(i, nnx, nkp)
524 : END DO
525 1680 : IF (ABS(SQRT(bb1) - SQRT(bbn)) > 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 256 : DO nkp = 1, num_kpts
538 976 : DO i = 1, 3
539 3120 : DO j = 1, 3
540 2160 : ddelta = 0.0_dp
541 2160 : nnx = 0
542 4320 : DO ndnnx = 1, num_shells
543 2160 : ndnn = shell_list(ndnnx)
544 17280 : DO nnsh = 1, nnshell(1, ndnn)
545 12960 : nnx = nnx + 1
546 15120 : ddelta = ddelta + wb_local(nnx)*bk_local(i, nnx, nkp)*bk_local(j, nnx, nkp)
547 : END DO
548 : END DO
549 2160 : IF ((i == j) .AND. (ABS(ddelta - 1.0_dp) > 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 2880 : IF ((i /= j) .AND. (ABS(ddelta) > 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 16 : WRITE (stdout, '(1x,a)') '| Completeness relation is fully satisfied [Eq. (B1), PRB 56, 12847 (1997)] |'
562 16 : WRITE (stdout, '(1x,"+",76("-"),"+")')
563 :
564 : !
565 16 : wbtot = 0.0_dp
566 16 : nnx = 0
567 32 : DO ndnnx = 1, num_shells
568 16 : ndnn = shell_list(ndnnx)
569 128 : DO nnsh = 1, nnshell(1, ndnn)
570 96 : nnx = nnx + 1
571 112 : wbtot = wbtot + wb_local(nnx)
572 : END DO
573 : END DO
574 :
575 16 : 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 16 : na = 0
579 112 : DO nn = 1, nntot
580 96 : ifound = 0
581 96 : IF (na /= 0) THEN
582 272 : DO nap = 1, na
583 192 : CALL utility_compar(bka(1, nap), bk_local(1, nn, 1), ispos, isneg)
584 272 : IF (isneg) ifound = 1
585 : END DO
586 : END IF
587 96 : IF (ifound == 0) THEN
588 : ! found new vector to add to set
589 48 : na = na + 1
590 48 : bka(1, na) = bk_local(1, nn, 1)
591 48 : bka(2, na) = bk_local(2, nn, 1)
592 48 : bka(3, na) = bk_local(3, nn, 1)
593 : END IF
594 : END DO
595 16 : IF (na /= nnh) CPABORT('Did not find right number of bk directions')
596 :
597 16 : WRITE (stdout, '(1x,a)') '| b_k Vectors (Bohr^-1) and Weights (Bohr^2) |'
598 16 : WRITE (stdout, '(1x,a)') '| ------------------------------------------ |'
599 16 : WRITE (stdout, '(1x,a)') '| No. b_k(x) b_k(y) b_k(z) w_b |'
600 16 : WRITE (stdout, '(1x,a)') '| --- -------------------------------- -------- |'
601 112 : DO i = 1, nntot
602 : WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
603 400 : i, (bk_local(j, i, 1)/lenconfac, j=1, 3), wb_local(i)*lenconfac**2
604 : END DO
605 16 : WRITE (stdout, '(1x,"+",76("-"),"+")')
606 16 : WRITE (stdout, '(1x,a)') '| b_k Directions (Bohr^-1) |'
607 16 : WRITE (stdout, '(1x,a)') '| ------------------------ |'
608 16 : WRITE (stdout, '(1x,a)') '| No. x y z |'
609 16 : WRITE (stdout, '(1x,a)') '| --- -------------------------------- |'
610 64 : DO i = 1, nnh
611 208 : WRITE (stdout, '(1x,"|",11x,i3,5x,3f12.6,21x,"|")') i, (bka(j, i)/lenconfac, j=1, 3)
612 : END DO
613 16 : WRITE (stdout, '(1x,"+",76("-"),"+")')
614 16 : WRITE (stdout, *) ' '
615 :
616 : ! find index array
617 256 : DO nkp = 1, num_kpts
618 976 : DO na = 1, nnh
619 : ! first, zero the index array so we can check it gets filled
620 720 : neigh(nkp, na) = 0
621 : ! now search through list of neighbours of this k-point
622 5040 : DO nn = 1, nntot
623 4320 : CALL utility_compar(bka(1, na), bk_local(1, nn, nkp), ispos, isneg)
624 5040 : IF (ispos) neigh(nkp, na) = nn
625 : END DO
626 : ! check found
627 960 : IF (neigh(nkp, na) == 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 112 : DO loop = 1, nntot
636 112 : wb(loop) = wb_local(loop)
637 : END DO
638 :
639 256 : DO loop_s = 1, num_kpts
640 1696 : DO loop = 1, nntot
641 6000 : bk(:, loop, loop_s) = bk_local(:, loop, loop_s)
642 : END DO
643 : END DO
644 :
645 16 : END SUBROUTINE w90_kmesh_get
646 :
647 : ! **************************************************************************************************
648 : !> \brief ...
649 : ! **************************************************************************************************
650 16 : 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 16 : counter = 1
666 64 : lmn(:, counter) = 0
667 16 : dist(counter) = 0.0_dp
668 192 : DO l = -nsupcell, nsupcell
669 2128 : DO m = -nsupcell, nsupcell
670 23408 : DO n = -nsupcell, nsupcell
671 21296 : IF (l == 0 .AND. m == 0 .AND. n == 0) CYCLE
672 21280 : counter = counter + 1
673 21280 : lmn(1, counter) = l; lmn(2, counter) = m; lmn(3, counter) = n
674 340480 : pos = MATMUL(lmn(:, counter), recip_lattice)
675 87056 : dist(counter) = SQRT(DOT_PRODUCT(pos, pos))
676 : END DO
677 : END DO
678 : END DO
679 :
680 21312 : DO loop = (2*nsupcell + 1)**3, 1, -1
681 42592 : indx = internal_maxloc(dist)
682 21296 : dist_cp(loop) = dist(indx(1))
683 85184 : lmn_cp(:, loop) = lmn(:, indx(1))
684 21312 : dist(indx(1)) = -1.0_dp
685 : END DO
686 :
687 16 : lmn = lmn_cp
688 : dist = dist_cp
689 :
690 16 : END SUBROUTINE w90_kmesh_supercell_sort
691 :
692 : ! **************************************************************************************************
693 : !> \brief ...
694 : !> \param dist ...
695 : !> \return ...
696 : ! **************************************************************************************************
697 21296 : 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 21296 : list = 0
711 21296 : counter = 1
712 :
713 28408864 : guess = MAXLOC(dist)
714 21296 : list(1) = guess(1)
715 : ! look for any degenerate values
716 28366272 : DO loop = 1, (2*nsupcell + 1)**3
717 28344976 : IF (loop == guess(1)) CYCLE
718 28344976 : IF (ABS(dist(loop) - dist(guess(1))) < eps8) THEN
719 416912 : counter = counter + 1
720 416912 : list(counter) = loop
721 : END IF
722 : END DO
723 : ! and always return the lowest index
724 459504 : internal_maxloc = MINVAL(list(1:counter))
725 :
726 21296 : END FUNCTION internal_maxloc
727 :
728 : ! **************************************************************************************************
729 : !> \brief ...
730 : !> \param multi ...
731 : !> \param dnn ...
732 : !> \param bweight ...
733 : ! **************************************************************************************************
734 16 : 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 16 : REAL(kind=dp), ALLOCATABLE :: bvector(:, :, :)
757 16 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: singv
758 16 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, smat, umat, vmat
759 :
760 256 : ALLOCATE (bvector(3, MAXVAL(multi), max_shells))
761 16 : bvector = 0.0_dp; bweight = 0.0_dp
762 :
763 16 : WRITE (stdout, '(1x,a)') '| The b-vectors are chosen automatically |'
764 :
765 16 : b1sat = .FALSE.
766 16 : DO shell = 1, search_shells
767 16 : cur_shell = num_shells + 1
768 :
769 : ! get the b vectors for the new shell
770 16 : 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 16 : lpar = .FALSE.
774 16 : 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 16 : num_shells = num_shells + 1
795 16 : shell_list(num_shells) = shell
796 :
797 48 : ALLOCATE (amat(max_shells, num_shells))
798 16 : ALLOCATE (umat(max_shells, max_shells))
799 64 : ALLOCATE (vmat(num_shells, num_shells))
800 32 : ALLOCATE (smat(num_shells, max_shells))
801 48 : ALLOCATE (singv(num_shells))
802 16 : amat = 0.0_dp; umat = 0.0_dp; vmat = 0.0_dp; smat = 0.0_dp; singv = 0.0_dp
803 :
804 16 : amat = 0.0_dp
805 32 : DO loop_s = 1, num_shells
806 128 : DO loop_b = 1, multi(shell_list(loop_s))
807 96 : amat(1, loop_s) = amat(1, loop_s) + bvector(1, loop_b, loop_s)*bvector(1, loop_b, loop_s)
808 96 : amat(2, loop_s) = amat(2, loop_s) + bvector(2, loop_b, loop_s)*bvector(2, loop_b, loop_s)
809 96 : amat(3, loop_s) = amat(3, loop_s) + bvector(3, loop_b, loop_s)*bvector(3, loop_b, loop_s)
810 96 : amat(4, loop_s) = amat(4, loop_s) + bvector(1, loop_b, loop_s)*bvector(2, loop_b, loop_s)
811 96 : amat(5, loop_s) = amat(5, loop_s) + bvector(2, loop_b, loop_s)*bvector(3, loop_b, loop_s)
812 112 : 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 16 : info = 0
817 : CALL dgesvd('A', 'A', max_shells, num_shells, amat, max_shells, singv, umat, &
818 16 : max_shells, vmat, num_shells, work, lwork, info)
819 16 : 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 16 : ELSE IF (info > 0) THEN
823 0 : CPABORT('kmesh_shell_automatic: Singular Value Decomposition did not converge')
824 : END IF
825 :
826 32 : 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 16 : smat = 0.0_dp
840 32 : DO loop_s = 1, num_shells
841 32 : smat(loop_s, loop_s) = 1/singv(loop_s)
842 : END DO
843 :
844 240 : bweight(1:num_shells) = MATMUL(TRANSPOSE(vmat), MATMUL(smat, MATMUL(TRANSPOSE(umat), TARGET)))
845 16 : 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 64 : DO loop_i = 1, 3
855 160 : DO loop_j = loop_i, 3
856 96 : delta = 0.0_dp
857 192 : DO loop_s = 1, num_shells
858 768 : DO loop_b = 1, multi(shell_list(loop_s))
859 672 : 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 96 : IF (loop_i == loop_j) THEN
863 48 : IF (ABS(delta - 1.0_dp) > kmesh_tol) b1sat = .FALSE.
864 : END IF
865 144 : IF (loop_i /= loop_j) THEN
866 48 : IF (ABS(delta) > kmesh_tol) b1sat = .FALSE.
867 : END IF
868 : END DO
869 : END DO
870 :
871 16 : 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 16 : DEALLOCATE (amat, umat, vmat, smat, singv)
885 :
886 16 : IF (b1sat) EXIT
887 :
888 : END DO
889 :
890 16 : 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 16 : END SUBROUTINE kmesh_shell_automatic
900 :
901 : ! **************************************************************************************************
902 : !> \brief ...
903 : !> \param multi ...
904 : !> \param kpt ...
905 : !> \param shell_dist ...
906 : !> \param bvector ...
907 : ! **************************************************************************************************
908 16 : 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 400 : bvector = 0.0_dp
923 :
924 : num_bvec = 0
925 21312 : ok: DO loop = 1, (2*nsupcell + 1)**3
926 340736 : vkpp2 = MATMUL(lmn(:, loop), recip_lattice)
927 22904 : DO nkp2 = 1, num_kpts
928 91168 : vkpp = vkpp2 + kpt_cart(:, nkp2)
929 : dist = SQRT((kpt_cart(1, kpt) - vkpp(1))**2 &
930 22792 : + (kpt_cart(2, kpt) - vkpp(2))**2 + (kpt_cart(3, kpt) - vkpp(3))**2)
931 22792 : IF ((dist >= shell_dist*(1.0_dp - kmesh_tol)) .AND. dist <= shell_dist*(1.0_dp + kmesh_tol)) THEN
932 96 : num_bvec = num_bvec + 1
933 384 : bvector(:, num_bvec) = vkpp(:) - kpt_cart(:, kpt)
934 : END IF
935 : !if we have the right number of neighbours we can exit
936 22888 : IF (num_bvec == multi) CYCLE ok
937 : END DO
938 : END DO ok
939 :
940 16 : IF (num_bvec < multi) CPABORT('kmesh_get_bvector: Not enough bvectors found')
941 :
942 16 : 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 4512 : 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 18048 : ispos = SUM((a - b)**2) < eps8
956 18048 : isneg = SUM((a + b)**2) < eps8
957 4512 : END SUBROUTINE utility_compar
958 :
959 : ! **************************************************************************************************
960 :
961 16 : END MODULE wannier90
|
