Merge pull request #472 from choderalab/use-nonbonded-offset

Use new 'offset' OpenMM feature for nonbonded electrostatics

devtools/conda-recipe/meta.yaml

@@ -18,7 +18,8 @@ requirements:
     - numexpr
     - autograd
     - pymbar
-    - openmm >=7.2.0
+    - cuda92
+    - openmm ==7.3.0
     - parmed
     - openmoltools
     - alchemy >=1.2.3
@@ -43,7 +44,8 @@ requirements:
     - numexpr
     - autograd
     - pymbar
-    - openmm >=7.2.0
+    - cuda92
+    - openmm ==7.3.0
     - parmed
     - openmoltools
     - alchemy >=1.2.3

examples/cdk2-example/cdk2_setup.yaml

@@ -13,8 +13,8 @@ new_ligand_index: 15
 #be provided with a full path
 forcefield_files:
     - gaff.xml
-    - amber99sbildn.xml
-    - tip3p.xml
+    - amber14/protein.ff14SB.xml
+    - amber14/tip3p.xml
 
 #the temperature and pressure of the simulation, as well as how much solvent paddding to add
 #units:
@@ -31,7 +31,9 @@ save_setup_pickle_as: fesetup_hbonds.pkl
 
 #the type of free energy calculation.
 #currently, this could be either nonequilibrium or sams
-fe_type: nonequilibrium
+fe_type: sams
+
+n_states: 100
 
 #the forward switching functions. The reverse ones will be computed from this. Only valid for nonequilibrium switching
 forward_functions:
@@ -42,17 +44,17 @@ forward_functions:
     lambda_torsions: lambda
 
 #the number of equilibration iterations:
-n_equilibration_iterations: 100
+n_equilibration_iterations: 10
 
 #The number of equilibrium steps to take between nonequilibrium switching events (only valid for nonequiibrium switching)
-n_equilibrium_steps_per_iteration: 100
+n_equilibrium_steps_per_iteration: 1000
 
 #The length of the ncmc protocol (only valid for nonequilibrium switching)
 n_steps_ncmc_protocol: 50
 
 #The number of NCMC steps per move application. This controls the output frequency
 #1 step/move application means writing out the work at every step. (only valid for nonequilibrium switching)
-n_steps_per_move_application: 10
+n_steps_per_move_application: 2500
 
 #where to put the trajectories
 trajectory_directory: cdk2_neq_hbonds
@@ -66,6 +68,7 @@ atom_selection: not water
 #which phases do we want to run (solvent, complex, or both solvent and complex in the list are acceptable):
 phases:
     - solvent
+    - complex
 
 #timestep in fs
 timestep: 1.0
@@ -83,5 +86,5 @@ measure_shadow_work: True
 scheduler_address: null
 
 #how many iterations to run (n_cycles*n_iterations_per_cycle) (only valid for nonequilibrium switching)
-n_cycles: 5
+n_cycles: 10000
 n_iterations_per_cycle: 1

perses/dispersed/relative_setup.py

@@ -187,7 +187,8 @@ def __init__(self, ligand_file, old_ligand_index, new_ligand_index, forcefield_f
                 barostat = None
             self._system_generator = SystemGenerator(forcefield_files, barostat=barostat,
                                                      forcefield_kwargs={'nonbondedMethod': self._nonbonded_method,
-                                                                        'constraints': app.HBonds})
+                                                                        'constraints': app.HBonds,
+                                                                        'hydrogenMass': 4 * unit.amus})
         else:
             self._system_generator = SystemGenerator(forcefield_files, forcefield_kwargs={'constraints': app.HBonds})
 
@@ -878,11 +879,11 @@ def __init__(self, *args, hybrid_factory=None, **kwargs):
         super(HybridSAMSSampler, self).__init__(*args, hybrid_factory=hybrid_factory, **kwargs)
         self._factory = hybrid_factory
 
-    def setup(self, n_states, temperature, storage_file, checkpoint_interval):
+    def setup(self, n_states, temperature, storage_file):
         hybrid_system = self._factory.hybrid_system
         initial_hybrid_positions = self._factory.hybrid_positions
         lambda_zero_alchemical_state = alchemy.AlchemicalState.from_system(hybrid_system)
-        lambda_zero_alchemical_state.set_alchemical_parameters(1.0)
+        #lambda_zero_alchemical_state.set_alchemical_parameters(1.0)
 
         thermostate = states.ThermodynamicState(hybrid_system, temperature=temperature)
         compound_thermodynamic_state = states.CompoundThermodynamicState(thermostate, composable_states=[lambda_zero_alchemical_state])
