From 7abc2edcfac2a1d1af038d0ac86682e3a1087416 Mon Sep 17 00:00:00 2001
From: nixnmtm <nixnmtm@gmail.com>
Date: Mon, 25 Nov 2019 20:59:35 +0800
Subject: [PATCH] Multiple Fixes

  1. Bug in rmax and rmin property setter fixed
  2. Skeleton of protein bonlists are maintained when selecting bondlist
     with rmax cutoff.
---
 src/fluctmatch/commands/cmd_convert.py |  4 ++--
 src/fluctmatch/models/enm.py           | 11 +++++++----
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/fluctmatch/commands/cmd_convert.py b/src/fluctmatch/commands/cmd_convert.py
index df24e8a..b00eb62 100644
--- a/src/fluctmatch/commands/cmd_convert.py
+++ b/src/fluctmatch/commands/cmd_convert.py
@@ -230,14 +230,14 @@ def cli(
         return
 
     kwargs = dict()
+    kwargs.update(dict(rmin=rmin, rmax=rmax,))
+
     universe = modeller(topology, trajectory, com=com, model=model, **kwargs)
 
     kwargs.update(
         dict(
             outdir=outdir,
             prefix=prefix,
-            rmin=rmin,
-            rmax=rmax,
             charmm_version=charmm_version,
             extended=extended,
             resid=not resid,
diff --git a/src/fluctmatch/models/enm.py b/src/fluctmatch/models/enm.py
index 37bbf02..13b6c6e 100644
--- a/src/fluctmatch/models/enm.py
+++ b/src/fluctmatch/models/enm.py
@@ -46,9 +46,9 @@ class Enm(ModelBase):
     describe = "Elastic network model"
 
     def __init__(self, *args, **kwargs):
+        self._rmin = kwargs.pop("rmin", 0.)
+        self._rmax = kwargs.pop("rmax", 10.)
         super().__init__(*args, **kwargs)
-        self._rmin = kwargs.get("rmin", 0.)
-        self._rmax = kwargs.get("rmax", 10.)
         self._initialize(*args, **kwargs)
 
     def __repr__(self):
@@ -91,8 +91,11 @@ def _add_bonds(self):
                               (distmat <= self._rmax))
         else:
             a0, a1 = np.where((distmat > self._rmin) & (distmat <= self._rmax))
-        bonds = topologyattrs.Bonds(
-            set([(x, y) for x, y in zip(a0, a1) if y > x]))
+        # backbone skeleton must be maintained.
+        skeleton_bonds = set([(skeleton[0].ix, skeleton[1].ix) for skeleton in self.atu.bonds])
+        rmax_bonds = set([(x, y) for x, y in zip(a0, a1) if y > x])
+        all_bonds = skeleton_bonds.union(rmax_bonds)
+        bonds = topologyattrs.Bonds(all_bonds)
         self._topology.add_TopologyAttr(bonds)
         self._generate_from_topology()