Skip to content

Commit

Permalink
Convert all FHI-aims stresses to be 3x3 instead of Voigt notation (#3476
Browse files Browse the repository at this point in the history 
)

* 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 <tpurcell90@users.noreply.github.com>

* 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 <tpurcell90@users.noreply.github.com>
Co-authored-by: Janosh Riebesell <janosh.riebesell@gmail.com>
Co-authored-by: Thomas Purcell <purcellt@arizona.edu>
  • Loading branch information
@tpurcell90 @janosh
3 people authored Nov 21, 2023
1 parent c8154cd commit ec750ca
Show file tree
Hide file tree
Showing 2 changed files with 28 additions and 55 deletions.
21 changes: 10 additions & 11 deletions pymatgen/io/aims/parsers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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"],
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.")
Expand All @@ -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.")
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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:
Expand All @@ -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",
Expand All @@ -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:
Expand Down
62 changes: 18 additions & 44 deletions tests/io/aims/test_aims_parsers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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],
]
)

Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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],
]
),
)
Expand All @@ -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):
Expand Down Expand Up @@ -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(
Expand Down

0 comments on commit ec750ca

Please sign in to comment.