Python molecule.Molecule类代码示例

    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)

    def testProteinLigand(self):
        from htmd.builder.solvate import solvate
        from htmd.parameterization.fftype import fftype, FFTypeMethod
        from htmd.parameterization.writers import writeFRCMOD

        # Test protein ligand building with parametrized ligand
        refdir = home(dataDir=join('test-amber-build', 'protLig'))
        tmpdir = os.path.join(self.testDir, 'protLig')
        os.makedirs(tmpdir)

        mol = Molecule(join(refdir, '3ptb_mod.pdb'))
        lig = Molecule(join(refdir, 'benzamidine.pdb'), guess=('bonds', 'angles', 'dihedrals'))
        prm, lig = fftype(lig, method=FFTypeMethod.GAFF2)
        writeFRCMOD(lig, prm, join(tmpdir, 'mol.frcmod'))
        lig.segid[:] = 'L'

        # params =
        newmol = Molecule()
        newmol.append(lig)
        newmol.append(mol)
        smol = solvate(newmol)

        params = defaultParam() + [join(tmpdir, 'mol.frcmod'),]

        _ = build(smol, outdir=tmpdir, param=params, ionize=False)

        refdir = home(dataDir=join('test-amber-build', 'protLig', 'results'))
        TestAmberBuild._compareResultFolders(refdir, tmpdir, '3PTB')

def _processSim(sim, projectionlist, uqmol, skip):
    pieces = sim.trajectory
    try:
        if uqmol is not None:
            mol = uqmol.copy()
        else:
            mol = Molecule(sim.molfile)
        logger.debug(pieces[0])

        mol.read(pieces, skip=skip)
        #Gianni testing
        #_highfreqFilter(mol,10)
 
        data = []
        for p in projectionlist:
            pj = _project(p, mol)
            if pj.ndim == 1:
                pj = np.atleast_2d(pj).T
            data.append(pj)
        data = np.hstack(data)
        if data.dtype == np.float64:
            data = data.astype(np.float32)
    except Exception as e:
        logger.warning('Error in simulation with id: ' + str(sim.simid) + '. "' + e.__str__() + '"')
        return None, None, None, True

    return data, _calcRef(pieces, mol.fileloc), mol.fstep, False

    def test_disulfideWithInsertion(self):
        from htmd.molecule.molecule import Molecule
        from htmd.builder.solvate import solvate
        from htmd.home import home
        from htmd.util import tempname, assertSameAsReferenceDir
        import os
        import numpy as np

        # Use pre-prepared files so we can tell whether the error is in prepare or in build
        # Inputs are reference outputs of proteinprepare.
        preparedInputDir = home(dataDir='test-proteinprepare')

        pdb = '3PTB'

        print('Building {}'.format(pdb))
        inFile = os.path.join(preparedInputDir, pdb, "{}-prepared.pdb".format(pdb))
        mol = Molecule(inFile)
        mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing

        np.random.seed(1)  # Needed for ions
        smol = solvate(mol)
        topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
        params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']

        smol.insertion[smol.resid == 42] = 'A'  # Adding an insertion to test that disulfide bonds with insertions work
        tmpdir = tempname()
        _ = build(smol, topo=topos, param=params, outdir=tmpdir)
        compareDir = home(dataDir=os.path.join('test-charmm-build', '3PTB_insertion'))
        assertSameAsReferenceDir(compareDir, tmpdir)

    def _prep_and_run(self, mol, rtf, prm, outdir, solvated):
        from htmd.builder.solvate import solvate
        from htmd.apps.acemdlocal import AcemdLocal
        # Do a simple solvation then run for 50ns
        ionize = True

        mol = Molecule(mol)
        mol.center()
        mol.set("segid", "L")
        d = maxDistance(mol, 'all') + 6

        if solvated:
            mol = solvate(mol, minmax=[[-d, -d, -d], [d, d, d]])

        if not solvated:
            ionize = False

        build_dir = os.path.join(outdir, "build")
        equil_dir = os.path.join(outdir, "equil")

        rtfs = ['top/top_water_ions.rtf', rtf]
        prms = ['par/par_water_ions.prm', prm]
        charmm.build(mol, topo=rtfs, param=prms, outdir=build_dir, ionize=ionize)
        md = Equilibration()
        md.runtime = 50
        md.timeunits = 'ns'
        md.temperature = 300
        md.write(build_dir, equil_dir)
        mdx = AcemdLocal()
        mdx.submit(equil_dir)
        mdx.wait()

    def test_build(self):
        from htmd.molecule.molecule import Molecule
        from htmd.builder.solvate import solvate
        from htmd.home import home
        from htmd.util import tempname, assertSameAsReferenceDir
        import os
        import numpy as np

        # Use pre-prepared files so we can tell whether the error is in prepare or in build
        # Inputs are reference outputs of proteinprepare.
        preparedInputDir = home(dataDir='test-proteinprepare')

        pdbids = ['3PTB', '1A25', '1GZM', '1U5U']
        for pdb in pdbids:
            with self.subTest(pdb=pdb):
                print('Building {}'.format(pdb))
                inFile = os.path.join(preparedInputDir, pdb, "{}-prepared.pdb".format(pdb))
                mol = Molecule(inFile)
                mol.filter('protein')  # Fix for bad proteinPrepare hydrogen placing

                np.random.seed(1)  # Needed for ions
                smol = solvate(mol)
                topos = ['top/top_all36_prot.rtf', 'top/top_water_ions.rtf']
                params = ['par/par_all36_prot_mod.prm', 'par/par_water_ions.prm']
                tmpdir = tempname()
                _ = build(smol, topo=topos, param=params, outdir=tmpdir)

                compareDir = home(dataDir=os.path.join('test-charmm-build', pdb))
                assertSameAsReferenceDir(compareDir, tmpdir)

    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)

