对于在 Python 中快速计算二面角，人们有什么建议？
在图中，phi 是二面角:
在计算 0 到 pi 范围内的角度时，您最擅长的是什么？那么0到2pi呢？
这里的“最佳”意味着快速和数值稳定的混合。首选在 0 到 2pi 的整个范围内返回值的方法，但如果您有一种非常快速的方法来计算 0 到 pi 的二面角，那么也可以分享。
这是我的 3 项最佳努力。只有第二个返回 0 到 2pi 之间的角度。它也是最慢的。
关于我的方法的一般评论:
Numpy 中的arccos() 看起来很稳定，但由于人们提出这个问题，我可能只是不完全理解它。
einsum 的使用来自这里。 Why is numpy's einsum faster than numpy's built in functions?
图表和一些灵感来自这里。 How do I calculate a dihedral angle given Cartesian coordinates?
带注释的 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 怪癖(einsum 神秘地比逻辑等效代码快)，并且更易于阅读。
这里不仅仅是黑客攻击，数学也不同。问题的 dihedral2 中使用的公式使用 3 个平方根和 1 个叉积，Wikipedia 上的公式使用 1 个平方根和 3 个叉积，但下面函数中使用的公式仅使用 1 个平方根和1个交叉产品。这可能就像数学可以得到的一样简单。
来自问题的 2pi 范围函数、用于比较的维基百科公式和新函数的函数:
二面体.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 的两倍。
...您还可以看到，用于此答案的硬件比问题中使用的硬件强大得多(dihedral2 为 3.74 与 1.64)；-P
如果你想变得更有侵略性，你可以使用 pypy。在撰写本文时 pypy 不支持 numpy.cross 但您可以只使用在 python 中实现的交叉产品。对于 3 向量叉积，C pypy 生成的可能至少与 numpy 使用的一样好。这样做可以将我的时间缩短到 0.60，但此时我们正在陷入愚蠢的骗局。
与问题中使用的基准相同但硬件相同:
old_dihedral2: 3.0171279907226562
wiki_dihedral: 3.415065050125122
new_dihedral: 2.086946964263916
关于python - Python中笛卡尔坐标中四个点的二面角/扭转角，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/20305272/