Is there any Python library for returning the atoms of a periodic material within a cube?

Question

Is there any Python package that allows to cut defined by the points (0, 0, 0) and (x, y, z)? That is, it should return all the atoms of the unit cell and periodic images inside the shaded region.

(For sake of simplicity the image is for 2D).

An example

I have attached a real example from the ZnMOF-74.cif. I want the substructure defined by the red box. And a link for the cif file.

Answer 1

You can use a python library called ASE (Atomic Simulation Environment) to achieve your requirement.

Here is a sample code:

from ase import Atoms
from ase.build import make_supercell
import numpy as np
from pprint import pprint


def atoms_in_plane(unit_cell_atoms, plane_point1, plane_point2):
    # Create a supercell to ensure periodic boundary conditions
    supercell = make_supercell(unit_cell_atoms, [[2, 0, 0], [0, 2, 0], [0, 0, 2]])

    # Define the plane using the points (0, 0, 0) and (x, y, z)
    plane_normal = np.array(plane_point2) - np.array(plane_point1)


    # Collect atoms within the plane
    atoms_in_plane = []
    for atom in supercell:
        position_vector = atom.position
        # Calculate the vector from the first point of the plane to the atom
        vector_to_atom = position_vector - plane_point1
        # Calculate the dot product between the vector to the atom and the plane normal
        dot_product = plane_normal.dot(vector_to_atom)
        # If the dot product is close to zero, the atom lies within the plane
        if abs(dot_product) < 1e-8:
            atoms_in_plane.append(atom)

    return atoms_in_plane

# Example usage:
# Define your unit cell atoms
unit_cell_atoms = Atoms(symbols=['Zn', 'O', 'O', 'O', 'C', 'C', 'C', 'C', 'H'],
                        positions=[
                         [0.64159, 0.68514,0.47189], 
                         [0.96943,0.29776,0.36038], 
                         [0.92763,0.23048,0.6102], 
                         [0.25494,0.94205,0.67383], 
                         [0.92667,0.2458,0.42703], 
                         [0.87918,0.20711,0.29083], 
                         [0.21254,0.88926,0.7549], 
                         [0.1663,0.84759,0.63558] ,
                         [0.16529,0.85803,0.47918]
                         ])

# Define the two points that define the plane
plane_point1 = [0,0,0]
plane_point2 = [1e-8, 1e-9, 1e-8]


# Get the atoms within the plane
atoms_in_defined_plane = atoms_in_plane(unit_cell_atoms, plane_point1, plane_point2)

print("Atoms within the defined plane:")
pprint(atoms_in_defined_plane)

Output:

Atoms within the defined plane:
[Atom('C', [0.1663, 0.84759, 0.63558], index=7),
 Atom('H', [0.16529, 0.85803, 0.47918], index=8),
 Atom('C', [0.1663, 0.84759, 0.63558], index=16),
 Atom('H', [0.16529, 0.85803, 0.47918], index=17),
 Atom('C', [0.1663, 0.84759, 0.63558], index=25),
 Atom('H', [0.16529, 0.85803, 0.47918], index=26),
 Atom('C', [0.1663, 0.84759, 0.63558], index=34),
 Atom('H', [0.16529, 0.85803, 0.47918], index=35),
 Atom('C', [0.1663, 0.84759, 0.63558], index=43),
 Atom('H', [0.16529, 0.85803, 0.47918], index=44),
 Atom('C', [0.1663, 0.84759, 0.63558], index=52),
 Atom('H', [0.16529, 0.85803, 0.47918], index=53),
 Atom('C', [0.1663, 0.84759, 0.63558], index=61),
 Atom('H', [0.16529, 0.85803, 0.47918], index=62),
 Atom('C', [0.1663, 0.84759, 0.63558], index=70),
 Atom('H', [0.16529, 0.85803, 0.47918], index=71)]

Answer 2

You can also use sisl to achieve this.

Best shown with an example:

from sisl import geom, Cuboid

graphene = geom.graphene()
gr2x2 = graphene.tile(2, 0).tile(2, 1) # double along along 1st and 2nd lattice vectors

# a cube of side-lengths 4 Ang, centerred around 0.5, 0.5, 0.5 Ang
cube = Cuboid(4, center=[0.5, 0.5, 0.5])

idx, xyz = gr2x2.within(cube, ret_xyz=True)
print(xyz)

Note here that the geometry gr2x2 may have an associated auxiliary supercell, in which case atoms in that supercell are also considered. If this is not desired one can easily change that behaviour.
sisl allows many more shapes, and also combinations of shapes. For instance:

cubeandcube = cube + Cuboid(2, center=[2, 2, 2])
gr2x2.within(cubeandcube)

Disclaimer: I am the author of sisl.