Python Molecule.atomselect方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def test_tmscore(self):
        from htmd.molecule.molecule import Molecule
        expectedTMscore = np.array([0.21418524, 0.2367377, 0.23433833, 0.21362964, 0.20935164,
                                    0.20279461, 0.27012895, 0.22675238, 0.21230793, 0.2372011])
        expectedRMSD = np.array([3.70322128, 3.43637027, 3.188193, 3.84455877, 3.53053882,
                                 3.46781854, 2.93777629, 2.97978692, 2.70792428, 2.63051318])

        mol = Molecule(os.path.join(home(dataDir='tmscore'), 'filtered.pdb'))
        mol.read(os.path.join(home(dataDir='tmscore'), 'traj.xtc'))
        ref = Molecule(os.path.join(home(dataDir='tmscore'), 'ntl9_2hbb.pdb'))
        tmscore, rmsd = molTMscore(mol, ref, mol.atomselect('protein'), ref.atomselect('protein'))

        self.assertTrue(np.allclose(tmscore, expectedTMscore))
        self.assertTrue(np.allclose(rmsd, expectedRMSD))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _listXTCs(path.join(outFolder, name))
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

    try:
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning(e.strerror + ', skipping trajectory')
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _listXTCs(path.join(outFolder, name))
    #bar.progress()
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def __init__(self, sims, refmol, trajalnstr, refalnstr, atomsel, centerstr):
        self._refmol = refmol
        self._refalnsel = self._refmol.atomselect(refalnstr)
        self._trajalnsel = trajalnstr
        self._centersel = centerstr
        self._atomsel = atomsel
        self._pc_trajalnsel = None  # pc = Pre-calculated
        self._pc_atomsel = None
        self._pc_centersel = None

        (single, molfile) = _singleMolfile(sims)
        if single:
            mol = Molecule(molfile)
            self._pc_trajalnsel = mol.atomselect(trajalnstr)
            self._pc_atomsel = mol.atomselect(atomsel)
            self._pc_centersel = mol.atomselect(centerstr)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def __init__(self, sims, sel, simple):
        self._pc_sel = None
        self._sel = sel
        self._simple = simple

        (single, molfile) = _singleMolfile(sims)
        if single:
            mol = Molecule(molfile)
            self._pc_sel = mol.atomselect(sel)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def __init__(self, sims, protsel, dih=None, sincos=True):
        self._protsel = protsel
        self._sincos = sincos
        self._dih = dih  # TODO: Calculate the dihedral
        self._pc_dih = None

        (single, molfile) = _singleMolfile(sims)
        if single:
            mol = Molecule(molfile)
            self._pc_dih = self._dihedralPrecalc(mol, mol.atomselect(protsel))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
        numframes = _getNumFrames(sims[i], ftrajectory)
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

    try:
        from htmd.molecule.molecule import Molecule
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning('{}, skipping trajectory'.format(e))
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
    numframes = _getNumFrames(sims[i], ftrajectory)
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def write(self, inputdir, outputdir):
        """ Writes the production protocol and files into a folder.

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a equilibration process.
        outputdir : str
            Directory where to write the production setup files.
        """
        from htmd.molecule.molecule import Molecule

        # Do version consistency check
        if (self._version == 2 and not isinstance(self.acemd, Acemd2)) and \
                (self._version == 3 and not isinstance(self.acemd, Acemd)):
            raise RuntimeError('Acemd object version ({}) inconsistent with protocol version at instantiation '
                               '({})'.format(type(self.acemd), self._version))

        self._findFiles(inputdir)
        self._amberFixes()

        if self._version == 2:
            self.acemd.temperature = str(self.temperature)
            self.acemd.langevintemp = str(self.temperature)
        elif self._version == 3:
            self.acemd.temperature = self.temperature
            self.acemd.thermostattemp = self.temperature

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        if self._version == 3:
            self.acemd.run = str(numsteps)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self._version == 2:
            if self.restraints:
                raise RuntimeWarning('restraints are only available on {}(_version=3)'.format(self.__class__.__name__))
            if self.fb_k > 0:  # use TCL only for flatbottom
                self.acemd.tclforces = 'on'
                if isinstance(self.acemd.TCL, tuple):
                    tcl = list(self.acemd.TCL)
                    tcl[0] = tcl[0].format(NUMSTEPS=numsteps, TEMPERATURE=self.temperature, KCONST=self.fb_k,
                                           REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                           SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                           BOX=' '.join(map(str, self.fb_box)))
                    self.acemd.TCL = tcl[0] + tcl[1]
                else:
                    logger.warning('{} default TCL was already formatted.'.format(self.__class__.__name__))
            else:
                self.acemd.TCL = 'set numsteps {NUMSTEPS}\n'.format(NUMSTEPS=numsteps)
            if self.useconstraints:
                # Turn on constraints
                self.acemd.constraints = 'on'
                self.acemd.constraintscaling = '1.0'
            else:
                if len(self.constraints) != 0:
                    logger.warning('You have setup constraints to {} but constraints are turned off. '
                                   'If you want to use constraints, define '
                                   'useconstraints=True'.format(self.constraints))
        elif self._version == 3:
            if self.restraints is not None:
                logger.info('Using user-provided restraints and ignoring constraints and fb_potential')
                self.acemd.restraints = self.restraints
            else:
                restraints = list()
                if self.fb_k > 0:
                    logger.warning('Converting fb_potential to restraints. This is a convenience '
                                   'functional conversion. We recommend start using restraints with '
                                   '{}(_version=3)'.format(self.__class__.__name__))
                    restraints += self._fb_potential2restraints(inputdir)
                if self.useconstraints:
                    logger.warning('Converting constraints to restraints. This is a convenience '
                                   'functional conversion. We recommend start using restraints with '
                                   '{}(_version=3)'.format(self.__class__.__name__))
                    restraints += self._constraints2restraints()
                else:
                    if len(self.constraints) != 0:
                        logger.warning('You have setup constraints to {} but constraints are turned off. '
                                       'If you want to use constraints, define '
                                       'useconstraints=True'.format(self.constraints))
                if len(restraints) != 0:
                    self.acemd.restraints = restraints

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        if self.adaptive:
            self.acemd.binvelocities = None

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        if self._version == 2:
            # Adding constraints by writing them to the consref file
            if self.useconstraints:
