|
| 1 | +#!/usr/bin/env python |
| 2 | +# Copyright 2014-2026 The PySCF Developers. All Rights Reserved. |
| 3 | +# |
| 4 | +# Licensed under the Apache License, Version 2.0 (the "License"); |
| 5 | +# you may not use this file except in compliance with the License. |
| 6 | +# You may obtain a copy of the License at |
| 7 | +# |
| 8 | +# http://www.apache.org/licenses/LICENSE-2.0 |
| 9 | +# |
| 10 | +# Unless required by applicable law or agreed to in writing, software |
| 11 | +# distributed under the License is distributed on an "AS IS" BASIS, |
| 12 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 13 | +# See the License for the specific language governing permissions and |
| 14 | +# limitations under the License. |
| 15 | + |
| 16 | +''' |
| 17 | +B-vector shells for finite-difference k-derivatives on a uniform mesh. |
| 18 | +
|
| 19 | +We need a set of nearest-neighbor offsets {b} whose weights {w_b} satisfy |
| 20 | +
|
| 21 | + sum_b w_b b_alpha b_beta = delta_{alpha, beta} (MV "B1") |
| 22 | +
|
| 23 | +Shells are added one at a time by distance and the weights solved by |
| 24 | +least squares until the residual drops below tolerance. Appendix B of |
| 25 | +Marzari & Vanderbilt, Phys. Rev. B 56, 12847 (1997). |
| 26 | +''' |
| 27 | + |
| 28 | +import numpy |
| 29 | + |
| 30 | + |
| 31 | +def find_bvectors(cell, kpts, *, kmesh=None, |
| 32 | + kmesh_tol=1e-6, b1_tol=1e-10, max_shells=12, search_range=3): |
| 33 | + '''Find the b-vector shells that satisfy B1 for this k-mesh. |
| 34 | +
|
| 35 | + Args: |
| 36 | + cell : pbc.gto.Cell |
| 37 | + kpts : (nk, 3) ndarray |
| 38 | + Cartesian k-points (1/Bohr). Must come from a uniform MP grid. |
| 39 | +
|
| 40 | + Kwargs: |
| 41 | + kmesh : (3,) or None |
| 42 | + MP grid size. Inferred from kpts when None. |
| 43 | + kmesh_tol : float |
| 44 | + Used for two related things: grouping b-vectors into shells of |
| 45 | + equal |b|, and matching k+b back to the k-point list. |
| 46 | + b1_tol : float |
| 47 | + How close the B1 residual has to get before we accept the shell set. |
| 48 | + max_shells, search_range : int |
| 49 | + Safety limits. search_range is the half-width of the integer |
| 50 | + lattice of candidate b-vectors (in units of the k-mesh spacing). |
| 51 | +
|
| 52 | + Returns: |
| 53 | + bvecs : (n_nn, 3) Cartesian b-vectors, 1/Bohr. |
| 54 | + weights : (n_nn,) per-b finite-difference weights, Bohr**2. |
| 55 | + kpb_idx : (nk, n_nn) int. The k-point index that k+b wraps to. |
| 56 | + kpb_g : (nk, n_nn, 3) int. The reciprocal-lattice shift G such that |
| 57 | + kpts[k] + bvecs[b] = kpts[kpb_idx[k, b]] + G @ recip. |
| 58 | + ''' |
| 59 | + kpts = numpy.asarray(kpts, dtype=float) |
| 60 | + if kpts.ndim != 2 or kpts.shape[1] != 3: |
| 61 | + raise ValueError(f'kpts must have shape (nk, 3); got {kpts.shape}') |
| 62 | + if kpts.shape[0] < 2: |
| 63 | + raise ValueError('find_bvectors requires at least 2 k-points') |
| 64 | + |
| 65 | + if kmesh is None: |
| 66 | + kmesh = _infer_kmesh(cell, kpts, kmesh_tol) |
| 67 | + else: |
| 68 | + kmesh = numpy.asarray(kmesh, dtype=int) |
| 69 | + if kmesh.shape != (3,): |
| 70 | + raise ValueError(f'kmesh must be a length-3 sequence; got shape {kmesh.shape}') |
| 71 | + if int(numpy.prod(kmesh)) != kpts.shape[0]: |
| 72 | + raise ValueError( |
| 73 | + f'kmesh={tuple(int(m) for m in kmesh)} has product ' |
| 74 | + f'{int(numpy.prod(kmesh))} but {kpts.shape[0]} k-points were given') |
| 75 | + |
| 76 | + recip = cell.reciprocal_vectors() |
| 77 | + shells = _candidate_shells(recip, kmesh, kmesh_tol, max_shells, search_range) |
| 78 | + |
| 79 | + bvecs, weights = _select_shells(shells, b1_tol) |
| 80 | + kpb_idx, kpb_g = _build_connectivity(recip, kpts, bvecs, kmesh_tol) |
| 81 | + return bvecs, weights, kpb_idx, kpb_g |
| 82 | + |
| 83 | + |
| 84 | +def _infer_kmesh(cell, kpts, tol): |
| 85 | + '''Read back the MP grid from a list of k-points.''' |
| 86 | + frac = cell.get_scaled_kpts(kpts) |
| 87 | + # Wrap into [0, 1). The +tol nudge keeps tiny negative values |
| 88 | + # (round-off below zero) from wrapping to ~1. |
| 89 | + frac = frac - numpy.floor(frac + tol) |
| 90 | + mp = numpy.zeros(3, dtype=int) |
| 91 | + for alpha in range(3): |
| 92 | + vals = numpy.sort(frac[:, alpha]) |
| 93 | + count = 1 |
| 94 | + for i in range(1, len(vals)): |
| 95 | + if vals[i] - vals[i - 1] > tol: |
| 96 | + count += 1 |
| 97 | + mp[alpha] = count |
| 98 | + if int(numpy.prod(mp)) != kpts.shape[0]: |
| 99 | + raise ValueError( |
| 100 | + f'Cannot infer uniform MP grid: guessed {tuple(int(m) for m in mp)} but ' |
| 101 | + f'have {kpts.shape[0]} k-points. Pass kmesh= explicitly.') |
| 102 | + return mp |
| 103 | + |
| 104 | + |
| 105 | +def _candidate_shells(recip, kmesh, tol, max_shells, search_range): |
| 106 | + '''Enumerate candidate b-vectors, grouped into shells of equal |b|.''' |
| 107 | + db = recip / kmesh[:, None] # spacing of the k-point sublattice in 1/Bohr |
| 108 | + grid = numpy.arange(-search_range, search_range + 1) |
| 109 | + ijk = numpy.array(numpy.meshgrid(grid, grid, grid, indexing='ij')).reshape(3, -1).T |
| 110 | + ijk = ijk[numpy.any(ijk != 0, axis=1)] |
| 111 | + bvecs = ijk @ db |
| 112 | + norms = numpy.linalg.norm(bvecs, axis=1) |
| 113 | + |
| 114 | + order = numpy.argsort(norms) |
| 115 | + bvecs = bvecs[order] |
| 116 | + norms = norms[order] |
| 117 | + |
| 118 | + shells = [] |
| 119 | + i = 0 |
| 120 | + n = len(norms) |
| 121 | + while i < n and len(shells) < max_shells: |
| 122 | + j = i |
| 123 | + while j < n and norms[j] - norms[i] < tol: |
| 124 | + j += 1 |
| 125 | + shells.append(bvecs[i:j].copy()) |
| 126 | + i = j |
| 127 | + return shells |
| 128 | + |
| 129 | + |
| 130 | +def _select_shells(shells, tol): |
| 131 | + '''Add shells one at a time and solve for weights until B1 holds.''' |
| 132 | + target = numpy.array([1.0, 1.0, 1.0, 0.0, 0.0, 0.0]) |
| 133 | + for n in range(1, len(shells) + 1): |
| 134 | + A = numpy.zeros((n, 6)) |
| 135 | + for s in range(n): |
| 136 | + shell = shells[s] |
| 137 | + A[s, 0] = (shell[:, 0] ** 2).sum() |
| 138 | + A[s, 1] = (shell[:, 1] ** 2).sum() |
| 139 | + A[s, 2] = (shell[:, 2] ** 2).sum() |
| 140 | + A[s, 3] = (shell[:, 0] * shell[:, 1]).sum() |
| 141 | + A[s, 4] = (shell[:, 0] * shell[:, 2]).sum() |
| 142 | + A[s, 5] = (shell[:, 1] * shell[:, 2]).sum() |
| 143 | + w, *_ = numpy.linalg.lstsq(A.T, target, rcond=None) |
| 144 | + if numpy.max(numpy.abs(A.T @ w - target)) < tol: |
| 145 | + bvecs = numpy.concatenate(shells[:n], axis=0) |
| 146 | + weights = numpy.concatenate([ |
| 147 | + numpy.full(len(shells[s]), w[s]) for s in range(n) |
| 148 | + ]) |
| 149 | + return bvecs, weights |
| 150 | + raise RuntimeError( |
| 151 | + f'B1 completeness relation not satisfied within {len(shells)} shells. ' |
| 152 | + f'Increase max_shells or search_range.') |
| 153 | + |
| 154 | + |
| 155 | +def _build_connectivity(recip, kpts, bvecs, tol): |
| 156 | + '''For each (k, b), find which k-point k+b wraps to and the G shift.''' |
| 157 | + inv_recip = numpy.linalg.inv(recip) |
| 158 | + frac_kpts = kpts @ inv_recip |
| 159 | + frac_kpts = frac_kpts - numpy.floor(frac_kpts + tol) |
| 160 | + frac_b = bvecs @ inv_recip |
| 161 | + |
| 162 | + nk = len(kpts) |
| 163 | + nb = len(bvecs) |
| 164 | + kpb_idx = numpy.zeros((nk, nb), dtype=int) |
| 165 | + kpb_g = numpy.zeros((nk, nb, 3), dtype=int) |
| 166 | + |
| 167 | + for i in range(nk): |
| 168 | + for b in range(nb): |
| 169 | + diff = (frac_kpts[i] + frac_b[b])[None, :] - frac_kpts |
| 170 | + g = numpy.round(diff) |
| 171 | + residual = numpy.linalg.norm(diff - g, axis=1) |
| 172 | + j = int(numpy.argmin(residual)) |
| 173 | + if residual[j] > tol: |
| 174 | + raise ValueError( |
| 175 | + f'k+b not matched to any k-point for k={i}, b={b} ' |
| 176 | + f'(min residual {residual[j]:.3e} > tol {tol:.3e})') |
| 177 | + kpb_idx[i, b] = j |
| 178 | + kpb_g[i, b] = g[j] |
| 179 | + return kpb_idx, kpb_g |
0 commit comments