Charmm36 nonbonded forces with Alchemy

[AJZein]

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)

[Lnaden]

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 or type2 is in the AlchemicalRegion, 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.

[AJZein]

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.

[Lnaden]

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 and 4eσ^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.

[Lnaden]

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.

[Lnaden]

Last thought for now: Also have to think about how to properly softcore a tabulated function where the 4e term is folded in...

[jchodera]

@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 and lambda_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 with scaling*({energy_expression})
• To ensure we make an appropriate softcore substitution, we replace r with effective separation r_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 nor epsilon, we would pre-multiply it by lambda_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.

[Lnaden]

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.

[AJZein]

@jchodera No problem, I can understand that covid has put a lot of people in difficult situations.