Update experiments branch to use new SystemGenerator

.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

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,

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

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)

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()

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():