PBC added to GITIM, centering improved, docs updated

.travis.yml

@@ -1,5 +1,9 @@
 language: python
 
+branches:
+  except:
+  - dev
+
 matrix:
     include:
         - os: linux
@@ -14,7 +18,7 @@ before_install:
     - if [[ "$TRAVIS_OS_NAME" == "osx"   ]] ; then pip install --upgrade pip               ; fi
     - if [[ "$TRAVIS_OS_NAME" == "linux" ]] ; then sudo apt-get install liblapack-dev -y   ; fi
 
-install: 
+install:
 - pip install Sphinx>=1.4.3
 - pip install matplotlib
 - pip install cython

docs/source/center.rst

@@ -30,37 +30,40 @@ If you want to save a local copy of the file, you just need to use this code fra
     u         = mda.Universe(pytim.datafiles.WATER_GRO)
     u.atoms.write('centering1.pdb')
 
-**Default, centering at the origin**
 
-When initializing/computing the interface, Pytim internally centers the interface in different ways, depending whether :class:`~pytim.itim.ITIM` or :class:`~pytim.gitim.GITIM` is used. For planar interfaces, by default the middle of the liquid phase is moved to the origin along the surface normal, while the positions along other directions are placed in the box. In our case, the result is:
+When initializing/computing the interface, Pytim internally centers the interface in different ways, depending whether :class:`~pytim.itim.ITIM` or :class:`~pytim.gitim.GITIM` is used. When the configuration with the surface layers information is written to file using :meth:`~pytim.itim.writepdb`, the system can be shifted so that the liquid phase is placed at different positions.
+
+**Default, no centering**
+
+If the (default) option `centered='no'` is used, then all atomic positions are kept the same as in the input file.
+This is like the initial case, however with information on the surface molecules added to the pdb.
 
 .. code-block:: python
 
-    interface = pytim.ITIM(u)
-    interface.writepdb('centering2.pdb')
+    interface.writepdb('centering3.pdb',centered='no')
 
 
-.. image:: centering2.png
+
+.. image:: centering3.png
    :width: 35%
    :align: center
 
+**Centering at the origin**
+
+The system can be shifted so that the liquid phase is placed across the origin of the normal axis, using the option `centered=origin`.
 This comes handy to quickly discriminate between upper and lower surface atoms, to spot immediately the middle of the liquid phase, or to verify that the automatic centering method behaved as expected.
 
-**No centering**
-
-If the option `centered='no'` is passed to :meth:`~pytim.itim.ITIM..writepdb`, then the first atom of the system is kept at its original position (i.e., no shift is applied) although the system is always put into the basic cell. 
-
 .. code-block:: python
 
-    interface.writepdb('centering3.pdb',centered='no')
-
+    interface = pytim.ITIM(u)
+    interface.writepdb('centering2.pdb',centered='origin')
 
-This is like the initial case, however with information on the surface molecules added to the pdb.
 
-.. image:: centering3.png
+.. image:: centering2.png
    :width: 35%
    :align: center
 
+
 **Centering in the middle**
 
 If the option `centered='middle'` is passed  to :meth:`~pytim.itim.ITIM.writepdb`, instead, then the liquid part of the system is placed in the middle of the box along the surface normal:

pytim/__init__.py

@@ -33,7 +33,7 @@ class PYTIM(object):
     __metaclass__ = ABCMeta
 
     directions_dict={0:'x',1:'y',2:'z','x':'x','y':'y','z':'z','X':'x','Y':'y','Z:':'z'}
-    symmetry_dict={'cylindrical':'cylindrical','spherical':'spherical'}
+    symmetry_dict={'cylindrical':'cylindrical','spherical':'spherical','planar':'planar'}
 
     ALPHA_NEGATIVE = "parameter alpha must be positive"
     ALPHA_LARGE= "parameter alpha must be smaller than the smaller box side"
@@ -50,7 +50,7 @@ class PYTIM(object):
     WRONG_DIRECTION="Wrong direction supplied. Use 'x','y','z' , 'X', 'Y', 'Z' or 0, 1, 2"
     CENTERING_FAILURE="Cannot center the group in the box. Wrong direction supplied?"
 
-    def writepdb(self,filename='layers.pdb',centered='origin',multiframe=True):
+    def writepdb(self,filename='layers.pdb',centered='no',multiframe=True):
         """ Write the frame to a pdb file, marking the atoms belonging
             to the layers with different beta factor.
 
@@ -70,41 +70,39 @@ def writepdb(self,filename='layers.pdb',centered='origin',multiframe=True):
             >>> interface.writepdb('layers.pdb',centered='no')
 
         """
