Merge pull request #23 from nomad-coe/18-add-shear-flow-md-parameters

18 add shear flow md parameters

simulationworkflowschema/molecular_dynamics.py

@@ -1564,6 +1564,68 @@ class BarostatParameters(ArchiveSection):
     )
 
 
+class ShearParameters(ArchiveSection):
+    """
+    Section containing the parameters pertaining to the shear flow for a molecular dynamics run.
+    """
+
+    m_def = Section(validate=False)
+
+    shear_type = Quantity(
+        type=MEnum('lees_edwards', 'trozzi_ciccotti', 'ashurst_hoover'),
+        shape=[],
+        description="""
+        The name of the method used to implement the effect of shear flow within the simulation.
+
+        Allowed values are:
+
+        | Shear Method          | Description                               |
+
+        | ---------------------- | ----------------------------------------- |
+
+        | `""`                   | No thermostat               |
+
+        | `"lees_edwards"`          | A.W. Lees and S.F. Edwards,
+        [J. Phys. C **5** (1972) 1921](https://doi.org/10.1088/0022-3719/5/15/006)|
+
+        | `"trozzi_ciccotti"`          | A.W. Lees and S.F. Edwards,
+        [Phys. Rev. A **29** (1984) 916](https://doi.org/10.1103/PhysRevA.29.916)|
+
+        | `"ashurst_hoover"`          | W. T. Ashurst and W. G. Hoover,
+        [Phys. Rev. A **11** (1975) 658](https://doi.org/10.1103/PhysRevA.11.658)|
+        """,
+    )
+
+    shear_rate = Quantity(
+        type=np.float64,
+        shape=[3, 3],
+        unit='ps^-1',
+        description="""
+        The external stress tensor include normal (diagonal elements; which are zero in shear simulations)
+        and shear stress' rates (off-diagonal elements).
+        Its elements are: [[σ_x, τ_yx, τ_zx], [τ_xy, σ_y, τ_zy], [τ_xz, τ_yz, σ_z]],
+		where σ and τ are the normal and shear stress' rates.
+        The first and second letters in the index correspond to the normal vector to the shear plane and the direction of shearing, respectively.
+        """,
+    )
+
+    step_start = Quantity(
+        type=int,
+        shape=[],
+        description="""
+        Trajectory step where this shearing starts.
+        """,
+    )
+
+    step_end = Quantity(
+        type=int,
+        shape=[],
+        description="""
+        Trajectory step number where this shearing ends.
+        """,
+    )
+
+
 class Lambdas(ArchiveSection):
     """
     Section for storing all lambda parameters for free energy perturbation
@@ -1834,6 +1896,8 @@ class MolecularDynamicsMethod(SimulationWorkflowMethod):
 
     barostat_parameters = SubSection(sub_section=BarostatParameters.m_def, repeats=True)
 
+    shear_parameters = SubSection(sub_section=ShearParameters.m_def, repeats=True)
+
     free_energy_calculation_parameters = SubSection(
         sub_section=FreeEnergyCalculationParameters.m_def, repeats=True
     )
@@ -2160,7 +2224,7 @@ def normalize(self, archive, logger):
             self.label = self._rg_results.get('label')
             # TODO Fix this assignment fails with TypeError
             try:
-                self.atomsgroup_ref = self._rg_results.get('atomsgroup_ref')
+                self.atomsgroup_ref = [self._rg_results.get('atomsgroup_ref')]
             except Exception:
                 pass
             self.n_frames = self._rg_results.get('n_frames')
@@ -2203,18 +2267,10 @@ class FreeEnergyCalculations(TrajectoryProperty):
         """,
     )
 
-    value_unit = Quantity(
-        type=str,
-        shape=[],
-        description="""
-        Unit of the property, using UnitRegistry() notation.
-        In this case, the unit corresponds to all `value` properties stored within this section.
-        """,
-    )
-
     value_total_energy_magnitude = Quantity(
         type=HDF5Dataset,
         shape=[],
+        unit='joule',
         description="""
         Value of the total energy for the present lambda state. The expected dimensions are ["n_frames"].
         This quantity is a reference to the data (file+path), which is stored in an HDF5 file for efficiency.
@@ -2224,6 +2280,7 @@ class FreeEnergyCalculations(TrajectoryProperty):
     value_PV_energy_magnitude = Quantity(
         type=HDF5Dataset,
         shape=[],
+        unit='joule',
         description="""
         Value of the pressure-volume energy (i.e., P*V) for the present lambda state. The expected dimensions are ["n_frames"].
         This quantity is a reference to the data (file+path), which is stored in an HDF5 file for efficiency.
@@ -2233,6 +2290,7 @@ class FreeEnergyCalculations(TrajectoryProperty):
     value_total_energy_differences_magnitude = Quantity(
         type=HDF5Dataset,
         shape=[],
+        unit='joule',
         description="""
         Values correspond to the difference in total energy between each specified lambda state
         and the reference state, which corresponds to the value of lambda of the current simulation.
@@ -2244,6 +2302,7 @@ class FreeEnergyCalculations(TrajectoryProperty):
     value_total_energy_derivative_magnitude = Quantity(
         type=HDF5Dataset,
         shape=[],
+        unit='joule',  # TODO check this unit
         description="""
         Value of the derivative of the total energy with respect to lambda, evaluated for the current
         lambda state. The expected dimensions are ["n_frames"].
@@ -2479,7 +2538,14 @@ class MolecularDynamicsResults(ThermodynamicsResults):
     def normalize(self, archive, logger):
         super().normalize(archive, logger)
 
-        universe = archive_to_universe(archive)
+        try:
+            universe = archive_to_universe(archive)
+        except Exception:
+            universe = None
+            logger.warning(
+                'Could not convert archive to MDAnalysis Universe, skipping MD results normalization.'
+            )
+
         if universe is None:
             return
 
@@ -2570,9 +2636,10 @@ def normalize(self, archive, logger):
                     sec_rg_values = sec_rgs_calc.m_create(
                         RadiusOfGyrationValuesCalculation
                     )
+
                     # TODO Fix this assignment fails with TypeError
                     try:
-                        sec_rg_values.atomsgroup_ref = rg.get('atomsgroup_ref')
+                        sec_rg_values.atomsgroup_ref = [rg.get('atomsgroup_ref')]
                     except Exception:
                         pass
                     sec_rg_values.label = rg.get('label')