Skip to content

atm #

Automated BFE calculations using the alchemical transfer method

Classes:

  • ATMSoftCore

    Configuration for the ATM soft-core potential.

  • ATMAlignmentRestraint

    Configuration for an ATM alignment restraint.

  • ATMRestraints

    Configure the restraints that will be applied during the ATM calculations.

  • ATMReferenceSelection

    Configure how receptor binding sites and ligand alignment reference atoms are

  • ATMSetupStage

    Configure how the complex will be solvated and restrained prior to

  • ATMStates

    Configure the lambda schedules.

  • ATMEquilibrateStage

    Configure how the system will be equilibrated prior to replica exchange.

  • ATMSamplingStage

    Configure how the system will be sampled using Hamiltonian replica exchange.

  • ATMConfig

    Configuration the stages of the ATM calculation.

  • ATMEdge

    Defines an ATM specific edge in a free energy network.

  • ATMNetwork

    Defines an ATM specific free energy network.

Functions:

  • load_config

    Load a configuration from a YAML file.

  • create_state_dicts

    Map the lambda states specified in the configuration to a list of dictionaries.

  • compute_ddg

    Computes the total binding free energy.

  • equilibrate_states

    Equilibrate the system at each lambda window.

  • run_hremd

    Perform replica exchange sampling for a system prepared for ATM calculations.

  • setup_system

    Prepares a system ready for running the ATM method.

  • select_displacement

    Attempts to automatically select a displacement vector for the ligands.

  • run_workflow

    Run the setup, equilibration, and sampling phases.

  • submit_network

    Submits a set of ATM calculations to the SLURM queueing manager.

Attributes:

DEFAULT_LAMBDA_1 module-attribute #

DEFAULT_LAMBDA_1 = (
    0.0,
    0.0,
    0.0,
    0.0,
    0.0,
    0.0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.5,
    0.4,
    0.3,
    0.2,
    0.1,
    0.0,
    0.0,
    0.0,
    0.0,
    0.0,
    0.0,
)

The default lambda 1 schedule.

DEFAULT_LAMBDA_2 module-attribute #

DEFAULT_LAMBDA_2 = (
    0.0,
    0.1,
    0.2,
    0.3,
    0.4,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.5,
    0.4,
    0.3,
    0.2,
    0.1,
    0.0,
)

The default lambda 2 schedule.

DEFAULT_DIRECTION module-attribute #

DEFAULT_DIRECTION = (
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
    -1,
)

The default direction schedule.

DEFAULT_ALPHA module-attribute #

DEFAULT_ALPHA = (
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
    0.1,
)

The default alpha schedule.

DEFAULT_U0 module-attribute #

DEFAULT_U0 = (
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
    110.0,
)

The default u0 schedule.

DEFAULT_W0 module-attribute #

DEFAULT_W0 = (
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
)

The default w0 schedule.

DEFAULT_MAX_REST_TEMPERATURE module-attribute #

DEFAULT_MAX_REST_TEMPERATURE = 900.0

The default maximum temperature to use during the REST2 calculations.

DEFAULT_BM_B0 module-attribute #

DEFAULT_BM_B0 = tuple(tolist())

The default beta scaling factors to use if running with REST2.

DEFAULT_EQUILIBRATE_INTEGRATOR module-attribute #

DEFAULT_EQUILIBRATE_INTEGRATOR = LangevinIntegrator(
    timestep=2.0 * femtosecond, friction=1.0 / picosecond
)

The default integrator to use during equilibration.

DEFAULT_EQUILIBRATE_RESTRAINTS module-attribute #

DEFAULT_EQUILIBRATE_RESTRAINTS = {
    DEFAULT_RESTRAINT_MASK: FlatBottomRestraint(
        k=25.0 * _KCAL_PER_ANG_SQR, radius=1.5 * _ANGSTROM
    )
}

The default position restraints to apply during equilibration.

ATMSoftCore pydantic-model #

Bases: BaseModel

Configuration for the ATM soft-core potential.

Fields:

u_max pydantic-field #

u_max: OpenMMQuantity[_KCAL_PER_MOL]

The 'u max' [kcal/mol] parameter.

u0 pydantic-field #

u0: OpenMMQuantity[_KCAL_PER_MOL]

The 'u0' [kcal/mol] parameter.

a pydantic-field #

a: float

The 'a' parameter.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMAlignmentRestraint pydantic-model #

Bases: BaseModel

Configuration for an ATM alignment restraint.

Fields:

k_distance pydantic-field #

k_distance: OpenMMQuantity[_KCAL_PER_ANG_SQR]

Force constant [kcal/mol/Å^2] of the flat-bottom potential restraining the distance between the ligands.

k_angle pydantic-field #

k_angle: OpenMMQuantity[_KCAL_PER_MOL]

Force constant [kcal/mol] of the 1-cos(theta) potential restraining the angle between the ligands.

k_dihedral pydantic-field #

k_dihedral: OpenMMQuantity[_KCAL_PER_MOL]

