Failing membrane energy test

[jthorton]

In #1561 there is an energy test test_abfe_protocol.py:TestA2AMembraneDryRun:test_complex_dry_run which is failing with the following error when comparing the total energies of the reference and alchemical systems.

reference: 18354753.825662132 kJ/mol
alchemical: 18354744.076934926 kJ/mol
PME correction: 0.0 kJ/mol
delta : -9.748727206140757 kJ/mol

Changing the GLOBAL_ALCHEMY_PLATFORM in openmmtools.tests.test_alchemy.py to Reference fixes this issue, highlighting that this might be a precision issue with the platform running the test. I have confirmed that all bonded forces agree exactly and that the issue stems from the nonbonded force comparison.

The attached script compares the nonbonded energy components on each of the platforms I have available locally [CPU, Reference and OpenCL] and shows that using the OpenCL platform with single precision leads to a large energy disagreement between the reference and alchemical systems, which is removed when switching to mixed precision. OpenCL in single precision is the default on my machine. What's more this issue can be triggered in an existing test which currently passes by setting the solvent padding to 10 nanometers in the test_abfe_protocol.py:TestT4LysozymeDryRun:test_complex_dry_run test, indicating that this issue accumulates with system size and wasnt seen till we included large systems with a membrane.

I think the ideal solution would be to use the Reference platform for all energy comparison tests however this is much slower in local testing and not great for the CI in general as most tests will get slower. The alternative is to enforce mixed precision when using OpenCL and CUDA platforms which will fix the energy tests we may then need to loosen the strictness of some energy comparisons which dissect the nonbonded energy (see test_abfe_protocol.py:TestA2AMembraneDryRun:_test_energies ) however due to small precision issues.

Script output:

================================================================================
NONBONDED ENERGY COMPARISON
================================================================================

CPU:
Ref NonbondedForce: 18264821.896272 kJ/mol
Alchem NonbondedForce: 18264767.075031 kJ/mol
Alchem CustomNonbonded(s): -27.857044 kJ/mol
#0: -18.417321 kJ/mol
#1: -9.439723 kJ/mol
Alchem CustomBond(s): 82.679269 kJ/mol
Alchem NB total: 18264821.897257 kJ/mol
Difference: +0.000985 kJ/mol

OpenCL_mixed:
Ref NonbondedForce: 18265110.984181 kJ/mol
Alchem NonbondedForce: 18265056.182793 kJ/mol
Alchem CustomNonbonded(s): -27.859857 kJ/mol
#0: -18.420266 kJ/mol
#1: -9.439590 kJ/mol
Alchem CustomBond(s): 82.679572 kJ/mol
Alchem NB total: 18265111.002508 kJ/mol
Difference: +0.018328 kJ/mol

OpenCL_single:
Ref NonbondedForce: 18265108.974100 kJ/mol
Alchem NonbondedForce: 18265064.411875 kJ/mol
Alchem CustomNonbonded(s): -27.859867 kJ/mol
#0: -18.420275 kJ/mol
#1: -9.439592 kJ/mol
Alchem CustomBond(s): 82.679550 kJ/mol
Alchem NB total: 18265119.231558 kJ/mol
Difference: +10.257458 kJ/mol

Reference:
Ref NonbondedForce: 18264895.239638 kJ/mol
Alchem NonbondedForce: 18264840.419630 kJ/mol
Alchem CustomNonbonded(s): -27.858380 kJ/mol
#0: -18.418737 kJ/mol
#1: -9.439643 kJ/mol
Alchem CustomBond(s): 82.679269 kJ/mol
Alchem NB total: 18264895.240520 kJ/mol
Difference: +0.000881 kJ/mol

positions.txt

alchem_membrane_system.xml.zip

ref_membrane_system.xml.zip

compare_platforms.py

[mikemhenry]

IMHO, I think setting the precision to mixed seems like the right balance + loosen any tolerances as needed. I think also breaking the forces down into bonded and non-bonded, then checking with higher tolerances on the bonded forces also should help us from making mistakes.