choderalab / alchemy Goto Github PK
View Code? Open in Web Editor NEWAlchemical tools for OpenMM
Alchemical tools for OpenMM
We need to add a note to the README.md
directing users to openmmtools.alchemy
.
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 &'
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
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 1System
object for system 2How 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:
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)alpha
(Lennard-Jones) and beta
(Coulomb) softcore parameterslambda_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?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?The current concept is to follow closely the create_relative_transformation code from examol
.
All Force
s 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 to
Custom*Forces. Force constants (
K) and equilibrium values (
r0,
theta0) are linearly interpolated between their initial and final values according to
lambda_bondsand
lambda_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
alpha
and beta
controlling this.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.
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?
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
.
Looks like openmm/openmm#1311 may have our GBSA support? Looking into this...
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.
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)
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.
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
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.?
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
The tests are too slow to complete on travis; it times out instead.
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.
I can't quite figure out this error:
https://travis-ci.org/choderalab/alchemy/jobs/81441089#L580-L584
This should give the same result as alchemical_bonds
, alchemical_angles
, and alchemical_torsions
= None
.
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.
We should investigate whether reworking some of the alchemical states can speed up simulations.
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
.
copy.deepcopy
for OpenMM forces is broken right now:
openmm/openmm#1190
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
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?
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.
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.
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.
Let's be sure to add the MIT License here.
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));
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.
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.
This name change is less likely to cause accidental name collisions.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.