Coverage for gpaw/symmetry.py: 88%

3# Wits University

4# Please see the accompanying LICENSE file for further information.

5from typing import Tuple

7from ase.utils import gcd

8import numpy as np

10import gpaw.cgpaw as cgpaw

11import gpaw.mpi as mpi

14def frac(f: float,

15 n: int = 2 * 3 * 4 * 5,

16 tol: float = 1e-6) -> Tuple[int, int]:

17 """Convert to fraction.

19 >>> frac(0.5)

20 (1, 2)

21 """

22 if f == 0:

23 return 0, 1

24 x = n * f

25 if abs(x - round(x)) > n * tol:

26 raise ValueError

27 x = int(round(x))

28 d = gcd(x, n)

29 return x // d, n // d

32def sfrac(f: float) -> str:

33 """Format as fraction.

35 >>> sfrac(0.5)

36 '1/2'

37 >>> sfrac(2 / 3)

38 '2/3'

39 >>> sfrac(0)

40 '0'

41 """

42 if f == 0:

43 return '0'

44 return '%d/%d' % frac(f)

47class Symmetry:

48 """Interface class for determination of symmetry, point and space groups.

50 It also provides to apply symmetry operations to kpoint grids,

51 wavefunctions and forces.

52 """

53 def __init__(self, id_a, cell_cv, pbc_c=np.ones(3, bool), tolerance=1e-7,

54 point_group=True, time_reversal=True, symmorphic=True,

55 allow_invert_aperiodic_axes=True):

56 """Construct symmetry object.

58 Parameters:

60 id_a: list of int

61 Numbered atomic types

62 cell_cv: array(3,3), float

63 Cartesian lattice vectors

64 pbc_c: array(3), bool

65 Periodic boundary conditions.

66 tolerance: float

67 Tolerance for symmetry determination.

68 symmorphic: bool

69 Switch for the use of non-symmorphic symmetries aka: symmetries

70 with fractional translations. Default is to use only symmorphic

71 symmetries.

72 point_group: bool

73 Use point-group symmetries.

74 time_reversal: bool

75 Use time-reversal symmetry.

76 tolerance: float

77 Relative tolerance.

79 Attributes:

81 op_scc:

82 Array of rotation matrices

83 ft_sc:

84 Array of fractional translation vectors

85 a_sa:

86 Array of atomic indices after symmetry operation

87 has_inversion:

88 (bool) Have inversion

89 """

91 self.id_a = id_a

92 self.cell_cv = np.array(cell_cv, float)

93 assert self.cell_cv.shape == (3, 3)

94 self.pbc_c = np.array(pbc_c, bool)

95 self.tol = tolerance

96 self.symmorphic = symmorphic

97 self.point_group = point_group

98 self.time_reversal = time_reversal

100 self.op_scc = np.identity(3, int).reshape((1, 3, 3))

101 self.ft_sc = np.zeros((1, 3))

102 self.a_sa = np.arange(len(id_a)).reshape((1, -1))

103 self.has_inversion = False

104 self.gcd_c = np.ones(3, int)

105

106 # For reading old gpw-files:

107 self.allow_invert_aperiodic_axes = allow_invert_aperiodic_axes

108

109 def analyze(self, spos_ac):

110 """Determine list of symmetry operations.

111

112 First determine all symmetry operations of the cell. Then call

113 ``prune_symmetries`` to remove those symmetries that are not satisfied

114 by the atoms.

115

116 It is not mandatory to call this method. If not called, only

117 time reversal symmetry may be used.

118 """

119 if self.point_group:

120 self.find_lattice_symmetry()

121 self.prune_symmetries_atoms(spos_ac)

122

123 def find_lattice_symmetry(self):

124 """Determine list of symmetry operations."""

125 # Symmetry operations as matrices in 123 basis.

126 # Operation is a 3x3 matrix, with possible elements -1, 0, 1, thus

127 # there are 3**9 = 19683 possible matrices:

128 combinations = 1 - np.indices([3] * 9)

