diff --git a/AUTHORS.md b/AUTHORS.md new file mode 100644 index 00000000..57108d4e --- /dev/null +++ b/AUTHORS.md @@ -0,0 +1,7 @@ +Marcello-Sega https://api.github.com/users/Marcello-Sega +gyorgy-hantal https://api.github.com/users/gyorgy-hantal +balazsfabian https://api.github.com/users/balazsfabian +elija-feigl https://api.github.com/users/elija-feigl +JakobMichl https://api.github.com/users/JakobMichl +iuvenis https://api.github.com/users/iuvenis +bnovak1 https://api.github.com/users/bnovak1 diff --git a/README.md b/README.md index 371cd70f..2f429e7f 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,6 @@ So far the following interface/phase identification methods have been implemente * GITIM * SASA * Willard Chandler -* Chacon Tarazona * DBSCAN filtering ## Supported formats @@ -475,7 +474,5 @@ Depending on which algorithm you are using, you might also want to cite the foll [M. Sega and G. Hantal._Phys. Chem. Chem. Phys._ **29**, 18968-18974 (2017)](https://doi.org/10.1039/C7CP02918G) Phase and interface determination in computer simulations of liquid mixtures with high partial miscibility. -[E. Chacón, P. Tarazona, Phys. Rev. Lett. **91**, 166103 (2003)](http://dx.doi.org/10.1103/PhysRevLett.91.166103) Intrinsic profiles beyond the capillary wave theory: A Monte Carlo study. - [A. P. Willard, D. Chandler, J. Phys. Chem. B **114**,1954 (2010)](http://dx.doi.org/10.1021/jp909219k) Instantaneous Liquid Interfaces diff --git a/README.rst b/README.rst index cdf1ff56..ae4db070 100644 --- a/README.rst +++ b/README.rst @@ -12,8 +12,8 @@ So far the following methods have been implemented: * ITIM * GITIM +* SASA * Willard Chandler -* Chacon Tarazona * DBSCAN filtering ---- diff --git a/docs/source/ChaconTarazona.rst b/docs/source/ChaconTarazona.rst deleted file mode 100644 index 45d0cb12..00000000 --- a/docs/source/ChaconTarazona.rst +++ /dev/null @@ -1,16 +0,0 @@ -ChaconTarazona -************** - -.. automodule:: pytim.chacon_tarazona - :members: - :inherited-members: - :member-order: bysource - :exclude-members: LayerAtomGroup - -* :ref:`genindex` -* :ref:`modindex` -* :ref:`search` - -.. raw:: html - :file: analytics.html - diff --git a/docs/source/SASA.rst b/docs/source/SASA.rst new file mode 100644 index 00000000..b987bd57 --- /dev/null +++ b/docs/source/SASA.rst @@ -0,0 +1,18 @@ +SASA +**** + +.. include:: refs.rst + +.. automodule:: pytim.sasa + :members: SASA + :inherited-members: + :member-order: bysource + :exclude-members: LayerAtomGroup, alpha,cluster_cut, cluster_threshold_density,extra_cluster_groups,info,itim_group,max_layers,molecular,multiproc,radii_dict,save_pdb, surfaces + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` + +.. raw:: html + :file: analytics.html + diff --git a/docs/source/Tutorials.rst b/docs/source/Tutorials.rst deleted file mode 100644 index 836e63a9..00000000 --- a/docs/source/Tutorials.rst +++ /dev/null @@ -1,28 +0,0 @@ -:orphan: - -********* -Tutorials -********* - -.. toctree:: - :caption: Tutorials - :name: Tutorials - - quick - micelle - guessing_radii - - - - -.... - - - - - -.. toctree:: - -.. raw:: html - :file: analytics.html - diff --git a/docs/source/conf.py b/docs/source/conf.py index b409021a..e0ccbd7d 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -77,16 +77,16 @@ # built documents. # # The short X.Y version. -version = u'0.7' +version = u'1.0' # The full version, including alpha/beta/rc tags. -release = u'0.7.0' +release = u'1.0.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. # # This is also used if you do content translation via gettext catalogs. # Usually you set "language" from the command line for these cases. -language = None +language = 'en' # There are two options for replacing |today|: either, you set today to some # non-false value, then it is used: diff --git a/docs/source/index.rst b/docs/source/index.rst index 89dacddd..52ccb491 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -23,8 +23,8 @@ ITIM GITIM + SASA WillardChandler - ChaconTarazona SimpleInterface observables datafiles diff --git a/docs/source/install.rst b/docs/source/install.rst index 9a7f1f08..46dbeb8c 100644 --- a/docs/source/install.rst +++ b/docs/source/install.rst @@ -17,10 +17,6 @@ The package will download all the dependencies which are needed. Prerequisites f .. code-block:: bash - pip install setuptools --user --upgrade - pip install numpy --user --upgrade - pip install cython --user --upgrade - git clone https://github.com/Marcello-Sega/pytim.git cd pytim python setup.py install --user @@ -35,12 +31,22 @@ The package can also be installed directly from the Python Package Index, along .. code-block:: bash - pip install setuptools --user --upgrade - pip install cython --user --upgrade - pip install numpy --user --upgrade - + git clone https://github.com/Marcello-Sega/pytim.git + cd pytim pip install pytim --user --upgrade + +Alternatively, you can use pip (directly or as a module) to install pytim, possibly in developer/editable mode: + +.. code-block:: bash + + git clone https://github.com/Marcello-Sega/pytim.git + cd pytim + python -m pip install -e . + +This way the changes in the source will be reflected immediately after restarting a python shell, avoiding the need to reinstall every time. + + Using Anaconda -------------- diff --git a/docs/source/observables.rst b/docs/source/observables.rst index 7f2d5612..f75a846c 100644 --- a/docs/source/observables.rst +++ b/docs/source/observables.rst @@ -6,6 +6,8 @@ Radial Distribution Functions ############################# .. automodule:: pytim.observables.rdf :members: +.. automodule:: pytim.observables.rdf2d + :members: Profiles ######## @@ -17,11 +19,32 @@ Time Correlation Functions .. autoclass:: pytim.observables.Correlator :members: +Contact Angles +############## +.. autoclass:: pytim.observables.ContactAngle + :members: + +3D Distributions +################ +.. autoclass:: pytim.observables.DistributionFunction + :members: + +Free Volume +########### +.. autoclass:: pytim.observables.FreeVolume + :members: + +Local Reference Frame +##################### + + Misc #### +.. autoclass:: pytim.observables.IntrinsicDistance + :members: .. autoclass:: pytim.observables.LayerTriangulation :members: -.. autoclass:: pytim.observables.IntrinsicDistance +.. autoclass:: pytim.observables.LocalReferenceFrame :members: diff --git a/docs/source/quick.rst b/docs/source/quick.rst index 099d6d68..ae94af6f 100644 --- a/docs/source/quick.rst +++ b/docs/source/quick.rst @@ -1,15 +1,12 @@ Pytim: Quick Tour ***************** +.. include:: refs.rst + :doc:`Pytim ` is a package based on MDAnalysis_ for the identification and analysis of surface molecules in configuration files or in trajectories from molecular dynamics simulations. -Although MDAnalysis_ is needed to perform the interfacial analysis, you can also use :doc:`Pytim ` on top of MDTraj_ or directly online from a simulation performed using OpenMM_. For more information see the :doc:`Tutorials`. +Although MDAnalysis_ is needed to perform the interfacial analysis, you can also use :doc:`Pytim ` on top of MDTraj_ or directly online from a simulation performed using OpenMM_. -.. _MDAnalysis: http://www.mdanalysis.org/ -.. _MDTraj: http://www.mdtraj.org/ -.. _OpenMM: http://www.openmm.org/ -.. _Paraview: https://www.paraview.org/ -.. _Supported_Formats: https://pythonhosted.org/MDAnalysis/documentation_pages/coordinates/init.html#id1 Basic Example ============= @@ -63,7 +60,7 @@ where surface oxygen atoms are highlighted in blue. :width: 70% :align: center -This is a very basic example, and more are given below, in the :doc:`Tutorials`, and in the documentation of the modules. +This is a very basic example, and more are given below and in the documentation of the modules. Non-planar interfaces ===================== @@ -100,6 +97,37 @@ We make here the example of multiple solvation layers around glucose: :width: 40% :align: center +SASA +---- + +The Solvent-Accessible-Surface-Area method of Lee and Richards has been extensively used +to compute the exposed surface of proteins to calculate solvation and electrostatic effects. +However, the method can be also used to produce the list of surface-exposed atoms. In this sense, +it is a very similar to GITIM. In the present implementation, the identification of surface atoms +scales quasi-linear with the number of atoms. + +The above example for GITIM can be rewritten using :class:`~pytim.sasa.SASA` with very similar results: + +.. code-block:: python + + import MDAnalysis as mda + import pytim + from pytim.datafiles import GLUCOSE_PDB + + u = mda.Universe(GLUCOSE_PDB) + solvent = u.select_atoms('name OW') + glc = u.atoms - solvent.residues.atoms + + # alpha is the probe-sphere radius + inter = pytim.SASA(u, group=solvent, max_layers=3, alpha=2) + + for i in [0,1,2]: + print ("Layer "+str(i),repr(inter.layers[i])) + + Layer 0 + Layer 1 + Layer 2 + Willard-Chandler ---------------- diff --git a/docs/source/refs.rst b/docs/source/refs.rst new file mode 100644 index 00000000..1a342124 --- /dev/null +++ b/docs/source/refs.rst @@ -0,0 +1,8 @@ + +.. _MDAnalysis: http://www.mdanalysis.org/ +.. _MDTraj: http://www.mdtraj.org/ +.. _OpenMM: http://www.openmm.org/ +.. _Paraview: https://www.paraview.org/ +.. _Supported_Formats: https://pythonhosted.org/MDAnalysis/documentation_pages/coordinates/init.html#id1 + + diff --git a/pytim/__init__.py b/pytim/__init__.py index b9ad9216..fb5701b5 100644 --- a/pytim/__init__.py +++ b/pytim/__init__.py @@ -7,7 +7,6 @@ from .gitim import GITIM from .sasa import SASA from .willard_chandler import WillardChandler -from .chacon_tarazona import ChaconTarazona from . import observables, utilities, datafiles from .version import __version__ import warnings diff --git a/pytim/chacon_tarazona.py b/pytim/chacon_tarazona.py deleted file mode 100644 index e73f05ef..00000000 --- a/pytim/chacon_tarazona.py +++ /dev/null @@ -1,221 +0,0 @@ -#!/usr/bin/python -# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 -""" Module: chacon_tarazona - ======================= -""" -from __future__ import print_function -import numpy as np -from . import utilities -from .surface import SurfaceFlatInterface as Surface -from .sanity_check import SanityCheck - -from .interface import Interface -from .patches import patchTrajectory, patchOpenMM, patchMDTRAJ - - -class ChaconTarazona(Interface): - """ Identifies the dividing surface using the Chacon-Tarazona method - - *(Chacón, E.; Tarazona, P. Phys. Rev. Lett. 91, 166103, 2003)* - *(Tarazona, P.; Chacón, E. Phys. Rev. B 70, 235407, 2004)* - - :param Universe universe: The MDAnalysis_ universe - :param float alpha: Molecular scale cutoff - :param float tau: Particles within this distance form the - surface will be added during the - self-consistent procedure. - :param bool molecular: Switches between search of interfacial - molecules / atoms (default: True) - :param AtomGroup group: Compute the density using this group - :param dict radii_dict: Dictionary with the atomic radii of the - elements in the group. If None is supplied, - the default one (from GROMOS43a1) will be - used. - - Example: - - >> import MDAnalysis as mda - >> import numpy as np - >> import pytim - >> from pytim.datafiles import WATER_GRO - >> u = mda.Universe(WATER_GRO) - >> g = u.select_atoms('name OW') - >> interface = pytim.ChaconTarazona(u,alpha=2.,tau=1.5,group=g,info=False,molecular=False) - >> interface.writepdb('CT.pdb',centered=True) # doctest: +SKIP - >> print (repr(interface.layers)) - array([[], - []], dtype=object) - - .. _MDAnalysis: http://www.mdanalysis.org/ - """ - _surface = None - - @property - def layers(self): - return self._layers - - def __init__(self, - universe, - alpha=2.0, - tau=1.5, - group=None, - radii_dict=None, - max_layers=1, - normal='guess', - molecular=True, - info=True, - mesh=None, - centered=False, - warnings=False, - autoassign=True, - **kargs): - - self.autoassign = autoassign - self.symmetry = 'planar' - self.do_center = centered - - sanity = SanityCheck(self, warnings=warnings) - sanity.assign_universe(universe, group) - - self.target_mesh = mesh - if mesh is not None: - raise Exception("FFT-based version not implemented") - self.info = info - self.alpha = alpha - self.tau = tau - if (max_layers != 1): - raise Exception("max_layers !=1 not implemented yet!") - - self.max_layers = max_layers - self._layers = np.empty([2, max_layers], dtype=type(universe.atoms)) - self._surfaces = np.empty(max_layers, dtype=type(Surface)) - self.normal = None - self.PDB = {} - - self.molecular = molecular - - # TODO implement cluster group - sanity.assign_cluster_params(None, None, None) - sanity.assign_normal(normal) - sanity.assign_radii(radii_dict=radii_dict) - - self.sorted_indices = None - self.surf = None - self.modes = [None, None] - - patchTrajectory(self.universe.trajectory, self) - self._assign_layers() - - def _points_next_to_surface(self, surf, modes, pivot): - """ Searches for points within a distance self.tau from the - interface. - """ - pivot_pos = self.cluster_group[pivot].positions - z_max = np.max(pivot_pos[:, 2]) - z_min = np.min(pivot_pos[:, 2]) - z_max += self.alpha * 2 - z_min -= self.alpha * 2 - positions = self.cluster_group.positions[:] - # TODO other directions - z = positions[:, 2] - condition = np.logical_and(z > z_min, z < z_max) - candidates = np.argwhere(condition)[:, 0] - dists = surf.surface_from_modes(positions[candidates], modes) - dists = dists - z[candidates] - return candidates[dists * dists < self.tau**2] - - def _initial_pivots(self, sorted_ind): - """ Defines the initial pivots as a set of 9 particles, where - each particle is in a distinct sector formed by dividing - the macroscopic plane into 3x3 regions. - """ - sectors = np.zeros((3, 3), dtype=int) - pivot = [] - box = utilities.get_box(self.universe, normal=self.normal) - pos = utilities.get_pos(self.cluster_group, normal=self.normal) - for ind in sorted_ind: - part = pos[ind] - nx, ny = list(map(int, 2.999 * part[0:2] / box[0:2])) - if sectors[nx, ny] == 0: - pivot.append(ind) - sectors[nx, ny] = 1 - if np.sum(sectors) >= 9: - break - return pivot - - def _assign_one_side(self, side): - """ Calculate the interfacial molecules on one side of the box. - """ - # TODO add successive layers - box = self.universe.dimensions[:3] - - if side == 0: - sorted_ind = self.sorted_indices[::-1] - else: - sorted_ind = self.sorted_indices - - pivot = np.sort(self._initial_pivots(sorted_ind)) - modes = None - surf = Surface(self, options={'layer': 0, 'from_modes': True}) - surf._compute_q_vectors(box) - while True: - p = self.cluster_group[pivot].positions - modes = surf.surface_modes(p) - s = surf.surface_from_modes(p, modes) - d = p[:, 2] - s - if self.info is True: - # d is too large means the decomposition failed - print("side", side, "->", len(pivot), "pivots, msd=", - np.sqrt(np.sum(d * d) / len(d))) - # TODO handle failure - new_pivot = np.sort( - self._points_next_to_surface(surf, modes, pivot)) - # If convergence reached... - if np.all(new_pivot == pivot): - self.surf = surf - self.modes[side] = modes - _inlayer_group = self.cluster_group[pivot] - if self.molecular is True: - _inlayer_group = _inlayer_group.residues.atoms - return _inlayer_group - - else: - pivot = new_pivot - - def _assign_layers(self): - """ This function identifies the dividing surface and the atoms in the - layers. - - """ - self.reset_labels() - - # TODO parallelize - self.prepare_box() - - # groups have been checked already in _sanity_checks() - - self._define_cluster_group() - - self.centered_positions = None - # we always (internally) center in Chacon-Tarazona - self.center(planar_to_origin=True) - # first we label all atoms in group to be in the gas phase - self.label_group(self.analysis_group.atoms, beta=0.5) - # then all atoms in the largest group are labelled as liquid-like - self.label_group(self.cluster_group.atoms, beta=0.0) - - pos = self.cluster_group.positions - self.sorted_indices = np.argsort(pos[:, self.normal]) - for side in [0, 1]: - self._layers[side][0] = self._assign_one_side(side) - - self.label_planar_sides() - - if self.do_center is False: - self.universe.atoms.positions = self.original_positions - else: - self._shift_positions_to_middle() - - -# diff --git a/pytim/gitim.py b/pytim/gitim.py index 4cf32013..0756fcc6 100644 --- a/pytim/gitim.py +++ b/pytim/gitim.py @@ -26,8 +26,7 @@ class GITIM(Interface): """ Identifies interfacial molecules at curved interfaces. - *(Sega, M.; Kantorovich, S.; Jedlovszky, P.; Jorge, M., \ -J. Chem. Phys. 138, 044110, 2013)* + *(Sega, M.; Kantorovich, S.; Jedlovszky, P.; Jorge, M., J. Chem. Phys. 138, 044110, 2013)* :param Object universe: The MDAnalysis_ Universe, MDTraj_ trajectory or OpenMM_ Simulation objects. @@ -68,6 +67,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 +127,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 +159,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 +259,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 diff --git a/pytim/interface.py b/pytim/interface.py index e3a9bd0f..fa169c41 100644 --- a/pytim/interface.py +++ b/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") diff --git a/pytim/itim.py b/pytim/itim.py index 1e154d9a..d2535f17 100644 --- a/pytim/itim.py +++ b/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 diff --git a/pytim/observables/contactangle.py b/pytim/observables/contactangle.py index 15c2560e..f7fa0335 100644 --- a/pytim/observables/contactangle.py +++ b/pytim/observables/contactangle.py @@ -15,22 +15,22 @@ class ContactAngle(object): - """ ContactAngle class implementation that uses interfacial atoms to compute + """ ContactAngle class implementation that uses interfacial atoms to compute\ droplet profiles and contact angles using different approaches. - :param PYTIM interface : Compute the contact angle for this interface - :param AtomGroup droplet : Atom group representative of the droplet - :param AtomGroup substrate : Atom group representative of the substrate - :param int periodic : direction along which the system is periodic. Default: None, not periodic - If None, the code performs best fit to ellipsoids, otherwise, to ellipses - If not None, selects the direction of the axis of macroscopic translational + :param PYTIM interface: Compute the contact angle for this interface + :param AtomGroup droplet: Atom group representative of the droplet + :param AtomGroup substrate: Atom group representative of the substrate + :param int periodic: direction along which the system is periodic. Default: None, not periodic\ + If None, the code performs best fit to ellipsoids, otherwise, to ellipses\ + If not None, selects the direction of the axis of macroscopic translational\ invariance: 0:x, 1:y, 2:z - :param float hcut : don't consider contributions from atoms closer than this to the substrate + :param float hcut: don't consider contributions from atoms closer than this to the substrate\ (used to disregard the microscopic contact angle) - :param float hcut_upper : don't consider contributions from atoms above this distance from the substrate + :param float hcut_upper: don't consider contributions from atoms above this distance from the substrate\ default: None - :param int bins : bins used for sampling the profile - :param int or array_like removeCOM : remove the COM motion along this direction(s). Default: None, does not remove COM motion + :param int bins: bins used for sampling the profile + :param int or array_like removeCOM: remove the COM motion along this direction(s). Default: None, does not remove COM motion Example: @@ -289,14 +289,15 @@ def rmsd_ellipse_penalty(p, x, N, check_coeffs=True): @staticmethod def rmsd_ellipsoid(p, x, N, check_coeffs=True): - """ RMSD between the points x and the ellipsoid defined by the general + """ RMSD between the points x and the ellipsoid defined by the general\ parameters p of the associated polynomial. - :param p list : general coefficients [a,b,c,f,g,h,p,q,r,d] - :param x ndarray : points coordinates as a (*,3)-ndarray - :param N int : number of points on the ellipsoid that are - generated and used to compute the rmsd - :param check_coeffs bool : if true, perform additional checks + :param list p: general coefficients [a,b,c,f,g,h,p,q,r,d] + :param ndarray x: points coordinates as a (x,3)-ndarray + :param int N: number of points on the ellipsoid that are\ + generated and used to compute the rmsd + :param bool check_coeffs: if true, perform additional checks + """ cp = ContactAngle._ellipsoid_general_to_affine(p, check_coeffs) xx, yy, zz = ContactAngle.ellipsoid(cp, N) @@ -320,16 +321,18 @@ def rmsd_circle(p, x): @staticmethod def ellipse(parmsc, npts=100, tmin=0., tmax=2.*np.pi): - """ Return npts points on the ellipse described by the canonical parameters + """ Return npts points on the ellipse described by the canonical parameters\ x0, y0, ap, bp, e, phi for values of the paramter between tmin and tmax. - :param dict parmsc : dictionary with keys: x0,y0,a,b,phi - :param float tmin : minimum value of the parameter - :param float tmax : maximum value of the parameter - :param int npts : number of points to use + :param dict parmsc: dictionary with keys: x0,y0,a,b,phi + :param float tmin: minimum value of the parameter + :param float tmax: maximum value of the parameter + :param int npts: number of points to use + + :return: + + :tuple: (x,y): coordinates as numpy arrays - return: - :tuple : (x,y) of coordinates as numpy arrays """ t = np.linspace(tmin, tmax, npts) x = parmsc['x0'] + parmsc['a'] * np.cos(t) * np.cos( @@ -343,13 +346,15 @@ def circle(parmsc, npts=100, tmin=0., tmax=2.*np.pi): """ Return npts points on the circle described by the canonical parameters R, x0, y0 for values of the paramter between tmin and tmax. - :param dict parmsc : dictionary with keys: R, x0, y0 - :param float tmin : minimum value of the parameter - :param float tmax : maximum value of the parameter - :param int npts : number of points to use + :param dict parmsc: dictionary with keys: R, x0, y0 + :param float tmin: minimum value of the parameter + :param float tmax: maximum value of the parameter + :param int npts: number of points to use + + :return: + + :tuple: (x,y): coordinates as numpy arrays - return: - :tuple : (x,y) of coordinates as numpy arrays """ t = np.linspace(tmin, tmax, npts) x = parmsc['x0'] + parmsc['R'] * np.cos(t) @@ -358,14 +363,14 @@ def circle(parmsc, npts=100, tmin=0., tmax=2.*np.pi): @staticmethod def ellipsoid(parmsc, npts=1000): - """ Return npts points on the ellipsoid described by the affine parameters - T, v + """ Compute npts points on the ellipsoid described by the affine parameters T, v - :param dict parmsc : dictionary with keys: T (3x3 matrix), v (3x1 vector) - :param int npts : number of points to use + :param dict parmsc: dictionary with keys: T (3x3 matrix), v (3x1 vector) + :param int npts: number of points to use + + :return: + :tuple: (x,y,z) : coordinates of points on the ellipsoid as ndarrays - return: - :tuple : (x,y,z) coordinates of points on the ellipsoid as ndarrays """ phi = np.arccos(2 * np.random.rand(npts) - 1) theta = 2 * np.pi * np.random.rand(npts) @@ -377,12 +382,13 @@ def ellipsoid(parmsc, npts=1000): @staticmethod def _fit_circle(hr, hh, nonlinear=True): """ fit an arc through the profile h(r) sampled by the class - :param list hr : list of arrays with the radial coordinates - :param list hh : list of arrays with the elevations - :param bool nonlinear : use the more accurate minimization of the rmsd instead of the algebraic distance - return: - :list : a list with the tuple (radius, base radius, cos(theta), center, rmsd) + :param list hr: list of arrays with the radial coordinates + :param list hh: list of arrays with the elevations + :param bool nonlinear: use the more accurate minimization of the rmsd instead of the algebraic distance + + :return: + :list: : a list with the tuple (radius, base radius, cos(theta), center, rmsd) for each bin. If only one bin is present, return just the tuple. """ parms = [] @@ -423,15 +429,15 @@ def _fit_circle(hr, hh, nonlinear=True): def _contact_angles_from_ellipse(p, off=0.0): """ compute the contact angle from the parameters of the polynomial representation - :parms array_like p : a sequence (list, tuple, ndarray) of parameters of the ellipse equation - in general form: + :parms array_like p: a sequence (list, tuple, ndarray) of parameters of the ellipse equation\ + in general form:\ p[0] x^2 + p[1] x y + p[2] y^2 + p[3] x + p[4] y + p[5] = 0 - :param float off : elevation from the substrate surface, where the - contact angle should be evaluated. Default; 0.0, corresponding + :param float off: elevation from the substrate surface, where the\ + contact angle should be evaluated. Default; 0.0, corresponding\ to the atomic center of the highest atom of the substrate - return: - :tuple : a tuple (theta1,theta2) with the left and right (internal) - contact angles + + :return: + :tuple : a tuple (theta1,theta2) with the left and right (internal) contact angles We require the solution of the line passing through the point to have two coincident solutions for the system ax2 + bxy + cy2 + dx + ey +f = 0 and y-y0 = m (x-x0) @@ -518,18 +524,19 @@ def _fit_ellipse_flusser(x, y): @staticmethod def _fit_ellipsoid_lstsq_cox(x, y, z): """ - Return the polynomial coefficients that perform a least square - fit to a set of points, using a modified version of: - Turner, D. A., I. J. Anderson, J. C. Mason, and M. G. Cox. - "An algorithm for fitting an ellipsoid to data." + Return the polynomial coefficients that perform a least square\ + fit to a set of points, using a modified version of:\ + Turner, D. A., I. J. Anderson, J. C. Mason, and M. G. Cox.\ + "An algorithm for fitting an ellipsoid to data."\ National Physical Laboratory, UK (1999). - :param array_like x : list or ndarray array with coordinate x - :param array_like y : list or ndarray array with coordinate y - :param array_like z : list or ndarray array with coordinate z + :param array_like x: list or ndarray array with coordinate x + :param array_like y: list or ndarray array with coordinate y + :param array_like z: list or ndarray array with coordinate z - return: - :ndarray : a (10,)-ndarray with the polynomial coefficients + :return: + + :ndarray: a (10,)-ndarray with the polynomial coefficients """ D = np.vstack([x**2+y**2-2*z**2, x**2 - 2*y**2 + z**2, 4*x*y, 2*x*z, 2*y*z, x, y, z, np.ones(len(x))]).T @@ -561,7 +568,7 @@ def _fit_ellipse(x, y, nonlinear=True, off=0.0, points_density=25, verbose=False to the atomic center of the highest atom of the substrate :param int points_density: number of points per Angstrom on the ellipse that are used to compute the rmsd - return: + :return: :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins If only one bin is present, return just the tuple @@ -634,7 +641,8 @@ def _fit_ellipsoid(x, y, z, nonlinear=True, off=0.0, points_density=25): to the atomic center of the highest atom of the substrate :param int points_density: number of points per Angstrom on the ellipse that are used to compute the rmsd - return: + + :return: :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins If only one bin is present, return just the tuple @@ -710,7 +718,7 @@ def _ellipse_general_to_canonical(coeffs, check_coeffs=True): during the minimization the constraints might be violated. - return: + :return: : dict : a dictionary with the canonical coefficients x0,y0,a,b,phi,e from https://scipython.com/blog/direct-linear-least-squares-fitting-of-an-ellipse/ @@ -826,7 +834,7 @@ def _ellipsoid_contact_line_and_angles(coeffs, off, npts=1000): the contact line :param int npts : number of points to sample - return: + :return: :tuple : a tuple containing the contact line on the (x,y) plane as a (npts,2)-ndarray and the contact angles as a (npts,)-ndarray. @@ -851,7 +859,7 @@ def _ellipsoid_general_to_affine(coeffs, check_coeffs=True): during the minimization the constraints might be violated. - return: + :return: : dict : a dictionary with the affine matrix and vector, T and v If a point on the unit sphere centered in the origin is s, then the correspoinding point r=(x,y,z) on the ellipsoid is @@ -1044,36 +1052,38 @@ def _select_coords(self, use, bins): def fit_circle(self, use='frame', nonlinear=True, bins=1): """ fit an arc through the profile h(r) sampled by the class - :param str use : 'frame' : use the positions of the current frame only (default) - 'histogram': use the binned values sampled so far - 'stored' : use the stored surface atoms positions, if the option store=True was passed at initialization - :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) + :param str use: 'frame' : use the positions of the current frame only (default)\ + 'histogram': use the binned values sampled so far\ + 'stored' : use the stored surface atoms positions, if the option store=True was passed at initialization + :param int bins: the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) - return: + :return: a list including, for each of the bins: - : tuple : radius, base radius, cos(theta), center + + :tuple: radius, base radius, cos(theta), center """ r, h = self._select_coords(use, bins=bins) return self._fit_circle(r, h, nonlinear=nonlinear) def fit_ellipse(self, use='frame', nonlinear=True, bins=1): - """ fit an ellipse through the points sampled by the class. See implementation details in _fit_ellipsoid() + """ fit an ellipse through the points sampled by the class. See implementation details in _fit_ellipsoid() - :param str use : 'frame' : use the positions of the current frame only (default) - 'histogram': use the binned values sampled so far - 'stored' : use the stored surface atoms positions, if the option store=True - was passed at initialization - :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) + :param str use: 'frame' : use the positions of the current frame only (default)\ + 'histogram': use the binned values sampled so far\ + 'stored' : use the stored surface atoms positions, if the option store=True\ + was passed at initialization + :param int bins: the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) :return: a list including, for each of the bins, a tuple with elements: - parms : parameters of the ellipse polynomial in general form: - a[0] x^2 + a[1] x y + a[2] y^2 + a[3] x + a[4] y = 0 - parmsc : dictionary of parameters in canoncial form: (a,b, x0,y0,phi, e) - with a,b the major and minor semiaxes, x0,y0 the center, phi the angle - (in rad) between x axis and major axis, and e the eccentricity. - theta : [left contact angle, right contact angle] + + :list: parms: parameters of the ellipse polynomial in general form:\ + a[0] x^2 + a[1] x y + a[2] y^2 + a[3] x + a[4] y = 0 + :dict: parmsc: dictionary of parameters in canoncial form: (a,b, x0,y0,phi, e)\ + with a,b the major and minor semiaxes, x0,y0 the center, phi the angle\ + (in rad) between x axis and major axis, and e the eccentricity. + :list: theta: [left contact angle, right contact angle] """ r, h = self._select_coords(use, bins=bins) @@ -1085,25 +1095,27 @@ def fit_ellipse(self, use='frame', nonlinear=True, bins=1): def fit_ellipsoid(self, use='frame', nonlinear=True, bins=1): """ fit an ellipsoid through the points sampled by the class. See implementation details in _fit_ellipsoid() - :param str use : 'frame' : use the positions of the current frame only (default) - 'histogram': use the binned values sampled so far - 'stored' : use the stored surface atoms positions, if the option store=True + :param str use : 'frame' : use the positions of the current frame only (default)\ + 'histogram': use the binned values sampled so far\ + 'stored' : use the stored surface atoms positions, if the option store=True\ was passed at initialization :param int bins : the number of bins to use along the symmetry direction (cylinder axis, azimuthal angle) :return: - :list : a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins - If only one bin is present, return just the tuple - parms : array of parameters of the ellipsoid equation in general form: - a[0] x^2 + a[1] y^2 + a[2] z^2 + a[3] yz + a[4] xz + + + a list with a tuple (parms,parmsc,theta, rmsd) for each of the bins;\ + If only one bin is present, return just the tuple. + + :array: parms : array of parameters of the ellipsoid equation in general form:\ + a[0] x^2 + a[1] y^2 + a[2] z^2 + a[3] yz + a[4] xz + \ a[5] xy + a[6] x + a[7] y + a[8] z + a[9] = 0, otherwise. - parmsc : dictionary with parameters of the ellipsoid affine form (T,v), - such that the ellipsoid points are - r = T s + v, + :dict: parmsc : dictionary with parameters of the ellipsoid affine form (T,v),\ + such that the ellipsoid points are\ + r = T s + v,\ if s are points from the unit sphere centered in the origin. - theta : the contact angle as a function of the azimuthal angle phi from phi=0, aligned + :float: theta : the contact angle as a function of the azimuthal angle phi from phi=0, aligned\ with (-1,0,0) to 2 pi. - rmsd : the rmsd to the best fit (linear or nonlinear) ellipsoid + :float: rmsd : the rmsd to the best fit (linear or nonlinear) ellipsoid """ x, y, z = self._select_coords(use, bins=bins) # TODO FIXME diff --git a/pytim/observables/free_volume.py b/pytim/observables/free_volume.py index 5f1059d4..484b0930 100644 --- a/pytim/observables/free_volume.py +++ b/pytim/observables/free_volume.py @@ -18,21 +18,23 @@ class FreeVolume(object): :param Universe universe: the universe :param int npoints: number of Monte Carlo sampling points (default: 10x the number of atoms in the universe) - Returns: - :(free volume, error) tuple : A tuple with the free volume and the estimated error + :return: + + :tuple: free volume, error : A tuple with the free volume and the estimated error Examples: + >>> import MDAnalysis as mda >>> import numpy as np >>> import pytim >>> from pytim.datafiles import CCL4_WATER_GRO, _TEST_BCC_GRO - + >>> >>> u = mda.Universe(CCL4_WATER_GRO) >>> inter = pytim.ITIM(u) # just to make sure all radii are set >>> np.random.seed(1) # ensure reproducibility of test >>> FV = pytim.observables.FreeVolume(u) >>> bins, prof ,err = FV.compute_profile() - + >>> >>> free, err = FV.compute() >>> print ('{:0.3f} +/- {:0.3f}'.format(free,err)) 0.431 +/- 0.001 @@ -95,7 +97,9 @@ def compute_profile(self, inp=None, nbins=30, direction=2): :param int nbins: number of bins, by default 30 :param int direction: direction along wich to compute the the profile, z (2) by default - :returns bins,fraction,error: the left limit of the bins, the free volume fraction in each bin, the associated std deviation + :return: + + :tuple: bins,fraction,error: the left limit of the bins, the free volume fraction in each bin, the associated std deviation """ box = self.u.dimensions[:3].copy() @@ -127,7 +131,9 @@ def compute(self, inp=None): :param AtomGroup inp: compute the volume fraction of this group, None selects the complete universe :param int nbins: number of bins, by default 30 - :returns fraction, error: the free volume fraction and associated error + :return: + + :tuple: fraction, error: the free volume fraction and associated error """ _, free, err = self.compute_profile(inp, nbins=1) diff --git a/pytim/observables/rdf.py b/pytim/observables/rdf.py index 565b3af6..92c4c270 100644 --- a/pytim/observables/rdf.py +++ b/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) diff --git a/pytim/sasa.py b/pytim/sasa.py index 30c5feac..1491c547 100644 --- a/pytim/sasa.py +++ b/pytim/sasa.py @@ -22,9 +22,9 @@ class SASA(GITIM): - """ Identifies interfacial molecules at curved interfaces using the - Lee-Richards SASA algorithm (Lee, B; Richards, FM. J Mol Biol. - 55, 379–400, 1971) + """ Identifies interfacial molecules at curved interfaces using the Lee-Richards SASA algorithm. + + *(Lee, B; Richards, FM. J Mol Biol. 55, 379–400, 1971)* :param Object universe: The MDAnalysis_ Universe, MDTraj_ trajectory or OpenMM_ Simulation objects. @@ -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 diff --git a/pytim/simple_interface.py b/pytim/simple_interface.py index cedb7e1a..d47a2a45 100644 --- a/pytim/simple_interface.py +++ b/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,6 +85,7 @@ def __init__(self, alpha=1.5, symmetry='generic', normal='z', + include_zero_radius=False, upper=None, lower=None): @@ -89,6 +93,7 @@ def __init__(self, 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: diff --git a/pytim/version.py b/pytim/version.py index ecc0ae4b..1f356cc5 100644 --- a/pytim/version.py +++ b/pytim/version.py @@ -1 +1 @@ -__version__ = '0.9.3' +__version__ = '1.0.0' diff --git a/pytim/willard_chandler.py b/pytim/willard_chandler.py index 71bff17b..9f433db1 100644 --- a/pytim/willard_chandler.py +++ b/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( diff --git a/setup.py b/setup.py index 649cdd6f..ba605a2c 100644 --- a/setup.py +++ b/setup.py @@ -91,7 +91,7 @@ def run_tests(self): # 3 - Alpha # 4 - Beta # 5 - Production/Stable - 'Development Status :: 4 - Beta', + 'Development Status :: 5 - Production/Stable', # Indicate who your project is intended for 'Intended Audience :: Science/Research',