Python Universe.selectAtoms方法代碼示例

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def count_interactions(grof, xtcf, btime, cutoff, debug):
    u = Universe(grof, xtcf)
    un_query = ('(resname PRO and (name CB or name CG or name CD)) or'
                '(resname VAL and (name CG1 or name CG2)) or'
                '(resname GLY and name CA) or'
                '(resname ALA and name CB)')
    vp_query = ('name OW')
    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs the unit is nm
    un_atoms = u.selectAtoms(un_query)
    for ts in u.trajectory:
        if ts.time >= btime:
            numcount = 0
            tropo_vp_atoms = u.selectAtoms(
                '({0}) and around 8 ({1})'.format(vp_query, un_query))
            # different from when calculating unun, there is no overlap atom
            # between un_atoms & tropo_vp_atoms
            for ai in un_atoms:
                for aj in tropo_vp_atoms:
                    d = np.linalg.norm(ai.pos - aj.pos)
                    if d <= cutoff:
                        numcount += 1
            yield '{0:10.0f}{1:8d}\n'.format(ts.time,  numcount)
        # per 100 frames, num of frames changes with the size of xtc file, for debugging
        if debug and ts.frame % 2 == 0: 
            print "time: {0:10.0f}; step: {1:10d}; frame: {2:10d}".format(ts.time, ts.step, ts.frame)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
 def test_write_selection(self):
     ref = Universe(mol2_molecule)
     gr0 = ref.selectAtoms("name C*")
     gr0.write(self.outfile)
     u = Universe(self.outfile)
     gr1 = u.selectAtoms("name C*")
     assert_equal(len(gr0), len(gr1))

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def calc_rama(grof, xtcf, btime, etime):
    u = Universe(grof, xtcf)

    resname_query = 'resname GLY or resname VAL or resname PRO'
    atoms = u.selectAtoms(resname_query)
    resname = atoms.resnames()[0] # [0] because .resnames() returns a list of one element
    resid = atoms.resids()[0] # [0] because .resnames() returns a list of one element

    phi_query = ('(resname ACE and name C) or '
                 '(resname GLY or resname VAL or resname PRO and '
                 '(name N or name CA or name C))')

    psi_query = ('(resname GLY or resname VAL or resname PRO and (name N or name CA or name C or name NT)) or '
                 '(resname NH2 and name N)')

    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs the unit is nm
    phi = u.selectAtoms(phi_query)
    psi = u.selectAtoms(psi_query)

    for _ in phi.atoms:
        print _

    for _ in psi.atoms:
        print _


    for ts in u.trajectory:
        if btime > ts.time:
            continue
        if etime > 0 and etime < ts.time:
            break

        yield '{0:.3f}  {1:.3f}  {2}-{3}\n'.format(
            phi.dihedral(), psi.dihedral(), resname, resid)
        U.print_progress(ts)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def count_interactions(A):
    logger.debug('loading {0}'.format(A.grof))
    univ = Universe(A.grof)
    logger.debug('loaded {0}'.format(A.grof))

    pro_atoms = univ.selectAtoms('protein and not resname ACE and not resname NH2')
    pl = pro_atoms.residues.numberOfResidues()
    # +1: for missing resname ACE, such that it's easier to proceed in the next
    # step

    logger.debug('loading {0}, {1}'.format(A.grof, A.xtcf))
    u = Universe(A.grof, A.xtcf)
    logger.debug('loaded {0}, {1}'.format(A.grof, A.xtcf))

    # Just for reference to the content of query when then code was first
    # written and used
    # query = ('(resname PRO and (name CB or name CG or name CD)) or'
    #          '(resname VAL and (name CG1 or name CG2)) or'
    #          '(resname GLY and name CA) or'
    #          '(resname ALA and name CB)')

    query = A.query
    atoms = u.selectAtoms(query)
    logger.info('Number of atoms selected: {0}'.format(atoms.numberOfAtoms()))

    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs
    # the unit is nm
    cutoff = A.cutoff * 10
    nres_away = A.nres_away
    btime = A.btime
    etime = A.etime
    nframe = 0
    unun_map = None
    for ts in u.trajectory:
        if btime > ts.time:
            continue
        if etime > 0 and etime < ts.time:
            break

        nframe += 1
        map_ = np.zeros((pl+1, pl+1))                   # map for a single frame
        for i, ai in enumerate(atoms):
            ai_resid = ai.resid
            for j, aj in enumerate(atoms):
                aj_resid = aj.resid
                # to avoid counting the same pair twices,
                # the 2 resid cannot be neigbors
                if i < j and aj_resid - ai_resid >= nres_away:
                    d = np.linalg.norm(ai.pos - aj.pos)
                    if d <= cutoff:
                        # -1: resid in MDAnalysis starts from 1
                        map_[ai_resid-1][aj_resid-1] += 1
        if unun_map is None:
            unun_map = map_
        else:
            unun_map = unun_map + map_
        utils.print_progress(ts)
    sys.stdout.write("\n")
    return unun_map / float(nframe)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def gen_hbond_map(xpm, ndx, grof):
    xpm = objs.XPM(xpm)
    hbndx = objs.HBNdx(ndx)

    univ = Universe(grof)
    pro_atoms = univ.selectAtoms('protein and not resname ACE and not resname NH2')
    hbonds_by_resid = hbndx.map_id2resid(pro_atoms)

    # pl: peptide length
    pl = pro_atoms.residues.numberOfResidues()

    hblist = []
    for i, j in zip(hbonds_by_resid, xpm.color_count):
        # j[1] is the probability of hbonds, while j[0] = 1 - j[1]
        # format: [resid of donor, resid of acceptor]
        # -1 is because resid in MDAnalysis starts from 1, minus so as to fit
        # -into hb_map initialized by hb_map
        hblist.append([i[0]-1, i[1]-1, j[1]])

    # +1: for missing resname ACE, such that it's easier to proceed in the next
    # step
    pl1 = pl + 1
    hb_map = np.zeros((pl1, pl1))
    for _ in hblist:
        hb_map[_[0]][_[1]] = _[2]

    return hb_map

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def main():
    arg_parser = argparse.ArgumentParser(description='通過給定殘基名稱，殘基內原子數目，兩個原子在殘基內的索引(從0開始)，計算所有殘基內這兩個原子之間的直線距離。')
    arg_parser.add_argument('resname', action='store', help='殘基名稱')
    arg_parser.add_argument('atoms_num', type=int, action='store', help='殘基內原子數目')
    arg_parser.add_argument('index1', type=int, action='store', help='第一個原子的索引，索引從0開始')
    arg_parser.add_argument('index2', type=int, action='store', help='第二個原子的索引，索引從0開始')
    arg_parser.add_argument('topology_file', action='store', help='拓撲文件，例如gro, pdb')
    args = arg_parser.parse_args()

    resname, atoms_num, index1, index2 = args.resname, args.atoms_num, args.index1, args.index2

    universe = Universe(args.topology_file)
    atom_groups = universe.selectAtoms("resname " + resname)
    if len(atom_groups) % atoms_num != 0:
        print("拓撲文件內對應殘基原子總數不是所給原子數目的整數倍，請給予正確的原子數目。")
        exit(1)

    atoms1 = []
    atoms2 = []
    for i in range(0, len(atom_groups), atoms_num):
        atoms1.append(atom_groups[i:i + atoms_num][index1])
        atoms2.append(atom_groups[i:i + atoms_num][index2])

    dists = dist(AtomGroup(atoms1), AtomGroup(atoms2))
    print("The distance between atoms %s and %s is:" % (index1, index2))
    for i in dists[2]:
        print(i)
    print("The average distance between atoms %s and %s is:" % (index1, index2))
    print(np.average(dists[2]))

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def count_interactions(grof, xtcf, btime, etime, cutoff):
    cutoff = cutoff * 10 # * 10: convert from nm to angstrom to work with MDAnalysis
    u = Universe(grof, xtcf)
    query = ('(resname PRO and (name CB or name CG or name CD)) or'
             '(resname VAL and (name CG1 or name CG2)) or'
             '(resname GLY and name CA) or'
             '(resname ALA and name CB)')
    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs the unit is nm
    atoms = u.selectAtoms(query)
    for ts in u.trajectory:
        if btime > ts.time:
            continue
        if etime > 0 and etime < ts.time:
            break

        numcount = 0
        for i, ai in enumerate(atoms):
            for j, aj in enumerate(atoms):
                # to avoid counting the same pair twices,
                # the 2 resid cannot be neigbors
                if i < j and abs(ai.resid - aj.resid) >= 2: 
                    d = np.linalg.norm(ai.pos - aj.pos)
                    if d <= cutoff:
                        numcount += 1
        yield '{0:10.0f}{1:8d}\n'.format(ts.time,  numcount)
        utils.print_progress(ts)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def main(struct):
    u = Universe(struct)

    phi = u.selectAtoms(PHI_SEL)
    psi = u.selectAtoms(PSI_SEL)
    
    print u.filename
    print 'phi: {0:8.2f}'.format(phi.dihedral())
    print 'psi: {0:8.2f}'.format(psi.dihedral())
    print 

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def calc_rg(grof, xtcf, btime, debug):
    u = Universe(grof, xtcf)
    query = 'name CA'
    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs the unit is nm
    atoms = u.selectAtoms(query)
    natoms = atoms.numberOfAtoms()
    for ts in u.trajectory:
        if ts.time >= btime:
            com = atoms.centerOfMass()                                # center of mass
            _sum = sum((sum(i**2 for i in (a.pos - com)) for a in atoms))
            rg = np.sqrt(_sum / natoms)
            yield '{0:10.0f}{1:15.6f}\n'.format(ts.time, rg)
        # per 100 frames, num of frames changes with the size of xtc file, for debugging
        if debug and ts.frame % 2 == 0: 
            print "time: {0:10.0f}; step: {1:10d}; frame: {2:10d}".format(ts.time, ts.step, ts.frame)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def sequence_spacing(grof, xtcf, btime, etime, peptide_length, atom_sel):
    u = Universe(grof, xtcf)
    # this selection part should be better customized
    # here, only have been backbone atoms are used, u.selectAtoms doesn't
    # include Hydrogen atoms
    # REMMEMBER: ARGS verification should be done in main ONLY!
    # range works like this:

    # in MDAnalysis, resid starts from 1, in sequence_spacing.py, we don't count
    # the C- and N- termini, so it's from 2 to peptide_len+2
    residues = [u.selectAtoms(atom_sel.format(i)) for i in range(2, peptide_length + 2)]
    ijdist_dict = {}
    for ts in u.trajectory:
        # btime, etime defaults to 0, if etime is 0, loop till the end of the
        # trajectory
        if btime > ts.time:
            continue
        if etime > 0 and etime < ts.time:
            break

        # the good stuff
        for i, resi in enumerate(residues):
            for j, resj in enumerate(residues):
                # to remove duplicate since resi & resj are within the same peptide 
                if i < j:
                    dij = abs(i - j)
                    d_atomi_atomj = []
                    # loop through every atom in both residues
                    for atomi in resi:
                        for atomj in resj:
                            d_atomi_atomj.append(
                                np.linalg.norm(atomi.pos - atomj.pos))
                # add the result to the dictionary
                    ij_dist = np.average(d_atomi_atomj)   # distance between i and j
                    if dij not in ijdist_dict.keys():
                        ijdist_dict[dij] = [ij_dist]
                    else:
                        ijdist_dict[dij].append(ij_dist)
        utils.print_progress(ts)

    return ijdist_dict

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def getWaterCoorWithH(self,centre,psf,dcd,outputFile):
    rho=Universe(psf,dcd)
    H2OCoordinate=[]
    no=0
    title='resname'+'    '+'atomid'+'    '+'resnumber'+'    X    Y     Z   '+'   '+'segname'+'  '+'frameNo'+'   '+'centreNo'+'\n'
    outputFile.write(title)
    for oxygenInforSet in self:
        H2OCoordinateSet=[]
        
        
        print 'There were',len(oxygenInforSet),'waters in the'
        for oxygenInfor in oxygenInforSet:
