Subdividing triclinic cell in Python

Question

I am currently analyzing molecular dynamics data, and would like to bin the supercell I am using in NxNxN parts. The cell shape is triclinic, or in other words it is a parallelepiped where the 3 main side lengths and 3 main angles are all different. Following this post here, I was able to figure out how to construct an image of my original cell. However, I am struggling to properly subdivide the original cell. My main issue actually is that there are small translations of the subcells that cause them to be misaligned from each other as well as misaligned relative to the original cell.

I have almost made everything correct with this following code:

def subdivide_cell(origin, a, b, c, n):
    subdivisions = []
    for i in range(n):
        for j in range(n):
            for k in range(n):
                sub_origin = origin + np.array([i * a / n, j * b / n, k * c / n])
                sub_a = a / n
                sub_b = b / n
                sub_c = c / n
                subdivisions.append((sub_origin, sub_a, sub_b, sub_c))
    return subdivisions

I also will attach the output figure I get in a couple orientations so it is clear how the subcells are misaligned. Any insight on how to figure the rest of this problem out is welcome!

Edit: I am testing the code with a 2x2x2 subdivision, but ideally I would like 4x4x4 subdivision. I expect that once the 2x2x2 works, the 4x4x4 should be pretty functional as well

Orientation 1 Orientation 2

Answer 1

[The image was rotated manually to have a good view of the subdivisions]

Below the code used to produce the triclinic cell and its subdivisions (this code should be more efficient that the code you can still find in the edit history).

I think it's quite straightforward, but please ask in comments if something is not clear.

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Line3DCollection as lc3d
import numpy as np
from itertools import product

# a solid parallelogram is defined in terms of an origin and 3 vectors,
# when no couple of vectors are orthogonal you have a tricline
o, a, b, c = np.array((
    ( 0.0,  0.0,  0.0),
    ( 6.0,  1.0, -1.0),
    ( 0.0,  7.0,  1.0),
    ( 2.0,  0.0,  7.0),
    ))

# let's get a 3d Axes
ax = plt.figure().add_subplot(projection='3d')

# first the inner cells, then the outer, single one, darker and thicker
draw_n_cells(o, a, b, c, 3, lw=0.5, color='gray')
draw_n_cells(o, a, b, c, 1, lw=1.5, color='k')
plt.show()

def draw_n_cells(o, a, b, c, n, ax=None, lw=1.5, color='blue'):
    ax = ax if ax else plt.gca()
    # the idea is to use two vectors to form a grid of points on one of the
    # six faces of the parallelogram, and draw the lines // to the third
    # vector that arrive to the opposite face, one line per point
    
    # helper vector, equispaced point on the interval [0, 1] spaced 1/n
    d = np.linspace(0, 1, n+1)[:,None] ####> d.shape is (n+1, 1)

    # we make a cycle to have the 2 vectors, r & s, that determine a face and
    # the vector t that determines the line direction and its extension
    for r, s, t in ((a, b, c), (b, c, a), (c, a, b)):
        # arrays containing equispaced 3d points on r and s respectively
        rn, sn = r*d, s*d
        # each point on the grid is defined by two vectors, say v0 and v1, that
        # must be summed to find the position of a point
        # 1. we compute the couples using itertools.product
        v0v1 =  product(rn, sn)
        # 2. we compute the points on the grid, spanning one face, and also the
        # points on the opposite face
        p0_p1 = [[o+p, o+p+t]for p in (v0+v1 for v0, v1 in v0v1)]
        # the list [[p0, p1]] is exactly in the format that's needed for LineCollection3D
        lc = lc3d(p0_p1, lw=lw, color=color)
        # or, more compact
        # lc = lc3d([[o+p, o+p+t]for p in (v0+v1 for v0, v1 in product(r*d, s*d))])
        # the line collection must be placed in the 3D Axes
        ax.add_artist(lc)
        ax.autoscale()