Диэдра / угол кручения из четырех точек в декартовых координатах в Python

Какие предложения имеют люди для быстрого вычисления двугранных углов в Python?

На диаграммах phi – двугранный угол:

Двугранный угол1

Двугранный угол2

Что лучше всего для расчета углов в диапазоне от 0 до pi? Как насчет от 0 до 2pi?

«Лучший» здесь означает некоторое сочетание быстрого и численно стабильного. Методы, возвращающие значения во всем диапазоне от 0 до 2pi, являются предпочтительными, но если у вас есть невероятно быстрый способ вычисления диэдра в пределах от 0 до pi, это тоже.

Вот мои 3 лучших усилия. Только второй возвращает углы между 0 и 2pi. Это также самый медленный.

Общие комментарии о моих подходах:

arccos () в Numpy кажется достаточно стабильным, но поскольку люди поднимают эту проблему, я могу просто не полностью ее понять.

Из этого вышло использование einsum. Почему функция numpy einsum работает быстрее, чем встроенные функции numpy?

Отсюда и диаграммы, и вдохновение. Как рассчитать двугранный угол, заданный декартовыми координатами?

3 подхода с комментариями:

import numpy as np from time import time # This approach tries to minimize magnitude and sqrt calculations def dihedral1(p): # Calculate vectors between points, b1, b2, and b3 in the diagram b = p[:-1] - p[1:] # "Flip" the first vector so that eclipsing vectors have dihedral=0 b[0] *= -1 # Use dot product to find the components of b1 and b3 that are not # perpendicular to b2. Subtract those components. The resulting vectors # lie in parallel planes. v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] ) # Use the relationship between cos and dot product to find the desired angle. return np.degrees(np.arccos( v[0].dot(v[1])/(np.linalg.norm(v[0]) * np.linalg.norm(v[1])))) # This is the straightforward approach as outlined in the answers to # "How do I calculate a dihedral angle given Cartesian coordinates?" def dihedral2(p): b = p[:-1] - p[1:] b[0] *= -1 v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] ) # Normalize vectors v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1) b1 = b[1] / np.linalg.norm(b[1]) x = np.dot(v[0], v[1]) m = np.cross(v[0], b1) y = np.dot(m, v[1]) return np.degrees(np.arctan2( y, x )) # This one starts with two cross products to get a vector perpendicular to # b2 and b1 and another perpendicular to b2 and b3. The angle between those vectors # is the dihedral angle. def dihedral3(p): b = p[:-1] - p[1:] b[0] *= -1 v = np.array( [np.cross(v,b[1]) for v in [b[0], b[2]] ] ) # Normalize vectors v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1) return np.degrees(np.arccos( v[0].dot(v[1]) )) dihedrals = [ ("dihedral1", dihedral1), ("dihedral2", dihedral2), ("dihedral3", dihedral3) ] 

Сравнительный анализ:

 # Testing arccos near 0 # Answer is 0.000057 p1 = np.array([ [ 1, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 1 ], [ 0.999999, 0.000001, 1 ] ]) # +x,+y p2 = np.array([ [ 1, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 1 ], [ 0.1, 0.6, 1 ] ]) # -x,+y p3 = np.array([ [ 1, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 1 ], [-0.3, 0.6, 1 ] ]) # -x,-y p4 = np.array([ [ 1, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 1 ], [-0.3, -0.6, 1 ] ]) # +x,-y p5 = np.array([ [ 1, 0, 0 ], [ 0, 0, 0 ], [ 0, 0, 1 ], [ 0.6, -0.6, 1 ] ]) for d in dihedrals: name = d[0] f = d[1] print "%s: %12.6f %12.6f %12.6f %12.6f %12.6f" \ % (name, f(p1), f(p2), f(p3), f(p4), f(p5)) print def profileDihedrals(f): t0 = time() for i in range(20000): p = np.random.random( (4,3) ) f(p) p = np.random.randn( 4,3 ) f(p) return(time() - t0) print "dihedral1: ", profileDihedrals(dihedral1) print "dihedral2: ", profileDihedrals(dihedral2) print "dihedral3: ", profileDihedrals(dihedral3) 

