Save Q vectors in the output file

MDANSE/Src/MDANSE/Framework/Configurators/CorrelationFramesConfigurator.py

@@ -19,15 +19,21 @@
 
 
 class CorrelationFramesConfigurator(FramesConfigurator):
+    """Parses the input of trajectory frames.
+
+    Configures the time frame range to be used in the calculations
+    together with a movable window used for correlations.
+    """
+
     def configure(self, value: tuple[int, int, int, int]):
-        """Configure the correlation and set the number of correlation
-        frames to use.
+        """Set the number of correlation frames to use.
 
         Parameters
         ----------
         value : tuple[int, int, int, int]
             The frames setting plus the number of frames used for the
             correlations.
+
         """
         trajConfig = self._configurable[self._dependencies["trajectory"]]
         n_steps = trajConfig["length"]
@@ -38,6 +44,7 @@ def configure(self, value: tuple[int, int, int, int]):
 
         first, last, step, c_frames = value
         super().configure((first, last, step))
+        self._original_input = value
 
         if c_frames > self["n_frames"]:
             self.error_status = (

MDANSE/Src/MDANSE/Framework/Configurators/QVectorsConfigurator.py

@@ -14,56 +14,55 @@
 #    along with this program.  If not, see <https://www.gnu.org/licenses/>.
 #
 
+from typing import Any
 
 from MDANSE.Framework.Configurators.IConfigurator import IConfigurator
 from MDANSE.Framework.QVectors.IQVectors import IQVectors
 
 
 class QVectorsConfigurator(IConfigurator):
-    """
-    This Configurator allows to set reciprocal vectors.
-
-    Reciprocal vectors are used in MDANSE for the analysis related to scattering experiments such as dynamic coherent structure
-    or elastic incoherent structure factor analysis. In MDANSE, properties that depends on Q vectors are always scalar regarding
-    Q vectors in the sense that the values of these properties will be computed for a given norm of Q vectors and not for a given Q vectors.
-    Hence, the Q vectors generator supported by MDANSE always generates Q vectors on Q-shells, each shell containing a set of Q vectors whose
+    """Creates and configures a q-vector generator.
+
+    Reciprocal vectors are used in MDANSE for analysis related to
+    scattering experiments, such as dynamic coherent structure
+    or elastic incoherent structure factor analysis. In MDANSE, properties
+    that depend on Q vectors are always scalar regarding Q vectors
+    in the sense that the values of these properties will be computed
+    for a given norm of Q vectors and not for a given Q vector.
+    Hence, the Q vectors generator supported by MDANSE always generates
+    Q vectors on Q-shells, each shell containing a set of Q vectors whose
     norm match the Q shell value within a given tolerance.
 
-    Depending on the generator selected, Q vectors can be generated isotropically or anistropically, on a lattice or randomly.
-
-    Q vectors can be saved to a user definition and, as such, can be further reused in another MDANSE session.
+    Depending on the generator selected, Q vectors can be generated
+    isotropically or anistropically, on a lattice or randomly.
 
-    To define a new Q vectors generator, you must inherit from MDANSE.Framework.QVectors.QVectors.QVector interface.
-
-    :note: this configurator depends on 'trajectory' configurator to be configured.
     """
 
     _default = (
         "SphericalLatticeQVectors",
         {"shells": (0.1, 5, 0.1), "width": 0.1, "n_vectors": 50, "seed": 0},
     )
 
-    def configure(self, value):
-        """
-        Configure a Q vectors generator.
+    def configure(self, value: tuple[str, dict[str, Any]]):
+        """Create a vector generator with given parameters.
+
+        Parameters
+        ----------
+        value : tuple[str, dict[str, Any]]
+            Class name and dictionary of input parameters
 
-        :param configuration: the current configuration.
-        :type configuration: a MDANSE.Framework.Configurable.Configurable object
-        :param value: the Q vectors generator definition. It can be a 2-tuple, whose 1st element is the name of the Q vector generator \
-        and 2nd element the parameters for this configurator or a string that matches a Q vectors user definition.
-        :type value: 2-tuple or str
         """
         self._original_input = value
 
         trajConfig = self._configurable[self._dependencies["trajectory"]]
         if isinstance(value, tuple):
             try:
-                generator, parameters = value
+                generator_name, parameters = value
             except ValueError:
                 self.error_status = f"Invalid q vectors settings {value}"
                 return
             generator = IQVectors.create(
-                generator, trajConfig["instance"].configuration(0)
+                generator_name, trajConfig["instance"].configuration(0)
             )
             try:
                 generator.setup(parameters)
@@ -84,7 +83,7 @@ def configure(self, value):
             if "q_vectors" not in generator.configuration:
                 self.error_status = "Wrong inputs for q-vector generation. At the moment there are no valid Q points."
                 return
-            elif not generator.configuration["q_vectors"]:
+            if not generator.configuration["q_vectors"]:
                 self.error_status = "no Q vectors could be generated"
                 return
 
@@ -99,23 +98,33 @@ def configure(self, value):
         self["shells"] = list(self["q_vectors"].keys())
         self["n_shells"] = len(self["q_vectors"])
         self["value"] = self["q_vectors"]
+        self["generator"] = generator
         self.error_status = "OK"
 
     def preview_output_axis(self):
+        """Output the values of |Q| from current parameters.
+
+        Returns
+        -------
+        list[float], str
+            Values of |Q| together with their physical unit
+
+        """
         if not self.is_configured():
             return None, None
         if not self._valid:
             return None, None
         return self["shells"], "1/nm"
 
-    def get_information(self):
-        """
-        Returns string information about this configurator.
+    def get_information(self) -> str:
+        """Return human-readable information about the generated vectors.
 
-        :return: the information about this configurator.
-        :rtype: str
-        """
+        Returns
+        -------
+        str
+            Summary of generated vectors.
 
