Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
161 changes: 161 additions & 0 deletions dpdata/formats/xyz/_unit_convert.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
"""Unit conversion helpers for extended XYZ (extxyz) format.

This module provides a table-driven approach to convert energy, force, and
stress/pressure values from various unit systems commonly found in extxyz
files into dpdata's internal units:

- Energy: eV
- Force: eV/angstrom
- Stress/Pressure: eV/angstrom^3 (before volume multiplication to get virial)
"""

from __future__ import annotations

from dpdata.unit import EnergyConversion, ForceConversion, PressureConversion

# ---------------------------------------------------------------------------
# Unit alias mapping tables
# Keys are LOWERCASE; values are canonical names recognized by dpdata.unit
# ---------------------------------------------------------------------------

_ENERGY_UNIT_MAP: dict[str, str] = {

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please keep common atomic-unit aliases in this table as well. au / a.u. are frequent spellings for Hartree in chemistry datasets, and this helper currently rejects energy-unit=au even though the PR is meant to handle unit synonyms. If force-unit=au / a.u. is intended to mean Hartree/Bohr, that should be handled explicitly too, with regression tests.

"ev": "eV",
"hartree": "hartree",
"ha": "hartree",
"au": "hartree",
"a.u.": "hartree",
"ry": "rydberg",
"rydberg": "rydberg",
"kcal/mol": "kcal_mol",
"kcal_mol": "kcal_mol",
"kj/mol": "kJ_mol",
"kj_mol": "kJ_mol",
}

_LENGTH_UNIT_MAP: dict[str, str] = {
"angstrom": "angstrom",
"ang": "angstrom",
"ang.": "angstrom",
"a": "angstrom",
"bohr": "bohr",
"nm": "nm",
}

_PRESSURE_UNIT_MAP: dict[str, str] = {
"gpa": "GPa",
"kpa": "kPa",
"pa": "Pa",
"kbar": "kbar",
"bar": "bar",
"ev/angstrom^3": "eV/angstrom^3",
"ev/ang^3": "eV/angstrom^3",
"ev/a^3": "eV/angstrom^3",
"ha/bohr^3": "hartree/bohr^3",
"hartree/bohr^3": "hartree/bohr^3",
"au/bohr^3": "hartree/bohr^3",
"a.u./bohr^3": "hartree/bohr^3",
}

# dpdata internal unit strings
_INTERNAL_ENERGY = "eV"
_INTERNAL_FORCE = "eV/angstrom"
_INTERNAL_PRESSURE = "eV/angstrom^3"


def _parse_force_unit(raw: str) -> tuple[str, str]:
"""Split a composite force unit string into (energy_part, length_part).

Examples
--------
>>> _parse_force_unit("kcal/mol/angstrom")
('kcal/mol', 'angstrom')
>>> _parse_force_unit("hartree/bohr")
('hartree', 'bohr')
>>> _parse_force_unit("ev/ang")
('ev', 'ang')
"""
# Handle atomic-unit shorthand: "au" or "a.u." means hartree/bohr for force
_FORCE_UNIT_ALIASES = {"au": ("hartree", "bohr"), "a.u.": ("hartree", "bohr")}
if raw in _FORCE_UNIT_ALIASES:
return _FORCE_UNIT_ALIASES[raw]

# Try matching known energy prefixes (longest first) to handle
# composite names like "kcal/mol" that themselves contain "/".
for e_key in sorted(_ENERGY_UNIT_MAP.keys(), key=len, reverse=True):
prefix = e_key + "/"
if raw.startswith(prefix):
l_part = raw[len(prefix) :]
if l_part:
return e_key, l_part
# Fallback: split on last "/"
parts = raw.rsplit("/", 1)
if len(parts) == 2 and parts[0] and parts[1]:
return parts[0], parts[1]
raise ValueError(f"Cannot parse force unit string: '{raw}'")


def _get_unit_factor(unit_str: str | None, quantity: str) -> float:
"""Return the multiplicative factor to convert from the given unit to dpdata internals.

Parameters
----------
unit_str : str or None
The unit string read from the extxyz header (e.g. "hartree", "kcal/mol/angstrom").
If None, returns 1.0 (assumes data is already in internal units).
quantity : str
One of "energy", "force", or "stress".

Returns
-------
float
Conversion factor such that ``value_internal = value_file * factor``.

Raises
------
ValueError
If the unit string is not recognized or the quantity type is invalid.
"""
if unit_str is None:
return 1.0

key = unit_str.lower().strip()

if quantity == "energy":
canonical = _ENERGY_UNIT_MAP.get(key)
if canonical is None:
raise ValueError(
f"Unsupported energy unit: '{unit_str}'. "
f"Supported: {list(_ENERGY_UNIT_MAP.keys())}"
)
return EnergyConversion(canonical, _INTERNAL_ENERGY).value()

elif quantity == "force":
e_part, l_part = _parse_force_unit(key)
e_canonical = _ENERGY_UNIT_MAP.get(e_part)
l_canonical = _LENGTH_UNIT_MAP.get(l_part)
if e_canonical is None:
raise ValueError(
f"Unsupported energy part in force unit: '{e_part}' "
f"(from '{unit_str}'). Supported: {list(_ENERGY_UNIT_MAP.keys())}"
)
if l_canonical is None:
raise ValueError(
f"Unsupported length part in force unit: '{l_part}' "
f"(from '{unit_str}'). Supported: {list(_LENGTH_UNIT_MAP.keys())}"
)
src_unit = f"{e_canonical}/{l_canonical}"
return ForceConversion(src_unit, _INTERNAL_FORCE).value()