-        center_options=['no','middle','origin']
+        center_options=['no','middle','origin',False,True]
         if centered not in center_options:
-            centered='origin'
+            centered='no'
+        if centered == False: centered='no'
+        if centered == True : centered='middle' 
         try:
             if centered=='no':
-                original_pos = self.universe.atoms.positions[:]
-                translation = self.reference_position - self.universe.atoms[0].position[:]
-                self.universe.atoms.translate(translation)
-                self.universe.atoms.pack_into_box(self.universe.dimensions[:3])
-
-            if centered=='middle' and self.symmetry=='planar':
-                original_pos = self.universe.atoms.positions[:]
-                translation = [0,0,0]
-                translation[self.normal] = self.universe.dimensions[self.normal]/2.
-                self.universe.atoms.translate(translation)
-                self.universe.atoms.pack_into_box(self.universe.dimensions[:3])
-
+                self.universe.atoms.positions=self.original_positions
+
+            if centered=='middle': 
+                # NOTE: this assumes that all method relying on 'planar' symmetry must center the interface along the normal
+                if self.symmetry=='planar':
+                    translation = [0,0,0]
+                    translation[self.normal] = self.universe.dimensions[self.normal]/2.
+                    self.universe.atoms.positions+=np.array(translation)
+                    self.universe.atoms.pack_into_box(self.universe.dimensions[:3])
+
             PDB=MDAnalysis.Writer(filename, multiframe=True, bonds=False,
                             n_atoms=self.universe.atoms.n_atoms)
-
+    
             PDB.write(self.universe.atoms)
-
-            if centered=='no' or centered=='middle':
-                self.universe.atoms.positions=original_pos
+
         except:
             print("Error writing pdb file")
 
-    def savepdb(self,filename='layers.pdb',centered=True,multiframe=True):
+    def savepdb(self,filename='layers.pdb',centered='no',multiframe=True):
         """ An alias to :func:`writepdb`
         """
         self.writepdb(filename,centered,multiframe)
 
     def assign_radii(self,radii_dict):
         try:
-            _groups = self.extra_cluster_groups[:] # deep copy
+            _groups = np.copy(self.extra_cluster_groups[:])
         except:
             _groups = []
         _groups.append(self.itim_group)
@@ -134,12 +132,22 @@ def assign_radii(self,radii_dict):
                         print("!! Pass a dictionary of radii (in Angstrom) with the option radii_dict")
                         print("!! for example: r={'"+_atype+"':1.2,...} ; inter=pytim.ITIM(u,radii_dict=r)")
 
-                _g.radii=_radii[:] #deep copy
+                _g.radii=np.copy(_radii[:])
 
                 assert not np.any(np.equal(_g.radii,None)) , self.UNDEFINED_RADIUS
                 del _radii
                 del _types
 
+    def _assign_normal(self,normal):
+        assert self.symmetry=='planar', "Error: wrong symmetry for normal assignement"
+        assert self.itim_group is not None, self.UNDEFINED_ITIM_GROUP
+        if normal=='guess':
+            self.normal=utilities.guess_normal(self.universe,self.itim_group)
+        else:
+            assert normal in self.directions_dict, self.WRONG_DIRECTION
+            self.normal = self.directions_dict[normal]
+
+
     def center(self, group, direction=None, halfbox_shift=True):
         """
         Centers the liquid slab in the simulation box.

pytim/examples/example_gitim_flat.py

@@ -0,0 +1,18 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
+import MDAnalysis as mda
+import pytim
+from   pytim.datafiles import *
+
+u         = mda.Universe(WATER_GRO)
+
+radii     = pytim_data.vdwradii(G43A1_TOP)
+
+interface = pytim.GITIM(u,molecular=True,symmetry='planar',alpha=2.5)
+
+layers     = interface.layers[:]
+print layers
+
+interface.writepdb('gitim_flat.pdb',centered='middle')
+
+

pytim/examples/example_micelle.py

@@ -11,7 +11,7 @@
 
 interface = pytim.GITIM(u,itim_group=g,molecular=False,symmetry='spherical',alpha=2.5,)
 
-layer     = interface.layers(1)
+layer     = interface.layers[0]
 
 interface.writepdb('gitim.pdb',centered=False)
 

pytim/gitim.py

@@ -11,6 +11,7 @@
 from   scipy.spatial     import Delaunay
 from   scipy.spatial     import distance
 from   scipy.interpolate import LinearNDInterpolator
+import itertools 
 from __builtin__         import zip as builtin_zip
 from   pytim             import utilities
 import pytim
@@ -46,8 +47,8 @@ class GITIM(pytim.PYTIM):
 
     """
 
