1"""
2    Saltbridge extension
3
4    Find all salt bridges as determined by the cutoff distance below.
5    Uses PDB2PQR to determine atom identities and distances, and write
6    out all located salt bridges to stdout.
7
8    NOTE:  A bond may be labeled BOTH hbond and salt-bridge if you use both
9           options in one pdb2pqr run.  Look out for double counting.
10
11    NOTE:  This extension currently does not support salt bridges with chain termini.
12
13    Author:  Mike Bradley (heavily copied from Todd Dolinsky's hbond extension)
14"""
15
16__date__ = "25 August 2006"
17__author__ = "Mike Bradley"
18
19from src.utilities import distance
20from src.routines import Cells
21
22DIST_CUTOFF = 4.0         # maximum cation to anion atom distance in angstroms
23
24def usage():
25    return 'Print a list of salt bridges to {output-path}.salt'
26
27def run_extension(routines, outroot, options):
28    """
29        Print a list of salt bridges.
30
31        Parameters
32            routines:  A link to the routines object
33            outroot:   The root of the output name
34            options:   options object
35    """
36    outname = outroot + ".salt"
37    outfile = open(outname, "w")
38
39    routines.write("Printing salt bridge list...\n")
40
41    # Define the potential salt bridge atoms here (the current lists are for the AMBER
42    # forcefield and are not necessarily exhaustive).
43    posresList = ["LYS","ARG","HIP"]
44    negresList = ["GLU","ASP","CYM"]
45    posatomList = ["NE","NH1","NH2","NZ","ND1","NE2",]
46    negatomList = ["SG","OE1","OE2","OD1","OD2"]
47
48    # Initialize - set nearby cells
49    # The cell size adds one for the salt bridge distance, and rounds up
50
51    cellsize = int(DIST_CUTOFF + 1.0 + 1.0)
52    protein = routines.protein
53    routines.cells = Cells(cellsize)
54    routines.cells.assignCells(protein)
55
56    # Loop over all the atoms
57    for cation in protein.getAtoms():
58        # check that we've found a cation
59        if cation.residue.name == "NMET":
60            print "YES NMET"
61        if cation.residue.name not in posresList:
62            continue
63        elif cation.name not in posatomList:
64                continue
65        # For each cation, grab all potential anions in nearby cells
66        closeatoms = routines.cells.getNearCells(cation)
67        for anion in closeatoms:
68            if cation.residue.name == anion.residue.name:
69                continue
70            if anion.residue.name not in negresList:
71                continue
72            elif anion.name not in negatomList:
73                    continue
74            # Do distance check
75            dist = distance(cation.getCoords(), anion.getCoords())
76            if dist > DIST_CUTOFF:
77                continue
78            #routines.write("Cation: %s %s\tAnion: %s %s\tsaltdist: %.2f\n" % \
79            #          (cation.residue, cation.name, anion.residue, anion.name, dist))
80            outfile.write("Cation: %s %s\tAnion: %s %s\tsaltdist: %.2f\n" % \
81                      (cation.residue, cation.name, anion.residue, anion.name, dist))
82    #routines.write("\n")
83    outfile.close()
84