Force constant [kcal/mol] of the 1-cos(phi) potential restraining the dihedral angle between the ligands.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMRestraints pydantic-model #

Bases: BaseModel

Configure the restraints that will be applied during the ATM calculations.

Fields:

com pydantic-field #

com: FlatBottomRestraint = FlatBottomRestraint(
    k=25.0 * _KCAL_PER_ANG_SQR, radius=5.0 * _ANGSTROM
)

The potential that restrains the ligands to the binding site.

alignment pydantic-field #

alignment: ATMAlignmentRestraint | None = (
    ATMAlignmentRestraint(
        k_distance=2.5 * _KCAL_PER_ANG_SQR,
        k_angle=25.0 * _KCAL_PER_MOL,
        k_dihedral=25.0 * _KCAL_PER_MOL,
    )
)

The potential that restrains the orientation of the two ligands during an RBFE calculation.

receptor pydantic-field #

receptor: FlatBottomRestraint = FlatBottomRestraint(
    k=25.0 * _KCAL_PER_ANG_SQR, radius=1.5 * _ANGSTROM
)

The potential that restrains specified receptor atoms to their initial coordinates.

receptor_query pydantic-field #

receptor_query: str = '@CA'

An Amber query used to identify which receptor atoms to restrain.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMReferenceSelection pydantic-model #

Bases: BaseModel

Configure how receptor binding sites and ligand alignment reference atoms are selected if the user does not explicitly provide them.

Fields:

receptor_cutoff pydantic-field #

receptor_cutoff: OpenMMQuantity[_ANGSTROM] = 5.0 * _ANGSTROM

The minimum distance between a residues' alpha carbon and a ligand atom to be considered part of the binding site.

ligand_method pydantic-field #

ligand_method: LigandReferenceMethod = 'chen'

The default method to use to select ligand reference atoms during RBFE calculations if the user does not explicitly provide them.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMSetupStage pydantic-model #

Bases: BaseModel

Configure how the complex will be solvated and restrained prior to equilibration

Fields:

Validators:

displacement pydantic-field #

displacement: (
    OpenMMQuantity[_ANGSTROM]
    | tuple[
        OpenMMQuantity[_ANGSTROM],
        OpenMMQuantity[_ANGSTROM],
        OpenMMQuantity[_ANGSTROM],
    ]
) = (
    38.0 * _ANGSTROM
)

The distance to displace ligands from the binding site along an automatically selected displacement vector, or the vector to displace the ligands by.

solvent pydantic-field #

solvent: Solvent = Solvent()

Control how the system should be solvated.

reference pydantic-field #

Selection of receptor and ligand reference atoms.

restraints pydantic-field #

restraints: ATMRestraints = ATMRestraints()

Control how the system should be restrained.

apply_hmr pydantic-field #

apply_hmr: bool = True

Whether to aply hydrogen mass repartitioning to the system.

hydrogen_mass pydantic-field #

hydrogen_mass: OpenMMQuantity[amu] = 1.5 * amu

The mass to assign to hydrogen atoms when applying HMR.

apply_rest pydantic-field #

apply_rest: bool = False

Whether to prepare the system for REST sampling.

rest_config pydantic-field #

rest_config: REST | None = REST(
    scale_nonbonded=True,
    scale_torsions=True,
    scale_angles=False,
    scale_bonds=False,
)

The REST configuration to use if apply_rest is True.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMStates pydantic-model #

Bases: BaseModel

Configure the lambda schedules.

Fields:

  • lambda_1 (list[float])
  • lambda_2 (list[float])
  • direction (list[Literal[-1, 1]])
  • alpha (list[float])
  • alpha_unit (ClassVar)
  • u0 (list[float])
  • u0_unit (ClassVar)
  • w0 (list[float])
  • w0_unit (ClassVar)
  • bm_b0 (list[float] | None)

Validators:

  • _validate_lambda_lengths

bm_b0 pydantic-field #

bm_b0: list[float] | None = None

The REST2 beta scaling factors (beta_m / beta_0) to use. apply_rest must be set to true in the setup config if using this.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMEquilibrateStage pydantic-model #

Bases: BaseModel

Configure how the system will be equilibrated prior to replica exchange.

Fields:

report_interval pydantic-field #

report_interval: int = 5000

The number of steps to report energy, volume, etc after.

soft_core pydantic-field #

soft_core: ATMSoftCore = ATMSoftCore(
    u_max=1000 * _KCAL_PER_MOL,
    u0=500 * _KCAL_PER_MOL,
    a=1.0 / 16.0,
)

The ATM soft-core potential parameters to use during equilibration.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMSamplingStage pydantic-model #

Bases: HREMD

Configure how the system will be sampled using Hamiltonian replica exchange.

Fields:

temperature pydantic-field #

temperature: OpenMMQuantity[kelvin] = DEFAULT_TEMPERATURE

The temperature to sample at.

n_warmup_steps pydantic-field #

n_warmup_steps: int = 150000

The number of steps to run each replica for before starting hremd trials. All energies gathered during this period will be discarded.

