Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def test_load_pdb_no_standard_names(get_fn):
# Minimal test. Standard_names=False will force load_pdb.py
# to NOT replace any non-standard atom or residue names in the topology
md.load(get_fn('native2.pdb'), standard_names=False, no_boxchk=True)
md.load_pdb(get_fn('native2.pdb'), standard_names=False, no_boxchk=True)
for assoc_term in assoc_terms:
# make the unit AssociationTypes
assoc_idx = sEH_TPPU_SystemType.make_unit_association_type(assoc_term)
association_type = sEH_TPPU_SystemType.association_types[assoc_idx]
# make HydrogenBondType interaction classes for this association
# in the inx_space
inx_space.add_association_subspace(association_type, HydrogenBondType)
print("Building profiler..")
# make a Profiler for the inx space
profiler = masticprof.InxSpaceProfiler(inx_space)
pdbpath = 'data/sEH_TPPU_system.pdb'
pdb = mdj.load_pdb(pdbpath)
lig_idxs = pdb.topology.select('resname "2RV"')
protein_idxs = np.array([atom.index for atom in pdb.topology.atoms if atom.residue.is_protein])
print("Setting up distance metric..")
hbond_distance = OpenMMHBondDistance(ligand_idxs=lig_idxs,
protein_idxs=protein_idxs,
profiler=profiler,
sys_type=sEH_TPPU_SystemType)
print("Loading walkers..")
walker_pkl_path = 'data/sEH_walkers.pkl'
with open(walker_pkl_path, 'rb') as rf:
sEH_walkers = pickle.load(rf)
print("Computing distance..")
dmat = hbond_distance.distance(sEH_walkers)
import mdtraj as mdj
from wepy.sim_manager import Manager
from wepy.resampling.clone_merge import CloneMergeDecision, CLONE_MERGE_INSTRUCTION_DTYPES
from wepy.resampling.clone_merge import RandomCloneMergeResampler
from wepy.openmm import OpenMMRunner, OpenMMWalker
from wepy.reporter.hdf5 import WepyHDF5Reporter
from wepy.resampling.clone_merge import clone_parent_table, clone_parent_panel
if __name__ == "__main__":
# write out an hdf5 storage of the system
mdj_pdb = mdj.load_pdb('sEH_TPPU_system.pdb')
mdj_pdb.save_hdf5('sEH_TPPU_system.h5')
# set up the system which will be passed to the arguments of the mapped functions
# the topology from the PSF
psf = omma.CharmmPsfFile('sEH_TPPU_system.psf')
# load the coordinates
omm_pdb = omma.PDBFile('sEH_TPPU_system.pdb')
# to use charmm forcefields get your parameters
params = omma.CharmmParameterSet('all36_cgenff.rtf',
'all36_cgenff.prm',
'all36_prot.rtf',
'all36_prot.prm',
'tppu.str',
'toppar_water_ions.str')
@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
def _construct_traj(self):
logger.debug('Loading Trajectory object for model {0} ({1}/{2})'.format(self.df.templateid.iloc[0], 0, len(self.df.model_filepath)))
traj = mdtraj.load_pdb(self.df.model_filepath[0])
remove_disulfide_bonds_from_topology(traj.topology)
self.traj = traj
for m, model_filepath in enumerate(self.df.model_filepath[1:]):
logger.debug('Loading Trajectory object for model {0} ({1}/{2})'.format(self.df.templateid.iloc[m+1], m+1, len(self.df.model_filepath)))
traj = mdtraj.load_pdb(model_filepath)
remove_disulfide_bonds_from_topology(traj.topology)
self.traj += traj
def _align_structure(self, structure_filepath, alignment_index):
traj = mdtraj.load_pdb(structure_filepath)
traj.superpose(self._ref_traj, atom_indices=self._heavy_atoms, parallel=False)
aligned_structure_filepath = os.path.join(self._tmpdir, 'aligned%d.pdb' % alignment_index)
traj.save(aligned_structure_filepath)
def _get_models(self):
self.model = {}
root, dirnames, filenames = next(os.walk(self.models_target_dir))
for dirname in dirnames:
if 'implicit' in self.model and 'explicit' in self.model:
break
if 'implicit' not in self.model:
implicit_model_filename = os.path.join(self.models_target_dir, dirname, 'implicit-refined.pdb.gz')
if os.path.exists(implicit_model_filename):
self.model['implicit'] = mdtraj.load_pdb(implicit_model_filename)
if 'explicit' not in self.model:
explicit_model_filename = os.path.join(self.models_target_dir, dirname, 'explicit-refined.pdb.gz')
if os.path.exists(explicit_model_filename):
self.model['explicit'] = mdtraj.load_pdb(explicit_model_filename)
current_wat = water_id_frame_list[max_index]
current_wat_coords = md.utils.in_units_of(coords[current_wat[0], current_wat[1], :],
"nanometers", "angstroms")
near_flag = 0
if len(cluster_list) != 0:
for clust in cluster_list:
clust_coords = coords[clust[0], clust[1], :]
dist = np.linalg.norm(current_wat_coords - clust_coords)
if dist < 1.20:
near_flag += 1
if near_flag == 0:
cluster_iter += 1
cluster_list.append(water_id_frame_list[max_index])
init_cluster_coords = [coords[cluster[0], cluster[1], :] for cluster in cluster_list]
else:
clusters_pdb_file = md.load_pdb(clustercenter_file, no_boxchk=True)
init_cluster_coords = clusters_pdb_file.xyz[0, :, :]
# Read full trajectory
print("Reading trajectory to obtain water molecules for each cluster.")
with md.open(self.trajectory) as f:
f.seek(self.start_frame)
if self.num_frames is None:
trj = f.read_as_traj(topology, stride=1,
atom_indices=np.concatenate((binding_site_atom_indices, self.wat_oxygen_atom_ids)))
self.num_frames = trj.n_frames
else:
trj = f.read_as_traj(topology, n_frames=self.num_frames, stride=1,
atom_indices=np.concatenate((binding_site_atom_indices, self.wat_oxygen_atom_ids)))
if trj.n_frames < self.num_frames:
print(("Warning: {0:d} frames found in the trajectory, resetting self.num_frames.".format(
trj.n_frames)))
def write_renumbered_structure(source_filepath, dest_filepath, renumbered_resnums):
traj = mdtraj.load_pdb(source_filepath)
for r, residue in enumerate(traj.top.residues):
residue.resSeq = renumbered_resnums[r]
traj.save_pdb(dest_filepath)
in the corresponding lists.
self.hsa_region_flat_ids:
Same as above except that indices are not atom indices from the topology
but in a sequence from 0 to N, where N is the total number of water oxygen atoms found in the
HSA region throughout the simulation.
self.hsa_region_water_coords:
An N x 3 numpy array is initialized, where N is the total number of water water oxygen atoms found in the
HSA region throughout the simulation. The array gets populated during individual frame processing.
"""
sphere_radius = md.utils.in_units_of(1.0, "angstroms", "nanometers")
topology = md.load_topology(self.topology_file)
if self.non_water_atom_ids.shape[0] == 0:
raise Exception(ValueError,
"Clustering is supported only for solute-solvent systems, no solute atoms found.")
ligand = md.load_pdb(ligand_file, no_boxchk=True)
ligand_coords = ligand.xyz[0, :, :]
binding_site_atom_indices = np.asarray(list(range(ligand_coords.shape[0])))
init_cluster_coords = None
# Step 1: Initial Clustering if user didn't provide cluster centers
if clustercenter_file is None:
clustering_stride = 10
print("Reading trajectory for clustering.")
with md.open(self.trajectory) as f:
f.seek(self.start_frame)
# read all frames if no frames specified by user
if self.num_frames is None:
trj_short = f.read_as_traj(topology, stride=clustering_stride,
atom_indices=np.concatenate(
(binding_site_atom_indices, self.wat_oxygen_atom_ids)))
else:
trj_short = f.read_as_traj(topology, n_frames=self.num_frames,