Python MMTK.Utility类代码示例

 def __init__(self, atoms, constraints):
     self.atoms = atoms
     natoms = len(self.atoms)
     nconst = sum([len(c) for c in constraints])
     b = N.zeros((nconst, natoms), N.Float)
     c = N.zeros((nconst,), N.Float)
     i = 0
     for cons in constraints:
         cons.setCoefficients(self.atoms, b, c, i)
         i = i + len(cons)
     u, s, vt = LA.singular_value_decomposition(b)
     self.rank = 0
     for i in range(min(natoms, nconst)):
         if s[i] > 0.:
             self.rank = self.rank + 1
     self.b = b
     self.bi = LA.generalized_inverse(b)
     self.p = N.identity(natoms)-N.dot(self.bi, self.b)
     self.c = c
     self.bi_c = N.dot(self.bi, c)
     c_test = N.dot(self.b, self.bi_c)
     if N.add.reduce((c_test-c)**2)/nconst > 1.e-12:
         Utility.warning("The charge constraints are inconsistent."
                         " They will be applied as a least-squares"
                         " condition.")

 def __init__(self, b1, b2, ca):
     self.b1 = b1 # bond 1
     self.b2 = b2 # bond 2
     self.ca = ca # common atom
     if Utility.uniqueID(self.b2) < Utility.uniqueID(self.b1):
         self.b1, self.b2 = self.b2, self.b1
     self.a1 = b1.otherAtom(ca)
     self.a2 = b2.otherAtom(ca)
     Utility.uniqueID.registerObject(self)

    def write(self, object, configuration = None, tag = None):
        """
        Write  an object to the file

        :param object: the object to be written
        :type object: :class:`~MMTK.Collections.GroupOfAtoms`
        :param configuration: the configuration from which the coordinates 
                              are taken (default: current configuration)
        :type configuration: :class:`~MMTK.ParticleProperties.Configuration`
        """
        if not ChemicalObjects.isChemicalObject(object):
            for o in object:
                self.write(o, configuration)
        else:
            toplevel = tag is None
            if toplevel:
                tag = Utility.uniqueAttribute()
            if hasattr(object, 'pdbmap'):
                for residue in object.pdbmap:
                    self.file.nextResidue(residue[0], )
                    sorted_atoms = residue[1].items()
                    sorted_atoms.sort(lambda x, y:
                                      cmp(x[1].number, y[1].number))
                    for atom_name, atom in sorted_atoms:
                        atom = object.getAtom(atom)
                        p = atom.position(configuration)
                        if Utility.isDefinedPosition(p):
                            try: occ = atom.occupancy
                            except AttributeError: occ = 0.
                            try: temp = atom.temperature_factor
                            except AttributeError: temp = 0.
                            self.file.writeAtom(atom_name, p/Units.Ang,
                                                occ, temp, atom.type.symbol)
                            self.atom_sequence.append(atom)
                        else:
                            self.warning = True
                        setattr(atom, tag, None)
            else:
                if hasattr(object, 'is_protein'):
                    for chain in object:                    
                        self.write(chain, configuration, tag)
                elif hasattr(object, 'is_chain'):
                        self.file.nextChain(None, object.name)
                        for residue in object:
                            self.write(residue, configuration, tag)
                        self.file.terminateChain()
                elif hasattr(object, 'molecules'):
                    for m in object.molecules:
                        self.write(m, configuration, tag)
                elif hasattr(object, 'groups'):
                    for g in object.groups:
                        self.write(g, configuration, tag)
            if toplevel:
                for a in object.atomList():
                    if not hasattr(a, tag):
                        self.write(a, configuration, tag)
                    delattr(a, tag)

 def close(self):
     """
     Closes the file. Must be called in order to prevent data loss.
     """
     if self.model_number is not None:
         self.file.writeLine('ENDMDL', '')
     self.file.close()
     if self.warning:
         Utility.warning('Some atoms are missing in the output file ' + \
                         'because their positions are undefined.')
         self.warning = False

