diff --git a/.travis.yml b/.travis.yml index 2a15f2a2..73fe3adb 100644 --- a/.travis.yml +++ b/.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 diff --git a/docs/source/center.rst b/docs/source/center.rst index 534e4ef6..4ecc4f16 100644 --- a/docs/source/center.rst +++ b/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: diff --git a/pytim/__init__.py b/pytim/__init__.py index 4bf09bec..47862dcc 100644 --- a/pytim/__init__.py +++ b/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. diff --git a/pytim/examples/example_gitim_flat.py b/pytim/examples/example_gitim_flat.py new file mode 100644 index 00000000..d8e6783b --- /dev/null +++ b/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') + + diff --git a/pytim/examples/example_micelle.py b/pytim/examples/example_micelle.py index 16c97069..082fd459 100644 --- a/pytim/examples/example_micelle.py +++ b/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) diff --git a/pytim/gitim.py b/pytim/gitim.py index 5d8dca7e..722d614a 100755 --- a/pytim/gitim.py +++ b/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=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=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 diff --git a/pytim/itim.py b/pytim/itim.py index 6b9d84ac..306331ae 100755 --- a/pytim/itim.py +++ b/pytim/itim.py @@ -18,16 +18,19 @@ class ITIM(pytim.PYTIM): """ Identifies the interfacial molecules at macroscopically flat interfaces. - :param Universe universe: the MDAnalysis universe - :param float mesh: the grid spacing used for the testlines - :param float alpha: the probe sphere radius - :param str normal: the macroscopic interface normal direction 'x','y' or 'z' (by default is 'guess') - :param AtomGroup itim_group: identify the interfacial molecules from this group - :param dict radii_dict: dictionary with the atomic radii of the elements in the itim_group. - If None is supplied, the default one (from GROMOS 43a1) will be used. - :param int max_layers: the number of layers to be identified - :param bool info: print additional info - :param bool multiproc: parallel version (default: True. Switch off for debugging) + :param Universe universe: the MDAnalysis universe + :param float mesh: the grid spacing used for the testlines + :param float alpha: the probe sphere radius + :param str normal: the macroscopic interface normal direction 'x','y' or 'z' (by default is 'guess') + :param AtomGroup itim_group: identify the interfacial molecules from this group + :param dict radii_dict: dictionary with the atomic radii of the elements in the itim_group. + If None is supplied, the default one (from MDAnalysis) will be used. + :param int max_layers: the number of layers to be identified + :param float cluster_cut: cutoff used for neighbors or density-based cluster search (default: None disables the cluster analysis) + :param float cluster_threshold_density: Number density threshold for the density-based cluster search. 'auto' determines the threshold automatically. Default: None uses simple neighbors cluster search, if cluster_cut is not None + :param bool molecular: Switches between search of interfacial molecules / atoms (default: True) + :param bool info: print additional info + :param bool multiproc: parallel version (default: True. Switch off for debugging) Example: @@ -102,6 +105,7 @@ def __init__(self,universe,mesh=0.4,alpha=2.0,normal='guess',itim_group=None,rad self._define_groups() self._assign_normal(normal) + self.assign_radii(radii_dict) self._sanity_checks() @@ -112,14 +116,6 @@ def __init__(self,universe,mesh=0.4,alpha=2.0,normal='guess',itim_group=None,rad self._assign_layers() - def _assign_normal(self,normal): - - 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 _assign_mesh(self): """ determine a mesh size for the testlines that is compatible with the simulation box @@ -259,7 +255,7 @@ def _assign_layers(self): dtype=int); # 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()