Python Utility.pairs方法代碼示例

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
 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

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
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))

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
 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)

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
 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

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
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]

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
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]

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
 def __init__(self, universe, rigid_bodies):
     """
     :param universe: the universe for which the subspace is created
     :type universe: :class:`~MMTK.Universe.Universe`
     :param rigid_bodies: a list or set of rigid bodies
                          with some common atoms
     """
     ex_ey_ez = [Vector(1.,0.,0.), Vector(0.,1.,0.), Vector(0.,0.,1.)]
     # Constructs
     # 1) a list of vectors describing the rigid-body motions of each
     #    rigid body as if it were independent.
     # 2) a list of pair-distance constraint vectors for all pairs of
     #    atoms inside a rigid body.
     # The LRB subspace is constructed from the projections of the
     # first set of vectors onto the orthogonal complement of the
     # subspace generated by the second set of vectors.
     vectors = []
     c_vectors = []
     for rb in rigid_bodies:
         atoms = rb.atomList()
         for d in ex_ey_ez:
             v = ParticleProperties.ParticleVector(universe)
             for a in atoms:
                 v[a] = d
             vectors.append(v)
         if len(atoms) > 1:
             center = rb.centerOfMass()
             iv = len(vectors)-3
             for d in ex_ey_ez:
                 v = ParticleProperties.ParticleVector(universe)
                 for a in atoms:
                     v[a] = d.cross(a.position()-center)
                 for vt in vectors[iv:]:
                     v -= v.dotProduct(vt)*vt
                 if v.dotProduct(v) > 0.:
                     vectors.append(v)
         for a1, a2 in Utility.pairs(atoms):
             distance = universe.distanceVector(a1.position(),
                                                a2.position())
             v = ParticleProperties.ParticleVector(universe)
             v[a1] = distance
             v[a2] = -distance
             c_vectors.append(v)
     if c_vectors:
         constraints = Subspace(universe, c_vectors)
         vectors = [constraints.projectionComplementOf(v)
                    for v in vectors]
     Subspace.__init__(self, universe, vectors)

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
def forceConstantTest(universe, atoms = None, delta = 0.0001):
    e0, grad0, fc = universe.energyGradientsAndForceConstants()
    if atoms is None:
	atoms = universe.atomList()
    for a1, a2 in map(lambda a: (a,a), atoms) + Utility.pairs(atoms):
	print a1, a2
	print fc[a1, a2]
	num_fc = []
	for v in [VectorModule.ex, VectorModule.ey, VectorModule.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 Numeric.array(map(lambda a: a.array, num_fc))

# 需要導入模塊: from MMTK import Utility [as 別名]
# 或者: from MMTK.Utility import pairs [as 別名]
 def excludedPairs(self, subset1, subset2, global_data):
     if 'excluded_pairs' not in global_data.get('initialized'):
         excluded_pairs = global_data.get('excluded_pairs')
         if subset1 is not None:
             for s1, s2 in [(subset1, subset2), (subset2, subset1)]:
                 set = {}
                 for a in s1.atomList():
                     set[a.index] = None
                 for a in s2.atomList():
                     try:
                         del set[a.index]
                     except KeyError: pass
                 excluded_pairs = excluded_pairs + \
                                  Utility.pairs(set.keys())
             set = {}
             for a in subset1.atomList():
                 set[a.index] = None
             for a in subset2.atomList():
                 set[a.index] = None
             atom_subset = set.keys()
             atom_subset.sort()
         else:
             atom_subset = None
         global_data.set('atom_subset', atom_subset)
         excluded_pairs = map(_normalizePair, excluded_pairs)
         excluded_pairs.sort(_cmpPair)
         _makeUnique(excluded_pairs)
         global_data.set('excluded_pairs', excluded_pairs)
         one_four_pairs = map(_normalizePair,
                              global_data.get('1_4_pairs'))
         one_four_pairs.sort(_cmpPair)
         _makeUnique(one_four_pairs)
         _eliminateExcluded(one_four_pairs, excluded_pairs)
         global_data.set('1_4_pairs', 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')