Python Molecule.copy方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
    def test_project_align(self):
        from htmd.molecule.molecule import Molecule
        from htmd.home import home
        from os import path
        mol = Molecule(path.join(home(), 'data', 'metricdistance', 'filtered.pdb'))
        mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))
        ref = mol.copy()
        ref.coords = np.atleast_3d(ref.coords[:, :, 0])
        metr = MetricCoordinate('protein and name CA', ref)
        data = metr.project(mol)

        lastcoors = np.array([6.79283285, 5.55226946, 4.49387407, 2.94484425,
                              5.36937141, 3.18590879, 5.75874281, 5.48864174,
                              1.69625032, 1.58790839, 0.57877392, -2.66498065,
                              -3.70919156, -3.33702421, -5.38465405, -8.43286991,
                              -8.15859032, -7.85062265, -10.92551327, -13.70733166], dtype=np.float32)
        assert np.all(np.abs(data[-1, -20:] - lastcoors) < 0.001), 'Coordinates calculation is broken'

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
    def test_project(self):
        from htmd.molecule.molecule import Molecule
        from htmd.home import home
        from os import path
        mol = Molecule(path.join(home(), 'data', 'metricdistance', 'filtered.pdb'))
        mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))
        ref = mol.copy()
        ref.coords = np.atleast_3d(ref.coords[:, :, 0])
        metr = MetricCoordinate('protein and name CA')
        data = metr.project(mol)

        lastcoors = np.array([-24.77000237, -27.76000023, -30.44000244, -33.65000153,
                              -33.40999985, -36.32000351, -36.02000427, -36.38000107,
                              -39.61000061, -41.01000214, -43.80000305, -45.56000137,
                              -45.36000061, -47.13000488, -49.54000473, -50.6000061 ,
                              -50.11000061, -52.15999985, -55.1400032 , -55.73000336], dtype=np.float32)
        assert np.all(np.abs(data[-1, -20:] - lastcoors) < 0.001), 'Coordinates calculation is broken'

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
            raise NameError('Atom selection resulted in 0 atoms.')
        if np.sum(centersel) == 0:
            raise NameError('Centering selection resulted in 0 atoms.')
        return trajalnsel, atomsel, centersel

    def _getMapping(self, mol):
        (xxx, atomsel, yyy) = self._getSelections(mol)
        premap = np.where(atomsel)[0]
        map = np.zeros(len(premap) * 3, dtype=int)
        map[0*len(premap):1*len(premap)] = premap
        map[1*len(premap):2*len(premap)] = premap
        map[2*len(premap):3*len(premap)] = premap


