新手求助：python求二面角代码应该怎么写

Jerryluo

请问用python语言读取txt文件中的四个原子的xyz坐标，并求二面角代码应该怎么写，谢谢。比如txt文件中乙烷的二面角：H3-C1-C2-H4
6
C2H6
H1                 -1.51326409   -0.45762891   -0.02447848
H2                 -1.51324568    1.05557928   -0.89812999
H3                 -2.93991852    0.55119428   -0.02447848
C1                 -1.86991852    0.55118109   -0.02447848
C2                 -1.35657630    1.27713737    1.23292649
H4                 -1.71163265    2.28651046    1.23194897
H5                 -1.71484541    0.77386893    2.10657638
H6                 -0.28657811    1.27543068    1.23390524

zjxitcc

可参考实例
https://dasher.wustl.edu/tinker/distribution/source/geometry.f

Mikasa

供参考：

1. #!/usr/bin/env python
   
2. 
   
3. """Dihedral angle calculation module.
   
4. K. Rother
   
5. Two functions are provided (see documentation of each function for
   
6. a more detailed description):
   
7. angle (v1,v2) - returns the angle between two 2D or 3D numpy arrays
   
8.     (in radians)
   
9. dihedral (v1,v2,v3,v4) - returns the dihedral angle between four 3D vectors
   
10.     (in degrees)
    
11. The vectors that dihedral uses can be lists, tuples or numpy arrays.
    
12. Scientific.Geometry.Vector objects behave differently on Win and Linux,
    
13. and they are therefore not supported (but may work anyway).
    
14. """
    
15. 
    
16. __author__ = "Kristian Rother"
    
17. __copyright__ = "Copyright 2007-2016, The Cogent Project"
    
18. __contributors__ = ["Kristian Rother", "Sandra Smit"]
    
19. __credits__ = ["Janusz Bujnicki", "Daniel McDonald"]
    
20. __license__ = "GPL"
    
21. __version__ = "1.9"
    
22. __maintainer__ = "Kristian Rother"
    
23. __email__ = "krother@rubor.de"
    
24. __status__ = "Production"
    
25. 
    
26. from numpy import array, cross, pi, cos, sqrt, sum, arccos as acos
    
27. 
    
28. 
    
29. class DihedralGeometryError(Exception):
    
30.     pass
    
31. 
    
32. 
    
33. class AngleGeometryError(Exception):
    
34.     pass
    
35. 
    
36. 
    
37. def fix_rounding_error(x):
    
38.     ROUND_ERROR = 1e-14    # fp rounding error: causes some tests to fail
    
39.     # will round to 0 if smaller in magnitude than this
    
40.     """
    
41.     If x is almost in the range 0-1, fixes it.
    
42.     Specifically, if x is between -ROUND_ERROR and 0, returns 0.
    
43.     If x is between 1 and 1+ROUND_ERROR, returns 1.
    
44.     """
    
45.     if -ROUND_ERROR < x < 0:
    
46.         return 0
    
47.     elif 1 < x < 1 + ROUND_ERROR:
    
48.         return 1
    
49.     else:
    
50.         return x
    
51. 
    
52. 
    
53. def norm(a):
    
54.     """Returns the norm of a matrix or vector
    
55. 
    
56.     Calculates the Euclidean norm of a vector.
    
57.     Applies the Frobenius norm function to a matrix
    
58.     (a.k.a. Euclidian matrix norm)
    
59. 
    
60.     a = numpy array
    
61.     """
    
62.     return sqrt(sum((a * a).flat))
    
63. 
    
64. 
    
65. def scalar(v1, v2):
    
66.     """
    
67.     calculates the scalar product of two vectors
    
68.     v1 and v2 are numpy.array objects.
    
69.     returns a float for a one-dimensional array.
    
70.     """
    
71.     return sum(v1 * v2)
    
72. 
    
73. 
    
74. def angle(v1, v2):
    
75.     """
    
76.     calculates the angle between two vectors.
    
77.     v1 and v2 are numpy.array objects.
    
78.     returns a float containing the angle in radians.
    
79.     """
    
