openmmtools
openmmtools copied to clipboard
Charmm36 nonbonded forces with Alchemy
Hello,
I am trying to perform alchemy using the charmm36 forcefield but I've noticed that there are remaining forces when lambda is set to 0. In the LysozymeImplicit test system for example, the pxylene_atoms wont separate from its binding pocket when lambda is set to 0 in charmm36, but it will when using the test system's forcefield. I've added a minimal example for both cases at the bottom.
I've tried looking into the problem and I believe it might be due to the LJ Parameters being implemented in a CustomNonbondedForce when using charmm36, which allows for adding nonbonded fix (NBFIX) as a lookup table. However I don't know how to get alchemy to work with this.
Thanks for any help, Zein
Working example
from openmmtools import alchemy, testsystems, integrators
import simtk.openmm.app as app
from sys import stdout
# Create the reference OpenMM System that will be alchemically modified.
lysozyme_pxylene = testsystems.LysozymeImplicit()
t4_system = lysozyme_pxylene.system
# Define the region of the System to be alchemically modified.
pxylene_atoms = lysozyme_pxylene.mdtraj_topology.select('resname TMP')
alchemical_region = alchemy.AlchemicalRegion(alchemical_atoms=pxylene_atoms)
# Set up simulation
factory = alchemy.AbsoluteAlchemicalFactory()
alchemical_system = factory.create_alchemical_system(t4_system, alchemical_region)
integrator = integrators.LangevinIntegrator()
simulation = app.Simulation(lysozyme_pxylene.topology, alchemical_system, integrator)
simulation.context.setPositions(lysozyme_pxylene.positions)
# Add reporting
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
app.PDBFile.writeFile(simulation.topology, lysozyme_pxylene.positions, open(f'working.pdb', 'w'))
simulation.reporters.append(app.DCDReporter(f'working.dcd', 1000))
# Remove interactions by setting lambda=0
alchemical_state = alchemy.AlchemicalState.from_system(alchemical_system)
alchemical_state.lambda_electrostatics = 0
alchemical_state.lambda_sterics = 0
alchemical_state.apply_to_context(simulation.context)
simulation.step(20000)
Failing example with charmm36
from openmmtools import alchemy, testsystems, integrators
import simtk.openmm.app as app
from sys import stdout
# Create the reference OpenMM System that will be alchemically modified.
lysozyme_pxylene = testsystems.LysozymeImplicit()
forcefield = app.ForceField('charmm36.xml')
t4_system = forcefield.createSystem(lysozyme_pxylene.topology, constraints=app.HBonds)
# Define the region of the System to be alchemically modified.
pxylene_atoms = lysozyme_pxylene.mdtraj_topology.select('resname TMP')
alchemical_region = alchemy.AlchemicalRegion(alchemical_atoms=pxylene_atoms)
# Set up simulation
factory = alchemy.AbsoluteAlchemicalFactory()
alchemical_system = factory.create_alchemical_system(t4_system, alchemical_region)
integrator = integrators.LangevinIntegrator()
simulation = app.Simulation(lysozyme_pxylene.topology, alchemical_system, integrator)
simulation.context.setPositions(lysozyme_pxylene.positions)
# Add reporting
simulation.reporters.append(app.StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True))
app.PDBFile.writeFile(simulation.topology, lysozyme_pxylene.positions, open(f'not_working.pdb', 'w'))
simulation.reporters.append(app.DCDReporter(f'not_working.dcd', 1000))
# Remove interactions by setting lambda=0
alchemical_state = alchemy.AlchemicalState.from_system(alchemical_system)
alchemical_state.lambda_electrostatics = 0
alchemical_state.lambda_sterics = 0
alchemical_state.apply_to_context(simulation.context)
simulation.step(20000)
I've done a bit of digging and I think I can better describe the source of the problem.
I think you are correct: This is a bug due to OpenMMTool's alchemical factory never being taught to recognize OpenMM's implementation of the CHARMM forcefield.
The t4_system
object you create from the charm36.xml
file and forcefield
object adds in a native CustomNonbondedForce
to the OpenMM System
. This force is a pure LJ force (no electrostatics) computed from two Discrete2DFunction
's:
'acoef(type1, type2)/r^12 - bcoef(type1, type2)/r^6;'
where acoef
and bcoef
are the tabulated functions based on the paired interaction. OpenMM's native forcefield.py
code makes this here.
When OpenMMTools creates its alchemical forces, it goes through the native forces it expects OpenMM to make, and then modifies parameters based on the AlchemicalRegion
and adds forces as needed. So the fix here is for OpenMMTools to recognize this corner case, then compensate for it accordingly.
I can probably write up a code you can use in the short term to get your system working, but I don't I personally will have time to implement a proper fix in the code base.
The short fix would be to modify the energy expression of the native CustomNonbondedForce
object to recognize if the particle is in the AlchemicalRegion
and then null that interaction out. Caveats:
- This might not give the correct energies at lambda=1 without me further studying how the particles in the
NonbondedForce
object in that system are parameterized. - You could rework that whole function to properly be softcore, but that would be much more work.
- Instead of zeroing out the interaction and changing the energy expression string, you could reparameterize the discrete function to have an outlier which drops to 0 if
type
ortype2
is in theAlchemicalRegion
, much more risky though, might have some odd derivatives if you're too close to the normal spline the function defines, and that might be more work than the proposed option.
Thanks for looking into this so quickly. I'll try explore some other options but in the end I may have to just rewrite the whole function.
I wanted to touch on this issue a bit more. I've done some looking at the forces which come out of the CHARMM forcefield.
It looks like the Lennard-Jones interaction is purely handled by the CustomNonbondedForce
object and none of the NonbondedForce
. The NonbondedForce
has every epsilon
set to 0, so the LJ interactions are never handled by that force. As such, most of the changes the AlchemicalFactory wont actually do anything. I think the fix would be to add another function to detect the form of the CHARMM LJ implementation in this CustomNonbondedForce
and modify accordingly.
That said, I don't know what the fix should be:
- The CHARMM-implementation
CustomNonbondedForce
has to be directly modified to handle alchemical particles, checking each pair of particles to either be alchemical or not so the non-alchemical/non-alchemical. Because its a tabulated function, zeroing out the interaction is harder. - The Alchemical/alchemical interactions have to be accounted for somehow, also through the tabulated function. Not too hard, just copy the 2D functions and use it for the
4eσ^12
and4eσ^6
terms, make them softcore as normal. - The Alchemical/non-alchemical interactions have to be accounted for in a similar way.
Although the CHARMM forcefield in OpenMM has this problem, I think any forcefield which calls the LennardJonesGenerator
class from forcefield.py
will have this problem.
The proper fix to this is complicated by a few other factors from alchemy.py
's design:
- Alchemy (it) assumes the nonbonded parameters are all contained in the Nonbonded Force
- It assumes the mixing rules are represented by simple Lorentz-Berthelot and NOT by tabulated function.
- It compiles all the alchemical custom forces (i.e. electrostatics and LJ) from basic Nonbonded force
- It assumes no other forces are contributing to the nonbonded interactions from the base system.
These assumptions make it simpler for the alchemical force factory to flow down the list and cast each native force through a single dispatch, getting out all the needed forces in one pass. However, to properly include these tabulated functions, all the Nonbonded forces would have to be compiled together and passed to the function which makes the alchemical forces, or we can treat both native forces as separate, which will result in more CustomNonbondedForces
Adding more CustomNonbondedForces will make the simulation a bit slower, but without doing an overhaul of the code base to handle multiple forces providing the native interactions, i don't know how this is avoidable.
Last thought for now: Also have to think about how to properly softcore a tabulated function where the 4e
term is folded in...
@ajzein: Apologies for the long delay in addressing this---we're shorthanded on developers during COVID, but should have a new developer aboard to start working through the backlog of issues later in June.
@Lnaden : Thanks for diagnosing this! #511 is a good temporary safety improvement.
We should be able to implement a simple solution to handle CustomNonbondedForce
terms generally, following a similar strategy to how we deal with alchemically modifying CustomGBForce
.
We could alter the CustomNonbondedForce
using the following rules:
-
lambda_sterics
andlambda_electrostatics
are added as global parameters -
is_alchemical
is added as a per-particle parameter. All atoms in the alchemical group have this parameter set to 1, else 0 - We define
is_alchemical = step(is_alchemical1 + is_alchemical2)
scaling = lambda_X * is_alchemical + (1-is_alchemical);
where lambda_X
is lambda_electrostatics
if the expression contains charge
, or lambda_sterics
if it contains sigma
and epsilon
. If neither, we default to lambda_sterics
since this is usually turned off last.
- We replace
energy_expression
withscaling*({energy_expression})
- To ensure we make an appropriate softcore substitution, we replace
r
with effective separationr_eff
defined by
r_eff = sigma*((softcore_beta*(1-scaling)^softcore_e + (r/sigma)^softcore_f))^(1/softcore_f)
where softcore_X
are appropriate softcore class parameters that can be tuned.
More recently, we've been experimenting with 4th dimensional lifting, which has a simpler form:
r_eff = (r^2 + ((1-scaling)*softcore_fraction*r_cutoff)^2)^(1/2)
where we are just "lifting" the interaction out into the 4th dimension a fraction of the cutoff as lambda
approaches 0
- If one of the parameters is named
epsilon
, we would substitute it with(lambda_sterics*epsilon)
- If it neither contains
charge
norepsilon
, we would pre-multiply it bylambda_sterics
to make sure it ends up in some group
This should safely work for all kinds of CustomNonbondedForce
implementations, including the lookup-table based LJ and reaction field electrostatics. The only problematic ones are likely to be the ones that mix sterics and electrostatics into a single CustomNonbondedForce
, but that might still be "safe" if we treat them all with the sterics.
Future fancier versions of this could make use of multiple CustomNonbondedForce
terms and particle groups to split out the alchemical interactions from the non-alchemical for better performance.
I like the solution for the alchemical side of things in the energy expression. I do wonder if we need to somehow incorporate the existing floats/constant of alpha, beta, a, b, c into this to be consistent with how the standard NonbondedForce
is cast.
The other engineering problem is how to detect this and not have an overlap of force objects generated from the Factory. Right now the factory detects the NonbondedForce
and splits out all the CustomNonbondedForce
objects to handle electrostatics and Lennard-Jones for all alchemical/nonalchemical permutations. The current logic would not be able to correctly and reliably handle the case where the electrostatics are contained in the NonbondedForce
and the Lennard-Jones are in this CustomNonbondedForce
. This is a solvable problem, but will take some engineering.
@jchodera No problem, I can understand that covid has put a lot of people in difficult situations.