n_steps_per_cycle pydantic-field #

n_steps_per_cycle: int = 1000

The number of steps to propagate the system by before attempting an exchange.

n_cycles pydantic-field #

n_cycles: int = 2500

The number of cycles of 'propagate the system' -> 'exchange replicas' to run.

max_step_retries pydantic-field #

max_step_retries: int = 5

The maximum number of times to attempt to step if a NaN is encountered before raising an exception

swap_mode pydantic-field #

swap_mode: HREMDSwapModeLiteral | None = ALL.value

The mode in which to propose state swaps between replicas. This can either be: 'neighbours', only try and swap adjacent states or ii. 'all', try and swap all states stochastically. If None, no replica exchanges will be attempted.

max_swaps pydantic-field #

max_swaps: int | None = None

The maximum number of swap proposals to make if running in 'all' mode. This variable does nothing when running in 'neighbours' mode.

trajectory_interval pydantic-field #

trajectory_interval: int | None = None

The number of cycles to run before saving the current replica states to DCD trajectory files. If None, no trajectories will be saved.

trajectory_enforce_pbc pydantic-field #

trajectory_enforce_pbc: bool = False

Whether to apply periodic boundary conditions when retrieving coordinates for writing to trajectory files.

checkpoint_interval pydantic-field #

checkpoint_interval: int | None = None

The number of cycles to run before saving the current replica states to checkpoint files. If None, no checkpoints will be saved.

integrator pydantic-field #

integrator: LangevinIntegrator = LangevinIntegrator(
    timestep=4.0 * femtosecond, friction=1.0 / picosecond
)

The MD integrator to use.

soft_core pydantic-field #

soft_core: ATMSoftCore = ATMSoftCore(
    u_max=200 * _KCAL_PER_MOL,
    u0=100 * _KCAL_PER_MOL,
    a=1.0 / 16.0,
)

The ATM soft-core potential parameters to use.

analysis_interval pydantic-field #

analysis_interval: int | None = None

The interval (in number of cycles) between estimating and reporting the free energy. If None, no analysis will be performed.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMConfig pydantic-model #

Bases: BaseModel

Configuration the stages of the ATM calculation.

Fields:

setup pydantic-field #

Prepare the system for equilibration.

states pydantic-field #

states: ATMStates = ATMStates()

Configure the lambda schedules.

equilibrate pydantic-field #

Equilibrate the system.

sample pydantic-field #

Sample across lambda windows using replica exchange.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMEdge pydantic-model #

Bases: Edge

Defines an ATM specific edge in a free energy network.

Fields:

ligand_1 pydantic-field #

ligand_1: str

The name of the first ligand.

ligand_2 pydantic-field #

ligand_2: str | None

The name of the second ligand. This should be None if running an ABFE calculation.

ligand_1_ref_atoms pydantic-field #

ligand_1_ref_atoms: tuple[str, str, str] | None = None

Three (optional) AMBER style queries that select the atoms of the first ligand to align during an RBFE calculation.

ligand_2_ref_atoms pydantic-field #

ligand_2_ref_atoms: tuple[str, str, str] | None = None

Three (optional) AMBER style queries that select the atoms of the second ligand to align during an RBFE calculation.

ligand_1_metadata pydantic-field #

ligand_1_metadata: dict[str, Any]

Any additional metadata about ligand 1.

ligand_2_metadata pydantic-field #

ligand_2_metadata: dict[str, Any]

Any additional metadata about ligand 2.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

ATMNetwork pydantic-model #

Bases: Network

Defines an ATM specific free energy network.

Fields:

receptor pydantic-field #

receptor: str | None = None

The name of the receptor. If None, the receptor will be identified from the input directory structure

receptor_ref_query pydantic-field #

receptor_ref_query: str | None = None

An (optional) AMBER style query to manually select the receptor atoms that define the binding site. If unspecified, they will be determined automatically based on the config.

edges pydantic-field #

edges: list[ATMEdge]

The edges in the free energy network.

receptor_metadata pydantic-field #

receptor_metadata: dict[str, Any]

Any additional metadata about the receptor.

model_dump_yaml #

model_dump_yaml(
    output_path: Path | None = None, **kwargs
) -> str

Dump the model to a YAML representation.

Parameters:

  • output_path (Path | None, default: None ) –

    The (optional) path to save the YAML representation to.

Returns:

  • str

    The YAML representation.

Source code in femto/md/utils/models.py
def model_dump_yaml(self, output_path: pathlib.Path | None = None, **kwargs) -> str:
    """Dump the model to a YAML representation.

    Args:
        output_path: The (optional) path to save the YAML representation to.

    Returns:
        The YAML representation.
    """

    model_yaml = yaml.safe_dump(self.model_dump(), **kwargs)

    if output_path is not None:
        output_path.parent.mkdir(exist_ok=True, parents=True)
        output_path.write_text(model_yaml)

    return model_yaml

load_config #

load_config(path: Path) -> ATMConfig

Load a configuration from a YAML file.

Parameters:

  • path (Path) –

    The path to the YAML configuration.

