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

1import numpy as np 

2 

3 

4def dagger(matrix): 

5 return np.conj(matrix.T) 

6 

7 

8def rotate_matrix(h, u): 

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

10 

11 

12def get_subspace(matrix, index): 

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

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

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

16 

17 

18def normalize(matrix, S=None): 

19 """Normalize column vectors. 

20 

21 :: 

22 

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

24 

25 """ 

26 for col in matrix.T: 

27 if S is None: 

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

29 else: 

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

31 

32 

33def subdiagonalize(h_ii, s_ii, index_j): 

34 nb = h_ii.shape[0] 

35 nb_sub = len(index_j) 

36 h_sub_jj = get_subspace(h_ii, index_j) 

37 s_sub_jj = get_subspace(s_ii, index_j) 

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

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

40 permute_list = np.argsort(e_j.real) 

41 e_j = np.take(e_j, permute_list) 

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

43 

44 # Setup transformation matrix 

45 c_ii = np.identity(nb, complex) 

46 for i in range(nb_sub): 

47 for j in range(nb_sub): 

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

49 

50 h1_ii = rotate_matrix(h_ii, c_ii) 

51 s1_ii = rotate_matrix(s_ii, c_ii) 

52 

53 return h1_ii, s1_ii, c_ii, e_j 

54 

55 

56def cutcoupling(h, s, index_n): 

57 for i in index_n: 

58 s[:, i] = 0.0 

59 s[i, :] = 0.0 

60 s[i, i] = 1.0 

61 Ei = h[i, i] 

62 h[:, i] = 0.0 

63 h[i, :] = 0.0 

64 h[i, i] = Ei 

65 

66 

67def fermidistribution(energy, kt): 

68 # fermi level is fixed to zero 

69 # energy can be a single number or a list 

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

71 

72 if kt == 0: 

73 if isinstance(energy, float): 

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

75 else: 

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

77 else: 

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