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

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

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

Dihedral angle1

Dihedral angle2

Что лучше всего для расчета углов в диапазоне от 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.

Ответы

Ответ 1

Здесь реализация угла кручения по всему диапазону 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 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 для меня, но при этом мы пробираемся в глупый hax.

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

old_dihedral2: 3.0171279907226562
wiki_dihedral: 3.415065050125122
new_dihedral:  2.086946964263916

Ответ 2

Мой подход:

Составьте векторы 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, используя формулу вытеснения, но это действительно не помогает.

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

Ответ 3

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

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

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

First published: 16 May 2005