129 U_scc = combinations.reshape((3, 3, 3**9)).transpose((2, 0, 1))

130

131 # The metric of the cell should be conserved after applying

132 # the operation:

133 metric_cc = self.cell_cv.dot(self.cell_cv.T)

134 metric_scc = np.einsum('sij, jk, slk -> sil',

135 U_scc, metric_cc, U_scc,

136 optimize=True)

137 mask_s = abs(metric_scc - metric_cc).sum(2).sum(1) <= self.tol

138 U_scc = U_scc[mask_s]

139

140 # Operation must not swap axes that don't have same PBC:

141 pbc_cc = np.logical_xor.outer(self.pbc_c, self.pbc_c)

142 mask_s = ~U_scc[:, pbc_cc].any(axis=1)

143 U_scc = U_scc[mask_s]

144

145 if not self.allow_invert_aperiodic_axes:

146 # Operation must not invert axes that are not periodic:

147 mask_s = (U_scc[:, np.diag(~self.pbc_c)] == 1).all(axis=1)

148 U_scc = U_scc[mask_s]

149

150 self.op_scc = U_scc

151 self.ft_sc = np.zeros((len(self.op_scc), 3))

152

153 def prune_symmetries_atoms(self, spos_ac):

154 """Remove symmetries that are not satisfied by the atoms."""

155

156 if len(spos_ac) == 0:

157 self.a_sa = np.zeros((len(self.op_scc), 0), int)

158 return

159

160 # Build lists of atom numbers for each type of atom - one

161 # list for each combination of atomic number, setup type,

162 # magnetic moment and basis set:

163 a_ij = {}

164 for a, id in enumerate(self.id_a):

165 if id in a_ij:

166 a_ij[id].append(a)

167 else:

168 a_ij[id] = [a]

169

170 a_j = a_ij[self.id_a[0]] # just pick the first species

171

172 # if supercell disable fractional translations:

173 if not self.symmorphic:

174 op_cc = np.identity(3, int)

175 ftrans_sc = spos_ac[a_j[1:]] - spos_ac[a_j[0]]

176 ftrans_sc -= np.rint(ftrans_sc)

177 for ft_c in ftrans_sc:

178 a_a = self.check_one_symmetry(spos_ac, op_cc, ft_c, a_ij)

179 if a_a is not None:

180 self.symmorphic = True

181 break

182

183 symmetries = []

184 ftsymmetries = []

185

186 # go through all possible symmetry operations

187 for op_cc in self.op_scc:

188 # first ignore fractional translations

189 a_a = self.check_one_symmetry(spos_ac, op_cc, [0, 0, 0], a_ij)

190 if a_a is not None:

191 symmetries.append((op_cc, [0, 0, 0], a_a))

192 elif not self.symmorphic:

193 # check fractional translations

194 sposrot_ac = np.dot(spos_ac, op_cc)

195 ftrans_jc = sposrot_ac[a_j] - spos_ac[a_j[0]]

196 ftrans_jc -= np.rint(ftrans_jc)

197 for ft_c in ftrans_jc:

198 try:

199 nom_c, denom_c = np.array([frac(ft, tol=self.tol)

200 for ft in ft_c]).T

201 except ValueError:

202 continue

203 ft_c = nom_c / denom_c

204 a_a = self.check_one_symmetry(spos_ac, op_cc, ft_c, a_ij)

205 if a_a is not None:

206 ftsymmetries.append((op_cc, ft_c, a_a))

207 for c, d in enumerate(denom_c):

208 if self.gcd_c[c] % d != 0:

209 self.gcd_c[c] *= d

210

211 # Add symmetry operations with fractional translations at the end:

212 symmetries.extend(ftsymmetries)

213 self.op_scc = np.array([sym[0] for sym in symmetries])

214 self.ft_sc = np.array([sym[1] for sym in symmetries])

215 self.a_sa = np.array([sym[2] for sym in symmetries])

216

217 inv_cc = -np.eye(3, dtype=int)

218 self.has_inversion = (self.op_scc == inv_cc).all(2).all(1).any()

219

