[WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching #431

Olllom · 2019-07-20T05:22:01Z

Description

This pull request contains two features that are largely independent:

An extension of the alchemical facility to work with custom forces, especially NBFIXes coming from OpenMM's PSF parser.
Custom van-der-Waals switching functions, in particular, CHARMM's VSWITCH and VFSWITCH.

If you would like to have one or both in openmmtools, feel free to take it. I am happy to amend it where needed.

1) Alchemical Functions

added methods _alchemically_modify_CustomNonbondedForce and _alchemically_modify_CustomBondForce to AlchemicalFactory. These methods identify Lennard-Jones forces based on their parameters (sigma, epsilon) and tabulated functions (acoeff, bcoeff) names. If LJ forces are identified, they are automatically augmented by a softcore potential. (If you find that check too loose, I could also check for specific energy functions.) By this change, the AlchemicalFactory will handle all LJ forces coming from openmm.app.CharmmPsfFile.createSystem as well as ones that have been augmented by switching functions (see 2) ).
Added functionality to facilitate alchemical manipulation of other custom forces: a function modify_custom_force that can add lambda variables and softcore potentials to an existing energy string, or simply replace the energy string by a user-defined alchemical transformation. AlchemicalState has been modified to pick up user-defined lambda variables lambda_...
Added a testsystem from CHARMM files (small molecule in salt water).

2) Van-der-Waals Switching Functions

Important, for example, for membrane simulations using the C36 force field.

CHARMM-GUI provides a VFSWITCH for OpenMM, but it may not be bad to have a similar implementation in openmmtools, too.
Support for custom, user-defined switching functions.

Todos

Profiling of VSWITCH and VFSWITCH.
Compare VSWITCH and VFSWITCH energies with CHARMM output.
Add documentation and example of modifying custom forces to AlchemicalFactory, AlchemicalState.
Update documentation as needed
Update changelogNotable points that this PR has either accomplished or will accomplish.

Status

Ready to go

jchodera · 2019-07-21T21:13:57Z

@Olllom : Thanks so much for this contribution! At first glance, it looks very much like how I'd implement this scheme.

@andrrizzi : Can you take a look on Monday and give some feedback?

Relatedly, I've been thinking about how we might modularize the choice of softcore functional forms in the future, especially given new efficient forms that have been proposed in addition to the generalized form from Michael Shirts we currently support:

This linear version from Helmut Grubmüller
An exciting improvement from Bert de Groot

We can always perform this generalization later, after this PR is merged.

andrrizzi · 2019-07-21T22:05:29Z

Thank you! Will take a look in the next few days.

I've been thinking about how we might modularize the choice of softcore functional forms

I've also worked on incapsulating some logic into Force objects to simplify the factory. Right now I have the nonbonded forces implemented, but I need to clean up docstrings and tests before pushing it.

…witching

Olllom · 2019-12-07T17:30:30Z

Andrea, Jaime, and I talked about this PR briefly earlier this week. We were wondering if somebody other than @andrrizzi could offer a (preliminary) review. I could also break this into smaller chunks if that helps. The new CharmmSolvated testsystem is straightforward. The switching functions should be easy, too.

The latest commits include

a test with respect to CHARMM energies for the VSWITCH and VFSWITCH LJ switching functions
a test with respect to CHARMM energies for the alchemical endstates
a rebase to the current master

The CHARMM reference calculations are hosted here.

Since we are now mentioning this PR in a paper it would be good if we could move towards merging it into master.

Best,
Andreas

andrrizzi · 2019-12-10T14:48:37Z

Hey Andreas! Thanks for updating the PR with master. If nobody gets to this, I'll spare a couple of hours during the weekend to review it.

andrrizzi

Looks great! Thank you for docs and the tests!

A few general comments:

Could you make the utility functions in alchemy.py private for now? I have a work in progress branch where we're transitioning towards Force objects for alchemy rather than using functional programming, and I think some of this logic will end up in forces.py after that. If we make the private we won't have to break the API when that'll happen.
Could you group the tests in test_alchemy.py into a TestClass similar to other tests? That really helps with testing code maintenance.
Could you take a pass at the tests and add a brief description of what tests are doing in the docstrings? Also for maintenance.

andrrizzi · 2019-12-16T09:18:38Z

openmmtools/alchemy.py

@@ -231,6 +233,9 @@ def from_system(cls, system, *args, **kwargs):
            If specified, the state will search for a modified
            version of the alchemical parameters with the name
            ``parameter_name + '_' + parameters_name_suffix``.
+        custom_lambdas : list of str, optional


I'd prefer not to expose this argument. The reason is simply that dynamically-created classes do not behave nicely in many cases (serialization above all), which we'll have to support if we expose this feature. The encouraged way to extend this is through (static) class extension

class MyAlchemicalState(AlchemicalState): custom_param1 = AlchemicalState._LambdaParameter(...) custom_param2 = AlchemicalState._LambdaParameter(...)

If your application requires the dynamic instantiation, could you wrap these lines in a function in your code? I don't think that a MultiStateSampler would currently be able to resume correctly if you'll pass one of these classes to it though.

andrrizzi · 2019-12-16T09:19:05Z

openmmtools/alchemy.py

