Python Molecule.get方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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 self.constraintsteps is None:
            constrsteps = int(numsteps / 2)
        else:
            constrsteps = int(self.constraintsteps)

        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]

        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
        inmol.set('occupancy', 0)
        inmol.set('beta', 0)
        for sel in self.constraints:
            inmol.set('beta', self.constraints[sel], sel)
        outfile = os.path.join(outputdir, self.acemd.coordinates)
        inmol.write(outfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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)
        self.acemd.temperature = self.temperature
        self.acemd.thermostattemp = self.temperature
        self.acemd.run = str(numsteps)

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

        # Adding the default restraints of the equilibration
        restraints = list()
        restraints.append(AtomRestraint('protein and noh and not name CA', 0, [(0.1, 0), (0, constrsteps)]))
        restraints.append(AtomRestraint('protein and name CA', 0, [(1, 0), (0, constrsteps)]))
        if self.restraints is not None:
            restraints += self.restraints
        self.acemd.restraints = restraints

        if self.acemd.celldimension is None and self.acemd.extendedsystem is None:
            inmol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            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)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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()

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

        self.acemd.TCL = self.acemd.TCL.replace('NUMSTEPS', str(self.numsteps))
        self.acemd.TCL = self.acemd.TCL.replace('TEMPERATURE', str(self.temperature))
        self.acemd.TCL = self.acemd.TCL.replace('KCONST', str(self.k))
        self.acemd.TCL = self.acemd.TCL.replace('REFINDEX', ' '.join(map(str, inmol.get('index', self.reference))))
        self.acemd.TCL = self.acemd.TCL.replace('SELINDEX', ' '.join(map(str, inmol.get('index', self.selection))))
        self.acemd.TCL = self.acemd.TCL.replace('BOX', ' '.join(map(str, self.box)))
        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
        inmol.set('occupancy', 0)
        inmol.set('beta', 0)
        for sel in self.constraints:
            inmol.set('beta', self.constraints[sel], sel)
        outfile = os.path.join(outputdir, self.acemd.coordinates)
        inmol.write(outfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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 = str(self.temperature)
        self.acemd.langevintemp = str(self.temperature)
        if self.k > 0: #use TCL only for flatbottom
            mol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            self.acemd.tclforces = 'on'
            TCL = self._TCL
            TCL = TCL.replace('KCONST', str(self.k))
            TCL = TCL.replace('REFINDEX', ' '.join(map(str, mol.get('index', self.reference))))
            TCL = TCL.replace('SELINDEX', ' '.join(map(str, mol.get('index', self.selection))))
            TCL = TCL.replace('BOX', ' '.join(map(str, self.box)))
            self.acemd.TCL = TCL
        else:
            self.acemd.TCL = ''
        self.acemd.setup(inputdir, outputdir, overwrite=True)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        self.acemd.temperature = str(self.temperature)
        self.acemd.langevintemp = str(self.temperature)
        if self.fb_k > 0: #use TCL only for flatbottom
            mol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            self.acemd.tclforces = 'on'
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, mol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, mol.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 {}\n'.format(numsteps)
        self.acemd.setup(inputdir, outputdir, overwrite=True)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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()

        from htmd.units import convert
        numsteps = convert(self.timeunits, 'timesteps', self.runtime, timestep=self.acemd.timestep)
        self.acemd.temperature = str(self.temperature)
        self.acemd.langevintemp = str(self.temperature)
        if self.fb_k > 0: #use TCL only for flatbottom
            mol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            self.acemd.tclforces = 'on'
            tcl = list(self.acemd.TCL)
            tcl[0] = tcl[0].format(NUMSTEPS=numsteps, KCONST=self.fb_k,
                                   REFINDEX=' '.join(map(str, mol.get('index', self.fb_reference))),
                                   SELINDEX=' '.join(map(str, mol.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 {}\n'.format(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))
        self.acemd.setup(inputdir, outputdir, overwrite=True)

        # Adding constraints
        if self.useconstraints:
            inmol = Molecule(os.path.join(inputdir, self.acemd.coordinates))
            inmol.set('occupancy', 0)
            inmol.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:
                    inmol.set('beta', self.constraints[sel], sel)
                outfile = os.path.join(outputdir, self.acemd.coordinates)
                inmol.write(outfile)
                self.acemd.consref = self.acemd.coordinates

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [as 别名]
def equilibrateSystem(pdbfile, psffile, outpdb, numsteps=30000, minimplatform='CPU', equilplatform='CUDA', device=0,
                      temp=300, minimize=500000, minimizetol=100, charmmfolder=None):
    pdb = Molecule(pdbfile)
    watcoo = pdb.get('coords', 'water')
    celld = unit.Quantity((watcoo.max(axis=0) - watcoo.min(axis=0)).squeeze(), unit=unit.angstrom)

    psf = app.CharmmPsfFile(psffile)
    pdb = app.PDBFile(pdbfile)
    psf.setBox(celld[0], celld[1], celld[2])

    params = readCharmmParameters(charmmfolder, defaultCharmmFiles)

    system = psf.createSystem(params, nonbondedMethod=app.PME,
                              nonbondedCutoff=1*unit.nanometer,
                              constraints=app.HBonds)
    system.addForce(mm.MonteCarloBarostat(1*unit.atmospheres, temp*unit.kelvin, 25))

    platform, prop = getPlatform(minimplatform, device)
    integrator = mm.LangevinIntegrator(temp*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds)
    simulation = app.Simulation(psf.topology, system, integrator, platform, prop)
    simulation.context.setPositions(pdb.positions)

    # Perform minimization
    printEnergies(simulation, 'Energy before minimization')
    simulation.minimizeEnergy(tolerance=minimizetol*unit.kilojoule/unit.mole, maxIterations=minimize)
    printEnergies(simulation, 'Energy after minimization')

    # Copy coordinates from miminization context to equilibration context
    state = simulation.context.getState(getPositions=True)
    pos = state.getPositions()

    # Set up the equilibration simulation
    platform, prop = getPlatform(equilplatform, device)
    integrator = mm.LangevinIntegrator(temp*unit.kelvin, 1/unit.picosecond, 0.002*unit.picoseconds)
    simulation = app.Simulation(psf.topology, system, integrator, platform, prop)
    simulation.context.setPositions(pos)

    # Generate initial velocities
    simulation.context.setVelocitiesToTemperature(temp)

    from htmd.membranebuilder.pdbreporter import PDBReporter
    from htmd.membranebuilder.dcdreporter import DCDReporter
    simulation.reporters.append(PDBReporter(outpdb, numsteps, enforcePeriodicBox=False))
    simulation.reporters.append(app.StateDataReporter(stdout, int(numsteps/10), step=True,
        potentialEnergy=True, kineticEnergy=True, totalEnergy=True, temperature=True,
        progress=True, speed=True, remainingTime=True, totalSteps=numsteps, separator='\t'))
    simulation.step(numsteps)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [as 别名]
    def _fb_potential2restraints(self, inputdir):
        from htmd.molecule.molecule import Molecule
        restraints = list()

        fb_box = np.array(self.fb_box)
        # convert fb_box to width
        width = list(np.concatenate(np.diff(np.array([fb_box[::2], fb_box[1::2]]), axis=0)))

        # If fb_box is not symmetrical
        if not np.all(fb_box[::2] == -fb_box[1::2]):
            # convert fb_box and fb_reference to fbcentre and width
            mol = Molecule(os.path.join(inputdir, self.acemd.structure))
            mol.read(os.path.join(inputdir, self.acemd.coordinates))
            fb_refcentre = mol.get('coords', sel=self.fb_reference).mean(axis=0).squeeze()

            fbcentre = list(np.around(np.mean(np.array([fb_box[::2], fb_box[1::2]]), axis=0) + fb_refcentre, 3))
            restraints.append(GroupRestraint(self.fb_selection, width, [(self.fb_k, 0)], fbcentre=fbcentre))
        else:
            restraints.append(GroupRestraint(self.fb_selection, width, [(self.fb_k, 0)], fbcentresel=self.fb_reference))

        return restraints

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [as 别名]
        # Cartesian to polar coordinates. Rho is the rotation angle
        rho = np.sqrt(x ** 2 + y ** 2)
        phi = np.arctan2(y, x)
        return rho, phi

    angle, _ = cart2pol(xa, ya)
    return angle + np.radians(45)


if __name__ == "__main__":
    from htmd import *
    from os import path

    p = Molecule(path.join(home(), 'data', 'building-protein-membrane', '4dkl.pdb'))
    p.filter('(chain B and protein) or water')
    p = autoSegment(p, 'protein', 'P')
    m = Molecule(path.join(home(), 'data', 'building-protein-membrane', 'membrane.pdb'))
    a = embed(p, m)
    print(np.unique(m.get('segid')))

    mol = Molecule(path.join(home(), 'data', 'building-protein-membrane', '1ITG_clean.pdb'))
    ref = Molecule(path.join(home(), 'data', 'building-protein-membrane', '1ITG.pdb'))
    mol = autoSegment(mol, sel='protein')
    assert np.all(mol.segid == ref.segid)

    mol = Molecule(path.join(home(), 'data', 'building-protein-membrane', '3PTB_clean.pdb'))
    ref = Molecule(path.join(home(), 'data', 'building-protein-membrane', '3PTB.pdb'))
    mol = autoSegment(mol, sel='protein')
    assert np.all(mol.segid == ref.segid)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import get [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 get [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 get [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 get [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)
#.........这里部分代码省略.........