Результат тестирования:

 dihedral1: 0.000057 80.537678 116.565051 116.565051 45.000000 dihedral2: 0.000057 80.537678 116.565051 -116.565051 -45.000000 dihedral3: 0.000057 80.537678 116.565051 116.565051 45.000000 dihedral1: 2.79781794548 dihedral2: 3.74271392822 dihedral3: 2.49604296684 

Как вы можете видеть в бенчмаркинге, последний имеет тенденцию быть самым быстрым, а второй – единственным, который возвращает углы из полного диапазона от 0 до 2pi, поскольку он использует arctan2.

Вот реализация угла кручения по всему диапазону 2pi, который немного быстрее, не прибегает к numpy quirks (einsum таинственно быстрее, чем логически эквивалентный код), и его легче читать.

Там даже немного больше, чем просто хаки, происходящие здесь – математика тоже отличается. Формула, используемая в dihedral2 вопроса, использует 3 квадратных корня и 1 кросс-продукт, формула в Википедии использует 1 квадратный корень и 3 кросс-произведения, но в формуле, используемой в нижеприведенной функции, используется только 1 квадратный корень и 1 кросс-продукт. Это, вероятно, так же просто, как математика.

Функции с функцией диапазона 2pi из вопроса, формула Википедии для сравнения и новая функция:

dihedrals.py

 #!/usr/bin/env python # -*- coding: utf-8 -*- import numpy as np def old_dihedral2(p): """http://stackoverflow.com/q/20305272/1128289""" b = p[:-1] - p[1:] b[0] *= -1 v = np.array( [ v - (v.dot(b[1])/b[1].dot(b[1])) * b[1] for v in [b[0], b[2]] ] ) # Normalize vectors v /= np.sqrt(np.einsum('...i,...i', v, v)).reshape(-1,1) b1 = b[1] / np.linalg.norm(b[1]) x = np.dot(v[0], v[1]) m = np.cross(v[0], b1) y = np.dot(m, v[1]) return np.degrees(np.arctan2( y, x )) def wiki_dihedral(p): """formula from Wikipedia article on "Dihedral angle"; formula was removed from the most recent version of article (no idea why, the article is a mess at the moment) but the formula can be found in at this permalink to an old version of the article: https://en.wikipedia.org/w/index.php?title=Dihedral_angle&oldid=689165217#Angle_between_three_vectors uses 1 sqrt, 3 cross products""" p0 = p[0] p1 = p[1] p2 = p[2] p3 = p[3] b0 = -1.0*(p1 - p0) b1 = p2 - p1 b2 = p3 - p2 b0xb1 = np.cross(b0, b1) b1xb2 = np.cross(b2, b1) b0xb1_x_b1xb2 = np.cross(b0xb1, b1xb2) y = np.dot(b0xb1_x_b1xb2, b1)*(1.0/np.linalg.norm(b1)) x = np.dot(b0xb1, b1xb2) return np.degrees(np.arctan2(y, x)) def new_dihedral(p): """Praxeolitic formula 1 sqrt, 1 cross product""" p0 = p[0] p1 = p[1] p2 = p[2] p3 = p[3] b0 = -1.0*(p1 - p0) b1 = p2 - p1 b2 = p3 - p2 # normalize b1 so that it does not influence magnitude of vector # rejections that come next b1 /= np.linalg.norm(b1) # vector rejections # v = projection of b0 onto plane perpendicular to b1 # = b0 minus component that aligns with b1 # w = projection of b2 onto plane perpendicular to b1 # = b2 minus component that aligns with b1 v = b0 - np.dot(b0, b1)*b1 w = b2 - np.dot(b2, b1)*b1 # angle between v and w in a plane is the torsion angle # v and w may not be normalized but that's fine since tan is y/x x = np.dot(v, w) y = np.dot(np.cross(b1, v), w) return np.degrees(np.arctan2(y, x)) 

Новая функция, вероятно, будет более удобна для вызова с 4 отдельными аргументами, но для соответствия сигнатуре в исходном вопросе она сразу же распаковывает аргумент.

Код для тестирования:

