Merge branch 'protos' into chi/documentation-cleanup

MDANSE/Src/MDANSE/Framework/Configurators/OptionalXYZFileConfigurator.py

@@ -0,0 +1,35 @@
+#    This file is part of MDANSE.
+#
+#    MDANSE is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    This program is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+#
+from .XYZFileConfigurator import XYZFileConfigurator
+
+
+class OptionalXYZFileConfigurator(XYZFileConfigurator):
+    def configure(self, filepath: str) -> None:
+        """Configure the XYZ file if the filepath is not empty.
+
+        Parameters
+        ----------
+        filepath : str
+            THe filepath of the xyz file.
+        """
+        if not filepath:
+            self._original_input = filepath
+            self["value"] = filepath
+            self["filename"] = filepath
+            self.error_status = "OK"
+            return
+
+        super().configure(filepath)

MDANSE/Src/MDANSE/Framework/Configurators/WeightsConfigurator.py

@@ -13,6 +13,7 @@
 #    You should have received a copy of the GNU General Public License
 #    along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #
+from typing import Optional
 from collections import defaultdict
 import itertools
 
@@ -109,7 +110,23 @@ def configure(self, value):
         self["property"] = value
         self.error_status = "OK"
 
-    def get_weights(self):
+    def get_weights(self, prop: Optional[str] = None):
+        """Generate a dictionary of weights.
+
+        Parameters
+        ----------
+        prop : str or None, optional
+            The property to generate the weights from, if None then the
+            property set in this configurator will be used.
+
+        Returns
+        -------
+        dict[str, float]
+            The dictionary of the weights.
+        """
+        if not prop:
+            prop = self["property"]
+
         atom_selection_configurator = self._configurable[
             self._dependencies["atom_selection"]
         ]
@@ -123,7 +140,7 @@ def get_weights(self):
             atom_selection_configurator["selection_length"],
         ):
             weights[name] += sum(
-                self._trajectory.get_atom_property(element, self["property"])
+                self._trajectory.get_atom_property(element, prop)
                 for element in elements
             )
 

MDANSE/Src/MDANSE/Framework/Converters/CP2K.py

@@ -122,14 +122,21 @@ class CP2K(Converter):
         },
     )
     settings["vel_file"] = (
-        "XYZFileConfigurator",
+        "OptionalXYZFileConfigurator",
         {
             "wildcard": "XYZ files (*.xyz);;All files (*)",
-            "default": "INPUT_FILENAME.xyz",
-            "optional": True,
+            "default": "",
             "label": "Velocity file (XYZ, optional)",
         },
     )
