Multiple alchemical regions

[adw62]

Provides feature for multiple alchemical regions, now with tests.

Todo:
Tests do not cover the interaction between alchemical regions.
Using reaction field causes errors in testing.

[andrrizzi]

Looks fantastic! Thank you!

I've only found one possible bug, but mostly minor things. I'm also about to push a few minor fixes. The only non-minor thing I changed was moving the copy of the reference_force above the for loops fixing the 0.0 sigmas so that the original system would not be modified.

I only went through alchemy.py so far and I need to switch, but I'll try to go through the tests tomorrow.

[andrrizzi]

It may be possible to avoid hardcoding these by having check_parameters finding the region name with something like

if parameter_name.startswith(parameter):
    region_name = parameter_name[len(parameter)+1:] if len(parameter_name) > len(parameter) else ''
    return True, region_name

[adw62]

Do you think this can remove all the hard coding in _find_force_components or just the lines you highlighted?

[andrrizzi]

I think it could remove the hard coding in _find_force_components altogether. You'd have to substitute this pattern

if check_energy_expression(force, 'lambda_sterics'):
    if check_energy_expression(force, 'lambda_sterics_zero'):
        sterics_bond_forces['zero'].append([force_index, force])
    ...

with something like

if check_energy_expression(force, 'lambda_sterics'):
    # region_name is None if the parameter is not found in the force and '' if it has no suffix.
    region_name = check_parameter(force, 'lambda_sterics')
    sterics_bond_forces[region_name].append([force_index, force])

unless there's something I don't see.

[andrrizzi]

Ah! Sorry, I remember now that I used check_energy_expression for the alchemical-alchemical custom forces that are decoupled instead of annihilated (and have no lambda_sterics parameter). To make it work we'll have to modify check_energy_expression too, not just check_parameter. I'm thinking using regex.

[andrrizzi]

Something like this works for me (see https://stackoverflow.com/questions/15340582/python-extract-pattern-matches)

>>> a = 'asdfa sdfalambda_sterics_zero*asdf'
>>> import re
>>> p = re.compile(r'lambda_sterics_([A-Za-z0-9]+)')
>>> p.search(a).group(0)
'lambda_sterics_zero'
>>> p.search(a).group(1)
'zero'

I think you can do a simple parameter in custom_force.getEnergyFunction() as it is done now, and then try that regex to see if it has a region associated.

[adw62]

I've reworked _find_force_components to remove the hard coded region names. It's not final and is causing HostGuestExplicit with exact PME to fail currently. But is this rework similar to what you were thinking?

[adw62]

Finished off _find_force_components think it's working as intended. Let me know if there are any more changes that need to be made.

[andrrizzi]

Done! Tests look great too.

I've removed the need for hardcoded 'zero','one','two' region names in is_alchemical_pme_treatment_exact in my last commit.

[adw62]

I think this simplification looks good. If i remember I had if len(lambda_var_suffixes) > 1: as the outer loop so that this did not have to be checked for every exception in the nonbonded force. If you think the performance will not be impacted I prefer how you have written it. Although I have changed it a little to remove one of the if statements. Let me know what you think.

[andrrizzi]

I see what you mean. Yeah, I think I'd prefer to simplify the code in this case as the cost of the if statement will be negligible w.r.t. the rest of the operations inside the foor loop.

Although I have changed it a little to remove one of the if statements

Much better now, thanks!

[andrrizzi]

Just looked at this briefly, but it looks like what I had in mind 👍

[andrrizzi]

Awesome! Thanks so much!

I think we are ready to merge this into master if there are no objections (@hannahbrucemacdonald @ajsilveira @dominicrufa).

[andrrizzi]

Merging! 🎉

[adw62]

Great. Thank you for all your help with this!

[andrrizzi]

No problem!