Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 15 commits into
base: main
Choose a base branch
from

Conversation

Olllom
Copy link
Contributor

@Olllom Olllom commented Jul 20, 2019

Description

This pull request contains two features that are largely independent:

  1. An extension of the alchemical facility to work with custom forces, especially NBFIXes coming from OpenMM's PSF parser.
  2. 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

@Olllom Olllom changed the title WIP Alchemical Custom Forces + Custom Lennard-Jones Switching [WIP] Alchemical Custom Forces + Custom Lennard-Jones Switching Jul 20, 2019
@jchodera
Copy link
Member

@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:

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

@andrrizzi
Copy link
Contributor

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.

@Olllom
Copy link
Contributor Author

Olllom commented Dec 7, 2019

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

  1. a test with respect to CHARMM energies for the VSWITCH and VFSWITCH LJ switching functions
  2. a test with respect to CHARMM energies for the alchemical endstates
  3. 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
Copy link
Contributor

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.

Copy link
Contributor

@andrrizzi andrrizzi left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

@@ -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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

sigma_string="sigma", bonded_force=True
)
terms = softcore_energy_string.split(";")
softcore_energy_string = ";".join(terms)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.

# No fitting modification found
return {'': [copy.deepcopy(reference_force)]}

if is_custom_lj_force(reference_force) or is_nbfix_force(reference_force):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

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.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

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 "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the () should be {}.

return False


def get_forces_that_define_vdw(system):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe rename this find_vdw_forces for consistency with forces.find_forces?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

assert len(nbfix_forces) == 1


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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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


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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants