Python FUNC.distance2方法代码示例

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
 def conect(self):
     # Return pairs of numbers that should be CONECTed
     # First extract the backbone IDs
     cg = self.cg()
     bb = [i+1 for i, j in zip(range(len(cg)), cg) if j[0] == "BB"]
     bb = zip(bb, bb[1:]+[len(bb)])
     # Set the backbone CONECTs (check whether the distance is consistent with binding)
     conect = [(i, j) for i, j in bb[:-1] if FUNC.distance2(cg[i-1][4:7], cg[j-1][4:7]) < 14]
     # Now add CONECTs for sidechains
     for i, j in bb:
         nsc = j-i-1

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
 def bbGetBond(self,r,ca,ss):
     # Retrieve parameters for each residue from tables defined above
     # Check is it DNA residue
     if r[0] in MAP.dnares3:
         return ca in self.dnaBbBondDictC.keys() and self.dnaBbBondDictC[ca] or None
     # RNA is not implemented properly yet
     elif r[0] in MAP.rnares3:
         return ca in self.rnaBbBondDictC.keys() and self.rnaBbBondDictC[ca] or None
     # If it's protein
     else:
         import FUNC 
         import math
         # The 150000 forceconstant gave an error message, turning to constraints would be better.
         return ( math.sqrt(FUNC.distance2(ca[0],ca[1]))/10., None )

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
def rubberBands(atomList, lowerBound, upperBound, decayFactor, decayPower, forceConstant, minimumForce):
    out = []
    u2  = upperBound**2
    while len(atomList) > 3:
        bi, xi = atomList.pop(0)
        for bj, xj in atomList[2:]:
            # Mind the nm/A conversion -- This has to be standardized! Global use of nm?
            d2 = FUNC.distance2(xi, xj)/100

            if d2 < u2:
                dij  = math.sqrt(d2)
                fscl = decayFunction(dij, lowerBound, decayFactor, decayPower)
                if fscl*forceConstant > minimumForce:
                    out.append({"atoms": (bi, bj), "parameters": (dij, "RUBBER_FC*%f" % fscl)})
    return out

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
def rubberBands(atomList,lowerBound,upperBound,decayFactor,decayPower,forceConstant,minimumForce):
    out = []
    u2  = upperBound**2
    while len(atomList) > 3:
        bi,xi = atomList.pop(0)
        # This is a bit weird (=wrong I think) way of doing the cutoff...
        #for bj,xj in atomList[2:]:
        for bj,xj in atomList:
            # Mind the nm/A conversion -- This has to be standardized! Global use of nm?
            d2 = FUNC.distance2(xi,xj)/100
            #if bi==73 and bj==79:
            #    print xi, xj, d2
            
            if d2 < u2:
                dij  = math.sqrt(d2)
                fscl = decayFunction(dij,lowerBound,decayFactor,decayPower)
                if fscl*forceConstant > minimumForce:
                    out.append({"atoms":(bi,bj),"parameters": (dij,"RUBBER_FC*%f"%fscl)})
    return out

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
 def bbGetBond(self,r,ca,ss):
     import FUNC 
     import math
     # The 150000 forceconstant gave an error message, turning to constraints would be better.
     return ( math.sqrt(FUNC.distance2(ca[0],ca[1]))/10., 150000 )

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
def check_merge(chains, m_list=[], l_list=[], ss_cutoff=0):
    chainIndex = range(len(chains))

    if 'all' in m_list:
        logging.info("All chains will be merged in a single moleculetype.")
        return chainIndex, [chainIndex]

    chainID = [chain.id for chain in chains]

    # Mark the combinations of chains that need to be merged
    merges = []
    if m_list:
        # Build a dictionary of chain IDs versus index
        # To give higher priority to top chains the lists are reversed
        # before building the dictionary
        chainIndex.reverse()
        chainID.reverse()
        dct = dict(zip(chainID, chainIndex))
        chainIndex.reverse()
        # Convert chains in the merge_list to numeric, if necessary
        # NOTE The internal numbering is zero-based, while the
        # command line chain indexing is one-based. We have to add
        # one to the number in the dictionary to bring it on par with
        # the numbering from the command line, but then from the
        # result we need to subtract one again to make indexing
        # zero-based
        merges = [[(i.isdigit() and int(i) or dct[i]+1)-1 for i in j] for j in m_list]
        for i in merges:
            i.sort()

    # Rearrange merge list to a list of pairs
    pairs = [(i[j], i[k]) for i in merges for j in range(len(i)-1) for k in range(j+1, len(i))]

    # Check each combination of chains for connections based on
    # ss-bridges, links and distance restraints
    for i in chainIndex[:-1]:
        for j in chainIndex[i+1:]:
            if (i, j) in pairs:
                continue
            # Check whether any link links these two groups
            for a, b in l_list:
                if ((a in chains[i] and b in chains[j]) or (a in chains[j] and b in chains[i])):
                    logging.info("Merging chains %d and %d to allow link %s" % (i+1, j+1, str((a, b))))
                    pairs.append(i < j and (i, j) or (j, i))
                    break
            if (i, j) in pairs:
                continue
            # Check whether any cystine bond given links these two groups
            #for a,b in s_list:
            #    if ((a in chains[i] and b in chains[j]) or
            #        (a in chains[j] and b in chains[i])):
            #        logging.info("Merging chains %d and %d to allow cystine bridge"%(i+1,j+1))
            #        pairs.append( i<j and (i,j) or (j,i) )
            #        break
            #if (i,j) in pairs:
            #    continue
            # Check for cystine bridges based on distance
            if not ss_cutoff:
                continue
            # Get SG atoms from cysteines from either chain
            # Check this pair of chains
            for cysA in chains[i]["CYS"]:
                for cysB in chains[j]["CYS"]:
                    d2 = FUNC.distance2(cysA["SG"][4:7], cysB["SG"][4:7])
                    if d2 <= ss_cutoff:
                        logging.info("Found SS contact linking chains %d and %d (%f nm)" % (i+1, j+1, math.sqrt(d2)/10))
                        pairs.append((i, j))
                    break
                if (i, j) in pairs:
                    break

    # Sort the combinations
    pairs.sort(reverse=True)

    merges = []
    while pairs:
        merges.append(set([pairs[-1][0]]))
        for i in range(len(pairs)-1, -1, -1):
            if pairs[i][0] in merges[-1]:
                merges[-1].add(pairs.pop(i)[1])
            elif pairs[i][1] in merges[-1]:
                merges[-1].add(pairs.pop(i)[0])
    merges = [list(i) for i in merges]
    for i in merges:
        i.sort()

    order = [j for i in merges for j in i]

    if merges:
        logging.warning("Merging chains.")
        logging.warning("This may change the order of atoms and will change the number of topology files.")
        logging.info("Merges: " + ", ".join([str([j+1 for j in i]) for i in merges]))

    if len(merges) == 1 and len(merges[0]) > 1 and set(merges[0]) == set(chainIndex):
        logging.info("All chains will be merged in a single moleculetype")

    # Determine the order for writing; merged chains go first
    merges.extend([[j] for j in chainIndex if j not in order])
    order.extend([j for j in chainIndex if j not in order])

