diff --git a/pyproject.toml b/pyproject.toml index dfb75c4..68d4200 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,6 @@ authors = [ version = "1.0.3" requires-python = ">=3.10" dependencies = [ - "ase", "numpy", "matplotlib" ] diff --git a/qse/__init__.py b/qse/__init__.py index 00c108a..0a54582 100644 --- a/qse/__init__.py +++ b/qse/__init__.py @@ -7,7 +7,7 @@ __all__ = [ "calc", - "cell", + "Cell", "draw", "lattices", "magnetic", @@ -21,6 +21,7 @@ __version__ = importlib.metadata.version("qse") +from qse.cell import Cell from qse.qbit import Qbit from qse.qbits import Qbits from qse.signal import Signal diff --git a/qse/cell.py b/qse/cell.py new file mode 100644 index 0000000..85573c6 --- /dev/null +++ b/qse/cell.py @@ -0,0 +1,74 @@ +import numpy as np + + +class Cell: + """ + A class to represent a lattice cell. + + Attributes + ---------- + lattice_vectors : numpy.ndarray + A 3x3 matrix representing the lattice vectors. + Each row is a lattice vector. + Defaults to all zeros. + """ + + def __init__(self, lattice_vectors=None): + if lattice_vectors is None: + lattice_vectors = np.zeros((3, 3)) + self.lattice_vectors = lattice_vectors + + def __repr__(self): + return self.to_str() + + @property + def lattice_vectors(self): + return self._lattice_vectors + + @lattice_vectors.setter + def lattice_vectors(self, value): + value = np.array(value) + assert value.shape == (3, 3) + self._lattice_vectors = value + + def rank(self): + """ + Compute the rank of the cell. + + Returns + ------- + int + The rank of the cell. + """ + return np.linalg.matrix_rank(self.lattice_vectors) + + def volume(self): + """ + Compute the volume of the cell. + + Returns + ------- + float + The cell volume. + """ + return np.linalg.det(self.lattice_vectors) + + def reciprocal(self): + """ + Compute the reciprocal lattice vectors. + + Returns + ------- + numpy.ndarray + The reciprocal cell, calculated as 2π * inverse of the + transposed cell matrix. + """ + return 2 * np.pi * np.linalg.inv(self.lattice_vectors.T) + + def to_str(self): + return "\n".join( + [ + " ".join(["{:10.5f}".format(val) for val in row]) + for row in self.lattice_vectors + ] + ) diff --git a/qse/qbits.py b/qse/qbits.py index a6eb2bf..35b83c0 100644 --- a/qse/qbits.py +++ b/qse/qbits.py @@ -6,8 +6,8 @@ from math import cos, sin import numpy as np -from ase.cell import Cell +from qse.cell import Cell from qse.qbit import Qbit from qse.visualise import draw as _draw @@ -33,15 +33,11 @@ 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)]. - First vector will lie in x-direction, second in xy-plane, - and the third one in z-positive subspace. - Default value: [0, 0, 0]. + cell: qse.Cell | 3x3 matrix + Unit cell vectors. + Either pass a matrix where each row corresponds to a lattice vector + or a qse.Cell object. + Default value is all zeros. pbc: one or three bool Periodic boundary conditions flags. Examples: True, False, 0, 1, (1, 1, 0), (True, False, False). Default @@ -68,8 +64,7 @@ class Qbits: >>> xd = np.array( ... [[0, 0, 0], ... [0.5, 0.5, 0.5]]) - >>> qdim = qse.Qbits(positions=xd) - >>> qdim.cell = [1,1,1] + >>> qdim = qse.Qbits(positions=xd, cell=np.eye(3)) >>> qdim.pbc = True >>> qlat = qdim.repeat([3,3,3]) @@ -130,18 +125,17 @@ 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) - - if scale_qbits: - M = np.linalg.solve(self.cell.complete(), cell.complete()) - self.positions[:] = np.dot(self.positions, M) - - self.cell[:] = cell - - 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() - - def get_reciprocal_cell(self): - """ - Get the three reciprocal lattice vectors as a 3x3 ndarray. - - Returns - ------- - np.ndarray - The reciprocal cell. - - Notes - ----- - 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) @@ -372,12 +289,11 @@ def copy(self): def todict(self): """For basic JSON (non-database) support.""" - # d = dict(self.arrays) d = {} 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 return d @@ -425,13 +341,8 @@ def __repr__(self): txt = "pbc={0}".format(self.pbc[0]) 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") + tokens.append("cell=\n{0}".format(self.cell.to_str()) + ",\n") + if self._calc is not None: tokens.append("calculator={0}".format(self._calc.__class__.__name__)) txt = "{0}(\n{1}{2})".format( @@ -560,7 +471,7 @@ def __imul__(self, m): if isinstance(m, int): m = (m, m, m) - for x, vec in zip(m, self.cell): + for x, vec in zip(m, self.cell.lattice_vectors): if x != 1 and not vec.any(): raise ValueError("Cannot repeat along undefined lattice " "vector") @@ -576,10 +487,12 @@ def __imul__(self, m): 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) + positions[i0:i1] += np.dot((m0, m1, m2), self.cell.lattice_vectors) i0 = i1 - self.cell = np.array([m[c] * self.cell[c] for c in range(3)]) + self.cell.lattice_vectors = np.array( + [m[c] * self.cell.lattice_vectors[c] for c in range(3)] + ) new_shape = tuple(self.shape * np.array(m)) self.shape = new_shape @@ -783,20 +696,18 @@ 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.lattice_vectors = ( + c * self.cell.lattice_vectors + - np.cross(self.cell.lattice_vectors, s * v) + + np.outer(np.dot(self.cell.lattice_vectors, v), (1.0 - c) * v) ) - self.set_cell(rotcell) def _centering_as_array(self, center): if isinstance(center, str): if center.lower() == "cop": center = self.get_centroid() elif center.lower() == "cou": - center = self.get_cell().sum(axis=0) / 2 + center = self.cell.lattice_vectors.sum(axis=0) / 2 else: raise ValueError("Cannot interpret center") else: @@ -1188,16 +1099,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"] @@ -1243,16 +1144,6 @@ def _set_labels(self, lbs): _get_labels, _set_labels, doc="Attribute for direct manipulation of labels" ) - @property - def cell(self): - """The :class:`ase.cell.Cell` for direct manipulation.""" - return self._cellobj - - @cell.setter - def cell(self, cell): - cell = Cell.ascell(cell) - self._cellobj[:] = cell - def write(self, filename, format=None, **kwargs): """Write qbits object to a file. diff --git a/tests/qse/cell_test.py b/tests/qse/cell_test.py new file mode 100644 index 0000000..6901cd2 --- /dev/null +++ b/tests/qse/cell_test.py @@ -0,0 +1,76 @@ +import numpy as np +import pytest + +import qse + + +@pytest.mark.parametrize( + "cell, expected", + [ + (4.232 * np.eye(3), np.eye(3) / 4.232), # square + ( + np.array( # rectangular + [ + [1.0, 0.0, 0.0], + [0.0, 2.0, 0.0], + [0.0, 0.0, 1.0], + ] + ), + np.array( + [ + [1.0, 0.0, 0.0], + [0.0, 1.0 / 2.0, 0.0], + [0.0, 0.0, 1.0], + ] + ), + ), + ( + np.array( # hex + [ + [1.0, 0.0, 0.0], + [0.5, 0.5 * np.sqrt(3), 0.0], + [0.0, 0.0, 1.0], + ] + ), + np.array( + [ + [1.0, -1.0 / np.sqrt(3), 0.0], + [0.0, 2.0 / np.sqrt(3), 0.0], + [0.0, 0.0, 1.0], + ] + ), + ), + ( + 0.5 + * np.array( # fcc + [ + [0.0, 1.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 0.0], + ] + ), + np.array( + [ + [-1.0, 1.0, 1.0], + [1.0, -1.0, 1.0], + [1.0, 1.0, -1.0], + ] + ), + ), + ], +) +def test_reciprocal_cell(cell, expected): + """Test the reciprocal method.""" + cell = qse.Cell(cell) + r_cell = cell.reciprocal() / (2 * np.pi) + assert np.allclose(r_cell, expected) + + # Check the two cells are orthonormal + assert np.allclose(cell.lattice_vectors @ r_cell.T, np.eye(3)) + + +@pytest.mark.parametrize("lattice_size", [1.0, 2.2, 3.11, 10.333]) +def test_volume(lattice_size): + """Test the volume method.""" + cell = qse.Cell(np.eye(3) * lattice_size) + assert np.isclose(cell.volume(), lattice_size**3) diff --git a/tests/qse/lattices_test.py b/tests/qse/lattices_test.py index a9e421b..e93e77a 100644 --- a/tests/qse/lattices_test.py +++ b/tests/qse/lattices_test.py @@ -30,7 +30,7 @@ def _lattice_checker( assert np.isclose(qbits.get_all_distances()[0][1:].min(), lattice_spacing) # Check the cell - assert np.allclose(qbits.cell, expected_cell) + assert np.allclose(qbits.cell.lattice_vectors, expected_cell) # Check the positions assert np.allclose(qbits.positions, expected_positions) diff --git a/tests/qse/qbits_methods_test.py b/tests/qse/qbits_methods_test.py index 3163b91..13c7e47 100644 --- a/tests/qse/qbits_methods_test.py +++ b/tests/qse/qbits_methods_test.py @@ -404,68 +404,3 @@ def test_repeat(nqbits, repeats): ps = np.array(ps) assert np.allclose(ps, qbits.positions) - - -@pytest.mark.parametrize( - "cell, expected", - [ - (4.232 * np.eye(3), np.eye(3) / 4.232), # square - ( - np.array( # rectangular - [ - [1.0, 0.0, 0.0], - [0.0, 2.0, 0.0], - [0.0, 0.0, 1.0], - ] - ), - np.array( - [ - [1.0, 0.0, 0.0], - [0.0, 1.0 / 2.0, 0.0], - [0.0, 0.0, 1.0], - ] - ), - ), - ( - np.array( # hex - [ - [1.0, 0.0, 0.0], - [0.5, 0.5 * np.sqrt(3), 0.0], - [0.0, 0.0, 1.0], - ] - ), - np.array( - [ - [1.0, -1.0 / np.sqrt(3), 0.0], - [0.0, 2.0 / np.sqrt(3), 0.0], - [0.0, 0.0, 1.0], - ] - ), - ), - ( - 0.5 - * np.array( # fcc - [ - [0.0, 1.0, 1.0], - [1.0, 0.0, 1.0], - [1.0, 1.0, 0.0], - ] - ), - np.array( - [ - [-1.0, 1.0, 1.0], - [1.0, -1.0, 1.0], - [1.0, 1.0, -1.0], - ] - ), - ), - ], -) -def test_reciprocal_cell(cell, expected): - """Test the get_reciprocal_cell method.""" - qbits = qse.Qbits(cell=cell) - r_cell = qbits.get_reciprocal_cell() - assert np.allclose(r_cell, expected) - - # Check the two cells are orthonormal - assert np.allclose(cell @ r_cell.T, np.eye(3))