Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add class method for Box class which create a box from Lx,Ly,Lz and angles #1234

Merged
merged 42 commits into from
Apr 25, 2024
Merged
Show file tree
Hide file tree
Changes from 19 commits
Commits
Show all changes
42 commits
Select commit Hold shift + click to select a range
0d6f975
add from lattice_vectors, and to and from lenghts and angles
DomFijan Feb 27, 2024
0717da1
Add methods to convert box dimensions and angles
DomFijan Feb 27, 2024
ae95e63
Merge branch 'main' into feat_from_to_lattice_vectors
DomFijan Mar 7, 2024
5a388bf
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 7, 2024
58a56f6
add type hints and a class method from box lengths and angles for uni…
DomFijan Mar 20, 2024
746b923
add type hints
DomFijan Mar 20, 2024
2f9adcf
Merge branch 'feat_from_to_lattice_vectors' of https://github.com/glo…
DomFijan Mar 20, 2024
c7ad750
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 20, 2024
c399ff3
fix format
DomFijan Mar 20, 2024
2a99532
more formatting fixes
DomFijan Mar 20, 2024
16d028a
fixes to logic and docs to all new functions
DomFijan Mar 20, 2024
dc74488
add tests
DomFijan Mar 20, 2024
7706c2f
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 20, 2024
8860961
update some docstrings for consistency
DomFijan Mar 20, 2024
12b4a1a
Merge branch 'feat_from_to_lattice_vectors' of https://github.com/glo…
DomFijan Mar 20, 2024
0e78e47
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 20, 2024
da9d7ea
linter fix
DomFijan Mar 20, 2024
aa4354a
added changelog and credits
DomFijan Mar 20, 2024
d536985
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 20, 2024
461485e
Test dimension detection for `from_lattice_vectors`
janbridley Mar 20, 2024
90dd30c
Update freud/box.pyx
DomFijan Mar 21, 2024
8e27b3c
apply requested changes
DomFijan Mar 21, 2024
863ac23
Merge branch 'feat_from_to_lattice_vectors' of https://github.com/glo…
DomFijan Mar 21, 2024
5bd2c08
revert unwanted changes and delete tests
DomFijan Mar 21, 2024
6a9b740
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 21, 2024
1a45db1
fix format
DomFijan Mar 21, 2024
ace7fac
final format fix
DomFijan Mar 21, 2024
6f6932e
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 21, 2024
ae96836
fix from_box_lengths_and_angles
DomFijan Mar 21, 2024
3aab7cb
fix changelog and credits
DomFijan Mar 21, 2024
a7c18e2
Update freud/box.pyx
tommy-waltmann Mar 21, 2024
c26e88e
Apply suggestions from code review
tommy-waltmann Mar 21, 2024
ae56eff
Update ChangeLog.md
DomFijan Mar 22, 2024
a3daf1b
update docs to state range of angles allowed
DomFijan Mar 22, 2024
dffb772
update test
tommy-waltmann Mar 22, 2024
a3548a1
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 22, 2024
a3449d8
Set nonzero atol for assert_allclose in test_box_Box.py
janbridley Mar 26, 2024
b32ed14
Merge branch 'main' into feat_from_to_lattice_vectors
janbridley Apr 11, 2024
ca710c3
update from_box_lengths_and_angles docstring to contain restrictions …
DomFijan Apr 11, 2024
ad8d8b3
raise error if box has negative volume
DomFijan Apr 11, 2024
b51d8d7
fix sqrt
DomFijan Apr 11, 2024
d78e08e
fix test to check for imposible boxes
DomFijan Apr 23, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions ChangeLog.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,10 @@ and this project adheres to

### Added
* New continuous coordination number compute `freud.order.ContinuousCoordination`.
* New class methods for construction of `freud.box.Box` and `freud.data.UnitCell` from
lattice vectors and box lengths and angles.
* New `freud.box.Box.to_box_lengths_and_angles` method for computing box vector lengths
and angles

### Removed
* `freud.order.Translational`.
Expand Down
7 changes: 7 additions & 0 deletions doc/source/reference/credits.rst
Original file line number Diff line number Diff line change
Expand Up @@ -369,6 +369,13 @@ Domagoj Fijan
* Contributed code, design, documentation, and testing for ``freud.locality.FilterSANN`` class.
* Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
* Added support for ``gsd.hoomd.Frame`` in ``NeighborQuery.from_system`` calls.
* Added support for ``freud.box.Box`` class methods for construction of boxes and unit
cells from lattice vectors (``freud.box.Box.from_lattice_vectors`` and
``freud.data.UnitCell.from_lattice_vectors``) and cell lengths and angles
(``freud.box.Box.from_box_lengths_and_angles`` and
``freud.data.UnitCell.from_box_lengths_and_angles``), as well as a method for
returning cell vector lengths and angles
(``freud.box.Box.to_box_lengths_and_angles``).

