Skip to content

minimum_image_displacement only valid for orthogonal cells #592

Description

@danielzuegner

The current implementation of minimum_image_displacement converts a displacement to fractional coordinates and rounds each periodic component independently:

dr_frac -= torch.where(pbc, torch.round(dr_frac), torch.zeros_like(dr_frac))

Component-wise rounding finds the minimum image for orthogonal cells, but not for a general triclinic cell. In a non-orthogonal lattice, minimizing each fractional component independently does not necessarily minimize the Cartesian norm.

Minimal example

import torch

from torch_sim.transforms import minimum_image_displacement

cell = torch.tensor(
    [
        [1.0, 0.9, 0.0],
        [0.0, 0.1, 0.0],
        [0.0, 0.0, 1.0],
    ],
    dtype=torch.float64,
)
frac = torch.tensor([0.49, 0.49, 0.0], dtype=torch.float64)
dr = cell @ frac

mic = minimum_image_displacement(
    dr=dr,
    cell=cell,
    pbc=torch.tensor([True, True, False]),
)

# The current implementation rounds frac to [0, 0, 0].
print(mic, torch.linalg.vector_norm(mic))
# tensor([0.9310, 0.0490, 0.0000]) tensor(0.9323)

# The image shifted by the second lattice vector is much closer.
closer_image = dr - cell[:, 1]
print(closer_image, torch.linalg.vector_norm(closer_image))
# tensor([ 0.0310, -0.0510, 0.0000]) tensor(0.0597)

The minimum-image problem for a general cell is instead the closest-vector problem

n* = argmin_{n in Z^3} ||dr - cell @ n||

over shifts along the periodic lattice directions.

Would it make sense to use lattice reduction followed by enumeration of nearby images, or another bounded closest-vector implementation? Checking only the 27 images adjacent in fractional coordinates is sufficient for suitable reduced cells, but is not guaranteed for an arbitrary unreduced, highly skewed basis. If support is intentionally limited to orthogonal cells, documenting and validating that restriction would also prevent silently incorrect distances.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions