merge update

prody/atomic/functions.py

@@ -424,8 +424,14 @@ def getClosest(a, B):
         ag = atoms.getAtomGroup()
     except AttributeError:
         ag = atoms
+
+    if hasattr(atoms, "getTitle"):
+        title = atoms.getTitle()
+    else:
+        title = str(atoms)
+
     atommap = AtomMap(ag, atom_indices, atoms.getACSIndex(),
-                      title=str(atoms), intarrays=True)
+                      title=title, intarrays=True)
 
     #LOGGER.report('Full atoms was extended in %2.fs.', label='_prody_extendAtoms')
 

prody/dynamics/__init__.py

@@ -152,6 +152,7 @@
   * :func:`.writeModes` - normal modes
   * :func:`.writeNMD` - normal mode, coordinate, and atomic data
   * :func:`.writeOverlapTable` - overlap between modes in a formatted table
+  * :func:`.writeBILD` - normal mode and coordinate data for ChimeraX
 
 Save/load models
 ================
@@ -315,6 +316,10 @@
 from .vmdfile import *
 __all__.extend(vmdfile.__all__)
 
+from . import bildfile
+from .bildfile import *
+__all__.extend(bildfile.__all__)
+
 #from . import bbenm
 #from .bbenm import *
 #__all__.extend(bbenm.__all__)

prody/dynamics/bildfile.py

@@ -0,0 +1,148 @@
+"""This module defines input and output functions for BILD format for ChimeraXs.
+
+.. _bild-format: https://rbvi.ucsf.edu/chimerax/docs/user/formats/bild.html
+
+"""
+
+
+__all__ = ['writeBILD']
+
+import os
+from os.path import join, split, splitext
+
+import numpy as np
+
+from prody import LOGGER
+from prody.utilities import openFile, getCoords
+
+from .nma import NMA
+from .anm import ANM
+from .gnm import GNM, ZERO
+from .pca import PCA
+from .mode import Vector, Mode
+from .modeset import ModeSet
+
+def writeBILD(filename, modes, atoms, scale=25.):
+    """Returns *filename* that contains *modes* and *atoms* data in BILD format
+    described at :ref:`bild-format`.  :file:`.bild` extension is appended to
+    filename, if it does not have an extension. 
+    
+    Arrows will be scaled by the square root of the variance times *scale*.
+
+    .. note::
+       #. This function skips modes with zero eigenvalues.
+       #. If a :class:`.Vector` instance is given, it will be normalized
+          before it is written. It's length before normalization will be
+          written as the scaling factor of the vector.
+       #. It is recommended to only use one mode to avoid having lots of arrows.
+    """
+
+    if not filename.lower().endswith(".bild"):
+        filename += '.bild'
+
+    if not isinstance(modes, (NMA, ModeSet, Mode, Vector)):
+        raise TypeError('modes must be NMA, ModeSet, Mode, or Vector, '
+                        'not {0}'.format(type(modes)))
+    if modes.numAtoms() != len(atoms):
+        raise Exception('number of atoms do not match')
+    out = openFile(filename, 'w')
+
+    name = modes.getTitle()
+    name = name.replace(' ', '_').replace('.', '_')
+    if not name.replace('_', '').isalnum() or len(name) > 30:
+        name = str(atoms)
+        name = name.replace(' ', '_').replace('.', '_')
+    if not name.replace('_', '').isalnum() or len(name) > 30:
+        name = splitext(split(filename)[1])[0]
+    out.write('.comment name {0}\n'.format(name))
+    try:
+        coords = getCoords(atoms)
+    except:
+        raise ValueError('coordinates could not be retrieved from atoms')
+    if coords is None:
+        raise ValueError('atom coordinates are not set')
+
+    try:
+        data = atoms.getNames()
+        if data is not None:
+            out.write('.comment atomnames {0}\n'.format(' '.join(data)))
+    except:
+        pass
+    try:
+        data = atoms.getResnames()
+        if data is not None:
+            out.write('.comment resnames {0}\n'.format(' '.join(data)))
+    except:
+        pass
+    try:
+        data = atoms.getResnums()
+        if data is not None:
+            out.write('.comment resids ')
+            data.tofile(out, ' ')
+            out.write('\n')
+    except:
+        pass
+    try:
+        data = atoms.getChids()
+        if data is not None:
+            out.write('.comment chainids {0}\n'.format(' '.join(data)))
+    except:
+        pass
+    try:
+        data = atoms.getSegnames()
+        if data is not None:
+            out.write('.comment segnames {0}\n'.format(' '.join(data)))
+    except:
+        pass
+
+    try:
+        data = atoms.getBetas()
+        if data is not None:
+            out.write('.comment bfactors ')
+            data.tofile(out, ' ', '%.2f')
+            out.write('\n')
+    except:
+        pass
+
+    format = '{0:.3f}'.format
+    out.write('.comment coordinates ')
+    coords.tofile(out, ' ', '%.3f')
+    out.write('\n')
+    count = 0
+    if isinstance(modes, Vector):
+        out.write('.comment mode 1 {0:.2f} \n'.format(abs(modes)))
+        normed = modes.getNormed()._getArray() * scale
+        for i in range(modes.numAtoms()):
+            out.write('.arrow ')
+            coords[3*i:3*i+1].tofile(out, ' ', '%.3f')
+            normed[3*i:3*i+1].tofile(out, ' ', '%.3f')
+            out.write('\n')
+        out.write('\n')
+        count += 1
+    else:
+        if isinstance(modes, Mode):
+            modes = [modes]
+        else:
+            LOGGER.warning('It is recommended to only provide one mode '
+                'to writeBILD as otherwise there will be a mess of arrows.')
+        for mode in modes:
+            if mode.getEigval() < ZERO:
+                continue
+            out.write('.comment mode {0} {1:.2f} \n'.format(
+                       mode.getIndex()+1, mode.getVariance()**0.5))
+
+            arr = mode._getArray() * mode.getVariance()**0.5 * scale
+            for i in range(mode.numAtoms()):
+                out.write('.arrow ')
+                coords[i].tofile(out, ' ', '%.3f')
+                out.write(' ')
+
+                arr_out = arr[3*i:3*i+3] + coords[i]
+                arr_out.tofile(out, ' ', '%.3f')
+                out.write('\n')
+            count += 1
+    if count == 0:
+        LOGGER.warning('No normal mode data was written. '
+                       'Given modes might have 0 eigenvalues.')
+    out.close()
+    return filename

prody/dynamics/gnm.py

@@ -75,46 +75,6 @@ def _getKirchhoff(self):
 
         return self._kirchhoff
 
-
-def checkENMParameters(cutoff, gamma):
-    """Check type and values of *cutoff* and *gamma*."""
-
-    if not isinstance(cutoff, (float, int)):
-        raise TypeError('cutoff must be a float or an integer')
-    elif cutoff < 4:
-        raise ValueError('cutoff must be greater or equal to 4')
-    if isinstance(gamma, Gamma):
-        gamma_func = gamma.gamma
-    elif isinstance(gamma, FunctionType):
-        gamma_func = gamma
-    else:
-        if not isinstance(gamma, (float, int)):
-            raise TypeError('gamma must be a float, an integer, derived '
-                            'from Gamma, or a function')
-        elif gamma <= 0:
-            raise ValueError('gamma must be greater than 0')
-        gamma = float(gamma)
-        gamma_func = lambda dist2, i, j: gamma
-    return cutoff, gamma, gamma_func
-
-
-class GNM(GNMBase):
-
-    """A class for Gaussian Network Model (GNM) analysis of proteins
-    ([IB97]_, [TH97]_).
-
-    See example :ref:`gnm`.
-
-    .. [IB97] Bahar I, Atilgan AR, Erman B. Direct evaluation of thermal
-       fluctuations in protein using a single parameter harmonic potential.
-       *Folding & Design* **1997** 2:173-181.
-
-    .. [TH97] Haliloglu T, Bahar I, Erman B. Gaussian dynamics of folded
-       proteins. *Phys. Rev. Lett.* **1997** 79:3090-3093."""
-
-    def __init__(self, name='Unknown'):
-        super(GNM, self).__init__(name)
-
     def setKirchhoff(self, kirchhoff):
         """Set Kirchhoff matrix."""
 
@@ -324,7 +284,48 @@ def getNormDistFluct(self, coords):
 
     def setEigens(self, vectors, values=None):
         self._clear()
-        super(GNM, self).setEigens(vectors, values)
+        super(GNMBase, self).setEigens(vectors, values)
+
+
+def checkENMParameters(cutoff, gamma):
+    """Check type and values of *cutoff* and *gamma*."""
+
+    if not isinstance(cutoff, (float, int)):
+        raise TypeError('cutoff must be a float or an integer')
+    elif cutoff < 4:
+        raise ValueError('cutoff must be greater or equal to 4')
+    if isinstance(gamma, Gamma):
+        gamma_func = gamma.gamma
+    elif isinstance(gamma, FunctionType):
+        gamma_func = gamma
+    else:
+        if not isinstance(gamma, (float, int)):
+            raise TypeError('gamma must be a float, an integer, derived '
+                            'from Gamma, or a function')
+        elif gamma <= 0:
+            raise ValueError('gamma must be greater than 0')
+        gamma = float(gamma)
+        gamma_func = lambda dist2, i, j: gamma
+    return cutoff, gamma, gamma_func
+
+
+class GNM(GNMBase):
+
+    """A class for Gaussian Network Model (GNM) analysis of proteins
+    ([IB97]_, [TH97]_).
+
+    See example :ref:`gnm`.
+
+    .. [IB97] Bahar I, Atilgan AR, Erman B. Direct evaluation of thermal
+       fluctuations in protein using a single parameter harmonic potential.
+       *Folding & Design* **1997** 2:173-181.
+
+    .. [TH97] Haliloglu T, Bahar I, Erman B. Gaussian dynamics of folded
+       proteins. *Phys. Rev. Lett.* **1997** 79:3090-3093."""
+
+    def __init__(self, name='Unknown'):
+        super(GNM, self).__init__(name)
+
 
 class MaskedGNM(GNM, MaskedNMA):
     def __init__(self, name='Unknown', mask=False, masked=True):

prody/dynamics/plotting.py

@@ -200,10 +200,11 @@ def showProjection(ensemble, modes, *args, **kwargs):
     :arg modes: up to three normal modes
     :type modes: :class:`.Mode`, :class:`.ModeSet`, :class:`.NMA`
 
-    :keyword by_time: whether to show a 1D projection by time (number of steps) 
-        on the x-axis, rather than making a population histogram. 
-        Default is **False** to maintain old behaviour.
-    :type by_time: bool
+    :keyword show_density: whether to show a density histogram or kernel density estimate
+        rather than points or a 1D projection by time (number of steps) 
+        on the x-axis. This option is not valid for 3D projections.
+        Default is **True** for 1D and **False** for 2D to maintain old behaviour.
+    :type show_density: bool
 
     :keyword color: a color name or a list of color names or values, 
         default is ``'blue'``
@@ -246,12 +247,13 @@ def showProjection(ensemble, modes, *args, **kwargs):
                                 kwargs.pop('norm', False))
 
     if projection.ndim == 1 or projection.shape[1] == 1:
-        by_time = kwargs.pop('by_time', False)
+        by_time = not kwargs.pop('show_density', True)
+        by_time = kwargs.pop('by_time', by_time)
         if by_time:
             show = plt.plot(range(len(projection)), projection.flatten(), *args, **kwargs)
             plt.ylabel('Mode {0} coordinate'.format(str(modes)))
             plt.xlabel('Conformation number')  
-        else:          
+        else:
             show = plt.hist(projection.flatten(), *args, **kwargs)
             plt.xlabel('Mode {0} coordinate'.format(str(modes)))
             plt.ylabel('Number of conformations')
@@ -311,8 +313,17 @@ def showProjection(ensemble, modes, *args, **kwargs):
         indict[opts].append(i)
 
     modes = [m for m in modes]
-    if len(modes) == 2: 
-        plot = plt.plot
+    if len(modes) == 2:
+        show_density = kwargs.pop("show_density", False)
+        if show_density:
+            try:
+                import seaborn as sns
+                plot = sns.kdeplot
+                kwargs["cmap"] = cmap
+            except ImportError:
+                raise ImportError('Please install seaborn to plot kernel density estimates')
+        else:
+            plot = plt.plot
         show = plt.gcf()
         text = plt.text
     else: 

prody/proteins/ciffile.py

@@ -69,9 +69,14 @@ def parseMMCIF(pdb, **kwargs):
 
     :arg chain: comma separated string or list-like of chain IDs
     :type chain: str, tuple, list, :class:`~numpy.ndarray`
+
+    :arg unite_chains: unite chains with the same segment name
+        Default is *False*
+    :type unite_chains: bool
     """
     chain = kwargs.pop('chain', None)
     title = kwargs.get('title', None)
+    unite_chains = kwargs.get('unite_chains', False)
     auto_bonds = SETTINGS.get('auto_bonds')
     get_bonds = kwargs.get('bonds', auto_bonds)
     if get_bonds:
@@ -116,6 +121,8 @@ def parseMMCIF(pdb, **kwargs):
     cif = openFile(pdb, 'rt')
     result = parseMMCIFStream(cif, chain=chain, **kwargs)
     cif.close()
+    if unite_chains:
+        result.setSegnames(result.getChids())
     return result
 
 

prody/proteins/functions.py

@@ -80,7 +80,11 @@ def view3D(*alist, **kwargs):
         data_list = wrapModes(data_list)
         n_data = len(data_list)
 
-    view = kwargs.get('view',py3Dmol.view(width=width, height=height, js=kwargs.get('js','https://3dmol.csb.pitt.edu/build/3Dmol-min.js')))
+    view = kwargs.get('view', None)
+    js = kwargs.get('js', 'https://3dmol.csb.pitt.edu/build/3Dmol-min.js')
+
+    if view is None:
+        view = py3Dmol.view(width=width, height=height, js=js)
 
     def _mapData(atoms, data):
         # construct map from residue to data property
@@ -162,13 +166,13 @@ def _mapData(atoms, data):
     # setting styles ...
     view.setBackgroundColor(bgcolor)
 
-    if n_modes:
+    if n_modes and anim:
         #create vibrations
         view.vibrate(frames, scale)
 
         animate = kwargs.get('animate', {'loop':'rock', 'interval':interval})
         view.animate(animate)     
-        
+
     if isinstance(style, dict):
         style = ({}, style)
     if isinstance(style, tuple):