#.........这里部分代码省略.........

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]

#.........这里部分代码省略.........
        param = []
    if caps is None:
        caps = _defaultCaps(mol)

    _missingSegID(mol)
    _checkMixedSegment(mol)

    logger.info('Converting CHARMM membranes to AMBER.')
    mol = _charmmLipid2Amber(mol)

    #_checkProteinGaps(mol)
    _applyCaps(mol, caps)

    f = open(path.join(outdir, 'tleap.in'), 'w')
    f.write('# tleap file generated by amber.build\n')

    # Printing out the forcefields
    for force in ff:
        f.write('source ' + force + '\n')
    f.write('\n')

    # Loading TIP3P water parameters
    f.write('# Loading ions and TIP3P water parameters\n')
    f.write('loadamberparams frcmod.ionsjc_tip3p\n\n')

    # Printing out topologies
    logger.info('Writing prepi files.')
    f.write('# Loading prepi topologies\n')
    for t in topo:
        shutil.copy(t, outdir)
        f.write('loadamberprep ' + path.basename(t) + '\n')
    f.write('\n')

    # Printing and loading the PDB file. AMBER can work with a single PDB file if the segments are separate by TER
    logger.info('Writing PDB file for input to tleap.')
    pdbname = path.join(outdir, 'input.pdb')
    mol.write(pdbname)
    if not os.path.isfile(pdbname):
        raise NameError('Could not write a PDB file out of the given Molecule.')
    f.write('# Loading the system\n')
    f.write('mol = loadpdb input.pdb\n\n')

    # Printing out patches for the disulfide bridges
    '''if disulfide is None and not ionize:
        logger.info('Detecting disulfide bonds.')
        disulfide = detectDisulfideBonds(mol)

    if not ionize and len(disulfide) != 0:  # Only make disu bonds after ionizing!
        f.write('# Adding disulfide bonds\n')
        for d in disulfide:
            # Convert to stupid amber residue numbering
            uqseqid = sequenceID(mol.resid)
            uqres1 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid1, d.resid1))]))
            uqres2 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid2, d.resid2))]))
            # Rename the CYS to CYX if there is a disulfide bond
            mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid1, d.resid1))
            mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid2, d.resid2))
            f.write('bond mol.{}.SG mol.{}.SG\n'.format(uqres1, uqres2))
        f.write('\n')'''

    f.write('# Writing out the results\n')
    f.write('savepdb mol ' + prefix + '.pdb\n')
    f.write('saveamberparm mol ' + prefix + '.prmtop ' + prefix + '.crd\n')
    f.write('quit')
    f.close()

    molbuilt = None
    if execute:
        logpath = os.path.abspath('{}/log.txt'.format(outdir))
        logger.info('Starting the build.')
        currdir = os.getcwd()
        os.chdir(outdir)
        f = open(logpath, 'w')
        try:
            call([tleap, '-f', './tleap.in'], stdout=f)
        except:
            raise NameError('tleap failed at execution')
        f.close()
        os.chdir(currdir)
        logger.info('Finished building.')

        if path.getsize(path.join(outdir, 'structure.pdb')) != 0 and path.getsize(path.join(outdir, 'structure.prmtop')) != 0:
            molbuilt = Molecule(path.join(outdir, 'structure.pdb'))
            molbuilt.read(path.join(outdir, 'structure.prmtop'))
            molbuilt.bonds = []  # Causes problems in ionization mol.remove and mol._removeBonds
        else:
            raise NameError('No structure pdb/prmtop file was generated. Check {} for errors in building.'.format(logpath))

        if ionize:
            shutil.move(path.join(outdir, 'structure.pdb'), path.join(outdir, 'structure.noions.pdb'))
            shutil.move(path.join(outdir, 'structure.crd'), path.join(outdir, 'structure.noions.crd'))
            shutil.move(path.join(outdir, 'structure.prmtop'), path.join(outdir, 'structure.noions.prmtop'))
            totalcharge = np.sum(molbuilt.charge)
            nwater = np.sum(molbuilt.atomselect('water and noh'))
            anion, cation, anionatom, cationatom, nanion, ncation = ionizef(totalcharge, nwater, saltconc=saltconc, ff='amber', anion=saltanion, cation=saltcation)
            newmol = ionizePlace(molbuilt, anion, cation, anionatom, cationatom, nanion, ncation)
            # Redo the whole build but now with ions included
            return build(newmol, ff=ff, topo=topo, param=param, prefix=prefix, outdir=outdir, caps={}, ionize=False,
                         execute=execute, saltconc=saltconc, disulfide=disulfide, tleap=tleap)
    return molbuilt

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def write(self, inputdir, outputdir):
        """ Write the equilibration protocol

        Writes the equilibration protocol and files into a folder for execution
        using files inside the inputdir directory

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a build process.
        outputdir : str
            Directory where to write the equilibration setup files.

        Examples
        --------
        >>> md = Equilibration()
        >>> md.write('./build','./equil')
        """
        self._findFiles(inputdir)
        self._amberFixes()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self.constraintsteps is None:
            constrsteps = int(numsteps / 2)
        else:
            constrsteps = int(self.constraintsteps)

        if isinstance(self.acemd.TCL, tuple):
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                   BOX=' '.join(map(str, self.fb_box)),
                                   NVTSTEPS=self.nvtsteps, CONSTRAINTSTEPS=constrsteps, TEMPERATURE=self.temperature)
            self.acemd.TCL = tcl[0] + tcl[1]
        else:
            logger.warning('{} default TCL was already formatted.'.format(self.__class__.__name__))

        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            coords = inmol.get('coords', sel='water')
            if coords.size == 0:  # It's a vacuum simulation
                coords = inmol.get('coords', sel='all')
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
                dim = dim + 12.
            else:
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
            self.acemd.celldimension = '{} {} {}'.format(dim[0], dim[1], dim[2])

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints by writing them to the consref file
        inconsreffile = os.path.join(inputdir, self.acemd.consref)
        consrefmol = Molecule(inconsreffile)
        consrefmol.set('occupancy', 0)
        consrefmol.set('beta', 0)
        if len(self.constraints) == 0:
            raise RuntimeError('You have not defined any constraints for the Equilibration (constraints={}).')
        else:
            for sel in self.constraints:
                consrefmol.set('beta', self.constraints[sel], sel)
        outconsreffile = os.path.join(outputdir, self.acemd.consref)
        consrefmol.write(outconsreffile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]