+        """
         try:
             info = [f"{self['n_shells']} Q shells generated\n"]
         except KeyError:

MDANSE/Src/MDANSE/Framework/Jobs/CurrentCorrelationFunction.py

@@ -370,6 +370,10 @@ def finalize(self):
         """Normalize, Fourier transform and write the results out to
         the MDA files.
         """
+        self.configuration["q_vectors"]["generator"].write_vectors_to_file(
+            self._outputData
+        )
+
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         for pair in self._elementsPairs:
             pair_str = "".join(map(str, pair))

MDANSE/Src/MDANSE/Framework/Jobs/DynamicCoherentStructureFactor.py

@@ -20,10 +20,14 @@
 import numpy as np
 from scipy.signal import correlate
 
+from MDANSE.MLogging import LOG
 from MDANSE.Core.Error import Error
 from MDANSE.Framework.Jobs.IJob import IJob
+from MDANSE.Framework.QVectors.LatticeQVectors import LatticeQVectors
 from MDANSE.Mathematics.Arithmetic import assign_weights, get_weights, weighted_sum
 from MDANSE.Mathematics.Signal import get_spectrum
+from MDANSE.Framework.QVectors.IQVectors import IQVectors
+from MDANSE.MolecularDynamics.UnitCell import UnitCell
 
 
 class DynamicCoherentStructureFactorError(Error):
@@ -99,6 +103,12 @@ def initialize(self):
         self.numberOfSteps = self.configuration["q_vectors"]["n_shells"]
 
         nQShells = self.configuration["q_vectors"]["n_shells"]
+        if not isinstance(
+            self.configuration["q_vectors"]["generator"], LatticeQVectors
+        ):
+            LOG.warning(
+                f"This task should be used with a lattice-based Q vector generator. You have picked {self.configuration['q_vectors']['generator'].__class__}. The results are likely to be incorrect."
+            )
 
         self._nFrames = self.configuration["frames"]["n_frames"]
 
@@ -182,6 +192,28 @@ def initialize(self):
             main_result=True,
         )
 
+        self._cell_std = 0.0
+        try:
+            all_cells = [
+                self.configuration["trajectory"]["instance"].unit_cell(frame)._unit_cell
+                for frame in self.configuration["frames"]["value"]
+            ]
+        except TypeError:
+            self._average_unit_cell = None
+        else:
+            self._average_unit_cell = UnitCell(
+                np.mean(
+                    all_cells,
+                    axis=0,
+                )
+            )
+            self._cell_std = UnitCell(
+                np.std(
+                    all_cells,
+                    axis=0,
+                )
+            )
+
     def run_step(self, index):
         """
         Runs a single step of the job.\n
@@ -212,9 +244,35 @@ def run_step(self, index):
                     dtype=np.complex64,
                 )
 
+            cell_present = True
+            cell_fixed = True
             # loop over the trajectory time steps
             for i, frame in enumerate(self.configuration["frames"]["value"]):
-                qVectors = self.configuration["q_vectors"]["value"][shell]["q_vectors"]
+                unit_cell = traj.unit_cell(frame)
+                if unit_cell is None:
+                    cell_present = False
+                elif not np.allclose(
+                    unit_cell._unit_cell, self._average_unit_cell._unit_cell
+                ):
+                    cell_fixed = False
+                if not cell_present:
+                    qVectors = self.configuration["q_vectors"]["value"][shell][
+                        "q_vectors"
+                    ]
+                else:
+                    try:
+                        hkls = self.configuration["q_vectors"]["value"][shell]["hkls"]
+                    except KeyError:
+                        qVectors = self.configuration["q_vectors"]["value"][shell][
+                            "q_vectors"
+                        ]
+                    else:
+                        if hkls is None:
+                            qVectors = self.configuration["q_vectors"]["value"][shell][
+                                "q_vectors"
+                            ]
+                        else:
+                            qVectors = IQVectors.hkl_to_qvectors(hkls, unit_cell)
 
                 coords = traj.configuration(frame)["coordinates"]
 
