Coverage for gpaw/bztools.py: 90%
183 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1from itertools import product
3import numpy as np
4from ase.dft.kpoints import monkhorst_pack
5from scipy.spatial import ConvexHull, Delaunay, Voronoi
6try:
7 from scipy.spatial import QhullError
8except ImportError: # scipy < 1.8
9 from scipy.spatial.qhull import QhullError
11import gpaw.mpi as mpi
12from gpaw import GPAW, restart
13from gpaw.kpt_descriptor import kpts2sizeandoffsets, to1bz
14from gpaw.mpi import world
15from gpaw.symmetry import Symmetry, aglomerate_points
18def get_lattice_symmetry(cell_cv, tolerance=1e-7):
19 """Return symmetry object of lattice group.
21 Parameters
22 ----------
23 cell_cv : ndarray
24 Unit cell.
26 Returns
27 -------
28 gpaw.symmetry object
30 """
31 # NB: Symmetry.find_lattice_symmetry() uses self.pbc_c, which defaults
32 # to pbc along all three dimensions. Hence, it seems that the lattice
33 # symmetry transformations produced by this method could be faulty if
34 # there are non-periodic dimensions in a system. XXX
35 latsym = Symmetry([0], cell_cv, tolerance=tolerance)
36 latsym.find_lattice_symmetry()
37 return latsym
40def find_high_symmetry_monkhorst_pack(calc, density):
41 """Make high symmetry Monkhorst Pack k-point grid.
43 Searches for and returns a Monkhorst Pack grid which
44 contains the corners of the irreducible BZ so that when the
45 number of kpoints are reduced the full irreducible brillouion
46 zone is spanned.
48 Parameters
49 ----------
50 calc : str
51 The path to a calculator object.
52 density : float
53 The required minimum density of the Monkhorst Pack grid.
55 Returns
56 -------
57 ndarray
58 Array of shape (nk, 3) containing the kpoints.
60 """
62 atoms, calc = restart(calc, txt=None)
63 pbc = atoms.pbc
64 minsize, offset = kpts2sizeandoffsets(density=density, even=True,
65 gamma=True, atoms=atoms)
67 # NB: get_bz() wants a pbc_c, but never gets it. This means that the
68 # pbc always will fall back to True along all dimensions. XXX
69 # NB: Why return latibzk_kc, if we never use it? XXX
70 bzk_kc, ibzk_kc, latibzk_kc = get_bz(calc)
72 maxsize = minsize + 10
73 minsize[~pbc] = 1
74 maxsize[~pbc] = 2
76 if mpi.rank == 0:
77 print('Brute force search for symmetry ' +
78 'complying MP-grid... please wait.')
80 for n1 in range(minsize[0], maxsize[0]):
81 for n2 in range(minsize[1], maxsize[1]):
82 for n3 in range(minsize[2], maxsize[2]):
83 size = [n1, n2, n3]
84 size, offset = kpts2sizeandoffsets(size=size, gamma=True,
85 atoms=atoms)
87 ints = ((ibzk_kc + 0.5 - offset) * size - 0.5)[:, pbc]
89 if (np.abs(ints - np.round(ints)) < 1e-5).all():
90 kpts_kc = monkhorst_pack(size) + offset
91 kpts_kc = to1bz(kpts_kc, calc.wfs.gd.cell_cv)
93 for ibzk_c in ibzk_kc:
94 diff_kc = np.abs(kpts_kc - ibzk_c)[:, pbc].round(6)
95 # NB: The second np.mod() statement seems redundant XXX
96 if not (np.mod(np.mod(diff_kc, 1), 1) <
97 1e-5).all(axis=1).any():
98 raise AssertionError('Did not find ' + str(ibzk_c))
99 if mpi.rank == 0:
100 print('Done. Monkhorst-Pack grid:', size, offset)
101 return kpts_kc
103 if mpi.rank == 0:
104 print('Did not find matching kpoints for the IBZ')
105 print(ibzk_kc.round(5))
107 raise RuntimeError
110def unfold_points(points, U_scc, tol=1e-8, mod=None):
111 """Unfold k-points using a given set of symmetry operators.
113 Parameters
114 ----------
115 points: ndarray
116 U_scc: ndarray
117 tol: float
118 Tolerance indicating when kpoints are considered to be
119 identical.
120 mod: integer 1 or None
121 Consider kpoints spaced by a full reciprocal lattice vector
122 to be identical.
124 Returns
125 -------
126 ndarray
127 Array of shape (nk, 3) containing the unfolded kpoints.
128 """
130 points = np.concatenate(np.dot(points, U_scc.transpose(0, 2, 1)))
131 return unique_rows(points, tol=tol, mod=mod)
134def unique_rows(ain, tol=1e-10, mod=None, aglomerate=True):
135 """Return unique rows of a 2D ndarray.
137 Parameters
138 ----------
139 ain : 2D ndarray
140 tol : float
141 Tolerance indicating when kpoints are considered to be
142 identical.
143 mod : integer 1 or None
144 Consider kpoints spaced by a full reciprocal lattice vector
145 to be identical.
146 aglomerate : bool
147 Aglomerate clusters of points before comparing.
149 Returns
150 -------
151 2D ndarray
152 Array containing only unique rows.
153 """
154 # Move to positive octant
155 a = ain - ain.min(0)
157 # First take modulus
158 if mod is not None:
159 a = np.mod(np.mod(a, mod), mod)
161 # Round and take modulus again
162 if aglomerate:
163 aglomerate_points(a, tol)
164 a = a.round(-np.log10(tol).astype(int))
165 if mod is not None:
166 a = np.mod(a, mod)
168 # Now perform ordering
169 order = np.lexsort(a.T)
170 a = a[order]
172 # Find unique rows
173 diff = np.diff(a, axis=0)
174 ui = np.ones(len(a), 'bool')
175 ui[1:] = (diff != 0).any(1)
177 return ain[order][ui]
180def get_smallest_Gvecs(cell_cv, n=5):
181 """Find smallest reciprocal lattice vectors.
183 Parameters
184 ----------
185 cell_cv : ndarray
186 Unit cell.
187 n : int
188 Sampling along each crystal axis.
190 Returns
191 -------
192 G_xv : ndarray
193 Reciprocal lattice vectors in cartesian coordinates.
194 N_xc : ndarray
195 Reciprocal lattice vectors in crystal coordinates.
197 """
198 B_cv = 2.0 * np.pi * np.linalg.inv(cell_cv).T
199 N_xc = np.indices((n, n, n)).reshape((3, n**3)).T - n // 2
200 G_xv = N_xc @ B_cv
202 return G_xv, N_xc
205def get_symmetry_operations(U_scc, time_reversal):
206 """Return point symmetry operations."""
208 if U_scc is None:
209 U_scc = np.array([np.eye(3)])
211 inv_cc = -np.eye(3, dtype=int)
212 has_inversion = (U_scc == inv_cc).all(2).all(1).any()
214 if has_inversion:
215 time_reversal = False
217 if time_reversal:
218 Utmp_scc = np.concatenate([U_scc, -U_scc])
219 else:
220 Utmp_scc = U_scc
222 return Utmp_scc
225def get_ibz_vertices(cell_cv, U_scc=None, time_reversal=None,
226 origin_c=None):
227 """Determine irreducible BZ.
229 Parameters
230 ----------
231 cell_cv : ndarray
232 Unit cell
233 U_scc : ndarray
234 Crystal symmetry operations.
235 time_reversal : bool
236 Use time reversal symmetry?
238 Returns
239 -------
240 ibzk_kc : ndarray
241 Vertices of the irreducible BZ.
242 """
243 # Choose an origin
244 if origin_c is None:
245 origin_c = np.array([0.12, 0.22, 0.21], float)
246 else:
247 assert (np.abs(origin_c) < 0.5).all()
249 if U_scc is None:
250 U_scc = np.array([np.eye(3)])
252 if time_reversal is None:
253 time_reversal = False
255 Utmp_scc = get_symmetry_operations(U_scc, time_reversal)
257 icell_cv = np.linalg.inv(cell_cv).T
258 B_cv = icell_cv * 2 * np.pi
259 A_cv = np.linalg.inv(B_cv).T
261 # Map a random point around
262 point_sc = np.dot(origin_c, Utmp_scc.transpose((0, 2, 1)))
263 assert len(point_sc) == len(unique_rows(point_sc))
264 point_sv = np.dot(point_sc, B_cv)
266 # Translate the points
267 n = 5
268 G_xv, N_xc = get_smallest_Gvecs(cell_cv, n=n)
269 G_xv = np.delete(G_xv, n**3 // 2, axis=0)
271 # Mirror points in plane
272 N_xv = G_xv / (((G_xv**2).sum(1))**0.5)[:, np.newaxis]
274 tp_sxv = (point_sv[:, np.newaxis] - G_xv[np.newaxis] / 2.)
275 delta_sxv = ((tp_sxv * N_xv[np.newaxis]).sum(2)[..., np.newaxis] *
276 N_xv[np.newaxis])
277 points_xv = (point_sv[:, np.newaxis] - 2 * delta_sxv).reshape((-1, 3))
278 points_xv = np.concatenate([point_sv, points_xv])
279 try:
280 voronoi = Voronoi(points_xv)
281 except QhullError:
282 return get_ibz_vertices(cell_cv, U_scc=U_scc,
283 time_reversal=time_reversal,
284 origin_c=origin_c + [0.01, -0.02, -0.01])
286 ibzregions = voronoi.point_region[0:len(point_sv)]
288 ibzregion = ibzregions[0]
289 ibzk_kv = voronoi.vertices[voronoi.regions[ibzregion]]
290 ibzk_kc = np.dot(ibzk_kv, A_cv.T)
292 return ibzk_kc
295def get_bz(calc, pbc_c=np.ones(3, bool)):
296 """Return the BZ and IBZ vertices.
298 Parameters
299 ----------
300 calc : str, GPAW calc instance
302 Returns
303 -------
304 bzk_kc : ndarray
305 Vertices of BZ in crystal coordinates
306 ibzk_kc : ndarray
307 Vertices of IBZ in crystal coordinates
309 """
311 if isinstance(calc, str):
312 calc = GPAW(calc, txt=None)
313 cell_cv = calc.wfs.gd.cell_cv
315 # Crystal symmetries
316 symmetry = calc.wfs.kd.symmetry
317 cU_scc = get_symmetry_operations(symmetry.op_scc,
318 symmetry.time_reversal)
320 return get_reduced_bz(cell_cv, cU_scc, False, pbc_c=pbc_c)
323def get_reduced_bz(cell_cv, cU_scc, time_reversal,
324 pbc_c=np.ones(3, bool), tolerance=1e-7):
326 """Reduce the BZ using the crystal symmetries to obtain the IBZ.
328 Parameters
329 ----------
330 cell_cv : ndarray
331 Unit cell.
332 cU_scc : ndarray
333 Crystal symmetry operations.
334 time_reversal : bool
335 Switch for time reversal.
336 pbc: bool or [bool, bool, bool]
337 Periodic bcs
338 """
340 if time_reversal:
341 # NB: The method never seems to be called with time_reversal=True,
342 # and hopefully get_bz() will generate the right symmetry operations
343 # always. So, can we remove this input? XXX
344 cU_scc = get_symmetry_operations(cU_scc, time_reversal)
346 # Lattice symmetries
347 latsym = get_lattice_symmetry(cell_cv, tolerance=tolerance)
348 lU_scc = get_symmetry_operations(latsym.op_scc,
349 latsym.time_reversal)
351 # Find Lattice IBZ
352 ibzk_kc = get_ibz_vertices(cell_cv,
353 U_scc=latsym.op_scc,
354 time_reversal=latsym.time_reversal)
355 latibzk_kc = ibzk_kc.copy()
357 # Expand lattice IBZ to crystal IBZ
358 ibzk_kc = expand_ibz(lU_scc, cU_scc, ibzk_kc, pbc_c=pbc_c)
360 # Fold out to full BZ
361 bzk_kc = unique_rows(np.concatenate(np.dot(ibzk_kc,
362 cU_scc.transpose(0, 2, 1))))
364 return bzk_kc, ibzk_kc, latibzk_kc
367def expand_ibz(lU_scc, cU_scc, ibzk_kc, pbc_c=np.ones(3, bool)):
368 """Expand IBZ from lattice group to crystal group.
370 Parameters
371 ----------
372 lU_scc : ndarray
373 Lattice symmetry operators.
374 cU_scc : ndarray
375 Crystal symmetry operators.
376 ibzk_kc : ndarray
377 Vertices of lattice IBZ.
379 Returns
380 -------
381 ibzk_kc : ndarray
382 Vertices of crystal IBZ.
384 """
386 # Find right cosets. The lattice group is partioned into right cosets of
387 # the crystal group. This can in practice be done by testing whether
388 # U1 U2^{-1} is in the crystal group as done below.
389 cosets = []
390 Utmp_scc = lU_scc.copy()
391 while len(Utmp_scc):
392 U1_cc = Utmp_scc[0].copy()
393 Utmp_scc = np.delete(Utmp_scc, 0, axis=0)
394 j = 0
395 new_coset = [U1_cc]
396 while j < len(Utmp_scc):
397 U2_cc = Utmp_scc[j]
398 U3_cc = np.dot(U1_cc, np.linalg.inv(U2_cc))
399 if (U3_cc == cU_scc).all(2).all(1).any():
400 new_coset.append(U2_cc)
401 Utmp_scc = np.delete(Utmp_scc, j, axis=0)
402 j -= 1
403 j += 1
404 cosets.append(new_coset)
406 volume = np.inf
407 nibzk_kc = ibzk_kc
408 U0_cc = cosets[0][0] # Origin
410 if np.any(~pbc_c):
411 nonpbcind = np.argwhere(~pbc_c)
413 # Once the coests are known the irreducible zone is given by picking one
414 # operation from each coset. To make sure that the IBZ produced is simply
415 # connected we compute the volume of the convex hull of the produced IBZ
416 # and pick (one of) the ones that have the smallest volume. This is done by
417 # brute force and can sometimes take a while, however, in most cases this
418 # is not a problem.
419 combs = list(product(*cosets[1:]))[world.rank::world.size]
420 for U_scc in combs:
421 if not len(U_scc):
422 continue
423 U_scc = np.concatenate([np.array(U_scc), [U0_cc]])
424 tmpk_kc = unfold_points(ibzk_kc, U_scc)
425 volumenew = convex_hull_volume(tmpk_kc)
427 if np.any(~pbc_c):
428 # Compute the area instead
429 volumenew /= (tmpk_kc[:, nonpbcind].max() -
430 tmpk_kc[:, nonpbcind].min())
432 if volumenew < volume:
433 nibzk_kc = tmpk_kc
434 volume = volumenew
436 ibzk_kc = unique_rows(nibzk_kc)
437 volume = np.array((volume,))
439 volumes = np.zeros(world.size, float)
440 world.all_gather(volume, volumes)
442 minrank = np.argmin(volumes)
443 minshape = np.array(ibzk_kc.shape)
444 world.broadcast(minshape, minrank)
446 if world.rank != minrank:
447 ibzk_kc = np.zeros(minshape, float)
448 world.broadcast(ibzk_kc, minrank)
450 return ibzk_kc
453def tetrahedron_volume(a, b, c, d):
454 """Calculate volume of tetrahedron.
456 Parameters
457 ----------
458 a, b, c, d : ndarray
459 Vertices of tetrahedron.
461 Returns
462 -------
463 float
464 Volume of tetrahedron.
466 """
467 return np.abs(np.einsum('ij,ij->i', a - d,
468 np.cross(b - d, c - d))) / 6
471def convex_hull_volume(pts):
472 """Calculate volume of the convex hull of a collection of points.
474 Parameters
475 ----------
476 pts : list, ndarray
477 A list of 3d points.
479 Returns
480 -------
481 float
482 Volume of convex hull.
484 """
485 hull = ConvexHull(pts)
486 dt = Delaunay(pts[hull.vertices])
487 tets = dt.points[dt.simplices]
488 vol = np.sum(tetrahedron_volume(tets[:, 0], tets[:, 1],
489 tets[:, 2], tets[:, 3]))
490 return vol