80.     length_product = norm(v1) * norm(v2)
    
81.     if length_product == 0:
    
82.         raise AngleGeometryError(
    
83.             "Cannot calculate angle for vectors with length zero")
    
84.     cosine = scalar(v1, v2) / length_product
    
85.     angle = acos(fix_rounding_error(cosine))
    
86.     return angle
    
87. 
    
88. 
    
89. def calc_angle(vec1, vec2, vec3):
    
90.     """Calculates a flat angle from three coordinates."""
    
91.     if len(vec1) == 3:
    
92.         v1, v2, v3 = map(create_vector, [vec1, vec2, vec3])
    
93.     else:
    
94.         v1, v2, v3 = map(create_vector2d, [vec1, vec2, vec3])
    
95.     v12 = v2 - v1
    
96.     v23 = v2 - v3
    
97.     angl = angle(v12, v23) * 180.0 / pi
    
98.     return angl
    
99. 
    
100. 
     
101. def create_vector2d(vec):
     
102.     """Returns a vector as a numpy array."""
     
103.     return array([vec[0], vec[1]])
     
104. 
     
105. 
     
106. def create_vector(vec):
     
107.     """Returns a vector as a numpy array."""
     
108.     return array([vec[0], vec[1], vec[2]])
     
109. 
     
110. 
     
111. def create_vectors(vec1, vec2, vec3, vec4):
     
112.     """Returns dihedral angle, takes four
     
113.     Scientific.Geometry.Vector objects
     
114.     (dihedral does not work for them because
     
115.     the Win and Linux libraries are not identical.
     
116.     """
     
117.     return map(create_vector, [vec1, vec2, vec3, vec4])
     
118. 
     
119. 
     
120. def dihedral(vec1, vec2, vec3, vec4):
     
121.     """
     
122.     Returns a float value for the dihedral angle between
     
123.     the four vectors. They define the bond for which the
     
124.     torsion is calculated (~) as:
     
125.     V1 - V2 ~ V3 - V4
     
126.     The vectors vec1 .. vec4 can be array objects, lists or tuples of length
     
127.     three containing floats.
     
128.     For Scientific.geometry.Vector objects the behavior is different
     
129.     on Windows and Linux. Therefore, the latter is not a featured input type
     
130.     even though it may work.
     
131. 
     
132.     If the dihedral angle cant be calculated (because vectors are collinear),
     
133.     the function raises a DihedralGeometryError
     
134.     """
     
135.     # create array instances.
     
136.     v1, v2, v3, v4 = create_vectors(vec1, vec2, vec3, vec4)
     
137.     all_vecs = [v1, v2, v3, v4]
     
138. 
     
139.     # rule out that two of the atoms are identical
     
140.     # except the first and last, which may be.
     
141.     for i in range(len(all_vecs) - 1):
     
142.         for j in range(i + 1, len(all_vecs)):
     
143.             if i > 0 or j < 3:  # exclude the (1,4) pair
     
144.                 equals = all_vecs[i] == all_vecs[j]
     
145.                 if equals.all():
     
146.                     raise DihedralGeometryError(
     
147.                         "Vectors #%i and #%i may not be identical!" % (i, j))
     
148. 
     
149.     # calculate vectors representing bonds
     
150.     v12 = v2 - v1
     
151.     v23 = v3 - v2
     
152.     v34 = v4 - v3
     
153. 
     
154.     # calculate vectors perpendicular to the bonds
     
155.     normal1 = cross(v12, v23)
     
156.     normal2 = cross(v23, v34)
     
157. 
     
158.     # check for linearity
     
159.     if norm(normal1) == 0 or norm(normal2) == 0:
     
160.         raise DihedralGeometryError(
     
161.             "Vectors are in one line; cannot calculate normals!")
     
162. 
     
163.     # normalize them to length 1.0
     
164.     normal1 = normal1 / norm(normal1)
     
165.     normal2 = normal2 / norm(normal2)
     
166. 
     
167.     # calculate torsion and convert to degrees
     
168.     torsion = angle(normal1, normal2) * 180.0 / pi
     