@@ -223,6 +281,14 @@ def run_step(self, index):
                     rho[element][i, :] = np.sum(
                         np.exp(1j * np.dot(selectedCoordinates, qVectors)), axis=0
                     )
+            if not cell_present:
+                LOG.warning(
+                    "You are running the DCSF calculation on a trajectory without periodic boundary conditions."
+                )
+            if not cell_fixed:
+                LOG.warning(
+                    f"The unit cell is VARIABLE with the standard deviation of {self._cell_std}. This analysis should not be used with NPT runs! PLEASE CHECK YOUR RESULTS CAREFULLY."
+                )
 
             return index, rho
 
@@ -249,6 +315,10 @@ def finalize(self):
         """
         Finalizes the calculations (e.g. averaging the total term, output files creations ...)
         """
+        self.configuration["q_vectors"]["generator"].write_vectors_to_file(
+            self._outputData
+        )
+
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         weights = self.configuration["weights"].get_weights()
         weight_dict = get_weights(weights, nAtomsPerElement, 2, conc_exp=0.5)

MDANSE/Src/MDANSE/Framework/Jobs/DynamicIncoherentStructureFactor.py

@@ -251,6 +251,9 @@ def finalize(self):
         """
         Finalizes the calculations (e.g. averaging the total term, output files creations ...)
         """
+        self.configuration["q_vectors"]["generator"].write_vectors_to_file(
+            self._outputData
+        )
 
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         weights = self.configuration["weights"].get_weights()

MDANSE/Src/MDANSE/Framework/Jobs/ElasticIncoherentStructureFactor.py

@@ -194,6 +194,10 @@ def finalize(self):
         Finalizes the calculations (e.g. averaging the total term, output files creations ...)
         """
 
+        self.configuration["q_vectors"]["generator"].write_vectors_to_file(
+            self._outputData
+        )
+
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         for element, number in list(nAtomsPerElement.items()):
             self._outputData[f"eisf_{element}"][:] /= number

MDANSE/Src/MDANSE/Framework/QVectors/ApproximateDispersionQVectors.py

@@ -24,7 +24,7 @@
 
 
 class ApproximateDispersionQVectors(LatticeQVectors):
-    """ """
+    """Generates Q vectors along a direction."""
 
     settings = collections.OrderedDict()
     settings["q_start"] = (
@@ -69,7 +69,7 @@ def _generate(self):
             np.array(qStart)[:, np.newaxis] + np.outer(n, np.arange(0, nSteps)) * qStep
         )
 
-        hkls = np.rint(np.dot(self._directUnitCell, vects))
+        hkls = self.qvectors_to_hkl(vects, self._unit_cell)
 
         dists = np.sqrt(np.sum(vects**2, axis=0))
         dists = list(zip(range(len(dists)), dists))
@@ -94,5 +94,4 @@ def _generate(self):
             if self._status is not None:
                 if self._status.is_stopped():
                     return
-                else:
-                    self._status.update()
+                self._status.update()

MDANSE/Src/MDANSE/Framework/QVectors/CircularLatticeQVectors.py

@@ -19,13 +19,21 @@
 
 import numpy as np
 
-from MDANSE.Mathematics.LinearAlgebra import Vector
-
 from MDANSE.Framework.QVectors.LatticeQVectors import LatticeQVectors
+from MDANSE.Mathematics.LinearAlgebra import Vector
 
 
 class CircularLatticeQVectors(LatticeQVectors):
-    """ """
+    """Generates Q vectors on a plane.
+
+    Vectors are grouped into annuli (called 'shells')
+    based on their length.
+
+    Only vectors commensurate with the reciprocal
+    space lattice vectors will be generated.
+    |Q| values for which no valid vectors can
+    be found are omitted in the output.
+    """
 
     settings = collections.OrderedDict()
     settings["seed"] = ("IntegerConfigurator", {"mini": 0, "default": 0})
@@ -61,7 +69,7 @@ def _generate(self):
             ]
         )
 
-        qVects = np.dot(self._inverseUnitCell, hkls)
+        qVects = self.hkl_to_qvectors(hkls, self._unit_cell)
 
         qMax = (
             self._configuration["shells"]["last"]
@@ -109,15 +117,11 @@ def _generate(self):
                 self._configuration["q_vectors"][q]["q_vectors"] = vects[:, hits]
                 self._configuration["q_vectors"][q]["n_q_vectors"] = n
                 self._configuration["q_vectors"][q]["q"] = q
-                self._configuration["q_vectors"][q]["hkls"] = np.rint(
-                    np.dot(
-                        self._directUnitCell,
-                        self._configuration["q_vectors"][q]["q_vectors"],
-                    )
+                self._configuration["q_vectors"][q]["hkls"] = self.qvectors_to_hkl(
+                    vects[:, hits], self._unit_cell
                 )
 
             if self._status is not None:
                 if self._status.is_stopped():
                     return
-                else:
-                    self._status.update()
+                self._status.update()