1import numpy as np 

2 

3 

4def rotation_matrix_from_points(m0, m1): 

5 """Returns a rigid transformation/rotation matrix that minimizes the 

6 RMSD between two set of points. 

7 

8 m0 and m1 should be (3, npoints) numpy arrays with 

9 coordinates as columns:: 

10 

11 (x1 x2 x3 ... xN 

12 y1 y2 y3 ... yN 

13 z1 z2 z3 ... zN) 

14 

15 The centeroids should be set to origin prior to 

16 computing the rotation matrix. 

17 

18 The rotation matrix is computed using quaternion 

19 algebra as detailed in:: 

20 

21 Melander et al. J. Chem. Theory Comput., 2015, 11,1055 

22 """ 

23 

24 v0 = np.copy(m0) 

25 v1 = np.copy(m1) 

26 

27 # compute the rotation quaternion 

28 

29 R11, R22, R33 = np.sum(v0 * v1, axis=1) 

30 R12, R23, R31 = np.sum(v0 * np.roll(v1, -1, axis=0), axis=1) 

31 R13, R21, R32 = np.sum(v0 * np.roll(v1, -2, axis=0), axis=1) 

32 

33 f = [[R11 + R22 + R33, R23 - R32, R31 - R13, R12 - R21], 

34 [R23 - R32, R11 - R22 - R33, R12 + R21, R13 + R31], 

35 [R31 - R13, R12 + R21, -R11 + R22 - R33, R23 + R32], 

36 [R12 - R21, R13 + R31, R23 + R32, -R11 - R22 + R33]] 

37 

38 F = np.array(f) 

39 

40 w, V = np.linalg.eigh(F) 

41 # eigenvector corresponding to the most 

42 # positive eigenvalue 

43 q = V[:, np.argmax(w)] 

44 

45 # Rotation matrix from the quaternion q 

46 

47 R = quaternion_to_matrix(q) 

48 

49 return R 

50 

51 

52def quaternion_to_matrix(q): 

53 """Returns a rotation matrix. 

54 

55 Computed from a unit quaternion Input as (4,) numpy array. 

56 """ 

57 

58 q0, q1, q2, q3 = q 

59 R_q = [[q0**2 + q1**2 - q2**2 - q3**2, 

60 2 * (q1 * q2 - q0 * q3), 

61 2 * (q1 * q3 + q0 * q2)], 

62 [2 * (q1 * q2 + q0 * q3), 

63 q0**2 - q1**2 + q2**2 - q3**2, 

64 2 * (q2 * q3 - q0 * q1)], 

65 [2 * (q1 * q3 - q0 * q2), 

66 2 * (q2 * q3 + q0 * q1), 

67 q0**2 - q1**2 - q2**2 + q3**2]] 

68 return np.array(R_q) 

69 

70 

71def minimize_rotation_and_translation(target, atoms): 

72 """Minimize RMSD between atoms and target. 

73 

74 Rotate and translate atoms to best match target. For more details, see:: 

75 

76 Melander et al. J. Chem. Theory Comput., 2015, 11,1055 

77 """ 

78 

79 p = atoms.get_positions() 

80 p0 = target.get_positions() 

81 

82 # centeroids to origin 

83 c = np.mean(p, axis=0) 

84 p -= c 

85 c0 = np.mean(p0, axis=0) 

86 p0 -= c0 

87 

88 # Compute rotation matrix 

89 R = rotation_matrix_from_points(p.T, p0.T) 

90 

91 atoms.set_positions(np.dot(p, R.T) + c0)