Coverage for gpaw/rotation.py: 65%

1# Copyright (C) 2003 CAMP 

2# Please see the accompanying LICENSE file for further information. 

3 

4from math import sqrt, cos, sin 

5 

6import numpy as np 

7 

8from gpaw.spherical_harmonics import Y 

9 

10 

11s = sqrt(0.5) 

12t = sqrt(3) / 3 

13# Points on the unit sphere: 

14sphere_lm = [ 

15 np.array([(1, 0, 0)]), # s 

16 np.array([(1, 0, 0), (0, 1, 0), (0, 0, 1)]), # p 

17 np.array([(s, s, 0), (0, s, s), (s, 0, s), (1, 0, 0), (0, 0, 1)]), # d 

18 np.array([(s, s, 0), (0, s, s), (s, 0, s), 

19 (1, 0, 0), (s, -s, 0), 

20 (t, -t, t), (t, t, -t)])] # f 

21 

22 

23def Y_matrix(l, U_vv): 

24 """YMatrix(l, symmetry) -> matrix. 

25 

26 For 2*l+1 points on the unit sphere (m1=0,...,2*l) calculate the 

27 value of Y_lm2 for m2=0,...,2*l. The points are those from the 

28 list sphere_lm[l] rotated as described by U_vv.""" 

29 

30 Y_mm = np.zeros((2 * l + 1, 2 * l + 1)) 

31 for m1, point in enumerate(sphere_lm[l]): 

32 x, y, z = np.dot(point, U_vv) 

33 for m2 in range(2 * l + 1): 

34 L = l**2 + m2 

35 Y_mm[m1, m2] = Y(L, x, y, z) 

36 return Y_mm 

37 

38 

39identity = np.eye(3) 

40iY_lmm = [np.linalg.inv(Y_matrix(l, identity)) for l in range(4)] 

41 

42 

43def rotation(l, U_vv): 

44 """Rotation(l, symmetry) -> transformation matrix. 

45 

46 Find the transformation from Y_lm1 to Y_lm2.""" 

47 

48 return np.dot(iY_lmm[l], Y_matrix(l, U_vv)) 

49 

50 

51def Y_rotation(l, angle): 

52 Y_mm = np.zeros((2 * l + 1, 2 * l + 1)) 

53 sn = sin(angle) 

54 cs = cos(angle) 

55 for m1, point in enumerate(sphere_lm[l]): 

56 x = point[0] 

57 y = cs * point[1] - sn * point[2] 

58 z = sn * point[1] + cs * point[2] 

59 for m2 in range(2 * l + 1): 

60 L = l**2 + m2 

61 Y_mm[m1, m2] = Y(L, x, y, z) 

62 return Y_mm 

63 

64 

65del s, identity