def setConfiguration(object, pdb_residues,
                     map = 'pdbmap', alt = 'pdb_alternative',
                     atom_map = None, toplevel = True):
    defined = 0
    if hasattr(object, 'is_protein'):
        i = 0
        for chain in object:
            l = len(chain)
            defined += setConfiguration(chain, pdb_residues[i:i+l],
                                        map, alt, atom_map, False)
            i = i + l
    elif hasattr(object, 'is_chain'):
        for i in range(len(object)):
            defined += setConfiguration(object[i], pdb_residues[i:i+1],
                                        map, alt, atom_map, False)
    elif hasattr(object, map):
        pdbmap = getattr(object, map)
        try: altmap = getattr(object, alt)
        except AttributeError: altmap = {}
        nres = len(pdb_residues)
        if len(pdbmap) != nres:
            raise IOError('PDB configuration does not match object ' +
                           object.fullName())
        for i in range(nres):
            defined += setResidueConfiguration(object, pdb_residues[i],
                                               pdbmap[i], altmap, atom_map)
    elif Collections.isCollection(object):
        nres = len(pdb_residues)
        if len(object) != nres:
            raise IOError('PDB configuration does not match object ' +
                           object.fullName())
        for i in range(nres):
            defined += setConfiguration(object[i], [pdb_residues[i]],
                                        map, alt, atom_map, False)
    else:
        try:
            name = object.fullName()
        except AttributeError:
            try:
                name = object.name
            except AttributeError:
                name = '???'
        raise IOError('PDB configuration does not match object ' + name)
              
    if toplevel and defined < object.numberOfAtoms():
        name = '[unnamed object]'
        try:
            name = object.fullName()
        except: pass
        if name: name = ' in ' + name
        Utility.warning(`object.numberOfAtoms()-defined` + ' atom(s)' + name +
                        ' were not assigned (new) positions.')
    return defined

 def pairsWithinCutoff(self, cutoff):
     """
     :param cutoff: a cutoff for pair distances
     :returns: a list containing all pairs of objects in the
               collection whose center-of-mass distance is less than
               the cutoff
     :rtype: list
     """
     pairs = []
     positions = {}
     for index, objects in self.partition.items():
         pos = map(lambda o: o.position(), objects)
         positions[index] = pos
         for o1, o2 in Utility.pairs(zip(objects, pos)):
             if (o2[1]-o1[1]).length() <= cutoff:
                 pairs.append((o1[0], o2[0]))
     partition_cutoff = int(N.floor((cutoff/self.partition_size)**2))
     ones = N.array([1,1,1])
     zeros = N.array([0,0,0])
     keys = self.partition.keys()
     for i in range(len(keys)):
         p1 = keys[i]
         for j in range(i+1, len(keys)):
             p2 = keys[j]
             d = N.maximum(abs(N.array(p2)-N.array(p1)) -
                           ones, zeros)
             if N.add.reduce(d*d) <= partition_cutoff:
                 for o1, pos1 in zip(self.partition[p1],
                                     positions[p1]):
                     for o2, pos2 in zip(self.partition[p2],
                                         positions[p2]):
                         if (pos2-pos1).length() <= cutoff:
                             pairs.append((o1, o2))
     return pairs

def forceConstantTest(universe, atoms = None, delta = 0.0001):
    """
    Test force constants by comparing to the numerical derivatives
    of the gradients.

    :param universe: the universe on which the test is performed
    :type universe: :class:`~MMTK.Universe.Universe`
    :param atoms: the atoms of the universe for which the gradient
                  is tested (default: all atoms)
    :type atoms: list
    :param delta: the step size used in calculating the numerical derivatives
    :type delta: float
    """
    e0, grad0, fc = universe.energyGradientsAndForceConstants()
    if atoms is None:
        atoms = universe.atomList()
    for a1, a2 in itertools.chain(itertools.izip(atoms, atoms),
                                  Utility.pairs(atoms)):
        print a1, a2
        print fc[a1, a2]
        num_fc = []
        for v in [ex, ey, ez]:
            x = a1.position()
            a1.setPosition(x+delta*v)
            e_plus, grad_plus = universe.energyAndGradients()
            a1.setPosition(x-delta*v)
            e_minus, grad_minus = universe.energyAndGradients()
            a1.setPosition(x)
            num_fc.append(0.5*(grad_plus[a2]-grad_minus[a2])/delta)
        print N.array(map(lambda a: a.array, num_fc))

 def dihedralAngles(self):
     """
     :returns: a list of all dihedral angles that can be formed from the
               bond angles in the list
     :rtype: :class:`~MMTK.Bonds.DihedralAngleList`
     """
     # find all bonds that are involved in more than one bond angle
     angles = {}
     bond_list = []
     for angle in self.data:
         try:
             al = angles[angle.b1]
         except KeyError:
             al = []
             angles[angle.b1] = al
             bond_list.append(angle.b1)
         al.append(angle)
         try:
             al = angles[angle.b2]
         except KeyError:
             al = []
             angles[angle.b2] = al
             bond_list.append(angle.b2)
         al.append(angle)
     dihedrals = []
     for bond in bond_list:
         # each pair of bond angles with a common bond defines a dihedral
         for p in Utility.pairs(angles[bond]):
             d = DihedralAngle(p[0], p[1], bond)
             if d.normalized:
                 dihedrals.append(d)
     return DihedralAngleList(dihedrals)

 def bondAngles(self):
     """
     :returns: a list of all bond angles that can be formed from the
               bonds in the list
     :rtype: :class:`~MMTK.Bonds.BondAngleList`
     """
     if self.bond_angles is None:
         # find all atoms that are involved in more than one bond
         bonds = {}
         atom_list = []
         for bond in self:
             try:
                 bl = bonds[bond.a1]
             except KeyError:
                 bl = []
                 bonds[bond.a1] = bl
                 atom_list.append(bond.a1)
             bl.append(bond)
             try:
                 bl = bonds[bond.a2]
             except KeyError:
                 bl = []
                 bonds[bond.a2] = bl
                 atom_list.append(bond.a2)
             bl.append(bond)
         angles = []
         for atom in atom_list:
             # each pair of bonds at the same atom defines a bond angle
             for p in Utility.pairs(bonds[atom]):
                 angles.append(BondAngle(p[0], p[1], atom))
         self.bond_angles = BondAngleList(angles)
     return self.bond_angles

 def setReferences(self):
     atom_refs = []
     for i in range(len(self.atoms)):
         atom_refs.append(AtomReference(i))
     for attr in vars(self).items():
         if attr[0] not in self.instance:
             setattr(self, attr[0],
                     Utility.substitute(getattr(self, attr[0]),
                                        self.atoms, atom_refs))

 def atomsWithDefinedPositions(self, conf = None):
     """
     :param conf: a configuration object, or None for the
                  current configuration
     :type conf: :class:`~MMTK.ParticleProperties.Configuration` or NoneType
     :returns: a collection of all atoms that have a position in the
               given configuration
     """
     return Collection([a for a in self.atomIterator()
                        if Utility.isDefinedPosition(a.position(conf))])

