zero radis handling

pytim/gitim.py

@@ -68,6 +68,10 @@ class GITIM(Interface):
                                     of the interface: 'generic' (default)
                                     or  'planar'
         :param bool centered:       Center the  :py:obj:`group`
+        :param bool include_zero_radius: if false (default) exclude atoms with zero radius
+                                    from the surface analysis (they are always included
+                                    in the cluster search, if present in the relevant
+                                    group) to avoid some artefacts.
         :param bool info:           Print additional info
         :param bool warnings:       Print warnings
         :param bool autoassign:     If true (default) detect the interface
@@ -124,6 +128,7 @@ def __init__(self,
                  max_layers=1,
                  radii_dict=None,
                  cluster_cut=None,
+                 include_zero_radius=False,
                  cluster_threshold_density=None,
                  extra_cluster_groups=None,
                  biggest_cluster_only=False,
@@ -155,6 +160,7 @@ def __init__(self,
         self.normal = None
         self.PDB = {}
         self.molecular = molecular
+        self.include_zero_radius = include_zero_radius
         sanity.assign_cluster_params(cluster_cut,
                                      cluster_threshold_density, extra_cluster_groups)
         sanity.check_multiple_layers_options()
@@ -254,7 +260,8 @@ def _assign_layers_setup(self):
         # then all atoms in the larges group are labelled as liquid-like
         self.label_group(self.cluster_group.atoms, beta=0.0)
 
-        alpha_group = self.cluster_group[:]
+        if self.include_zero_radius: alpha_group = self.cluster_group[:]
+        else: alpha_group = self.cluster_group[self.cluster_group.radii > 0.0]
 
         # TODO the successive layers analysis should be done by removing points
         # from the triangulation and updating the circumradius of the neighbors

pytim/interface.py

@@ -73,6 +73,9 @@ class Interface(object):
     autoassign, _autoassign =\
         _create_property('autoassign',
                          "(bool) assign layers every time a frame changes")
+    include_zero_radius, _include_zero_radius=\
+        _create_property('include_zero_radius',
+                         "(bool) include atoms with zero radius in the analysis (excluded by default)")
     cluster_threshold_density, _cluster_threshold_density =\
         _create_property('cluster_threshold_density',
                          "(float) threshold for the density-based filtering")

pytim/itim.py

@@ -56,6 +56,10 @@ class ITIM(Interface):
                                   search, if cluster_cut is not None
         :param Object extra_cluster_groups: Additional groups, to allow for
                                   mixed interfaces
+        :param bool include_zero_radius: if false (default) exclude atoms with zero radius
+                                  from the surface analysis (they are always included
+                                  in the cluster search, if present in the relevant
+                                  group) to avoid some artefacts.
         :param bool info:         Print additional info
         :param bool centered:     Center the  :py:obj:`group`
         :param bool warnings:     Print warnings
@@ -186,6 +190,7 @@ def __init__(self,
                  max_layers=1,
                  radii_dict=None,
                  cluster_cut=None,
+                 include_zero_radius=False,
                  cluster_threshold_density=None,
                  extra_cluster_groups=None,
                  info=False,
@@ -214,6 +219,7 @@ def __init__(self,
         self.normal = None
         self.PDB = {}
         self.molecular = molecular
+        self.include_zero_radius = include_zero_radius
 
         sanity.assign_cluster_params(cluster_cut,
                                      cluster_threshold_density, extra_cluster_groups)
@@ -294,6 +300,8 @@ def _assign_one_side(self,
             # atom here goes to 0 to #sorted_atoms, it is not a MDAnalysis
             # index/atom
             for atom in sorted_atoms:
+                if self.include_zero_radius == False and not (_radius[atom] > 0.0):
+                    continue
                 if self._seen[uplow][atom] != 0:
                     continue
 

pytim/observables/rdf.py

@@ -316,7 +316,7 @@ def sample(self, g1=None, g2=None, kargs1=None, kargs2=None):
         self.count += count
 
         box = self.universe.dimensions
-        self.volume += np.product(box[:3])
+        self.volume += np.prod(box[:3])
         self.nsamples += 1
         self.n_squared += len(self.g1) * len(self.g2)
 

pytim/sasa.py

@@ -56,6 +56,10 @@ class SASA(GITIM):
                                     those in the largest cluster. Need to
                                     specify also a :py:obj:`cluster_cut` value.
         :param bool centered:       Center the  :py:obj:`group`
+        :param bool include_zero_radius: if false (default) exclude atoms with zero radius
+                                    from the surface analysis (they are always included
+                                    in the cluster search, if present in the relevant
+                                    group) to avoid some artefacts.
         :param bool info:           Print additional info
         :param bool warnings:       Print warnings
         :param bool autoassign:     If true (default) detect the interface

pytim/simple_interface.py

@@ -29,7 +29,10 @@ class SimpleInterface(Interface):
                                   be also passed (upper and lower)
         :param AtomGroup upper:   The upper surface group (if symmetry='planar')
         :param AtomGroup lower:   The lower surface group (if symmetry='planar')
-
+        :param bool include_zero_radius: if false (default) exclude atoms with zero radius
+                                  from the surface analysis (they are always included
+                                  in the cluster search, if present in the relevant
+                                  group) to avoid some artefacts.
 
         Example: compute an intrinsic profile by interpolating the surface
         position from the P atoms of a DPPC bilayer
@@ -82,13 +85,15 @@ def __init__(self,
                  alpha=1.5,
                  symmetry='generic',
                  normal='z',
+                 include_zero_radius=False,
                  upper=None,
                  lower=None):
 
         self.symmetry = symmetry
         self.universe = universe
         self.group = group
         self.alpha = alpha
+        self.include_zero_radius = include_zero_radius
         self.upper = upper
         self.lower = lower
         emptyg = universe.atoms[0:0]
@@ -99,6 +104,9 @@ def __init__(self,
         sanity.assign_universe(universe, group)
         sanity.assign_alpha(alpha)
         sanity.assign_radii(radii_dict=None)
+
+        if not self.include_zero_radius: self.group = self.group[self.group.radii > 0.0]
+
         if normal in [0, 1, 2]:
             self.normal = normal
         else:

pytim/willard_chandler.py

@@ -54,6 +54,11 @@ class WillardChandler(Interface):
                                   search, if cluster_cut is not None
         :param Object extra_cluster_groups: Additional groups, to allow for
                                   mixed interfaces
+
+        :param bool include_zero_radius: if false (default) exclude atoms with zero radius
+                                  from the surface analysis (they are always included
+                                  in the cluster search, if present in the relevant
+                                  group) to avoid some artefacts.
         :param bool centered:     Center the  :py:obj:`group`
         :param bool warnings:     Print warnings
 
@@ -116,6 +121,7 @@ def __init__(self,
                  mesh=2.0,
                  symmetry='spherical',
                  cluster_cut=None,
+                 include_zero_radius=False,
                  cluster_threshold_density=None,
                  extra_cluster_groups=None,
                  centered=False,
@@ -124,6 +130,7 @@ def __init__(self,
                  **kargs):
 
         self.autoassign, self.do_center = autoassign, centered
+        self.include_zero_radius = include_zero_radius
         sanity = SanityCheck(self, warnings=warnings)
         sanity.assign_universe(universe, group)
         sanity.assign_alpha(alpha)
@@ -216,7 +223,10 @@ def _assign_layers(self):
         if self.do_center is True:
             self.center()
 
-        pos = self.cluster_group.positions
+        if self.include_zero_radius:
+            pos = self.cluster_group.positions
+        else:
+            pos = self.cluster_group.positions[self.cluster_group.radii > 0.0]
         box = self.universe.dimensions[:3]
 
         ngrid, spacing = utilities.compute_compatible_mesh_params(