1# Copyright (C) 2016 R. Warmbier Materials for Energy Research Group, 2# Wits University 3 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 9 10""" 11This file provides routines to create non-uniform k-point grids. We use 12locally uniform grids of different densities, practically grid refinement. 13 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. 16 17The user can further choose, whether the original symmetry of the system is 18enforced (careful!) or whether a reduced symmetry should be used. 19 20Optionally, the user can ask to add all k+q points, if not already included. 21 22Please cite https://doi.org/10.1016/j.cpc.2018.03.001 23 24Example (Graphene): 25 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 ) 32 33""" 34 35 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!') 43 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) 47 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)) 54 55 # Arguments needed for k-point descriptor construction 56 kwargs = {'nspins': nspins, 'atoms': atoms, 'symmetry': symmetry, 57 'comm': comm} 58 59 # Define coarse grid points 60 bzk_coarse_kc = bzkpts_kc 61 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) 66 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 76 77 # Create new descriptor with both sets of points 78 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) 85 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) 90 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") 103 104 return kd 105 106 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 113 114 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] 118 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) 126 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] 137 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 141 142 # return kd_new.copy() 143 return kd_new 144 145 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. """ 149 150 # Define Monkhorst-Pack grid as starting point 151 kd_coarse = create_kpoint_descriptor(bzk_coarse_kc, **kwargs) 152 153 # Determine how many (cubic) nearest neighbour shells to replace 154 nshells = size.shape[0] 155 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) 159 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) 171 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) 178 179 # Remove redundant entries, which come from symmetry 180 center_i, index_map = np.unique(center_i, return_index=True) 181 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) 187 188 # Assign the kpoint vectors corresponding to the indexes 189 centers_kc = bzk_coarse_kc[center_i] 190 191 # Get scaled Monkhorst-Pack for refinement 192 mh_kc = get_reduced_monkhorst(size[shell], kd_coarse.N_c) 193 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) 198 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) 204 205 weight_k = np.concatenate(weight_k) 206 bzk_fine_kc = np.concatenate(bzk_fine_kc) 207 208 assert np.abs(len(centers_i) - sum(weight_k)) < 1e-6 209 assert len(weight_k) == bzk_fine_kc.shape[0] 210 211 return centers_i, bzk_fine_kc, weight_k 212 213 214def find_equivalent_kpoints(point, kd): 215 """Find coordinates and BZ index for all refinement centers.""" 216 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!') 221 222 equiv_k = list(set(kd.bz2bz_ks[k])) 223 # equiv_kc = kd.bzk_kc[equiv_k] 224 225 return np.array(equiv_k) # , equiv_kc 226 227 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.""" 231 232 mh_kc = monkhorst_pack(size) 233 return mh_kc / N_c 234 235 236def construct_neighbours_by_shells(nshells, N_c): 237 """Construct a list of neighbours (translations from center for each 238 shell around the center.""" 239 240 neighbours = [] 241 242 # The innermost shell is always the center itself 243 neighbours.append(np.array([0, 0, 0], dtype=float, ndmin=2)) 244 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))) 254 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) 265 266 this_list = np.array(this_list, dtype=float, ndmin=2) / N_c 267 neighbours.append(this_list) 268 269 return np.array(neighbours) 270 271 272def prune_symmetries_kpoints(kd, symmetry): 273 """Remove the symmetries which are violated by a kpoint set.""" 274 275 new_symmetry = copy.copy(symmetry) 276 277 nsyms = len(symmetry.op_scc) 278 279 where = np.where(~np.any(kd.bz2bz_ks[:, 0:nsyms] == -1, 0)) 280 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] 284 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() 288 289 return new_symmetry 290 291 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 297 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))) 317 318 return buf 319 320 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 326 327 kd_new = create_new_descriptor_with_zero_points(kd, add_points_kc, kwargs) 328 assert -1 not in kd_new.bz2bz_ks 329 330 return kd_new 331 332 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.""" 336 337 # Add missing points to retrieve full symmetry. Might be redundant. 338 _kd = add_missing_points(kd, kwargs) 339 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) 353 354 if add_points_kc.shape[0] == 0: 355 return _kd 356 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 363 364 kd_new = create_new_descriptor_with_zero_points(kd, add_points_kc, kwargs) 365 return kd_new 366 367 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] 373 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] 377 378 # Construct the new KPointDescriptor 379 kd_new = create_kpoint_descriptor(bzk_kc, **kwargs) 380 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) 389 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 394 395 return kd_new 396 397 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 403 404 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 413 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 420 421 def set_unrefined_nbzkpts(self, mhnbzkpts): 422 self.mhnbzkpts = mhnbzkpts 423 424 def get_unrefined_nbzkpts(self): 425 return self.mhnbzkpts 426 427 def set_weight_k(self, weight_k): 428 self.weight_k = weight_k 429 430 def get_weight_k(self): 431 return self.weight_k 432 433 def set_label_k(self, label_k): 434 self.label_k = label_k 435 436 def get_label_k(self): 437 return self.label_k 438 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 446