Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def _run_amber_traj(traj, ext_ref):
# Test triclinic case where simple approach in Tuckerman text does not
# always work
distopt = md.compute_distances(traj, [[0, 9999]], opt=True)
distslw = md.compute_distances(traj, [[0, 9999]], opt=False)
dispopt = md.compute_displacements(traj, [[0, 9999]], opt=True)
dispslw = md.compute_displacements(traj, [[0, 9999]], opt=False)
eq(distopt, distslw, decimal=5)
eq(dispopt, dispslw, decimal=5)
assert_allclose(distopt.flatten(), ext_ref, atol=2e-5)
# Make sure distances from displacements are the same
eq(np.sqrt((dispopt.squeeze() ** 2).sum(axis=1)), distopt.squeeze())
eq(np.sqrt((dispslw.squeeze() ** 2).sum(axis=1)), distslw.squeeze())
eq(dispopt, dispslw, decimal=5)
def _run_amber_traj(traj, ext_ref):
# Test triclinic case where simple approach in Tuckerman text does not
# always work
distopt = md.compute_distances(traj, [[0, 9999]], opt=True)
distslw = md.compute_distances(traj, [[0, 9999]], opt=False)
dispopt = md.compute_displacements(traj, [[0, 9999]], opt=True)
dispslw = md.compute_displacements(traj, [[0, 9999]], opt=False)
eq(distopt, distslw, decimal=5)
eq(dispopt, dispslw, decimal=5)
assert_allclose(distopt.flatten(), ext_ref, atol=2e-5)
# Make sure distances from displacements are the same
eq(np.sqrt((dispopt.squeeze() ** 2).sum(axis=1)), distopt.squeeze())
eq(np.sqrt((dispslw.squeeze() ** 2).sum(axis=1)), distslw.squeeze())
eq(dispopt, dispslw, decimal=5)
"""
cell_lengths, cell_angles = box_vectors_to_lengths_angles(walker.state['box_vectors'])
t2 = time.time()
# make a traj out of it so we can calculate distances through
# the periodic boundary conditions
walker_traj = mdj.Trajectory(walker.state['positions'],
topology=self._mdj_top,
unitcell_lengths=cell_lengths,
unitcell_angles=cell_angles)
t3 = time.time()
# calculate the distances through periodic boundary conditions
# and get hte minimum distance
min_distance = np.min(mdj.compute_distances(walker_traj,
it.product(self.ligand_idxs,
self.receptor_idxs),
periodic=self._periodic)
)
t4 = time.time()
logging.info("Make a traj: {0}; Calc dists: {1}".format(t3-t2,t4-t3))
return min_distance
self.section("Bond Check (without PBCs)")
radii = []
pairs = []
for (a, b) in self.topology.bonds:
try:
radsum = COVALENT_RADII[a.element.symbol] + COVALENT_RADII[b.element.symbol]
except KeyError:
raise NotImplementedError("I don't have radii information for all of your atoms")
radii.append(radsum)
pairs.append((a.index, b.index))
radii = np.array(radii)
pairs = np.array(pairs)
distances = md.compute_distances(self.t, pairs, periodic=False)
low, high = self.bond_low * radii, self.bond_high * radii
extreme = np.logical_or(distances < low, distances > high)
if np.any(extreme):
frames, bonds = np.nonzero(extreme)
frame, bond = frames[0], bonds[0]
a1 = self.topology.atom(pairs[bond][0])
a2 = self.topology.atom(pairs[bond][0])
self.log('error: atoms (%s) and (%s) are bonded according to the topology ' % (a1, a2))
self.log('but they are a distance of %.3f nm apart in frame %d' % (distances[frame, bond], frame))
else:
self.log("All good.")
equivalent atoms, respectively.
min_dist : boolean, optional, default: True
Optionally reports the minimum distance between equivalent
atoms. If False, will report the maximum distance.
Returns
----------
dists : array, shape [trj_length, n_pairs]
Distances between specified pairs of atoms
"""
top = trj.topology
apairs = _convert_atom_names(top, apairs)
dists = []
for num in range(len(apairs)):
dpairs = itertools.product(apairs[num][0], apairs[num][1])
dists_tmp = md.compute_distances(trj, atom_pairs=dpairs)
if min_dist:
dists.append(dists_tmp.min(axis=1))
else:
dists.append(dists_tmp.max(axis=1))
return np.array(dists).T
def __init__(self, trajectory, atom_contacts, cutoff=0.45):
atom_contacts = _regularize_contact_input(atom_contacts, "atom")
atom_pairs = [[contact[0][0].index, contact[0][1].index]
for contact in atom_contacts]
labels = [str(contact[0]) for contact in atom_contacts]
distances = md.compute_distances(trajectory, atom_pairs=atom_pairs)
vector_f = np.vectorize(lambda d: d < cutoff)
# transpose because distances is ndarray shape (n_frames,
# n_contacts); values should be list shape (n_contacts, n_frames)
values = vector_f(distances).T.tolist()
super(AtomContactConcurrence, self).__init__(values=values,
labels=labels)
"""
filtered_candidate_pairs = {} # outer dictionary key = time, value = list of tuples (ca_1, ca_2)
ca_keys, ca_indices = [], []
for ca_atom in chain_top.atoms_by_name('CA'):
ca_keys.append(ca_atom)
ca_indices.append(ca_atom.index)
ca_pair_keys, ca_pair_indices = [], []
n = len(ca_indices)
if(n == 0): return
for i1 in range(n):
for i2 in range(i1 + 1, n):
ca_pair_keys.append((ca_keys[i1], ca_keys[i2]))
ca_pair_indices.append([ca_indices[i1], ca_indices[i2]])
ca_pair_indices = np.array(ca_pair_indices)
if(len(ca_pair_indices) == 0): return
pairDistances = md.compute_distances(traj, ca_pair_indices)
for time in range(len(pairDistances)):
t_distances = pairDistances[time]
potential_indices = [i for i in range(len(t_distances)) if t_distances[i] <= ALPHA_CARBON_DIST_CUTOFF]
filtered_candidate_pairs[time] = [ca_pair_keys[i] for i in potential_indices]
filtered_candidate_pairs = fillTimeGaps(range(len(pairDistances)), filtered_candidate_pairs)
return filtered_candidate_pairs
Perform soft distance cutoff
"""
filtered_candidate_pairs = {}
pairKeys = []
atomPairs = []
tempDict = {}
for k1 in cand_aromatic_dict.keys():
for k2 in cand_aromatic_dict.keys():
if(k1 != k2 and (k1, k2) not in tempDict and (k2, k1) not in tempDict):
tempDict[(k1, k2)] = 1
aromatic1 = cand_aromatic_dict[k1]
aromatic2 = cand_aromatic_dict[k2]
generatePairKeys(pairKeys, atomPairs, k1, aromatic1, k2, aromatic2)
# calculate atom pair distances
atomPairs = np.array(atomPairs)
pairDistances = md.compute_distances(traj, atomPairs)
for time in range(len(pairDistances)):
t_distances = pairDistances[time]
for i in range(0, len(t_distances), 9):
arom_arom_distances = t_distances[i:i+9]
if(min(arom_arom_distances) <= soft_cutoff_dist):
arom1_key = pairKeys[i][0][0]
arom2_key = pairKeys[i][1][0]
if(time not in filtered_candidate_pairs.keys()):
filtered_candidate_pairs[time] = {arom1_key:[arom2_key]}
else:
if(arom1_key not in filtered_candidate_pairs[time].keys()):
filtered_candidate_pairs[time][arom1_key] = [arom2_key]
else:
filtered_candidate_pairs[time][arom1_key].append(arom2_key)
filtered_candidate_pairs = fillTimeGaps(range(len(pairDistances)), filtered_candidate_pairs)
return filtered_candidate_pairs