diff --git a/.travis.yml b/.travis.yml index 9e23edb35..b07ed7616 100644 --- a/.travis.yml +++ b/.travis.yml @@ -58,6 +58,13 @@ script: - conda install --yes --use-local ${PACKAGENAME}-dev # Install testing dependencies - conda install --yes --quiet nose nose-timer + # Install openmm-forcefields via pip (until available via omnia) + - conda install --yes tinydb "openforcefield>=0.6.0" "openforcefields>=1.0.0" + - pip install git+https://github.com/openmm/openmm-forcefields.git + # Workaround for https://github.com/openmm/openmm/pull/2511 + # TODO: Remove this when OpenMM builds including https://github.com/openmm/openmm/pull/2511 are released + - echo "Overwriting OpenMM forcefield.py with fix from https://github.com/openmm/openmm/pull/2511" + - wget -q https://raw.githubusercontent.com/openmm/openmm/4f48402f1be3e0f049ae0e8595db638d297b0d75/wrappers/python/simtk/openmm/app/forcefield.py -O `python -c "from simtk.openmm.app import forcefield; print(forcefield.__file__)"` # Install desired OpenMM version - if [ "$OPENMM" == "latest" ]; then echo "Using latest release OpenMM."; conda install --yes -c omnia openmm; fi - if [ "$OPENMM" == "beta" ]; then echo "Using OpenMM beta"; conda install --yes -c omnia/label/beta openmm; fi diff --git a/perses/app/experiments.py b/perses/app/experiments.py index e10f7e534..e0efe4c0e 100644 --- a/perses/app/experiments.py +++ b/perses/app/experiments.py @@ -11,7 +11,7 @@ from perses.annihilation.relative import HybridTopologyFactory from perses.app.relative_setup import NonequilibriumSwitchingFEP, RelativeFEPSetup from perses.annihilation.lambda_protocol import LambdaProtocol -from perses.rjmc.topology_proposal import TopologyProposal, SystemGenerator,SmallMoleculeSetProposalEngine +from perses.rjmc.topology_proposal import TopologyProposal, SmallMoleculeSetProposalEngine from perses.rjmc.geometry import FFAllAngleGeometryEngine from perses.app.simulation import Simulation @@ -52,7 +52,9 @@ class BuildProposalNetwork(object): 'hmass': 4 * unit.amus, 'map_strength': 'default', 'phases': ['vacuum', 'solvent', 'complex'], - 'forcefield_files': ['gaff.xml', 'amber14/protein.ff14SB.xml', 'amber14/tip3p.xml'], + 'forcefield_files': ['amber/protein.ff14SB.xml', 'amber/tip3p_standard.xml', 'amber/tip3p_HFE_multivalent.xml'], + 'small_molecule_forcefield' : 'openff-1.0.0', + 'small_molecule_parameters_cache' : 'cache.json', 'neglect_angles': False, 'anneal_14s': False, 'water_model': 'tip3p', @@ -75,7 +77,7 @@ class BuildProposalNetwork(object): 'n_steps_per_move_application': 1, 'reassign_velocities': False, 'n_restart_attempts': 20, - 'splitting': "V R R R O R R R V", + 'splitting': "V R O R V", 'constraint_tolerance' : 1e-6, 'offline_freq': 10, @@ -631,19 +633,19 @@ def _parse_ligand_input(self): self.smiles_list = load_smi(self.ligand_input, index = self.ligand_indices) #create a ligand data list to hold all ligand oemols, systems, positions, topologies + from openforcefield.topology import Molecule + self.molecules = list() # openforcefield Molecule objects for smiles in self.smiles_list: _logger.debug(f"creating oemol, system, positions, and openmm.Topology for smiles: {smiles}...") + # TODO: Should we avoid createSystemFromSMILES? oemol, system, positions, topology = createSystemFromSMILES(smiles, title=smiles) self.ligand_oemol_pos_top.append([oemol, positions]) - - #pull all of the oemols (in order) to make an appropriate ffxml - mol_list = [_tuple[0] for _tuple in self.ligand_oemol_pos_top] - self.ligand_ffxml = forcefield_generators.generateForceFieldFromMolecules(mol_list) + self.molecules.append(Molecule.from_openeye(oemol)) #now make all of the oemol titles 'MOL' [self.ligand_oemol_pos_top[i][0].SetTitle("MOL") for i in range(len(self.ligand_oemol_pos_top))] - #the last thing to do is to make ligand topologies + #make ligand topologies ligand_topologies = [forcefield_generators.generateTopologyFromOEMol(data[0]) for data in self.ligand_oemol_pos_top] [self.ligand_oemol_pos_top[i].append(topology) for i, topology in enumerate(ligand_topologies)] @@ -656,6 +658,9 @@ def _parse_ligand_input(self): [oemol.SetTitle("MOL") for oemol in oemols] self.smiles_list = [ oechem.OECreateSmiString(oemol, oechem.OESMILESFlag_DEFAULT | oechem.OESMILESFlag_Hydrogens) for oemol in oemols] self.ligand_oemol_pos_top = [[oemol, position, forcefield_generators.generateTopologyFromOEMol(oemol)] for oemol, position in zip(oemols, positions)] + # Create openforcefield Molecule objects + from openforcefield.topology import Molecule + self.molecules = [ Molecule.from_openeye(oemol) for oemol in oemols ] else: raise Exception(f"the ligand input can only be a string pointing to an .sdf or .smi file. Aborting!") @@ -783,27 +788,26 @@ def _create_system_generator(self): """ Wrap the process for generating a dict of system generators for each phase. """ + from openmmforcefields.generators import SystemGenerator + + barostat = None if self.proposal_arguments['pressure'] is not None: if self.nonbonded_method == app.PME: barostat = openmm.MonteCarloBarostat(self.proposal_arguments['pressure'], self.proposal_arguments['temperature'], 50) - else: - barostat = None - self.system_generator = SystemGenerator(self.proposal_arguments['forcefield_files'], - barostat=barostat, - forcefield_kwargs={'removeCMMotion': False, - 'nonbondedMethod': self.nonbonded_method, - 'constraints' : app.HBonds, - 'hydrogenMass' : self.proposal_arguments['hmass']}) - else: - self.system_generator = SystemGenerator(forcefield_files, - forcefield_kwargs={'removeCMMotion': False, - 'nonbondedMethod': self.nonbonded_method, - 'constraints' : app.HBonds, - 'hydrogenMass' : self.proposal_arguments['hmass']}) - self.system_generator._forcefield.loadFile(StringIO(self.ligand_ffxml)) + forcefield_kwargs = {'removeCMMotion': False, + 'nonbondedMethod': self.nonbonded_method, + 'constraints' : app.HBonds, + 'hydrogenMass' : self.proposal_arguments['hmass']} + + self.system_generator = SystemGenerator(forcefields=self.proposal_arguments['forcefield_files'], + small_molecule_forcefield=self.proposal_arguments['small_molecule_forcefield'], + cache=self.proposal_arguments['small_molecule_parameters_cache'], + molecules=self.molecules, + barostat=barostat, + forcefield_kwargs=forcefield_kwargs) def _setup_complex_phase(self, ligand_oemol, ligand_positions, ligand_topology): """ @@ -878,10 +882,10 @@ def _solvate(self, topology, positions, model = 'tip3p', vacuum = False): modeller = app.Modeller(topology, positions) hs = [atom for atom in modeller.topology.atoms() if atom.element.symbol in ['H'] and atom.residue.name not in ['MOL','OLD','NEW']] modeller.delete(hs) - modeller.addHydrogens(forcefield = self.system_generator._forcefield) + modeller.addHydrogens(forcefield = self.system_generator.forcefield) if not vacuum: _logger.info(f"\tpreparing to add solvent") - modeller.addSolvent(self.system_generator._forcefield, + modeller.addSolvent(self.system_generator.forcefield, model=model, padding = self.proposal_arguments['solvent_padding'], ionicStrength = 0.15*unit.molar) @@ -891,7 +895,7 @@ def _solvate(self, topology, positions, model = 'tip3p', vacuum = False): solvated_positions = modeller.getPositions() solvated_positions = unit.quantity.Quantity(value = np.array([list(atom_pos) for atom_pos in solvated_positions.value_in_unit_system(unit.md_unit_system)]), unit = unit.nanometers) - solvated_system = self.system_generator.build_system(solvated_topology) + solvated_system = self.system_generator.create_system(solvated_topology) return solvated_topology, solvated_positions, solvated_system def _generate_solvent_topologies(self, topology_proposal, old_positions): @@ -948,7 +952,7 @@ def _generate_solvent_topologies(self, topology_proposal, old_positions): new_solvated_ligand_omm_topology.setPeriodicBoxVectors(old_solvated_topology.getPeriodicBoxVectors()) # create the new ligand system: - new_solvated_system = self.system_generator.build_system(new_solvated_ligand_omm_topology) + 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()} @@ -1011,8 +1015,8 @@ def _generate_vacuum_topologies(self, topology_proposal, old_positions, system_g new_ligand_topology = new_ligand_topology.to_openmm() # create the new ligand system: - old_ligand_system = system_generator.build_system(old_ligand_topology) - new_ligand_system = system_generator.build_system(new_ligand_topology) + old_ligand_system = system_generator.create_system(old_ligand_topology) + new_ligand_system = system_generator.create_system(new_ligand_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()} @@ -1154,11 +1158,12 @@ def _generate_proposals(self, current_oemol, proposed_oemol, current_positions, if 'vacuum' in self.proposal_arguments['phases']: _logger.debug(f"\t\tvacuum:") - vacuum_system_generator = SystemGenerator(self.proposal_arguments['forcefield_files'], - forcefield_kwargs={'removeCMMotion': False, - 'nonbondedMethod': app.NoCutoff, - 'constraints' : app.HBonds}) - vacuum_system_generator._forcefield.loadFile(StringIO(self.ligand_ffxml)) + # Create modified SystemGenerator for vacuum systems + # TODO: Should we instead allow SystemGenerator.create_system() to override kwargs? + # Or perhaps it should recognize whether the Topology object is periodic or not? + vacuum_system_generator = copy.deepcopy(self.system_generator) + vacuum_system_generator.forcefield_kwargs['nonbondedMethod'] = app.NoCutoff + if 'complex' not in self.proposal_arguments['phases'] and 'solvent' not in self.proposal_arguments['phases']: vacuum_ligand_topology, vacuum_ligand_positions, vacuum_ligand_system = self._solvate(current_topology, current_positions, diff --git a/perses/tests/test.smi b/perses/data/test.smi similarity index 100% rename from perses/tests/test.smi rename to perses/data/test.smi diff --git a/perses/dispersed/smc.py b/perses/dispersed/smc.py index fbd11e5a4..31174a4a0 100644 --- a/perses/dispersed/smc.py +++ b/perses/dispersed/smc.py @@ -182,7 +182,7 @@ def __init__(self, #instantiate nonequilibrium work dicts: the keys indicate from which equilibrium thermodynamic state the neq_switching is conducted FROM (as opposed to TO) self.cumulative_work = {'forward': [], 'reverse': []} - self._neq_timers = copy.deepcopy() + self._neq_timers = copy.deepcopy(self._eq_timers) self.incremental_work = copy.deepcopy(self.cumulative_work) self.shadow_work = copy.deepcopy(self.cumulative_work) self.total_jobs = {'forward': 0, 'reverse': 0} #the total number of jobs that have been run in either direction diff --git a/perses/rjmc/topology_proposal.py b/perses/rjmc/topology_proposal.py index b2f7c3b0a..a54e62d91 100644 --- a/perses/rjmc/topology_proposal.py +++ b/perses/rjmc/topology_proposal.py @@ -1077,8 +1077,11 @@ def propose(self, current_system, current_topology, current_metadata=None): # Copy periodic box vectors from current topology new_topology.setPeriodicBoxVectors(current_topology.getPeriodicBoxVectors()) - # Build system - new_system = self._system_generator.build_system(new_topology) + # TODO: Remove build_system() branch once we convert entirely to new openmm-forcefields SystemBuilder + if hasattr(self._system_generator, 'create_system'): + new_system = self._system_generator.create_system(new_topology) + else: + new_system = self._system_generator.build_system(new_topology) # Adjust logp_propose based on HIS presence his_residues = ['HID', 'HIE'] @@ -2439,7 +2442,11 @@ def propose(self, current_system, current_topology, current_mol=None, proposed_m # Generate an OpenMM System from the proposed Topology _logger.info(f"proceeding to build the new system from the new topology...") - new_system = self._system_generator.build_system(new_topology) + # TODO: Remove build_system() branch once we convert entirely to new openmm-forcefields SystemBuilder + if hasattr(self._system_generator, 'create_system'): + new_system = self._system_generator.create_system(new_topology) + else: + new_system = self._system_generator.build_system(new_topology) # Determine atom mapping between old and new molecules _logger.info(f"determining atom map between old and new molecules...") @@ -3130,7 +3137,11 @@ def propose(self, current_system, current_topology, current_smiles=None, propose new_mol_start_index, len_new_mol = self._find_mol_start_index(new_topology) # Generate an OpenMM System from the proposed Topology - new_system = self._system_generator.build_system(new_topology) + # TODO: Remove build_system() branch once we convert entirely to new openmm-forcefields SystemBuilder + if hasattr(self._system_generator, 'create_system'): + new_system = self._system_generator.create_system(new_topology) + else: + new_system = self._system_generator.build_system(new_topology) # Determine atom mapping between old and new molecules mol_atom_maps = self._atom_mapper.get_atom_maps(current_mol_smiles, proposed_mol_smiles) diff --git a/perses/tests/test_experiments.py b/perses/tests/test_experiments.py index 3b705fa07..bfa2062bc 100644 --- a/perses/tests/test_experiments.py +++ b/perses/tests/test_experiments.py @@ -57,21 +57,27 @@ def test_BuildProposalNetwork(): # resources = None, # proposal_parameters = None, # simulation_parameters = _simulation_parameters) - _parallelism = Parallelism() - _parallelism.activate_client(library = None, - num_processes = 1, - timeout = 1800, - processor = 'cpu') - network = BuildProposalNetwork(parallelism = _parallelism) - network.setup_engines(ligand_input = os.path.join(os.getcwd(), 'test.smi'), - ligand_indices = [0,1], - receptor_filename = None, - graph_connectivity = 'fully_connected', - proposal_parameters = {'phases': ['vacuum', 'solvent']}, - simulation_parameters = _simulation_parameters) - network.create_network() - print(vars(network)) - return network + + from perses.tests.utils import enter_temp_directory + with enter_temp_directory() as tmpdirname: + print(f'Running example in temporary directory: {tmpdirname}') + + _parallelism = Parallelism() + _parallelism.activate_client(library = None, + num_processes = 1, + timeout = 1800, + processor = 'cpu') + network = BuildProposalNetwork(parallelism = _parallelism) + from pkg_resources import resource_filename + smiles_filename = resource_filename("perses", os.path.join("data", "test.smi")) + network.setup_engines(ligand_input = smiles_filename, + ligand_indices = [0,1], + receptor_filename = None, + graph_connectivity = 'fully_connected', + proposal_parameters = {'phases': ['vacuum', 'solvent']}, + simulation_parameters = _simulation_parameters) + network.create_network() + print(vars(network)) if __name__ == "__main__": test_BuildProposalNetwork() diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 59009a192..b4b1e0441 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -1352,35 +1352,32 @@ def _get_capped_amino_acid(amino_acid='ALA'): saveamberparm system {amino_acid}.prmtop {amino_acid}.inpcrd quit """.format(amino_acid=amino_acid) - cwd = os.getcwd() - temp_dir = tempfile.mkdtemp() - os.chdir(temp_dir) - tleap_file = open('tleap_commands', 'w') - tleap_file.writelines(tleapstr) - tleap_file.close() - tleap_cmd_str = "tleap -f %s " % tleap_file.name - - #call tleap, log output to logger - output = getoutput(tleap_cmd_str) - logging.debug(output) - - prmtop = app.AmberPrmtopFile("{amino_acid}.prmtop".format(amino_acid=amino_acid)) - inpcrd = app.AmberInpcrdFile("{amino_acid}.inpcrd".format(amino_acid=amino_acid)) - topology = prmtop.topology - positions = inpcrd.positions - - debug = False - if debug: - system = prmtop.createSystem() - integrator = openmm.VerletIntegrator(1) - context = openmm.Context(system, integrator) - context.setPositions(positions) - openmm.LocalEnergyMinimizer.minimize(context) - state = context.getState(getEnergy=True) - print("%s energy: %s" % (amino_acid, str(state.getPotentialEnergy()))) + from perses.tests.utils import enter_temp_directory + with enter_temp_directory() as tmpdirname: + tleap_file = open('tleap_commands', 'w') + tleap_file.writelines(tleapstr) + tleap_file.close() + tleap_cmd_str = "tleap -f %s " % tleap_file.name + + #call tleap, log output to logger + output = getoutput(tleap_cmd_str) + logging.debug(output) + + prmtop = app.AmberPrmtopFile("{amino_acid}.prmtop".format(amino_acid=amino_acid)) + inpcrd = app.AmberInpcrdFile("{amino_acid}.inpcrd".format(amino_acid=amino_acid)) + topology = prmtop.topology + positions = inpcrd.positions + + debug = False + if debug: + system = prmtop.createSystem() + integrator = openmm.VerletIntegrator(1) + context = openmm.Context(system, integrator) + context.setPositions(positions) + openmm.LocalEnergyMinimizer.minimize(context) + state = context.getState(getEnergy=True) + print("%s energy: %s" % (amino_acid, str(state.getPotentialEnergy()))) - os.chdir(cwd) - shutil.rmtree(temp_dir) return topology, positions def _tleap_all(): diff --git a/perses/tests/test_smc.py b/perses/tests/test_smc.py index a12e56a4b..7d7a1d6fd 100644 --- a/perses/tests/test_smc.py +++ b/perses/tests/test_smc.py @@ -48,7 +48,9 @@ def sMC_setup(): """ function to setup local sMC """ - fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", + from pkg_resources import resource_filename + smiles_filename = resource_filename("perses", os.path.join("data", "test.smi")) + fe_setup = RelativeFEPSetup(ligand_input = smiles_filename, old_ligand_index = 0, new_ligand_index = 1, forcefield_files = ['gaff.xml'], @@ -229,7 +231,9 @@ def test_create_endstates(): """ test the creation of unsampled endstates """ - fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", + from pkg_resources import resource_filename + smiles_filename = resource_filename("perses", os.path.join("data", "test.smi")) + fe_setup = RelativeFEPSetup(ligand_input = smiles_filename, old_ligand_index = 0, new_ligand_index = 1, forcefield_files = ['gaff.xml'], diff --git a/perses/tests/utils.py b/perses/tests/utils.py index 3253d714b..b95b4aa75 100644 --- a/perses/tests/utils.py +++ b/perses/tests/utils.py @@ -27,6 +27,7 @@ from commands import getstatusoutput from openmmtools.constants import kB from openmmtools import alchemy, states +import contextlib ################################################################################ # CONSTANTS @@ -41,6 +42,19 @@ # UTILITIES ################################################################################] +@contextlib.contextmanager +def enter_temp_directory(): + """Create and enter a temporary directory; used as context manager.""" + import tempfile + temp_dir = tempfile.mkdtemp() + import os + cwd = os.getcwd() + os.chdir(temp_dir) + yield temp_dir + os.chdir(cwd) + import shutil + shutil.rmtree(temp_dir) + # TODO: Move some of these utility routines to openmoltools. class NaNException(Exception):