Returns:

Source code in femto/fe/atm/_config.py
def load_config(path: pathlib.Path) -> ATMConfig:
    """Load a configuration from a YAML file.

    Args:
        path: The path to the YAML configuration.

    Returns:
        The loaded configuration.
    """
    return ATMConfig(**yaml.safe_load(path.read_text()))

create_state_dicts #

create_state_dicts(
    states: ATMSamplingStage,
) -> list[dict[str, float]]

Map the lambda states specified in the configuration to a list of dictionaries.

Parameters:

Returns:

  • list[dict[str, float]]

    The dictionaries of lambda states.

Source code in femto/fe/atm/_utils.py
def create_state_dicts(
    states: "femto.fe.atm.ATMSamplingStage",
) -> list[dict[str, float]]:
    """Map the lambda states specified in the configuration to a list of dictionaries.

    Args:
        states: The lambda states.

    Returns:
        The dictionaries of lambda states.
    """
    return [
        {
            openmm.ATMForce.Lambda1(): states.lambda_1[i],
            openmm.ATMForce.Lambda2(): states.lambda_2[i],
            openmm.ATMForce.Direction(): states.direction[i],
            openmm.ATMForce.Alpha(): (
                states.alpha[i] * states.alpha_unit
            ).value_in_unit(openmm.unit.kilojoules_per_mole**-1),
            openmm.ATMForce.Uh(): (states.u0[i] * states.u0_unit).value_in_unit(
                openmm.unit.kilojoules_per_mole
            ),
            openmm.ATMForce.W0(): (states.w0[i] * states.w0_unit).value_in_unit(
                openmm.unit.kilojoules_per_mole
            ),
            **(
                {femto.md.rest.REST_CTX_PARAM: states.bm_b0[i]}
                if states.bm_b0 is not None
                else {}
            ),
        }
        for i in range(len(states.lambda_1))
    ]

compute_ddg #

compute_ddg(
    config: ATMSamplingStage,
    states: ATMStates,
    u_kn: ndarray,
    n_k: ndarray,
) -> DataFrame

Computes the total binding free energy.

Parameters:

  • config (ATMSamplingStage) –

    The sampling configuration.

  • states (ATMStates) –

    The sampled states.

  • u_kn (ndarray) –

    The samples.

  • n_k (ndarray) –

    The sample counts.

Returns:

  • DataFrame

    A pandas DataFrame containing the total binding free energy and its components.

Source code in femto/fe/atm/_analyze.py
def compute_ddg(
    config: "femto.fe.atm.ATMSamplingStage",
    states: "femto.fe.atm.ATMStates",
    u_kn: numpy.ndarray,
    n_k: numpy.ndarray,
) -> "pandas.DataFrame":
    """Computes the total binding free energy.

    Args:
        config: The sampling configuration.
        states: The sampled states.
        u_kn: The samples.
        n_k: The sample counts.

    Returns:
        A pandas DataFrame containing the total binding free energy and its components.
    """
    import pandas

    n_states = len(states.lambda_1)
    n_states_leg_1 = states.direction.index(-1)

    state_groups = [(n_states_leg_1, 1.0), (n_states - n_states_leg_1, 1.0)]

    estimated, _ = femto.fe.ddg.estimate_ddg(
        u_kn, n_k, config.temperature, state_groups
    )
    return pandas.DataFrame([estimated])

equilibrate_states #

equilibrate_states(
    system: System,
    topology: Structure,
    states: ATMStates,
    config: ATMEquilibrateStage,
    offset: Quantity,
    platform: OpenMMPlatform,
    reporter: Reporter | None = None,
) -> list[State]

Equilibrate the system at each lambda window.

Parameters:

  • system (System) –

    The system to simulate.

  • topology (Structure) –

    The topology of the system to simulate.

  • states (ATMStates) –

    The states of the system to simulate.

  • config (ATMEquilibrateStage) –

    Configuration settings.

  • offset (Quantity) –

    The vector to offset the ligand by using the ATM force.

  • platform (OpenMMPlatform) –

    The accelerator to use.

  • reporter (Reporter | None, default: None ) –

    The (optional) reporter to use to record system statistics such as volume and energy.

Returns:

  • list[State]

    The final equilibrated state.

