Source code for prody.dynamics.exanm

# -*- coding: utf-8 -*-
"""This module defines a class and a function for explicit membrane ANM calculations."""

import numpy as np

from prody import LOGGER, PY3K
from prody.atomic import Atomic, AtomGroup
from prody.utilities import importLA, checkCoords, copy
from numpy import sqrt, zeros, array, ceil, dot

from .anm import ANM
from .nma import MaskedNMA
from .gnm import checkENMParameters
from .editing import _reduceModel

LA = importLA()
inv = LA.inv
pinv = LA.pinv
norm = LA.norm

__all__ = ['exANM']

[docs]class exANM(ANM): """Class for explicit ANM (exANM) method ([FT00]_). Optional arguments build a membrane lattice permit analysis of membrane effect on elastic network models in *exANM* method described in [TL12]_. .. [TL12] Lezon TR, Bahar I, Constraints Imposed by the Membrane Selectively Guide the Alternating Access Dynamics of the Glutamate Transporter GltPh """ def __init__(self, name='Unknown'): super(exANM, self).__init__(name) self._membrane = None self._combined = None
[docs] def buildMembrane(self, coords, **kwargs): """Build membrane lattice around **coords**. :arg coords: a coordinate set or an object with ``getCoords`` method :type coords: :class:`numpy.ndarray` :arg membrane_high: the maximum z coordinate of the membrane. Default is **13.0** :type membrane_high: float :arg membrane_low: the minimum z coordinate of the membrane. Default is **-13.0** :type membrane_low: float :arg R: radius of all membrane in x-y direction. Default is **80** :type R: float :arg Ri: inner radius of the membrane in x-y direction if it needs to be hollow. Default is **0**, which is not hollow :type Ri: float :arg r: radius of each membrane node. Default is **3.1** :type r: float :arg lat: lattice type which could be **FCC** (face-centered-cubic, default), **SC** (simple cubic), **SH** (simple hexagonal) :type lat: str :arg exr: exclusive radius of each protein node. Default is **5.0** :type exr: float :arg hull: whether use convex hull to determine the protein's interior. Turn it off if protein is multimer. Default is **True** :type hull: bool :arg center: whether transform the structure to the origin (only x- and y-axis). Default is **True** :type center: bool """ atoms = coords try: coords = (coords._getCoords() if hasattr(coords, '_getCoords') else coords.getCoords()) except AttributeError: try: checkCoords(coords) except TypeError: raise TypeError('coords must be a Numpy array or an object ' 'with `getCoords` method') self._n_atoms = int(coords.shape[0]) LOGGER.timeit('_membrane') depth = kwargs.pop('depth', None) h = depth / 2 if depth is not None else None h = kwargs.pop('h', h) if h is not None: h = float(h) hu = h hl = -h else: hu = kwargs.pop('membrane_high', 13.0) hu = kwargs.pop('high', hu) hu = float(hu) hl = kwargs.pop('membrane_low', -13.0) hl = kwargs.pop('low', hl) hl = float(hl) R = float(kwargs.pop('R', 80.)) Ri = float(kwargs.pop('Ri', 0.)) r = float(kwargs.pop('r', 3.1)) lat = str(kwargs.pop('lat', 'FCC')) exr = float(kwargs.pop('exr', 5.)) use_hull = kwargs.pop('hull', True) centering = kwargs.pop('center', True) V = assign_lpvs(lat) if centering: c0 = coords.mean(axis=0) c0[-1] = 0. coords -= c0 # determine transmembrane part torf = np.logical_and(coords[:, -1] < hu, coords[:, -1] > hl) transmembrane = coords[torf, :] if not np.any(torf): raise ValueError('No region was identified as membrane. Please use a structure from opm/ppm.') if use_hull: from scipy.spatial import ConvexHull hull = ConvexHull(transmembrane) else: hull = transmembrane ## determine the bound for ijk imax = (R + V[0,2] * (hu - hl)/2.)/r jmax = (R + V[1,2] * (hu - hl)/2.)/r kmax = (R + V[2,2] * (hu - hl)/2.)/r imax = int(ceil(imax)) jmax = int(ceil(jmax)) kmax = int(ceil(kmax)) membrane = [] atm = 0 for i in range(-imax, imax): for j in range(-jmax, jmax): for k in range(-kmax, kmax): c = array([i, j, k]) xyz = 2.*r*dot(c, V) if xyz[2]>hl and xyz[2]<hu and \ xyz[0]>-R and xyz[0]<R and \ xyz[1]>-R and xyz[1]<R: dd = norm(xyz[:2]) if dd < R and dd > Ri: if checkClash(xyz, hull, radius=exr): membrane.append(xyz) atm = atm + 1 membrane = array(membrane) if len(membrane) == 0: self._membrane = None LOGGER.warn('no membrane is built. The protein should be transformed to the correct origin as in OPM') return coords else: self._membrane = AtomGroup(title="Membrane") self._membrane.setCoords(membrane) self._membrane.setResnums(range(atm)) self._membrane.setResnames(["NE1" for i in range(atm)]) self._membrane.setChids(["Q" for i in range(atm)]) self._membrane.setElements(["Q1" for i in range(atm)]) self._membrane.setNames(["Q1" for i in range(atm)]) LOGGER.report('Membrane was built in %2.fs.', label='_membrane') coords = self._combineMembraneProtein(atoms) return coords
[docs] def buildHessian(self, coords, cutoff=15., gamma=1., **kwargs): """Build Hessian matrix for given coordinate set. **kwargs** are passed to :method:`.buildMembrane`. :arg coords: a coordinate set or an object with ``getCoords`` method :type coords: :class:`numpy.ndarray` :arg cutoff: cutoff distance (Å) for pairwise interactions, default is 15.0 Å :type cutoff: float :arg gamma: spring constant, default is 1.0 :type gamma: float """ atoms = coords turbo = kwargs.pop('turbo', True) try: coords = (coords._getCoords() if hasattr(coords, '_getCoords') else coords.getCoords()) except AttributeError: try: checkCoords(coords) except TypeError: raise TypeError('coords must be a Numpy array or an object ' 'with `getCoords` method') n_atoms = int(coords.shape[0]) if self._membrane is None: coords = self.buildMembrane(atoms, **kwargs) else: coords = self._combined.getCoords() system = zeros(coords.shape[0], dtype=bool) system[:n_atoms] = True LOGGER.timeit('_exanm') if turbo: self._hessian = buildReducedHessian(coords, system, cutoff, gamma, **kwargs) else: super(exANM, self).buildHessian(coords, cutoff, gamma, **kwargs) system = np.repeat(system, 3) self._hessian = _reduceModel(self._hessian, system) LOGGER.report('Hessian was built in %.2fs.', label='_exanm') self._dof = self._hessian.shape[0] self._n_atoms = n_atoms
[docs] def calcModes(self, n_modes=20, zeros=False, turbo=True): """Calculate normal modes. This method uses :func:`scipy.linalg.eigh` function to diagonalize the Hessian matrix. When Scipy is not found, :func:`numpy.linalg.eigh` is used. :arg n_modes: number of non-zero eigenvalues/vectors to calculate. If **None** is given, all modes will be calculated. :type n_modes: int or None, default is 20 :arg zeros: If **True**, modes with zero eigenvalues will be kept. :type zeros: bool, default is **True** :arg turbo: Use a memory intensive, but faster way to calculate modes. :type turbo: bool, default is **True** """ super(exANM, self).calcModes(n_modes, zeros, turbo)
[docs] def getMembrane(self): """Returns a copy of the membrane coordinates.""" return self._membrane
[docs] def getCombined(self): """Returns a copy of the combined atoms or coordinates.""" return self._combined
def _combineMembraneProtein(self, coords): if self._membrane is not None: if isinstance(coords, Atomic): self._combined = coords.copy() + self._membrane coords = self._combined.getCoords() else: self._combined = coords = np.concatenate((coords, self._membrane.getCoords()), axis=0) else: self._combined = coords = copy(coords) return coords
def assign_lpvs(lat): """ Given lattice type return 3 lattice primitive vectors""" lpv = zeros((3,3)) if lat=='FCC': lpv[0,1]=1./sqrt(2) lpv[0,2]=1./sqrt(2) lpv[1,0]=1./sqrt(2) lpv[1,2]=1./sqrt(2) lpv[2,0]=1./sqrt(2) lpv[2,1]=1./sqrt(2) elif lat=='SC': lpv[0,0]=1 lpv[1,1]=1 lpv[2,2]=1 elif lat=='SH': lpv[0,0]=1./2 lpv[0,1]=-sqrt(3)/2 lpv[1,0]=1./2 lpv[1,1]=sqrt(3)/2 lpv[2,2]=1. return lpv def checkClash(node, hull, radius=5.): """ Check there is a clash between given coordinate and all pdb coordinates. **False** for clashing and **True** for not clashing.""" if isinstance(hull, np.ndarray): H = hull ishull = False else: H = hull.points ishull = True lb = H.min(axis=0) - radius ub = H.max(axis=0) + radius if np.all(node > ub) or np.all(node < lb): return True if ishull: in_hull = all([dot(eq[:-1], node) + eq[-1] <= 0 for eq in hull.equations]) if in_hull: return False for coord in H: if norm(node-coord)<radius: return False return True def peelr(coords, system, r0=20., dr=20.): n_sys_atoms = int(system.sum()) n_atoms = len(system) labels = np.zeros(n_atoms, dtype=int) # identify system beads sys_coords = coords[system, :2] sys_norms = norm(sys_coords, axis=1) sys_r = max(sys_norms) r0 += sys_r # label environment beads env_coords = coords[~system, :2] env_norms = norm(env_coords, axis=1) L = (env_norms - r0) // dr + 1 L = np.clip(L, 0, None) + 1 labels[n_sys_atoms:] = L uniq_labels = np.unique(labels) if len(uniq_labels) >= 3: uniq_labels.sort() lbl_last = uniq_labels[-1] lbl_2nd_last = uniq_labels[-2] n_last = np.sum(labels == lbl_last) n_2nd_last = np.sum(labels == lbl_2nd_last) if n_last < 0.2 * n_2nd_last: LOGGER.debug('edge nodes detected (%d/%d)'%(n_2nd_last, n_last)) labels[labels == lbl_last] = lbl_2nd_last if len(uniq_labels) >= 3: uniq_labels.sort() lbl_first = uniq_labels[1] lbl_2nd = uniq_labels[2] n_first = np.sum(labels == lbl_first) n_2nd = np.sum(labels == lbl_2nd) if n_first < 0.2 * n_2nd: LOGGER.debug('inner nodes detected (%d/%d)'%(n_2nd, n_first)) labels[labels == lbl_first] = lbl_2nd if not any(uniq_labels == 1): LOGGER.debug('no layer inside the system') for i in range(len(labels)): if labels[i] > 1: labels[i] -= 1 uniq_labels = np.unique(labels) for i, label in enumerate(uniq_labels): labels[labels==label] = i return labels def buildReducedHessian(coords, system, cutoff=15., gamma=1.0, **kwargs): r0 = kwargs.pop('r0', 20.) dr = kwargs.pop('dr', 20.) labels = peelr(coords, system, r0, dr) LOGGER.debug('layers: ' + str(np.unique(labels))) H = calcHessianRecursion(coords, labels, 0, cutoff=cutoff, gamma=gamma, **kwargs) return H def calcHessianRecursion(coords, layers, layer, cutoff=15., gamma=1.0, **kwargs): if layer == 0: LOGGER.debug('max layer: %d'%max(layers)) LOGGER.debug('layer: %d'%layer) Hss, Hse = buildLayerHessian(coords, layers, layer, cutoff=cutoff, gamma=gamma, **kwargs) if Hse is None: # last layer, Hee=Hss H = Hss else: Hee = calcHessianRecursion(coords, layers, layer+1, cutoff=cutoff, gamma=gamma, **kwargs) Cee = inv(Hee) #H = Hss - Hse.dot(Cee.dot(Hse.T)) #H = Hss - Hse @ Cee @ Hse.T if PY3K: H = Hss - Hse.__matmul__(Cee).__matmul__(Hse.T) else: H = Hss - Hse.dot(Cee.dot(Hse.T)) LOGGER.debug('layer: %d finished'%layer) return H def buildLayerHessian(coords, layers, layer, cutoff=15., gamma=1.0, **kwargs): gmem = kwargs.pop('gamma_memb', gamma) torf_inner = layers == (layer-1) torf_sys = layers == layer torf_env = layers == (layer+1) coords_inner = coords[torf_inner, :] # inner layer coords coords_sys = coords[torf_sys, :] # sys coords coords_env = coords[torf_env, :] # env coords n_sys_atoms = len(coords_sys) n_env_atoms = len(coords_env) n_inner_atoms = len(coords_inner) Hss = np.zeros((n_sys_atoms*3, n_sys_atoms*3)) Hse = np.zeros((n_sys_atoms*3, n_env_atoms*3)) if n_env_atoms else None cutoff2 = cutoff * cutoff for i in range(n_sys_atoms): coordi = coords_sys[i, :] I = slice(i*3, (i+1)*3) # sys-sys for j in range(i+1, n_sys_atoms): coordj = coords_sys[j, :] J = slice(j*3, (j+1)*3) v = coordi - coordj dist2 = np.inner(v, v) if dist2 < cutoff2: g = gamma if layer == 0 else gmem superelement = np.outer(v, v)*g/dist2 Hss[I, J] = Hss[J, I] = -superelement Hss[I, I] += superelement Hss[J, J] += superelement # sys-env for k in range(n_env_atoms): coordk = coords_env[k, :] K = slice(k*3, (k+1)*3) v = coordi - coordk dist2 = np.inner(v, v) if dist2 < cutoff2: g = gmem superelement = np.outer(v, v)*g/dist2 Hse[I, K] = -superelement Hss[I, I] += superelement # sys-inner for k in range(n_inner_atoms): coordk = coords_inner[k, :] K = slice(k*3, (k+1)*3) v = coordi - coordk dist2 = np.inner(v, v) if dist2 < cutoff2: g = gmem superelement = np.outer(v, v)*g/dist2 Hss[I, I] += superelement return Hss, Hse class MaskedExANM(exANM, MaskedNMA): def __init__(self, name='Unknown', mask=False, masked=True): exANM.__init__(self, name) MaskedNMA.__init__(self, name, mask, masked) def calcModes(self, n_modes=20, zeros=False, turbo=True): self._maskedarray = None super(MaskedExANM, self).calcModes(n_modes, zeros, turbo) def _reset(self): super(MaskedExANM, self)._reset() self._maskedarray = None def setHessian(self, hessian): """Set Hessian matrix. A symmetric matrix is expected, i.e. not a lower- or upper-triangular matrix.""" if not isinstance(hessian, np.ndarray): raise TypeError('hessian must be a Numpy array') elif hessian.ndim != 2 or hessian.shape[0] != hessian.shape[1]: raise ValueError('hessian must be square matrix') elif hessian.shape[0] % 3: raise ValueError('hessian.shape must be (3*n_atoms,3*n_atoms)') elif hessian.dtype != float: try: hessian = hessian.astype(float) except: raise ValueError('hessian.dtype must be float') mask = self.mask if not self.masked: if not np.isscalar(mask): if self.is3d(): mask = np.repeat(mask, 3) try: hessian = hessian[mask, :][:, mask] except IndexError: raise IndexError('size mismatch between Hessian (%d) and mask (%d).' 'Try set masked to False or reset mask'%(len(hessian), len(mask))) super(MaskedExANM, self).setHessian(hessian) def _getHessian(self): """Returns the Hessian matrix.""" hessian = self._hessian if hessian is None: return None if not self._isOriginal(): hessian = self._extend(hessian, axis=None) return hessian