Implementing Heating and Restraint Ladder

[beitom]

Hello Im looking for advice on implementing a custom heating and restraint ladder before the production MD simulation. The specific details for reference are:

For each edge simulation the hybrid system at λ = 0.5 using the protocol below:

• Energy minimization (50k steps).
• Heats 0→300 K over 200ps.
• 200 ps heavy-atom restraints (10 kcal /mol A harmonic restraint).
• 200 ps backbone restraints (10 kcal /mol A harmonic restraint).
• 200 ps unrestrained.

Does openfe support adding this kind of multistage protocol directly or would I be better off running to with openmm and configuring openfe to use the openmm results as the intial cooridfinates. For the ladder how would I generate the hybrid system λ = 0.5 to feed easily into openmm.

[IAlibay]

@beitom the brief answer is that unfortunately OpenFE doesn't support extended equilibration workflows at the moment.

It is on our to-do list, but not something we've been able to prioritise since the RBFE protocol usually works well if you come with a system somewhat ready for MD.

Re: equilibrating elsewhere and then bringing it to OpenFE. That is feasible, but for now it does require removing non-binding site waters (you can use a tool like MDAnalysis to remove waters more than a given cutoff from your ligand). We are in the process of adding support for ingesting a fully pre-solvated box (with box vectors, etc...), although it might not make it this side of the year.

[beitom]

Thank you @IAlibay. Do you have any suggestions for how to use openfe to prepare hybrid system but feed that into my own openmm pipeline. Can openfe create system.xml like files to load into openmm. Otherwise I'd be worried that if i created the hybrid system file using another tool like openmm alchemy it may not match exactly what openfe is doing and my equilibration stage could be off.

Thanks

[beitom]

Hello @IAlibay hoping for your insight on this being a valid way to extract the openmm hybrid system to run equiliriation

# assume transformation is defined in scope as an openfe.Transformation object that uses a protocol with 
# lambda_settings.lambda_windows = 1, simulation_settings.n_replicas = 1

protocol_dag = protocol.create(
        stateA=transformation.stateA,
        stateB=transformation.stateB,
        mapping=transformation.mapping,
)

protocol_unit = protocol_dag.protocol_units[0]

dry_run_output = protocol_unit.run(dry=True, verbose=False)  # type: ignore

sampler = dry_run_output["debug"]["sampler"]
factory = sampler._factory
hybrid_system: openmm.System = factory.hybrid_system
hybrid_topology: Topology = factory.hybrid_topology
hybrid_positions: openmm.unit.Quantity = factory.hybrid_positions

Then I build an openmm simulation for $\lambda=0.5$

def get_available_parameters(system):
    """Get all global parameters available in the system."""
    available_parameters = set()
    for i in range(system.getNumForces()):
        force = system.getForce(i)
        if hasattr(force, "getNumGlobalParameters"):
            for j in range(force.getNumGlobalParameters()):
                param_name = force.getGlobalParameterName(j)
                available_parameters.add(param_name)
    return available_parameters


integrator = VerletIntegrator(1.0)
simulation = Simulation(hybrid_topology, hybrid_system, integrator)

simulation.context.setPositions(hybrid_positions)

lambda_params = [
    param
    for param in get_available_parameters(hybrid_system)
    if param.startswith("lambda") and "_core" not in param
]

for param in lambda_params:
    simulation.context.setParameter(param, 0.5)

# example: minimize energy

LocalEnergyMinimizer.minimize(
            simulation.context
        )