Infrared calculation for bulk systems

Doc/pages/analysis.rst

@@ -8,8 +8,9 @@ Analysis: Other
 
 This section contains background theory for following plugins:
 
--  :ref:`infrared`
 -  :ref:`dipole-autocorrelation-function`
+-  :ref:`infraredbulk`
+-  :ref:`infraredmolecular`
 -  :ref:`density`
 -  :ref:`temperature`
 -  :ref:`center-of-masses-trajectory`
@@ -18,23 +19,6 @@ This section contains background theory for following plugins:
 Infrared
 ^^^^^^^^
 
-.. _infrared:
-
-Infrared
-''''''''
-Calculates the molecular infrared spectrum averaged over all molecules
-in the trajectory. The infrared spectrum is calculated from the Fourier
-transform of the autocorrelation of the time-derivative of the
-molecular dipole:
-
-.. math::
-   :label: ir1
-
-   I(\omega) \propto  \frac{1}{N_{m}}\sum_{m} \frac{1}{6\pi} \int \mathrm{d}t \,  \left\langle \dot{\vec{\mu}}_{m}(0) \cdot \dot{\vec{\mu}}_{m}(t) \right\rangle e^{-i\omega t}
-
-where :math:`N_{m}` is the number of molecules and :math:`\dot{\vec{\mu}}_{m}(t)` is
-the time-derivative of the molecular dipole moment of molecule :math:`m`.
-
 .. _dipole-autocorrelation-function:
 
 Dipole Autocorrelation Function
@@ -44,12 +28,70 @@ Calculates the molecular dipole autocorrelation function which is closely
 related to the molecular infrared spectrum
 
 .. math::
-   :label: ir2
+   :label: ir1
 
    \mathrm{DACF}(t) = \frac{1}{3 N_{m}}\sum_{m} \left\langle \vec{\mu}_{m}(0) \cdot \vec{\mu}_{m}(t) \right\rangle
 
 where :math:`N_{m}` is the number of molecules :math:`m` and :math:`\vec{\mu}(t)` is
