choderalab/alchemy

Alchemical tools for OpenMM

alchemy's Issues

alchemy-dev conda package is not being pushed

Assign MIT License

Let's be sure to add the MIT License here.

Eliminate ligand_atoms and receptor_atoms in favor of a single alchemical_atoms

We may later want the ability to have multiple alchemical atomsets, supporting a more complex set of alchemical parameters such as lambda_sterics_[name], etc.

Add support for atom-by-atom alchemical intermediate definition

Copied from issue in YANK

We can easily add an alternative alchemical intermediate generator based on this scheme:
http://dx.doi.org/10.1002/jcc.21829

Enable support for all OpenMM GBSA models

Looks like openmm/openmm#1311 may have our GBSA support? Looking into this...

Correct build not being uploaded to alchemy-dev conda package

I don't know how this is happening, but it seems like alchemy-dev on conda contains a 1.1.2 build, not the latest github build.

https://anaconda.org/omnia/alchemy-dev/files

Different Alchemical Transformations of Multiple Atom Groups?

Dear All,

In "AbsoluteAlchemicalFactory(object).init(...), I found the following ToDo:
"* Can we specify multiple groups of alchemically-modified atoms that have different alchemical parameters associated with them, like lambda_sterics_[groupname]?"

It would be very nice to have this feature in the near future to address important questions regarding multiple ligand-binding sites, mutations, etc.
I am sure you are very busy developing other parts of the code right now, but how would you envision the implementation of such feature? (I am also asking in hopes of contributing myself after having understood your code a little better).

I could imagine this to be implemented by generalizing the concept of "ligand_atoms" to two separate arrays of atom index sets ("array_of_deleted_index_sets" and "array_of_inserted_index_sets") and optional corresponding arrays of lambda dependencies (each defined as a list of string, e.g., ["lambda_sterics = lambda*2", "lambda_electrostatics = (lambda>0.3)(lambda-0.3)/0.7"]), so that users can
(i) have a simple lambda-protocol for multiple alchemical transformations (with .addExclusion(...) automatically set between all atom indices in the i-th set of "array_of_deleted_index_sets" and all atoms in the i-th set of "array_of_inserted_index_sets")
(ii) if the number of force groups <= 32, be able to manually define different lambda dependencies for each index set.

Please let me know about your thoughts,
Thank you, Best,
Sebastian

Rename context parameters from `lambda_[name]` to `alchemy_[name]` in release 2.0`

This name change is less likely to cause accidental name collisions.

There may be something wrong with alchemical angles

cc choderalab/perses#207

Add a note that this repo has been deprecated

We need to add a note to the README.md directing users to openmmtools.alchemy.

Merge alchemy and init?

Is there a specific reason why there is both an alchemy.py and __init__.py?

__init__.py just contains from .alchemy import * , so alchemy.py seems rather redundant.

Add API to allow bonds, angles, and torsions to be alchemically softened

See choderalab/perses#9

Add sanity checks to alchemical factory

Given the trouble we're having with large systems causing alchemically-modified systems to give NaN energies, it will be useful to add a number of sanity checks that help debug which force component is causing the NaN.

We can add an optional argument to AlchemicalFactor.__init__() to compute the total potential energy for a given set of positions and if NaN is detected, compute the energy (on the Reference platform) by force component to pin down the culprit.

Moved from yank repo to here.

Change ligand_atoms and receptor_atoms to just alchemical_atoms

Don't use openmm-dev with alchemy

copy.deepcopy for OpenMM forces is broken right now:
openmm/openmm#1190

Not sure why python 3.x is not building

I can't quite figure out this error:
https://travis-ci.org/choderalab/alchemy/jobs/81441089#L580-L584

Cache simulations generated by test_overlap

We can save a great deal of testing time if we enable travis caching and store the trajectory snapshots generated by check_overlap() in a $HOME/.cache/alchemy/ directory.

More information on caching can be found here:
https://docs.travis-ci.com/user/caching/#before_cache-phase

Improve PME overlap by computing reciprocal-space exceptions

We can improve the PME overlap by subtracting out the exceptions computed in reciprocal space.

The relevant OpenMM code that implements this is located here.

To implement this, we would add a CustomBondForce that computes all exceptions and exclusions as bond terms using the direct-space Ewald energy:

-(lambda_electrostatics^softcore_d)*ONE_4PI_EPS0*chargeprod*erf(alpha_ewald*reff_electrostatics)/reff_electrostatics

(note the minus sign, since we are subtracting these).

chargeprod would be computed as charge1*charge2, and would NOT match the chargeprod used in the exception/exclusion. This is because we are subtracting off the PME reciprocal space contribution, which uses the product of charges.

We would iterate over all exceptions and exclusions in the original NonbondedForce.

I think we can just use r instead of bothering to compute softcore reff_electrostatics.

Proposed API for DualTopologyFactory

We need the ability to create OpenMM System objects with context parameters that allow us to perform relative free energy calculations between two System endpoints.

We presume that the user starts with:

• System object for system 1
• System object for system 2
• a mapping between atoms in systems 1 and 2, where we also require that the positions of atoms listed in this map are identical in both sets of positions

Proposed API

How does this look for an API:

from alchemy.relative import DualTopologyFactory
# Create a factory that generates alchemical intermediates.
factory = DualTopologyFactory(system1, system2, atom_mapping)
# Retrieve the alchemically-modified system.
alchemical_system = factory.createPerturbedSystem()

Here atom_mapping is a dict that maps atoms in system1 to system2.

We can add some other optional arguments to the factory:

• a context_parameter_prefix that allows us to specify the prefix prepended to alchemical context parameters (with the default of lambda, e.g. lambda_sterics, lambda_electrostatics); this could avoid parameter namespace collisions if other alchemical components are in use (e.g. alchemical softening of a loop or residues near a binding site through an orthogonal parameter)
• additional optional factory arguments could control the alchemical softcore potential, such as the softness parameters alpha (Lennard-Jones) and beta (Coulomb) softcore parameters

Some questions:

• Do we want to allow the user full control over all of the individual alchemical context parameters (lambda_sterics, lambda_electrostatics, lambda_torsion, lambda_bonds, lambda_angles), or do we want to have a single lambda parameter and let the user specify a set of functions to describe how these various context parameters are slaved to the single controlling function?
• Some of the additional parameters controlling the softcore potential, such as softness parameters (alpha, beta) and various exponents, could either be hard-coded at System creation time, could be left as alchemical parameters the user could change, or could have functions specified that describe how they are slaved to a master lambda parameter. Which one would be more useful?

Implementation details

The current concept is to follow closely the create_relative_transformation code from examol.

All Forces must appear in the same order in each system. Failure to do so will trigger an exception.

Each Force term in is handled differently:

HarmonicBondForce and HarmonicAngleForce

Bonds and angles present in both system1 and system2are converted toCustom*Forces. Force constants (K) and equilibrium values (r0,theta0) are linearly interpolated between their initial and final values according tolambda_bondsandlambda_angles`.

