Serial simulations using same alchemical states as YANK

[daveminh]

Could you post some quick demonstration code on how to set up an OpenMM system using a specific value of the alchemical coupling parameter?

[jchodera]

We'll also want to point you toward how to use the new ThermodynamicState way of updating the alchemical parameters.

[andrrizzi]

Unless I misunderstood your use case, you should be able to do this using openmmtools only. This is an example.

import openmmtools as mmtools
from simtk import openmm, unit

# Create an alchemically-modified system from a reference OpenMM system.
test_system = mmtools.testsystems.AlanineDipeptideExplicit()
openmm_system, positions = test_system.system, test_system.positions
alchemical_factory = mmtools.alchemy.AlchemicalFactory()
alchemical_region = mmtools.alchemy.AlchemicalRegion(alchemical_atoms=range(22))
alchemical_system = alchemical_factory.create_alchemical_system(reference_system=openmm_system,
                                                                alchemical_regions=alchemical_region)

# You can simply use this system as a normal OpenMM system and
# modify the context parameters.
# NOTE: this won't work when we'll implement scaled charges.
integrator = openmm.VerletIntegrator(1.0*unit.femtosecond)
context = openmm.Context(alchemical_system, integrator)
context.setParameter('lambda_sterics', 0.5)
context.setParameter('lambda_electrostatics', 0.5)
integrator.step(1000)

# Or you can use the State classes that allow you control
# the thermodynamic state of a System more generally.
thermodynamic_state = mmtools.states.ThermodynamicState(system=alchemical_system,
                                                        temperature=300*unit.kelvin,
                                                        pressure=1.0*unit.atmosphere)
alchemical_state = mmtools.alchemy.AlchemicalState.from_system(alchemical_system)
compound_state = mmtools.states.CompoundThermodynamicState(thermodynamic_state=thermodynamic_state,
                                                           composable_states=[alchemical_state])
compound_state.lambda_sterics = 0.5
compound_state.lambda_electrostatics = 0.5
integrator = openmm.VerletIntegrator(1.0*unit.femtosecond)
context = compound_state.create_context(integrator)
integrator.step(1000)

# You can use AlchemicalFunction to enslave multiple
# lambda parameters to a single variable.
compound_state.set_alchemical_variable('lambda', 0.5)
compound_state.lambda_sterics = mmtools.alchemy.AlchemicalFunction('lambda')
compound_state.lambda_electrostatics = mmtools.alchemy.AlchemicalFunction('lambda')
compound_state.set_alchemical_variable('lambda', 0.0)  # This turns off both VdW and charges.

Note that we plan to provide relatively soon an implementation that directly scales the charges of the NonbondedForce instead of doing it through the global parameter lambda_electrostatics. This will have the advantage of including the reciprocal space in the alchemical states, but the first method (i.e. setting the context parameter) won't work with this.

[andrrizzi]

Closing this to keep the tracker clean. Feel free to reopen in case you still have questions!