Source code in femto/fe/atm/_equilibrate.py
def equilibrate_states(
    system: openmm.System,
    topology: parmed.Structure,
    states: "femto.fe.atm.ATMStates",
    config: "femto.fe.atm.ATMEquilibrateStage",
    offset: openmm.unit.Quantity,
    platform: femto.md.constants.OpenMMPlatform,
    reporter: femto.md.reporting.Reporter | None = None,
) -> list[openmm.State]:
    """Equilibrate the system at each lambda window.

    Args:
        system: The system to simulate.
        topology: The topology of the system to simulate.
        states: The states of the system to simulate.
        config: Configuration settings.
        offset: The vector to offset the ligand by using the ATM force.
        platform: The accelerator to use.
        reporter: The (optional) reporter to use to record system statistics such as
            volume and energy.

    Returns:
        The final equilibrated state.
    """
    import femto.fe.atm._utils

    reporter = femto.md.reporting.NullReporter() if reporter is None else reporter

    openmm_reporter = femto.md.reporting.openmm.OpenMMStateReporter(
        reporter, "equilibration", config.report_interval
    )

    state_dicts = femto.fe.atm._utils.create_state_dicts(states)

    system = copy.deepcopy(system)
    femto.fe.atm._utils.add_atm_force(system, topology, config.soft_core, offset)

    equilibrated_coords = femto.md.simulate.simulate_states(
        system, topology, state_dicts, config.stages, platform, openmm_reporter
    )

    box_vectors = [
        coords.getPeriodicBoxVectors(asNumpy=True).value_in_unit(openmm.unit.angstrom)
        for coords in equilibrated_coords
    ]
    box_vectors_max = numpy.max(numpy.stack(box_vectors), axis=0) * openmm.unit.angstrom

    for i, coords in enumerate(equilibrated_coords):
        context = openmm.Context(system, openmm.VerletIntegrator(0.00001))
        context.setState(coords)
        context.setPeriodicBoxVectors(*box_vectors_max)

        equilibrated_coords[i] = context.getState(getPositions=True)

    return equilibrated_coords

run_hremd #

run_hremd(
    system: System,
    topology: Structure,
    coords: list[State],
    states: ATMStates,
    config: ATMSamplingStage,
    offset: Quantity,
    platform: OpenMMPlatform,
    output_dir: Path,
    reporter: Reporter | None = None,
)

Perform replica exchange sampling for a system prepared for ATM calculations.

Parameters:

  • system (System) –

    The system to simulate. It should not already contain an ATM force.

  • topology (Structure) –

    The topology associated with the system.

  • coords (list[State]) –

    The starting coordinates for each state.

  • states (ATMStates) –

    The lambda states to sample.

  • config (ATMSamplingStage) –

    Configuration settings.

  • offset (Quantity) –

    The vector to offset the ligand by using the ATM force.

  • platform (OpenMMPlatform) –

    The platform to run on.

  • output_dir (Path) –

    The directory to store the sampled energies and statistics to, and any trajectory files if requested.

  • reporter (Reporter | None, default: None ) –

    The reporter to log statistics such as online estimates of the free energy to.

Source code in femto/fe/atm/_sample.py
def run_hremd(
    system: openmm.System,
    topology: parmed.Structure,
    coords: list[openmm.State],
    states: "femto.fe.atm.ATMStates",
    config: "femto.fe.atm.ATMSamplingStage",
    offset: openmm.unit.Quantity,
    platform: femto.md.constants.OpenMMPlatform,
    output_dir: pathlib.Path,
    reporter: femto.md.reporting.Reporter | None = None,
):
    """Perform replica exchange sampling for a system prepared for ATM calculations.

    Args:
        system: The system to simulate. It should *not* already contain an ATM force.
        topology: The topology associated with the system.
        coords: The starting coordinates for each state.
        states: The lambda states to sample.
        config: Configuration settings.
        offset: The vector to offset the ligand by using the ATM force.
        platform: The platform to run on.
        output_dir: The directory to store the sampled energies and statistics to, and
            any trajectory files if requested.
        reporter: The reporter to log statistics such as online estimates of the
            free energy to.
    """
    import femto.fe.atm._utils

    state_dicts = femto.fe.atm._utils.create_state_dicts(states)

    system = copy.deepcopy(system)
    femto.fe.atm._utils.add_atm_force(system, topology, config.soft_core, offset)

    integrator = femto.md.utils.openmm.create_integrator(
        config.integrator, config.temperature
    )
    simulation = femto.md.utils.openmm.create_simulation(
        system, topology, coords[0], integrator, state_dicts[0], platform
    )

    n_states = len(state_dicts)

    swap_mask = {
        (i, j)
        for i in range(n_states)
        for j in range(n_states)
        if states.direction[i] != states.direction[j]
    }

    analysis_fn = None

    if reporter is not None and config.analysis_interval is not None:
        analysis_fn = functools.partial(
            _analyze, config=config, states=states, reporter=reporter
        )

    return femto.md.hremd.run_hremd(
        simulation,
        state_dicts,
        config,
        output_dir,
        swap_mask,
        initial_coords=coords,
        analysis_fn=analysis_fn,
        analysis_interval=config.analysis_interval,
    )

setup_system #

setup_system(
    config: ATMSetupStage,
    receptor: AmberParm,
    ligand_1: AmberParm,
    ligand_2: AmberParm | None,
    displacement: Quantity,
    receptor_ref_query: str | None,
    ligand_1_ref_query: tuple[str, str, str] | None = None,
    ligand_2_ref_query: tuple[str, str, str] | None = None,
) -> tuple[Structure, System]

Prepares a system ready for running the ATM method.

Returns:

  • tuple[Structure, System]

    The prepared topology and OpenMM system object.

