From ab08e91634ff2bce6d7239d444991eab04ea635f Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 21 Nov 2024 01:05:21 +0100 Subject: [PATCH 1/4] Fix traj parsing --- ipi/utils/parsing.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/ipi/utils/parsing.py b/ipi/utils/parsing.py index ef9455bb1..1313323e8 100644 --- a/ipi/utils/parsing.py +++ b/ipi/utils/parsing.py @@ -139,7 +139,7 @@ def read_trajectory( file_handle = open(filename, "r") bohr2angstrom = unit_to_user("length", "angstrom", 1.0) - comment_regex = re.compile(r"([^)]+)\{([^}]+)\}") + comment_regex = re.compile(r"(\w+)\{([^}]+)\}") step_regex = re.compile(r"Step:\s+(\d+)") frames = [] @@ -207,10 +207,10 @@ def read_trajectory( # parse comment to get the property matches = comment_regex.findall(ret["comment"]) - + print("MATCHES ", matches) # get what we have found - if len(matches) >= 2: - what = matches[-2][0] + if len(matches) == 2: + what = matches[0][0] else: # defaults to reading positions what = "positions" @@ -219,13 +219,14 @@ def read_trajectory( if len(matches) >= 1: frame.info["step"] = int(matches[-1][0]) - # if we have forces, set positions to zero (that data is missing!) and set forces instead + # if we have forces, set positions to zero (that data is missing!) + # and set forces instead if what == "forces": # set forces and convert to eV/angstrom frame.positions *= 0 frame.arrays["forces"] = ret["atoms"].q.reshape( (-1, 3) - ) * unit_to_user("force", "ev/ang", 1.0) + ) * unit_to_user("force", "ase", 1.0) frames.append(frame) From 16dadfbb1619606e4fd6ca3d6f194a655886d1cc Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 21 Nov 2024 01:07:13 +0100 Subject: [PATCH 2/4] Remove debug print --- ipi/utils/parsing.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ipi/utils/parsing.py b/ipi/utils/parsing.py index 1313323e8..225713604 100644 --- a/ipi/utils/parsing.py +++ b/ipi/utils/parsing.py @@ -207,7 +207,6 @@ def read_trajectory( # parse comment to get the property matches = comment_regex.findall(ret["comment"]) - print("MATCHES ", matches) # get what we have found if len(matches) == 2: what = matches[0][0] From fda7aec239de1434d94878d1507c189b70e20606 Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Thu, 28 Nov 2024 13:27:10 +0100 Subject: [PATCH 3/4] Improved handling of other trajectory types --- ipi/utils/parsing.py | 36 ++++++++++++++++++++++++++++++------ 1 file changed, 30 insertions(+), 6 deletions(-) diff --git a/ipi/utils/parsing.py b/ipi/utils/parsing.py index 225713604..2b1bd4f9e 100644 --- a/ipi/utils/parsing.py +++ b/ipi/utils/parsing.py @@ -9,6 +9,7 @@ import re import numpy as np from ipi.utils.units import unit_to_user +from ipi.utils.messages import warning, verbosity try: import ase @@ -220,12 +221,35 @@ def read_trajectory( # if we have forces, set positions to zero (that data is missing!) # and set forces instead - if what == "forces": - # set forces and convert to eV/angstrom - frame.positions *= 0 - frame.arrays["forces"] = ret["atoms"].q.reshape( - (-1, 3) - ) * unit_to_user("force", "ase", 1.0) + if what in [ + "positions", + "x_centroid", + "x_centroid_even", + "x_centroid_odd", + ]: + pass # nothing to do here, positions is the right place to store + else: + frame.positions *= 0.0 # trajectory does NOT contain position data! + frame.arrays[what] = ret["atoms"].q.reshape((-1, 3)) + if what in [ + "forces", + "forces_spring", + "forces_component", + "forces_sc", + "Eforces", + "f_centroid", + "forces_component_raw", + ]: + frame.arrays[what] *= unit_to_user("force", "ase", 1.0) + elif what in ["velocities", "v_centroid"]: + frame.arrays[what] *= unit_to_user("velocity", "ase", 1.0) + elif what in ["momenta", "p_centroid"]: + frame.arrays[what] *= unit_to_user("momentum", "ase", 1.0) + else: + warning( + f"No units conversion implemented for trajectory type {what}", + verbosity.low, + ) frames.append(frame) From 3d89bacef8c7fa1c4247dde2d0a2dfa7d978f02f Mon Sep 17 00:00:00 2001 From: Michele Ceriotti Date: Fri, 29 Nov 2024 00:38:14 +0100 Subject: [PATCH 4/4] Fixed units conversion in parser in a smarter way --- ipi/utils/parsing.py | 54 +++++++++++++++++--------------------------- ipi/utils/units.py | 6 ++--- 2 files changed, 24 insertions(+), 36 deletions(-) diff --git a/ipi/utils/parsing.py b/ipi/utils/parsing.py index 2b1bd4f9e..ef4ea1a42 100644 --- a/ipi/utils/parsing.py +++ b/ipi/utils/parsing.py @@ -139,7 +139,6 @@ def read_trajectory( raise ValueError(f"Unrecognized file format: {format}") file_handle = open(filename, "r") - bohr2angstrom = unit_to_user("length", "angstrom", 1.0) comment_regex = re.compile(r"(\w+)\{([^}]+)\}") step_regex = re.compile(r"Step:\s+(\d+)") @@ -201,8 +200,9 @@ def read_trajectory( frame = ase.Atoms( ret["atoms"].names, - positions=ret["atoms"].q.reshape((-1, 3)) * bohr2angstrom, - cell=ret["cell"].h.T * bohr2angstrom, + # will apply conversion later! + positions=ret["atoms"].q.reshape((-1, 3)), + cell=ret["cell"].h.T * unit_to_user("length", "ase"), pbc=True, ) @@ -219,37 +219,25 @@ def read_trajectory( if len(matches) >= 1: frame.info["step"] = int(matches[-1][0]) - # if we have forces, set positions to zero (that data is missing!) - # and set forces instead - if what in [ - "positions", - "x_centroid", - "x_centroid_even", - "x_centroid_odd", - ]: - pass # nothing to do here, positions is the right place to store + # fetch the list of known traj types, cf. `engine/properties.py`` + from ipi.engine.properties import Trajectories # avoids circular import + + traj_types = Trajectories().traj_dict + if not what in traj_types: + warning( + f"{what} is not a known trajectory type. Will apply no units conversion", + verbosity.low, + ) + elif traj_types[what]["dimension"] == "length": + # positions is the right place to store, and we just need to convert + frame.positions *= unit_to_user("length", "ase") else: - frame.positions *= 0.0 # trajectory does NOT contain position data! - frame.arrays[what] = ret["atoms"].q.reshape((-1, 3)) - if what in [ - "forces", - "forces_spring", - "forces_component", - "forces_sc", - "Eforces", - "f_centroid", - "forces_component_raw", - ]: - frame.arrays[what] *= unit_to_user("force", "ase", 1.0) - elif what in ["velocities", "v_centroid"]: - frame.arrays[what] *= unit_to_user("velocity", "ase", 1.0) - elif what in ["momenta", "p_centroid"]: - frame.arrays[what] *= unit_to_user("momentum", "ase", 1.0) - else: - warning( - f"No units conversion implemented for trajectory type {what}", - verbosity.low, - ) + # if we have another type of value, set positions to zero + # (that data is missing!) and set an array instead + frame.positions *= 0.0 + frame.arrays[what] = ret["atoms"].q.reshape((-1, 3)) * unit_to_user( + traj_types[what]["dimension"], "ase" + ) frames.append(frame) diff --git a/ipi/utils/units.py b/ipi/utils/units.py index 4841a1950..6f621bbf2 100644 --- a/ipi/utils/units.py +++ b/ipi/utils/units.py @@ -365,7 +365,7 @@ def mass(cls, label): # -def unit_to_internal(family, unit, number): +def unit_to_internal(family, unit, number=1.0): """Converts a number of given dimensions and units into internal units. Args: @@ -411,7 +411,7 @@ def unit_to_internal(family, unit, number): return number * UnitMap[family][base.lower()] * UnitPrefix[prefix] -def unit_to_user(family, unit, number): +def unit_to_user(family, unit, number=1.0): """Converts a number of given dimensions from internal to user units. Args: @@ -423,4 +423,4 @@ def unit_to_user(family, unit, number): The number in the user specified units """ - return number / unit_to_internal(family, unit, 1.0) + return number / unit_to_internal(family, unit)