Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
for (index, atom) in enumerate(mol.GetAtoms()):
[x,y,z] = mol.GetCoords(atom)
positions[index,0] = x*unit.angstroms
positions[index,1] = y*unit.angstroms
positions[index,2] = z*unit.angstroms
#Initialize SMIRFF parameters
# using smarty ForceField object
smirnoff_xml_filename ="smirff99Frosst.ffxml"
forcefield = ForceField(smirnoff_xml_filename)
#Create OpenMM System
system = forcefield.createSystem(topology, [mol])
#Load System to ParmEd Structure
structure = pmd.openmm.load_topology(topology, system)
return structure
The new positions of the system
logp_proposal : float
The log probability of the forward-only proposal
"""
_logger.info("propose: performing forward proposal")
# Ensure positions have units compatible with nanometers
check_dimensionality(current_positions, unit.nanometers)
check_dimensionality(beta, unit.kilojoules_per_mole**(-1))
# TODO: Change this to use md_unit_system instead of hard-coding nanometers
if not top_proposal.unique_new_atoms:
_logger.info("propose: there are no unique new atoms; logp_proposal = 0.0.")
# If there are no unique new atoms, return new positions in correct order for new topology object and log probability of zero
# TODO: Carefully check this
import parmed
structure = parmed.openmm.load_topology(top_proposal.old_topology, top_proposal._old_system)
atoms_with_positions = [ structure.atoms[atom_idx] for atom_idx in top_proposal.new_to_old_atom_map.keys() ]
new_positions = self._copy_positions(atoms_with_positions, top_proposal, current_positions)
logp_proposal, rjmc_info, atoms_with_positions_reduced_potential, final_context_reduced_potential, neglected_angle_terms = 0.0, None, None, None, None
self.forward_final_growth_system = None
else:
_logger.info("propose: unique new atoms detected; proceeding to _logp_propose...")
logp_proposal, new_positions, rjmc_info, atoms_with_positions_reduced_potential, final_context_reduced_potential, neglected_angle_terms = self._logp_propose(top_proposal, current_positions, beta, direction='forward')
self.nproposed += 1
check_dimensionality(new_positions, unit.nanometers)
check_dimensionality(logp_proposal, float)
#define forward attributes
self.forward_rjmc_info = rjmc_info
self.forward_atoms_with_positions_reduced_potential, self.forward_final_context_reduced_potential = atoms_with_positions_reduced_potential, final_context_reduced_potential
self.forward_neglected_angle_terms = neglected_angle_terms
check_dimensionality(new_positions, unit.angstroms)
# Determine order in which atoms (and the torsions they are involved in) will be proposed
_logger.info("Computing proposal order with NetworkX...")
proposal_order_tool = NetworkXProposalOrder(top_proposal, direction=direction)
torsion_proposal_order, logp_choice = proposal_order_tool.determine_proposal_order()
atom_proposal_order = [ torsion[0] for torsion in torsion_proposal_order ]
_logger.info(f"number of atoms to be placed: {len(atom_proposal_order)}")
_logger.info(f"Atom index proposal order is {atom_proposal_order}")
growth_parameter_name = 'growth_stage'
if direction=="forward":
_logger.info("direction of proposal is forward; creating atoms_with_positions and new positions from old system/topology...")
# Find and copy known positions to match new topology
import parmed
structure = parmed.openmm.load_topology(top_proposal.new_topology, top_proposal.new_system)
atoms_with_positions = [structure.atoms[atom_idx] for atom_idx in top_proposal.new_to_old_atom_map.keys()]
new_positions = self._copy_positions(atoms_with_positions, top_proposal, old_positions)
self._new_posits = copy.deepcopy(new_positions)
# Create modified System object
_logger.info("creating growth system...")
growth_system_generator = GeometrySystemGenerator(top_proposal.new_system, torsion_proposal_order, global_parameter_name=growth_parameter_name, reference_topology=top_proposal.new_topology, use_sterics=self.use_sterics, neglect_angles = self.neglect_angles, use_14_nonbondeds = self._use_14_nonbondeds)
growth_system = growth_system_generator.get_modified_system()
elif direction=='reverse':
_logger.info("direction of proposal is reverse; creating atoms_with_positions from old system/topology")
if new_positions is None:
raise ValueError("For reverse proposals, new_positions must not be none.")
# Find and copy known positions to match old topology
import parmed
new_positions = self._copy_positions(atoms_with_positions, top_proposal, old_positions)
self._new_posits = copy.deepcopy(new_positions)
# Create modified System object
_logger.info("creating growth system...")
growth_system_generator = GeometrySystemGenerator(top_proposal.new_system, torsion_proposal_order, global_parameter_name=growth_parameter_name, reference_topology=top_proposal.new_topology, use_sterics=self.use_sterics, neglect_angles = self.neglect_angles, use_14_nonbondeds = self._use_14_nonbondeds)
growth_system = growth_system_generator.get_modified_system()
elif direction=='reverse':
_logger.info("direction of proposal is reverse; creating atoms_with_positions from old system/topology")
if new_positions is None:
raise ValueError("For reverse proposals, new_positions must not be none.")
# Find and copy known positions to match old topology
import parmed
structure = parmed.openmm.load_topology(top_proposal.old_topology, top_proposal.old_system)
atoms_with_positions = [structure.atoms[atom_idx] for atom_idx in top_proposal.old_to_new_atom_map.keys()]
# Create modified System object
_logger.info("creating growth system...")
growth_system_generator = GeometrySystemGenerator(top_proposal.old_system, torsion_proposal_order, global_parameter_name=growth_parameter_name, reference_topology=top_proposal.old_topology, use_sterics=self.use_sterics, neglect_angles = self.neglect_angles, use_14_nonbondeds = self._use_14_nonbondeds)
growth_system = growth_system_generator.get_modified_system()
else:
raise ValueError("Parameter 'direction' must be forward or reverse")
# Define a system for the core atoms before new atoms are placed
self.atoms_with_positions_system = growth_system_generator._atoms_with_positions_system
self.growth_system = growth_system
# Get the angle terms that are neglected from the growth system
neglected_angle_terms = growth_system_generator.neglected_angle_terms
_logger.info(f"neglected angle terms include {neglected_angle_terms}")