在Python中细分三斜晶格

Question

我目前正在分析分子动力学数据，并希望将我正在使用的超级单元放入 NxNxN 零件中。泡孔形状是三斜晶系，或者换句话说，它是一个平行六面体，其中 3 个主边长和 3 个主角都不同。在这篇文章here之后，我能够弄清楚如何构建我的原始细胞的图像。但是，我正在努力正确细分原始单元格。我的主要问题实际上是子单元格存在小的平移，导致它们彼此不对齐以及相对于原始单元格不对齐。

我用下面的代码几乎使一切都正确了：

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

我还将附上我在几个方向上获得的输出图，以便清楚地了解子单元是如何未对齐的。欢迎任何关于如何解决这个问题的其余部分的见解！

编辑：我正在使用 2x2x2 细分测试代码，但理想情况下我想要 4x4x4 细分。我预计一旦 2x2x2 工作，4x4x4 也应该能正常工作

方向1 方向2

Answer 1

[手动旋转图像以更好地查看细分]

下面是用于生成三斜晶格及其细分的代码（此代码应该比您仍然可以在编辑历史记录中找到的代码更有效）。

我认为这很简单，但如果有不清楚的地方请在评论中询问。

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)
        
        # draw an invisible line from the origin to the opposite vertex
        # to get a proper scaling of the plot
        ax.plot(*zip(o, o+a+b+c), lw=0)