169. 
     
170.     # take into account the determinant
     
171.     # (the determinant is a scalar value distinguishing
     
172.     # between clockwise and counter-clockwise torsion.
     
173.     if scalar(normal1, v34) >= 0:
     
174.         return torsion
     
175.     else:
     
176.         torsion = -torsion
     
177.         if torsion == 360:
     
178.             torsion = 0.0
     
179.         return torsion
     

复制代码

Daniel_Arndt

我之前写过一段awk代码，比较粗糙。你可以参考一下。

function Dihedral(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,    norm1,x5,y5,z5,norm2,x6,y6,z6)
{
    norm1=sqrt((y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))*(y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))+(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))*(z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))+(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))*(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2)))
    x5 = (y1*(z2-z3)+y2*(z3-z1)+y3*(z1-z2))/norm1
    y5 = (z1*(x2-x3)+z2*(x3-x1)+z3*(x1-x2))/norm1
    z5 = (x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/norm1
    norm2=sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2))
    x6 = (y5*(z3-z2)-(y3-y2)*z5)/norm2
    y6 = (z5*(x3-x2)-(z3-z2)*x5)/norm2
    z6 = (x5*(y3-y2)-(x3-x2)*y5)/norm2
    return -atan2(x6*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y6*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z6*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)),x5*(y2*(z3-z4)+y3*(z4-z2)+y4*(z2-z3))+y5*(z2*(x3-x4)+z3*(x4-x2)+z4*(x2-x3))+z5*(x2*(y3-y4)+x3*(y4-y2)+x4*(y2-y3)))
}

让你变成回忆

基本思想就是求出两个平面的法向量（定义平面内两个向量，然后用叉乘算出法向量），再求两个法向量之间的夹角。

k64_cc

def dihedral(a, b, c, d):
    """
    Calculate the dihedral angle between abc surface and bcd surface (- pi ~ pi).
    """
    na = np.cross(c - b, a - b)
    nb = np.cross(b - c, d - c)
    return np.arccos(np.abs(np.dot(na, nb)) / np.sqrt(np.power(na, 2).sum()) / np.sqrt(np.power(nb, 2).sum()))

没查过边界值，凑合用的。

highlight

搬运 https://stackoverflow.com/questions/20305272/dihedral-torsion-angle-from-four-points-in-cartesian-coordinates-in-python
里面用 wiki 定义的那个

1. 
   
2. #!/bin/env python3.6
   
3. 
   
4. import numpy as np
   
5. 
   
6. 
   
7. def read_file(filename):
   
8.     with open(filename, 'r') as f:
   
9.         natom = int(f.readline())
   
10.         f.readline()
    
11.         coordinates = []
    
12.         for _ in range(natom):
    
13.             line = f.readline()
    
14.             coordinate = line.split()[1:]
    
15.             coordinates.append(coordinate)
    
16.     return np.array(coordinates, dtype=float)
    
17. 
    
18. 
    
19. def dihedral(coordinates, a, b, c, d):
    
20.     i = coordinates[a]
    
21.     j = coordinates[b]
    
22.     k = coordinates[c]
    
23.     l = coordinates[d]
    
24.     f, g, h = i-j, j-k, l-k
    
25.     a = np.cross(f, g)
    
26.     b = np.cross(h, g)
    
27.     axb = np.cross(a, b)
    
28.     cos = np.dot(a, b)
    
29.     sin = np.dot(axb, g) /  np.linalg.norm(g)
    
30.     r = -np.arctan2(sin, cos)
    
31.     return np.degrees(r)
    
32. 
    
33. 
    
34. if __name__ == "__main__":
    
35.     c = read_file('C2H6.xyz')
    
36.     print(dihedral(c, 2, 3, 4, 5))
    

复制代码

Jerryluo

十分感谢大家！

我本是个娃娃

diangle.py (1.28 KB, 下载次数 Times of downloads: 41)

2019-7-15 13:01 上传 Uploaded

点击下载
Click to download

based on the blog: http://sobereva.com/302

you can modify according to the case you faced.