##            no1+=1
##            print no1
            frameNo=oxygenInfor[-2]
            frameNo=int(frameNo)-1
            segName=oxygenInfor[-3]
            resNumber=oxygenInfor[2]
            frame=rho.trajectory[frameNo]
            infor='segid '+segName+' and resid '+resNumber
            selected=rho.selectAtoms(infor)
            atomID=[]
            for atoms in selected.atoms:
                ID=str(atoms).split()[2][:-1]
                atomID.append(ID)
            selectedResId=selected.resids()
            selectedResNa=selected.resnames()
            coordsOH1H2=selected.coordinates()
            for i in range(3):
                atomInfor=str(selectedResNa[0])+'    '+str(atomID[i])+'    '+str(resNumber)+'    '+str(coordsOH1H2[i])[1:-1]+'   '+segName+'    '+str(frameNo)+'    '+str(no)+'\n'
                outputFile.write(atomInfor)
            H2OCoordinateSet.append(coordsOH1H2)

        no+=1

        H2OCoordinate.append(H2OCoordinateSet)
        print no,'is finished'
    outputFile.close()
    return H2OCoordinate

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def count_interactions(grof, xtcf, btime, cutoff, debug):
    u = Universe(grof, xtcf)
    query = ('(resname PRO and (name CB or name CG or name CD)) or'
             '(resname VAL and (name CG1 or name CG2)) or'
             '(resname GLY and name CA) or'
             '(resname ALA and name CB)')
    # MDAnalysis will convert the unit of length to angstrom, though in Gromacs the unit is nm
    atoms = u.selectAtoms(query)
    for ts in u.trajectory:
        if ts.time >= btime:
            numcount = 0
            for i, ai in enumerate(atoms):
                for j, aj in enumerate(atoms):
                    # to avoid counting the same pair twices,
                    # the 2 resid cannot be neigbors
                    if i < j and abs(ai.resid - aj.resid) >= 2: 
                        d = np.linalg.norm(ai.pos - aj.pos)
                        if d <= cutoff:
                            numcount += 1
            yield '{0:10.0f}{1:8d}\n'.format(ts.time,  numcount)
        # per 100 frames, num of frames changes with the size of xtc file, for debugging
        if debug and ts.frame % 2 == 0: 
            print "time: {0:10.0f}; step: {1:10d}; frame: {2:10d}".format(ts.time, ts.step, ts.frame)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def sequence_spacing(pf, grof, xtcf, peptide_length, atom_sel, output=None):
    u = Universe(grof, xtcf)
    # this selection part should be better customized
    # here, only have been backbone atoms are used, u.selectAtoms doesn't
    # include Hydrogen atoms
    # REMMEMBER: OPTIONS verification should be done in main ONLY!
    residues = [u.selectAtoms(atom_sel.format(i)) for i in range(2, peptide_length)]
    ijdist_dict = {}
    for ts in u.trajectory:
        for i, resi in enumerate(residues):
            for j, resj in enumerate(residues):
                if i < j:
                    resi_pos = resi.centerOfGeometry()                # residue i position
                    resj_pos = resj.centerOfGeometry()                # residue j position
                    ijdist = np.linalg.norm(resi_pos - resj_pos)
                    dij = j - i                                       # distance between i and j
                    if dij not in ijdist_dict.keys():
                        ijdist_dict[dij] = [dij]
                    else:
                        ijdist_dict[dij].append(ijdist)
        if ts.step % 2000000 == 0:                             # 2000ps
            print "time step: {0:d}".format(ts.step)
    return ijdist_dict

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
def main():
    arg_parser = argparse.ArgumentParser(description='通過給定殘基名稱，殘基內原子數目，原子在殘基內的索引(從0開始)，計算原子的坐標。')
    arg_parser.add_argument('resname', action='store', help='殘基名稱')
    arg_parser.add_argument('atoms_num', type=int, action='store', help='殘基內原子數目')
    arg_parser.add_argument('index', type=int, action='store', help='原子的索引，索引從0開始')
    arg_parser.add_argument('topology_file', action='store', help='拓撲文件，例如gro, pdb')
    args = arg_parser.parse_args()

    resname, atoms_num, index = args.resname, args.atoms_num, args.index

    universe = Universe(args.topology_file)
    atom_groups = universe.selectAtoms("resname " + resname)
    if len(atom_groups) % atoms_num != 0:
        print("拓撲文件內對應殘基原子總數不是所給原子數目的整數倍，請給予正確的原子數目。")
        exit(1)

    positions = []
    for i in range(0, len(atom_groups), atoms_num):
        positions.append(atom_groups[i:i + atoms_num][index].position)

    print("The positions of atoms %s is:" % (index))
    for i in positions:
        print(i)

# 需要導入模塊: from MDAnalysis import Universe [as 別名]
# 或者: from MDAnalysis.Universe import selectAtoms [as 別名]
import sys
sys.path.append('/home/x/xiansu/pfs/program/numpy/lib/python2.6/site-packages')
from MDAnalysis import Universe, Writer
from MDAnalysis.analysis.distances import distance_array
import MDAnalysis
import numpy
from Numeric import *


top='npt.gro'
traj='md_extract1.trr'


water=Universe(top,traj)

o=water.selectAtoms('name O*')

resid=o.resids()
print resid

#resnu=o.resnums()
#resna=o.resnames()

atomInf=[]
for i in o.atoms:
    atomid= str(i).split()[2]
    atomseg=str(i).split()[-1]
    atomidandseg=[]
    atomidandseg.append(atomid)
    atomidandseg.append(atomseg)
    atomInf.append(atomidandseg)