PeriodicTorsionForce

Because the number and periodicity of Fourier terms can differ between system1 and system2, it seems most appropriate to have lambda_torsions produce a linear interpolation of the energies of the two systems.

NonbondedForce

• Atoms appearing only in one system but not the other will end up as uncharged, noninteracting dummy atoms in those systems, using softcore electrostatics and sterics to create/delete those atoms.
• Any atoms shared between both systems but differing in Lennard-Jones parameters or charges will have both sterics and electrostatics use an intermediate softcore form, with softness parameters alpha and beta controlling this.
• Any atoms shared between both systems where both Lennard-Jones parameters and charges match between both systems will be treated by normal forces and will not be softcore.
• We are unable to support the reciprocal-space component of PME, requiring us to zero out the charges for any atoms whose charges are . This means that atoms being created, destroyed, or undergoing mutations will have their charges zeroed out in the PME calculation.

GBSAOBCForce

This will be converted to a CustomGBForce implementation.

CustomGBForce

This will work only of the energy function and parameter names are the same; otherwise, an Exception will be thrown.
Parameters will be interpreted between the two systems.
Exclusions will be combined from both systems.
Exclusions in the shared components of the system must match.

Other forces

Right now, we will throw an Exception, but it is possible we will want to treat other forces (MonteCarloBarostat, COMMotionRemover) by checking that all parameters are the same and adding it to the resulting alchemically-modified system

Anything I have forgotten?

Reduce test time on travis

The tests are too slow to complete on travis; it times out instead.

python 3.x issues

I've encountered a failure for the Abl:imatinib explicit solvent YANK test system in Python 3.4:

Traceback (most recent call last):
  File "soften-ligand.py", line 77, in <module>
    factory = AbsoluteAlchemicalFactory(reference_system, ligand_atoms=alchemical_atoms, annihilate_electrostatics=True, annihilate_sterics=False)
  File "/cbio/jclab/home/chodera/miniconda3/envs/sams/lib/python3.4/site-packages/alchemy/alchemy.py", line 356, in __init__
    [self.alchemically_modified_system, self.force_labels] = self._createAlchemicallyModifiedSystem(self.reference_system)
  File "/cbio/jclab/home/chodera/miniconda3/envs/sams/lib/python3.4/site-packages/alchemy/alchemy.py", line 1480, in _createAlchemicallyModifiedSystem
    self._alchemicallyModifyNonbondedForce(system, reference_force, force_labels)
  File "/cbio/jclab/home/chodera/miniconda3/envs/sams/lib/python3.4/site-packages/alchemy/alchemy.py", line 1197, in _alchemicallyModifyNonbondedForce
    na_sterics_custom_nonbonded_force.addInteractionGroup(nonalchemical_atomset, alchemical_atomset)
  File "/cbio/jclab/home/chodera/miniconda3/envs/sams/lib/python3.4/site-packages/simtk/openmm/openmm.py", line 9950, in addInteractionGroup
    return _openmm.CustomNonbondedForce_addInteractionGroup(self, set1, set2)
