Source code for prody.dynamics.exgnm

# -*- coding: utf-8 -*-
"""This module defines a class and a function for explicit membrane GNM 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 .gnm import GNM, checkENMParameters
from .editing import _reduceModel

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

__all__ = ['exGNM']

[docs]class exGNM(GNM): """Class for explicit GNM (exGNM) method ([FT00]_). Optional arguments build a membrane lattice permit analysis of membrane effect on elastic network models in *exGNM* 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(exGNM, self).__init__(name) self._membrane = None self._combined = None
[docs] def buildMembrane(self, coords, **kwargs): """Build Kirchhoff matrix for given coordinate set. :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 *.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 buildKirchhoff(self, coords, cutoff=15., gamma=1., **kwargs): """Build Kirchhoff 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('_exgnm') if turbo: self._kirchhoff = buildReducedKirchhoff(coords, system, cutoff, gamma, **kwargs) else: super(exGNM, self).buildKirchhoff(coords, cutoff, gamma, **kwargs) system = np.repeat(system, 3) self._kirchhoff = _reduceModel(self._kirchhoff, system) LOGGER.report('Kirchhoff was built in %.2fs.', label='_exgnm') self._dof = self._kirchhoff.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 Kirchhoff 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(exGNM, 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 buildReducedKirchhoff(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))) G = calcKirchhoffRecursion(coords, labels, 0, cutoff=cutoff, gamma=gamma, **kwargs) return G def calcKirchhoffRecursion(coords, layers, layer, cutoff=15., gamma=1.0, **kwargs): if layer == 0: LOGGER.debug('max layer: %d'%max(layers)) LOGGER.debug('layer: %d'%layer) Gss, Gse = buildLayerKirchhoff(coords, layers, layer, cutoff=cutoff, gamma=gamma, **kwargs) if Gse is None: # last layer, Gee=Gss G = Gss else: Gee = calcKirchhoffRecursion(coords, layers, layer+1, cutoff=cutoff, gamma=gamma, **kwargs) Cee = inv(Gee) #G = Gss - Gse.dot(Cee.dot(Gse.T)) #G = Gss - Gse @ Cee @ Gse.T if PY3K: G = Gss - Gse.__matmul__(Cee).__matmul__(Gse.T) else: G = Gss - Gse.dot(Cee.dot(Gse.T)) LOGGER.debug('layer: %d finished'%layer) return G def buildLayerKirchhoff(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) Gss = np.zeros((n_sys_atoms, n_sys_atoms)) Gse = np.zeros((n_sys_atoms, n_env_atoms)) if n_env_atoms else None cutoff2 = cutoff * cutoff for i in range(n_sys_atoms): coordi = coords_sys[i, :] I = slice(i, (i+1)) # sys-sys for j in range(i+1, n_sys_atoms): coordj = coords_sys[j, :] J = slice(j, (j+1)) v = coordi - coordj dist2 = np.inner(v, v) if dist2 < cutoff2: g = gamma if layer == 0 else gmem Gss[I, J] = Gss[J, I] = -g Gss[I, I] += g Gss[J, J] += g # sys-env for k in range(n_env_atoms): coordk = coords_env[k, :] K = slice(k, (k+1)) v = coordi - coordk dist2 = np.inner(v, v) if dist2 < cutoff2: g = gmem Gse[I, K] = -g Gss[I, I] += g # sys-inner for k in range(n_inner_atoms): coordk = coords_inner[k, :] K = slice(k, (k+1)) v = coordi - coordk dist2 = np.inner(v, v) if dist2 < cutoff2: g = gmem Gss[I, I] += g return Gss, Gse