test_dihedrals.ph

 from dihedrals import * # some atom coordinates for testing p0 = np.array([24.969, 13.428, 30.692]) # N p1 = np.array([24.044, 12.661, 29.808]) # CA p2 = np.array([22.785, 13.482, 29.543]) # C p3 = np.array([21.951, 13.670, 30.431]) # O p4 = np.array([23.672, 11.328, 30.466]) # CB p5 = np.array([22.881, 10.326, 29.620]) # CG p6 = np.array([23.691, 9.935, 28.389]) # CD1 p7 = np.array([22.557, 9.096, 30.459]) # CD2 # I guess these tests do leave 1 quadrant (-x, +y) untested, oh well... def test_old_dihedral2(): assert(abs(old_dihedral2(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4) assert(abs(old_dihedral2(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4) assert(abs(old_dihedral2(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4) assert(abs(old_dihedral2(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4) def test_new_dihedral1(): assert(abs(wiki_dihedral(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4) assert(abs(wiki_dihedral(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4) assert(abs(wiki_dihedral(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4) assert(abs(wiki_dihedral(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4) def test_new_dihedral2(): assert(abs(new_dihedral(np.array([p0, p1, p2, p3])) - (-71.21515)) < 1E-4) assert(abs(new_dihedral(np.array([p0, p1, p4, p5])) - (-171.94319)) < 1E-4) assert(abs(new_dihedral(np.array([p1, p4, p5, p6])) - (60.82226)) < 1E-4) assert(abs(new_dihedral(np.array([p1, p4, p5, p7])) - (-177.63641)) < 1E-4) 

Код для синхронизации:

time_dihedrals.py

 #!/usr/bin/env python # -*- coding: utf-8 -*- from dihedrals import * from time import time def profileDihedrals(f): t0 = time() for i in range(20000): p = np.random.random( (4,3) ) f(p) p = np.random.randn( 4,3 ) f(p) return(time() - t0) print("old_dihedral2: ", profileDihedrals(old_dihedral2)) print("wiki_dihedral: ", profileDihedrals(wiki_dihedral)) print("new_dihedral: ", profileDihedrals(new_dihedral)) 

Эти функции можно протестировать с помощью pytest как pytest ./test_dihedrals.py .

Сроки:

 ./time_dihedrals.py old_dihedral2: 1.6442952156066895 wiki_dihedral: 1.3895585536956787 new_dihedral: 0.8703620433807373 

new_dihedral примерно в два раза быстрее, чем old_dihedral2 .

… вы также можете видеть, что аппаратное обеспечение, используемое для этого ответа, намного сильнее, чем аппаратное обеспечение, используемое в вопросе (3.74 против 1.64 для dihedral2); -P

Если вы хотите стать еще более агрессивным, вы можете использовать pypy. Во время написания pypy не поддерживает numpy.cross но вы можете просто использовать кросс-продукт, реализованный в python. Для 3-векторного кросс-продукта, генерируемого Cpypy, вероятно, по меньшей мере так же хорошо, как то, что использует numpy. Для меня время доходит до 0.60 для меня, но при этом мы пробираемся в глупый хакс.

Тот же тест, но с тем же оборудованием, что и в вопросе:

 old_dihedral2: 3.0171279907226562 wiki_dihedral: 3.415065050125122 new_dihedral: 2.086946964263916 

Мой подход:

Составьте векторы b4 = b1 / \ b2 и b5 = b2 / \ b4. Они образуют ортогональный кадр с b2, а длина b5 – это b2, чем длина b4 (поскольку они ортогональны). Проецирование b3 на эти два вектора дает вам (масштабированные) 2D-координаты кончика b3, как на втором рисунке. Угол задается atan2 в диапазоне -Pi .. + Pi.

 b4= cross(b1, b2); b5= cross(b2, b4); Dihedral= atan2(dot(b3, b4), dot(b3, b5) * sqrt(dot(b2, b2))); 

Как и у вашего диэдра2. 12 добавляет, 21 умножает, 1 квадратный корень, 1 арктангенс. Вы можете переписать выражение для b5 с использованием формулы вытеснения, но это действительно не помогает.

ПРЕДОСТЕРЕЖЕНИЕ: Я не тщательно проверял проблемы с сигналами / квадрантами!

Вот окончательный ответ. У авторов есть версии python, а также версия C.

http://onlinelibrary.wiley.com/doi/10.1002/jcc.20237/abstract

Практическое преобразование из торсионного пространства в декартово пространство для синтеза силиконового белка

 First published: 16 May 2005