Coverage for gpaw/response/kpoints.py: 99%
93 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1from functools import cached_property
3import numpy as np
4from scipy.spatial import Delaunay, cKDTree
6from gpaw.bztools import get_reduced_bz, unique_rows
9class KPointFinder:
10 def __init__(self, bzk_kc):
11 self.kdtree = cKDTree(self._round(bzk_kc))
13 @staticmethod
14 def _round(bzk_kc):
15 return np.mod(np.mod(bzk_kc, 1).round(6), 1)
17 def find(self, kpt_c):
18 distance, k = self.kdtree.query(self._round(kpt_c))
19 if distance > 1.e-6:
20 raise ValueError('Requested k-point is not on the grid. '
21 'Please check that your q-points of interest '
22 'are commensurate with the k-point grid.')
24 return k
27class ResponseKPointGrid:
28 def __init__(self, kd):
29 self.kd = kd
31 @cached_property
32 def kptfinder(self):
33 return KPointFinder(self.kd.bzk_kc)
36class KPointDomain:
37 def __init__(self, k_kc, icell_cv):
38 self.k_kc = k_kc
39 self.icell_cv = icell_cv
41 def __len__(self):
42 return len(self.k_kc)
44 @cached_property
45 def k_kv(self):
46 return self.k_kc @ (2 * np.pi * self.icell_cv)
49class KPointDomainGenerator:
50 def __init__(self, symmetries, kpoints):
51 self.symmetries = symmetries
53 self.kd = kpoints.kd
54 self.kptfinder = kpoints.kptfinder
56 def how_many_symmetries(self):
57 # temporary backwards compatibility for external calls
58 return len(self.symmetries)
60 def get_infostring(self):
61 # Maybe we can avoid calling this somehow, we're only using
62 # it to print:
63 K_gK = self.group_kpoints()
64 ng = len(K_gK)
65 txt = f'{ng} groups of equivalent kpoints. '
66 percent = (1. - (ng + 0.) / self.kd.nbzkpts) * 100
67 txt += f'{percent}% reduction.\n'
68 return txt
70 def group_kpoints(self, K_k=None):
71 """Group kpoints according to the reduced symmetries"""
72 if K_k is None:
73 K_k = np.arange(self.kd.nbzkpts)
74 bz2bz_kS = self.kd.bz2bz_ks # on kd, s is the global symmetry index
75 nk = len(bz2bz_kS)
76 sbz2sbz_ks = bz2bz_kS[K_k][:, self.symmetries.S_s] # s: q-symmetries
77 # Avoid -1 (see documentation in gpaw.symmetry)
78 sbz2sbz_ks[sbz2sbz_ks == -1] = nk
80 smallestk_k = np.sort(sbz2sbz_ks)[:, 0]
81 k2g_g = np.unique(smallestk_k, return_index=True)[1]
83 K_gs = sbz2sbz_ks[k2g_g]
84 K_gK = [np.unique(K_s[K_s != nk]) for K_s in K_gs]
86 return K_gK
88 def get_kpt_domain(self):
89 k_kc = np.array([self.kd.bzk_kc[K_K[0]] for
90 K_K in self.group_kpoints()])
91 return k_kc
93 def get_tetrahedron_ikpts(self, *, pbc_c, cell_cv):
94 """Find irreducible k-points for tetrahedron integration."""
95 U_scc = np.array([ # little group of q
96 sign * U_cc for U_cc, sign, _ in self.symmetries])
98 # Determine the irreducible BZ
99 bzk_kc, ibzk_kc, _ = get_reduced_bz(cell_cv,
100 U_scc,
101 False,
102 pbc_c=pbc_c)
104 n = 3
105 N_xc = np.indices((n, n, n)).reshape((3, n**3)).T - n // 2
107 # Find the irreducible kpoints
108 tess = Delaunay(ibzk_kc)
109 ik_kc = []
110 for N_c in N_xc:
111 k_kc = self.kd.bzk_kc + N_c
112 k_kc = k_kc[tess.find_simplex(k_kc) >= 0]
113 if not len(ik_kc) and len(k_kc):
114 ik_kc = unique_rows(k_kc)
115 elif len(k_kc):
116 ik_kc = unique_rows(np.append(k_kc, ik_kc, axis=0))
118 return ik_kc
120 def get_tetrahedron_kpt_domain(self, *, pbc_c, cell_cv):
121 ik_kc = self.get_tetrahedron_ikpts(pbc_c=pbc_c, cell_cv=cell_cv)
122 if pbc_c.all():
123 k_kc = ik_kc
124 else:
125 k_kc = np.append(ik_kc,
126 ik_kc + (~pbc_c).astype(int),
127 axis=0)
128 return k_kc
130 def get_kpoint_weight(self, k_c):
131 K = self.kptfinder.find(k_c)
132 iK = self.kd.bz2ibz_k[K]
133 K_k = self.unfold_ibz_kpoint(iK)
134 K_gK = self.group_kpoints(K_k)
136 for K_k in K_gK:
137 if K in K_k:
138 return len(K_k)
140 def unfold_ibz_kpoint(self, ik):
141 """Return kpoints related to irreducible kpoint."""
142 kd = self.kd
143 K_k = np.unique(kd.bz2bz_ks[kd.ibz2bz_k[ik]])
144 K_k = K_k[K_k != -1]
145 return K_k