From fee82e3026f35e6e76df4e90ae429cbec001da6c Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Sat, 25 May 2024 17:52:20 +0200 Subject: [PATCH 1/8] QE out parser now also reads energy contributions --- mala/targets/target.py | 141 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 141 insertions(+) diff --git a/mala/targets/target.py b/mala/targets/target.py index 23212470b..ce0362f96 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -119,6 +119,12 @@ def __init__(self, params): self.band_energy_dft_calculation = None self.total_energy_dft_calculation = None self.entropy_contribution_dft_calculation = None + self.total_energy_contributions_dft_calculation = { + "one_electron_contribution": None, + "hartree_contribution": None, + "xc_contribution": None, + "ewald_contribution": None, + } self.atoms = None self.electrons_per_atom = None self.qe_input_data = { @@ -319,6 +325,12 @@ def read_additional_calculation_data(self, data, data_type=None): self.band_energy_dft_calculation = None self.total_energy_dft_calculation = None self.entropy_contribution_dft_calculation = None + self.total_energy_contributions_dft_calculation = { + "one_electron_contribution": None, + "hartree_contribution": None, + "xc_contribution": None, + "ewald_contribution": None, + } self.grid_dimensions = [0, 0, 0] self.atoms = None @@ -333,7 +345,14 @@ def read_additional_calculation_data(self, data, data_type=None): total_energy = None past_calculation_part = False bands_included = True + + # individual energy contributions. entropy_contribution = None + one_electron_contribution = None + hartree_contribution = None + xc_contribution = None + ewald_contribution = None + with open(data) as out: pseudolinefound = False lastpseudo = None @@ -390,6 +409,32 @@ def read_additional_calculation_data(self, data, data_type=None): entropy_contribution = float( (line.split("=")[1]).split("Ry")[0] ) + if ( + "one-electron contribution" in line + and past_calculation_part + ): + if one_electron_contribution is None: + one_electron_contribution = float( + (line.split("=")[1]).split("Ry")[0] + ) + if ( + "hartree contribution" in line + and past_calculation_part + ): + if hartree_contribution is None: + hartree_contribution = float( + (line.split("=")[1]).split("Ry")[0] + ) + if "xc contribution" in line and past_calculation_part: + if xc_contribution is None: + xc_contribution = float( + (line.split("=")[1]).split("Ry")[0] + ) + if "ewald contribution" in line and past_calculation_part: + if ewald_contribution is None: + ewald_contribution = float( + (line.split("=")[1]).split("Ry")[0] + ) if "set verbosity='high' to print them." in line: bands_included = False @@ -413,6 +458,25 @@ def read_additional_calculation_data(self, data, data_type=None): self.entropy_contribution_dft_calculation = ( entropy_contribution * Rydberg ) + if one_electron_contribution is not None: + self.total_energy_contributions_dft_calculation[ + "one_electron_contribution" + ] = (one_electron_contribution * Rydberg) + + if hartree_contribution is not None: + self.total_energy_contributions_dft_calculation[ + "hartree_contribution" + ] = (hartree_contribution * Rydberg) + + if xc_contribution is not None: + self.total_energy_contributions_dft_calculation[ + "xc_contribution" + ] = (xc_contribution * Rydberg) + + if ewald_contribution is not None: + self.total_energy_contributions_dft_calculation[ + "ewald_contribution" + ] = (ewald_contribution * Rydberg) # Calculate band energy, if the necessary data is included in # the output file. @@ -446,6 +510,12 @@ def read_additional_calculation_data(self, data, data_type=None): self.band_energy_dft_calculation = None self.total_energy_dft_calculation = None self.entropy_contribution_dft_calculation = None + self.total_energy_contributions_dft_calculation = { + "one_electron_contribution": None, + "hartree_contribution": None, + "xc_contribution": None, + "ewald_contribution": None, + } self.grid_dimensions = [0, 0, 0] self.atoms: ase.Atoms = data[0] @@ -490,6 +560,12 @@ def read_additional_calculation_data(self, data, data_type=None): self.voxel = None self.band_energy_dft_calculation = None self.total_energy_dft_calculation = None + self.total_energy_contributions_dft_calculation = { + "one_electron_contribution": None, + "hartree_contribution": None, + "xc_contribution": None, + "ewald_contribution": None, + } self.entropy_contribution_dft_calculation = None self.grid_dimensions = [0, 0, 0] self.atoms = None @@ -505,6 +581,21 @@ def read_additional_calculation_data(self, data, data_type=None): self.qe_input_data["degauss"] = json_dict["degauss"] self.qe_pseudopotentials = json_dict["pseudopotentials"] + # These attributes are only needed for debugging purposes. + # The interace should not break if they are not present in the + # json file. + energy_contribution_ids = [ + "one_electron_contribution", + "hartree_contribution", + "xc_contribution", + "ewald_contribution", + ] + for key in energy_contribution_ids: + if key in json_dict: + self.total_energy_contributions_dft_calculation[key] = ( + json_dict[key] + ) + else: raise Exception("Unsupported auxiliary file type.") @@ -524,6 +615,7 @@ def write_additional_calculation_data(self, filepath, return_string=False): If True, no file will be written, and instead a json dict will be returned. """ + additional_calculation_data = { "fermi_energy_dft": self.fermi_energy_dft, "temperature": self.temperature, @@ -539,6 +631,18 @@ def write_additional_calculation_data(self, filepath, return_string=False): "degauss": self.qe_input_data["degauss"], "pseudopotentials": self.qe_pseudopotentials, "entropy_contribution_dft_calculation": self.entropy_contribution_dft_calculation, + "one_electron_contribution": self.total_energy_contributions_dft_calculation[ + "one_electron_contribution" + ], + "hartree_contribution": self.total_energy_contributions_dft_calculation[ + "hartree_contribution" + ], + "xc_contribution": self.total_energy_contributions_dft_calculation[ + "xc_contribution" + ], + "ewald_contribution": self.total_energy_contributions_dft_calculation[ + "ewald_contribution" + ], } if self.voxel is not None: additional_calculation_data["voxel"] = self.voxel.todict() @@ -1447,6 +1551,43 @@ def _process_openpmd_attributes(self, series, iteration, mesh): iteration.get_attribute(attribute) ) + self.total_energy_contributions_dft_calculation[ + "one_electron_contribution" + ] = self._get_attribute_if_attribute_exists( + iteration, + "one_electron_contribution", + default_value=self.total_energy_contributions_dft_calculation[ + "one_electron_contribution" + ], + ) + self.total_energy_contributions_dft_calculation[ + "hartree_contribution" + ] = self._get_attribute_if_attribute_exists( + iteration, + "hartree_contribution", + default_value=self.total_energy_contributions_dft_calculation[ + "hartree_contribution" + ], + ) + self.total_energy_contributions_dft_calculation["xc_contribution"] = ( + self._get_attribute_if_attribute_exists( + iteration, + "xc_contribution", + default_value=self.total_energy_contributions_dft_calculation[ + "xc_contribution" + ], + ) + ) + self.total_energy_contributions_dft_calculation[ + "ewald_contribution" + ] = self._get_attribute_if_attribute_exists( + iteration, + "ewald_contribution", + default_value=self.total_energy_contributions_dft_calculation[ + "ewald_contribution" + ], + ) + def _set_geometry_info(self, mesh): # Geometry: Save the cell parameters and angles of the grid. if self.atoms is not None: From 789f9096ca5434369e2a6a1f354a858ac9d05100 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Tue, 11 Jun 2024 15:06:22 +0200 Subject: [PATCH 2/8] Also added the forces for good measure --- mala/targets/target.py | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index ce0362f96..31d962508 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -125,6 +125,7 @@ def __init__(self, params): "xc_contribution": None, "ewald_contribution": None, } + self.atomic_forces_dft = None self.atoms = None self.electrons_per_atom = None self.qe_input_data = { @@ -331,6 +332,7 @@ def read_additional_calculation_data(self, data, data_type=None): "xc_contribution": None, "ewald_contribution": None, } + self.atomic_forces_dft = None self.grid_dimensions = [0, 0, 0] self.atoms = None @@ -340,6 +342,14 @@ def read_additional_calculation_data(self, data, data_type=None): self.fermi_energy_dft = ( self.atoms.get_calculator().get_fermi_level() ) + # The forces may not have been computed. If they are indeed + # not computed, ASE will by default throw an PropertyNotImplementedError + # error + try: + self.atomic_forces_dft = self.atoms.get_forces() + except ase.calculators.calculator.PropertyNotImplementedError: + print("CAUGHT AN ERROR!") + pass # Parse the file for energy values. total_energy = None @@ -516,6 +526,7 @@ def read_additional_calculation_data(self, data, data_type=None): "xc_contribution": None, "ewald_contribution": None, } + self.atomic_forces_dft = None self.grid_dimensions = [0, 0, 0] self.atoms: ase.Atoms = data[0] @@ -566,6 +577,7 @@ def read_additional_calculation_data(self, data, data_type=None): "xc_contribution": None, "ewald_contribution": None, } + self.atomic_forces_dft = None self.entropy_contribution_dft_calculation = None self.grid_dimensions = [0, 0, 0] self.atoms = None @@ -581,9 +593,6 @@ def read_additional_calculation_data(self, data, data_type=None): self.qe_input_data["degauss"] = json_dict["degauss"] self.qe_pseudopotentials = json_dict["pseudopotentials"] - # These attributes are only needed for debugging purposes. - # The interace should not break if they are not present in the - # json file. energy_contribution_ids = [ "one_electron_contribution", "hartree_contribution", @@ -596,6 +605,12 @@ def read_additional_calculation_data(self, data, data_type=None): json_dict[key] ) + # Not always read from DFT files. + if "atomic_forces_dft" in json_dict: + self.atomic_forces_dft = np.array( + json_dict["atomic_forces_dft"] + ) + else: raise Exception("Unsupported auxiliary file type.") @@ -664,6 +679,11 @@ def write_additional_calculation_data(self, filepath, return_string=False): additional_calculation_data["atoms"]["pbc"] = ( additional_calculation_data["atoms"]["pbc"].tolist() ) + if self.atomic_forces_dft is not None: + additional_calculation_data["atomic_forces_dft"] = ( + self.atomic_forces_dft.tolist() + ) + if return_string is False: with open(filepath, "w", encoding="utf-8") as f: json.dump( From ba660fe5d803e9ec66239f84be64bbbed4c68586 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Wed, 17 Jul 2024 17:38:33 +0200 Subject: [PATCH 3/8] Got rid of superfluous print statement --- mala/targets/target.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index 31d962508..79c5222b5 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -348,7 +348,6 @@ def read_additional_calculation_data(self, data, data_type=None): try: self.atomic_forces_dft = self.atoms.get_forces() except ase.calculators.calculator.PropertyNotImplementedError: - print("CAUGHT AN ERROR!") pass # Parse the file for energy values. From e66879983e6735c0bc6d4abe67dd4774c10fb603 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Wed, 23 Oct 2024 16:30:47 +0200 Subject: [PATCH 4/8] Fixed docstring --- mala/targets/target.py | 1 - 1 file changed, 1 deletion(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index 79c5222b5..4c7782849 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -629,7 +629,6 @@ def write_additional_calculation_data(self, filepath, return_string=False): If True, no file will be written, and instead a json dict will be returned. """ - additional_calculation_data = { "fermi_energy_dft": self.fermi_energy_dft, "temperature": self.temperature, From 2de6782e808b2c7e1815bee6366e5f9d7ac393af Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Wed, 23 Oct 2024 17:38:18 +0200 Subject: [PATCH 5/8] Working on a fix for the OpenPMD interface --- mala/targets/target.py | 42 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index 4c7782849..f0272cb58 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -14,7 +14,7 @@ from scipy.integrate import simps from mala.common.parameters import Parameters, ParametersTargets -from mala.common.parallelizer import printout, parallel_warn +from mala.common.parallelizer import printout, parallel_warn, get_rank from mala.targets.calculation_helpers import fermi_function from mala.common.physical_data import PhysicalData from mala.descriptors.atomic_density import AtomicDensity @@ -1428,6 +1428,8 @@ def _process_additional_metadata(self, additional_metadata): ) def _set_openpmd_attribtues(self, iteration, mesh): + import openpmd_api as io + super(Target, self)._set_openpmd_attribtues(iteration, mesh) # If no atoms have been read, neither have any of the other @@ -1443,6 +1445,7 @@ def _set_openpmd_attribtues(self, iteration, mesh): and key is not None and key != "pseudopotentials" and additional_calculation_data[key] is not None + and key != "atomic_forces_dft" ): iteration.set_attribute(key, additional_calculation_data[key]) if key == "pseudopotentials": @@ -1456,6 +1459,43 @@ def _set_openpmd_attribtues(self, iteration, mesh): ], ) + # If the data contains atomic forces from a DFT calculation, we need + # to process it in much the same fashion as the atoms. + if "atomic_forces_dft" in additional_calculation_data: + atomic_forces = additional_calculation_data["atomic_forces_dft"] + atomic_forces_2 = self.atomic_forces_dft + if atomic_forces is not None: + # This data is equivalent across the ranks, so just write it once + atomic_forces_dft_openpmd = iteration.particles[ + "atomic_forces_dft" + ] + forces = io.Dataset( + # Need bugfix https://github.com/openPMD/openPMD-api/pull/1357 + ( + np.array(atomic_forces[0]).dtype + if io.__version__ >= "0.15.0" + else io.Datatype.DOUBLE + ), + np.array(atomic_forces[0]).shape, + ) + # atoms_openpmd["position"].time_offset = 0.0 + # atoms_openpmd["positionOffset"].time_offset = 0.0 + for atom in range(0, np.shape(atomic_forces_dft_openpmd)[0]): + atomic_forces_dft_openpmd["atomic_forces_dft"][ + str(atom) + ].reset_dataset(forces) + + individual_force = atomic_forces_dft_openpmd[ + "atomic_forces_dft" + ][str(atom)] + if get_rank() == 0: + individual_force.store_chunk(atomic_forces[atom]) + + # Positions are stored in Angstrom. + atomic_forces_dft_openpmd["position"][ + str(atom) + ].unit_SI = 1.0e-10 + def _process_openpmd_attributes(self, series, iteration, mesh): super(Target, self)._process_openpmd_attributes( series, iteration, mesh From 6d6df8a2c1bac98595a985eb3548eca1e111beba Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Thu, 24 Oct 2024 09:34:17 +0200 Subject: [PATCH 6/8] Writing seems to work, reading not yet --- mala/targets/target.py | 37 +++++++++++++++++++++++++++++-------- 1 file changed, 29 insertions(+), 8 deletions(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index f0272cb58..0c4d6b884 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -1463,7 +1463,6 @@ def _set_openpmd_attribtues(self, iteration, mesh): # to process it in much the same fashion as the atoms. if "atomic_forces_dft" in additional_calculation_data: atomic_forces = additional_calculation_data["atomic_forces_dft"] - atomic_forces_2 = self.atomic_forces_dft if atomic_forces is not None: # This data is equivalent across the ranks, so just write it once atomic_forces_dft_openpmd = iteration.particles[ @@ -1478,21 +1477,21 @@ def _set_openpmd_attribtues(self, iteration, mesh): ), np.array(atomic_forces[0]).shape, ) - # atoms_openpmd["position"].time_offset = 0.0 - # atoms_openpmd["positionOffset"].time_offset = 0.0 - for atom in range(0, np.shape(atomic_forces_dft_openpmd)[0]): - atomic_forces_dft_openpmd["atomic_forces_dft"][ + for atom in range(0, np.shape(atomic_forces)[0]): + atomic_forces_dft_openpmd["force_compopnents"][ str(atom) ].reset_dataset(forces) individual_force = atomic_forces_dft_openpmd[ - "atomic_forces_dft" + "force_compopnents" ][str(atom)] if get_rank() == 0: - individual_force.store_chunk(atomic_forces[atom]) + individual_force.store_chunk( + np.array(atomic_forces)[atom] + ) # Positions are stored in Angstrom. - atomic_forces_dft_openpmd["position"][ + atomic_forces_dft_openpmd["force_compopnents"][ str(atom) ].unit_SI = 1.0e-10 @@ -1536,6 +1535,28 @@ def _process_openpmd_attributes(self, series, iteration, mesh): "periodic_boundary_conditions_z" ) + # Forces may not necessarily have been read (and therefore written) + atomic_forces_dft = iteration.particles["atomic_forces_dft"] + nr_atoms = len(atomic_forces_dft["force_compopnents"]) + self.atomic_forces_dft = np.zeros((nr_atoms, 3)) + for i in range(0, nr_atoms): + self.atomic_forces_dft["force_compopnents"][str(i)].load_chunk( + self.atomic_forces_dft[i, :] + ) + series.flush() + + # try: + # atomic_forces_dft = iteration.particles["atomic_forces_dft"] + # nr_atoms = len(atomic_forces_dft["atomic_forces_dft"]) + # self.atomic_forces_dft = np.zeros((nr_atoms, 3)) + # for i in range(0, nr_atoms): + # self.atomic_forces_dft["atomic_forces_dft"][str(i)].load_chunk( + # self.atomic_forces_dft[i, :] + # ) + # series.flush() + # except IndexError: + # pass + # Process all the regular meta info. self.fermi_energy_dft = self._get_attribute_if_attribute_exists( iteration, "fermi_energy_dft", default_value=self.fermi_energy_dft From d9eb1048ac6a74d3247db4f89d0390cbcf8343c4 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Thu, 24 Oct 2024 10:41:00 +0200 Subject: [PATCH 7/8] Fixed OpenPMD interface --- mala/targets/target.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index 0c4d6b884..f775a9ad1 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -1540,7 +1540,7 @@ def _process_openpmd_attributes(self, series, iteration, mesh): nr_atoms = len(atomic_forces_dft["force_compopnents"]) self.atomic_forces_dft = np.zeros((nr_atoms, 3)) for i in range(0, nr_atoms): - self.atomic_forces_dft["force_compopnents"][str(i)].load_chunk( + atomic_forces_dft["force_compopnents"][str(i)].load_chunk( self.atomic_forces_dft[i, :] ) series.flush() From e5a51fa99e9f3c969445809ba345b28aa98067b4 Mon Sep 17 00:00:00 2001 From: Lenz Fiedler Date: Thu, 24 Oct 2024 10:44:21 +0200 Subject: [PATCH 8/8] Now actually fixed the interface --- mala/targets/target.py | 32 ++++++++++++-------------------- 1 file changed, 12 insertions(+), 20 deletions(-) diff --git a/mala/targets/target.py b/mala/targets/target.py index f775a9ad1..dac17dcf3 100644 --- a/mala/targets/target.py +++ b/mala/targets/target.py @@ -1536,26 +1536,18 @@ def _process_openpmd_attributes(self, series, iteration, mesh): ) # Forces may not necessarily have been read (and therefore written) - atomic_forces_dft = iteration.particles["atomic_forces_dft"] - nr_atoms = len(atomic_forces_dft["force_compopnents"]) - self.atomic_forces_dft = np.zeros((nr_atoms, 3)) - for i in range(0, nr_atoms): - atomic_forces_dft["force_compopnents"][str(i)].load_chunk( - self.atomic_forces_dft[i, :] - ) - series.flush() - - # try: - # atomic_forces_dft = iteration.particles["atomic_forces_dft"] - # nr_atoms = len(atomic_forces_dft["atomic_forces_dft"]) - # self.atomic_forces_dft = np.zeros((nr_atoms, 3)) - # for i in range(0, nr_atoms): - # self.atomic_forces_dft["atomic_forces_dft"][str(i)].load_chunk( - # self.atomic_forces_dft[i, :] - # ) - # series.flush() - # except IndexError: - # pass + + try: + atomic_forces_dft = iteration.particles["atomic_forces_dft"] + nr_atoms = len(atomic_forces_dft["force_compopnents"]) + self.atomic_forces_dft = np.zeros((nr_atoms, 3)) + for i in range(0, nr_atoms): + atomic_forces_dft["force_compopnents"][str(i)].load_chunk( + self.atomic_forces_dft[i, :] + ) + series.flush() + except IndexError: + pass # Process all the regular meta info. self.fermi_energy_dft = self._get_attribute_if_attribute_exists(