Giter Site home page Giter Site logo

alchemy's People

Contributors

andrrizzi avatar jchodera avatar lnaden avatar mikemhenry avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

alchemy's Issues

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 &'

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

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?

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.

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.?

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.

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.

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.

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.

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));

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.

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    ๐Ÿ–– Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. ๐Ÿ“Š๐Ÿ“ˆ๐ŸŽ‰

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google โค๏ธ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.