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