Coverage for /builds/ase/ase/ase/transport/tools.py: 86.96%

46 statements  

« prev     ^ index     » next       coverage.py v7.5.3, created at 2025-08-02 00:12 +0000

1# fmt: off 

2 

3import numpy as np 

4 

5 

6def dagger(matrix): 

7 return np.conj(matrix.T) 

8 

9 

10def rotate_matrix(h, u): 

11 return np.dot(u.T.conj(), np.dot(h, u)) 

12 

13 

14def get_subspace(matrix, index): 

15 """Get the subspace spanned by the basis function listed in index""" 

16 assert matrix.ndim == 2 and matrix.shape[0] == matrix.shape[1] 

17 return matrix.take(index, 0).take(index, 1) 

18 

19 

20def normalize(matrix, S=None): 

21 """Normalize column vectors. 

22 

23 :: 

24 

25 <matrix[:,i]| S |matrix[:,i]> = 1 

26 

27 """ 

28 for col in matrix.T: 

29 if S is None: 

30 col /= np.linalg.norm(col) 

31 else: 

32 col /= np.sqrt(np.dot(col.conj(), np.dot(S, col))) 

33 

34 

35def subdiagonalize(h_ii, s_ii, index_j): 

36 nb = h_ii.shape[0] 

37 nb_sub = len(index_j) 

38 h_sub_jj = get_subspace(h_ii, index_j) 

39 s_sub_jj = get_subspace(s_ii, index_j) 

40 e_j, v_jj = np.linalg.eig(np.linalg.solve(s_sub_jj, h_sub_jj)) 

41 normalize(v_jj, s_sub_jj) # normalize: <v_j|s|v_j> = 1 

42 permute_list = np.argsort(e_j.real) 

43 e_j = np.take(e_j, permute_list) 

44 v_jj = np.take(v_jj, permute_list, axis=1) 

45 

46 # Setup transformation matrix 

47 c_ii = np.identity(nb, complex) 

48 for i in range(nb_sub): 

49 for j in range(nb_sub): 

50 c_ii[index_j[i], index_j[j]] = v_jj[i, j] 

51 

52 h1_ii = rotate_matrix(h_ii, c_ii) 

53 s1_ii = rotate_matrix(s_ii, c_ii) 

54 

55 return h1_ii, s1_ii, c_ii, e_j 

56 

57 

58def cutcoupling(h, s, index_n): 

59 for i in index_n: 

60 s[:, i] = 0.0 

61 s[i, :] = 0.0 

62 s[i, i] = 1.0 

63 Ei = h[i, i] 

64 h[:, i] = 0.0 

65 h[i, :] = 0.0 

66 h[i, i] = Ei 

67 

68 

69def fermidistribution(energy, kt): 

70 # fermi level is fixed to zero 

71 # energy can be a single number or a list 

72 assert kt >= 0., 'Negative temperature encountered!' 

73 

74 if kt == 0: 

75 if isinstance(energy, float): 

76 return int(energy / 2. <= 0) 

77 else: 

78 return (energy / 2. <= 0).astype(int) 

79 else: 

80 return 1. / (1. + np.exp(energy / kt))