Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
conv_hatom = indices2.index(global_hatom)
typeHeavy = next((x.htype for x in heavyatoms if x.name == self.type_array[convindex1]),
AtomHBondType.none)
if typeHeavy == AtomHBondType.acc and (distarray[idx1, conv_hatom] <= hbondcutoff):
donorPosition = sel2.positions[idx2]
hydrogenPosition = sel2.positions[conv_hatom]
acceptorPosition = sel1.positions[idx1]
v1 = hydrogenPosition - acceptorPosition
v2 = hydrogenPosition - donorPosition
v1norm = np.linalg.norm(v1)
v2norm = np.linalg.norm(v2)
dot = np.dot(v1, v2)
angle = np.degrees(np.arccos(dot / (v1norm * v2norm)))
if angle >= hbondcutangle:
dist = distarray[idx1, conv_hatom]
new_hbond = HydrogenBond(convindex2, convindex1, global_hatom, dist, angle,
hbondcutoff,
hbondcutangle)
hydrogenBonds.append(new_hbond)
# print str(convindex1) + " " + str(convindex2)
# print "hbond found: %d,%d,%d"%(convindex2,global_hatom,convindex1)
# print angle
# finalize
newAtomContact = AtomContact(int(frame), float(distance), float(weight), int(convindex1),
int(convindex2),
hydrogenBonds)
currentFrameContacts.append(newAtomContact)
contactResults.append(currentFrameContacts)
stop = time.time()
# show trajectory information and selection information
print("trajectory with %d frames loaded" % len(u.trajectory))
typeHeavy = next((x.htype for x in heavyatoms if x.name == self.type_array[convindex2]),
AtomHBondType.none)
# print(typeHeavy)
if typeHeavy == AtomHBondType.acc and (distarray[conv_hatom, idx2] <= hbondcutoff):
donorPosition = sel1.positions[idx1]
hydrogenPosition = sel1.positions[conv_hatom]
acceptorPosition = sel2.positions[idx2]
v1 = hydrogenPosition - acceptorPosition
v2 = hydrogenPosition - donorPosition
v1norm = np.linalg.norm(v1)
v2norm = np.linalg.norm(v2)
dot = np.dot(v1, v2)
angle = np.degrees(np.arccos(dot / (v1norm * v2norm)))
if angle >= hbondcutangle:
dist = distarray[conv_hatom, idx2]
new_hbond = HydrogenBond(convindex1, convindex2, global_hatom, dist, angle,
hbondcutoff,
hbondcutangle)
hydrogenBonds.append(new_hbond)
# print(str(convindex1) + " " + str(convindex2)
# print("hbond found: %d,%d,%d"%(convindex1,global_hatom,convindex2))
# print(angle)
for global_hatom in hydrogenAtomsBoundToAtom2:
conv_hatom = indices2.index(global_hatom)
typeHeavy = next((x.htype for x in heavyatoms if x.name == self.type_array[convindex1]),
AtomHBondType.none)
if typeHeavy == AtomHBondType.acc and (distarray[idx1, conv_hatom] <= hbondcutoff):
donorPosition = sel2.positions[idx2]
hydrogenPosition = sel2.positions[conv_hatom]
acceptorPosition = sel1.positions[idx1]
v1 = hydrogenPosition - acceptorPosition
v2 = hydrogenPosition - donorPosition