@@ -1050,7 +1051,7 @@ def run_setup(setup_options):
             ne_fep[phase] = NonequilibriumSwitchingFEP(top_prop['%s_topology_proposal' % phase],
                                                        top_prop['%s_old_positions' % phase],
                                                        top_prop['%s_new_positions' % phase],
-                                                       forward_functions=forward_functions,
+                                                       forward_functions=forward_functions, 
                                                        n_equil_steps=n_equilibrium_steps_per_iteration,
                                                        ncmc_nsteps=n_steps_ncmc_protocol,
                                                        nsteps_per_iteration=n_steps_per_move_application,
@@ -1074,24 +1075,25 @@ def run_setup(setup_options):
         for phase in phases:
             htf[phase] = HybridTopologyFactory(top_prop['%s_topology_proposal' % phase],
                                                top_prop['%s_old_positions' % phase],
-                                               top_prop['%s_new_positions' % phase])
+                                               top_prop['%s_new_positions' % phase], softcore_method='amber')
 
             if atom_selection:
                 selection_indices = htf[phase].hybrid_topology.select(atom_selection)
             else:
                 selection_indices = None
 
             storage_name = "-".join([trajectory_prefix, '%s.nc' % phase])
-            reporter = MultiStateReporter(storage_name, analysis_particle_indices=selection_indices)
-
-            hss[phase] = HybridSAMSSampler(mcmc_moves=mcmc.LangevinDynamicsMove(timestep=2.0 * unit.femtosecond,
-                                                                                collision_rate=5.0 / unit.picosecond,
-                                                                                n_steps=n_steps_per_move_application,
-                                                                                reassign_velocities=False,
-                                                                                n_restart_attempts=6),
-                                           hybrid_factory=htf[phase])
-            hss[phase].setup(n_states=n_states, temperature=300.0 * unit.kelvin,
-                             storage_file=reporter,
-                             checkpoint_interval=10)
+            reporter = MultiStateReporter(storage_name, analysis_particle_indices=selection_indices,
+                                          checkpoint_interval=10)
+
+            hss[phase] = HybridSAMSSampler(mcmc_moves=mcmc.LangevinSplittingDynamicsMove(timestep=4.0 * unit.femtosecond,
+                                                                                         collision_rate=5.0 / unit.picosecond,
+                                                                                         n_steps=n_steps_per_move_application,
+                                                                                         reassign_velocities=False,
+                                                                                         n_restart_attempts=6,
+                                                                                         splitting="V R R R O R R R V"),
+                                           hybrid_factory=htf[phase], online_analysis_interval=10,
+                                           online_analysis_target_error=0.2, online_analysis_minimum_iterations=10)
+            hss[phase].setup(n_states=n_states, temperature=300.0 * unit.kelvin, storage_file=reporter)
 
         return {'topology_proposals': top_prop, 'hybrid_topology_factories': htf, 'hybrid_sams_samplers': hss}

perses/rjmc/geometry.py

@@ -373,6 +373,7 @@ def _oemol_from_residue(res, verbose=True):
         oemol : openeye.oechem.OEMol
             an oemol representation of the residue with topology indices
         """
+        # TODO: This seems to be broken. Can we fix it?
         from openmoltools.forcefield_generators import generateOEMolFromTopologyResidue
         external_bonds = list(res.external_bonds())
         for bond in external_bonds:

perses/rjmc/topology_proposal.py

@@ -1324,7 +1324,7 @@ def __init__(self, list_of_smiles, system_generator, residue_name='MOL',
 
         super(SmallMoleculeSetProposalEngine, self).__init__(system_generator, proposal_metadata=proposal_metadata, always_change=always_change)
 
-    def propose(self, current_system, current_topology, current_metadata=None):
+    def propose(self, current_system, current_topology, current_mol=None, proposed_mol=None, current_metadata=None):
         """
         Propose the next state, given the current state
 
@@ -1336,14 +1336,22 @@ def propose(self, current_system, current_topology, current_metadata=None):
             the topology of the current state
         current_metadata : dict
             dict containing current smiles as a key