#.........这里部分代码省略.........

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
def contacts(atoms, cutoff=5):
    rla = range(len(atoms))
    crd = [atom[4:] for atom in atoms]
    return [(i, j) for i in rla[:-1] for j in rla[i+1:]
            if FUNC.distance2(crd[i], crd[j]) < cutoff]

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]
def residueDistance2(r1, r2):
    return min([FUNC.distance2(i, j) for i in r1 for j in r2])

# 需要导入模块: import FUNC [as 别名]
# 或者: from FUNC import distance2 [as 别名]

#.........这里部分代码省略.........
            chain.set_ss(ssAver[:len(chain)])
            ssAver = ssAver[len(chain):]

        # Now the chains are complete, each consisting of a residuelist,
        # and a secondary structure designation if the chain is of type 'Protein'.
        # There may be mixed chains, there may be HETATM things.
        # Water has been discarded. Maybe this has to be changed at some point.
        # The order in the coarse grained files matches the order in the set of chains.
        #
        # If there are no merges to be done, i.e. no global Elnedyn network, no
        # disulphide bridges, no links, no distance restraints and no explicit merges,
        # then we can write out the topology, which will match the coarse grained file.
        #
        # If there are merges to be done, the order of things may be changed, in which
        # case the coarse grained structure will not match with the topology...

        # CYSTINE BRIDGES #
        # Extract the cysteine coordinates (for all frames) and the cysteine identifiers
        if options['CystineCheckBonds']:
            logging.info("Checking for cystine bridges, based on sulphur (SG) atoms lying closer than %.4f nm" % math.sqrt(options['CystineMaxDist2']/100))

            cyscoord  = zip(*[[j[4:7] for j in i] for i in cysteines])
            cysteines = [i[:4] for i in cysteines[0]]

            bl, kb    = options['ForceField'].special[(("SC1", "CYS"), ("SC1", "CYS"))]

            # Check the distances and add the cysteines to the link list if the
            # SG atoms have a distance smaller than the cutoff.
            rlc = range(len(cysteines))
            for i in rlc[:-1]:
                for j in rlc[i+1:]:
                    # Checking the minimum distance over all frames
                    # But we could also take the maximum, or the mean
                    d2 = min([FUNC.distance2(a, b) for a, b in zip(cyscoord[i], cyscoord[j])])
                    if d2 <= options['CystineMaxDist2']:
                        a, b = cysteines[i], cysteines[j]
                        options['linkListCG'].append((("SC1", "CYS", a[2], a[3]), ("SC1", "CYS", b[2], b[3]), bl, kb))
                        a, b = (a[0], a[1], a[2]-(32 << 20), a[3]), (b[0], b[1], b[2]-(32 << 20), b[3])
                        logging.info("Detected SS bridge between %s and %s (%f nm)" % (a, b, math.sqrt(d2)/10))

        # REAL ITP STUFF #
        # Check whether we have identical chains, in which case we
        # only write the ITP for one...
        # This means making a distinction between chains and
        # moleculetypes.

        molecules = [tuple([chains[i] for i in j]) for j in merge]

        # At this point we should have a list or dictionary of chains
        # Each chain should be given a unique name, based on the value
        # of options["-o"] combined with the chain identifier and possibly
        # a number if there are chains with identical identifiers.
        # For each chain we then write an ITP file using the name for
        # moleculetype and name + ".itp" for the topology include file.
        # In addition we write a master topology file, using the value of
        # options["-o"], with an added extension ".top" if not given.

        # XXX *NOTE*: This should probably be gathered in a 'Universe' class
        itp = 0
        moleculeTypes = {}
        for mi in range(len(molecules)):
            mol = molecules[mi]
            # Check if the moleculetype is already listed
            # If not, generate the topology from the chain definition
            if mol not in moleculeTypes or options['SeparateTop']:
                # Name of the moleculetype