From 7f25bfecad622109531c95be66113da02b4d380a Mon Sep 17 00:00:00 2001 From: John Chodera Date: Sun, 5 Jan 2020 19:53:14 -0800 Subject: [PATCH 01/21] Start of minimal changes necessary for relative free energy calculations to use new SystemGenerator --- .travis.yml | 5 +- devtools/conda-recipe/meta.yaml | 1 + perses/app/relative_setup.py | 56 +++++++----- perses/app/setup_relative_calculation.py | 16 +++- perses/data/cdk2-example/cdk2_setup.yaml | 89 ------------------- perses/data/cdk2-example/cdk2_setup_neq.yaml | 67 ++++++++++++++ .../data/cdk2-example/cdk2_setup_repex.yaml | 73 +++++++++++++++ perses/tests/test_relative_setup.py | 40 +++++++-- 8 files changed, 231 insertions(+), 116 deletions(-) delete mode 100644 perses/data/cdk2-example/cdk2_setup.yaml create mode 100644 perses/data/cdk2-example/cdk2_setup_neq.yaml create mode 100644 perses/data/cdk2-example/cdk2_setup_repex.yaml diff --git a/.travis.yml b/.travis.yml index 9e23edb35..60b5989fc 100644 --- a/.travis.yml +++ b/.travis.yml @@ -58,6 +58,10 @@ script: - conda install --yes --use-local ${PACKAGENAME}-dev # Install testing dependencies - conda install --yes --quiet nose nose-timer + # 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 @@ -69,7 +73,6 @@ script: env: global: - ORGNAME="omnia" - - OE_LICENSE="$HOME/oe_license.txt" - PACKAGENAME="perses" # Location of decrypted OpenEye license file - OE_LICENSE="$HOME/oe_license.txt" diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 8fd229d52..2d29a818e 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -40,6 +40,7 @@ requirements: - progressbar2 - parmed - tqdm + - openmm-forcefields about: home: https://github.com/choderalab/perses diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 52186d7cb..019af54cb 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -5,7 +5,7 @@ from perses.utils.data import load_smi from perses.annihilation.relative import HybridTopologyFactory from perses.annihilation.lambda_protocol import RelativeAlchemicalState, LambdaProtocol -from perses.rjmc.topology_proposal import TopologyProposal, TwoMoleculeSetProposalEngine, SystemGenerator,SmallMoleculeSetProposalEngine +from perses.rjmc.topology_proposal import TopologyProposal, TwoMoleculeSetProposalEngine, SmallMoleculeSetProposalEngine from perses.rjmc.geometry import FFAllAngleGeometryEngine from openmmtools.states import ThermodynamicState, CompoundThermodynamicState, SamplerState @@ -49,7 +49,8 @@ class RelativeFEPSetup(object): def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_files, phases, protein_pdb_filename=None,receptor_mol2_filename=None, pressure=1.0 * unit.atmosphere, temperature=300.0 * unit.kelvin, solvent_padding=9.0 * unit.angstroms, atom_map=None, - hmass=4*unit.amus, neglect_angles=False, map_strength='default', anneal_14s = False): + hmass=4*unit.amus, neglect_angles=False, map_strength='default', anneal_14s = False, + small_molecule_forcefield='gaff-2.11', small_molecule_parameters_cache=None): """ Initialize a NonequilibriumFEPSetup object @@ -78,6 +79,12 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ Whether to anneal 1,4 interactions over the protocol; if True, then geometry_engine takes the argument use_14_nonbondeds = False; if False, then geometry_engine takes the argument use_14_nonbondeds = True; + small_molecule_forcefield : str, optional, default='gaff-2.11' + Small molecule force field name. + Anything supported by SystemGenerator is supported, but we recommend one of + ['gaff-1.81', 'gaff-2.11', 'smirnoff99Frosst-1.1.0', 'openff-1.0.0'] + small_molecule_parameters_cache : str, optional, default=None + If specified, this filename will be used for a small molecule parameter cache by the SystemGenerator. """ self._pressure = pressure self._temperature = temperature @@ -198,20 +205,31 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ self._nonbonded_method = app.NoCutoff _logger.info(f"Detected vacuum phase: setting noCutoff nonbonded method.") + # Select barostat + # TODO: Does this need to omit this barostat for the vacuum phase? + from simtk.openmm.app import NoCutoff, CutoffNonPeriodic + NONPERIODIC_NONBONDED_METHODS = [NoCutoff, CutoffNonPeriodic] if pressure is not None: - if self._nonbonded_method == app.PME: + if self._nonbonded_method not in NONPERIODIC_NONBONDED_METHODS: barostat = openmm.MonteCarloBarostat(self._pressure, self._temperature, self._barostat_period) - _logger.info(f"set MonteCarloBarostat.") + _logger.info(f"set MonteCarloBarostat because pressure was specified as {pressure} atmospheres") else: barostat = None - _logger.info(f"omitted MonteCarloBarostat.") - self._system_generator = SystemGenerator(forcefield_files, barostat=barostat, - forcefield_kwargs={'removeCMMotion': False, 'ewaldErrorTolerance': self._pme_tol, 'nonbondedMethod': self._nonbonded_method,'constraints' : app.HBonds, 'hydrogenMass' : self._hmass}) + _logger.info(f"omitted MonteCarloBarostat because pressure was specified but system was not periodic") else: - self._system_generator = SystemGenerator(forcefield_files, forcefield_kwargs={'removeCMMotion': False, 'ewaldErrorTolerance': self._pme_tol, 'nonbondedMethod': self._nonbonded_method,'constraints' : app.HBonds, 'hydrogenMass' : self._hmass}) + barostat = None + _logger.info(f"omitted MonteCarloBarostat because pressure was not specified") - _logger.info("successfully called TopologyProposal.SystemGenerator to create ligand systems.") - self._system_generator._forcefield.loadFile(StringIO(ffxml)) + # Create openforcefield Molecule objects for old and new molecules + from openforcefield.topology import Molecule + molecules = [ Molecule.from_openeye(oemol) for oemol in [self._ligand_oemol_old, self._ligand_oemol_new] ] + + # Create SystemGenerator + from openmmforcefields.generators import SystemGenerator + forcefield_kwargs = {'removeCMMotion': False, 'ewaldErrorTolerance': self._pme_tol, 'nonbondedMethod': self._nonbonded_method,'constraints' : app.HBonds, 'hydrogenMass' : self._hmass} + self._system_generator = SystemGenerator(forcefields=forcefield_files, barostat=barostat, forcefield_kwargs=forcefield_kwargs, + small_molecule_forcefield=small_molecule_forcefield, molecules=molecules, cache=cache) + _logger.info("successfully created SystemGenerator to create ligand systems") _logger.info(f"executing SmallMoleculeSetProposalEngine...") self._proposal_engine = SmallMoleculeSetProposalEngine([self._ligand_smiles_old, self._ligand_smiles_new], self._system_generator,map_strength=self._map_strength, residue_name='MOL') @@ -311,10 +329,8 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ # need to change nonbonded cutoff and remove barostat for vacuum leg _logger.info(f"assgning noCutoff to nonbonded_method") self._nonbonded_method = app.NoCutoff - _logger.info(f"calling TopologyProposal.SystemGenerator to create ligand systems.") - self._system_generator = SystemGenerator(forcefield_files, forcefield_kwargs={'removeCMMotion': False, 'ewaldErrorTolerance': self._pme_tol, - 'nonbondedMethod': self._nonbonded_method,'constraints' : app.HBonds}) - self._system_generator._forcefield.loadFile(StringIO(ffxml)) + _logger.info(f"calling SystemGenerator to create ligand systems.") + if self._proposal_phase is None: _logger.info('No complex or solvent leg, so performing topology proposal for vacuum leg') self._vacuum_topology_old, self._vacuum_positions_old, self._vacuum_system_old = self._solvate_system(self._ligand_topology_old, @@ -448,7 +464,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()} @@ -507,8 +523,8 @@ def _generate_vacuum_topologies(self, topology_proposal, old_positions): new_ligand_topology = new_ligand_topology.to_openmm() # create the new ligand system: - old_ligand_system = self._system_generator.build_system(old_ligand_topology) - new_ligand_system = self._system_generator.build_system(new_ligand_topology) + old_ligand_system = self._system_generator.create_system(old_ligand_topology) + new_ligand_system = self._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()} @@ -551,10 +567,10 @@ def _solvate_system(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, model=model, padding=self._padding, ionicStrength=0.15*unit.molar) + modeller.addSolvent(self._system_generator.forcefield, model=model, padding=self._padding, ionicStrength=0.15*unit.molar) else: _logger.info(f"\tSkipping solvation of vacuum perturbation") solvated_topology = modeller.getTopology() @@ -563,7 +579,7 @@ def _solvate_system(self, topology, positions, model='tip3p',vacuum=False): # canonicalize the solvated positions: turn tuples into np.array 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) _logger.info(f"\tparameterizing...") - solvated_system = self._system_generator.build_system(solvated_topology) + solvated_system = self._system_generator.create_system(solvated_topology) _logger.info(f"\tSystem parameterized") return solvated_topology, solvated_positions, solvated_system diff --git a/perses/app/setup_relative_calculation.py b/perses/app/setup_relative_calculation.py index 8848592e8..257199f1e 100644 --- a/perses/app/setup_relative_calculation.py +++ b/perses/app/setup_relative_calculation.py @@ -53,6 +53,19 @@ def getSetupOptions(filename): if 'protocol-type' not in setup_options: setup_options['protocol-type'] = 'default' + if 'small_molecule_forcefield' not in setup_options: + setup_options['small_molecule_forcefield'] = 'gaff-2.11' + + if 'small_molecule_parameters_cache' not in setup_options: + setup_options['small_molecule_parameters_cache'] = None + + # Not sure why these are needed + # TODO: Revisit these? + if 'neglect_angles' not in setup_options: + setup_options['neglect_angles'] = False + if 'anneal_1,4s' not in setup_options: + setup_options['anneal_1,4s'] = False + if 'run_type' not in setup_options: _logger.info(f"\t\t\trun_type is not specified; default to None") setup_options['run_type'] = None @@ -318,7 +331,8 @@ def run_setup(setup_options): protein_pdb_filename=protein_pdb_filename, receptor_mol2_filename=receptor_mol2, pressure=pressure, temperature=temperature, solvent_padding=solvent_padding_angstroms, - atom_map=atom_map, neglect_angles = setup_options['neglect_angles'], anneal_14s = setup_options['anneal_1,4s']) + atom_map=atom_map, neglect_angles = setup_options['neglect_angles'], anneal_14s = setup_options['anneal_1,4s'], + small_molecule_forcefield=setup_options['small_molecule_forcefield'], small_molecule_parameters_cache=setup_options['small_molecule_parameters_cache']) _logger.info(f"\n\n\n") _logger.info(f"\twriting pickle output...") diff --git a/perses/data/cdk2-example/cdk2_setup.yaml b/perses/data/cdk2-example/cdk2_setup.yaml deleted file mode 100644 index c4108da84..000000000 --- a/perses/data/cdk2-example/cdk2_setup.yaml +++ /dev/null @@ -1,89 +0,0 @@ -#provide the full path of the protein PDB file -protein_pdb: CDK2_fixed_nohet.pdb - -#provide the path to the ligand file with coordinates -ligand_file: CDK2_ligands.sdf - -#The ligand file contains multiple ligands. Choose the indices of the ligands -#between which we should compute a relative free energy -old_ligand_index: 14 -new_ligand_index: 15 - -#Provide the list of forcefield files. Non-standard files should -#be provided with a full path -forcefield_files: - - gaff.xml - - amber99sbildn.xml - - tip3p.xml - -#the temperature and pressure of the simulation, as well as how much solvent paddding to add -#units: -#pressure: atm -#temperature: Kelvin -#padding: angstroms -pressure: 1.0 -temperature: 300.0 -solvent_padding: 9.0 - - -#The name of the pickle file where we will save the setup object -save_setup_pickle_as: fesetup_hbonds.pkl - -#the type of free energy calculation. -#currently, this could be either nonequilibrium or sams -fe_type: nonequilibrium - -#the forward switching functions. The reverse ones will be computed from this. Only valid for nonequilibrium switching -forward_functions: - lambda_sterics: lambda - lambda_electrostatics: lambda - lambda_bonds: lambda - lambda_angles: lambda - lambda_torsions: lambda - -#the number of equilibration iterations: -n_equilibration_iterations: 100 - -#The number of equilibrium steps to take between nonequilibrium switching events (only valid for nonequiibrium switching) -n_equilibrium_steps_per_iteration: 100 - -#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 - -#where to put the trajectories -trajectory_directory: cdk2_neq_hbonds - -#how to prefix the trajectory files (project-specific name) -trajectory_prefix: cdk2 - -#which atoms to save (MDTraj selection syntax) -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 - -#timestep in fs -timestep: 1.0 - -#splitting string for equilibrium integration (only valid for nonequilibrium switching): -eq_splitting: V R O R V - -#splitting string for nonequilibrium switching (only valid for nonequilibrium switching): -neq_splitting: V R O H R V - -#whether to measure shadow work (only valid for nonequilibrium switching): -measure_shadow_work: True - -#The location of the schduler. If it's null, a localhost scheduler is created (only valid for nonequilibrium switching) -scheduler_address: null - -#how many iterations to run (n_cycles*n_iterations_per_cycle) (only valid for nonequilibrium switching) -n_cycles: 5 -n_iterations_per_cycle: 1 - -neglect_angles: False diff --git a/perses/data/cdk2-example/cdk2_setup_neq.yaml b/perses/data/cdk2-example/cdk2_setup_neq.yaml new file mode 100644 index 000000000..b2754bb70 --- /dev/null +++ b/perses/data/cdk2-example/cdk2_setup_neq.yaml @@ -0,0 +1,67 @@ +# Relative path to protein PDB file (with no missing heavy atoms) +protein_pdb: CDK2_fixed_nohet.pdb + +# Relative path to ligand SDF file +ligand_file: CDK2_ligands.sdf + +# Indices of the ligands in the ligand_file between which we should compute a relative free energy +# (zero-indexed) +old_ligand_index: 14 +new_ligand_index: 15 + +# OpenMM ffxml force field files installed via the openmm-forcefields package +# for biopolymers and solvents. +# Note that small molecule force field files should NOT be included here. +forcefield_files: + - amber/protein.ff14SB.xml # ff14SB protein force field + - amber/tip3p_standard.xml # TIP3P and recommended monovalent ion parameters + - amber/tip3p_HFE_multivalent.xml # for divalent ions + +# Small molecule force field +# We recommnd one of ['gaff-1.81', 'gaff-2.11', 'smirnoff99Frosst-1.1.0', 'openff-1.0.0'] +small_molecule_forcefield: gaff-2.11 + +# Thermodynamic state and system setup information +pressure: 1.0 # atmospheres +temperature: 300.0 # kelvin +solvent_padding: 9.0 # angstroms + +# Pickle file to be generated to save the setup object +save_setup_pickle_as: fesetup_hbonds.pkl + +#number of equilibrium steps per move application +n_steps_per_move_application: 250 + +#The name of the pickle file where we will save the setup object +save_setup_pickle_as: fesetup_hbonds.pkl + +#the type of free energy calculation. +#currently, this could be either nonequilibrium or sams +fe_type: repex + +#checkpoint interval in iterations: +checkpoint_interval: 250 + +#number of iterations +n_cycles: 5000 + +#The number of SAMS states +n_states: 16 + +# Number of equilibration iterations: +n_equilibration_iterations: 0 + +# Destination for trajectories +trajectory_directory: output + +#how to prefix the trajectory files (project-specific name) +trajectory_prefix: out + +#which atoms to save (MDTraj selection syntax) +atom_selection: not water + +#which phases do we want to run. Any combination or number of complex, solvent or vacuum will be accepted +phases: + - complex + - solvent + - vacuum diff --git a/perses/data/cdk2-example/cdk2_setup_repex.yaml b/perses/data/cdk2-example/cdk2_setup_repex.yaml new file mode 100644 index 000000000..b4a8c9120 --- /dev/null +++ b/perses/data/cdk2-example/cdk2_setup_repex.yaml @@ -0,0 +1,73 @@ +# Relative path to protein PDB file (with no missing heavy atoms) +protein_pdb: CDK2_fixed_nohet.pdb + +# Relative path to ligand SDF file +ligand_file: CDK2_ligands.sdf + +# Indices of the ligands in the ligand_file between which we should compute a relative free energy +# (zero-indexed) +old_ligand_index: 14 +new_ligand_index: 15 + +# OpenMM ffxml force field files installed via the openmm-forcefields package +# for biopolymers and solvents. +# Note that small molecule force field files should NOT be included here. +forcefield_files: + - amber/protein.ff14SB.xml # ff14SB protein force field + - amber/tip3p_standard.xml # TIP3P and recommended monovalent ion parameters + - amber/tip3p_HFE_multivalent.xml # for divalent ions + +# Small molecule force field +# We recommnd one of ['gaff-1.81', 'gaff-2.11', 'smirnoff99Frosst-1.1.0', 'openff-1.0.0'] +small_molecule_forcefield: gaff-2.11 + +# Small molecule parameters cache +# If specified, this will cache small molecule parameters so they do not need to be regenerated. +# The cache can contain parameters for many different choices of 'small_molecule_forcefield' without fear of collisions +small_molecule_parameters_cache: cdk2-cache.json + +# Thermodynamic state and system setup information +pressure: 1.0 # atmospheres +temperature: 300.0 # kelvin +solvent_padding: 9.0 # angstroms + +# Pickle file to be generated to save the setup object +save_setup_pickle_as: fesetup_hbonds.pkl + +#number of equilibrium steps per move application +n_steps_per_move_application: 250 + +#The name of the pickle file where we will save the setup object +save_setup_pickle_as: fesetup_hbonds.pkl + +#the type of free energy calculation. +#currently, this could be either nonequilibrium or sams +fe_type: repex + +#checkpoint interval in iterations: +checkpoint_interval: 250 + +#number of iterations +n_cycles: 5000 + +#The number of SAMS states +n_states: 16 + +# Number of equilibration iterations: +n_equilibration_iterations: 0 + +# Destination for trajectories +trajectory_directory: output +overwrite_output: true + +#how to prefix the trajectory files (project-specific name) +trajectory_prefix: out + +#which atoms to save (MDTraj selection syntax) +atom_selection: not water + +#which phases do we want to run. Any combination or number of complex, solvent or vacuum will be accepted +phases: + - complex + - solvent + - vacuum diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index 78b843dfc..4c0446a94 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -70,15 +70,15 @@ def test_run_nonequilibrium_switching_move(): assert integrator.getGlobalVariableByName("lambda") == 1.0 @skipIf(os.environ.get("TRAVIS", None) == 'true', "Skip slow test on TRAVIS.") -def test_run_cdk2_iterations(): +def test_run_cdk2_iterations_neq(): """ - Ensure that we can instantiate and run the cdk2 ligands in vacuum + Ensure that we can instantiate and run a nonequilibrium relative free energy calculation for the cdk2 ligands in vacuum """ setup_directory = resource_filename("perses", "data/cdk2-example") os.chdir(setup_directory) n_iterations = 2 - yaml_filename = "cdk2_setup.yaml" + yaml_filename = "cdk2_setup_neq.yaml" yaml_file = open(yaml_filename, "r") setup_options = yaml.safe_load(yaml_file) yaml_file.close() @@ -114,6 +114,36 @@ def test_run_cdk2_iterations(): lambda_one_npy = np.stack([np.load(filename) for filename in lambda_one_filenames]) assert np.shape(lambda_one_npy) == (n_iterations, n_work_values_per_iteration+1) +@skipIf(os.environ.get("TRAVIS", None) == 'true', "Skip slow test on TRAVIS.") +def test_run_cdk2_iterations_repex(): + """ + Ensure that we can instantiate and run a repex relative free energy calculation the cdk2 ligands in vacuum + """ + setup_directory = resource_filename("perses", "data/cdk2-example") + os.chdir(setup_directory) + n_iterations = 2 + + yaml_filename = "cdk2_setup_repex.yaml" + from perses.app.setup_relative_calculation import getSetupOptions + setup_options = getSetupOptions(yaml_filename) + + if not os.path.exists(setup_options['trajectory_directory']): + os.makedirs(setup_options['trajectory_directory']) + + #setup_options['solvate'] = False + setup_options['scheduler_address'] = None + + #length_of_protocol = setup_options['n_steps_ncmc_protocol'] + #write_interval = setup_options['n_steps_per_move_application'] + + #n_work_values_per_iteration = length_of_protocol // write_interval + + setup_dict = setup_relative_calculation.run_setup(setup_options) + print(setup_dict) # DEBUG + setup_dict['hybrid_samplers']['solvent'].run(n_iterations=n_iterations) + + # TODO: Check output + + if __name__=="__main__": - test_run_cdk2_iterations() - test_run_nonequilibrium_switching_move() + test_run_cdk2_iterations_repex() From 783e203eda332e2368426493b9ecedfee222a6ea Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 7 Jan 2020 20:31:27 -0500 Subject: [PATCH 02/21] Clean up CDK2 repex example. --- perses/app/relative_setup.py | 2 +- perses/app/setup_relative_calculation.py | 5 +- .../data/cdk2-example/cdk2_setup_repex.yaml | 6 +- perses/tests/test_relative_setup.py | 58 +++++++++++-------- 4 files changed, 39 insertions(+), 32 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index a71b34c1d..198696599 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -226,7 +226,7 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ from openmmforcefields.generators import SystemGenerator forcefield_kwargs = {'removeCMMotion': False, 'ewaldErrorTolerance': self._pme_tol, 'nonbondedMethod': self._nonbonded_method,'constraints' : app.HBonds, 'hydrogenMass' : self._hmass} self._system_generator = SystemGenerator(forcefields=forcefield_files, barostat=barostat, forcefield_kwargs=forcefield_kwargs, - small_molecule_forcefield=small_molecule_forcefield, molecules=molecules, cache=cache) + small_molecule_forcefield=small_molecule_forcefield, molecules=molecules, cache=small_molecule_parameters_cache) _logger.info("successfully created SystemGenerator to create ligand systems") _logger.info(f"executing SmallMoleculeSetProposalEngine...") diff --git a/perses/app/setup_relative_calculation.py b/perses/app/setup_relative_calculation.py index b3345e07c..2c303ac48 100644 --- a/perses/app/setup_relative_calculation.py +++ b/perses/app/setup_relative_calculation.py @@ -152,8 +152,9 @@ def getSetupOptions(filename): setup_options['softcore_v2'] = False _logger.info(f"\t'softcore_v2' not specified: default to 'False'") - _logger.info(f"\ttrajectory_directory detected: {trajectory_directory}. making dir...") - assert (os.path.exists(trajectory_directory) == False), 'Output trajectory directory already exists. Refusing to overwrite' + _logger.info(f"\tCreating '{trajectory_directory}'...") + print(os.getcwd()) + assert (not os.path.exists(trajectory_directory)), f'Output trajectory directory "{trajectory_directory}" already exists. Refusing to overwrite' os.makedirs(trajectory_directory) return setup_options diff --git a/perses/data/cdk2-example/cdk2_setup_repex.yaml b/perses/data/cdk2-example/cdk2_setup_repex.yaml index b4a8c9120..3eabf4ef1 100644 --- a/perses/data/cdk2-example/cdk2_setup_repex.yaml +++ b/perses/data/cdk2-example/cdk2_setup_repex.yaml @@ -34,12 +34,9 @@ solvent_padding: 9.0 # angstroms # Pickle file to be generated to save the setup object save_setup_pickle_as: fesetup_hbonds.pkl -#number of equilibrium steps per move application +# Number of equilibrium steps per move application n_steps_per_move_application: 250 -#The name of the pickle file where we will save the setup object -save_setup_pickle_as: fesetup_hbonds.pkl - #the type of free energy calculation. #currently, this could be either nonequilibrium or sams fe_type: repex @@ -58,7 +55,6 @@ n_equilibration_iterations: 0 # Destination for trajectories trajectory_directory: output -overwrite_output: true #how to prefix the trajectory files (project-specific name) trajectory_prefix: out diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index 4c0446a94..ad4968a24 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -119,30 +119,40 @@ def test_run_cdk2_iterations_repex(): """ Ensure that we can instantiate and run a repex relative free energy calculation the cdk2 ligands in vacuum """ - setup_directory = resource_filename("perses", "data/cdk2-example") - os.chdir(setup_directory) - n_iterations = 2 - - yaml_filename = "cdk2_setup_repex.yaml" - from perses.app.setup_relative_calculation import getSetupOptions - setup_options = getSetupOptions(yaml_filename) - - if not os.path.exists(setup_options['trajectory_directory']): - os.makedirs(setup_options['trajectory_directory']) - - #setup_options['solvate'] = False - setup_options['scheduler_address'] = None - - #length_of_protocol = setup_options['n_steps_ncmc_protocol'] - #write_interval = setup_options['n_steps_per_move_application'] - - #n_work_values_per_iteration = length_of_protocol // write_interval - - setup_dict = setup_relative_calculation.run_setup(setup_options) - print(setup_dict) # DEBUG - setup_dict['hybrid_samplers']['solvent'].run(n_iterations=n_iterations) - - # TODO: Check output + # Enter a temporary directory + import tempfile + with tempfile.TemporaryDirectory() as tmpdirname: + # Move to temporary directory + os.chdir(tmpdirname) + print(f'Running example in temporary directory: {tmpdirname}') + + # Setup directory + setup_directory = resource_filename("perses", "data/cdk2-example") + + # Get options + from perses.app.setup_relative_calculation import getSetupOptions + yaml_filename = os.path.join(setup_directory, "cdk2_setup_repex.yaml") + setup_options = getSetupOptions(yaml_filename) + + # Update options + #setup_options['solvate'] = False + #setup_options['n_cycles'] = 2 + setup_options['scheduler_address'] = None + for parameter in ['protein_pdb', 'ligand_file', 'small_molecule_parameters_cache']: + setup_options[parameter] = os.path.join(setup_directory, setup_options[parameter]) + for parameter in ['trajectory_directory', 'trajectory_prefix', 'save_setup_pickle_as']: + setup_options[parameter] = os.path.join(tmpdirname, setup_options[parameter]) + print(setup_options) + + #length_of_protocol = setup_options['n_steps_ncmc_protocol'] + #write_interval = setup_options['n_steps_per_move_application'] + #n_work_values_per_iteration = length_of_protocol // write_interval + + # Run setup + setup_dict = setup_relative_calculation.run_setup(setup_options) + setup_dict['hybrid_samplers']['solvent'].run(n_iterations=n_iterations) + + # TODO: Check output if __name__=="__main__": From adb7d3ee97348202f7b1cacc6d533d23d15c3c96 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 7 Jan 2020 21:15:52 -0500 Subject: [PATCH 03/21] More debugging --- perses/app/relative_setup.py | 20 ++++++-------------- perses/tests/test_relative_setup.py | 11 ++++++++++- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 198696599..7cbfa83cd 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -58,7 +58,7 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ the name of the ligand file (any openeye supported format) this can either be an .sdf or list of .sdf files, or a list of SMILES strings forcefield_files : list of str - The list of ffxml files that contain the forcefields that will be used + The list of ffxml files that contain the forcefields that will be used (excepting for small molecules) phases : list of str The phases to simulate protein_pdb_filename : str, default None @@ -122,9 +122,6 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ mol_list.append(self._ligand_oemol_old) mol_list.append(self._ligand_oemol_new) - ffxml = forcefield_generators.generateForceFieldFromMolecules(mol_list) - _logger.info(f"\tsuccessfully generated ffxml from molecules.") - # forcefield_generators needs to be able to distinguish between the two ligands # while topology_proposal needs them to have the same residue name self._ligand_oemol_old.SetTitle("MOL") @@ -146,9 +143,6 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ self._ligand_positions_old = extractPositionsFromOEMol(self._ligand_oemol_old) _logger.info(f"\tsuccessfully extracted positions from OEMOL.") - ffxml = forcefield_generators.generateForceFieldFromMolecules(mol_list) - _logger.info(f"\tsuccessfully generated ffxml from molecules.") - self._ligand_oemol_old.SetTitle("MOL") self._ligand_oemol_new.SetTitle("MOL") _logger.info(f"\tsetting both molecule oemol titles to 'MOL'.") @@ -185,13 +179,6 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ mol_list.append(self._ligand_oemol_old) mol_list.append(self._ligand_oemol_new) - old_ligand_parameter_set = pm.openmm.OpenMMParameterSet.from_structure(old_ligand) - new_ligand_parameter_set = pm.openmm.OpenMMParameterSet.from_structure(new_ligand) - ffxml = StringIO() - old_ligand_parameter_set.write(ffxml) - new_ligand_parameter_set.write(ffxml) - ffxml = ffxml.getvalue() - self._ligand_md_topology_old = md.Topology.from_openmm(self._ligand_topology_old) self._ligand_md_topology_new = md.Topology.from_openmm(self._ligand_topology_new) _logger.info(f"Created mdtraj topologies for both ligands.") @@ -562,6 +549,11 @@ def _solvate_system(self, topology, positions, model='tip3p',vacuum=False): solvated_system : openmm.System The parameterized system, containing a barostat if one was specified. """ + # DEBUG: Write PDB file being fed into Modeller to check why MOL isn't being matched + from simtk.openmm.app import PDBFile + with open('/Users/choderaj/github/choderalab/perses/modeller.pdb', 'w') as outfile: + PDBFile.writeFile(topology, positions, outfile) + 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) diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index ad4968a24..80f0be48f 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -134,6 +134,16 @@ def test_run_cdk2_iterations_repex(): yaml_filename = os.path.join(setup_directory, "cdk2_setup_repex.yaml") setup_options = getSetupOptions(yaml_filename) + # DEBUG: Print traceback for any UserWarnings + import traceback + import warnings + _old_warn = warnings.warn + def warn(*args, **kwargs): + tb = traceback.extract_stack() + _old_warn(*args, **kwargs) + print("".join(traceback.format_list(tb)[:-1])) + warnings.warn = warn + # Update options #setup_options['solvate'] = False #setup_options['n_cycles'] = 2 @@ -142,7 +152,6 @@ def test_run_cdk2_iterations_repex(): setup_options[parameter] = os.path.join(setup_directory, setup_options[parameter]) for parameter in ['trajectory_directory', 'trajectory_prefix', 'save_setup_pickle_as']: setup_options[parameter] = os.path.join(tmpdirname, setup_options[parameter]) - print(setup_options) #length_of_protocol = setup_options['n_steps_ncmc_protocol'] #write_interval = setup_options['n_steps_per_move_application'] From 2256c3303f545468ca0970add09dfdd4724d60cf Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 7 Jan 2020 21:18:11 -0500 Subject: [PATCH 04/21] Use relative path for debugging --- perses/app/relative_setup.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 7cbfa83cd..9c5206874 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -551,7 +551,9 @@ def _solvate_system(self, topology, positions, model='tip3p',vacuum=False): """ # DEBUG: Write PDB file being fed into Modeller to check why MOL isn't being matched from simtk.openmm.app import PDBFile - with open('/Users/choderaj/github/choderalab/perses/modeller.pdb', 'w') as outfile: + import os + pdb_filename = os.path.join(os.environ['PWD'], 'modeller.pdb') + with open(pdb_filename, 'w') as outfile: PDBFile.writeFile(topology, positions, outfile) modeller = app.Modeller(topology, positions) From 0d7510b22fd330e282b952915133e0deaaf09e15 Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Wed, 8 Jan 2020 17:11:20 -0500 Subject: [PATCH 05/21] generating unique atom names for oemols --- perses/app/relative_setup.py | 6 ++++++ perses/utils/openeye.py | 37 ++++++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 9c5206874..fffb9d72e 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -114,9 +114,11 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ all_old_mol = createSystemFromSMILES(self._ligand_smiles_old, title='MOL') # should be stereospecific self._ligand_oemol_old, self._ligand_system_old, self._ligand_positions_old, self._ligand_topology_old = all_old_mol + self._ligand_oemol_old = generate_unique_atom_names(self._ligand_oemol_old) all_new_mol = createSystemFromSMILES(self._ligand_smiles_new, title='NEW') self._ligand_oemol_new, self._ligand_system_new, self._ligand_positions_new, self._ligand_topology_new = all_new_mol + self._ligand_oemol_new = generate_unique_atom_names(self._ligand_oemol_new) _logger.info(f"\tsuccessfully created old and new systems from smiles") mol_list.append(self._ligand_oemol_old) @@ -136,6 +138,8 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ _logger.info(f"Detected .sdf format. Proceeding...") #TODO: write checkpoints for sdf format self._ligand_oemol_old = createOEMolFromSDF(self._ligand_input, index=self._old_ligand_index) self._ligand_oemol_new = createOEMolFromSDF(self._ligand_input, index=self._new_ligand_index) + self._ligand_oemol_old = generate_unique_atom_names(self._ligand_oemol_old) + self._ligand_oemol_new = generate_unique_atom_names(self._ligand_oemol_new) mol_list.append(self._ligand_oemol_old) mol_list.append(self._ligand_oemol_new) @@ -166,6 +170,7 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ self._ligand_topology_old = old_ligand.topology self._ligand_positions_old = old_ligand.positions self._ligand_oemol_old = createOEMolFromSDF('%s.mol2' % self._ligand_input[0]) + self._ligand_oemol_old = generate_unique_atom_names(self._ligand_oemol_old) self._ligand_smiles_old = oechem.OECreateSmiString(self._ligand_oemol_old, oechem.OESMILESFlag_DEFAULT | oechem.OESMILESFlag_Hydrogens) @@ -173,6 +178,7 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ self._ligand_topology_new = new_ligand.topology self._ligand_positions_new = new_ligand.positions self._ligand_oemol_new = createOEMolFromSDF('%s.mol2' % self._ligand_input[1]) + self._ligand_oemol_new = generate_unique_atom_names(self._ligand_oemol_new) self._ligand_smiles_new = oechem.OECreateSmiString(self._ligand_oemol_new, oechem.OESMILESFlag_DEFAULT | oechem.OESMILESFlag_Hydrogens) diff --git a/perses/utils/openeye.py b/perses/utils/openeye.py index db2083fa1..789aee2ab 100644 --- a/perses/utils/openeye.py +++ b/perses/utils/openeye.py @@ -269,3 +269,40 @@ def createSMILESfromOEMol(molecule): smiles = oechem.OECreateSmiString(molecule, oechem.OESMILESFlag_DEFAULT | oechem.OESMILESFlag_Hydrogens) return smiles + + +def generate_unique_atom_names(molecule): + """ + Check if an oemol has unique atom names, and if not, then assigns them + + Parameters + ---------- + molecule : openeye.oechem.OEMol object + oemol object to check + + Returns + ------- + molecule : openeye.oechem.OEMol object + oemol, either unchanged if atom names are already unique, or newly generated atom names + """ + atom_names = [] + + atom_count = 0 + for atom in molecule.GetAtoms(): + atom_names.append(atom.GetName()) + atom_count += 1 + + if len(set(atom_names)) == atom_count: + # one name per atom therefore unique + return molecule + else: + # generating new atom names + from collections import defaultdict + from simtk.openmm.app.element import Element + element_counts = defaultdict(int) + for atom in molecule.GetAtoms(): + element = Element.getByAtomicNumber(atom.GetAtomicNum()) + element_counts[element._symbol] += 1 + name = element._symbol + str(element_counts[element._symbol]) + atom.SetName(name) + return molecule From 2a3704d6f68807178e4cdd3307084b30856f236a Mon Sep 17 00:00:00 2001 From: John Chodera Date: Wed, 8 Jan 2020 17:27:45 -0500 Subject: [PATCH 06/21] Fix build_system errors --- perses/rjmc/topology_proposal.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/topology_proposal.py b/perses/rjmc/topology_proposal.py index 1b642a477..998bf8c78 100644 --- a/perses/rjmc/topology_proposal.py +++ b/perses/rjmc/topology_proposal.py @@ -1078,7 +1078,11 @@ def propose(self, current_system, current_topology, current_metadata=None): 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 +2443,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 +3138,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) From 8812b50a6f465b0fc26b6e766bce300915cfbb6c Mon Sep 17 00:00:00 2001 From: John Chodera Date: Wed, 8 Jan 2020 17:32:04 -0500 Subject: [PATCH 07/21] Comment out stacktraces for warnings --- perses/tests/test_relative_setup.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/perses/tests/test_relative_setup.py b/perses/tests/test_relative_setup.py index 80f0be48f..c31a8ae98 100644 --- a/perses/tests/test_relative_setup.py +++ b/perses/tests/test_relative_setup.py @@ -135,14 +135,16 @@ def test_run_cdk2_iterations_repex(): setup_options = getSetupOptions(yaml_filename) # DEBUG: Print traceback for any UserWarnings - import traceback - import warnings - _old_warn = warnings.warn - def warn(*args, **kwargs): - tb = traceback.extract_stack() - _old_warn(*args, **kwargs) - print("".join(traceback.format_list(tb)[:-1])) - warnings.warn = warn + show_warning_stacktraces = False + if show_warning_stacktraces: + import traceback + import warnings + _old_warn = warnings.warn + def warn(*args, **kwargs): + tb = traceback.extract_stack() + _old_warn(*args, **kwargs) + print("".join(traceback.format_list(tb)[:-1])) + warnings.warn = warn # Update options #setup_options['solvate'] = False From f10076b4a4ccbf05f1c3787df7b397fcbc332d40 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Wed, 8 Jan 2020 17:35:54 -0500 Subject: [PATCH 08/21] Add JSON cache --- perses/data/cdk2-example/cdk2-cache.json | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 perses/data/cdk2-example/cdk2-cache.json diff --git a/perses/data/cdk2-example/cdk2-cache.json b/perses/data/cdk2-example/cdk2-cache.json new file mode 100644 index 000000000..f8e195031 --- /dev/null +++ b/perses/data/cdk2-example/cdk2-cache.json @@ -0,0 +1,15 @@ +{ + "_default": {}, + "gaff-2.11": { + "1": { + "ffxml": "\n \n 2020-01-08\n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "iupac": "~{N}-(3-bromophenyl)-6-(cyclohexylmethoxy)-9~{H}-purin-2-amine", + "smiles": "[H]c1c(c(c(c(c1[H])Br)[H])N([H])c2nc3c(c(n2)OC([H])([H])C4(C(C(C(C(C4([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])N=C(N3[H])[H])[H]" + }, + "2": { + "ffxml": "\n \n 2020-01-08\n \n \n \n \n \n \n \n \n \n \n \n \n\n", + "iupac": "6-(cyclohexylmethoxy)-~{N}-phenyl-9~{H}-purin-2-amine", + "smiles": "[H]c1c(c(c(c(c1[H])[H])N([H])c2nc3c(c(n2)OC([H])([H])C4(C(C(C(C(C4([H])[H])([H])[H])([H])[H])([H])[H])([H])[H])[H])N=C(N3[H])[H])[H])[H]" + } + } +} \ No newline at end of file From b27bf959c1b83cfa436cfa8e2d56924d7cbfe566 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 10 Jan 2020 12:26:07 -0500 Subject: [PATCH 09/21] Install openmm-forcefields via pip for now --- .travis.yml | 2 ++ devtools/conda-recipe/meta.yaml | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index f76579f16..8f52fdcd8 100644 --- a/.travis.yml +++ b/.travis.yml @@ -56,6 +56,8 @@ 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) + - 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" diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 2d29a818e..39e19d8a6 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -40,7 +40,7 @@ requirements: - progressbar2 - parmed - tqdm - - openmm-forcefields + #- openmm-forcefields # uncomment once openmm-forcefields package is pushed to conda about: home: https://github.com/choderalab/perses From 8b1f2bfa6476f23bf4ce2e274802edadf46ebdf9 Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Sat, 11 Jan 2020 07:24:13 -0500 Subject: [PATCH 10/21] adding logger comments --- perses/utils/openeye.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/perses/utils/openeye.py b/perses/utils/openeye.py index 789aee2ab..6e9bf1f14 100644 --- a/perses/utils/openeye.py +++ b/perses/utils/openeye.py @@ -11,6 +11,11 @@ from openmoltools.openeye import iupac_to_oemol, generate_conformers import simtk.unit as unit import numpy as np +import logging + +logging.basicConfig(level = logging.NOTSET) +_logger = logging.getLogger("utils.openeye") +_logger.setLevel(logging.INFO) def smiles_to_oemol(smiles, title='MOL',max_confs=1): """ @@ -294,11 +299,13 @@ def generate_unique_atom_names(molecule): if len(set(atom_names)) == atom_count: # one name per atom therefore unique + _logger.info(f'molecule {molecule.GetTitle()} has unique atom names already') return molecule else: # generating new atom names from collections import defaultdict from simtk.openmm.app.element import Element + _logger.info(f'molecule {molecule.GetTitle()} does not have unique atom names. Generating now...') element_counts = defaultdict(int) for atom in molecule.GetAtoms(): element = Element.getByAtomicNumber(atom.GetAtomicNum()) From f4cf4355e24290dbd27f6d65f714107f60a7f258 Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Sat, 11 Jan 2020 07:24:29 -0500 Subject: [PATCH 11/21] removing de- and re- protonation of protein --- perses/app/relative_setup.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index f8fea01f1..5ead03d26 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -566,9 +566,10 @@ def _solvate_system(self, topology, positions, model='tip3p',vacuum=False): PDBFile.writeFile(topology, positions, outfile) 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) + # retaining protein protonation from input files + #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) if not vacuum: _logger.info(f"\tpreparing to add solvent") modeller.addSolvent(self._system_generator.forcefield, model=model, padding=self._padding, ionicStrength=0.15*unit.molar) From 26b5bd16296277331352013d2a09e69383da9220 Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Sat, 11 Jan 2020 14:10:54 -0500 Subject: [PATCH 12/21] minor logger changes --- perses/app/setup_relative_calculation.py | 3 +-- perses/rjmc/topology_proposal.py | 1 - 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/perses/app/setup_relative_calculation.py b/perses/app/setup_relative_calculation.py index 047d9deea..128df105e 100644 --- a/perses/app/setup_relative_calculation.py +++ b/perses/app/setup_relative_calculation.py @@ -277,7 +277,7 @@ def run_setup(setup_options): if "timestep" in setup_options: timestep = setup_options['timestep'] * unit.femtoseconds - _logger.info(f"\ttimestep: {timestep}fs.") + _logger.info(f"\ttimestep: {timestep}.") else: timestep = 1.0 * unit.femtoseconds _logger.info(f"\tno timestep detected: setting default as 1.0fs.") @@ -333,7 +333,6 @@ def run_setup(setup_options): temperature=temperature, solvent_padding=solvent_padding_angstroms, atom_map=atom_map, neglect_angles = setup_options['neglect_angles'], anneal_14s = setup_options['anneal_1,4s'], small_molecule_forcefield=setup_options['small_molecule_forcefield'], small_molecule_parameters_cache=setup_options['small_molecule_parameters_cache']) - _logger.info(f"\n\n\n") _logger.info(f"\twriting pickle output...") with open(os.path.join(os.getcwd(), trajectory_directory, setup_pickle_file), 'wb') as f: diff --git a/perses/rjmc/topology_proposal.py b/perses/rjmc/topology_proposal.py index 6b03cd877..30287f3d8 100644 --- a/perses/rjmc/topology_proposal.py +++ b/perses/rjmc/topology_proposal.py @@ -2319,7 +2319,6 @@ def __init__(self, list_of_smiles, system_generator, residue_name='MOL', # Default atom and bond expressions for MCSS if atom_expr is None: _logger.info(f'Setting the atom expression to {map_strength}') - _logger.info(type(map_strength)) if map_strength == 'default': self.atom_expr = DEFAULT_ATOM_EXPRESSION elif map_strength == 'weak': From 3e219d324714eb90be3037dd9aa895c8b00405ab Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Sat, 11 Jan 2020 14:11:26 -0500 Subject: [PATCH 13/21] fixing barostat/cutoff for vacuum phase --- perses/app/relative_setup.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 5ead03d26..5449fec0a 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -194,12 +194,11 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ if 'complex' in phases or 'solvent' in phases: self._nonbonded_method = app.PME _logger.info(f"Detected complex or solvent phases: setting PME nonbonded method.") - elif 'vacuum' in phases: + else: self._nonbonded_method = app.NoCutoff _logger.info(f"Detected vacuum phase: setting noCutoff nonbonded method.") # Select barostat - # TODO: Does this need to omit this barostat for the vacuum phase? from simtk.openmm.app import NoCutoff, CutoffNonPeriodic NONPERIODIC_NONBONDED_METHODS = [NoCutoff, CutoffNonPeriodic] if pressure is not None: @@ -323,6 +322,11 @@ def __init__(self, ligand_input, old_ligand_index, new_ligand_index, forcefield_ _logger.info(f"assgning noCutoff to nonbonded_method") self._nonbonded_method = app.NoCutoff + self._system_generator.barostat = None + _logger.info(f'Removing barostat for vacuum phase') + self._system_generator.forcefield_kwargs['nonbondedMethod'] = self._nonbonded_method + _logger.info(f'Setting nonbondedMethod to NoCutoff for vacuum phase') + _logger.info(f"calling SystemGenerator to create ligand systems.") if self._proposal_phase is None: From 994b983a45d6b6fe3b04a4057f6e9bd95d5e6f8c Mon Sep 17 00:00:00 2001 From: John Chodera Date: Sun, 12 Jan 2020 12:29:29 -0800 Subject: [PATCH 14/21] Fix travis by installing prerequisites for openmm-forcefields --- .travis.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.travis.yml b/.travis.yml index 9c412b8d3..ee5bb7fe6 100644 --- a/.travis.yml +++ b/.travis.yml @@ -59,6 +59,7 @@ script: # 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 From 26b5aceac00e0c159521ca9dd0a6996ce1f31e39 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Sun, 12 Jan 2020 15:27:40 -0800 Subject: [PATCH 15/21] Fix test failure --- perses/tests/test_smc.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/perses/tests/test_smc.py b/perses/tests/test_smc.py index a12e56a4b..b6fcdbfb7 100644 --- a/perses/tests/test_smc.py +++ b/perses/tests/test_smc.py @@ -51,7 +51,8 @@ def sMC_setup(): fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", old_ligand_index = 0, new_ligand_index = 1, - forcefield_files = ['gaff.xml'], + forcefield_files = [], + small_molecule_forcefield = 'gaff-2.11', phases = ['vacuum']) hybrid_factory = HybridTopologyFactory(topology_proposal = fe_setup._vacuum_topology_proposal, @@ -232,7 +233,8 @@ def test_create_endstates(): fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", old_ligand_index = 0, new_ligand_index = 1, - forcefield_files = ['gaff.xml'], + forcefield_files = [], + small_molecule_forcefield = 'gaff-2.11', phases = ['vacuum']) hybrid_factory = HybridTopologyFactory(topology_proposal = fe_setup._vacuum_topology_proposal, From caea5234e9fb2a69efaed6c125bebb30d4e3fe4c Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 14 Jan 2020 14:23:03 -0700 Subject: [PATCH 16/21] Fix test.smi path --- perses/tests/test_smc.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/perses/tests/test_smc.py b/perses/tests/test_smc.py index b6fcdbfb7..8d9a8c918 100644 --- a/perses/tests/test_smc.py +++ b/perses/tests/test_smc.py @@ -48,7 +48,8 @@ def sMC_setup(): """ function to setup local sMC """ - fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", + 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 = [], @@ -230,7 +231,8 @@ def test_create_endstates(): """ test the creation of unsampled endstates """ - fe_setup = RelativeFEPSetup(ligand_input = f"{os.getcwd()}/test.smi", + 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 = [], From 58c56eddd7f0584da1539c5f1c8d00ca44179574 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 14 Jan 2020 17:23:13 -0700 Subject: [PATCH 17/21] Fix temporary directory issues --- perses/tests/test_geometry_engine.py | 53 +++++++++++++--------------- perses/tests/utils.py | 11 ++++++ 2 files changed, 36 insertions(+), 28 deletions(-) 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/utils.py b/perses/tests/utils.py index 3253d714b..60aa6cc18 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,16 @@ # UTILITIES ################################################################################] +@contextlib.contextmanager +def enter_temp_directory(): + """Create and enter a temporary directory; used as context manager.""" + temp_dir = tempfile.mkdtemp() + cwd = os.getcwd() + os.chdir(temp_dir) + yield temp_dir + os.chdir(cwd) + shutil.rmtree(temp_dir) + # TODO: Move some of these utility routines to openmoltools. class NaNException(Exception): From 2a65fb0abdd3428d8fe1bb4701281ebcc63ffd77 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 14 Jan 2020 21:57:08 -0700 Subject: [PATCH 18/21] Update enter_temp_directory --- perses/tests/utils.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/perses/tests/utils.py b/perses/tests/utils.py index 60aa6cc18..b95b4aa75 100644 --- a/perses/tests/utils.py +++ b/perses/tests/utils.py @@ -45,11 +45,14 @@ @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. From 96663ab6f9b17deaada7aa27d93254fbaf3a97fc Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 14 Jan 2020 22:11:03 -0700 Subject: [PATCH 19/21] Add test.smi --- perses/{tests => data}/test.smi | 0 perses/tests/test_smc.py | 1 + 2 files changed, 1 insertion(+) rename perses/{tests => data}/test.smi (100%) 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/tests/test_smc.py b/perses/tests/test_smc.py index 8d9a8c918..2d074dde0 100644 --- a/perses/tests/test_smc.py +++ b/perses/tests/test_smc.py @@ -48,6 +48,7 @@ def sMC_setup(): """ function to setup local sMC """ + 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, From 6985bbbd9ed9fff1917ab4cfd5b9da9da30a0176 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Tue, 14 Jan 2020 22:48:17 -0700 Subject: [PATCH 20/21] Fix failing tests --- perses/tests/test_smc.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/perses/tests/test_smc.py b/perses/tests/test_smc.py index 2d074dde0..b8ca22b13 100644 --- a/perses/tests/test_smc.py +++ b/perses/tests/test_smc.py @@ -48,7 +48,7 @@ def sMC_setup(): """ function to setup local sMC """ - from pkg_resources import resource_filename + 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, @@ -232,6 +232,7 @@ def test_create_endstates(): """ test the creation of unsampled endstates """ + 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, From 27caad32c77ff1e022ad67068c1929edd12c4ff9 Mon Sep 17 00:00:00 2001 From: Hannah Bruce Macdonald Date: Wed, 15 Jan 2020 09:26:17 -0500 Subject: [PATCH 21/21] commenting out printing of intermediate pdb file for testing --- perses/app/relative_setup.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/app/relative_setup.py b/perses/app/relative_setup.py index 5449fec0a..5a67f9cfc 100644 --- a/perses/app/relative_setup.py +++ b/perses/app/relative_setup.py @@ -565,9 +565,9 @@ def _solvate_system(self, topology, positions, model='tip3p',vacuum=False): # DEBUG: Write PDB file being fed into Modeller to check why MOL isn't being matched from simtk.openmm.app import PDBFile import os - pdb_filename = os.path.join(os.environ['PWD'], 'modeller.pdb') - with open(pdb_filename, 'w') as outfile: - PDBFile.writeFile(topology, positions, outfile) + #pdb_filename = os.path.join(os.environ['PWD'], 'modeller.pdb') + #with open(pdb_filename, 'w') as outfile: + # PDBFile.writeFile(topology, positions, outfile) modeller = app.Modeller(topology, positions) # retaining protein protonation from input files