#.........这里部分代码省略.........
    f.write('\n')

    # Loading TIP3P water parameters
    f.write('# Loading ions and TIP3P water parameters\n')
    f.write('loadamberparams frcmod.ionsjc_tip3p\n\n')

    # Loading user parameters
    f.write('# Loading parameter files\n')
    for p in param:
        try:
            shutil.copy(p, outdir)
            f.write('loadamberparams ' + path.basename(p) + '\n')
        except:
            f.write('loadamberparams ' + p + '\n')
            logger.info("Path {:s} not found, assuming a standard AmberTools file.".
                        format(p))
    f.write('\n')

    # Printing out topologies
    f.write('# Loading prepi topologies\n')
    for t in topo:
        shutil.copy(t, outdir)
        f.write('loadamberprep ' + path.basename(t) + '\n')
    f.write('\n')

    # Detect disulfide bonds
    if disulfide is None and not ionize:
        logger.info('Detecting disulfide bonds.')
        disulfide = detectDisulfideBonds(mol)
        if len(disulfide) != 0:
            for d in disulfide:
                # Convert to stupid amber residue numbering
                uqseqid = sequenceID((mol.resid, mol.insertion, mol.segid)) + mol.resid[0] - 1
                uqres1 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid1, d.resid1))]))
                uqres2 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid2, d.resid2))]))
                # Rename the CYS to CYX if there is a disulfide bond
                mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid1, d.resid1))
                mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid2, d.resid2))
                # Remove (eventual) HG hydrogens on these CYS (from proteinPrepare)
                mol.remove('name HG and segid {} and resid {}'.format(d.segid1, d.resid1), _logger=False)
                mol.remove('name HG and segid {} and resid {}'.format(d.segid2, d.resid2), _logger=False)

    # Printing and loading the PDB file. AMBER can work with a single PDB file if the segments are separate by TER
    logger.info('Writing PDB file for input to tleap.')
    pdbname = path.join(outdir, 'input.pdb')

    # mol2 files have atomtype, here we only write parts not coming from mol2
    mol.write(pdbname, mol.atomtype == '')
    if not os.path.isfile(pdbname):
        raise NameError('Could not write a PDB file out of the given Molecule.')
    f.write('# Loading the system\n')
    f.write('mol = loadpdb input.pdb\n\n')

    if np.sum(mol.atomtype != '') != 0:
        logger.info('Writing mol2 files for input to tleap.')
        segs = np.unique(mol.segid[mol.atomtype != ''])
        combstr = 'mol = combine {mol'
        for s in segs:
            name = 'segment{}'.format(s)
            mol2name = path.join(outdir, '{}.mol2'.format(name))
            mol.write(mol2name, (mol.atomtype != '') & (mol.segid == s))
            if not os.path.isfile(mol2name):
                raise NameError('Could not write a mol2 file out of the given Molecule.')
            f.write('# Loading the rest of the system\n')
            f.write('{} = loadmol2 {}.mol2\n\n'.format(name, name))
            combstr += ' {}'.format(name)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [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 atomselect [as 别名]
    def write(self, inputdir, outputdir):
        """ Writes the production protocol and files into a folder.

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a equilibration process.
        outputdir : str
            Directory where to write the production setup files.
        """
        self._findFiles(inputdir)
        self._amberFixes()

        self.acemd.temperature = self.temperature
        self.acemd.langevintemp = self.temperature

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self.fb_k > 0:  # use TCL only for flatbottom
            self.acemd.tclforces = 'on'
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, TEMPERATURE=self.temperature, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                   BOX=' '.join(map(str, self.fb_box)))
            self.acemd.TCL = tcl[0] + tcl[1]
        else:
            self.acemd.TCL = 'set numsteps {NUMSTEPS}\n' \
                             'set temperature {TEMPERATURE}\n'.format(NUMSTEPS=numsteps, TEMPERATURE=self.temperature)

        if self.useconstraints:
            # Turn on constraints
            self.acemd.constraints = 'on'
            self.acemd.constraintscaling = 1.0
        else:
            if len(self.constraints) != 0:
                logger.warning('You have setup constraints to {} but constraints are turned off. '
                               'If you want to use constraints, define useconstraints=True'.format(self.constraints))

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        if self.adaptive:
            self.acemd.binvelocities = None

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints by writing them to the consref file
        if self.useconstraints:
            inconsreffile = os.path.join(inputdir, self.acemd.consref)
            consrefmol = Molecule(inconsreffile)
            consrefmol.set('occupancy', 0)
            consrefmol.set('beta', 0)
            if len(self.constraints) == 0:
                raise RuntimeError('You have set the production to use constraints (useconstraints=True), but have not '
                                   'defined any constraints (constraints={}).')
            else:
                for sel in self.constraints:
                    consrefmol.set('beta', self.constraints[sel], sel)
            outconsreffile = os.path.join(outputdir, self.acemd.consref)
            consrefmol.write(outconsreffile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    return angles, dihedrals


# A test method
if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    from htmd.home import home
    import numpy as np
    from os import path
    import doctest
    #doctest.testmod(extraglobs={"mol" : Molecule("3PTB")})

    expectedTMscore = np.array([ 0.21418524,  0.2367377 ,  0.23433833,  0.21362964,  0.20935164,
        0.20279461,  0.27012895,  0.22675238,  0.21230793,  0.2372011 ])
    expectedRMSD = np.array([ 3.70322128,  3.43637027,  3.188193  ,  3.84455877,  3.53053882,
        3.46781854,  2.93777629,  2.97978692,  2.70792428,  2.63051318])

    mol = Molecule(path.join(home(), 'data', 'tmscore', 'filtered.pdb'))
    mol.read(path.join(home(), 'data', 'tmscore', 'traj.xtc'))
    ref = Molecule(path.join(home(), 'data', 'tmscore', 'ntl9_2hbb.pdb'))
    tmscore, rmsd = molTMscore(mol, ref, mol.atomselect('protein'), ref.atomselect('protein'))

    assert np.allclose(tmscore, expectedTMscore)
    assert np.allclose(rmsd, expectedRMSD)

    #
    # rhodopsin = Molecule('1F88')
    # d3r = Molecule('3PBL')
    # alnmol = sequenceStructureAlignment(rhodopsin, d3r)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]

#.........这里部分代码省略.........
    #_checkProteinGaps(mol)
    if patches is None:
        patches = []
    if isinstance(patches, str):
        patches = [patches]
    allpatches = []
    allpatches += patches
    # Find protonated residues and add patches for them
    allpatches += _protonationPatches(mol)

    f = open(path.join(outdir, 'build.vmd'), 'w')
    f.write('# psfgen file generated by charmm.build\n')
    f.write('package require psfgen;\n')
    f.write('psfcontext reset;\n\n')

    # Copying and printing out the topologies
    if not path.exists(path.join(outdir, 'topologies')):
        os.makedirs(path.join(outdir, 'topologies'))
    for i in range(len(alltopo)):
        if alltopo[i][0] != '.' and path.isfile(path.join(charmmdir, alltopo[i])):
            alltopo[i] = path.join(charmmdir, alltopo[i])
        localname = '{}.'.format(i) + path.basename(alltopo[i])
        shutil.copy(alltopo[i], path.join(outdir, 'topologies', localname))
        f.write('topology ' + path.join('topologies', localname) + '\n')
    f.write('\n')

    _printAliases(f)

    # Printing out segments
    if not path.exists(path.join(outdir, 'segments')):
        os.makedirs(path.join(outdir, 'segments'))
    logger.info('Writing out segments.')
    segments = _getSegments(mol)
    wateratoms = mol.atomselect('water')
    for seg in segments:
        pdbname = 'segment' + seg + '.pdb'
        segatoms = mol.segid == seg
        mol.write(path.join(outdir, 'segments', pdbname), sel=segatoms)

        segwater = wateratoms & segatoms
        f.write('segment ' + seg + ' {\n')
        if np.all(segatoms == segwater):  # If segment only contains waters, set: auto none
            f.write('\tauto none\n')
        f.write('\tpdb ' + path.join('segments', pdbname) + '\n')
        if caps is not None and seg in caps:
            for c in caps[seg]:
                f.write('\t' + c + '\n')
        f.write('}\n')
        f.write('coordpdb ' + path.join('segments', pdbname) + ' ' + seg + '\n\n')

    # Printing out patches for the disulfide bridges
    if disulfide is None:
        disulfide = detectDisulfideBonds(mol)

    if len(disulfide) != 0:
        for d in disulfide:
            f.write('patch DISU {}:{} {}:{}\n'.format(d.segid1, d.resid1, d.segid2, d.resid2))
        f.write('\n')

    noregenpatches = [p for p in allpatches if p.split()[1] in noregen]
    regenpatches = [p for p in allpatches if p.split()[1] not in noregen]

    # Printing regenerable patches
    if len(regenpatches) != 0:
        for p in regenpatches:
            f.write(p + '\n')

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import atomselect [as 别名]
    def write(self, inputdir, outputdir):
        """ Write the equilibration protocol

        Writes the equilibration protocol and files into a folder for execution
        using files inside the inputdir directory

        Parameters
        ----------
        inputdir : str
            Path to a directory containing the files produced by a build process.
        outputdir : str
            Directory where to write the equilibration setup files.

        Examples
        --------
        >>> md = Equilibration()
        >>> md.write('./build','./equil')
        """

        from htmd.molecule.molecule import Molecule

        # Do version consistency check
        if (self._version == 2 and not isinstance(self.acemd, Acemd2)) and \
                (self._version == 3 and not isinstance(self.acemd, Acemd)):
            raise RuntimeError('Acemd object version ({}) inconsistent with protocol version at instantiation '
                               '({})'.format(type(self.acemd), self._version))

        self._findFiles(inputdir)
        self._amberFixes()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        if self._version == 3:
            self.acemd.temperature = self.temperature
            self.acemd.thermostattemp = self.temperature
            self.acemd.run = str(numsteps)

        pdbfile = os.path.join(inputdir, self.acemd.coordinates)
        inmol = Molecule(pdbfile)

        if np.any(inmol.atomselect('lipids')) and not self.useconstantratio:
            logger.warning('Lipids detected in input structure. We highly recommend setting useconstantratio=True '
                           'for membrane simulations.')

        if self.constraintsteps is None:
            constrsteps = int(numsteps / 2)
        else:
            constrsteps = int(self.constraintsteps)

        if self._version == 2:
            if self.restraints:
                raise RuntimeWarning('restraints are only available on {}(_version=3)'.format(self.__class__.__name__))
            if isinstance(self.acemd.TCL, tuple):
                tcl = list(self.acemd.TCL)
                tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                       REFINDEX=' '.join(map(str, inmol.get('index', self.fb_reference))),
                                       SELINDEX=' '.join(map(str, inmol.get('index', self.fb_selection))),
                                       BOX=' '.join(map(str, self.fb_box)),
                                       NVTSTEPS=self.nvtsteps, CONSTRAINTSTEPS=constrsteps, TEMPERATURE=self.temperature)
                self.acemd.TCL = tcl[0] + tcl[1]
            else:
                logger.warning('{} default TCL was already formatted.'.format(self.__class__.__name__))
        elif self._version == 3:
            if self.restraints is not None:
                logger.info('Using user-provided restraints and ignoring constraints and fb_potential')
                self.acemd.restraints = self.restraints
            else:
                logger.warning('Converting constraints and fb_potential to restraints. This is a convenience '
                               'functional conversion. We recommend start using restraints with '
                               '{}(_version=3)'.format(self.__class__.__name__))
                restraints = list()
                # convert constraints to restraints and add them
                if self.constraints is not None:
                    restraints += self._constraints2restraints(constrsteps)
                # convert fb_potential to restraints and add them
                if self.fb_k > 0:
                    restraints += self._fb_potential2restraints(inputdir)
                self.acemd.restraints = restraints

        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            coords = inmol.get('coords', sel='water')
            if coords.size == 0:  # It's a vacuum simulation
                coords = inmol.get('coords', sel='all')
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
                dim += 12.
            else:
                dim = np.max(coords, axis=0) - np.min(coords, axis=0)
            self.acemd.celldimension = '{} {} {}'.format(dim[0], dim[1], dim[2])

        if self.useconstantratio:
            self.acemd.useconstantratio = 'on'

        self.acemd.setup(inputdir, outputdir, overwrite=True)

        if self._version == 2:
            # Adding constraints by writing them to the consref file
            inconsreffile = os.path.join(inputdir, self.acemd.consref)
            consrefmol = Molecule(inconsreffile)
            consrefmol.set('occupancy', 0)
            consrefmol.set('beta', 0)
#.........这里部分代码省略.........