Coverage for calorine/tools/stiffness.py: 100%
41 statements
« prev ^ index » next coverage.py v7.2.5, created at 2024-07-26 09:34 +0000
« prev ^ index » next coverage.py v7.2.5, created at 2024-07-26 09:34 +0000
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
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>`_.
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
31 Returns
32 -------
33 Stiffness tensor in units of GPa
34 """
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
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)
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)
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)]
78 return Cv