Andrew Kerr

Expand Down
88 changes: 88 additions & 0 deletions freud/box.pyx
tommy-waltmann marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
Expand Up @@ -689,6 +689,25 @@ cdef class Box:
[0, self.Ly, self.yz * self.Lz],
[0, 0, self.Lz]])

def to_box_lengths_and_angles(self):
r"""Return the box lengths and angles.

Returns:
tuple: The box vector lengths and angles in radians
DomFijan marked this conversation as resolved.
Show resolved Hide resolved
"""
alpha = np.arccos(
(self.xy * self.xz + self.yz)
/ (np.sqrt(1 + self.xy**2) * np.sqrt(1 + self.xz**2 + self.yz**2))
)
beta = np.arccos(self.xz/np.sqrt(1+self.xz**2+self.yz**2))
gamma = np.arccos(self.xy/np.sqrt(1+self.xy**2))
L1 = self.Lx
a2 = [self.Ly*self.xy, self.Ly, 0]
a3 = [self.Lz*self.xz, self.Lz*self.yz, self.Lz]
L2 = np.linalg.norm(a2)
L3 = np.linalg.norm(a3)
return (L1, L2, L3, alpha, beta, gamma)

def __repr__(self):
return ("freud.box.{cls}(Lx={Lx}, Ly={Ly}, Lz={Lz}, "
"xy={xy}, xz={xz}, yz={yz}, "
Expand Down Expand Up @@ -921,6 +940,75 @@ cdef class Box:
"positional argument: L")
return cls(Lx=L, Ly=L, Lz=0, xy=0, xz=0, yz=0, is2D=True)

@classmethod
def from_lattice_vectors(cls, lattice_vectors, dimensions: int = None):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How is this different than from_matrix? Box.from_matrix makes a box from a matrix where the columns are the lattice vectors. This seems very similar to that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure how this got past me but with this in mind this opens a whole new debate. from_matrix in freud basically uses lattice vectors as input, but this is not stated explicitly. Instead a link to hoomd docs is given which further complicates the problem. In hoomd from_matrix wants upper triangular matrix that is already a box like object as defined in hoomd docs. What we should do is to have freuds from_matrix take upper triangular matrix as input, such that its functionality is equivalent to hoomd's and add from_lattice_vectors which has same function as current from_matrix?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be best to add to freud's documentation and explicitly state that the columns of the input to from_matrix are the lattice vectors of the box. Any other improvements to the documentation of from_matrix to make it less confusing would also be welcome.

Looking at the diff of the code between from_lattice_vectors and from_matrix, they do the exact same thing so I don't see a need to have both factory methods hanging around.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the functionality of freud's from_matrix equivalent to hoomd's from_matrix?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It appears that hoomd's from_matrix requires an upper triangular matrix, whereas freud's is more general and does not.

"""Create a unit cell from lattice vectors.

Args:
lattice_vectors (array-like):
The lattice vectors. The dimensions of the object should be 3x3. Lattice
vector a1 is lattice_vectors[0], etc.
dimensions (int): The number of dimensions (Default value = :code:`None`)

Returns:
:class:`~.UnitCell`: A unit cell with the given lattice vectors.
"""
lattice_matrix = np.asarray(lattice_vectors, dtype=np.float32)
v0 = lattice_matrix[0]
v1 = lattice_matrix[1]
v2 = lattice_matrix[2]
Lx = np.sqrt(np.dot(v0, v0))
a2x = np.dot(v0, v1) / Lx
Ly = np.sqrt(np.dot(v1, v1) - a2x * a2x)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Ly = np.sqrt(np.dot(v1, v1) - a2x * a2x)
Ly = np.sqrt(np.dot(v1, v1) - a2x**2)

I think a square here better conveys the intent, as below we have a very similar function with two different variables (a2x * a3x)

xy = a2x / Ly
v0xv1 = np.cross(v0, v1)
v0xv1mag = np.sqrt(np.dot(v0xv1, v0xv1))
Lz = np.dot(v2, v0xv1) / v0xv1mag
if Lz != 0:
a3x = np.dot(v0, v2) / Lx
xz = a3x / Lz
yz = (np.dot(v1, v2) - a2x * a3x) / (Ly * Lz)
else:
xz = yz = 0
if dimensions is None:
dimensions = 2 if Lz == 0 else 3
return cls.from_box([Lx, Ly, Lz, xy, xz, yz], dimensions=dimensions)

