From ec750ca15d02cdd51b0c0a7a4408af8e0d259223 Mon Sep 17 00:00:00 2001 From: Thomas Purcell Date: Tue, 21 Nov 2023 10:42:53 -0600 Subject: [PATCH] Convert all FHI-aims stresses to be 3x3 instead of Voigt notation (#3476) * Convert all FHI-aims stresses to be 3x3 Previously used voigt * Update type hinting for parsers.py Use np array instead of Sequence Signed-off-by: Thomas Purcell * fix typo, snake case var name * Use Tensor.from_voigt instead of my own function This should now be more consistent --------- Signed-off-by: Thomas Purcell Co-authored-by: Janosh Riebesell Co-authored-by: Thomas Purcell --- pymatgen/io/aims/parsers.py | 21 +++++----- tests/io/aims/test_aims_parsers.py | 62 +++++++++--------------------- 2 files changed, 28 insertions(+), 55 deletions(-) diff --git a/pymatgen/io/aims/parsers.py b/pymatgen/io/aims/parsers.py index 7f577e0aa9b..daa3c686c7d 100644 --- a/pymatgen/io/aims/parsers.py +++ b/pymatgen/io/aims/parsers.py @@ -9,6 +9,7 @@ import numpy as np from pymatgen.core import Lattice, Molecule, Structure +from pymatgen.core.tensors import Tensor if TYPE_CHECKING: from collections.abc import Generator, Sequence @@ -40,7 +41,7 @@ def __init__(self, message: str) -> None: # Read aims.out files -scalar_property_to_line_key = { +SCALAR_PROPERTY_TO_LINE_KEY = { "free_energy": ["| Electronic free energy"], "number_of_iterations": ["| Number of self-consistency cycles"], "magnetic_moment": ["N_up - N_down"], @@ -110,7 +111,7 @@ def parse_scalar(self, property: str) -> float | None: Returns: The scalar value of the property or None if not found """ - line_start = self.reverse_search_for(scalar_property_to_line_key[property]) + line_start = self.reverse_search_for(SCALAR_PROPERTY_TO_LINE_KEY[property]) if line_start == LINE_NOT_FOUND: return None @@ -360,7 +361,7 @@ def n_atoms(self) -> int: @property def n_bands(self) -> int | None: """The number of Kohn-Sham states for the chunk.""" - line_start = self.reverse_search_for(scalar_property_to_line_key["n_bands"]) + line_start = self.reverse_search_for(SCALAR_PROPERTY_TO_LINE_KEY["n_bands"]) if line_start == LINE_NOT_FOUND: raise AimsParseError("No information about the number of Kohn-Sham states in the header.") @@ -374,7 +375,7 @@ def n_bands(self) -> int | None: @property def n_electrons(self) -> int | None: """The number of electrons for the chunk.""" - line_start = self.reverse_search_for(scalar_property_to_line_key["n_electrons"]) + line_start = self.reverse_search_for(SCALAR_PROPERTY_TO_LINE_KEY["n_electrons"]) if line_start == LINE_NOT_FOUND: raise AimsParseError("No information about the number of electrons in the header.") @@ -402,7 +403,7 @@ def n_spins(self) -> int | None: @property def electronic_temperature(self) -> float: """The electronic temperature for the chunk.""" - line_start = self.reverse_search_for(scalar_property_to_line_key["electronic_temp"]) + line_start = self.reverse_search_for(SCALAR_PROPERTY_TO_LINE_KEY["electronic_temp"]) # TARP: Default FHI-aims value if line_start == LINE_NOT_FOUND: return 0.00 @@ -620,7 +621,7 @@ def lattice(self) -> Lattice: return self._cache["lattice"] @property - def forces(self) -> list[Vector3D] | None: + def forces(self) -> np.array[Vector3D] | None: """The forces from the aims.out file.""" line_start = self.reverse_search_for(["Total atomic forces"]) if line_start == LINE_NOT_FOUND: @@ -633,7 +634,7 @@ def forces(self) -> list[Vector3D] | None: ) @property - def stresses(self) -> list[Matrix3D] | None: + def stresses(self) -> np.array[Matrix3D] | None: """The stresses from the aims.out file and convert to kbar.""" line_start = self.reverse_search_for(["Per atom stress (eV) used for heat flux calculation"]) if line_start == LINE_NOT_FOUND: @@ -642,15 +643,13 @@ def stresses(self) -> list[Matrix3D] | None: stresses = [] for line in self.lines[line_start : line_start + self.n_atoms]: xx, yy, zz, xy, xz, yz = (float(d) for d in line.split()[2:8]) - stresses.append([xx, yy, zz, yz, xz, xy]) + stresses.append(Tensor.from_voigt([xx, yy, zz, yz, xz, xy])) return np.array(stresses) * EV_PER_A3_TO_KBAR @property def stress(self) -> Matrix3D | None: """The stress from the aims.out file and convert to kbar.""" - from ase.stress import full_3x3_to_voigt_6_stress - line_start = self.reverse_search_for( [ "Analytical stress tensor - Symmetrized", @@ -661,7 +660,7 @@ def stress(self) -> Matrix3D | None: return None stress = [[float(inp) for inp in line.split()[2:5]] for line in self.lines[line_start + 5 : line_start + 8]] - return full_3x3_to_voigt_6_stress(stress) * EV_PER_A3_TO_KBAR + return np.array(stress) * EV_PER_A3_TO_KBAR @property def is_metallic(self) -> bool: diff --git a/tests/io/aims/test_aims_parsers.py b/tests/io/aims/test_aims_parsers.py index 9916eca1feb..6cecc3e4ff9 100644 --- a/tests/io/aims/test_aims_parsers.py +++ b/tests/io/aims/test_aims_parsers.py @@ -8,14 +8,14 @@ EV_PER_A3_TO_KBAR, LINE_NOT_FOUND, ) +from pymatgen.core.tensors import Tensor import gzip from pathlib import Path -from ase.stress import full_3x3_to_voigt_6_stress import pytest -eps_hp = 1e-15 # The espsilon value used to compare numbers that are high-precision -eps_lp = 1e-7 # The espsilon value used to compare numbers that are low-precision +eps_hp = 1e-15 # The epsilon value used to compare numbers that are high-precision +eps_lp = 1e-7 # The epsilon value used to compare numbers that are low-precision parser_file_dir = Path(__file__).parent / "parser_checks" @@ -52,10 +52,10 @@ def empty_header_chunk(): return AimsOutHeaderChunk([]) -@pytest.mark.parametrize("attrname", ["n_atoms", "n_bands", "n_electrons", "n_spins", "initial_structure"]) -def test_missing_parameter(attrname, empty_header_chunk): +@pytest.mark.parametrize("attr_name", ["n_atoms", "n_bands", "n_electrons", "n_spins", "initial_structure"]) +def test_missing_parameter(attr_name, empty_header_chunk): with pytest.raises(AimsParseError, match="No information about"): - getattr(empty_header_chunk, attrname) + getattr(empty_header_chunk, attr_name) def test_default_header_electronic_temperature(empty_header_chunk): @@ -90,9 +90,9 @@ def test_default_header_k_point_weights(empty_header_chunk): def initial_lattice(): return np.array( [ - [1.00000000, 2.70300000, 3.70300000], - [4.70300000, 2.00000000, 6.70300000], - [8.70300000, 7.70300000, 3.00000000], + [1, 2.70300000, 3.70300000], + [4.70300000, 2, 6.70300000], + [8.70300000, 7.70300000, 3], ] ) @@ -329,8 +329,8 @@ def test_calc_forces(calc_chunk): def test_calc_stresses(calc_chunk): stresses = EV_PER_A3_TO_KBAR * np.array( [ - [-10.0, -20.0, -30.0, -60.0, -50.0, -40.0], - [10.0, 20.0, 30.0, 60.0, 50.0, 40.0], + Tensor.from_voigt([-10.0, -20.0, -30.0, -60.0, -50.0, -40.0]), + Tensor.from_voigt([10.0, 20.0, 30.0, 60.0, 50.0, 40.0]), ] ) assert np.allclose(calc_chunk.stresses, stresses) @@ -339,30 +339,14 @@ def test_calc_stresses(calc_chunk): def test_calc_stress(calc_chunk): - stress = EV_PER_A3_TO_KBAR * full_3x3_to_voigt_6_stress( - np.array( - [ - [1.00000000, 2.00000000, 3.00000000], - [2.00000000, 5.00000000, 6.00000000], - [3.00000000, 6.00000000, 7.00000000], - ] - ) - ) + stress = EV_PER_A3_TO_KBAR * np.array([[1.0, 2.0, 3.0], [2.0, 5.0, 6.0], [3.0, 6.0, 7.0]]) assert np.allclose(calc_chunk.stress, stress) assert np.allclose(calc_chunk.structure.properties["stress"], stress) assert np.allclose(calc_chunk.results["stress"], stress) def test_calc_num_stress(numerical_stress_chunk): - stress = EV_PER_A3_TO_KBAR * full_3x3_to_voigt_6_stress( - np.array( - [ - [1.00000000, 2.00000000, 3.00000000], - [2.00000000, 5.00000000, 6.00000000], - [3.00000000, 6.00000000, 7.00000000], - ] - ) - ) + stress = EV_PER_A3_TO_KBAR * np.array([[1.0, 2.0, 3.0], [2.0, 5.0, 6.0], [3.0, 6.0, 7.0]]) assert np.allclose(numerical_stress_chunk.stress, stress) assert np.allclose(numerical_stress_chunk.structure.properties["stress"], stress) assert np.allclose(numerical_stress_chunk.results["stress"], stress) @@ -470,9 +454,9 @@ def test_molecular_header_initial_structure(molecular_header_chunk, molecular_po molecular_header_chunk.initial_structure.cart_coords, np.array( [ - [0.00000000, 0.00000000, 0.00000000], - [0.95840000, 0.00000000, 0.00000000], - [-0.24000000, 0.92790000, 0.00000000], + [0, 0, 0], + [0.95840000, 0, 0], + [-0.24000000, 0.92790000, 0], ] ), ) @@ -490,13 +474,7 @@ def molecular_calc_chunk(molecular_header_chunk): @pytest.fixture def molecular_positions(): - return np.array( - [ - [-0.00191785, -0.00243279, 0.00000000], - [0.97071531, -0.00756333, 0.00000000], - [-0.25039746, 0.93789612, 0.00000000], - ] - ) + return np.array([[-0.00191785, -0.00243279, 0], [0.97071531, -0.00756333, 0], [-0.25039746, 0.93789612, 0]]) def test_molecular_calc_atoms(molecular_calc_chunk, molecular_positions): @@ -572,11 +550,7 @@ def test_molecular_calc_hirshfeld_volumes(molecular_calc_chunk): def test_molecular_calc_hirshfeld_atomic_dipoles(molecular_calc_chunk): hirshfeld_atomic_dipoles = np.array( - [ - [0.04249319, 0.05486053, 0.00000000], - [0.13710134, -0.00105126, 0.00000000], - [-0.03534982, 0.13248706, 0.00000000], - ] + [[0.04249319, 0.05486053, 0], [0.13710134, -0.00105126, 0], [-0.03534982, 0.13248706, 0]] ) assert np.allclose(molecular_calc_chunk.hirshfeld_atomic_dipoles, hirshfeld_atomic_dipoles) assert np.allclose(