def addToForceConstants(total_fc, indices, small_fc):
    indices = zip(indices, range(len(indices)))
    for i1, i2 in itertools.chain(itertools.izip(indices, indices),
                                  Utility.pairs(indices)):
        ii1, jj1 = i1
        ii2, jj2 = i2
        if ii1 > ii2:
            ii1, ii2 = ii2, ii1
            jj1, jj2 = jj2, jj1
        jj1 = 3*jj1
        jj2 = 3*jj2
        total_fc[ii1,:,ii2,:] += small_fc[jj1:jj1+3, jj2:jj2+3]

def addToForceConstants(total_fc, indices, small_fc):
    indices = map(None, indices, range(len(indices)))
    for i1, i2 in map(lambda i: (i,i), indices) + Utility.pairs(indices):
	ii1, jj1 = i1
	ii2, jj2 = i2
	if ii1 > ii2:
	    ii1, ii2 = ii2, ii1
	    jj1, jj2 = jj2, jj1
	jj1 = 3*jj1
	jj2 = 3*jj2
	total_fc[ii1,:,ii2,:] = total_fc[ii1,:,ii2,:] + \
				small_fc[jj1:jj1+3, jj2:jj2+3]

 def excludedPairs(self, subset1, subset2, global_data):
     if 'excluded_pairs' not in global_data.get('initialized'):
         excluded_pairs = set(global_data.get('excluded_pairs'))
         if subset1 is not None:
             set1 = set(a.index for a in subset1.atomList())
             set2 = set(a.index for a in subset2.atomList())
             excluded_pairs |= set(Utility.orderedPairs(list(set1-set2)))
             excluded_pairs |= set(Utility.orderedPairs(list(set2-set1)))
             atom_subset = list(set1 | set2)
             atom_subset.sort()
         else:
             atom_subset = None
         global_data.set('atom_subset', atom_subset)
         global_data.set('excluded_pairs', list(excluded_pairs))
         one_four_pairs = set(global_data.get('1_4_pairs')) \
                            - excluded_pairs
         global_data.set('1_4_pairs', list(one_four_pairs))
         global_data.add('initialized', 'excluded_pairs')
     return global_data.get('excluded_pairs'), \
            global_data.get('1_4_pairs'), \
            global_data.get('atom_subset')

def databasePath(filename, directory, try_direct = False):
    if Utility.isURL(filename):
        return filename
    filename = os.path.expanduser(filename)
    if try_direct and os.path.exists(filename):
        return os.path.normcase(filename)
    entries = []
    if os.path.split(filename)[0] == '':
        for p in path:
            if Utility.isURL(p):
                url = Utility.joinURL(p, directory+'/'+filename)
                if Utility.checkURL(url):
                    entries.append(url)
            else:
                full_name = os.path.join(os.path.join(p, directory), filename)
                if os.path.exists(full_name):
                    entries.append(os.path.normcase(full_name))
    if len(entries) == 0:
        raise IOError("Database entry %s/%s not found" % (directory, filename))
    else:
        if len(entries) > 1:
            Utility.warning("multiple database entries for %s/%s, using first one"
                            % (directory, filename))
            for e in entries:
                sys.stderr.write(e+'\n')
        return entries[0]