elif quantity == "stress":
canonical = _PRESSURE_UNIT_MAP.get(key)
if canonical is None:
raise ValueError(
f"Unsupported stress/pressure unit: '{unit_str}'. "
f"Supported: {list(_PRESSURE_UNIT_MAP.keys())}"
)
return PressureConversion(canonical, _INTERNAL_PRESSURE).value()

else:
raise ValueError(
f"Unknown quantity type: '{quantity}'. Must be 'energy', 'force', or 'stress'."
)
55 changes: 47 additions & 8 deletions dpdata/formats/xyz/quip_gap_xyz.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,14 @@
# %%
from __future__ import annotations

import logging
import re
import warnings
from collections import OrderedDict

import numpy as np

from dpdata.formats.xyz._unit_convert import _get_unit_factor
from dpdata.periodic_table import Element

# Possible keys for the energy field in the extxyz comment line,
Expand All @@ -24,7 +26,7 @@
_STRESS_KEYS = ("stress", "stresses")


def _parse_stress_to_virials(stress_str, cell, stress_sign=-1):
def _parse_stress_to_virials(stress_str, cell, stress_sign=-1, stress_factor=1.0):
"""Convert a stress field string to virial tensor.

Parameters
Expand All @@ -38,6 +40,10 @@ def _parse_stress_to_virials(stress_str, cell, stress_sign=-1):
Sign convention for ``virial = stress_sign * volume * stress``.
Default ``-1`` follows the ASE convention where
``virial = -V * stress`` (stress in eV/angstrom^3).
stress_factor : float
Multiplicative factor to convert stress from file units to
eV/angstrom^3 before computing the virial. Default 1.0
(assumes stress is already in eV/angstrom^3).

Returns
-------
Expand All @@ -57,7 +63,7 @@ def _parse_stress_to_virials(stress_str, cell, stress_sign=-1):
f"stress field must have 6 (Voigt) or 9 (3x3) values, got {len(vals)}"
)
volume = abs(np.linalg.det(cell))
virials = stress_sign * volume * stress
virials = stress_sign * volume * stress * stress_factor
return np.array([virials])


Expand Down Expand Up @@ -249,13 +255,46 @@ def handle_single_xyz_frame(lines, stress_sign=-1, **kwargs):
stress_raw = field_dict[skey]
break

# --- unit conversion factors ---
# Read optional unit metadata from extxyz header.
# When absent, factor is 1.0 (assumes dpdata internal units).
e_unit_str = field_dict.get("energy-unit") or field_dict.get("energy_unit")
f_unit_str = field_dict.get("force-unit") or field_dict.get("force_unit")
s_unit_str = field_dict.get("stress-unit") or field_dict.get("stress_unit")

e_factor = _get_unit_factor(e_unit_str, "energy")
f_factor = _get_unit_factor(f_unit_str, "force")
s_factor = _get_unit_factor(s_unit_str, "stress")

if e_factor != 1.0:
logging.info(
"extxyz: converting energy from '%s' to eV (factor=%.10g)",
e_unit_str,
e_factor,
)
if f_factor != 1.0:
logging.info(
"extxyz: converting forces from '%s' to eV/angstrom (factor=%.10g)",
f_unit_str,
f_factor,
)
if s_factor != 1.0:
logging.info(
"extxyz: converting stress from '%s' to eV/angstrom^3 (factor=%.10g)",
s_unit_str,
s_factor,
)

if virial_raw is not None:
virials = np.array(
[np.array(list(filter(bool, virial_raw.split(" ")))).reshape(3, 3)]
).astype(np.float64)
virials = (
np.array(
[np.array(list(filter(bool, virial_raw.split(" ")))).reshape(3, 3)]
).astype(np.float64)
* e_factor # virial has energy dimensions; convert if energy-unit is set
)
elif stress_raw is not None:
virials = _parse_stress_to_virials(
stress_raw, cells, stress_sign=stress_sign
stress_raw, cells, stress_sign=stress_sign, stress_factor=s_factor
)

# --- energy (try several common keys) ---
Expand All @@ -275,8 +314,8 @@ def handle_single_xyz_frame(lines, stress_sign=-1, **kwargs):
info_dict["atom_numbs"] = list(type_num_array[:, 1].astype(int))
info_dict["atom_types"] = np.array(atom_type_list).astype(int)
info_dict["coords"] = np.array([coords_array]).astype(np.float64)
info_dict["energies"] = np.array([energy_value]).astype(np.float64)
info_dict["forces"] = np.array([force_array]).astype(np.float64)
info_dict["energies"] = np.array([energy_value], dtype=np.float64) * e_factor
info_dict["forces"] = np.array([force_array], dtype=np.float64) * f_factor
if virials is not None:
info_dict["virials"] = virials
info_dict["orig"] = np.zeros(3)
Expand Down
Loading
Loading