Source code in femto/fe/atm/_setup.py
def setup_system(
    config: "femto.fe.atm.ATMSetupStage",
    receptor: parmed.amber.AmberParm,
    ligand_1: parmed.amber.AmberParm,
    ligand_2: parmed.amber.AmberParm | None,
    displacement: openmm.unit.Quantity,
    receptor_ref_query: str | None,
    ligand_1_ref_query: tuple[str, str, str] | None = None,
    ligand_2_ref_query: tuple[str, str, str] | None = None,
) -> tuple[parmed.Structure, openmm.System]:
    """Prepares a system ready for running the ATM method.

    Returns:
        The prepared topology and OpenMM system object.
    """

    _LOGGER.info(f"setting up an {'ABFE' if ligand_2 is None else 'RBFE'} calculation")

    if receptor_ref_query is None:
        # we need to select the receptor cavity atoms before offsetting any ligands
        # as the query is distance based

        _LOGGER.info("selecting receptor reference atoms")
        receptor_ref_query = femto.fe.reference.select_protein_cavity_atoms(
            receptor,
            [ligand_1] + ([] if ligand_2 is None else [ligand_2]),
            config.reference.receptor_cutoff,
        )

    ligand_1_ref_idxs, ligand_2_ref_idxs = None, None

    # we carve out a 'cavity' where the first ligand will be displaced into during the
    # ATM calculations. this should make equilibration at all states easier.
    cavity_formers = [_offset_ligand(ligand_1, displacement)]

    if ligand_2 is not None:
        # we make sure that when placing solvent molecules we don't accidentally place
        # any on top of the ligands in the cavity itself
        cavity_formers.append(ligand_2)

        (
            ligand_1_ref_idxs,
            ligand_2_ref_idxs,
        ) = femto.fe.reference.select_ligand_idxs(
            ligand_1,
            ligand_2,
            config.reference.ligand_method,
            ligand_1_ref_query,
            ligand_2_ref_query,
        )
        assert ligand_2_ref_idxs is not None, "ligand 2 ref atoms were not selected"
        ligand_2_ref_idxs = tuple(i + len(ligand_1.atoms) for i in ligand_2_ref_idxs)

        ligand_2 = _offset_ligand(ligand_2, displacement)

    _LOGGER.info("solvating system")
    topology = femto.md.solvate.solvate_system(
        receptor,
        ligand_1,
        ligand_2,
        config.solvent,
        displacement,
        cavity_formers=cavity_formers,
    )

    _LOGGER.info("creating OpenMM system")
    system = topology.createSystem(
        nonbondedMethod=openmm.app.PME,
        nonbondedCutoff=0.9 * openmm.unit.nanometer,
        constraints=openmm.app.HBonds,
        rigidWater=True,
    )

    if config.apply_hmr:
        _LOGGER.info("applying HMR.")

        hydrogen_mass = config.hydrogen_mass
        femto.md.system.apply_hmr(system, topology, hydrogen_mass)

    ligand_1_idxs = list(range(len(ligand_1.atoms)))
    ligand_2_idxs = None

    if ligand_2 is not None:
        ligand_2_idxs = [i + len(ligand_1_idxs) for i in range(len(ligand_2.atoms))]

    if config.apply_rest:
        _LOGGER.info("applying REST2.")

        solute_idxs = ligand_1_idxs + ([] if ligand_2_idxs is None else ligand_2_idxs)
        femto.md.rest.apply_rest(system, set(solute_idxs), config.rest_config)

    _LOGGER.info("applying restraints.")
    ligands = [ligand_1] + ([] if ligand_2 is None else [ligand_2])
    idx_offset = sum(len(ligand.atoms) for ligand in ligands)

    receptor_ref_mask = parmed.amber.AmberMask(receptor, receptor_ref_query).Selection()
    receptor_ref_idxs = [i + idx_offset for i, m in enumerate(receptor_ref_mask) if m]
    _LOGGER.info(f"receptor ref idxs={receptor_ref_idxs}")

    _apply_atm_restraints(
        system,
        config.restraints,
        ligand_1_com_idxs=ligand_1_idxs,
        ligand_1_ref_idxs=ligand_1_ref_idxs,
        ligand_2_com_idxs=ligand_2_idxs,
        ligand_2_ref_idxs=ligand_2_ref_idxs,
        receptor_ref_idxs=receptor_ref_idxs,
        offset=displacement,
    )

    restraint_query = config.restraints.receptor_query
    restraint_mask = parmed.amber.AmberMask(receptor, restraint_query).Selection()
    restraint_idxs = [i + idx_offset for i, match in enumerate(restraint_mask) if match]

    _apply_receptor_restraints(
        system, config.restraints, {i: topology.positions[i] for i in restraint_idxs}
    )
    femto.md.utils.openmm.assign_force_groups(system)

    return topology, system

select_displacement #

select_displacement(
    receptor: AmberParm,
    ligand_1: AmberParm,
    ligand_2: AmberParm | None,
    distance: Quantity,
) -> Quantity

Attempts to automatically select a displacement vector for the ligands.