@classmethod
def from_box_lengths_and_angles(
cls,
L1: float,
L2: float,
L3: float,
alpha: float,
beta: float,
gamma: float,
dimensions: int = None,
DomFijan marked this conversation as resolved.
Show resolved Hide resolved
):
r"""Construct a box from lengths and angles.

Args:
L1 (float): The length of the first lattice vector
L2 (float): The length of the second lattice vector
L3 (float): The length of the third lattice vector
alpha (float): The angle between the y and z axes in radians
beta (float): The angle between the x and z axes in radians
gamma (float): The angle between the x and y axes in radians
DomFijan marked this conversation as resolved.
Show resolved Hide resolved
dimensions (int): The number of dimensions (Default value = :code:`None`)

Returns:
:class:`freud.box.Box`: The resulting box object.
"""
a1 = np.array([L1, 0, 0])
a2 = np.array([L2 * np.cos(gamma), L2 * np.sin(gamma), 0])
a3x = np.cos(beta)
a3y = (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) / np.sin(gamma)
a3z = np.sqrt(1 - a3x**2 - a3y**2)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This line appears to take the square root of a negative number at times, resulting in NaNs.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some combinations of angles are not possible-there is likely a constraint on their sum or something like that. Is this taken into account?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not in how I wrote the test case, which would explain the failures. The code should throw an error if the user gives inputs which violate the constraints.

a3 = np.array([L3 * a3x, L3 * a3y, L3 * a3z])
if dimensions is None:
dimensions = 2 if L3 == 0 else 3
return cls.from_lattice_vectors(np.array([a1, a2, a3]), dimensions=dimensions)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
return cls.from_lattice_vectors(np.array([a1, a2, a3]), dimensions=dimensions)
return cls.from_matrix(np.array([a1, a2, a3]).T, dimensions=dimensions)

from_lattice_vectors was removed



cdef BoxFromCPP(const freud._box.Box & cppbox):
b = Box(cppbox.getLx(), cppbox.getLy(), cppbox.getLz(),
Expand Down
50 changes: 50 additions & 0 deletions freud/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -233,6 +233,56 @@ def hex(cls):
fractions = np.array([[0, 0, 0], [0.5, 0.5, 0]])
return cls([1, np.sqrt(3)], fractions)

@classmethod
def from_lattice_vectors(
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since this method (and from_box_lengths_and_angles) just replaces two lines of user script, I don't think it's needed. I am open to discussion though.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's just a convenience function. I love these. If you strongly feel they are not needed we can remove them. Up to you.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My personal logic is that if we have a convenience function for making UnitCell's in the data module for Box.from_lattice_vectors and Box.from_box_lengths_and_angles, then we should also have a convenience function for every factory method in the Box class (like cube, from_box, from_matrix, square, etc.). There are quite a few of those and it may get cumbersome to have so many convenience functions hanging around.

If you like having that convenience function, add one to your python scripts. It only takes two lines to write it, after all :)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While I see what you're saying, there are a significant number of users/use cases for this specific functionality: making the translation between freud boxes and data from crystallography files more transparent is a good thing in my opinion. I think adding this particular convenience function provides enough to users that we can probably do without the rest mentioned above. If this truly doesn't fit though, I can live without it.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What data does a crystallography file typically contain? If its a common use case to make a freud UnitCell from crystallography file data, then it makes sense to have a convenience function for it. I personally don't work a lot with those, so I'm not entirely sure what would be most helpful for others who work with it.

I think this topic is worth further offline discussion as well.

Copy link
Contributor

@janbridley janbridley Mar 21, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CIF files (the most common file type, at least for computational/experimental researchers) have unit cell information encoded as a/b/c/α/β/γ. For DFT peeps, VASP and ABINIT provide lattice vectors so those would work without this function.

Crystallography Open Database

_cell_angle_alpha                90.000
_cell_angle_beta                 102.104(8)
_cell_angle_gamma                90.000
_cell_formula_units_Z            4
_cell_length_a                   7.9661(7)
_cell_length_b                   9.205(1)
_cell_length_c                   7.3198(8)
_cell_volume                     524.8(9)

AFlow:

_cell_length_a  5.4030000000
_cell_length_b  3.4330000000
_cell_length_c  5.0790000000
_cell_angle_alpha  90.0000000000
_cell_angle_beta  132.3200000000
_cell_angle_gamma  90.0000000000

cls, lattice_vectors: np.ndarray, unique_positions: np.ndarray
):
"""Create a unit cell from lattice vectors.

Args:
lattice_vectors (:math:`(3, 3)` :class:`np.ndarray
The lattice vectors. Lattice vector a1 is lattice_vectors[0], etc.
unique_positions (:math:`(N_{points}, 3)` :class:`np.ndarray`):
The basis positions.

Returns:
:class:`~.UnitCell`: A unit cell with the given lattice vectors.
"""
return cls(
freud.box.Box.from_lattice_vectors(lattice_vectors), unique_positions
)

