Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
for residue in self.residues:
if not isinstance(residue, aa.Amino):
continue
residue.dihedrals = []
refangles = residue.reference.dihedrals
for dihed in refangles:
coords = []
atoms = dihed.split()
for i in range(4):
atomname = atoms[i]
if residue.has_atom(atomname):
coords.append(residue.get_atom(atomname).coords)
if len(coords) == 4:
angle = util.dihedral(coords[0], coords[1], coords[2], coords[3])
else:
angle = None
residue.add_dihedral_angle(angle)
if self.has_atom(bond) and bond.startswith("H"):
hatoms.append(self.get_atom(bond))
# If this is more than two something is wrong
if len(hatoms) != 2:
return False
# Use the existing hydrogen and rotate about the bond
self.rotate_tetrahedral(nextatom, bondatom, 120)
newcoords1 = hatoms[0].coords
self.rotate_tetrahedral(nextatom, bondatom, 120)
newcoords2 = hatoms[0].coords
self.rotate_tetrahedral(nextatom, bondatom, 120)
# Determine which one hatoms[1] is not in
if util.distance(hatoms[1].coords, newcoords1) > 0.1:
self.create_atom(atomname, newcoords1)
else:
self.create_atom(atomname, newcoords2)
return True
return False
for i in range(len(chain.residues) - 1):
res1 = chain.residues[i]
res2 = chain.residues[i + 1]
if not isinstance(res1, aa.Amino) or not isinstance(res2, aa.Amino):
continue
atom1 = res1.get_atom("C")
atom2 = res2.get_atom("N")
if atom1 is not None:
res2.peptide_c = atom1
if atom2 is not None:
res1.peptide_n = atom2
if atom1 is None or atom2 is None:
continue
if util.distance(atom1.coords, atom2.coords) > PEPTIDE_DIST:
text = "Gap in backbone detected between %s and %s!" % \
(res1, res2)
_LOGGER.warning(text)
res2.peptide_c = None
res1.peptide_n = None
options: The name of the forcefield (string)
Returns
header: The PQR file header (string)
lines: The PQR file atoms (list)
missedligandresidues: A list of ligand residue names whose charges could
not be assigned (ligand)
protein: The protein object
"""
pkaname = ""
lines = []
Lig = None
atomcount = 0 # Count the number of ATOM records in pdb
outroot = utilities.getPQRBaseFileName(options.output_pqr)
if options.pka_method == 'propka':
pkaname = outroot + ".propka"
#TODO: What? Shouldn't it be up to propka on how to handle this?
if os.path.isfile(pkaname):
os.remove(pkaname)
start = time.time()
_LOGGER.info("Beginning PDB2PQR...")
myDefinition = definitions.Definition()
_LOGGER.info("Parsed Amino Acid definition file.")
if options.drop_water:
# Remove the waters
pdblist_new = []
def rotate_tetrahedral(cls, atom1, atom2, angle):
"""Rotate about the atom1-atom2 bond by a given angle All atoms connected
to atom2 will rotate.
Args:
atom1: The first atom of the bond to rotate about (atom)
atom2: The second atom of the bond to rotate about (atom)
angle: The number of degrees to rotate (float)
"""
moveatoms = []
movecoords = []
initcoords = util.subtract(atom2.coords, atom1.coords)
# Determine which atoms to rotate
for atom in atom2.bonds:
if atom == atom1:
continue
moveatoms.append(atom)
movecoords.append(util.subtract(atom.coords, atom1.coords))
newcoords = quat.qchichange(initcoords, movecoords, angle)
for iatom, atom in enumerate(moveatoms):
x = newcoords[iatom][0] + atom1.x
y = newcoords[iatom][1] + atom1.y
z = newcoords[iatom][2] + atom1.z
atom.x = x
atom.y = y
atom.z = z
y = (newcoords[iatom][1] + coordlist[1][1])
z = (newcoords[iatom][2] + coordlist[1][2])
atom.x = x
atom.y = y
atom.z = z
self.cells.add_cell(atom)
# Set the new angle
coordlist = []
for atomname in atomnames:
if residue.has_atom(atomname):
coordlist.append(residue.get_atom(atomname).coords)
else:
raise ValueError("Error occurred while trying to debump!")
dihed = util.dihedral(coordlist[0], coordlist[1], coordlist[2], coordlist[3])
residue.dihedrals[anglenum] = dihed
closeatoms = self.debumper.cells.get_near_cells(atom)
for closeatom in closeatoms:
# Conditions for continuing
if atom.residue == closeatom.residue:
continue
if not (closeatom.hacceptor or closeatom.hdonor):
continue
if atom.hdonor and not atom.hacceptor:
if not closeatom.hacceptor:
continue
if atom.hacceptor:
if not atom.hdonor and not closeatom.hdonor:
continue
dist = util.distance(atom.coords, closeatom.coords)
if dist < 4.3:
residue = atom.residue
hbond = structures.PotentialBond(atom, closeatom, dist)
# Store the potential bond
obj.hbonds.append(hbond)
# Keep track of connectivity
if closeatom in self.atomlist:
closeobj = self.resmap[closeatom.residue]
if closeobj not in connectivity[obj]:
connectivity[obj].append(closeobj)
progress += increment
while progress >= 0.0499:
progress -= 0.05
if not isinstance(closeresidue, (aa.Amino, aa.WAT)):
continue
if isinstance(residue, aa.CYS) and residue.ss_bonded_partner == closeatom:
continue
# Also ignore if this is a donor/acceptor pair
if (atom.is_hydrogen and len(atom.bonds) != 0 and \
atom.bonds[0].hdonor and closeatom.hacceptor):
continue
if (closeatom.is_hydrogen and len(closeatom.bonds) != 0 and \
closeatom.bonds[0].hdonor and atom.hacceptor):
continue
dist = util.distance(atom.coords, closeatom.coords)
other_size = BUMP_HYDROGEN_SIZE if closeatom.is_hydrogen else BUMP_HEAVY_SIZE
cutoff = atom_size + other_size
if dist < cutoff:
nearatoms[closeatom] = cutoff - dist
return nearatoms
return 0
elif residue.has_atom("LP1"):
newname = "LP2"
else:
newname = "LP1"
_LOGGER.debug("Working on %s %s (acceptor) to %s %s (donor)",
acc.residue, acc.name, donor.residue, donor.name)
# Act depending on the number of bonds
if len(acc.bonds) == 0:
if self.is_hbond(donor, acc):
# Find the best donor hydrogen and use that
bestdist = util.distance(acc.coords, donor.coords)
for donorh in donor.bonds:
dist = util.distance(acc.coords, donorh.coords)
if dist < bestdist:
bestdist = dist
# Point the LP to the best H
# TODO - this looks like a bug...
# donorh ends up being the last item in donor.bonds
# This may be fixed by setting a best_donorh to go with bestdist
# and using best_donorh in the function below
self.make_atom_with_no_bonds(acc, donorh, newname)
_LOGGER.warning("The best donorH was not picked (BUG?).")
_LOGGER.debug("Added %s to %s", newname, acc.residue)
return 1
return 0
if not residue.has_atom("H1"):
addname = "H1"
else:
addname = "H2"
_LOGGER.debug("Finalizing %s by adding %s (%i current O bonds)",
residue, addname, len(atom.bonds))
if len(atom.bonds) == 0:
newcoords = []
# Build hydrogen away from closest atom
closeatom = self.routines.get_closest_atom(atom)
if closeatom != None:
vec = util.subtract(atom.coords, closeatom.coords)
dist = util.distance(atom.coords, closeatom.coords)
for i in range(3):
newcoords.append(vec[i]/dist + atom.coords[i])
else:
newcoords = util.add(atom.coords, [1.0, 0.0, 0.0])
residue.create_atom(addname, newcoords)
self.routines.cells.add_cell(residue.get_atom(addname))
self.finalize()
elif len(atom.bonds) == 1:
# Initialize variables