220 def check_one_symmetry(self, spos_ac, op_cc, ft_c, a_ij):

221 """Checks whether atoms satisfy one given symmetry operation."""

222

223 a_a = np.zeros(len(spos_ac), int)

224 for a_j in a_ij.values():

225 spos_jc = spos_ac[a_j]

226 for a in a_j:

227 spos_c = np.dot(spos_ac[a], op_cc)

228 sdiff_jc = spos_c - spos_jc - ft_c

229 sdiff_jc -= sdiff_jc.round()

230 indices = np.where(abs(sdiff_jc).max(1) < self.tol)[0]

231 if len(indices) == 1:

232 j = indices[0]

233 a_a[a] = a_j[j]

234 else:

235 assert len(indices) == 0

236 return

237

238 return a_a

239

240 def check(self, spos_ac):

241 """Check if positions satisfy symmetry operations."""

242

243 nsymold = len(self.op_scc)

244 self.prune_symmetries_atoms(spos_ac)

245 if len(self.op_scc) < nsymold:

246 raise RuntimeError('Broken symmetry!')

247

248 def reduce(self, bzk_kc, comm=None):

249 """Reduce k-points to irreducible part of the BZ.

250

251 Returns the irreducible k-points and the weights and other stuff.

252

253 """

254 return reduce_kpts(bzk_kc,

255 self.op_scc,

256 self.time_reversal and not self.has_inversion,

257 comm,

258 self.tol)

259

260 def check_grid(self, N_c) -> bool:

261 """Check that symmetries are commensurate with grid."""

262 for s, (U_cc, ft_c) in enumerate(zip(self.op_scc, self.ft_sc)):

263 t_c = ft_c * N_c

264 # Make sure all grid-points map onto another grid-point:

265 if (((N_c * U_cc).T % N_c).any() or

266 not np.allclose(t_c, t_c.round())):

267 return False

268 return True

269

270 def symmetrize(self, a, gd):

271 """Symmetrize array."""

272 gd.symmetrize(a, self.op_scc, self.ft_sc)

273

274 def symmetrize_positions(self, spos_ac):

275 """Symmetrizes the atomic positions."""

276 spos_tmp_ac = np.zeros_like(spos_ac)

277 spos_new_ac = np.zeros_like(spos_ac)

278 for i, op_cc in enumerate(self.op_scc):

279 spos_tmp_ac[:] = 0.

280 for a in range(len(spos_ac)):

281 spos_c = np.dot(spos_ac[a], op_cc) - self.ft_sc[i]

282 # Bring back the negative ones:

283 spos_c = spos_c - np.floor(spos_c + 1e-5)

284 spos_tmp_ac[self.a_sa[i][a]] += spos_c

285 spos_new_ac += spos_tmp_ac

286

287 spos_new_ac /= len(self.op_scc)

288 return spos_new_ac

289

290 def symmetrize_wavefunction(self, a_g, kibz_c, kbz_c, op_cc,

291 time_reversal):

292 """Generate Bloch function from symmetry related function in the IBZ.

293

294 a_g: ndarray

295 Array with Bloch function from the irreducible BZ.

296 kibz_c: ndarray

297 Corresponding k-point coordinates.

298 kbz_c: ndarray

299 K-point coordinates of the symmetry related k-point.

300 op_cc: ndarray

301 Point group operation connecting the two k-points.

302 time-reversal: bool

303 Time-reversal symmetry required in addition to the point group

304 symmetry to connect the two k-points.

305 """

306

307 # Identity

308 if (np.abs(op_cc - np.eye(3, dtype=int)) < 1e-10).all():

309 if time_reversal:

310 return a_g.conj()

311 else:

312 return a_g

313 # Inversion symmetry

314 elif (np.abs(op_cc + np.eye(3, dtype=int)) < 1e-10).all():

315 return a_g.conj()

316 # General point group symmetry

317 else:

318 import gpaw.cgpaw as cgpaw

319 b_g = np.zeros_like(a_g)

320 if time_reversal:

321 # assert abs(np.dot(op_cc, kibz_c) - -kbz_c) < tol

