Coverage for /builds/ase/ase/ase/stress.py: 93.02%

1# fmt: off 

2 

3import numpy as np 

4 

5# The indices of the full stiffness matrix of (orthorhombic) interest 

6voigt_notation = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)] 

7 

8 

9def get_elasticity_tensor(atoms, h=0.001, verbose=False): 

10 """ 

11 

12 1 dE dσ_ij 

13 C = - ----------- = ----- 

14 ijkl V dε_ij dε_kl dε_kl 

15 

16 """ 

17 cell0 = atoms.cell.copy() 

18 C_ijkl = np.zeros((3, 3, 3, 3)) 

19 f = voigt_6_to_full_3x3_stress 

20 for k in range(3): 

21 for l in range(3): 

22 strain = np.eye(3) 

23 strain[k, l] += h 

24 atoms.set_cell(cell0 @ strain, scale_atoms=True) 

25 stressp_ij = f(atoms.get_stress()) 

26 strain[k, l] -= 2 * h 

27 atoms.set_cell(cell0 @ strain, scale_atoms=True) 

28 stressm_ij = f(atoms.get_stress()) 

29 C_ijkl[k, l] = (stressp_ij - stressm_ij) / (2 * h) 

30 

31 if verbose: 

32 for i in range(3): 

33 for j in range(3): 

34 print(f'C_ijkl[{i}, {j}] =') 

35 for k in range(3): 

36 for l in range(3): 

37 print(round(C_ijkl[i, j, k, l], 2), end=' ') 

38 print() 

39 print() 

40 print() 

41 

42 return C_ijkl 

43 

44 

45def full_3x3_to_voigt_6_index(i, j): 

46 if i == j: 

47 return i 

48 return 6 - i - j 

49 

50 

51def voigt_6_to_full_3x3_strain(strain_vector): 

52 """ 

53 Form a 3x3 strain matrix from a 6 component vector in Voigt notation 

54 """ 

55 e1, e2, e3, e4, e5, e6 = np.transpose(strain_vector) 

56 return np.transpose([[1.0 + e1, 0.5 * e6, 0.5 * e5], 

57 [0.5 * e6, 1.0 + e2, 0.5 * e4], 

58 [0.5 * e5, 0.5 * e4, 1.0 + e3]]) 

59 

60 

61def voigt_6_to_full_3x3_stress(stress_vector): 

62 """ 

63 Form a 3x3 stress matrix from a 6 component vector in Voigt notation 

64 """ 

65 s1, s2, s3, s4, s5, s6 = np.transpose(stress_vector) 

66 return np.transpose([[s1, s6, s5], 

67 [s6, s2, s4], 

68 [s5, s4, s3]]) 

69 

70 

71def full_3x3_to_voigt_6_strain(strain_matrix): 

72 """ 

73 Form a 6 component strain vector in Voigt notation from a 3x3 matrix 

74 """ 

75 strain_matrix = np.asarray(strain_matrix) 

76 return np.transpose([strain_matrix[..., 0, 0] - 1.0, 

77 strain_matrix[..., 1, 1] - 1.0, 

78 strain_matrix[..., 2, 2] - 1.0, 

79 strain_matrix[..., 1, 2] + strain_matrix[..., 2, 1], 

80 strain_matrix[..., 0, 2] + strain_matrix[..., 2, 0], 

81 strain_matrix[..., 0, 1] + strain_matrix[..., 1, 0]]) 

82 

83 

84def full_3x3_to_voigt_6_stress(stress_matrix): 

85 """ 

86 Form a 6 component stress vector in Voigt notation from a 3x3 matrix 

87 """ 

88 stress_matrix = np.asarray(stress_matrix) 

89 return np.transpose([stress_matrix[..., 0, 0], 

90 stress_matrix[..., 1, 1], 

91 stress_matrix[..., 2, 2], 

92 (stress_matrix[..., 1, 2] + 

93 stress_matrix[..., 2, 1]) / 2, 

94 (stress_matrix[..., 0, 2] + 

95 stress_matrix[..., 2, 0]) / 2, 

96 (stress_matrix[..., 0, 1] + 

97 stress_matrix[..., 1, 0]) / 2])