Merge pull request #175 from CCPBioSim/172-coordinates-and-forces-in-separate-files

Add Support for Calculating Entropy When Forces and Coordinates Are Within Seperate Files

Author: skfegan

CodeEntropy/config/arg_config_manager.py

@@ -12,7 +12,22 @@
     "top_traj_file": {
         "type": str,
         "nargs": "+",
-        "help": "Path to Structure/topology file followed by Trajectory file(s)",
+        "help": "Path to structure/topology file followed by trajectory file",
+    },
+    "force_file": {
+        "type": str,
+        "default": None,
+        "help": "Optional path to force file if forces are not in trajectory file",
+    },
+    "file_format": {
+        "type": str,
+        "default": None,
+        "help": "String for file format as recognised by MDAnalysis",
+    },
+    "kcal_force_units": {
+        "type": bool,
+        "default": False,
+        "help": "Set this to True if you have a separate force file with nonSI units.",
     },
     "selection_string": {
         "type": str,

CodeEntropy/entropy.py

@@ -86,7 +86,7 @@ def execute(self):
         )
 
         reduced_atom, number_molecules, levels, groups = self._initialize_molecules()
-
+        logger.debug(f"Universe 3: {reduced_atom}")
         water_atoms = self._universe.select_atoms("water")
         water_resids = set(res.resid for res in water_atoms.residues)
 

CodeEntropy/run.py

@@ -247,8 +247,43 @@ def run_entropy_workflow(self):
                 # Load MDAnalysis Universe
                 tprfile = args.top_traj_file[0]
                 trrfile = args.top_traj_file[1:]
+                fileformat = args.file_format
                 logger.debug(f"Loading Universe with {tprfile} and {trrfile}")
-                u = mda.Universe(tprfile, trrfile)
+                u = mda.Universe(tprfile, trrfile, format=fileformat)
+
+                # If forces are in separate file merge them with the
+                # coordinates from the trajectory file
+                forcefile = args.force_file
+                if forcefile is not None:
+                    logger.debug(f"Loading Universe with {forcefile}")
+                    u_force = mda.Universe(tprfile, forcefile, format=fileformat)
+                    select_atom = u.select_atoms("all")
+                    select_atom_force = u_force.select_atoms("all")
+
+                    coordinates = (
+                        AnalysisFromFunction(
+                            lambda ag: ag.positions.copy(), select_atom
+                        )
+                        .run()
+                        .results["timeseries"]
+                    )
+                    forces = (
+                        AnalysisFromFunction(
+                            lambda ag: ag.positions.copy(), select_atom_force
+                        )
+                        .run()
+                        .results["timeseries"]
+                    )
+
+                    if args.kcal_force_units:
+                        # Convert from kcal to kJ
+                        forces *= 4.184
+
+                    logger.debug("Merging forces with coordinates universe.")
+                    new_universe = mda.Merge(select_atom)
+                    new_universe.load_new(coordinates, forces=forces)
+
+                    u = new_universe
 
                 self._config_manager.input_parameters_validation(u, args)
 
@@ -311,15 +346,8 @@ def new_U_select_frame(self, u, start=None, end=None, step=1):
             .run()
             .results["timeseries"][start:end:step]
         )
-        dimensions = (
-            AnalysisFromFunction(lambda ag: ag.dimensions.copy(), select_atom)
-            .run()
-            .results["timeseries"][start:end:step]
-        )
         u2 = mda.Merge(select_atom)
-        u2.load_new(
-            coordinates, format=MemoryReader, forces=forces, dimensions=dimensions
-        )
+        u2.load_new(coordinates, format=MemoryReader, forces=forces)
         logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
         return u2
 
@@ -351,15 +379,8 @@ def new_U_select_atom(self, u, select_string="all"):
             .run()
             .results["timeseries"]
         )
-        dimensions = (
-            AnalysisFromFunction(lambda ag: ag.dimensions.copy(), select_atom)
-            .run()
-            .results["timeseries"]
-        )
         u2 = mda.Merge(select_atom)
-        u2.load_new(
-            coordinates, format=MemoryReader, forces=forces, dimensions=dimensions
-        )
+        u2.load_new(coordinates, format=MemoryReader, forces=forces)
         logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
         return u2
 

docs/config.yaml

@@ -2,6 +2,8 @@
 
 run1:
   top_traj_file: ["1AKI_prod.tpr", "1AKI_prod.trr"]
+  force_file:
+  file_formate: 
   selection_string: 'all'
   start: 0
   end: 500

docs/faq.rst

@@ -4,7 +4,15 @@ Frequently asked questions
 Why do I get a ``WARNING`` about invalid eigenvalues?
 -----------------------------------------------------
 
-Insufficient sampling might introduce noise and cause matrix elements to deviate to values that would not reflect the uncorrelated nature of force-force covariance of distantly positioned residues.Try increasing the sampling time. This is especially true at the residue level. 
+Insufficient sampling might introduce noise and cause matrix elements to deviate to values that would not reflect the uncorrelated nature of force-force covariance of distantly positioned residues. 
+Try increasing the sampling time. 
+This is especially true at the residue level. 
 
 For example in a lysozyme system, the residue level contains the largest force and torque covariance matrices because at this level we have the largest number of beads (which is equal to the number of residues in a protein) compared to the molecule level (3 beads) and united-atom level (~10 beads per amino acid). 
 
+What do I do if there is an error from MDAnalysis not recognising the file type?
+-------------------------------------------------------------------------------
+
+Use the file_format option. 
+The MDAnalysis documentation has a list of acceptable formats: https://userguide.mdanalysis.org/1.1.1/formats/index.html#id1.
+The first column of their table gives you the string you need (case sensitive).

docs/getting_started.rst

@@ -81,10 +81,18 @@ The top_traj_file argument is necessary to identify your simulation data, the ot
      - Description
      - Default
      - Type
-   * - ``--top_traj_file`` 
-     - Path to Structure/topology file followed by Trajectory file(s). Any MDAnalysis readable files should work  (for example ``GROMACS TPR and TRR`` or ``AMBER PRMTOP and NETCDF``). You will need to output the **coordinates** and **forces** to the **same file** . 
+   * - ``--top_traj_file``
+     - Path to Structure/topology file followed by Trajectory file. Any MDAnalysis readable files should work  (for example ``GROMACS TPR and TRR`` or ``AMBER PRMTOP and NETCDF``). 
      - Required, no default value
      - list of ``str`` 
+   * - ``--force_file``
+     - Path to a file with forces. This option should be used if the forces are not in the same file as the coordinates. It is expected that the force file has the same number of atoms and frames as the trajectory file. Any MDAnalysis readable files should work  (for example ``AMBER NETCDF`` or ``LAMMPS DCD``). 
+     - None
+     - ``str`` 
+   * - ``--file_format``
+     - Use to tell MDAnalysis the format if the trajectory or force file does not have the standard extension recognised by MDAnalysis.
+     - None
+     - ``str`` 
    * - ``--selection_string``
      - Selection string for CodeEntropy such as protein or resid, refer to ``MDAnalysis.select_atoms`` for more information.
      - ``"all"``: select all atom in trajectory