322 cgpaw.symmetrize_wavefunction(a_g, b_g, op_cc.T.copy(),

323 kibz_c, -kbz_c)

324 return b_g.conj()

325 else:

326 # assert abs(np.dot(op_cc, kibz_c) - kbz_c) < tol

327 cgpaw.symmetrize_wavefunction(a_g, b_g, op_cc.T.copy(),

328 kibz_c, kbz_c)

329 return b_g

330

331 def symmetrize_forces(self, F0_av):

332 """Symmetrize forces."""

333 F_av = np.zeros_like(F0_av)

334 for map_a, op_cc in zip(self.a_sa, self.op_scc):

335 op_vv = np.dot(np.linalg.inv(self.cell_cv),

336 np.dot(op_cc, self.cell_cv))

337 for a1, a2 in enumerate(map_a):

338 F_av[a2] += np.dot(F0_av[a1], op_vv)

339 return F_av / len(self.op_scc)

340

341 def __str__(self):

342 n = len(self.op_scc)

343 nft = self.ft_sc.any(1).sum()

344 lines = [f'Symmetries present (total): {n}']

345 if not self.symmorphic:

346 lines.append(

347 f'Symmetries with fractional translations: {nft}')

348

349 # X-Y grid of symmetry matrices:

350

351 lines.append('')

352 nx = 6 if self.symmorphic else 3

353 ns = len(self.op_scc)

354 y = 0

