What suggestions do people have for quickly calculating dihedral angles in Python?
In the diagrams, phi is the dihedral angle:
What's your best for calculating angles in the range 0 to pi? What about 0 to 2pi?
"Best" here means some mix of fast and numerically stable. Methods that return values over the full range 0 to 2pi are preferred but if you have an incredibly fast way of calculating the dihedral over 0 to pi share that too.
Here are my 3 best efforts. Only the 2nd one returns angles between 0 and 2pi. It's also the slowest.
General comments about my approaches:
arccos() in Numpy seems plenty stable but since people raise this issue I may just not fully understand it.
The use of einsum came from here. Why is numpy's einsum faster than numpy's built in functions?
The diagrams and some inspiration came from here. How do I calculate a dihedral angle given Cartesian coordinates?
The 3 approaches with comments:
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) ]
Benchmarking:
# 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)
Benchmarking output:
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
As you can see in the benchmarking, the last one tends to be the fastest while the second one is the only one that returns angles from the full range of 0 to 2pi since it uses arctan2.