Parameters:

  • receptor (AmberParm) –

    The receptor.

  • ligand_1 (AmberParm) –

    The first ligand positioned in the binding site.

  • ligand_2 (AmberParm | None) –

    The second ligand positioned in the binding site.

  • distance (Quantity) –

    The distance to translate ligands along the displacement vector by.

Returns:

  • Quantity

    The displacement vector.

Source code in femto/fe/atm/_setup.py
def select_displacement(
    receptor: parmed.amber.AmberParm,
    ligand_1: parmed.amber.AmberParm,
    ligand_2: parmed.amber.AmberParm | None,
    distance: openmm.unit.Quantity,
) -> openmm.unit.Quantity:
    """Attempts to automatically select a displacement vector for the ligands.

    Args:
        receptor: The receptor.
        ligand_1: The first ligand positioned in the binding site.
        ligand_2: The second ligand positioned in the binding site.
        distance: The distance to translate ligands along the displacement vector by.

    Returns:
        The displacement vector.
    """

    ligand_coords = numpy.vstack(
        [ligand_1.coordinates] + ([] if ligand_2 is None else [ligand_2.coordinates])
    )
    receptor_coords = receptor.coordinates

    directions = numpy.array(
        [
            [-1.0, -1.0, -1.0],
            [+1.0, -1.0, -1.0],
            [+1.0, +1.0, -1.0],
            [-1.0, +1.0, -1.0],
            [-1.0, -1.0, +1.0],
            [+1.0, -1.0, +1.0],
            [+1.0, +1.0, +1.0],
            [-1.0, +1.0, +1.0],
        ]
    )
    directions /= numpy.linalg.norm(directions, axis=1, keepdims=True)

    closest_distances = []

    for direction in directions:
        displacement = direction * distance.value_in_unit(openmm.unit.angstrom)

        offset_coords = ligand_coords + displacement

        distances = scipy.spatial.distance.cdist(offset_coords, receptor_coords)
        closest_distances.append(distances.min())

    direction = directions[numpy.argmax(closest_distances)]
    return direction.flatten() * distance

run_workflow #

run_workflow(
    config: ATMConfig,
    ligand_1_coords: Path,
    ligand_1_params: Path,
    ligand_2_coords: Path | None,
    ligand_2_params: Path | None,
    receptor_coords: Path,
    receptor_params: Path | None,
    output_dir: Path,
    report_dir: Path | None = None,
    displacement: Quantity | None = None,
    ligand_1_ref_atoms: tuple[str, str, str] | None = None,
    ligand_2_ref_atoms: tuple[str, str, str] | None = None,
    receptor_ref_atoms: str | None = None,
)

Run the setup, equilibration, and sampling phases.

Parameters:

  • config (ATMConfig) –

    The configuration.

  • ligand_1_coords (Path) –

    The path to the first ligand coordinates.

  • ligand_1_params (Path) –

    The path to the first ligand parameters.

  • ligand_2_coords (Path | None) –

    The path to the second ligand coordinates.

  • ligand_2_params (Path | None) –

    The path to the second ligand parameters.

  • receptor_coords (Path) –

    The path to the receptor coordinates.

  • receptor_params (Path | None) –

    The path to the receptor parameters.

  • report_dir (Path | None, default: None ) –

    The directory to write any statistics to.

  • output_dir (Path) –

    The directory to store all outputs in.

  • displacement (Quantity | None, default: None ) –

    The displacement to offset the ligands by.

  • ligand_1_ref_atoms (tuple[str, str, str] | None, default: None ) –

    The AMBER style query masks that select the first ligands' reference atoms.

  • ligand_2_ref_atoms (tuple[str, str, str] | None, default: None ) –

    The AMBER style query masks that select the second ligands' reference atoms.

  • receptor_ref_atoms (str | None, default: None ) –

    The AMBER style query mask that selects the receptor atoms that form the binding site.