+        current_mol : OEMol, optional, default=None
+            If specified, use this OEMol instead of converting from topology
+        proposed_mol : OEMol, optional, default=None
+            If specified, use this OEMol instead of converting from topology
 
         Returns
         -------
         proposal : TopologyProposal object
            topology proposal object
         """
         # Determine SMILES string for current small molecule
-        current_mol_smiles, current_mol = self._topology_to_smiles(current_topology)
+        if current_mol is None:
+            current_mol_smiles, current_mol = self._topology_to_smiles(current_topology)
+        else:
+            # TODO: Make sure we're using canonical mol to smiles conversion
+            current_mol_smiles = oechem.OEMolToSmiles(current_mol)
 
         # Remove the small molecule from the current Topology object
         current_receptor_topology = self._remove_small_molecule(current_topology)
@@ -1355,7 +1363,12 @@ def propose(self, current_system, current_topology, current_metadata=None):
         old_alchemical_atoms = range(old_mol_start_index, len_old_mol)
 
         # Select the next molecule SMILES given proposal probabilities
-        proposed_mol_smiles, proposed_mol, logp_proposal = self._propose_molecule(current_system, current_topology, current_mol_smiles)
+        if proposed_mol is None:
+            proposed_mol_smiles, proposed_mol, logp_proposal = self._propose_molecule(current_system, current_topology, current_mol_smiles)
+        else:
+            # TODO: Make sure we're using canonical mol to smiles conversion
+            proposed_mol_smiles = oechem.OEMolToSmiles(current_mol)
+            logp_proposal = 0.0
 
         # Build the new Topology object, including the proposed molecule
         new_topology = self._build_new_topology(current_receptor_topology, proposed_mol)

perses/tests/test_hybrid_builder.py

@@ -2,6 +2,12 @@
 from simtk import unit, openmm
 import numpy as np
 import os
+
+try:
+    from StringIO import StringIO
+except ImportError:
+    from io import StringIO
+
 from perses.annihilation.new_relative import HybridTopologyFactory
 from perses.rjmc.geometry import FFAllAngleGeometryEngine
 from perses.rjmc.topology_proposal import SmallMoleculeSetProposalEngine, SystemGenerator, TopologyProposal
@@ -123,18 +129,18 @@
 
 forcefield = app.ForceField('amber99sbildn.xml')
 
-def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benzene"):
+def generate_vacuum_topology_proposal(current_mol_name="benzene", proposed_mol_name="toluene"):
     """
     Generate a test vacuum topology proposal, current positions, and new positions triplet
     from two IUPAC molecule names.
 
     Parameters
     ----------
-    mol_name : str, optional
+    current_mol_name : str, optional
         name of the first molecule
-    ref_mol_name : str, optional
+    proposed_mol_name : str, optional
         name of the second molecule
-    
+
     Returns
     -------
     topology_proposal : perses.rjmc.topology_proposal
@@ -148,11 +154,11 @@ def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benz
 
     from perses.tests.utils import createOEMolFromIUPAC, createSystemFromIUPAC, get_data_filename
 
-    m, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(mol_name)
-    refmol = createOEMolFromIUPAC(ref_mol_name)
+    current_mol, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(current_mol_name)
+    proposed_mol = createOEMolFromIUPAC(proposed_mol_name)
 
-    initial_smiles = oechem.OEMolToSmiles(m)
-    final_smiles = oechem.OEMolToSmiles(refmol)
+    initial_smiles = oechem.OEMolToSmiles(current_mol)
+    final_smiles = oechem.OEMolToSmiles(proposed_mol)
 
     gaff_xml_filename = get_data_filename("data/gaff.xml")
     forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml')
@@ -161,31 +167,31 @@ def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benz
     solvated_system = forcefield.createSystem(top_old, removeCMMotion=False)
 
     gaff_filename = get_data_filename('data/gaff.xml')
-    system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml'])
+    system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml'], forcefield_kwargs={'removeCMMotion': False, 'nonbondedMethod': app.NoCutoff})
     geometry_engine = FFAllAngleGeometryEngine()
     proposal_engine = SmallMoleculeSetProposalEngine(
-        [initial_smiles, final_smiles], system_generator, residue_name=mol_name)
+        [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name)
 
     #generate topology proposal
-    topology_proposal = proposal_engine.propose(solvated_system, top_old)
+    topology_proposal = proposal_engine.propose(solvated_system, top_old, current_mol=current_mol, proposed_mol=proposed_mol)
 
     #generate new positions with geometry engine
     new_positions, _ = geometry_engine.propose(topology_proposal, pos_old, beta)
 
     return topology_proposal, pos_old, new_positions
 