+
+        # this is a quick fix.
+        # TODO: extend to multiple alchemical regions
+        assert len(alchemical_regions) == 1


Could you raise a NotImplementedError here rather than an assert explaining this is not currently supported?

andrrizzi · 2019-12-16T09:19:33Z

openmmtools/alchemy.py

+                sigma_string="sigma", bonded_force=True
+            )
+            terms = softcore_energy_string.split(";")
+            softcore_energy_string = ";".join(terms)


Does this have any effect since terms = softcore_energy_string.split(";")?

andrrizzi · 2019-12-16T09:19:52Z

openmmtools/alchemy.py

+
+        # this is a quick fix
+        # TODO: extend to multiple alchemical regions
+        assert len(alchemical_regions) == 1


andrrizzi · 2019-12-16T09:20:14Z

openmmtools/alchemy.py

+            # No fitting modification found
+            return {'': [copy.deepcopy(reference_force)]}
+
+        if is_custom_lj_force(reference_force) or is_nbfix_force(reference_force):


I think this if-clause is always True, so you could remove one level of indentation here.

andrrizzi · 2019-12-16T09:21:24Z

openmmtools/alchemy.py

+        passed as scale_string or effective_radius_string.
+    kwargs:
+        Any additional global parameters that the force depends on. Parameters that begin with "lambda_" will
+        be added to the alchemical state upon creation.


I don't think this happens automatically currently. Maybe add that "will have to be added"?

andrrizzi · 2019-12-16T09:21:46Z

openmmtools/alchemy.py

+            force.addGlobalParameter(arg, kwargs[arg])
+            if not arg.startswith("lambda_"):
+                if hasattr(alchemical_region, arg) and getattr(alchemical_region, arg) != kwargs[arg]:
+                        warnings.warn("Argument {} in alchemical region was () but is getting overwritten "


I think the () should be {}.

andrrizzi · 2019-12-16T09:22:55Z

openmmtools/forcefactories.py

+    return False
+
+
+def get_forces_that_define_vdw(system):


Maybe rename this find_vdw_forces for consistency with forces.find_forces?

It looks like this function encapsulates some logic that is also exposed in alchemy.is_nbfix_force and is_custom_lj_force. Does it make sense to remove the redundant logic?

andrrizzi · 2019-12-16T09:24:20Z

openmmtools/tests/test_alchemy.py

+    assert len(nbfix_forces) == 1
+
+
+def endstate_test(use_vswitch=False, pme_tol=1e-5, platform=openmm.Platform.getPlatformByName("Reference")):


We usually use check_endstates for this kind of non-test functions.

andrrizzi · 2019-12-16T09:25:31Z

openmmtools/tests/test_alchemy.py

+
+    # utility functions
+    def calc_energy(testsystem, platform=platform):
+        dummy_integrator = openmm.VerletIntegrator(1.0*unit.femtosecond)


Could you implement this using the existing test_alchemy.compute_energy function in here?

Olllom added 8 commits February 5, 2019 22:09

Added custom Lennard-Jones switching functions

76bde74

Clean up implementation and documentation of forcefactories

883c193

Added a testsystem for solvation using charmm input files

44dc396

CharmmSolvated testsystem optionally without vdw/charges

51ea8fe

alchemical softcore for custom forces

259caf5

bugfix in charmm testsystem and vdw factory

30e1ecb

Modifying custom forces.

b6aecf8

Merge branch 'master' into vdw_switching

606c339

Olllom changed the title ~~WIP Alchemical Custom Forces + Custom Lennard-Jones Switching~~ [WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching Jul 20, 2019

Olllom added 2 commits September 23, 2019 15:09

Merge branch 'master' of github.com:choderalab/openmmtools into vdw_s…

c3af895

…witching

merge master into vdw_switching

a163240

Olllom mentioned this pull request Oct 29, 2019

ThermodynamicState support for more barostats #437

Merged

5 tasks

Olllom added 5 commits December 6, 2019 21:09

added energy comparison to charmm for lj switch

bbe0e74

bugfix

7d77d72

CHARMM reference for uncharged and annihilated energies

6cc4192

rename get_energy to calculate_energy in LJPairBox

b9ec52e

Merge branch 'master' into vdw_switching

af68f39

andrrizzi reviewed Dec 16, 2019

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching #431

[WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching #431

Olllom commented Jul 20, 2019 •

edited

jchodera commented Jul 21, 2019

andrrizzi commented Jul 21, 2019

Olllom commented Dec 7, 2019

andrrizzi commented Dec 10, 2019

andrrizzi left a comment

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

andrrizzi Dec 16, 2019

		assert len(nbfix_forces) == 1


		def endstate_test(use_vswitch=False, pme_tol=1e-5, platform=openmm.Platform.getPlatformByName("Reference")):

[WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching #431

Are you sure you want to change the base?

[WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching #431

Conversation

Olllom commented Jul 20, 2019 • edited

Description

1) Alchemical Functions

2) Van-der-Waals Switching Functions

Todos

Status

jchodera commented Jul 21, 2019

andrrizzi commented Jul 21, 2019

Olllom commented Dec 7, 2019

andrrizzi commented Dec 10, 2019

andrrizzi left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Olllom commented Jul 20, 2019 •

edited