Python Molecule.insertion方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import insertion [as 别名]
def _fillMolecule(name, resname, chain, resid, insertion, coords, segid, elements):
    numAtoms = len(name)
    mol = Molecule()
    mol.empty(numAtoms)

    mol.name = np.array(name, dtype=mol._pdb_fields['name'])
    mol.resname = np.array(resname, dtype=mol._pdb_fields['resname'])
    mol.chain = np.array(chain, dtype=mol._pdb_fields['chain'])
    mol.resid = np.array(resid, dtype=mol._pdb_fields['resid'])
    mol.insertion = np.array(insertion, dtype=mol._pdb_fields['insertion'])
    mol.coords = np.array(np.atleast_3d(np.vstack(coords)), dtype=mol._pdb_fields['coords'])
    mol.segid = np.array(segid, dtype=mol._pdb_fields['segid'])
    mol.element = np.array(elements, dtype=mol._pdb_fields['element'])
    return mol

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import insertion [as 别名]
def _fillMolecule(name, resname, chain, resid, insertion, coords, segid, element,
                  occupancy, beta, charge, record):
    numAtoms = len(name)
    mol = Molecule()
    mol.empty(numAtoms)

    mol.name = np.array(name, dtype=mol._dtypes['name'])
    mol.resname = np.array(resname, dtype=mol._dtypes['resname'])
    mol.chain = np.array(chain, dtype=mol._dtypes['chain'])
    mol.resid = np.array(resid, dtype=mol._dtypes['resid'])
    mol.insertion = np.array(insertion, dtype=mol._dtypes['insertion'])
    mol.coords = np.array(np.atleast_3d(np.vstack(coords)), dtype=mol._dtypes['coords'])
    mol.segid = np.array(segid, dtype=mol._dtypes['segid'])
    mol.element = np.array(element, dtype=mol._dtypes['element'])
    mol.occupancy = np.array(occupancy, dtype=mol._dtypes['occupancy'])
    mol.beta = np.array(beta, dtype=mol._dtypes['beta'])
    # mol.charge = np.array(charge, dtype=mol._dtypes['charge'])
    # mol.record = np.array(record, dtype=mol._dtypes['record'])
    return mol

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import insertion [as 别名]
def ionizePlace(mol, anion, cation, anionatom, cationatom, nanion, ncation, dfrom=5, dbetween=5, segname=None):
    newmol = mol.copy()

    logger.info('Min distance of ions from molecule: ' + str(dfrom) + 'A')
    logger.info('Min distance between ions: ' + str(dbetween) + 'A')
    logger.info('Placing ' + str(nanion+ncation) + ' ions.')

    if (nanion + ncation) == 0:
        return newmol

    segname = _getSegname(newmol, segname)
    nions = nanion + ncation

    betabackup = newmol.beta
    newmol.set('beta', sequenceID((newmol.resid, newmol.segid)))

    # Find water oxygens to replace with ions
    ntries = 0
    maxtries = 10
    while True:
        ionlist = np.empty(0, dtype=int)
        watindex = newmol.atomselect('noh and water and not (within ' + str(dfrom) + ' of not water)', indexes=True)
        watsize = len(watindex)

        if watsize == 0:
            raise NameError('No waters could be found further than ' + str(dfrom) + ' from other molecules to be replaced by ions. You might need to solvate with a bigger box.')

        while len(ionlist) < nions:
            if len(watindex) == 0:
                break
            randwat = np.random.randint(len(watindex))
            thision = watindex[randwat]
            addit = True
            if len(ionlist) != 0:  # Check for distance from precious ions
                ionspos = newmol.get('coords', sel=ionlist)
                thispos = newmol.get('coords', sel=thision)
                dists = distance.cdist(np.atleast_2d(ionspos), np.atleast_2d(thispos), metric='euclidean')

                if np.any(dists < dbetween):
                    addit = False
            if addit:
                ionlist = np.append(ionlist, thision)
                watindex = np.delete(watindex, randwat)
        if len(ionlist) == nions:
            break

        ntries += 1
        if ntries == maxtries:
            raise NameError('Failed to add ions after ' + str(maxtries) + ' attempts. Try decreasing the ''from'' and ''between'' parameters, decreasing ion concentration or making a larger water box.')

    # Delete waters but keep their coordinates
    waterpos = np.atleast_2d(newmol.get('coords', ionlist))
    stringsel = 'beta'
    for x in newmol.beta[ionlist]:
        stringsel += ' ' + str(int(x))
    atmrem = np.sum(newmol.atomselect(stringsel))
    atmput = 3 * len(ionlist)
    # assert atmrem == atmput, 'Removing {} atoms instead of {}. Report this bug.'.format(atmrem, atmput)
    sel = newmol.atomselect(stringsel, indexes=True)
    newmol.remove(sel, _logger=False)
    # assert np.size(sel) == atmput, 'Removed {} atoms instead of {}. Report this bug.'.format(np.size(sel), atmput)
    betabackup = np.delete(betabackup, sel)

    # Add the ions
    randidx = np.random.permutation(np.size(waterpos, 0))
    atom = Molecule()
    atom.record = 'ATOM'
    atom.chain = 'I'
    atom.segid = 'I'
    atom.occupancy = 0
    atom.beta = 0
    atom.insertion = ''
    atom.element = ''
    atom.altloc = ''

    for i in range(nanion):
        atom.name = anionatom
        atom.resname = anion
        atom.resid = newmol.resid[-1] + 1
        atom.coords = waterpos[randidx[i], :]
        newmol.insert(atom, len(newmol.name))
    for i in range(ncation):
        atom.name = cationatom
        atom.resname = cation
        atom.resid = newmol.resid[-1] + 1
        atom.coords = waterpos[randidx[i+nanion], :]
        newmol.insert(atom, len(newmol.name))

    # Restoring the original betas
    sel = np.ones(len(betabackup) + nions, dtype=bool)
    sel[len(betabackup)::] = False
    newmol.set('beta', betabackup, sel=sel)
    return newmol