+    settings["force_file"] = (
+        "OptionalXYZFileConfigurator",
+        {
+            "wildcard": "XYZ files (*.xyz);;All files (*)",
+            "default": "",
+            "label": "Force file (XYZ, optional)",
+        },
+    )
     settings["cell_file"] = (
         "InputFileConfigurator",
         {
@@ -169,7 +176,7 @@ def initialize(self):
         self._atomicAliases = self.configuration["atom_aliases"]["value"]
         self._xyzFile = self.configuration["pos_file"]
 
-        if self.configuration["vel_file"]:
+        if self.configuration["vel_file"]["value"]:
             self._velFile = self.configuration["vel_file"]
             if abs(self._xyzFile["time_step"] - self._velFile["time_step"]) > 1.0e-09:
                 raise CP2KConverterError(
@@ -181,6 +188,18 @@ def initialize(self):
                     "Inconsistent number of frames between pos and vel files"
                 )
 
+        if self.configuration["force_file"]["value"]:
+            self._forceFile = self.configuration["force_file"]
+            if abs(self._xyzFile["time_step"] - self._forceFile["time_step"]) > 1.0e-09:
+                raise CP2KConverterError(
+                    "Inconsistent time step between pos and force files"
+                )
+
+            if self._xyzFile["n_frames"] != self._forceFile["n_frames"]:
+                raise CP2KConverterError(
+                    "Inconsistent number of frames between pos and force files"
+                )
+
         self._cellFile = CellFile(self.configuration["cell_file"]["filename"])
 
         if abs(self._cellFile["time_step"] - self._xyzFile["time_step"]) > 1.0e-09:
@@ -216,8 +235,10 @@ def initialize(self):
         )
 
         data_to_be_written = ["configuration", "time"]
-        if self.configuration["vel_file"]:
+        if self.configuration["vel_file"]["value"]:
             data_to_be_written.append("velocities")
+        if self.configuration["force_file"]["value"]:
+            data_to_be_written.append("forces")
 
     def run_step(self, index):
         """Runs a single step of the job.
@@ -236,10 +257,14 @@ def run_step(self, index):
         )
 
         variables = {}
-        if self.configuration["vel_file"]:
+        if self.configuration["vel_file"]["value"]:
             variables["velocities"] = self._velFile.read_step(index) * measure(
                 1.0, iunit="ang/fs"
             ).toval("nm/ps")
+        if self.configuration["force_file"]["value"]:
+            variables["forces"] = self._forceFile.read_step(index) * measure(
+                1.0, iunit="uma ang / fs2"
+            ).toval("uma nm / ps2")
 
         real_conf = PeriodicRealConfiguration(
             self._trajectory.chemical_system, coords, unitcell, **variables
@@ -284,9 +309,12 @@ def finalize(self):
 
         self._xyzFile.close()
 
-        if self.configuration["vel_file"]:
+        if self.configuration["vel_file"]["value"]:
             self._velFile.close()
 
+        if self.configuration["force_file"]["value"]:
+            self._forceFile.close()
+
         self._cellFile.close()
 
         # Close the output trajectory.

MDANSE/Src/MDANSE/Framework/Jobs/GaussianDynamicIncoherentStructureFactor.py

@@ -166,6 +166,15 @@ def initialize(self):
                 main_result=True,
                 partial_result=True,
             )
+            self._outputData.add(
+                f"msd_{element}",
+                "LineOutputVariable",
+                (self._nFrames,),
+                axis="time",
+                units="nm2",
+                main_result=True,
+                partial_result=True,
+            )
 
         self._outputData.add(
             "f(q,t)_total",
@@ -182,20 +191,32 @@ def initialize(self):
             units="nm2/ps",
             main_result=True,
         )
+        self._outputData.add(
+            "msd_total",
+            "LineOutputVariable",
+            (self._nFrames,),
+            axis="time",
+            units="nm2",
+            main_result=True,
+        )
 
         self._atoms = self.configuration["trajectory"][
             "instance"
         ].chemical_system.atom_list
 
-    def run_step(self, index):
-        """
-        Runs a single step of the job.\n
+    def run_step(self, index: int):
+        """Calculates the GDISF and MSD of an atom.
 
-        :Parameters:
-            #. index (int): The index of the step.
-        :Returns:
-            #. index (int): The index of the step.
-            #. atomicSF (np.array): The atomic structure factor
+        Parameters
+        ----------
+        index : int
+            The index of the atom that the calculation will be run over.
+
+        Returns
+        -------
+        tuple[int, tuple[np.ndarray, np.ndarray]]
+            A tuple which contains the job index and a tuple of the
+            GDISF and MSD of an atom.
         """
 
         # get atom index
@@ -218,29 +239,29 @@ def run_step(self, index):
 
         for i, q2 in enumerate(self._kSquare):
             gaussian = np.exp(-msd * q2 / 6.0)
-
             atomicSF[i, :] += gaussian
 
-        return index, atomicSF
+        return index, (atomicSF, msd)
 
-    def combine(self, index, x):
-        """
-        Combines returned results of run_step.\n
-        :Parameters:
-            #. index (int): The index of the step.\n
-            #. x (any): The returned result(s) of run_step
-        """
+    def combine(self, index: int, x: tuple[np.ndarray, np.ndarray]):
+        """Add the results to the output files.
 
-        # The symbol of the atom.
+        Parameters
+        ----------
+        index : int
+            The atom index that the calculation was run over.
+        x : tuple[np.ndarray, np.ndarray]
+            A tuple of the GDISF and MSD of an atom.
+        """
         element = self.configuration["atom_selection"]["names"][index]
-
-        self._outputData[f"f(q,t)_{element}"] += x
+        atomicSF, msd = x
+        self._outputData[f"f(q,t)_{element}"] += atomicSF
+        self._outputData[f"msd_{element}"] += msd
 
     def finalize(self):
         """
         Finalizes the calculations (e.g. averaging the total term, output files creations ...)
         """
-
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         for element, number in nAtomsPerElement.items():
             self._outputData[f"f(q,t)_{element}"][:] /= number
@@ -250,6 +271,8 @@ def finalize(self):
                 self.configuration["instrument_resolution"]["time_step"],
                 axis=1,
             )