Source code in femto/fe/atm/_runner.py
def run_workflow(
    config: "femto.fe.atm.ATMConfig",
    ligand_1_coords: pathlib.Path,
    ligand_1_params: pathlib.Path,
    ligand_2_coords: pathlib.Path | None,
    ligand_2_params: pathlib.Path | None,
    receptor_coords: pathlib.Path,
    receptor_params: pathlib.Path | None,
    output_dir: pathlib.Path,
    report_dir: pathlib.Path | None = None,
    displacement: openmm.unit.Quantity | None = None,
    ligand_1_ref_atoms: tuple[str, str, str] | None = None,
    ligand_2_ref_atoms: tuple[str, str, str] | None = None,
    receptor_ref_atoms: str | None = None,
):
    """Run the setup, equilibration, and sampling phases.

    Args:
        config: The configuration.
        ligand_1_coords: The path to the first ligand coordinates.
        ligand_1_params: The path to the first ligand parameters.
        ligand_2_coords: The path to the second ligand coordinates.
        ligand_2_params: The path to the second ligand parameters.
        receptor_coords: The path to the receptor coordinates.
        receptor_params: The path to the receptor parameters.
        report_dir: The directory to write any statistics to.
        output_dir: The directory to store all outputs in.
        displacement: The displacement to offset the ligands by.
        ligand_1_ref_atoms: The AMBER style query masks that select the first ligands'
            reference atoms.
        ligand_2_ref_atoms: The AMBER style query masks that select the second ligands'
            reference atoms.
        receptor_ref_atoms: The AMBER style query mask that selects the receptor atoms
            that form the binding site.
    """
    import femto.fe.atm._equilibrate
    import femto.fe.atm._sample

    reporter = (
        femto.md.reporting.NullReporter()
        if report_dir is None or not femto.md.utils.mpi.is_rank_zero()
        else femto.md.reporting.TensorboardReporter(report_dir)
    )

    topology, system, displacement = _prepare_system(
        config.setup,
        ligand_1_coords,
        ligand_1_params,
        ligand_2_coords,
        ligand_2_params,
        receptor_coords,
        receptor_params,
        displacement,
        ligand_1_ref_atoms,
        ligand_2_ref_atoms,
        receptor_ref_atoms,
        output_dir / "_setup",
    )
    topology.symmetry = None  # needed as attr is lost after pickling by MPI

    equilibrate_dir = output_dir / "_equilibrate"
    equilibrate_dir.mkdir(exist_ok=True, parents=True)

    coord_paths = [
        equilibrate_dir / f"state_{state_idx}.xml"
        for state_idx in range(len(config.states.lambda_1))
    ]

    if any(not path.exists() for path in coord_paths):
        coords = femto.fe.atm._equilibrate.equilibrate_states(
            system,
            topology,
            config.states,
            config.equilibrate,
            displacement,
            femto.md.constants.OpenMMPlatform.CUDA,
            reporter,
        )
        _cache_equilibrate_outputs(coords, coord_paths)
    else:
        coords = [
            openmm.XmlSerializer.deserialize(path.read_text()) for path in coord_paths
        ]

    sample_dir = output_dir / "_sample"
    result_path = output_dir / "ddg.csv"

    if not result_path.exists():
        femto.fe.atm._sample.run_hremd(
            system,
            topology,
            coords,
            config.states,
            config.sample,
            displacement,
            femto.md.constants.OpenMMPlatform.CUDA,
            sample_dir,
            reporter,
        )
        _analyze_results(config, sample_dir, result_path)

submit_network #

submit_network(
    config: ATMConfig,
    network: Network,
    output_dir: Path,
    queue_options: SLURMOptions,
    mpi_command: list[str] | None = None,
) -> list[str]

Submits a set of ATM calculations to the SLURM queueing manager.

Parameters:

  • config (ATMConfig) –

    The configuration.

  • network (Network) –

    The network of edges to run.

  • output_dir (Path) –

    The directory to store any outputs in.

  • queue_options (SLURMOptions) –

    The options to use when submitting the jobs.

  • mpi_command (list[str] | None, default: None ) –

    The mpi runner command to use. The default is "srun --mpi=pmix".

Returns:

  • list[str]

    The ids of the submitted jobs.

Source code in femto/fe/atm/_runner.py
def submit_network(
    config: "femto.fe.atm.ATMConfig",
    network: femto.fe.inputs.Network,
    output_dir: pathlib.Path,
    queue_options: femto.fe.utils.queue.SLURMOptions,
    mpi_command: list[str] | None = None,
) -> list[str]:
    """Submits a set of ATM calculations to the SLURM queueing manager.

    Args:
        config: The configuration.
        network: The network of edges to run.
        output_dir: The directory to store any outputs in.
        queue_options: The options to use when submitting the jobs.
        mpi_command: The mpi runner command to use. The default is
            ``"srun --mpi=pmix"``.

    Returns:
        The ids of the submitted jobs.
    """

    mpi_command = mpi_command if mpi_command is not None else ["srun", "--mpi=pmix"]

    output_dir.mkdir(exist_ok=True, parents=True)

    date_str = datetime.datetime.now().strftime("%Y-%m-%d--%H-%M-%S")
    config_path = output_dir / f"config-{date_str}.yaml"
    config_path.write_text(config.model_dump_yaml(sort_keys=False))

    femto_command = [
        "femto",
        "atm",
        f"--config={config_path}",
        "run-workflow",
    ]

    slurm_job_ids = []

    for edge in network.edges:
        edge_dir = output_dir / f"{edge.ligand_1.name}~{edge.ligand_2.name}"

        job_id = femto.fe.utils.queue.submit_slurm_job(
            [
                *mpi_command,
                *femto_command,
                *_create_run_flags(network.receptor, "receptor"),
                *_create_run_flags(edge.ligand_1, "ligand-1"),
                *_create_run_flags(edge.ligand_2, "ligand-2"),
                f"--output-dir={edge_dir}",
                f"--report-dir={edge_dir}",
            ],
            queue_options,
            edge_dir / f"run-{date_str}.out",
        )

        slurm_job_ids.append(job_id)

    return slurm_job_ids