Разделение триклинной ячейки в Python
В настоящее время я анализирую данные молекулярной динамики и хотел бы разделить суперячейку, которую я использую, на части NxNxN. Форма ячейки триклинная, или, другими словами, это параллелепипед, у которого все три основные длины сторон и три главных угла различны. После этого поста здесь я смог понять, как построить изображение моей исходной клетки. Однако я изо всех сил пытаюсь правильно разделить исходную ячейку. На самом деле моя основная проблема заключается в том, что существуют небольшие перемещения подячеек, которые приводят к их смещению друг относительно друга, а также смещению относительно исходной ячейки.
Я почти все сделал правильно с помощью следующего кода:
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
[Изображение было повернуто вручную, чтобы лучше видеть подразделения]
Ниже приведен код, используемый для создания триклинной ячейки и ее подразделений (этот код должен быть более эффективным, чем код, который вы все еще можете найти в истории редактирования).
Я думаю, что это довольно просто, но если что-то неясно, спрашивайте в комментариях.
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()
В этом случае все было бы в порядке, если бы у вас уже были векторы ребер, верно?
@BabaBooey Вы определяете триклинную ячейку через начало координат O и 3 вектора V1, V2 и V3 или через 4 точки P0, P1, P2, P3. В первом случае у вас уже есть векторы ребер, во втором 3 вектора можно вычислить как V1=P1-P0 и т. д.
Ах хорошо. Мне нужно будет вычислить края вектора, чтобы это можно было правильно использовать, но это не должно быть слишком сложно. Теперь, когда я не занят, я протестирую это и посмотрю, как это работает. Я отчитаюсь в ближайшее время
Кажется, это работает очень хорошо, спасибо!
Здесь вы упускаете некоторую информацию. Вы не можете разделить ячейку, зная только начало координат и (я предполагаю; вы не сказали) три измерения и количество подячеек. Как a,b,c могут определить, правильная это ячейка или перекошенная? (и, конечно, насколько перекошено)? Думаю, вы это знаете, поскольку рисуете перекошенные подячейки. Но вы этого не предоставили