-the molecular dipole moment of molecule :math:`m`.
+the molecular dipole moment of molecule :math:`m`. The electric dipole of a
+given molecule is calculated relative to the molecules COM
+
+.. math::
+   :label: ir2
+
+   \vec{\mu}_{m}(t) = \sum_{i \in m} q_{i}(t)(\mathbf{r}_{i}(t) - \mathbf{r}_{\mathrm{COM},m}).
+
+.. _infraredbulk:
+
+InfraredBulk
+''''''''''''
+Calculates the infrared spectrum, should usually be used for systems where
+molecules are not defined and the system is amorphous. The infrared spectrum
+is calculated following the equation
+
+.. math::
+   :label: ir3
+
+   I(\omega) = \frac{1}{N}\sum_{i \in \mathrm{u.c.}} \sum_{\substack{j \\ \vert \mathbf{r}_j-\mathbf{r}_i \vert < r}} \frac{1}{6\pi} \int \mathrm{d}t \,  \left\langle \dot{\vec{\mu}}_{i}(0) \cdot \dot{\vec{\mu}}_{j}(t) \right\rangle e^{-i\omega t}
+
+where
+
+.. math::
+   :label: ir4
+
+   \dot{\vec{\mu}}_{i}(t) = \frac{\mathrm{d}}{\mathrm{d}t} [q(t)\mathbf{r}(t)]
+
+and the summation over :math:`i` goes over all atoms in the unit cell
+and the summation over :math:`j` goes over all atoms which have approached
+atom :math:`i`. We define atom :math:`j` to have approached :math:`i` when
+it has reached a distance less than the user defined cutoff, :math:`r`, at
+any point during the trajectory. We therefore assume that the correlation function
+
+.. math::
+   :label: ir5
+
+   \Big\langle \dot{\vec{\mu}}_{i}(0) \cdot  \dot{\vec{\mu}}_{j}(t) \Big\rangle
+
+is small when atoms :math:`i` and :math:`j` are separeted by large distances
+which may not be true for periodic systems where the motion of atoms
+separeted by large distances may still be correlated.
+
+.. _infraredmolecular:
+
+InfraredMolecular
+'''''''''''''''''
+Calculates the molecular infrared spectrum averaged over all molecules
+in the trajectory. The infrared spectrum is calculated from the Fourier
+transform of the autocorrelation of the time-derivative of the
+molecular dipole:
+
+.. math::
+   :label: ir6
+
+   I(\omega) = \frac{1}{N_{m}}\sum_{m} \frac{1}{6\pi} \int \mathrm{d}t \,  \left\langle \dot{\vec{\mu}}_{m}(0) \cdot \dot{\vec{\mu}}_{m}(t) \right\rangle e^{-i\omega t}
+
+where :math:`N_{m}` is the number of molecules and :math:`\dot{\vec{\mu}}_{m}(t)` is
+the time-derivative of the molecular dipole moment of molecule :math:`m`.
 
 
 Thermodynamics

MDANSE/Src/MDANSE/Framework/Configurators/DistCutoffConfigurator.py

@@ -0,0 +1,89 @@
+#    This file is part of MDANSE.
+#
+#    MDANSE_GUI 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 __future__ import annotations
+
+import numpy as np
+
+from .FloatConfigurator import FloatConfigurator
+
+
+def get_largest_cutoff(traj_config) -> float:
+    """Get the largest cutoff value for the given trajectories
+    unit cells.
+
+    Returns
+    -------
+    traj_config
+        The trajectory configuration object
+    """
+    try:
+        trajectory_array = np.array(
+            [
+                traj_config.unit_cell(frame)._unit_cell
+                for frame in range(len(traj_config))
+            ]
+        )
+    except Exception:
+        return np.linalg.norm(traj_config.min_span)
+
+    if np.allclose(trajectory_array, 0.0):
+        return np.linalg.norm(traj_config.min_span)
+
+    # calculated the radius of the largest sphere that can
+    # fit into the unit cell
+    min_d = np.min(trajectory_array, axis=0)
+    vec_a, vec_b, vec_c = min_d
+
+    cross_bc = np.cross(vec_b, vec_c)
+    cross_ca = np.cross(vec_c, vec_a)
+    cross_ab = np.cross(vec_a, vec_b)
+
+    if any(np.allclose(vec, 0.0) for vec in (cross_bc, cross_ca, cross_ab)):
+        raise ValueError("Trajectory contains invalid unit cell.")
+
+    h_1 = abs(np.dot(vec_a, cross_bc)) / np.linalg.norm(cross_bc)
+    h_2 = abs(np.dot(vec_b, cross_ca)) / np.linalg.norm(cross_ca)
+    h_3 = abs(np.dot(vec_c, cross_ab)) / np.linalg.norm(cross_ab)
+
+    return 0.5 * min(h_1, h_2, h_3)
+
+
+class DistCutoffConfigurator(FloatConfigurator):
+    """.
+
+    It does not allow distances large enough to include
+    the periodic image of any atom in the system.
+    """
+
+    def configure(self, value):
+        """Configure the distance histogram cutoff configurator.
+
+        Parameters
+        ----------
+        value : tuple
+            A tuple of the range parameters.
+        """
+        super().configure(value)
+
+        if float(value) > round(self.get_max_cutoff(), 2):
+            self.error_status = (
+                "The cutoff distance goes into the simulation box periodic images."
+            )
+            return
+
+    def get_max_cutoff(self):
+        traj_config = self.configurable[self.dependencies["trajectory"]]["instance"]
+        return get_largest_cutoff(traj_config)

MDANSE/Src/MDANSE/Framework/Configurators/DistHistCutoffConfigurator.py

@@ -1,9 +1,21 @@
+#    This file is part of MDANSE.
+#
+#    MDANSE_GUI 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 __future__ import annotations
 
-from math import floor
-
-import numpy as np
-
+from .DistCutoffConfigurator import get_largest_cutoff
 from .IConfigurator import PredictionSettings
 from .RangeConfigurator import RangeConfigurator
 
@@ -34,55 +46,14 @@ def configure(self, value):
         if not self.update_needed(value):
             return
 
-        if self._max_value and value[1] > floor(self.get_largest_cutoff() * 100) / 100:
+        if self._max_value and value[1] > round(self.get_max_cutoff(), 2):
             self.error_status = (
                 "The cutoff distance goes into the simulation box periodic images."
             )
             return
 
         super().configure(value)
 
-    def get_largest_cutoff(self) -> float:
-        """Get the largest cutoff value for the given trajectories
-        unit cells.
-
-        Returns
-        -------
-        float
-            The maximum cutoff for the distance histogram job.
-        """
+    def get_max_cutoff(self):
         traj_config = self.configurable[self.dependencies["trajectory"]]["instance"]
-        try:
-            trajectory_array = np.array(
-                [
-                    traj_config.unit_cell(frame)._unit_cell
-                    for frame in range(len(traj_config))
-                ]
-            )
-        except Exception:
-            return np.linalg.norm(traj_config.min_span)
-        else:
-            if np.allclose(trajectory_array, 0.0):
-                return np.linalg.norm(traj_config.min_span)
-            else:
-                # calculated the radius of the largest sphere that can
-                # fit into the unit cell
-                min_d = np.min(trajectory_array, axis=0)
-                vec_a, vec_b, vec_c = min_d
-
-                cross_bc = np.cross(vec_b, vec_c)
-                cross_ca = np.cross(vec_c, vec_a)
-                cross_ab = np.cross(vec_a, vec_b)
-
-                if (
-                    np.allclose(cross_bc, 0.0)
-                    or np.allclose(cross_ca, 0.0)
-                    or np.allclose(cross_ab, 0.0)
-                ):
-                    raise ValueError("Trajectory contains invalid unit cell.")
-
-                h_1 = abs(np.dot(vec_a, cross_bc)) / np.linalg.norm(cross_bc)
-                h_2 = abs(np.dot(vec_b, cross_ca)) / np.linalg.norm(cross_ca)
-                h_3 = abs(np.dot(vec_c, cross_ab)) / np.linalg.norm(cross_ab)
-
-                return 0.5 * min(h_1, h_2, h_3)
+        return get_largest_cutoff(traj_config)

MDANSE/Src/MDANSE/Framework/Configurators/__init__.py

@@ -35,6 +35,9 @@
 from .DerivativeOrderConfigurator import (
     DerivativeOrderConfigurator as DerivativeOrderConfigurator,
 )
+from .DistCutoffConfigurator import (
+    DistCutoffConfigurator as DistCutoffConfigurator,
+)
 from .DistHistCutoffConfigurator import (
     DistHistCutoffConfigurator as DistHistCutoffConfigurator,
 )