Coverage for calorine/tools/stiffness.py: 100%

1from itertools import combinations_with_replacement, product 

2import numpy as np 

3from ase import Atoms 

4from ase.units import GPa 

5from .structures import relax_structure 

6 

7 

8def get_elastic_stiffness_tensor(structure: Atoms, 

9 clamped: bool = False, 

10 epsilon: float = 1e-3, 

11 **kwargs) -> np.ndarray: 

12 """Calculate and return the elastic stiffness tensor in units of GPa for the 

13 given structure in `Voigt form <https://en.wikipedia.org/wiki/Voigt_notation>`_. 

14 

15 Parameters 

16 ---------- 

17 structure 

18 input structure; should be fully relaxed 

19 clamped 

20 if ``False`` (default) return the *relaxed* elastic stiffness tensor; 

21 if ``True`` return the *clamped ion* elastic stiffness tensor 

22 epsilon 

23 magnitude of the applied strain 

24 kwargs 

25 keyword arguments forwarded to the :func:`relax_structure 

26 <calorine.tools.relax_structure>` function used for relaxing the 

27 structure when computing the relaxed stiffness tensor; it should not be 

28 necessary to change the default for the vast majority of use cases; use 

29 with care 

30 

31 Returns 

32 ------- 

33 Stiffness tensor in units of GPa 

34 """ 

35 

36 # set up of deformations 

37 deformations = [] 

38 for i, j in combinations_with_replacement(range(9), r=2): 

39 for s1, s2 in product([-1, 1], repeat=2): 

40 S = np.zeros((3, 3)) 

41 S.flat[i] = s1 

42 S.flat[j] = s2 

43 deformations.append(S) 

44 deformations = np.array(deformations) 

45 deformations *= epsilon 

46 

47 # compute strain energies 

48 reference_energy = structure.get_potential_energy() 

49 energies = [] 

50 for S in deformations: 

51 cell = structure.get_cell() 

52 cell += cell @ S.T 

53 deformed_structure = structure.copy() 

54 deformed_structure.calc = structure.calc 

55 deformed_structure.set_cell(cell, scale_atoms=True) 

56 if not clamped: 

57 relax_structure(deformed_structure, constant_cell=True, **kwargs) 

58 energy = deformed_structure.get_potential_energy() 

59 energies.append(energy - reference_energy) 

60 energies = np.array(energies) 

61 

62 # extract stiffness tensor (full rank) 

63 SS = np.einsum('nij,nkl->nijkl', deformations, deformations) 

64 M = SS.reshape(len(SS), -1) 

65 M *= 0.5 

66 C, *_ = np.linalg.lstsq(M, energies, rcond=None) 

67 C = C.reshape(3, 3, 3, 3) 

68 C /= (structure.cell.volume * GPa) 

69 

70 # convert stiffness tensor to Voigt form 

71 voigts = np.array([1, 1, 2, 2, 3, 3, 2, 3, 3, 1, 1, 2]).reshape(-1, 2) - 1 

72 Cv = np.zeros((6, 6)) 

73 for i, j in product(range(6), repeat=2): 

74 v1 = voigts[i] 

75 v2 = voigts[j] 

76 Cv[i, j] = C[(*v1, *v2)] 

77 

78 return Cv