Merge branch 'dev'

Bumped version to 0.0.6 and updated changelog

README.md

@@ -31,9 +31,8 @@ These you can either preinstall, or let either pip or conda install for you in t
 1. importlib_resources (only with with python<3.10) 
 1. [MDAnalysis](https://docs.mdanalysis.org/stable/index.html)
 1. [scikit-learn](https://scikit-learn.org/stable/)
-1. Optional: [NetworkX](https://networkx.org/)
-
-NetworkX is only needed for ancestry graph building and is not installed by default. You can run the tool just fine, until you use the `fst analysis ancestry` command, which will throw an error. You can then just install the package any time and ancestry graphing will start working.
+1. [NetworkX](https://networkx.org/)
+1. [pygraphviz](https://networkx.org/)
 
 
 ## Installation
@@ -56,7 +55,7 @@ In the recommended way we will use conda to make a new envirnonment. This way wo
 First, we will make a new conda environment dedicated just for the tool and activate it (if you have [mamba](https://mamba.readthedocs.io/en/latest/) installed, use it for the first command to speed up the process):
 
 ```sh
-conda create -c conda-forge -n fst_env python=3.10 numpy matplotlib mdanalysis scikit-learn networkx
+conda create -c conda-forge -n fst_env python=3.10 numpy matplotlib mdanalysis scikit-learn networkx pygraphviz
 conda activate fst_env
 ```
 
@@ -224,10 +223,12 @@ These options change the on-the-fly transformations that are done to the traject
 
 | Variable | Description | Default value |
 | --- | - | - |
+| `precenter` | Shift the atoms in before unwrapping or putting mols in box, such that a single atom (defined by `precenter_atom`) is centered. If you need to stop multichain system from being split into separate periodic images when unwrapping and wrapping, you can precentre an atom close to the centre of the structure, such that the complex will be away from the edges and the subunits will always be put next to each other. | `False` |
+| `precenter_atom` | The selection string (or index group) which defined a selection of one atom to use for precentering. By default it is None, in which case the atom that is closest to the box centre in the starting structure will be used. | `None` |
 | `unwrap_mols` | Whether to unwrap molecules (make them whole if broken over the PBC). | `False` |
 | `unwrap_sel` | The selection string (or index group) used for unwrapping. Only atoms within this selection will be considered, so if parts of a chain are missing, only those parts that have unbroken bonded graphs will be made whole, each of them separately. In general you should make a selection with at least the backbone connecting whichever parts you want whole. This should be fast enough even with large proteins, but including the water can effect performance badly. The default selection of "protein" should work in most cases.| `"protein"` |
 | `unwrap_starters` | `None` or the selection string (or index group) used as "starters". These atoms will be the starting points of making the molecules whole. In effect, they are guaranteed to be inside the box after the process. If `None`, or if a molecule does not have any atoms in the group, the atom with the smallest index will be used. | `None` |
-| `mols_in_box` | Whether to put centre of mass of molecules back in box after unwrapping. Only the coordinates in `unwrap_sel` are moved and considered for the centre of mass, **however**, unlike for unwrapping, molecules are moved as a whole, even if they are missing parts in between. This option is ignored if `unwrap_mols=False`.| `False` |
+| `mols_in_box` | Whether to put centre of mass of molecules back in box after unwrapping. Only the coordinates in `unwrap_sel` are moved and considered for the centre of mass, **however**, unlike for unwrapping, molecules are moved as a whole, even if they are missing parts in between. In such a case it is best to include at least the backbone of the whole protein in `unwrap_sel`. | `False` |
 
 ##### Clustering
 

changelog.txt

@@ -1,3 +1,8 @@
+v0.0.6
+- Change in behaviour: Wrapping can be done without unwrapping
+- New feature: Added precentering cordinates with a single atom before unwrapping
+- Updated requirements: NetworkX is no longer optional and is installed automatically
+
 v0.0.5
 - BIG CHANGE: clustering will be done if at least epochs_pre_clust epochs have been run (used to be strictly more than)
 - bugfix: current_dir behaviour fixed when the function changed, but xtc was not reread.

setup.cfg

@@ -22,6 +22,8 @@ dependencies =
     matplotlib
     mdanalysis
     scikit-learn
+    networkx
+    pygraphviz
 
 [options]
 include_package_data = True
@@ -31,7 +33,9 @@ install_requires =
     matplotlib>=3.2
     numpy>=1.21
     scipy>=1.1.0
-    MDAnalysis>=2.0
+    MDAnalysis>=2.4
+    networkx>=2.5
+    pygraphviz>=1.5
     importlib_resources>=4.0;python_version<'3.10'
 
 [options.entry_points]

src/__init__.py

@@ -1 +1 @@
-__version__ = "0.0.5"
+__version__ = "0.0.6"

src/analysis/extract.py

@@ -6,6 +6,7 @@
 from .. import inout
 from .. import utils
 from .. import transformations
+from ..exceptions import WrongSelectionSizeError
 
 # Type hints
 from typing import Any, Tuple, Optional, List, Dict
@@ -29,26 +30,49 @@ def load_struct(args: argparse.Namespace) -> Tuple[mda.Universe, AtomGroup, List
         args.index = args.cfg.index_file
     indexes = utils.read_ndx(args.index)
 
-    args.selection = get_cfg_sel(args.selection,    args.cfg)
-    args.sel_unwrap = get_cfg_sel(args.sel_unwrap,   args.cfg, args.selection)
-    args.sel_superpos = get_cfg_sel(
-        args.sel_superpos, args.cfg, args.selection)
+    args.selection = get_cfg_sel(args.selection, args.cfg)
+    args.sel_unwrap = get_cfg_sel(args.sel_unwrap, args.cfg, args.selection)
+    args.sel_superpos = get_cfg_sel(args.sel_superpos, args.cfg,
+                                    args.selection)
+    args.precenter_atom = get_cfg_sel(args.precenter_atom, args.cfg)
 
     if (args.unwrap_starters == "unwrap_starters"):
         args.unwrap_starters = args.cfg.unwrap_starters
 
     print("Loading structure")
-    u = mda.Universe(str(args.cfg.initial_struct[0]))
+    u = mda.Universe(args.cfg.initial_struct[0])
     sel = utils.load_sel(args.selection, u, indexes)
     print("Selected %d atoms for extraction" % len(sel))
     sel_superpos = utils.load_sel(args.sel_superpos, u, indexes)
 
+    if (args.unwrap or args.wrap or args.precenter):
+        unwrap_sel = utils.load_sel(args.sel_unwrap, u, indexes)
     if (args.unwrap or args.wrap):
         # Preparing molecule unwrapper
         bonded_struct = mda.Universe(
-            "epoch01/rep01/mdrun.tpr", str(args.cfg.initial_struct[0]))
-        unwrap_sel = utils.load_sel(args.sel_unwrap, u, indexes)
+            "epoch01/rep01/mdrun.tpr",
+            args.cfg.initial_struct[0]
+        )
         unwrap_sel = bonded_struct.atoms[unwrap_sel.indices]
+
+    traj_transforms = []
+    if (args.precenter):
+        if (args.precenter_atom is not None):
+            centre_atom_group = utils.load_sel(args.precenter_atom,
+                                               u, indexes)
+            if (len(centre_atom_group) != 1):
+                raise WrongSelectionSizeError(f"Selection precenter_atom={repr(args.precenter_atom)} resulted "
+                                              f"in {len(centre_atom_group)} atoms. Should be exactly 1!")
+            centre_atom = centre_atom_group[0]
+        else:
+            centre_atom = None
+        traj_transforms.append(
+            transformations.Precenter(unwrap_sel, centre_atom))
+        ca = u.atoms[traj_transforms[-1].centre_atom]
+        print(f"Precentering using atom index {ca.index}",
+              f"({ca.name} of {ca.resname}:{ca.resid})")
+
+    if (args.unwrap):
         print("Selected %d atoms for unwrapping" % len(unwrap_sel))
         if (args.unwrap_starters is None):
             unwrap_starters = []
@@ -57,9 +81,6 @@ def load_struct(args: argparse.Namespace) -> Tuple[mda.Universe, AtomGroup, List
                 args.unwrap_starters, unwrap_sel, indexes)
             print("Selected %d atoms as unwrap starters" %
                   len(unwrap_starters))
-
-    traj_transforms = []
-    if (args.unwrap):
         traj_transforms.append(
             transformations.Unwrapper(unwrap_sel, unwrap_starters))
 

src/analysis/parsers.py

@@ -1,18 +1,13 @@
 import argparse
 import pathlib
 
-from ..exceptions import NoNetworkxError, handle_errors
+from ..exceptions import handle_errors
 from ..cmd_line import config_importer
 
 
 @handle_errors
 def ancestry_graphing(args: argparse.Namespace):
-    try:
-        from .ancestry import ancestry
-    except ImportError as e:
-        raise NoNetworkxError("Failed to import functional_sampling_tool.analysis.extract.\n"
-                              "Make sure you have networkx installed before making ancestry graphs.\n"
-                              f"Original error message: {e}")
+    from .ancestry import ancestry
     ancestry(args)
 
 
@@ -62,6 +57,12 @@ def extract_subparser(subparsers: "argparse._SubParsersAction[argparse.ArgumentP
     extr_parser.add_argument("--beginning", "-b", metavar="N",
                              help="Ignore this many frames from the beginning of each simulation. By default \"ignore_from_start\" from config.",
                              default=None,  type=int)
+    extr_parser.add_argument("--precenter",   action="store_true",
+                             help="Precenter the --sel-unwrap by shifting it such that --precenter-atom is in the centre of the box  (default: No)")
+    extr_parser.add_argument("--precenter-atom", metavar="str",
+                             help="The selection to unwrap. Same syntax as --selection, but should result in a single atom being selected. "
+                             "By default None, which results in using whichever atom in --sel-unwrap is closest to box center initially.",
+                             default=None)
     extr_parser.add_argument("--unwrap",   action="store_true",
                              help="Unwrap molecules (default: No)")
     extr_parser.add_argument("--sel-unwrap", metavar="str",

src/default_config.py

@@ -33,6 +33,10 @@
 
 ############################## Trajectrory transformations #####################
 
+# whether to centre by a single atom before making anything whole
+precenter = False
+# The atom to use for precentering. By default use the one closest to box centre.
+precenter_atom = None
 # whether to make molecules whole from being broken over pbc
 unwrap_mols = False
 # The selection to make whole

src/epoch_starting.py

@@ -32,7 +32,7 @@ def init_rep(i: int,
             f"{d} already exists, will not try to overwrite. Try running \"fst clean\" to remove the latest epoch directory.")
     d.mkdir()
 
-    atoms.write(str(d / "start.gro"))
+    atoms.write(d / "start.gro")
 
     # Make a note of the origin
     with (d / "origin.txt").open("w") as f:
@@ -91,10 +91,10 @@ def next_rep(i: int,
     d_old = pathlib.Path("epoch%02d" % oldepoch) / ("rep%02d" % rep)
     next_xtc = d_old / "mdrun.xtc"
     if (next_xtc != cfg.struct.filename):
-        cfg.struct.load_new(str(next_xtc))
+        cfg.struct.load_new(next_xtc)
     cfg.struct.trajectory[frm]
 
-    cfg.struct.atoms.write(str(d / "start.gro"))
+    cfg.struct.atoms.write(d / "start.gro")
     print("Wrote structure to %s, starting to grompp..." % (d / "start.gro"))
 
     # Make a note of the origin
@@ -159,7 +159,7 @@ def start_epoch(nextepoch: int, cfg: Any,
         res = []
         if (nextepoch == 1):
             # Initial structures and first epoch
-            struct = mda.Universe(*[str(f) for f in cfg.initial_struct])
+            struct = mda.Universe(*cfg.initial_struct)
             num_frames = len(struct.trajectory)
             print("Found %d starting structures" % num_frames)
             if (num_frames > cfg.N):

src/exceptions.py

@@ -7,8 +7,8 @@
 
 class FSTException(Exception):
     """
-    A base exception differnetiate between fst-internal exception,
-    to be caught nicely, and some other code's exceptions, which should not be caught.
+    A base exception to differentiate between fst-internal exception (to be caught
+    nicely) and some other code's exceptions (which should not be caught).
     """
 
     def __init__(self, *args: Any, code=1) -> None:
@@ -17,6 +17,18 @@ def __init__(self, *args: Any, code=1) -> None:
 
 
 def handle_errors(func: Callable[[argparse.Namespace], None]) -> Callable[[argparse.Namespace], None]:
+    """
+    A function decorator that will handle FSTException and KeyboardInterrupt nicely
+    while not catching any other Exceptions.
+
+    parameters:
+        func: Callable(Namespace) -> None
+            The function to wrap with the error handling.
+
+    return:
+        wrapped_func: Callable(Namespace) -> None
+            The wrapped function.
+    """
     def wrapped_func(args: argparse.Namespace) -> None:
         try:
             return func(args)
@@ -62,6 +74,13 @@ class NotEnoughDataError(FSTException, ValueError):
     pass
 
 
+class WrongSelectionSizeError(FSTException, ValueError):
+    """
+    An exception raised, when a selection should be a certain size, but is not.
+    """
+    pass
+
+
 class NoEpochsFoundError(FSTException, FileNotFoundError):
     """
     An exception raised, when no Epoch data is found, even though the command requires it.
@@ -92,29 +111,22 @@ class NonzeroReturnError(FSTException, RuntimeError):
     pass
 
 
-class NoConfigError(FSTException, FileNotFoundError):
-    """
-    An exception raised, when the config file does not exists.
-    """
-    pass
-
-
 class RequiredFileMissingError(FSTException, FileNotFoundError):
     """
-    An exception raised, when a erquired file does not exists.
+    An exception raised, when a required file does not exists.
     """
     pass
 
 
-class NoNetworkxError(FSTException, ImportError):
+class NoConfigError(RequiredFileMissingError):
     """
     An exception raised, when the config file does not exists.
     """
     pass
 
 
-class NoSbatchLaunchError(FSTException, FileNotFoundError):
+class NoSbatchLaunchError(RequiredFileMissingError):
     """
-    An exception raised, when the config file does not exists.
+    An exception raised, when the sbatch file does not exists.
     """
     pass