From d3caaf08ff4f0b95c7543d882e42b4730bad32d8 Mon Sep 17 00:00:00 2001 From: James Nelson Date: Fri, 22 Aug 2025 08:34:06 +0100 Subject: [PATCH 1/2] rm cell --- qse/qbits.py | 223 ++++--------------------------------- tests/qse/lattices_test.py | 22 ++-- 2 files changed, 30 insertions(+), 215 deletions(-) diff --git a/qse/qbits.py b/qse/qbits.py index 63441581..e48c5047 100644 --- a/qse/qbits.py +++ b/qse/qbits.py @@ -7,12 +7,10 @@ from math import cos, sin import numpy as np -from ase.cell import Cell from ase.geometry import ( find_mic, get_dihedrals, ) -from ase.utils import deprecated from qse.qbit import Qbit from qse.visualise import draw as _draw @@ -39,15 +37,10 @@ class Qbits: scaled_positions: list of scaled-positions Like positions, but given in units of the unit cell. Can not be set at the same time as positions. - cell: 3x3 matrix or length 3 or 6 vector - Unit cell vectors. Can also be given as just three - numbers for orthorhombic cells, or 6 numbers, where - first three are lengths of unit cell vectors, and the - other three are angles between them (in degrees), in following order: - [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)]. + cell: 3x3 matrix, optional + The unit cell vectors. First vector will lie in x-direction, second in xy-plane, and the third one in z-positive subspace. - Default value: [0, 0, 0]. celldisp: Vector Unit cell displacement vector. To visualize a displaced cell around the center of mass of a Systems of qbits. Default value @@ -95,8 +88,7 @@ class Qbits: ... [[0, 0, 0], ... [0.5, 0.5, 0.5]]) >>> qdim = qse.Qbits(positions=xd) - >>> qdim.cell = [1,1,1] - >>> qdim.pbc = True + >>> qdim.cell = np.eye(3) >>> qlat = qdim.repeat([3,3,3]) The qdim will have shape = (2,1,1) and qlat will have shape = (6, 3, 3) @@ -159,10 +151,8 @@ def __init__( self.new_array("labels", labels, ">> qbits = Qbits('He') - >>> a, b, c = 7, 7.5, 8 - >>> qbits.set_cell([a, b, c]) - >>> qbits.set_cell([(a, 0, 0), (0, b, 0), (0, 0, c)]) - - FCC unit cell: - - >>> qbits.set_cell([(0, b, b), (b, 0, b), (b, b, 0)]) - - Hexagonal unit cell: - - >>> qbits.set_cell([a, a, c, 90, 90, 120]) - - Rhombohedral unit cell: - - >>> alpha = 77 - >>> qbits.set_cell([a, a, a, alpha, alpha, alpha]) - """ - - # Override pbcs if and only if given a Cell object: - # print('here', cell) - cell = Cell.new(cell) - - # XXX not working well during initialize due to missing _constraints - if apply_constraint and hasattr(self, "_constraints"): - for constraint in self.constraints: - if hasattr(constraint, "adjust_cell"): - constraint.adjust_cell(self, cell) - - if scale_qbits: - M = np.linalg.solve(self.cell.complete(), cell.complete()) - self.positions[:] = np.dot(self.positions, M) - - self.cell[:] = cell - def set_celldisp(self, celldisp): """Set the unit cell displacement vectors.""" celldisp = np.array(celldisp, float) @@ -326,37 +258,6 @@ def get_celldisp(self): """Get the unit cell displacement vectors.""" return self._celldisp.copy() - def get_cell(self, complete=False): - """Get the three unit cell vectors as a `class`:ase.cell.Cell` object. - - The Cell object resembles a 3x3 ndarray, and cell[i, j] - is the jth Cartesian coordinate of the ith cell vector.""" - if complete: - return self.cell.complete() - return self.cell.copy() - - @deprecated("Please use qbits.cell.cellpar() instead") - def get_cell_lengths_and_angles(self): - """Get unit cell parameters. Sequence of 6 numbers. - - First three are unit cell vector lengths and second three - are angles between them:: - - [len(a), len(b), len(c), angle(b,c), angle(a,c), angle(a,b)] - - in degrees. - """ - return self.cell.cellpar() - - @deprecated("Please use qbits.cell.reciprocal()") - def get_reciprocal_cell(self): - """Get the three reciprocal lattice vectors as a 3x3 ndarray. - - Note that the commonly used factor of 2 pi for Fourier - transforms is not included here.""" - - return self.cell.reciprocal() - @property def nqbits(self): return len(self) @@ -482,7 +383,7 @@ def todict(self): d["labels"] = self.arrays["labels"] d["positions"] = self.arrays["positions"] d["states"] = self.arrays["states"] - d["cell"] = self.cell # np.asarray(self.cell) + d["cell"] = self.cell d["pbc"] = self.pbc if self._celldisp.any(): d["celldisp"] = self._celldisp @@ -549,12 +450,8 @@ def __repr__(self): tokens.append(txt + ",\n") cell = self.cell - if cell: - if cell.orthorhombic: - cell = cell.lengths().tolist() - else: - cell = cell.tolist() - tokens.append("cell={0}".format(cell) + ",\n") + if cell is not None: + tokens.append("cell={0}".format(list(cell)) + ",\n") if self._calc is not None: tokens.append("calculator={0}".format(self._calc.__class__.__name__)) txt = "{0}(\n{1}{2})".format( @@ -688,9 +585,12 @@ def __imul__(self, m): if isinstance(m, int): m = (m, m, m) + if self.cell is None: + raise ValueError("Cannot repeat without a defined unit cell.") + for x, vec in zip(m, self.cell): if x != 1 and not vec.any(): - raise ValueError("Cannot repeat along undefined lattice " "vector") + raise ValueError("Cannot repeat without a defined unit cell.") M = np.prod(m) n = len(self) @@ -775,87 +675,6 @@ def translate(self, displacement): """ self.positions += np.array(displacement) - def center_in_unit_cell(self, vacuum=None, axis=(0, 1, 2), about=None): - """ - Center qbits in unit cell. - - Centers the qbits in the unit cell, so there is the same - amount of vacuum on all sides. - - vacuum: float (default: None) - If specified adjust the amount of vacuum when centering. - If vacuum=10.0 there will thus be 10 Angstrom of vacuum - on each side. - axis: int or sequence of ints - Axis or axes to act on. Default: Act on all axes. - about: float or array (default: None) - If specified, center the qbits about . - I.e., about=(0., 0., 0.) (or just "about=0.", interpreted - identically), to center about the origin. - """ - - # Find the orientations of the faces of the unit cell - cell = self.cell.complete() - dirs = np.zeros_like(cell) - - lengths = cell.lengths() - for i in range(3): - dirs[i] = np.cross(cell[i - 1], cell[i - 2]) - dirs[i] /= np.linalg.norm(dirs[i]) - if dirs[i] @ cell[i] < 0.0: - dirs[i] *= -1 - - if isinstance(axis, int): - axes = (axis,) - else: - axes = axis - - # Now, decide how much each basis vector should be made longer - pos = self.positions - longer = np.zeros(3) - shift = np.zeros(3) - for i in axes: - if len(pos): - scalarprod = pos @ dirs[i] - p0 = scalarprod.min() - p1 = scalarprod.max() - else: - p0 = 0 - p1 = 0 - height = cell[i] @ dirs[i] - if vacuum is not None: - lng = (p1 - p0 + 2 * vacuum) - height - else: - lng = 0.0 # Do not change unit cell size! - top = lng + height - p1 - shf = 0.5 * (top - p0) - cosphi = cell[i] @ dirs[i] / lengths[i] - longer[i] = lng / cosphi - shift[i] = shf / cosphi - - # Now, do it! - translation = np.zeros(3) - for i in axes: - nowlen = lengths[i] - if vacuum is not None: - self.cell[i] = cell[i] * (1 + longer[i] / nowlen) - translation += shift[i] * cell[i] / nowlen - - # We calculated translations using the completed cell, - # so directions without cell vectors will have been centered - # along a "fake" vector of length 1. - # Therefore, we adjust by -0.5: - if not any(self.cell[i]): - translation[i] -= 0.5 - - # Optionally, translate to center about a point in space. - if about is not None: - for vector in self.cell: - translation -= vector / 2.0 - translation += about - - self.positions += translation - def get_centroid(self): r""" Get the centroid of the positions. @@ -984,13 +803,11 @@ def rotate(self, a, v="z", center=(0, 0, 0), rotate_cell=False): c * p - np.cross(p, s * v) + np.outer(np.dot(p, v), (1.0 - c) * v) + center ) if rotate_cell: - rotcell = self.get_cell() - rotcell[:] = ( - c * rotcell - - np.cross(rotcell, s * v) - + np.outer(np.dot(rotcell, v), (1.0 - c) * v) + self.cell[:] = ( + c * self.cell() + - np.cross(self.cell(), s * v) + + np.outer(np.dot(self.cell(), v), (1.0 - c) * v) ) - self.set_cell(rotcell) def _centering_as_array(self, center): if isinstance(center, str): @@ -1575,8 +1392,12 @@ def cell(self): @cell.setter def cell(self, cell): - cell = Cell.ascell(cell) - self._cellobj[:] = cell + if cell is None: + self._cellobj = None + else: + if (cell.shape[0] != 3) or (cell.shape[1] != 3): + raise Exception("The cell must be a 3x3 matrix.") + self._cellobj = cell def write(self, filename, format=None, **kwargs): """Write qbits object to a file. diff --git a/tests/qse/lattices_test.py b/tests/qse/lattices_test.py index a3547beb..b59f10e1 100644 --- a/tests/qse/lattices_test.py +++ b/tests/qse/lattices_test.py @@ -15,7 +15,7 @@ _repeats = [2, 3] -def _lattice_checker(qbits, expected_qbits, lattice_spacing, expected_cellpar): +def _lattice_checker(qbits, expected_qbits, lattice_spacing): assert isinstance(qbits, qse.Qbits) assert qbits.nqbits == expected_qbits @@ -23,8 +23,9 @@ def _lattice_checker(qbits, expected_qbits, lattice_spacing, expected_cellpar): # is the lattice spacing. assert np.isclose(qbits.get_all_distances()[0][1:].min(), lattice_spacing) - # Check the cellpar corresponds to what we expect. - assert np.allclose(qbits.cell.cellpar(), expected_cellpar) + # Check 2D + assert np.allclose(qbits.cell[2], np.zeros(3)) + assert np.allclose(qbits.positions[:, 2], np.zeros(qbits.nqbits)) # Check the labels assert all(qbits.labels[i] == str(i) for i in range(qbits.nqbits)) @@ -34,7 +35,7 @@ def _lattice_checker(qbits, expected_qbits, lattice_spacing, expected_cellpar): @pytest.mark.parametrize("N", _repeats) def test_linear(lattice_spacing, N): qbits = chain(lattice_spacing, N) - _lattice_checker(qbits, N, lattice_spacing, [lattice_spacing * N, 0, 0, 90, 90, 90]) + _lattice_checker(qbits, N, lattice_spacing) def test_linear_fail(): @@ -48,19 +49,15 @@ def test_linear_fail(): @pytest.mark.parametrize("N1", _repeats) @pytest.mark.parametrize("N2", _repeats) @pytest.mark.parametrize( - "lattice_func, angles", - [ - (square, [90, 90, 90]), - (triangular, [90, 90, 60]), - ], + "lattice_func", + [square, triangular], ) -def test_square_and_triangular(lattice_func, lattice_spacing, N1, N2, angles): +def test_square_and_triangular(lattice_func, lattice_spacing, N1, N2): qbits = lattice_func(lattice_spacing, N1, N2) _lattice_checker( qbits, N1 * N2, lattice_spacing, - [lattice_spacing * N1, lattice_spacing * N2, 0] + angles, ) @@ -73,8 +70,6 @@ def test_hexagonal(lattice_spacing, N1, N2): qbits, N1 * N2 * 2, lattice_spacing, - [lattice_spacing * N1 * np.sqrt(3), lattice_spacing * N2 * np.sqrt(3), 0] - + [90, 90, 60], ) @@ -87,7 +82,6 @@ def test_kagome(lattice_spacing, N1, N2): qbits, N1 * N2 * 3, lattice_spacing, - [lattice_spacing * N1 * 2, lattice_spacing * N2 * 2, 0] + [90, 90, 60], ) From add90ae4012b4a81bf9f44dca741c73e19f843c4 Mon Sep 17 00:00:00 2001 From: James Nelson Date: Wed, 27 Aug 2025 14:32:52 +0100 Subject: [PATCH 2/2] vis --- qse/qbits.py | 34 ++++++++-------------------------- qse/visualise.py | 18 +++++++++++------- 2 files changed, 19 insertions(+), 33 deletions(-) diff --git a/qse/qbits.py b/qse/qbits.py index e48c5047..291f0f36 100644 --- a/qse/qbits.py +++ b/qse/qbits.py @@ -37,10 +37,10 @@ class Qbits: scaled_positions: list of scaled-positions Like positions, but given in units of the unit cell. Can not be set at the same time as positions. - cell: 3x3 matrix, optional + cell: np.ndarray, optional The unit cell vectors. - First vector will lie in x-direction, second in xy-plane, - and the third one in z-positive subspace. + The first vector corresponds to the x-direction, second to the + xy-plane and the third to the z-positive subspace. celldisp: Vector Unit cell displacement vector. To visualize a displaced cell around the center of mass of a Systems of qbits. Default value @@ -159,7 +159,6 @@ def __init__( if scaled_positions is None: positions = np.zeros((len(self.arrays["labels"]), 3)) else: - assert self.cell.rank == 3 positions = np.dot(scaled_positions, self.cell) self.new_array("positions", positions, float, (3,)) @@ -408,8 +407,6 @@ def fromdict(cls, dct): constraints = [dict2constraint(d) for d in constraints] - # labels = dct.pop('labels', None) - info = dct.pop("info", None) qbits = cls( @@ -582,15 +579,14 @@ def pop(self, i=-1): def __imul__(self, m): """In-place repeat of qbits.""" - if isinstance(m, int): - m = (m, m, m) - if self.cell is None: raise ValueError("Cannot repeat without a defined unit cell.") for x, vec in zip(m, self.cell): if x != 1 and not vec.any(): raise ValueError("Cannot repeat without a defined unit cell.") + if isinstance(m, int): + m = (m, m, m) M = np.prod(m) n = len(self) @@ -599,13 +595,12 @@ def __imul__(self, m): self.arrays[name] = np.tile(a, (M,) + (1,) * (len(a.shape) - 1)) positions = self.arrays["positions"] - i0 = 0 + i = 0 for m0 in range(m[0]): for m1 in range(m[1]): for m2 in range(m[2]): - i1 = i0 + n - positions[i0:i1] += np.dot((m0, m1, m2), self.cell) - i0 = i1 + positions[i : i + n] += np.dot((m0, m1, m2), self.cell) + i += n if self.constraints is not None: self.constraints = [c.repeat(m, n) for c in self.constraints] @@ -1330,16 +1325,6 @@ def __ne__(self, other): else: return not eq - def get_volume(self): - """Get volume of unit cell.""" - if self.cell.rank != 3: - raise ValueError( - "You have {0} lattice vectors: volume not defined".format( - self.cell.rank - ) - ) - return self.cell.volume - def _get_positions(self): """Return reference to positions-array for in-place manipulations.""" return self.arrays["positions"] @@ -1365,14 +1350,12 @@ def _set_states(self, sts): : ] = sts # (sts.T / np.linalg.norm(sts, axis=1)).T # need to be normalized - # below is equivalent to defining @property states states = property( _get_states, _set_states, doc="Attribute for direct " + "manipulation of the states.", ) - # Rajarshi: Write method to get labels def _get_labels(self): """Return array of labels""" return self.arrays["labels"] @@ -1387,7 +1370,6 @@ def _set_labels(self, lbs): @property def cell(self): - """The :class:`ase.cell.Cell` for direct manipulation.""" return self._cellobj @cell.setter diff --git a/qse/visualise.py b/qse/visualise.py index cb5ab04d..f9800dc5 100644 --- a/qse/visualise.py +++ b/qse/visualise.py @@ -33,15 +33,19 @@ def draw(qbits, radius=None, show_labels=False, colouring=None, units=None): raise Exception("The length of colouring must equal the number of Qubits.") colouring = [int(i) for i in colouring] - cell_rank = qbits.cell.rank position_rank = np.linalg.matrix_rank(qbits.positions) - # if cell_rank is more than position_rank, it means that a higher dimensional cell - # is present, and actual structure requires repettition of cells to - # clearly visualize. if position_rank is more than cell_rank, it means the repettion - # happens in lower dimension of a local geometry which visualized in - # higher dimension. Either way, we need to see things in higher dimension. - rank = max(cell_rank, position_rank) + if qbits.cell is None: + rank = position_rank + else: + cell_rank = np.linalg.matrix_rank(qbits.cell) + rank = max(cell_rank, position_rank) + # if cell_rank is more than position_rank, it means that a higher + # dimensional cell is present, and actual structure requires + # repettition of cells to clearly visualize. + # If position_rank is more than cell_rank, it means the repettion + # happens in lower dimension of a local geometry which visualized in + # higher dimension. Either way, we need to see things in higher dimension. draw_bonds = False if radius is None else True rij = None