diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml new file mode 100644 index 0000000..e41b20e --- /dev/null +++ b/.github/workflows/ci.yaml @@ -0,0 +1,84 @@ +name: "CI" +on: + pull_request: + branches: + - main + push: + branches: + - main + schedule: + # At 07:00 UTC on Monday and Thursday. + - cron: "0 7 * * *" + release: + types: + - published + +concurrency: + group: "${{ github.workflow }}-${{ github.ref }}" + cancel-in-progress: true + +defaults: + run: + shell: bash -l {0} + +jobs: + tests: + runs-on: ${{ matrix.os }}-latest + name: "💻-${{matrix.os }} 🐍-${{ matrix.python-version }}" + strategy: + fail-fast: false + matrix: + os: ["ubuntu"] + python-version: + - "3.9" + - "3.10" + - "3.11" + include: + - os: "macos" + python-version: "3.11" + + steps: + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: "Setup Micromamba" + uses: mamba-org/setup-micromamba@v1 + with: + environment-file: devtools/conda-envs/test_env.yaml + environment-name: feflow-env + cache-environment: true + cache-downloads: true + create-args: >- + python=${{ matrix.python-version }} + init-shell: bash + + - name: "Install" + run: python -m pip install --no-deps -e . + + - name: "Test imports" + run: | + # if we add more to this, consider changing to for + env vars + python -Ic "import feflow; print(feflow.__version__)" + + - name: "Environment Information" + run: | + micromamba info + micromamba list + + - name: "Run tests" + env: + # Set the OFE_SLOW_TESTS to True if running a Cron job + OFE_SLOW_TESTS: ${{ fromJSON('{"false":"false","true":"true"}')[github.event_name != 'pull_request'] }} + run: | + pytest -n auto -v --cov=feflow --cov-report=xml --durations=10 + + - name: codecov + if: ${{ github.repository == 'choderalab/feflow' + && github.event_name == 'pull_request' }} + uses: codecov/codecov-action@v3 + with: + token: ${{ secrets.CODECOV_TOKEN }} + file: coverage.xml + fail_ci_if_error: False + verbose: True diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 0849b26..5456c86 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -2,26 +2,18 @@ name: feflow-test channels: - conda-forge - defaults - - openeye # TODO: Remove once we don't depend on openeye dependencies: # Base depends + - gufe - openfe # TODO: Remove once we don't depend on openfe - python - pip - # Perses depends (TODO: Remove once we don't depend on perses) - - openmoltools - - cloudpathlib - - dask - - openeye-toolkits - # Testing - pytest - pytest-cov + - pytest-xdist - codecov - # Pip-only installs - - pip: - - git+https://github.com/choderalab/perses.git@protocol-neqcyc # TODO: Remove once it only depends on openmmtools/openfe # - codecov diff --git a/feflow/__init__.py b/feflow/__init__.py index 9aef5c1..3add673 100644 --- a/feflow/__init__.py +++ b/feflow/__init__.py @@ -1,6 +1,5 @@ -"""Recipes and protocols for molecular free energy calculations using the openmmtools/perses and Open Free Energy toolkits""" - -# Add imports here -from .feflow import * +"""Recipes and protocols for molecular free energy calculations using the +openmmtools/perses and Open Free Energy toolkits""" from importlib.metadata import version +__version__ = version("openfe") diff --git a/feflow/protocols/nonequilibrium_cycling.py b/feflow/protocols/nonequilibrium_cycling.py index 3316a46..2e48107 100644 --- a/feflow/protocols/nonequilibrium_cycling.py +++ b/feflow/protocols/nonequilibrium_cycling.py @@ -6,7 +6,7 @@ import logging import time -from gufe.settings import Settings, OpenMMSystemGeneratorFFSettings, ThermoSettings +from gufe.settings import Settings from gufe.chemicalsystem import ChemicalSystem from gufe.mapping import ComponentMapping from gufe.protocols import ( @@ -16,12 +16,12 @@ ProtocolDAGResult, ) -from perses.app.relative_setup import RelativeFEPSetup -from perses.app.setup_relative_calculation import get_openmm_platform -from perses.annihilation.relative import HybridTopologyFactory +# TODO: Remove/change when things get migrated to openmmtools or feflow +from openfe.protocols.openmm_utils import system_creation +from openfe.protocols.openmm_rfe._rfe_utils.compute import get_openmm_platform from openff.units import unit -from openff.units.openmm import to_openmm +from openff.units.openmm import to_openmm, from_openmm from ..utils.data import serialize @@ -94,6 +94,58 @@ def _detect_phase(state_a, state_b): return detected_phase + @staticmethod + def extract_positions(context, hybrid_topology_factory, atom_selection_exp="not water"): + """ + Extract positions from initial and final systems based from the hybrid topology. + + Parameters + ---------- + context: openmm.Context + Current simulation context where from extract positions. + hybrid_topology_factory: perses.annihilation.relative.HybridTopologyFactory + Hybrid topology factory where to extract positions and mapping information + atom_selection_exp: str, optional + Atom selection expression using mdtraj syntax. Defaults to "not water" + + Returns + ------- + + Notes + ----- + It achieves this by taking the positions and indices from the initial and final states of + the transformation, and computing the overlap of these with the indices of the complete + hybrid topology, filtered by some mdtraj selection expression. + + 1. Get positions from context + 2. Get topology from HTF (already mdtraj topology) + 3. Merge that information into mdtraj.Trajectory + 4. Filter positions for initial/final according to selection string + """ + # TODO: Maybe we want this as a helper/utils function in perses. We also need tests for this. + import mdtraj as md + import numpy as np + + # Get positions from current openmm context + positions = context.getState(getPositions=True).getPositions(asNumpy=True) + + # Get topology from HTF - indices for initial and final topologies in hybrid topology + initial_indices = np.asarray(hybrid_topology_factory.initial_atom_indices) + final_indices = np.asarray(hybrid_topology_factory.final_atom_indices) + hybrid_topology = hybrid_topology_factory.hybrid_topology + selection = atom_selection_exp + md_trajectory = md.Trajectory(xyz=positions, topology=hybrid_topology) + selection_indices = md_trajectory.topology.select(selection) + + # Now we have to find the intersection/overlap between selected indices in the hybrid + # topology and the initial/final positions, respectively + initial_selected_indices = np.intersect1d(initial_indices, selection_indices) + final_selected_indices = np.intersect1d(final_indices, selection_indices) + initial_selected_positions = md_trajectory.xyz[0, initial_selected_indices, :] + final_selected_positions = md_trajectory.xyz[0, final_selected_indices, :] + + return initial_selected_positions, final_selected_positions + def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): """ Execute the simulation part of the Nonequilibrium switching protocol using GUFE objects. @@ -126,8 +178,11 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): import numpy as np import openmm import openmm.unit as openmm_unit + from openff.units.openmm import ensure_quantity from openmmtools.integrators import PeriodicNonequilibriumIntegrator - from perses.utils.openeye import generate_unique_atom_names + from gufe.components import SmallMoleculeComponent + from openfe.protocols.openmm_rfe import _rfe_utils + from feflow.utils.hybrid_topology import HybridTopologyFactoryModded as HybridTopologyFactory # Check compatibility between states (same receptor and solvent) self._check_states_compatibility(state_a, state_b) @@ -140,70 +195,127 @@ def _execute(self, ctx, *, state_a, state_b, mapping, settings, **inputs): solvent_a = state_a.components.get("solvent") # solvent_b = state_b.components.get("solvent") # Should not be needed - # Check first state for receptor if not get receptor from second one - if receptor_a: - receptor_top = receptor_a.to_openmm_topology() - receptor_pos = receptor_a.to_openmm_positions() - else: - receptor_top, receptor_pos = None, None - - # Get ligands cheminformatics molecules - ligand_a = ligand_a.to_openff().to_openeye() - ligand_b = ligand_b.to_openff().to_openeye() - # Generating unique atom names for ligands -- openmmforcefields needs them - ligand_a = generate_unique_atom_names(ligand_a) - ligand_b = generate_unique_atom_names(ligand_b) - - # Get solvent parameters from component - if solvent_a: - ion_concentration = solvent_a.ion_concentration.to_openmm() - positive_ion = solvent_a.positive_ion - negative_ion = solvent_a.negative_ion + # Get settings for system generator + forcefield_settings = settings.forcefield_settings + thermodynamic_settings = settings.thermo_settings + system_settings = settings.system_settings + + # handle cache for system generator + if settings.forcefield_cache is not None: + ffcache = ctx.shared / settings.forcefield_cache else: - ion_concentration, positive_ion, negative_ion = None, None, None + ffcache = None + + system_generator = system_creation.get_system_generator( + forcefield_settings=forcefield_settings, + thermo_settings=thermodynamic_settings, + system_settings=system_settings, + cache=ffcache, + has_solvent=solvent_a is not None, + ) - # Get settings - thermodynamic_settings = settings.thermo_settings - phase = self._detect_phase(state_a, state_b) - traj_save_frequency = settings.traj_save_frequency - work_save_frequency = settings.work_save_frequency # Note: this is divisor of traj save freq. - selection_expression = settings.atom_selection_expression + # b. force the creation of parameters + # This is necessary because we need to have the FF generated ahead of + # solvating the system. + # Note: by default this is cached to ctx.shared/db.json so shouldn't + # incur too large a cost + self.logger.info("Parameterizing molecules") + small_mols_a = [] + for comp in state_a.components.values(): + if isinstance(comp, SmallMoleculeComponent): + small_mols_a.append(comp) + + for comp in small_mols_a: + offmol = comp.to_openff() + system_generator.create_system(offmol.to_topology().to_openmm(), + molecules=[offmol]) + if comp == ligand_a: + mol_b = ligand_b.to_openff() + system_generator.create_system(mol_b.to_topology().to_openmm(), + molecules=[mol_b]) + + # c. get OpenMM Modeller + a dictionary of resids for each component + solvation_settings = settings.solvation_settings + state_a_modeller, comp_resids = system_creation.get_omm_modeller( + protein_comp=receptor_a, + solvent_comp=solvent_a, + small_mols=small_mols_a, + omm_forcefield=system_generator.forcefield, + solvent_settings=solvation_settings, + ) + + # d. get topology & positions + # Note: roundtrip positions to remove vec3 issues + state_a_topology = state_a_modeller.getTopology() + state_a_positions = to_openmm( + from_openmm(state_a_modeller.getPositions()) + ) - # Get the ligand mapping from ComponentMapping object - # NOTE: perses to date has a different "directionality" sense in terms of the mapping, - # see perses.rjmc.topology_proposal.propose docstring for detailed information. - ligand_mapping = mapping['ligand'].componentB_to_componentA - - # Setup relative FE calculation - fe_setup = RelativeFEPSetup( - old_ligand=ligand_a, - new_ligand=ligand_b, - receptor=receptor_top, - receptor_positions=receptor_pos, - forcefield_files=settings.forcefield_settings.forcefields, - small_molecule_forcefield=settings.forcefield_settings.small_molecule_forcefield, - phases=[phase], - transformation_atom_map=ligand_mapping, # Handle atom mapping between systems - ionic_strength=ion_concentration, - positive_ion=positive_ion, - negative_ion=negative_ion, + # e. create the stateA System + state_a_system = system_generator.create_system( + state_a_modeller.topology, + molecules=[s.to_openff() for s in small_mols_a], ) - topology_proposals = fe_setup.topology_proposals - old_positions = fe_setup.old_positions - new_positions = fe_setup.new_positions - - # Generate Hybrid Topology Factory - Generic HTF - htf = HybridTopologyFactory( - topology_proposal=topology_proposals[phase], - current_positions=old_positions[phase], - new_positions=new_positions[phase], - softcore_LJ_v2=settings.softcore_LJ_v2, - interpolate_old_and_new_14s=settings.interpolate_old_and_new_14s, + # 2. Get stateB system + # a. get the topology + state_b_topology, state_b_alchem_resids = _rfe_utils.topologyhelpers.combined_topology( + state_a_topology, + ligand_b.to_openff().to_topology().to_openmm(), + exclude_resids=comp_resids[ligand_a], ) - system = htf.hybrid_system - positions = htf.hybrid_positions + # b. get a list of small molecules for stateB + off_mols_state_b = [ligand_b.to_openff(), ] + for comp in small_mols_a: + if comp != ligand_a: + off_mols_state_b.append(comp.to_openff()) + + state_b_system = system_generator.create_system( + state_b_topology, + molecules=off_mols_state_b, + ) + + # c. Define correspondence mappings between the two systems + ligand_mappings = _rfe_utils.topologyhelpers.get_system_mappings( + mapping["ligand"].componentA_to_componentB, + state_a_system, state_a_topology, comp_resids[ligand_a], + state_b_system, state_b_topology, state_b_alchem_resids, + # These are non-optional settings for this method + fix_constraints=True, + ) + + # d. Finally get the positions + state_b_positions = _rfe_utils.topologyhelpers.set_and_check_new_positions( + ligand_mappings, state_a_topology, state_b_topology, + old_positions=ensure_quantity(state_a_positions, 'openmm'), + insert_positions=ensure_quantity(ligand_b.to_openff().conformers[0], 'openmm'), + ) + + # Get alchemical settings + alchemical_settings = settings.alchemical_settings + + # Now we can create the HTF from the previous objects + hybrid_factory = HybridTopologyFactory( + state_a_system, state_a_positions, state_a_topology, + state_b_system, state_b_positions, state_b_topology, + old_to_new_atom_map=ligand_mappings['old_to_new_atom_map'], + old_to_new_core_atom_map=ligand_mappings['old_to_new_core_atom_map'], + use_dispersion_correction=alchemical_settings.use_dispersion_correction, + softcore_alpha=alchemical_settings.softcore_alpha, + softcore_LJ_v2=alchemical_settings.softcore_LJ_v2, + softcore_LJ_v2_alpha=alchemical_settings.softcore_alpha, + interpolate_old_and_new_14s=alchemical_settings.interpolate_old_and_new_14s, + flatten_torsions=alchemical_settings.flatten_torsions, + ) + ####### END OF SETUP ######### + + traj_save_frequency = settings.traj_save_frequency + work_save_frequency = settings.work_save_frequency # Note: this is divisor of traj save freq. + selection_expression = settings.atom_selection_expression + + system = hybrid_factory.hybrid_system + positions = hybrid_factory.hybrid_positions # Set up integrator temperature = to_openmm(thermodynamic_settings.temperature) @@ -376,13 +488,14 @@ def _execute(self, ctx, *, setup, settings, **inputs): # Save positions if step % traj_save_frequency == 0: file_logger.debug(f"coarse step: {step}: saving trajectory (freq {traj_save_frequency})") - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) # Make sure trajectories are stored at the end of the eq loop file_logger.debug(f"coarse step: {step}: saving trajectory (freq {traj_save_frequency})") - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) @@ -399,12 +512,13 @@ def _execute(self, ctx, *, setup, settings, **inputs): integrator.step(work_save_frequency) forward_works.append(integrator.get_protocol_work(dimensionless=True)) if fwd_step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_neq_old.append(initial_positions) forward_neq_new.append(final_positions) # Make sure trajectories are stored at the end of the neq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_neq_old.append(initial_positions) forward_neq_new.append(final_positions) @@ -417,12 +531,13 @@ def _execute(self, ctx, *, setup, settings, **inputs): for step in range(coarse_eq_steps): integrator.step(work_save_frequency) if step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_eq_new.append(initial_positions) # TODO: Maybe better naming not old/new but initial/final reverse_eq_old.append(final_positions) # Make sure trajectories are stored at the end of the eq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_eq_old.append(initial_positions) reverse_eq_new.append(final_positions) @@ -437,12 +552,13 @@ def _execute(self, ctx, *, setup, settings, **inputs): integrator.step(work_save_frequency) reverse_works.append(integrator.get_protocol_work(dimensionless=True)) if rev_step % traj_save_frequency == 0: - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, + hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) reverse_neq_old.append(initial_positions) reverse_neq_new.append(final_positions) # Make sure trajectories are stored at the end of the neq loop - initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=htf, + initial_positions, final_positions = self.extract_positions(context, hybrid_topology_factory=hybrid_factory, atom_selection_exp=selection_expression) forward_eq_old.append(initial_positions) forward_eq_new.append(final_positions) @@ -462,6 +578,7 @@ def _execute(self, ctx, *, setup, settings, **inputs): file_logger.info(f"replicate_{self.name} Estimated performance: {estimated_performance} ns/day") # Serialize works + phase = self._detect_phase(state_a, state_b) # infer phase from systems and components forward_work_path = ctx.shared / f"forward_{phase}_{self.name}.npy" reverse_work_path = ctx.shared / f"reverse_{phase}_{self.name}.npy" with open(forward_work_path, 'wb') as out_file: @@ -679,11 +796,17 @@ def __init__(self, settings: Settings): @classmethod def _default_settings(cls): - from perses.protocols.settings import NonEquilibriumCyclingSettings + from feflow.settings.nonequilibrium_cycling import NonEquilibriumCyclingSettings + from gufe.settings import OpenMMSystemGeneratorFFSettings, ThermoSettings + from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings + from openfe.protocols.openmm_rfe.equil_rfe_settings import AlchemicalSettings return NonEquilibriumCyclingSettings( - forcefield_settings=OpenMMSystemGeneratorFFSettings(), - thermo_settings=ThermoSettings(temperature=300 * unit.kelvin), - ) + forcefield_settings=OpenMMSystemGeneratorFFSettings(), + thermo_settings=ThermoSettings(temperature=300 * unit.kelvin), + system_settings=SystemSettings(), + solvation_settings=SolvationSettings(), + alchemical_settings=AlchemicalSettings(), + ) # NOTE: create method should be really fast, since it would be running in the work units not the clients!! def _create( diff --git a/feflow/settings/neq_cycling.py b/feflow/settings/nonequilibrium_cycling.py similarity index 81% rename from feflow/settings/neq_cycling.py rename to feflow/settings/nonequilibrium_cycling.py index 0b547b6..e2c97d3 100644 --- a/feflow/settings/neq_cycling.py +++ b/feflow/settings/nonequilibrium_cycling.py @@ -4,10 +4,12 @@ This module implements the objects that will be needed to run relative binding free energy calculations using perses. """ - +from typing import Optional from gufe.settings import Settings from openff.units import unit from pydantic import root_validator +from openfe.protocols.openmm_utils.omm_settings import SystemSettings, SolvationSettings +from openfe.protocols.openmm_rfe.equil_rfe_settings import AlchemicalSettings # Default settings for the lambda functions x = 'lambda' @@ -44,12 +46,18 @@ class NonEquilibriumCyclingSettings(Settings): class Config: arbitrary_types_allowed = True + # System Settings (from openfe) + system_settings: SystemSettings + forcefield_cache: Optional[str] = "db.json" # TODO: Remove once it has been integrated with openfe settings + + # Solvation settings + solvation_settings: SolvationSettings + # Lambda settings lambda_functions = DEFAULT_ALCHEMICAL_FUNCTIONS # alchemical settings - softcore_LJ_v2 = True - interpolate_old_and_new_14s = False + alchemical_settings: AlchemicalSettings # NEQ integration settings timestep = 4.0 * unit.femtoseconds @@ -63,7 +71,7 @@ class Config: work_save_frequency: int = 500 atom_selection_expression: str = "not water" - # Number of cycles to run + # Number of replicates to run (1 cycle/replicate) num_replicates: int = 1 @root_validator diff --git a/feflow/tests/conftest.py b/feflow/tests/conftest.py index 0f10b4a..d9473c2 100644 --- a/feflow/tests/conftest.py +++ b/feflow/tests/conftest.py @@ -69,7 +69,7 @@ def toluene_solvent_system(toluene, solvent_comp): @pytest.fixture def short_settings(): from openff.units import unit - from perses.protocols import NonEquilibriumCyclingProtocol + from feflow.protocols import NonEquilibriumCyclingProtocol settings = NonEquilibriumCyclingProtocol.default_settings() @@ -94,7 +94,7 @@ def short_settings_gpu(short_settings): @pytest.fixture def short_settings_multiple_cycles(): from openff.units import unit - from perses.protocols import NonEquilibriumCyclingProtocol + from feflow.protocols import NonEquilibriumCyclingProtocol settings = NonEquilibriumCyclingProtocol.default_settings() diff --git a/feflow/utils/hybrid_topology.py b/feflow/utils/hybrid_topology.py new file mode 100644 index 0000000..ee73a22 --- /dev/null +++ b/feflow/utils/hybrid_topology.py @@ -0,0 +1,104 @@ +from openfe.protocols.openmm_rfe._rfe_utils.relative import HybridTopologyFactory + + +# TODO: This is an utility class. To be deleted when we migrate/extend the base HybridTopologyFactory +class HybridTopologyFactoryModded(HybridTopologyFactory): + """ + Utility class that extends the base HybridTopologyFactory class with properties for + getting the indices from initial and final states. + """ + + def __init__(self, + old_system, old_positions, old_topology, + new_system, new_positions, new_topology, + old_to_new_atom_map, old_to_new_core_atom_map, + use_dispersion_correction=False, + softcore_alpha=0.5, + softcore_LJ_v2=True, + softcore_LJ_v2_alpha=0.85, + interpolate_old_and_new_14s=False, + flatten_torsions=False, + **kwargs): + """ + Initialize the Hybrid topology factory. + + Parameters + ---------- + old_system : openmm.System + OpenMM system defining the "old" (i.e. starting) state. + old_positions : [n,3] np.ndarray of float + The positions of the "old system". + old_topology : openmm.Topology + OpenMM topology defining the "old" state. + new_system: opemm.System + OpenMM system defining the "new" (i.e. end) state. + new_positions : [m,3] np.ndarray of float + The positions of the "new system" + new_topology : openmm.Topology + OpenMM topology defining the "new" state. + old_to_new_atom_map : dict of int : int + Dictionary of corresponding atoms between the old and new systems. + Unique atoms are not included in this atom map. + old_to_new_core_atom_map : dict of int : int + Dictionary of corresponding atoms between the alchemical "core + atoms" (i.e. residues which are changing) between the old and + new systems. + use_dispersion_correction : bool, default False + Whether to use the long range correction in the custom sterics + force. This can be very expensive for NCMC. + softcore_alpha: float, default None + "alpha" parameter of softcore sterics, default 0.5. + softcore_LJ_v2 : bool, default True + Implement the softcore LJ as defined by Gapsys et al. JCTC 2012. + softcore_LJ_v2_alpha : float, default 0.85 + Softcore alpha parameter for LJ v2 + interpolate_old_and_new_14s : bool, default False + Whether to turn off interactions for new exceptions (not just + 1,4s) at lambda = 0 and old exceptions at lambda = 1; if False, + they are present in the nonbonded force. + flatten_torsions : bool, default False + If True, torsion terms involving `unique_new_atoms` will be + scaled such that at lambda=0,1, the torsion term is turned off/on + respectively. The opposite is true for `unique_old_atoms`. + """ + super().__init__(old_system, old_positions, old_topology, + new_system, new_positions, new_topology, + old_to_new_atom_map, old_to_new_core_atom_map, + use_dispersion_correction=use_dispersion_correction, + softcore_alpha=softcore_alpha, + softcore_LJ_v2=softcore_LJ_v2, + softcore_LJ_v2_alpha=softcore_LJ_v2_alpha, + interpolate_old_and_new_14s=interpolate_old_and_new_14s, + flatten_torsions=flatten_torsions, + **kwargs) + + # TODO: We need to refactor for the init to use these properties and have an attribute with the indices + @property + def initial_atom_indices(self): + """ + Indices of atoms in hybrid topology for atoms in the old topology. + + Returns + ------- + hybrid_indices: list[int] + Indices of old atoms in hybrid topology + + """ + n_atoms_old = self._old_system.getNumParticles() + hybrid_indices = [self._old_to_hybrid_map[idx] for idx in range(n_atoms_old)] + return hybrid_indices + + @property + def final_atom_indices(self): + """ + Indices of atoms in hybrid topology for atoms in the new topology. + + Returns + ------- + hybrid_indices: list[int] + Indices of new atoms in hybrid topology + + """ + n_atoms_new = self._new_system.getNumParticles() + hybrid_indices = [self._new_to_hybrid_map[idx] for idx in range(n_atoms_new)] + return hybrid_indices diff --git a/pyproject.toml b/pyproject.toml index 683fa9c..c81476a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,7 +10,7 @@ description = "Recipes and protocols for molecular free energy calculations usin dynamic = ["version"] readme = "README.md" authors = [ - { name = "Iván Pulido", email = "ivan.pulido@choderalab.org" } + { name = "ChoderaLab", email = "ivan.pulido@choderalab.org" } ] license = { text = "MIT" } # See https://pypi.org/classifiers/ @@ -67,7 +67,7 @@ feflow = [ ] [tool.versioningit] -default-version = "1+unknown" +default-version = "0.1.0+unknown" [tool.versioningit.format] distance = "{base_version}+{distance}.{vcs}{rev}" @@ -79,7 +79,7 @@ distance-dirty = "{base_version}+{distance}.{vcs}{rev}.dirty" method = "git" # <- The method name # Parameters to pass to the method: match = ["*"] -default-tag = "1.0.0" +default-tag = "0.1.0" [tool.versioningit.write] file = "feflow/_version.py"