Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
sampler_state: openmmtools.SamplerState
sampler state with positions and box vectors if applicable
"""
#pull a random index
_logger.debug(f"\tpulling a decorrelated trajectory snapshot...")
index = random.choice(self._eq_dict[f"{endstate}_decorrelated"])
_logger.debug(f"\t\tpulled decorrelated index label {index}")
files = [key for key in self._eq_files_dict[endstate].keys() if index in self._eq_files_dict[endstate][key]]
_logger.debug(f"\t\t files corresponding to index {index}: {files}")
assert len(files) == 1, f"files: {files} doesn't have one entry; index: {index}, eq_files_dict: {self._eq_files_dict[endstate]}"
file = files[0]
file_index = self._eq_files_dict[endstate][file].index(index)
_logger.debug(f"\t\tfile_index: {file_index}")
#now we load file as a traj and create a sampler state with it
traj = md.load_frame(file, file_index)
positions = traj.openmm_positions(0)
box_vectors = traj.openmm_boxes(0)
sampler_state = SamplerState(positions, box_vectors = box_vectors)
return sampler_state
def onLoadTrajectory(self):
#self.trajectory_file = str(QtGui.QFileDialog.getOpenFileName(self, 'Open Trajectory-File', '.h5', 'H5-Trajectory-Files (*.h5)'))
filenames = mfm.widgets.open_files('Open Trajectory-File', 'H5-Trajectory-Files (*.h5)')
self.filenames = filenames
self.trajectory_file = filenames[0]
frame0 = md.load_frame(self.trajectory_file, 0)
_, tmp = tempfile.mkstemp(
suffix=".pdb"
)
frame0.save(tmp)
self.topology_file = tmp
self.d1.atoms = self.pdb
self.d2.atoms = self.pdb
self.a1.atoms = self.pdb
self.a2.atoms = self.pdb
@function_timer
def generate_clusters(self, ligand_file):
# Obtain binding site solute atoms using ligand atom coordinates
ligand = md.load_pdb(ligand_file)
ligand_coords = md.utils.in_units_of(ligand.xyz[0, :, :], "nanometers", "angstroms", inplace=True)
first_frame = md.load_frame(self.trajectory, 0, top=self.topology_file)
solute_pos = md.utils.in_units_of(first_frame.xyz[0, self.non_water_atom_ids, :], "nanometers", "angstroms")
search_space = NeighborSearch(solute_pos, 5.0)
near_indices = search_space.query_nbrs_multiple_points(ligand_coords)
binding_site_atom_indices = [self.non_water_atom_ids[nbr_index] for nbr_index in near_indices]
# Obtain water molecules solvating the binding site
stride = 10
print "Reading in trajectory for clustering."
trj = md.load(self.trajectory, top=self.topology)
trj_short = trj[self.start_frame:self.start_frame + trj.n_frames:stride]
print "Obtaining a superconfiguration of all water molecules found in the binding site throught the trajectory."
binding_site_waters = md.compute_neighbors(trj_short, 0.50, binding_site_atom_indices, haystack_indices=self.wat_oxygen_atom_ids)
# generate a list of all waters with their frame ids
water_id_frame_list = [(i, nbr) for i in range(len(binding_site_waters)) for nbr in binding_site_waters[i]]
# Set up clustering loop
print "Performing clustering on the superconfiguration."
cutoff = trj_short.n_frames * 2 * 0.1401
def calculate_angular_structure(self, dist_cutoff=6.0, site_indices=[]):
'''
Returns energetic quantities for each hydration site
'''
pbar = ProgressBar(widgets=[Percentage(), Bar(), ETA()], maxval=self.start_frame + self.num_frames).start()
angular_strucure_data = [[] for site in self.hsa_data]
for i in xrange(self.start_frame, self.start_frame + self.num_frames):
# print "Processing frame: ", i+1, "..."
frame = md.load_frame(self.trajectory, i, top=self.topology)
pos = frame.xyz[0, :, :] * 10.0
pbc = frame.unitcell_lengths[0] * 10.0
#pos = self.trj[i].xyz[0,:,:]*10.0
# obtain coords of O-atoms
oxygen_pos = pos[self.wat_oxygen_atom_ids]
cluster_search_space = NeighborSearch(oxygen_pos, 1.0)
water_search_space = NeighborSearch(oxygen_pos, dist_cutoff)
for site_i, site_data in enumerate(self.hsa_data):
#cluster_center_coords = (site_data["x"], site_data["y"], site_data["z"])
cluster_center_coords = (site_data[1][0], site_data[1][1], site_data[1][2])
nbr_indices = cluster_search_space.query_nbrs_single_point(cluster_center_coords)
cluster_wat_oxygens = [self.wat_oxygen_atom_ids[nbr_index] for nbr_index in nbr_indices]
# begin iterating over water oxygens found in this cluster in
# current frame
process = lambda x, frame: md.load_frame(x, frame, top=top)
else:
f.seek(0, 0)
df = pd.read_csv(f)
if not all(e in df.columns for e in ('filename', 'index', 'state')):
self.error('CSV file not read properly')
for k in np.unique(df['state']):
fn = self.outfn(k)
if os.path.exists(fn):
self.error('IOError: file exists: %s' % fn)
frames = defaultdict(lambda: [])
for fn, group in df.groupby('filename'):
for _, row in group.sort('index').iterrows():
frames[row['state']].append(
md.load_frame(fn, row['index'], top=self.top))
for state, samples in list(frames.items()):
traj = samples[0].join(samples[1:])
print('saving %s...' % self.outfn(state))
traj.save(self.outfn(state), force_overwrite=False)
print('done')
def onJoinTrajectories(self):
target_filename = str(QtWidgets.QFileDialog.getSaveFileName(None, 'Save H5-Model file', '', 'H5-files (*.h5)'))[0]
fn1 = self.trajectory_filename_1
fn2 = self.trajectory_filename_2
r1 = self.reverse_traj_1
r2 = self.reverse_traj_2
traj_1 = md.load_frame(fn1, index=0)
traj_2 = md.load_frame(fn2, index=0)
# Create empty trajectory
if self.join_mode == 'time':
traj_join = traj_1.join(traj_2)
axis = 0
elif self.join_mode == 'atoms':
traj_join = traj_1.stack(traj_2)
axis = 1
target_traj = md.Trajectory(xyz=np.empty((0, traj_join.n_atoms, 3)), topology=traj_join.topology)
target_traj.save(target_filename)
chunk_size = self.chunk_size
table = tables.open_file(target_filename, 'a')
for i, (c1, c2) in enumerate(izip(md.iterload(fn1, chunk=chunk_size), md.iterload(fn2, chunk=chunk_size))):
def onSaveTrajectory(self, target_filename=None):
if target_filename is None:
target_filename = str(QtWidgets.QFileDialog.getSaveFileName(None, 'Save H5-Model file', '', 'H5-files (*.h5)'))[0]
translation_vector = self.translation_vector
rotation_matrix = self.rotation_matrix
stride = self.stride
if self.verbose:
print("Stride: %s" % stride)
print("\nRotation Matrix")
print(rotation_matrix)
print("\nTranslation vector")
print(translation_vector)
first_frame = mdtraj.load_frame(self.trajectory_filename, 0)
traj_new = mdtraj.Trajectory(xyz=np.empty((1, first_frame.n_atoms, 3)), topology=first_frame.topology)
traj_new.save(target_filename)
chunk_size = 1000
table = tables.open_file(target_filename, 'a')
for i, chunk in enumerate(
mdtraj.iterload(
self.trajectory_filename,
chunk=chunk_size,
stride=stride
)
):
xyz = chunk.xyz.copy()
rotate(xyz, rotation_matrix)
translate(xyz, translation_vector)
table.root.xyz.append(xyz)
@function_timer
def calculate_site_quantities(self, energy=True, hbonds=True, entropy=True):
'''
TODO: replace TIP3P nbr count constant with variable that depends on water model
'''
pbar = ProgressBar(widgets=[Percentage(), Bar(), ETA()], maxval=self.start_frame + self.num_frames).start()
for frame_i in xrange(self.start_frame, self.start_frame + self.num_frames):
frame = md.load_frame(self.trajectory, frame_i, top=self.topology)
pos = md.utils.in_units_of(frame.xyz, "nanometers", "angstroms")
pbc = md.utils.in_units_of(frame.unitcell_lengths, "nanometers", "angstroms")
oxygen_pos = pos[0, self.wat_oxygen_atom_ids, :]
cluster_search_space = NeighborSearch(oxygen_pos, 1.0)
water_search_space = NeighborSearch(oxygen_pos, 3.5)
for site_i in xrange(self.hsa_data.shape[0]):
wat_O = None
if self.site_waters is not None:
if len(self.site_waters[site_i]) != 0:
if self.site_waters[site_i][0][0] == frame_i:
wat_O = self.site_waters[site_i].pop(0)[1]
self.hsa_data[site_i, 4] += 1
else:
cluster_center_coords = (self.hsa_data[site_i, 1], self.hsa_data[site_i, 2], self.hsa_data[site_i, 3])
nbr_indices = cluster_search_space.query_nbrs_single_point(cluster_center_coords)
cluster_wat_oxygens = [self.wat_oxygen_atom_ids[nbr_index] for nbr_index in nbr_indices]