def _filterPDBPSF(sim, outfolder, filtsel):
    try:
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise NameError('simFilter: ' + e.strerror + ' Cannot create filtered.pdb due to problematic pdb: ' + sim.molfile)

    if not path.isfile(path.join(outfolder, 'filtered.pdb')):
        mol.write(path.join(outfolder, 'filtered.pdb'), filtsel)

    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)

def _filterPDBPSF(sim, outfolder, filtsel):
    try:
        mol = Molecule(sim.molfile)
    except IOError as e:
        raise NameError('simFilter: {}. Cannot create filtered.pdb due to problematic pdb: {}'.format(e, sim.molfile))

    if not path.isfile(path.join(outfolder, 'filtered.pdb')):
        if mol.coords.size == 0:  # If we read for example psf or prmtop which have no coords, just add 0s everywhere
            mol.coords = np.zeros((mol.numAtoms, 3, 1), dtype=np.float32)
        mol.write(path.join(outfolder, 'filtered.pdb'), filtsel)

    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))

def _loadMolecules(lipids, files):
    from htmd.rotationmatrix import rotationMatrix
    # Create Molecules
    for l in lipids:
        randidx = np.random.randint(len(files[l.resname]))
        mol = Molecule(files[l.resname][randidx])
        mol.filter('not water', _logger=False)
        if l.xyz[2] < 0:
            mol.rotateBy(rotationMatrix([1, 0, 0], np.deg2rad(180)))  # Rotate the lower leaflet lipids upside down
        l.mol = mol
        l.rot = np.random.random() * 360 - 180  # Random starting rotation

    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)

def tileMembrane(memb, xmin, ymin, xmax, ymax, buffer=1.5):
    """ Tile a membrane in the X and Y dimensions to reach a specific size.

    Parameters
    ----------
    memb
    xmin
    ymin
    xmax
    ymax
    buffer

    Returns
    -------
    megamemb :
        A big membrane Molecule
    """
    from htmd.progress.progress import ProgressBar
    memb = memb.copy()
    memb.resid = sequenceID(memb.resid)

    minmemb = np.min(memb.get('coords', 'water'), axis=0).flatten()

    size = np.max(memb.get('coords', 'water'), axis=0) - np.min(memb.get('coords', 'water'), axis=0)
    size = size.flatten()
    xreps = int(np.ceil((xmax - xmin) / size[0]))
    yreps = int(np.ceil((ymax - ymin) / size[1]))

    logger.info('Replicating Membrane {}x{}'.format(xreps, yreps))

    from htmd.molecule.molecule import Molecule
    megamemb = Molecule()
    bar = ProgressBar(xreps * yreps, description='Replicating Membrane')
    k = 0
    for x in range(xreps):
        for y in range(yreps):
            tmpmemb = memb.copy()
            xpos = xmin + x * (size[0] + buffer)
            ypos = ymin + y * (size[1] + buffer)

            tmpmemb.moveBy([-float(minmemb[0]) + xpos, -float(minmemb[1]) + ypos, 0])
            tmpmemb.remove('same resid as (x > {} or y > {})'.format(xmax, ymax), _logger=False)
            tmpmemb.set('segid', 'M{}'.format(k))

            megamemb.append(tmpmemb)
            k += 1
            bar.progress()
    bar.stop()
    # Membranes don't tile perfectly. Need to remove waters that clash with lipids of other tiles
    # Some clashes will still occur between periodic images however
    megamemb.remove('same fragment as water and within 1.5 of not water', _logger=False)
    return megamemb

def opm(pdb):
    """Download a molecule from the OPM.

    Removes DUM atoms.

    Parameters
    ----------
    pdb: str
        The 4-letter PDB code

    Returns
    -------
    mol: Molecule
        The oriented molecule

    thickness: float
        The bilayer thickness (both layers)

    Examples
    --------
    >>> mol, thickness = opm("1z98")
    >>> mol.numAtoms
    7902
    >>> thickness
    28.2

    """

    from htmd.molecule.support import string_to_tempfile
    from htmd.molecule.molecule import Molecule
    # http://opm.phar.umich.edu/pdb/1z98.pdb
    r = requests.get("http://opm.phar.umich.edu/pdb/{:s}.pdb".format(pdb.lower()))

    if r.status_code != 200:
        raise NameError('PDB code not found in the OPM database')

    tempfile = string_to_tempfile(r.content.decode('ascii'), "pdb")
    mol = Molecule(tempfile)
    mol.filter("not resname DUM")

    # Assuming the half-thickness is the last word in the first line
    # REMARK      1/2 of bilayer thickness:   14.1
    f = open(tempfile)
    h = f.readline()
    f.close()
    os.unlink(tempfile)

    hs = h.split()
    thickness = 2.0 * float(hs[-1])

    return mol, thickness