+            self._outputData[f"msd_{element}"][:] /= number
+
         weights = self.configuration["weights"].get_weights()
         weight_dict = get_weights(weights, nAtomsPerElement, 1)
         assign_weights(self._outputData, weight_dict, "f(q,t)_%s")
@@ -259,13 +282,23 @@ def finalize(self):
             weight_dict,
             "f(q,t)_%s",
         )
-
         self._outputData["s(q,f)_total"][:] = weighted_sum(
             self._outputData,
             weight_dict,
             "s(q,f)_%s",
         )
 
+        # since GDISF ~ exp(-msd * q2 / 6.0) the MSD isn't weighted in
+        # the exp lets save the MSD with equal weights
+        weights = self.configuration["weights"].get_weights("equal")
+        weight_dict = get_weights(weights, nAtomsPerElement, 1)
+        assign_weights(self._outputData, weight_dict, "msd_%s")
+        self._outputData["msd_total"][:] = weighted_sum(
+            self._outputData,
+            weight_dict,
+            "msd_%s",
+        )
+
         self._outputData.write(
             self.configuration["output_files"]["root"],
             self.configuration["output_files"]["formats"],

MDANSE/Src/MDANSE/Mathematics/Arithmetic.py

@@ -98,7 +98,6 @@ def assign_weights(
     dim = key.count("%s")
 
     for k in values.keys() & matches:
-
         if symmetric:
             permutations = set(itertools.permutations(matches[k], r=dim))
             w = sum(weights[p] for p in permutations)

MDANSE/Tests/UnitTests/Analysis/test_scattering.py

@@ -291,12 +291,12 @@ def test_gdisf(traj_info):
     result_file = os.path.join(result_dir, f"gdisf_{traj_info[0]}.mda")
     with h5py.File(temp_name + ".mda") as actual, h5py.File(result_file) as desired:
         keys = [
-            i for i in desired.keys() if any([j in i for j in ["f(q,t)", "s(q,f)"]])
+            i for i in desired.keys() if any([j in i for j in ["f(q,t)", "s(q,f)", "msd"]])
         ]
         for key in keys:
             np.testing.assert_array_almost_equal(
                 actual[f"/{key}"] * actual[f"/{key}"].attrs["scaling_factor"],
-                desired[f"/{key}"],
+                desired[f"/{key}"] * desired[f"/{key}"].attrs["scaling_factor"],
             )
 
     os.remove(temp_name + ".mda")

MDANSE/Tests/UnitTests/Data/SrTiO3_MD-cell-1.cell

@@ -0,0 +1,11 @@
+#   Step   Time [fs]       Ax [Angstrom]       Ay [Angstrom]       Az [Angstrom]       Bx [Angstrom]       By [Angstrom]       Bz [Angstrom]       Cx [Angstrom]       Cy [Angstrom]       Cz [Angstrom]      Volume [Angstrom^3]
+2 0.6 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+3 1.2 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+4 1.8 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+5 2.4 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+6 3.0 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+7 3.6 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+8 4.2 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+9 4.8 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+10 5.4 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875
+11 6.0 11.715 0.0 0.0 0.0 11.715 0.0 0.0 0.0 11.715 1607.780950875