SCHWED-1740: document multi-loop in MmSystemCreator

doc/tests/scripts/loop_mm_sys_creation.py

@@ -12,24 +12,20 @@ num_residues = len(res_list)
 
 # get all atom positions for full protein
 all_atoms = loop.AllAtomPositions(res_list)
-# get res. indices for loop followed by rest
-# -> here we add everything to rest but in practice this could be
-#    just a subset of residues relevant to the loop of interest
-loop_start = 10
-loop_length = 6
-res_indices = list(range(loop_start, loop_start + loop_length)) \
-            + list(range(0, loop_start)) \
-            + list(range(loop_start + loop_length, num_residues))
+# here full structure in res_indices but in practice this could
+# be just a subset of residues relevant to the loop of interest
+res_indices = range(num_residues)
+# define two loops (indices into res_indices list)
+loop_start_indices = [10, 20]
+loop_lengths = [6, 4]
 # define which of the res_indices is terminal
-is_n_ter = [False] * num_residues
-is_c_ter = [False] * num_residues
-is_n_ter[res_indices.index(0)] = True
-is_c_ter[res_indices.index(num_residues-1)] = True
+is_n_ter = [True] + [False] * (num_residues - 1)
+is_c_ter = [False] * (num_residues - 1) + [True]
 # get disulfid bridges
 disulfid_bridges = mm_sys.GetDisulfidBridges(all_atoms, res_indices)
 # setup MM system
-mm_sys.SetupSystem(all_atoms, res_indices, loop_length,
-                   is_n_ter, is_c_ter, disulfid_bridges)
+mm_sys.SetupSystem(all_atoms, res_indices, loop_start_indices,
+                   loop_lengths, is_n_ter, is_c_ter, disulfid_bridges)
 
 # run simulation
 sim = mm_sys.GetSimulation()

loop/doc/mm_system_creation.rst

@@ -11,6 +11,15 @@ The example below showcases the creation and use of an MM system:
 
 .. literalinclude:: ../../../tests/doc/scripts/loop_mm_sys_creation.py
 
+.. note::
+
+  To simplify use when sidechains may be missing and the region of interest for
+  the MM system has to be determined, it might be better to use the simpler
+  :class:`promod3.modelling.SidechainReconstructor` and
+  :class:`promod3.modelling.AllAtomRelaxer` classes. Even if all the sidechains
+  are available, those classes will be helpful since the overhead to check
+  sidechains without reconstructing them is minimal.
+
 
 Create MM systems for loops
 --------------------------------------------------------------------------------
@@ -61,6 +70,8 @@ Create MM systems for loops
 
   .. method:: SetupSystem(all_pos, res_indices, loop_length, is_n_ter, \
                           is_c_ter, disulfid_bridges)
+              SetupSystem(all_pos, res_indices, loop_start_indices, \
+                          loop_lengths, is_n_ter, is_c_ter, disulfid_bridges)
 
     Setup a new system for the loop / surrounding defined by *all_pos* and
     *res_indices*. Positions are taken from *all_pos[res_indices[i]]* and each
@@ -72,11 +83,16 @@ Create MM systems for loops
     approximation. Similarly, we do not generate OXT atoms for residues followed
     by a gap unless they are specifically labelled as C-terminal.
 