355 for y in range((ns + nx - 1) // nx):

356 for c in range(3):

357 line = ''

358 for x in range(nx):

359 s = x + y * nx

360 if s == ns:

361 break

362 op_c = self.op_scc[s, c]

363 ft = self.ft_sc[s, c]

364 line += ' (%2d %2d %2d)' % tuple(op_c)

365 if not self.symmorphic:

366 line += ' + (%4s)' % sfrac(ft)

367 lines.append(line)

368 lines.append('')

369 return '\n'.join(lines)

370

371

372def map_k_points(bzk_kc, U_scc, time_reversal, comm=None, tol=1e-11):

373 """Find symmetry relations between k-points.

374

375 This is a Python-wrapper for a C-function that does the hard work

376 which is distributed over comm.

377

378 The map bz2bz_ks is returned. If there is a k2 for which::

379

380 = _ _ _

381 U q = q + N,

382 s k1 k2

383

384 where N is a vector of integers, then bz2bz_ks[k1, s] = k2, otherwise

385 if there is a k2 for which::

386

387 = _ _ _

388 U q = -q + N,

389 s k1 k2

390

391 then bz2bz_ks[k1, s + nsym] = k2, where nsym = len(U_scc). Otherwise

392 bz2bz_ks[k1, s] = -1.

393 """

394

395 if comm is None or isinstance(comm, mpi.DryRunCommunicator):

396 comm = mpi.serial_comm

397

398 nbzkpts = len(bzk_kc)

399 ka = nbzkpts * comm.rank // comm.size

400 kb = nbzkpts * (comm.rank + 1) // comm.size

401 assert comm.sum_scalar(kb - ka) == nbzkpts

402

403 if time_reversal:

404 U_scc = np.concatenate([U_scc, -U_scc])

405

406 bz2bz_ks = np.zeros((nbzkpts, len(U_scc)), int)

407 bz2bz_ks[ka:kb] = -1

408 cgpaw.map_k_points(np.ascontiguousarray(bzk_kc),

409 np.ascontiguousarray(U_scc), tol, bz2bz_ks, ka, kb)

410 comm.sum(bz2bz_ks)

411 return bz2bz_ks

412

413

414def reduce_kpts(bzk_kc,

415 U_scc,

416 use_time_reversal,

417 comm,

418 tol):

419 nbzkpts = len(bzk_kc)

420 nsym = len(U_scc)

421

422 bz2bz_ks = map_k_points_fast(bzk_kc, U_scc, use_time_reversal,

423 comm, tol)

424

425 bz2bz_k = -np.ones(nbzkpts + 1, int)

426 ibz2bz_k = []

427 for k in range(nbzkpts - 1, -1, -1):

428 # Reverse order looks more natural

429 if bz2bz_k[k] == -1:

430 bz2bz_k[bz2bz_ks[k]] = k

431 ibz2bz_k.append(k)

432 ibz2bz_k = np.array(ibz2bz_k[::-1])

433 bz2bz_k = bz2bz_k[:-1].copy()

434

435 bz2ibz_k = np.empty(nbzkpts, int)

436 bz2ibz_k[ibz2bz_k] = np.arange(len(ibz2bz_k))

437 bz2ibz_k = bz2ibz_k[bz2bz_k]

438

439 weight_k = np.bincount(bz2ibz_k) * (1.0 / nbzkpts)

440

441 # Symmetry operation mapping IBZ to BZ:

442 sym_k = np.empty(nbzkpts, int)

443 for k in range(nbzkpts):

444 # We pick the first one found:

445 try:

446 sym_k[k] = np.where(bz2bz_ks[bz2bz_k[k]] == k)[0][0]

447 except IndexError:

448 print(nbzkpts)

449 print(k)

450 print(bz2bz_k)

451 print(bz2bz_ks[bz2bz_k[k]])

452 print(np.shape(np.where(bz2bz_ks[bz2bz_k[k]] == k)))

453 print(bz2bz_k[k])

454 print(bz2bz_ks[bz2bz_k[k]] == k)

455 raise

456

457 # Time-reversal symmetry used on top of the point group operation:

458 if use_time_reversal:

459 time_reversal_k = sym_k >= nsym

460 sym_k %= nsym

461 else:

462 time_reversal_k = np.zeros(nbzkpts, bool)

463

464 assert (ibz2bz_k[bz2ibz_k] == bz2bz_k).all()

465 for k in range(nbzkpts):

466 sign = 1 - 2 * time_reversal_k[k]

467 dq_c = (np.dot(U_scc[sym_k[k]], bzk_kc[bz2bz_k[k]]) -

468 sign * bzk_kc[k])

469 dq_c -= dq_c.round()

470 assert abs(dq_c).max() < 1e-10

471

472 return (bzk_kc[ibz2bz_k], weight_k,

473 sym_k, time_reversal_k, bz2ibz_k, ibz2bz_k, bz2bz_ks)

474

475

476def map_k_points_fast(bzk_kc, U_scc, time_reversal, comm=None, tol=1e-7):

477 """Find symmetry relations between k-points.

478

479 Performs the same task as map_k_points(), but much faster.

480 This is achieved by finding the symmetry related kpoints using

481 lexical sorting instead of brute force searching.

482

483 bzk_kc: ndarray

484 kpoint coordinates.

485 U_scc: ndarray

486 Symmetry operations

487 time_reversal: Bool

488 Use time reversal symmetry in mapping.

489 comm:

490 Communicator

491 tol: float

492 When kpoint are closer than tol, they are

493 considered to be identical.

494 """

495

496 nbzkpts = len(bzk_kc)

497

498 if time_reversal:

499 U_scc = np.concatenate([U_scc, -U_scc])

500

501 bz2bz_ks = np.zeros((nbzkpts, len(U_scc)), int)

502 bz2bz_ks[:] = -1

503

504 for s, U_cc in enumerate(U_scc):

505 # Find mapped kpoints

506 Ubzk_kc = np.dot(bzk_kc, U_cc.T)

507

508 # Do some work on the input

509 k_kc = np.concatenate([bzk_kc, Ubzk_kc])

510 k_kc = np.mod(np.mod(k_kc, 1), 1)

511 aglomerate_points(k_kc, tol)

512 k_kc = k_kc.round(-np.log10(tol).astype(int))

513 k_kc = np.mod(k_kc, 1)

514

515 # Find the lexicographical order

516 order = np.lexsort(k_kc.T)

517 k_kc = k_kc[order]

518 diff_kc = np.diff(k_kc, axis=0)

519 equivalentpairs_k = np.array((diff_kc == 0).all(1),

520 bool)

521

522 # Mapping array.

523 orders = np.array([order[:-1][equivalentpairs_k],

524 order[1:][equivalentpairs_k]])

525

526 # This has to be true.

527 assert (orders[0] < nbzkpts).all()

528 assert (orders[1] >= nbzkpts).all()

529 bz2bz_ks[orders[1] - nbzkpts, s] = orders[0]

530

531 return bz2bz_ks

532

533

534def aglomerate_points(k_kc, tol):

535 nd = k_kc.shape[1]

536 nbzkpts = len(k_kc)

537 inds_kc = np.argsort(k_kc, axis=0)

538 for c in range(nd):

539 sk_k = k_kc[inds_kc[:, c], c]

540 dk_k = np.diff(sk_k)

541

542 # Partition the kpoints into groups

543 pt_K = np.argwhere(dk_k > tol)[:, 0]

544 pt_K = np.append(np.append(0, pt_K + 1), 2 * nbzkpts)

545 for i in range(len(pt_K) - 1):

546 k_kc[inds_kc[pt_K[i]:pt_K[i + 1], c],

547 c] = k_kc[inds_kc[pt_K[i], c], c]

548

549

550def atoms2symmetry(atoms, id_a=None, tolerance=1e-7):

551 """Create symmetry object from atoms object."""

552 if id_a is None:

553 id_a = atoms.get_atomic_numbers()

554 symmetry = Symmetry(id_a, atoms.cell, atoms.pbc,

555 symmorphic=False,

556 time_reversal=False,

557 tolerance=tolerance)

558 symmetry.analyze(atoms.get_scaled_positions())

559 return symmetry

560

561

562class CLICommand:

563 """Analyse symmetry (and show IBZ k-points).

564

565 Example:

566

567 $ ase build -x bcc -a 3.5 Li | gpaw symmetry -k "{density:3,gamma:1}"

568 symmetry:

569 number of symmetries: 48

570 number of symmetries with translation: 0

571

572 bz sampling:

573 number of bz points: 512

574 number of ibz points: 29

575 monkhorst-pack size: [8, 8, 8]

576 monkhorst-pack shift: [0.0625, 0.0625, 0.0625]

577 """

578

579 @staticmethod

580 def add_arguments(parser):

581 parser.add_argument('-t', '--tolerance', type=float, default=1e-7,

582 help='Tolerance used for identifying symmetries.')

583 parser.add_argument(

584 '-k', '--k-points',

585 help='Use symmetries to reduce number of k-points. '

586 'Exapmples: "4,4,4", "{density:3.5,gamma:True}".')

587 parser.add_argument('-v', '--verbose', action='store_true',

588 help='Show symmetry operations (and k-points).')

589 parser.add_argument('-s', '--symmorphic', action='store_true',

590 help='Only find symmorphic symmetries.')

591 parser.add_argument('filename', nargs='?', default='-',

592 help='Filename to read structure from. '

593 'Use "-" for reading from stdin. '

594 'Default is "-".')

595

596 @staticmethod

597 def run(args):

598 import sys

599 from gpaw.new.symmetry import create_symmetries_object

600 from gpaw.dft import MonkhorstPack

601 from ase.cli.run import str2dict

602 from ase.db import connect

603 from ase.io import read

604

605 if args.filename == '-':

606 atoms = next(connect(sys.stdin).select()).toatoms()

607 else:

608 atoms = read(args.filename)

609 symmetries = create_symmetries_object(

610 atoms,

611 tolerance=args.tolerance,

612 symmorphic=args.symmorphic)

613 txt = str(symmetries)

614 if not args.verbose:

615 txt = txt.split(' rotations', 1)[0]

616 print(txt)

617 if args.k_points:

618 k = str2dict('kpts=' + args.k_points)['kpts']

619 bz = MonkhorstPack.from_param(k).build(atoms)

620 ibz = bz.reduce(symmetries)

621 txt = str(ibz)

622 if not args.verbose:

623 txt = txt.split(' points and weights:', 1)[0]

624 print(txt)