Merge d1b0367 into f19873e

pyiron_atomistics/atomistics/master/elastic.py

@@ -147,11 +147,9 @@ def _fit_coeffs_with_energies(
     energy = np.tile(energy, len(rotations))
     volume = np.tile(volume, len(rotations))
     # Create symmetric tensor for elastic tensor
-    strain = 0.5 * np.einsum("ni,nj->nij", strain_voigt, strain_voigt)
-    # Set lower triangle to 0 (which is the same as the upper triangle)
-    strain = np.triu(strain).reshape(-1, 36)
+    strain = 0.5 * np.einsum("ni,nj->nij", *2 * [strain_voigt])
     # Remove lower triangle
-    strain = strain[:, np.sum(strain, axis=0) != 0]
+    strain = np.array([s[np.triu_indices(6)] for s in strain])
     if higher_strains is not None:
         strain = np.concatenate((strain, higher_strains), axis=-1)
     if fit_first_order:
@@ -160,9 +158,9 @@ def _fit_coeffs_with_energies(
     reg = LinearRegression().fit(strain, energy)
     score = reg.score(strain, energy)
     # Create base tensor for elastic tensor
-    coeff = np.triu(np.ones((6, 6))).flatten()
+    coeff = np.zeros((6, 6))
     # Multiply upper triangle with upper triangle coeffs (v.s.)
-    coeff[coeff != 0] *= reg.coef_[:21] * eV_div_A3_to_GPa
+    coeff[np.triu_indices(6)] = reg.coef_[:21] * eV_div_A3_to_GPa
     coeff = coeff.reshape(6, 6)
     coeff = 0.5 * (coeff + coeff.T)
     return coeff, score

tests/atomistics/master/test_elastic.py

@@ -6,8 +6,13 @@
 import unittest
 from pyiron_atomistics.atomistics.structure.atoms import CrystalStructure
 from pyiron_base import Project
-from pyiron_atomistics.atomistics.master.elastic import (calc_elastic_tensor, _get_higher_order_strains,
-    calc_elastic_constants, get_elastic_tensor_by_orientation)
+from pyiron_atomistics.atomistics.master.elastic import (
+    calc_elastic_tensor,
+    _get_higher_order_strains,
+    calc_elastic_constants,
+    get_elastic_tensor_by_orientation,
+    _convert_to_voigt,
+)
 import numpy as np
 
 
@@ -88,6 +93,11 @@ def test_calc_elastic_constants(self):
             np.linalg.norm(C-get_elastic_tensor_by_orientation([[1, 0, 0], [0, 0, -1], [0, 1, 0]], C)), 0
         )
 
+    def test_convert_to_voigt(self):
+        voigt = _convert_to_voigt(np.arange(3**4).reshape(*4 * (3,)))
+        self.assertTrue(np.array_equal(voigt[0], [ 0,  4,  8,  5,  2,  1]))
+        self.assertTrue(np.all(np.diff(voigt, axis=0) == 9))
+
 
 if __name__ == "__main__":
     unittest.main()