-def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name="benzene"):
+def generate_solvated_hybrid_test_topology(current_mol_name="naphthalene", proposed_mol_name="benzene"):
     """
     Generate a test solvated topology proposal, current positions, and new positions triplet
     from two IUPAC molecule names.
 
     Parameters
     ----------
-    mol_name : str, optional
+    current_mol_name : str, optional
         name of the first molecule
-    ref_mol_name : str, optional
+    proposed_mol_name : str, optional
         name of the second molecule
-    
+
     Returns
     -------
     topology_proposal : perses.rjmc.topology_proposal
@@ -200,11 +206,11 @@ def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name=
 
     from perses.tests.utils import createOEMolFromIUPAC, createSystemFromIUPAC, get_data_filename
 
-    m, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(mol_name)
-    refmol = createOEMolFromIUPAC(ref_mol_name)
+    current_mol, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(current_mol_name)
+    proposed_mol = createOEMolFromIUPAC(proposed_mol_name)
 
-    initial_smiles = oechem.OEMolToSmiles(m)
-    final_smiles = oechem.OEMolToSmiles(refmol)
+    initial_smiles = oechem.OEMolToSmiles(current_mol)
+    final_smiles = oechem.OEMolToSmiles(proposed_mol)
 
     gaff_xml_filename = get_data_filename("data/gaff.xml")
     forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml')
@@ -219,13 +225,12 @@ def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name=
 
     solvated_system.addForce(barostat)
 
-
     gaff_filename = get_data_filename('data/gaff.xml')
 
     system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml'], barostat=barostat, forcefield_kwargs={'removeCMMotion': False, 'nonbondedMethod': app.PME})
     geometry_engine = FFAllAngleGeometryEngine()
     proposal_engine = SmallMoleculeSetProposalEngine(
-        [initial_smiles, final_smiles], system_generator, residue_name=mol_name)
+        [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name)
 
     #generate topology proposal
     topology_proposal = proposal_engine.propose(solvated_system, solvated_topology)
@@ -306,7 +311,7 @@ def check_result(results, threshold=3.0, neffmin=10):
 
 def test_simple_overlap():
     """Test that the variance of the endpoint->nonalchemical perturbation is sufficiently small for pentane->butane in vacuum"""
-    topology_proposal, current_positions, new_positions = generate_vacuum_topology_proposal()
+    topology_proposal, current_positions, new_positions = generate_vacuum_topology_proposal(current_mol_name='imatinib', proposed_mol_name='nilotinib')
     results = run_hybrid_endpoint_overlap(topology_proposal, current_positions, new_positions)
 
     for idx, lambda_result in enumerate(results):
@@ -320,7 +325,23 @@ def test_simple_overlap():
 @skipIf(istravis, "Skip expensive test on travis")
 def test_difficult_overlap():
     """Test that the variance of the endpoint->nonalchemical perturbation is sufficiently small for imatinib->nilotinib in solvent"""
-    topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(mol_name='imatinib', ref_mol_name='nilotinib')
+    name1 = 'imatinib'
+    name2 = 'nilotinib'
+
+    print(name1, name2)
+    topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(current_mol_name=name1, proposed_mol_name=name2)
+    results = run_hybrid_endpoint_overlap(topology_proposal, solvated_positions, new_positions)
+
+    for idx, lambda_result in enumerate(results):
+        try:
+            check_result(lambda_result)
+        except Exception as e:
+            message = "solvated imatinib->nilotinib failed at lambda %d \n" % idx
+            message += str(e)
+            raise Exception(message)
+
+    print(name2, name1)
+    topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(current_mol_name=name2, proposed_mol_name=name1)
     results = run_hybrid_endpoint_overlap(topology_proposal, solvated_positions, new_positions)
 
     for idx, lambda_result in enumerate(results):

scripts/setup_relative_calculation.py

@@ -32,8 +32,9 @@
 
     setup_dict = relative_setup.run_setup(setup_options)
     print("setup complete")
+    print('setup_dict keys: {}'.format(setup_dict.keys()))
 
-    n_equilibration_iterations = setup_dict['n_equilibration_iterations']
+    n_equilibration_iterations = setup_options['n_equilibration_iterations']
 
     trajectory_prefix = setup_options['trajectory_prefix']
     #write out topology proposals
@@ -97,13 +98,17 @@
                 setup_dict['hybrid_topology_factories'])
 
         hss = setup_dict['hybrid_sams_samplers']
+        logZ = dict()
         free_energies = dict()
         for phase in phases:
             hss_run = hss[phase]
             hss_run.minimize()
             hss_run.equilibrate(n_equilibration_iterations)
-            hss_run.extend(1000)
-            free_energies[phase] = hss_run._logZ[-1] - hss_run._logZ[0]
-            print("Finished phase %s with dG estimated as %.4f kT" % (phase, free_energies[phase]))
+            hss_run.extend(setup_options['n_cycles'])
+            logZ[phase] = hss_run._logZ[-1] - hss_run._logZ[0]
+            free_energies[phase] = hss_run._last_mbar_f_k[-1] - hss_run._last_mbar_f_k[0]
+            print("Finished phase %s with dG estimated as %.4f +/- %.4f kT" % (
+            phase, free_energies[phase], hss_run._last_err_free_energy))
+            print("Finished phase %s with logZ dG estimated as %.4f kT" % (phase, logZ[phase]))
 
         print("Total ddG is estimated as %.4f kT" % (free_energies['complex'] - free_energies['solvent']))