Coverage for gpaw/kpt_refine.py: 72%
236 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1# Copyright (C) 2016 R. Warmbier Materials for Energy Research Group,
2# Wits University
4import copy
5import numpy as np
6from ase.dft.kpoints import monkhorst_pack, get_monkhorst_pack_size_and_offset
7from gpaw import KPointError
8from gpaw.kpt_descriptor import KPointDescriptor
10"""
11This file provides routines to create non-uniform k-point grids. We use
12locally uniform grids of different densities, practically grid refinement.
14One starts from a standard Monkhorst-Pack grid and chooses a number of grid
15points, which shall be replaced with a finer grid of a given size.
17The user can further choose, whether the original symmetry of the system is
18enforced (careful!) or whether a reduced symmetry should be used.
20Optionally, the user can ask to add all k+q points, if not already included.
22Please cite https://doi.org/10.1016/j.cpc.2018.03.001
24Example (Graphene):
26calc = GPAW(mode=PW(ecut=450),
27 kpts={"size":[15,15,1], "gamma":True},
28 kpt_refine={"center":[1./3,1./3,0.], "size":[3,3,1],
29 "reduce_symmetry":False,
30 "q":[1./15,1./15,0.]},
31 )
33"""
36def create_kpoint_descriptor_with_refinement(refine, bzkpts_kc, nspins, atoms,
37 symmetry, comm, timer):
38 """Main routine to build refined k-point grids."""
39 if 'center' not in refine:
40 raise RuntimeError('Center for refinement not given!')
41 if 'size' not in refine:
42 raise RuntimeError('Grid size for refinement not given!')
44 center_ic = np.array(refine.get('center'), dtype=float, ndmin=2)
45 size = np.array(refine.get('size'), ndmin=2)
46 reduce_symmetry = refine.get('reduce_symmetry', True)
48 # Check that all sizes are odd. That's not so much an issue really.
49 # But even Monkhorst-Pack grids have points on the boundary,
50 # which just would require more special casing, which I want to avoid.
51 if (np.array(size) % 2 == 0).any():
52 raise RuntimeError('Grid size for refinement must be odd! Is: {}'.
53 format(size))
55 # Arguments needed for k-point descriptor construction
56 kwargs = {'nspins': nspins, 'atoms': atoms, 'symmetry': symmetry,
57 'comm': comm}
59 # Define coarse grid points
60 bzk_coarse_kc = bzkpts_kc
62 # Define fine grid points
63 centers_i, bzk_fine_kc, weight_fine_k = get_fine_bzkpts(center_ic, size,
64 bzk_coarse_kc,
65 kwargs)
67 if reduce_symmetry:
68 # Define new symmetry object ignoring symmetries violated by the
69 # refined kpoints
70 kd_fine = create_kpoint_descriptor(bzk_fine_kc, **kwargs)
71 symm = prune_symmetries_kpoints(kd_fine, symmetry)
72 del kd_fine
73 else:
74 symm = copy.copy(symmetry)
75 kwargs['symmetry'] = symm
77 # Create new descriptor with both sets of points
79 with timer('Create mixed descriptor'):
80 kd = create_mixed_kpoint_descriptor(bzk_coarse_kc,
81 bzk_fine_kc,
82 centers_i,
83 weight_fine_k,
84 kwargs)
86 # Add missing k-points to fulfill group properties with
87 # zero-weighted points
88 with timer('Add_missing_points'):
89 kd = add_missing_points(kd, kwargs)
91 # Add additional +q k-points, if necessary
92 if 'q' in refine:
93 timer.start("+q")
94 N_coarse_c = get_monkhorst_pack_size_and_offset(bzk_coarse_kc)[0]
95 bla = N_coarse_c * refine['q']
96 if not max(abs(bla - np.rint(bla))) < 1e-8:
97 kd.refine_info.almostoptical = True
98 kd = add_plusq_points(kd, refine['q'], kwargs)
99 symm = kd.symmetry
100 kwargs['symmetry'] = symm
101 kd = add_missing_points(kd, kwargs)
102 timer.stop("+q")
104 return kd
107def create_kpoint_descriptor(bzkpts_kc, nspins, atoms, symmetry, comm):
108 kd = KPointDescriptor(bzkpts_kc, nspins)
109 # self.timer.start('Set symmetry')
110 kd.set_symmetry(atoms, symmetry, comm=comm)
111 # self.timer.stop('Set symmetry')
112 return kd
115def create_mixed_kpoint_descriptor(bzk_coarse_kc, bzk_fine_kc, centers_i,
116 weight_fine_k, kwargs):
117 assert len(weight_fine_k) == bzk_fine_kc.shape[0]
119 # Get old Monkhorst-Pack grid points and delete the centers from them
120 nbzkpts_coarse = bzk_coarse_kc.shape[0]
121 bzk_new_kc = np.delete(bzk_coarse_kc, centers_i, axis=0)
122 nbzkpts_coarse_new = bzk_new_kc.shape[0]
123 # Add refined points
124 nbzkpts_fine = bzk_fine_kc.shape[0]
125 bzk_new_kc = np.append(bzk_new_kc, bzk_fine_kc, axis=0)
127 # Construct the new KPointDescriptor
128 kd_new = create_kpoint_descriptor(bzk_new_kc, **kwargs)
129 refine_info = KRefinement()
130 refine_info.set_unrefined_nbzkpts(nbzkpts_coarse)
131 label_k = np.array(nbzkpts_coarse_new * ['mh'] + nbzkpts_fine * ['refine'])
132 refine_info.set_label_k(label_k)
133 weight_k = np.append(np.ones(nbzkpts_coarse_new), weight_fine_k)
134 refine_info.set_weight_k(weight_k)
135 kd_new.refine_info = refine_info
136 assert len(weight_k) == bzk_new_kc.shape[0]
138 # This new Descriptor is good, except the weights are completely wrong now.
139 kd_new.weight_k = np.bincount(kd_new.bz2ibz_k, weight_k) / nbzkpts_coarse
140 assert np.abs(np.sum(kd_new.weight_k) - 1) < 1e-6
142 # return kd_new.copy()
143 return kd_new
146def get_fine_bzkpts(center_ic, size, bzk_coarse_kc, kwargs):
147 """Return list of reducible refined kpoints and the indexes for the points
148 on the coarse kpoint grid. """
150 # Define Monkhorst-Pack grid as starting point
151 kd_coarse = create_kpoint_descriptor(bzk_coarse_kc, **kwargs)
153 # Determine how many (cubic) nearest neighbour shells to replace
154 nshells = size.shape[0]
156 # Get coordinates of neighbour k-points. First index is the 'shell' aka
157 # order of neighbour. Coordinates always centred around zero.
158 neighbours_kc = construct_neighbours_by_shells(nshells, kd_coarse.N_c)
160 # loop over all the different shells
161 centers_i = np.empty((0), dtype=int)
162 bzk_fine_kc = []
163 weight_k = []
164 for shell in range(nshells):
165 # Determine all k-points in shell, which will be refinement centers
166 shell_kc = []
167 for i, center_c in enumerate(center_ic):
168 shell_c = neighbours_kc[shell] + center_c
169 shell_kc.append(shell_c)
170 shell_kc = np.concatenate(shell_kc)
172 # Find all equivalent centers using full symmetry
173 center_i = []
174 for k, shell_c in enumerate(shell_kc):
175 equiv_i = find_equivalent_kpoints(shell_c, kd_coarse)
176 center_i.append(equiv_i)
177 center_i = np.concatenate(center_i)
179 # Remove redundant entries, which come from symmetry
180 center_i, index_map = np.unique(center_i, return_index=True)
182 # Append to list for all shells
183 # Remove also points which occur also in a previous shell. This can
184 # happen due to symmetry or overlap of centers' shells
185 center_i = np.setdiff1d(center_i, centers_i)
186 centers_i = np.append(centers_i, center_i)
188 # Assign the kpoint vectors corresponding to the indexes
189 centers_kc = bzk_coarse_kc[center_i]
191 # Get scaled Monkhorst-Pack for refinement
192 mh_kc = get_reduced_monkhorst(size[shell], kd_coarse.N_c)
194 # Gather the refined points of the grid for all centers
195 for k, centers_c in enumerate(centers_kc):
196 this_kc = mh_kc + centers_c
197 bzk_fine_kc.append(this_kc)
199 # Determine weight of refined points
200 weight = 1. / np.prod(size[shell])
201 nkpts = len(center_i) * np.prod(size[shell])
202 weight *= np.ones(nkpts)
203 weight_k.append(weight)
205 weight_k = np.concatenate(weight_k)
206 bzk_fine_kc = np.concatenate(bzk_fine_kc)
208 assert np.abs(len(centers_i) - sum(weight_k)) < 1e-6
209 assert len(weight_k) == bzk_fine_kc.shape[0]
211 return centers_i, bzk_fine_kc, weight_k
214def find_equivalent_kpoints(point, kd):
215 """Find coordinates and BZ index for all refinement centers."""
217 # check whether point in bzkpts_kc
218 min_dist, k = minimal_point_distance(point, kd.bzk_kc)
219 if min_dist > 1e-8:
220 raise RuntimeError('Could not find k-point in kd list!')
222 equiv_k = list(set(kd.bz2bz_ks[k]))
223 # equiv_kc = kd.bzk_kc[equiv_k]
225 return np.array(equiv_k) # , equiv_kc
228def get_reduced_monkhorst(size, N_c):
229 """Find monkhorst_pack according to size of refined grid and shrink it into
230 the volume of one original kpoint."""
232 mh_kc = monkhorst_pack(size)
233 return mh_kc / N_c
236def construct_neighbours_by_shells(nshells, N_c):
237 """Construct a list of neighbours (translations from center for each
238 shell around the center."""
240 neighbours = []
242 # The innermost shell is always the center itself
243 neighbours.append(np.array([0, 0, 0], dtype=float, ndmin=2))
245 # Loop through the other shells
246 for shell in range(1, nshells):
247 # Construct the displacements/elements for each dimension.
248 elements = []
249 for i in range(3):
250 if N_c[i] == 1:
251 elements.append([0, ])
252 else:
253 elements.append(list(range(-shell, shell + 1)))
255 # Construct the vectors
256 # For each valid point, at least one component must have the value
257 # of the shell index (+ or -).
258 this_list = []
259 for cx in elements[0]:
260 for cy in elements[1]:
261 for cz in elements[2]:
262 candidate = [cx, cy, cz]
263 if (np.abs(np.array(candidate)) == shell).any():
264 this_list.append(candidate)
266 this_list = np.array(this_list, dtype=float, ndmin=2) / N_c
267 neighbours.append(this_list)
269 return np.array(neighbours)
272def prune_symmetries_kpoints(kd, symmetry):
273 """Remove the symmetries which are violated by a kpoint set."""
275 new_symmetry = copy.copy(symmetry)
277 nsyms = len(symmetry.op_scc)
279 where = np.where(~np.any(kd.bz2bz_ks[:, 0:nsyms] == -1, 0))
281 new_symmetry.op_scc = new_symmetry.op_scc[where]
282 new_symmetry.ft_sc = new_symmetry.ft_sc[where]
283 new_symmetry.a_sa = new_symmetry.a_sa[where]
285 inv_cc = -np.eye(3, dtype=int)
286 new_symmetry.has_inversion = (new_symmetry.op_scc ==
287 inv_cc).all(2).all(1).any()
289 return new_symmetry
292def find_missing_points(kd):
293 """Find points in k-point descriptor, which miss in the set to fill the
294 group."""
295 if -1 not in kd.bz2bz_ks:
296 return None
298 # Find unaccounted points
299 buf = np.empty((0, 3))
300 for k in range(kd.bz2bz_ks.shape[0]):
301 for s in range(kd.bz2bz_ks.shape[1]):
302 if kd.bz2bz_ks[k, s] == -1:
303 sign = 1.0
304 i = s
305 if kd.symmetry.time_reversal:
306 if s / kd.bz2bz_ks.shape[1] >= 0.5:
307 i = s - int(kd.bz2bz_ks.shape[1] / 2)
308 sign = -1.
309 k_c = sign * np.dot(kd.symmetry.op_scc[i], kd.bzk_kc[k])
310 k_c -= np.rint(k_c)
311 # Check that k_c hasn't been added yet
312 if buf.shape[0] > 0:
313 min_dist, index = minimal_point_distance(k_c, buf)
314 if min_dist < 1e-6:
315 continue
316 buf = np.concatenate((buf, np.array(k_c, ndmin=2)))
318 return buf
321def add_missing_points(kd, kwargs):
322 """Add points to k-point descriptor, which miss to fill the group."""
323 add_points_kc = find_missing_points(kd)
324 if add_points_kc is None:
325 return kd
327 kd_new = create_new_descriptor_with_zero_points(kd, add_points_kc, kwargs)
328 assert -1 not in kd_new.bz2bz_ks
330 return kd_new
333def add_plusq_points(kd, q_c, kwargs):
334 """Add +q points to k-point descriptor, if missing. Also, reduce the
335 symmetry of the system as necessary."""
337 # Add missing points to retrieve full symmetry. Might be redundant.
338 _kd = add_missing_points(kd, kwargs)
340 # Find missing q
341 add_points_kc = []
342 for k in range(_kd.nbzkpts):
343 # if q_c is small, use q_c = 0.0 for mh points, else they don't need
344 # extra points anyway
345 if _kd.refine_info.label_k[k] == 'mh':
346 continue
347 try:
348 _kd.find_k_plus_q(q_c, [k])
349 except KPointError:
350 k_c = _kd.bzk_kc[k] + q_c
351 add_points_kc.append(np.array(k_c, ndmin=2))
352 add_points_kc = np.concatenate(add_points_kc)
354 if add_points_kc.shape[0] == 0:
355 return _kd
357 # Find reduced +q symmetry
358 bzk_kc = np.append(_kd.bzk_kc, add_points_kc, axis=0)
359 _kd_new = create_kpoint_descriptor(bzk_kc, **kwargs)
360 symm_new = prune_symmetries_kpoints(_kd_new, kwargs.get('symmetry'))
361 kwargs['symmetry'] = symm_new
362 del _kd_new, _kd
364 kd_new = create_new_descriptor_with_zero_points(kd, add_points_kc, kwargs)
365 return kd_new
368def create_new_descriptor_with_zero_points(kd, add_points_kc, kwargs):
369 """Create a new k-point descriptor including some additional zero-weighted
370 points."""
371 bzk_kc = np.append(kd.bzk_kc, add_points_kc, axis=0)
372 nbzkpts_add = add_points_kc.shape[0]
374 # Make sure all points are unique
375 assert bzk_kc.shape[0] == np.vstack([tuple(row)
376 for row in bzk_kc]).shape[0]
378 # Construct the new KPointDescriptor
379 kd_new = create_kpoint_descriptor(bzk_kc, **kwargs)
381 # Update refine_info
382 kd_new.refine_info = kd.refine_info.copy()
383 label_k = np.append(kd.refine_info.label_k,
384 np.array(nbzkpts_add * ['zero']))
385 kd_new.refine_info.set_label_k(label_k)
386 # Avoid exact zero, as that would screw up the occupations calculation
387 weight_k = np.append(kd.refine_info.weight_k, 1e-10 * np.ones(nbzkpts_add))
388 kd_new.refine_info.set_weight_k(weight_k)
390 # Correct ibz weights
391 kd_new.weight_k = np.bincount(kd_new.bz2ibz_k, weight_k)
392 kd_new.weight_k *= 1.0 / kd_new.refine_info.get_unrefined_nbzkpts()
393 assert np.abs(np.sum(kd_new.weight_k) - 1) < 1e-6
395 return kd_new
398def minimal_point_distance(point_c, bzk_kc):
399 d_kc = point_c - bzk_kc
400 d_k = abs(d_kc - d_kc.round()).sum(1)
401 k = d_k.argmin()
402 return d_k[k], k
405class KRefinement:
406 """Additional information for refined kpoint grids.
407 """
408 def __init__(self):
409 self.mhnbzkpts = None
410 self.label_k = None
411 self.weight_k = None
412 self.almostoptical = None
414 def __str__(self):
415 s = "Using k-point grid refinement"
416 if self.almostoptical:
417 s += " with almostoptical approximation"
418 s += "\n"
419 return s
421 def set_unrefined_nbzkpts(self, mhnbzkpts):
422 self.mhnbzkpts = mhnbzkpts
424 def get_unrefined_nbzkpts(self):
425 return self.mhnbzkpts
427 def set_weight_k(self, weight_k):
428 self.weight_k = weight_k
430 def get_weight_k(self):
431 return self.weight_k
433 def set_label_k(self, label_k):
434 self.label_k = label_k
436 def get_label_k(self):
437 return self.label_k
439 def copy(self):
440 refine_info = KRefinement()
441 refine_info.mhnbzkpts = self.mhnbzkpts
442 refine_info.label_k = self.label_k
443 refine_info.weight_k = self.weight_k
444 refine_info.almostoptical = self.almostoptical
445 return refine_info