-    def __init__(self,universe,alpha=2.0,symmetry='spherical',itim_group=None,radii_dict=None,
-                 max_layers=1,cluster_cut=None,molecular=True,extra_cluster_groups=None,
+    def __init__(self,universe,alpha=2.0,symmetry='spherical',normal='guess',itim_group=None,radii_dict=None,
+                 max_layers=1,cluster_cut=None,cluster_threshold_density=None,molecular=True,extra_cluster_groups=None,
                  info=False,multiproc=True):
 
         #TODO add type checking for each MDA class passed
@@ -56,6 +57,7 @@ def __init__(self,universe,alpha=2.0,symmetry='spherical',itim_group=None,radii_
         pytim.PatchTrajectory(universe.trajectory,self)
 
         self.universe=universe
+        self.cluster_threshold_density = cluster_threshold_density
         self.alpha=alpha
         self.max_layers=max_layers
         self._layers=np.empty([max_layers],dtype=type(universe.atoms))
@@ -74,6 +76,10 @@ def __init__(self,universe,alpha=2.0,symmetry='spherical',itim_group=None,radii_
         self._define_groups()
 
         self._assign_symmetry(symmetry)
+
+        if(self.symmetry=='planar'):
+            self._assign_normal(normal)
+
         self.assign_radii(radii_dict)
         self._sanity_checks()
 
@@ -138,71 +144,56 @@ def circumradius(self,simplex):
         positiveR = R[R>=0]
         return np.min(positiveR) if positiveR.size == 1 else 0
 
-##       check configuration  ##
-##            print r_i
-##            print R_i
-##            from mpl_toolkits.mplot3d import Axes3D
-##            import matplotlib.pyplot as plt
-##
-##            def drawSphere(xCenter, yCenter, zCenter, r):
-##                #draw sphere
-##                u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
-##                x=np.cos(u)*np.sin(v)
-##                y=np.sin(u)*np.sin(v)
-##                z=np.cos(v)
-##                # shift and scale sphere
-##                x = r*x + xCenter
-##                y = r*y + yCenter
-##                z = r*z + zCenter
-##                return (x,y,z)
-##
-##            fig = plt.figure()
-##            ax = fig.add_subplot(111, projection='3d')
-##
-##            # draw a sphere for each data point
-##            for (xi,yi,zi,ri) in zip(r_i[:,0],r_i[:,1],r_i[:,2],R_i):
-##                (xs,ys,zs) = drawSphere(xi,yi,zi,ri)
-##                ax.plot_wireframe(xs, ys, zs, color="r")
-##            plt.show()
-##            exit()
-
-
-
     def alpha_shape(self,alpha):
+        #print  utilities.lap()
         box    = self.universe.dimensions[:3]
         delta  = 2.*self.alpha+1e-6
         points = self.cluster_group.positions[:]
-        extrapoints  = []
+        extrapoints = np.copy(points)
+        nrealpoints = len(points)
+        extraids    = np.arange(len(points),dtype=np.int)
+        tmpPoints   = []
+        tmpIDs      = []
+        for shift in np.array(list(itertools.product([1,-1,0],repeat=3))):
+            if(np.sum(shift*shift)): # avoid [0,0,0]
+                # this needs some explanation:
+                # if shift ==0  -> the condition is always true
+                # if shift ==1  -> the condition is x > box - delta
+                # if shift ==-1 -> the condition is -x > 0 - delta -> x < delta
+                # Requiring np.all() to be true makes the logical and returns (axis=1) True for all indices whose atoms satisfy the condition
+                selection = np.all(shift * points >= shift*shift*( (box + shift * box)/2. - delta)   ,axis=1)
+                # add the new points at the border of the box
+                extrapoints=np.append(extrapoints,points[selection]-shift*box,axis=0)
+                # we keep track of the original ids.
+                extraids=np.append(extraids,np.where(selection)[0])
+
         # add points at the vertices of the expanded (by 2 alpha) box
         for dim in range(8):
                 # [0,0,0],[0,0,1],[0,1,0],...,[1,1,1]
                 tmp = np.array(np.array(list(np.binary_repr(dim,width=3)),dtype=np.int8),dtype=np.float)
                 tmp *=(box+delta)
                 tmp[tmp<box/2.]-=delta
-                extrapoints.append(tmp)
-        points = np.append(points,extrapoints,axis=0)
-
-#       time=dict();utilities.lap()
+                tmp=np.reshape(tmp,(1,3))
+                extrapoints=np.append(extrapoints,tmp,axis=0)
+                extraids=np.append(extraids,-1)
 
-        self.triangulation = Delaunay(points)
-        self.triangulation.radii = np.append(self.cluster_group.radii[:],np.zeros(8))
-
-#       time['triangulate']=utilities.lap()
+        #print utilities.lap()
+        self.triangulation = Delaunay(extrapoints)
+        self.triangulation.radii = np.append(self.cluster_group.radii[extraids[extraids>=0]],np.zeros(8))
+        #print utilities.lap()
 
         prefiltered = self.alpha_prefilter(self.triangulation, alpha)
-#       time['prefilter']=utilities.lap()
 
-        alpha_shape = prefiltered [ [ self.circumradius(simplex) >=self.alpha for simplex in prefiltered      ]   ]
-#       time['alpha']=utilities.lap()
-#       print time
-#       print len(self.triangulation.simplices),len(prefiltered),len(alpha_shape),len(_ids)
+        #print utilities.lap()
+        a_shape = prefiltered [ [ self.circumradius(simplex) >=self.alpha for simplex in prefiltered      ]   ]
 
-        _ids = np.unique(alpha_shape)
+        #print utilities.lap()
+        _ids = np.unique(a_shape.flatten())
 
-        # remove the indices corresponding to the 8 additional points
-        #ids =_ids
-        ids = _ids[_ids<len(points)-8]
+        # remove the indices corresponding to the 8 additional points, which have extraid==-1
+        ids = _ids[np.logical_and(_ids>=0 , _ids < nrealpoints)]
 
+        #print utilities.lap()
         return ids
 
     def _assign_layers(self):
@@ -211,12 +202,11 @@ def _assign_layers(self):
 
         """
         # this can be used later to shift back to the original shift
-        self.reference_position=self.universe.atoms[0].position[:]
-
+        self.original_positions=np.copy(self.universe.atoms.positions[:])
         self.universe.atoms.pack_into_box()
 
         if(self.cluster_cut is not None): # groups have been checked already in _sanity_checks()
-            labels,counts,n_neigh = utilities.do_cluster_analysis_DBSCAN(self.itim_group,self.cluster_cut[0],self.universe.dimensions[:6],self.molecular)
+            labels,counts,n_neigh = utilities.do_cluster_analysis_DBSCAN(self.itim_group,self.cluster_cut[0],self.universe.dimensions[:6],self.cluster_threshold_density,self.molecular)
             labels = np.array(labels)
             label_max = np.argmax(counts) # the label of atoms in the largest cluster
             ids_max   = np.where(labels==label_max)[0]  # the indices (within the group) of the
@@ -225,10 +215,18 @@ def _assign_layers(self):
 
         else:
             self.cluster_group=self.itim_group ;
+
+        if self.symmetry=='planar':
+            utilities.centerbox(self.universe,center_direction=self.normal)
+            self.center(self.cluster_group,self.normal)
+            utilities.centerbox(self.universe,center_direction=self.normal)
         if self.symmetry=='spherical':
             self.center(self.cluster_group,'x',halfbox_shift=False)
             self.center(self.cluster_group,'y',halfbox_shift=False)
             self.center(self.cluster_group,'z',halfbox_shift=False)
+            self.universe.atoms.pack_into_box(self.universe.dimensions[:3])
+
+
 
         # first we label all atoms in itim_group to be in the gas phase
         self.itim_group.atoms.bfactors = 0.5
@@ -242,7 +240,10 @@ def _assign_layers(self):
         alpha_ids = self.alpha_shape(self.alpha)
 
         # only the 1st layer is implemented in gitim so far
-        self._layers[0] = self.cluster_group[alpha_ids]
+        if self.molecular == True:
+            self._layers[0] = self.cluster_group[alpha_ids].residues.atoms
+        else:
+            self._layers[0] = self.cluster_group[alpha_ids]
 
         for layer in self._layers:
                 layer.bfactors = 1