TypeError: in method 'CustomNonbondedForce_addInteractionGroup', argument 3 of type 'std::set< int,std::less< int >,std::allocator< int > > const &'

Unpin conda-build from 2.0.1 to fix travis

Currently, devtools/travis-ci/install.sh pins conda-build==2.0.1, which causes travis to fail.

@Lnaden : Can you create a PR to un-pin?

Allow alchemical_bonds, angles, and torsions to be specified as False

This should give the same result as alchemical_bonds, alchemical_angles, and alchemical_torsions = None.

Fix reaction field alchemical treatment

By default, OpenMM uses a form of reaction field electrostatics that adds an additive term that depends on the charges to ensure charge-charge interaction energies go to zero at the cutoff. Unfortunately, the charge-dependent additive shift (c_rf) means that computing the charging free energy of a dipole in a cavity with an exterior dielectric has an arbitrary constant added to it that cannot be eliminated.

Until OpenMM permits us to set c_rf = 0, we will need to implement reaction field ourselves without this c_rf in order to permit alchemical free energy calculations with CutoffPeriodic.

Direct-space PME contribution incorrect if alpha_ewald = 0

If the reference System object does not manually specify alpha_ewald for PME, it appears that the direct-space PME contribution is incorrect. The logic in this block does not work as intended, instead skipping the if alpha_ewald == 0.0 clause (since alpha_ewald has units of 1/nm) and using a value of 0 in the calculation of the electrostatic energy expression.

If alpha_ewald is not manually set, it will be determined automatically from the PME tolerance upon Context creation. This is thankfully standardized in NonbondedForceImpl:

    double tol = force.getEwaldErrorTolerance();
    alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));

Bond force expects 3 particles in bond

Looks like the AbsoluteAlchemicalFactory tries to get 3 particles out of a HarmonicBondForce term. I can fix this unless someone else is itching to do it.

Fix alchemical pathway

Multiple users have noted the current AbsoluteAlchemicalFactory has very poor performance in terms of easily bridging lambda = 0 and lambda = 1 with a pathway of minimal thermodynamic length. I'm working on investigating this problem, but it seems that the overlap is very poor near lambda = 1. We may need to replace this with another functional form, or at least improve the default parameters used to define the functional form.

Citeable releases

Would it be possible to add doi links to alchemy releases using something like zenodo?
https://zenodo.org/

I think it'll be important and useful in any papers that will depend on alchemy.

Proposed API changes for alchemy 2.0

In alchemy 2.0, I would like to support the definition of multiple alchemical regions in a single system, all accessible via independent OpenMM context parameters.

We would need an API to support this new feature. There are a few schemes I can propose, none of which are backwards-compatible, necessitating a major version number change.

Current API:

The alchemical region is specified in the constructor:

factory = AbsoluteAlchemicalFactory(system, ligand_atoms=[0,1], alchemical_torsions=True, alchemical_angles=True, annihilate_sterics=True, annihilate_electrostatics=True)

Scheme 1 : Listomania

Everything becomes lists!

factory = AbsoluteAlchemicalFactory(system, ligand_atoms=[[0,1], [1,2,3]], alchemical_torsions=[False,False], alchemical_angles=[True,False], annihilate_sterics=[False,True], annihilate_electrostatics=[True,True])

Kinda clunky.

Scheme 2 : AlchemicalRegion class

We create a new object do describe the alchemical region, modeled after the old interface, and make a dict or list of that:

region1 = AlchemicalRegion(alchemical_atoms=[0,1], alchemical_torsions=True, alchemical_angles=True, annihilate_sterics=True, annihilate_electrostatics=True) 
region2 = AlchemicalRegion(alchemical_atoms=[1,2,3], alchemical_torsions=False, alchemical_angles=False, annihilate_sterics=True, annihilate_electrostatics=False) 
regions = { 'region1' : region1, 'region2' : region2 }
factory = AbsoluteAlchemicalFactory(system, alchemical_regions=regions)

The region names are then used in the alchemical parameters

Scheme 3 : API extensions

We create more API methods for AbsoluteAlchemicalFactory to define the regions:

factory = AbsoluteAlchemicalFactory(system)
factory.addAlchemicalRegion(name='region1', alchemical_atoms=[0,1], alchemical_torsions=True, alchemical_angles=True, annihilate_sterics=True, annihilate_electrostatics=True)
factory.addAlchemicalRegion(name='region2', alchemical_atoms=[1,2,3], alchemical_torsions=False, alchemical_angles=False, annihilate_sterics=True, annihilate_electrostatics=False) 

This has much of the same benefit of Scheme 2 without requiring an additional class.

Any comments/suggestions from @andrrizzi, @pgrinaway, etc.?

Switch to versioneer

Look into performance regression introduced by exceptions overhaul in 1.2.1

We should investigate whether reworking some of the alchemical states can speed up simulations.

cc: choderalab/yank#577 (comment)