Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
import mdtraj as md
import pandas as pd
except ImportError as e:
print("Error: generate_dummy_trajectory() requires mdtraj and pandas!")
raise(e)
n_atoms = len(xyz)
data = []
for i in range(n_atoms):
data.append(dict(serial=i, name="H", element="H", resSeq=i + 1, resName="UNK", chainID=0))
data = pd.DataFrame(data)
unitcell_lengths = box * np.ones((1, 3))
unitcell_angles = 90 * np.ones((1, 3))
top = md.Topology.from_dataframe(data, np.zeros((0, 2), dtype='int'))
traj = md.Trajectory(xyz, top, unitcell_lengths=unitcell_lengths, unitcell_angles=unitcell_angles)
return traj
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
init_sim_state = context.getState(**get_state_kwargs)
thermostat = omm.AndersenThermostat(300.0 * unit.kelvin, 1/unit.picosecond)
barostat = omm.MonteCarloBarostat(1.0*unit.atmosphere, 300.0*unit.kelvin, 50)
runner = OpenMMRunner(test_sys.system, test_sys.topology, integrator, platform='CUDA')
num_walkers = 10
init_weight = 1.0 / num_walkers
init_walkers = [Walker(OpenMMState(init_sim_state), init_weight) for i in range(num_walkers)]
# the mdtraj here is needed for the distance function
mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)
# make a distance object which can be used to compute the distance
# between two walkers, for our scorer class
distance = PairDistance()
# make a WExplore2 resampler with default parameters and our
# distance metric
resampler = WExplore1Resampler(max_region_sizes=(0.5, 0.5, 0.5, 0.5),
max_n_regions=(10, 10, 10, 10),
distance=distance,
pmin=1e-12, pmax=0.5)
ubc = UnbindingBC(cutoff_distance=2.0,
initial_state=init_walkers[0].state,
topology=mdtraj_topology,
ligand_idxs=np.array(test_sys.ligand_indices),
protein_pdbfile = open(protein_filename, 'r')
protein_pdb = app.PDBFile(protein_pdbfile)
protein_pdbfile.close()
protein_positions, protein_topology, protein_md_topology = protein_pdb.positions, protein_pdb.topology, md.Topology.from_openmm(protein_pdb.topology)
protein_topology = protein_md_topology.to_openmm()
protein_n_atoms = protein_md_topology.n_atoms
# Load the ligand, if present
molecules = []
if ligand_file:
if ligand_file.endswith('.sdf'): # small molecule
ligand_mol = createOEMolFromSDF(ligand_file, index=ligand_index)
ligand_mol = generate_unique_atom_names(ligand_mol)
molecules.append(Molecule.from_openeye(ligand_mol, allow_undefined_stereo=False))
ligand_positions, ligand_topology = extractPositionsFromOEMol(ligand_mol), forcefield_generators.generateTopologyFromOEMol(ligand_mol)
ligand_md_topology = md.Topology.from_openmm(ligand_topology)
ligand_n_atoms = ligand_md_topology.n_atoms
elif ligand_file.endswith('pdb'): # protein
ligand_pdbfile = open(ligand_file, 'r')
ligand_pdb = app.PDBFile(ligand_pdbfile)
ligand_pdbfile.close()
ligand_positions, ligand_topology, ligand_md_topology = ligand_pdb.positions, ligand_pdb.topology, md.Topology.from_openmm(
ligand_pdb.topology)
ligand_n_atoms = ligand_md_topology.n_atoms
else:
_logger.warning(f'ligand filetype not recognised. Please provide a path to a .pdb or .sdf file')
return
# Now create a complex
complex_md_topology = protein_md_topology.join(ligand_md_topology)
hdf5_filename = 'tmp.wepy.h5'
wepy_h5 = WepyHDF5(hdf5_filename, mode='w', topology=top_json)
# make the test system from openmmtools
test_sys = LennardJonesPair()
# integrator
TEMPERATURE= 300.0*unit.kelvin
FRICTION_COEFFICIENT = 1/unit.picosecond
STEP_SIZE = 0.002*unit.picoseconds
integrator = omm.LangevinIntegrator(TEMPERATURE, FRICTION_COEFFICIENT, STEP_SIZE)
# the mdtraj here is needed for unbininding BC
mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)
# the initial state for the BC
context = omm.Context(test_sys.system, copy(integrator))
context.setPositions(test_sys.positions)
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
init_sim_state = context.getState(**get_state_kwargs)
init_state = OpenMMState(init_sim_state)
# initialize the unbinding boundary conditions
ubc = UnbindingBC(cutoff_distance=0.5,
initial_state=init_state,
topology=mdtraj_topology,
ligand_idxs=np.array(test_sys.ligand_indices),
binding_site_idxs=np.array(test_sys.receptor_indices))
def restrain_atoms_by_dsl(thermodynamic_state, sampler_state, topology, atoms_dsl, **kwargs):
# Make sure the topology is an MDTraj topology.
if isinstance(topology, mdtraj.Topology):
mdtraj_topology = topology
else:
mdtraj_topology = mdtraj.Topology.from_openmm(topology)
# Determine indices of the atoms to restrain.
restrained_atoms = mdtraj_topology.select(atoms_dsl).tolist()
restrain_atoms(thermodynamic_state, sampler_state, restrained_atoms, **kwargs)
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
init_sim_state = context.getState(**get_state_kwargs)
thermostat = omm.AndersenThermostat(300.0 * unit.kelvin, 1/unit.picosecond)
barostat = omm.MonteCarloBarostat(1.0*unit.atmosphere, 300.0*unit.kelvin, 50)
runner = OpenMMRunner(test_sys.system, test_sys.topology, integrator, platform='Reference')
num_walkers = 10
init_weight = 1.0 / num_walkers
init_walkers = [Walker(OpenMMState(init_sim_state), init_weight) for i in range(num_walkers)]
# the mdtraj here is needed for the distance function
mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)
# make a distance object which can be used to compute the distance
# between two walkers, for our scorer class
distance = PairDistance()
# we need a scorer class to perform the all-to-all distances
# using our distance class
scorer = AllToAllScorer(distance=distance)
# make a WExplore2 resampler with default parameters and our
# distance metric
resampler = WExplore2Resampler(scorer=scorer,
pmax=0.5)
ubc = UnbindingBC(cutoff_distance=1.5,
initial_state=init_walkers[0].state,
def _topology_from_arrays(AtomID, AtomNames, ChainID, ResidueID, ResidueNames):
"""Build topology object from the arrays stored in the lh5 file"""
# Delayed import due to wacky recursive imports in compatibilty
from mdtraj import Topology
topology = Topology()
# assert that the ChainID is just an array of empty strings, which appears
# to be the case in our test systems for this legacy format
if not np.all(chainid == '' for chainid in ChainID):
raise NotImplementedError('Im not prepared to parse multiple chains')
chain0 = topology.add_chain()
# register the residues
registered_residues = {}
for i in np.argsort(ResidueID):
residue_name = ResidueNames[i]
if not isinstance(residue_name, basestring):
residue_name = residue_name.decode()
if ResidueID[i] not in registered_residues:
res = topology.add_residue(residue_name, chain0)
registered_residues[ResidueID[i]] = res
def topologyfile(self, topfile):
self._topologyfile = topfile
if isinstance(topfile, str):
self.topology = load_topology_cached(topfile)
self._topologyfile = topfile
elif isinstance(topfile, mdtraj.Topology):
self.topology = topfile
elif isinstance(topfile, mdtraj.Trajectory):
self.topology = topfile.topology
else:
raise ValueError("no valid topfile arg: type was %s, "
"but only string or mdtraj.Topology allowed." % type(topfile))
old_ligand_topology = old_complex.subset(old_complex.select("resname == 'MOL' "))
new_ligand_topology = new_complex.subset(new_complex.select("resname == 'MOL' "))
# solvate the old ligand topology:
old_solvated_topology, old_solvated_positions, old_solvated_system = self._solvate(old_ligand_topology.to_openmm(), old_ligand_positions)
old_solvated_md_topology = md.Topology.from_openmm(old_solvated_topology)
# now remove the old ligand, leaving only the solvent
solvent_only_topology = old_solvated_md_topology.subset(old_solvated_md_topology.select("not resname MOL"))
# append the solvent to the new ligand-only topology:
new_solvated_ligand_md_topology = new_ligand_topology.join(solvent_only_topology)
nsl, b = new_solvated_ligand_md_topology.to_dataframe()
# dirty hack because new_solvated_ligand_md_topology.to_openmm() was throwing bond topology error
new_solvated_ligand_md_topology = md.Topology.from_dataframe(nsl, b)
new_solvated_ligand_omm_topology = new_solvated_ligand_md_topology.to_openmm()
new_solvated_ligand_omm_topology.setPeriodicBoxVectors(old_solvated_topology.getPeriodicBoxVectors())
# create the new ligand system:
new_solvated_system = self.system_generator.create_system(new_solvated_ligand_omm_topology)
new_to_old_atom_map = {atom_map[x] - new_mol_start_index: x - old_mol_start_index for x in
old_complex.select("resname == 'MOL' ") if x in atom_map.keys()}
# adjust the atom map to account for the presence of solvent degrees of freedom:
# By design, all atoms after the ligands are water, and should be mapped.
n_water_atoms = solvent_only_topology.to_openmm().getNumAtoms()
for i in range(n_water_atoms):
new_to_old_atom_map[new_mol_len + i] = old_mol_len + i