1"""
2    Contact extension
3
4    Find all hydrogen bonds as determined by the DISTANCE cutoff below.
5    Uses PDB2PQR to determine donors and acceptors, and displays
6    all available bonds to stdout in a WHATIF-like format.
7
8    Author:  Julie C. Mitchell
9"""
10
11__date__ = "April 2007"
12__author__ = "Julie C. Mitchell"
13
14from src.utilities import distance
15from src.routines import Cells
16
17DIST_CUTOFF = 3.5         # max distance
18
19def usage():
20    return 'Print a list of contacts to {output-path}.con\n'
21
22def run_extension(routines, outroot, options):
23    """
24        Print a list of contacts.
25
26        Parameters
27            routines:  A link to the routines object
28            outroot:   The root of the output name
29            options:   options object
30    """
31    outname = outroot + ".con"
32    outfile = open(outname, "w")
33
34    # Initialize - set nearby cells, donors/acceptors
35
36    cellsize = int(DIST_CUTOFF + 1.0 + 1.0)
37    protein = routines.protein
38    routines.setDonorsAndAcceptors()
39    routines.cells = Cells(cellsize)
40    routines.cells.assignCells(protein)
41
42    for thisatom in protein.getAtoms():
43        # Grab the list of thisatoms
44        if not thisatom.hdonor:
45            continue
46        thisatomhs = []
47        for bond in thisatom.bonds:
48            if bond.isHydrogen():
49                thisatomhs.append(bond)
50        if thisatomhs == []:
51            continue
52
53        # For each thisatom, grab all thatatomeptors
54
55        count = 0
56        closeatoms = routines.cells.getNearCells(thisatom)
57        for thatatom in closeatoms:
58            if (thisatom.residue == thatatom.residue):
59                continue  # comment this out to include interresidue contacts
60            if (thatatom.isHydrogen()):
61                continue
62            thisdist = distance(thisatom.getCoords(), thatatom.getCoords())
63            if (thisdist <= DIST_CUTOFF):
64                count = count+1
65                thisBstring='S'
66                thatBstring='S'
67                hscore= 0.0
68                if (thisatom.hdonor & thatatom.hacceptor):
69                    hscore = 1.0
70                if (thisatom.hacceptor & thatatom.hdonor):
71                    hscore = 1.0
72                if (thisatom.isBackbone()):
73                    thisBstring='B'
74                if (thatatom.isBackbone()):
75                    thatBstring='B'
76                outfile.write("%4d %4d %-4s (%4d  ) %s     %-4s<>%4d %-4s (%4d  ) %s     %-4s D=%6.2f  H-ene=%6.2f  Sym=  (%s-%s)\n" % \
77                  (count, thisatom.residue.resSeq,thisatom.residue.name,thisatom.residue.resSeq, thisatom.residue.chainID,thisatom.name,thatatom.residue.resSeq,thatatom.residue.name,thatatom.residue.resSeq, thatatom.residue.chainID,thatatom.name, thisdist, hscore, thisBstring, thatBstring))
78
79    routines.write("\n")
80    outfile.close()
81