if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    from htmd.home import home
    from os import path
    mol = Molecule(path.join(home(), 'data', 'metricdistance', 'filtered.pdb'))
    mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))
    ref = mol.copy()
    ref.coords = np.atleast_3d(ref.coords[:, :, 0])
    metr = MetricCoordinate(ref, 'protein and name CA')
    data = metr.project(mol)

    lastcoors = np.array([6.79283285,   5.55226946,   4.49387407,   2.94484425,
                          5.36937141,   3.18590879,   5.75874281,   5.48864174,
                          1.69625032,   1.58790839,   0.57877392,  -2.66498065,
                          -3.70919156,  -3.33702421,  -5.38465405,  -8.43286991,
                          -8.15859032,  -7.85062265, -10.92551327, -13.70733166], dtype=np.float32)
    assert np.all(np.abs(data[-1, -20:] - lastcoors) < 0.001), 'Coordinates calculation is broken'

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
            A Molecule object which will be used to calculate the descriptions of the projected dimensions.

        Returns
        -------
        map : :class:`DataFrame <pandas.core.frame.DataFrame>` object
            A DataFrame containing the descriptions of each dimension
        """
        trajalnsel, targetcom, refcom, centersel = self._getSelections(mol)
        targetatomidx = np.where(targetcom)[0]
        refatomidx = np.where(refcom)[0]
        from pandas import DataFrame
        types = ['r', 'theta', 'phi']
        indexes = [[targetatomidx, refatomidx]] * 3
        description = ['r', 'theta', 'phi']
        return DataFrame({'type': types, 'atomIndexes': indexes, 'description': description})


if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    from htmd.home import home
    from os import path
    ref = Molecule(path.join(home(dataDir='metricdistance'), 'filtered.pdb'))
    mol = ref.copy()
    mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))

    res = MetricSphericalCoordinate(ref, 'resname MOL', 'within 8 of resid 98').project(mol)
    _ = MetricSphericalCoordinate(ref, 'resname MOL', 'within 8 of resid 98').getMapping(mol)

    ref_array = np.load(path.join(home(dataDir='test-metrics'), 'metricsphericalcoordinate', 'res.npy'))
    assert np.allclose(res, ref_array, rtol=0, atol=1e-04)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
class MutualInformation:
    def __init__(self, model, mol=None, fstep=0.1, skip=1):
        """ Class that calculates the mutual information of protein residues.

        Parameters
        ----------
        model : :class:`Model <htmd.model.Model>` object
            A Model object with a calculated MSM
        mol : :class:`Molecule <htmd.molecule.molecule.Molecule>` object
            A reference molecule from which to obtain structural information. By default model.data.simlist[0].molfile
            will be used.
        fstep : float
            The frame step of the simulations
        skip : int
            Frame skipping

        Examples
        --------
        >>> from htmd.mutualinformation import MutualInformation
        >>> from htmd.ui import *
        >>>
        >>> sims = simlist(glob('./filtered/*/'), './filtered/filtered.pdb')
        >>> mol = Molecule('./filtered/filtered.pdb')
        >>>
        >>> metr = Metric(sims)
        >>> metr.set(MetricDihedral())
        >>> datadih = metr.project()
        >>> datadih.fstep = 0.1
        >>> datadih.dropTraj()
        >>>
        >>> metr = Metric(datadih.simlist)
        >>> metr.set(MetricSelfDistance('protein and name CA', metric='contacts'))
        >>> dataco = metr.project()
        >>> dataco.fstep = 0.1
        >>> dataco.dropTraj()
        >>>
        >>> tica = TICA(datadih, 20, units='ns')
        >>> datatica = tica.project(3)
        >>> datatica.cluster(MiniBatchKMeans(n_clusters=1500))
        >>> model = Model(datatica)
        >>> model.markovModel(12, 4, units='ns')
        >>>
        >>> mu = MutualInformation(model)
        >>> mu.calculate()
        >>> mu.saveMI('./mi_matrix.npy')
        >>> mu.weightGraph(dataco, 0.005)
        >>> mu.save_graphml('./weightgraph.graphml')
        """
        self.model = model

        self.mol = mol
        if mol is None:
            self.mol = Molecule(self.model.data.simlist[0].molfile)

        self._computeChiDihedrals(fstep=fstep, skip=skip)

        self.bindihcat = self._histogram()  # Lasts two minutes
        self.resids = self.mol.get('resid', 'name CA')
        self.residmap = np.ones(self.mol.resid.max() + 1, dtype=int) * -1
        self.residmap[self.resids] = np.arange(len(self.resids))

        self.mi_matrix = None
        self.graph_array = None
        self.graph = None


    def calculate(self):
        from htmd.config import _config
        from htmd.parallelprogress import ParallelExecutor
        numchi = self.chi.numDimensions
        statdist = self.model.msm.stationary_distribution
        stconcat = np.concatenate(self.model.data.St)
        microcat = self.model.micro_ofcluster[stconcat]

        aprun = ParallelExecutor(n_jobs=_config['ncpus'])
        res = aprun(total=numchi, desc='Calculating MI')(delayed(self._parallelAll)(numchi, dih1, 4, self.model.micronum, self.bindihcat, microcat, statdist) for dih1 in range(numchi))
        MI_all = np.zeros((len(self.resids), len(self.resids)))
        for r in res:
            dihcounts = r[0]
            pairs = r[1]
            for dihc, p in zip(dihcounts, pairs):
                dih1, dih2 = p
                if dih1 == dih2:
                    continue
                resid1 = self.residmap[self.mol.resid[self.chi.description.atomIndexes[dih1][0]]]
                resid2 = self.residmap[self.mol.resid[self.chi.description.atomIndexes[dih2][0]]]
                MI_all[resid1][resid2] = self._calcMutualInfo(dihc)
        self.mi_matrix = self._cleanautocorrelations(MI_all)

    def _computeChiDihedrals(self, fstep=0.1, skip=1):
        chis = []
        protmol = self.mol.copy()
        protmol.filter('protein')
        caidx = self.mol.atomselect('protein and name CA')
        resids = self.mol.resid[caidx]
        resnames = self.mol.resname[caidx]
        for residue, resname in zip(resids, resnames):
            ch = Dihedral.chi1(protmol, residue)
            if ch is not None:
                chis.append(ch)
#.........这里部分代码省略.........

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import copy [as 别名]
def solvate(mol, pad=None, minmax=None, negx=0, posx=0, negy=0, posy=0, negz=0, posz=0, buffer=2.4, watsize=65.4195,
            prefix='WT', keysel='name OH2', rotate=False, rotsel='all', rotinc=36, spdb=None,
            spsf=None, stop=None):
    """ Solvates the system in a water box


    Parameters
    ----------
    mol : :class:`Molecule <htmd.molecule.molecule.Molecule>` object
        The molecule object we want to solvate
    pad : float
        The padding to add to the minmax in all dimensions. You can specify different padding in each dimension using
        the negx, negy, negz, posx, posy, posz options. This option will override any values in the neg and pos options.
    minmax : list
        Min and max dimensions. Should be a 2D matrix of the form [[minx, miny, minz], [maxx, maxy, maxz]]. If none is
        given, it is calculated from the minimum and maximum coordinates in the mol.
    negx : float
        The padding in the -x dimension
    posx : float
        The padding in the +x dimension
    negy : float
        The padding in the -y dimension
    posy : float
        The padding in the +y dimension
    negz : float
        The padding in the -z dimension
    posz : float
        The padding in the +z dimension
    buffer : float
        How much buffer space to leave empty between waters and other molecules
    watsize : float
        The size of the water box
    prefix : str
        The prefix used for water segments
    keysel : str
        The key selection for water atoms
    rotate : bool
        Enable automated rotation of molecule to fit best in box
    rotsel : str
        The selection of atoms to rotate
    rotinc : float
        The increment in degrees to rotate
    spdb : str
        The path to the water pdb file
    spsf : str
        The path to the water psf file
    stop : str
        The path to the water topology file

    Returns
    -------
    mol : :class:`Molecule <htmd.molecule.molecule.Molecule>` object
        A solvated molecule

    Examples
    --------
    >>> smol = solvate(mol, pad=10)
    >>> smol = solvate(mol, minmax=[[-20, -20, -20],[20, 20, 20]])
    """
    mol = mol.copy()
    if mol.numFrames > 1:
        logger.warning('Multiple frames in Molecule. Solvate keeps only frame 0 and discards the rest.')
        mol.coords = np.atleast_3d(mol.coords[:, :, 0])

    if spdb is None:
        spdb = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'wat.pdb')

    if os.path.isfile(spdb):
        logger.info('Using water pdb file at: ' + spdb)
        water = Molecule(spdb)
    else:
        raise NameError('No solvent pdb file found in ' + spdb)

    if pad is not None:
        negx = pad; posx = pad; negy = pad; posy = pad; negz = pad; posz = pad

    if rotate:
        raise NameError('Rotation not implemented yet')

    # Calculate min max coordinates from molecule
    if mol.numAtoms > 0:
        minmol = np.min(mol.get('coords'), axis=0)
        maxmol = np.max(mol.get('coords'), axis=0)
    else:
        minmol = [np.inf, np.inf, np.inf]
        maxmol = [-np.inf, -np.inf, -np.inf]

    if minmax is None:
        minc = minmol
        maxc = maxmol
    else:
        if isinstance(minmax, list):
            minmax = np.array(minmax)
        minc = minmax[0, :]
        maxc = minmax[1, :]

    xmin = float(minc[0] - negx)
    xmax = float(maxc[0] + posx)
    ymin = float(minc[1] - negy)
    ymax = float(maxc[1] + posy)
#.........这里部分代码省略.........