-    Residues *res_indices[i]* for *i* < *loop_length* are loop, rest is
-    surrounding. Residues *res_indices[0]* / *res_indices[loop_length-1]* are
-    assumed to be N- / C-stems of loop. The stems of the loop are kept fixed (N,
+    Each loop is defined by an entry in *loop_start_indices* and *loop_lengths*.
+    For loop *i_loop*, *res_indices[loop_start_indices[i_loop]]* is the N-stem
+    and *res_indices[loop_start_indices[i_loop] + loop_lengths[i_loop] - 1]* is
+    the C-stem of the loop. The loop indices are expected to be contiguous in
+    *res_indices* between the stems. The stems of the loop are kept fixed (N,
     CA, CB for N-stem / CA, CB, C, O for C-stem) unless they are terminal
-    residues.
+    residues. Overlapping loops are merged and 0-length loops are removed.
+
+    If *loop_length* is given, a single loop with start index 0 and the given
+    loop length is defined.
 
     If a system was setup previously, it is overwritten completely.
 
@@ -88,8 +104,12 @@ Create MM systems for loops
     :type all_pos:  :class:`AllAtomPositions`
     :param res_indices: Residue indices into *all_pos*.
     :type res_indices:  :class:`list` of :class:`int`
-    :param loop_length: Length of loop.
+    :param loop_length: Length of single loop (incl. stems).
     :type loop_length:  :class:`int`
+    :param loop_start_indices: Start indices of loops.
+    :type loop_start_indices:  :class:`list` of :class:`int`
+    :param loop_lengths: Lengths of loops (incl. stems).
+    :type loop_lengths:  :class:`list` of :class:`int`
     :param is_n_ter: For each *res_indices[i]*, *is_n_ter[i]* defines whether
                      that residue is N-terminal.
     :type is_n_ter:  :class:`list` of :class:`bool`
@@ -102,10 +122,11 @@ Create MM systems for loops
     :type disulfid_bridges:  :class:`list` of :class:`tuple` with two
                              :class:`int`
 
-    :raises: :exc:`~exceptions.RuntimeError` if *loop_length* is longer than
-             *res_indices*, *res_indices* / *is_n_ter* / *is_c_ter* have
-             inconsistent lengths or if any *all_pos[res_indices[i]]* is invalid
-             or has unset heavy atoms.
+    :raises: :exc:`~exceptions.RuntimeError` if loops out of bounds with respect
+             to *res_indices*, all loops have length 0, *loop_start_indices* /
+             *loop_lengths* or *res_indices* / *is_n_ter* / *is_c_ter* have
+             inconsistent lengths, or if any *all_pos[res_indices[i]]* is
+             invalid or has unset heavy atoms.
 
   .. method:: UpdatePositions(all_pos, res_indices)
 
@@ -130,14 +151,21 @@ Create MM systems for loops
     :meth:`SetupSystem`. This function then extracts these positions back into a
     compatible output storage.
 
-    :param loop_pos: Storage for loop positions. Residues *loop_pos[i]* and
-                     *all_pos[res_indices[i]]* (of :meth:`SetupSystem`) must
-                     have matching amino acid types for *i* < *loop_length*.
+    If *loop_pos* is passed, only the loops are stored. The loop residues are
+    stored contiguously and the loops are stored in the order given by
+    :meth:`GetLoopStartIndices` / :meth:`GetLoopLengths`. The *loop_pos* storage
+    must have at least :meth:`GetNumLoopResidues` residues and must have
+    matching amino acid types with respect to the loop residues passed as
+    *all_pos* of :meth:`SetupSystem`.
+
+    If *out_pos* and *res_indices* are passed, residues must have matching amino
+    acid types for *out_pos[res_indices[i]]* and *all_pos[res_indices[i]]* of
+    :meth:`SetupSystem`. Only loop residues with *i* in the ranges defined by
+    :meth:`GetLoopStartIndices` / :meth:`GetLoopLengths` are touched.
+
+    :param loop_pos: Reduced storage only covering loop positions.
     :type loop_pos:  :class:`AllAtomPositions`
-    :param out_pos: Storage for loop positions. Residues
-                    *out_pos[res_indices[i]]* and *all_pos[res_indices[i]]* (of
-                    :meth:`SetupSystem`) must have matching amino acid types for
-                    *i* < *loop_length*.
+    :param out_pos: Storage for loop positions linked to *res_indices*.
     :type out_pos:  :class:`AllAtomPositions`
     :param res_indices: Residue indices into *out_pos*.
     :type res_indices:  :class:`list` of :class:`int`
@@ -151,6 +179,26 @@ Create MM systems for loops
              MM simulations.
     :rtype:  :class:`ost.mol.mm.Simulation`
 
+  .. method:: GetNumResidues
+
+    :return: Number of residues of current simulation setup.
+    :rtype:  :class:`int`
+
+  .. method:: GetNumLoopResidues
+
+    :return: Number of loop residues (incl. stems) of current simulation setup.
+    :rtype:  :class:`int`
+
+  .. method:: GetLoopStartIndices
+
+    :return: Start indices of loops (see :meth:`SetupSystem`).
+    :rtype:  :class:`list` of :class:`int`
+
+  .. method:: GetLoopLengths
+
+    :return: Lengths of loops (see :meth:`SetupSystem`).
+    :rtype:  :class:`list` of :class:`int`
+
   .. method:: GetForcefieldAminoAcids
 
     :return: Forcefield-specific amino acid type for each residue of last

loop/src/mm_system_creator.cc

@@ -197,7 +197,8 @@ MmSystemCreator::MmSystemCreator(ForcefieldLookupPtr ff_lookup,
                                  bool kill_electrostatics,
                                  Real nonbonded_cutoff,
                                  bool inaccurate_pot_energy)
-                                 : aa_lookup_(AminoAcidLookup::GetInstance()) {
+                                 : num_loop_residues_(0)
+                                 , aa_lookup_(AminoAcidLookup::GetInstance()) {
   // store settings and check for CPU support
   ff_lookup_ = ff_lookup;
   fix_surrounding_hydrogens_ = fix_surrounding_hydrogens;