Python中笛卡尔坐标中四点的二面角/扭转角

Question

人们对在Python中快速计算二面角有什么建议？

在图中,phi是二面角:

在0到pi范围内计算角度的最佳方法是什么？0到2pi怎么样？

这里的"最佳"意味着快速和数字稳定的混合.在0到2pi的整个范围内返回值的方法是首选,但如果您有一种非常快速的方法来计算0到pi的二面角,那么也是如此.

这是我最好的3个努力.只有第二个返回0到2pi之间的角度.它也是最慢的.

关于我的方法的一般评论:

Numpy的arccos()似乎很稳定但是由于人们提出这个问题,我可能还没有完全理解它.

使用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.

Answer 1

这是一个在整个2pi范围内扭转角度的实现,速度更快,不会诉诸于numpy怪癖(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 vs 1.64 for dihedral2)更强大;-P

如果你想要更加激进,你可以使用pypy.在撰写本文时,pypy不支持,numpy.cross但您可以使用python中实现的交叉产品.对于3载体交叉产物,C pypy产生的可能性至少与numpy使用的产品一样好.这样做可以让我的时间降到0.60,但在此我们正在涉及愚蠢的hax.

相同的基准,但使用与问题中使用的相同的硬件:

old_dihedral2: 3.0171279907226562
wiki_dihedral: 3.415065050125122
new_dihedral:  2.086946964263916