@classmethod
def from_box_lengths_and_angles(
cls,
L1: float,
L2: float,
L3: float,
alpha: float,
beta: float,
gamma: float,
unique_positions: np.ndarray,
):
"""Create a unit cell from box lengths and angles.

Args:
L1 (float): The length of the first box vector.
L2 (float): The length of the second box vector.
L3 (float): The length of the third box vector.
alpha (float): The angle between the x and y lattice vectors.
beta (float): The angle between the x and z lattice vectors.
gamma (float): The angle between the y and z lattice vectors.
unique_positions (:math:`(N_{points}, 3)` :class:`np.ndarray
The basis positions.

Returns:
:class:`~.UnitCell`: A unit cell with the given box lengths and angles.
"""
return cls(
freud.box.Box.from_box_lengths_and_angles(L1, L2, L3, alpha, beta, gamma),
unique_positions,
)


def make_random_system(box_size, num_points, is2D=False, seed=None):
r"""Helper function to make random points with a cubic or square box.
Expand Down
57 changes: 57 additions & 0 deletions tests/test_box_Box.py
Original file line number Diff line number Diff line change
Expand Up @@ -485,6 +485,63 @@ def test_from_box(self):
box7 = freud.box.Box.from_matrix(box.to_matrix())
assert np.isclose(box.to_matrix(), box7.to_matrix()).all()

def test_standard_orthogonal_box(self):
box = freud.box.Box.from_box((1, 2, 3, 0, 0, 0))
Lx, Ly, Lz, alpha, beta, gamma = box.to_box_lengths_and_angles()
npt.assert_allclose(
(Lx, Ly, Lz, alpha, beta, gamma), (1, 2, 3, np.pi / 2, np.pi / 2, np.pi / 2)
)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a great test to add would be to start with 3 random lengths and 3 random angles, then call from_box_lengths_and_angles to get a box, then call to_box_lengths_and_angles on that box and assert that you recover the original random lengths and angles.

def test_to_and_from_box_lengths_and_angles(self):
original_box_lengths_and_angles = (
1,
2,
3,
np.pi / 3,
np.pi / 2.25,
np.pi / 2.35,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Each of the six numbers should be the output of a different call to np.random.rand. That way, each time the test is run it gets run with a different random value.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am currently working on adding this.

)
box = freud.box.Box.from_box_lengths_and_angles(
*original_box_lengths_and_angles
)
lengths_and_angles_computed = box.to_box_lengths_and_angles()
assert np.allclose(lengths_and_angles_computed, original_box_lengths_and_angles)

def test_from_lattice_vectors_cubic_lattice(self):
DomFijan marked this conversation as resolved.
Show resolved Hide resolved
lattice_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
unit_cell = freud.box.Box.from_lattice_vectors(lattice_vectors)
assert unit_cell.Lx == 1
assert unit_cell.Ly == 1
assert unit_cell.Lz == 1
assert unit_cell.xy == 0
assert unit_cell.xz == 0
assert unit_cell.yz == 0

def test_from_lattice_vectors_non_cubic_lattice(self):
box = freud.box.Box.from_lattice_vectors([[1, 0, 0], [0, 2, 0], [0, 0, 3]])
assert box.Lx == 1
assert box.Ly == 2
assert box.Lz == 3
assert box.xy == 0
assert box.xz == 0
assert box.yz == 0

def test_tilted_lattice(self):
lattice_vectors = np.array([[1, 0, 0], [0.5, np.sqrt(3) / 2, 0], [0, 0, 1]])
unit_cell = freud.box.Box.from_lattice_vectors(lattice_vectors)
assert np.isclose(unit_cell.Lx, 1.0)
assert np.isclose(unit_cell.Ly, np.sqrt(3) / 2)
assert np.isclose(unit_cell.Lz, 1.0)
assert np.isclose(unit_cell.xy, np.sqrt(3) / 3)
assert np.isclose(unit_cell.xz, 0.0)
assert np.isclose(unit_cell.yz, 0.0)

def test_dimensions(self):
lattice_vectors = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
unit_cell = freud.box.Box.from_lattice_vectors(lattice_vectors, dimensions=2)
assert unit_cell.dimensions == 2
assert np.isclose(unit_cell.Lz, 0)

def test_matrix(self):
box = freud.box.Box(2, 2, 2, 1, 0.5, 0.1)
box2 = freud.box.Box.from_matrix(box.to_matrix())
Expand Down