[WIP] Solvent-solute splitting #439

maxentile · 2019-09-26T16:59:14Z

Description

Intention: Efficient implementations of force factories for splitting solvent-solvent interactions from other interactions, for use in multiple-timestepping schemes.

Status

Incomplete. Currently the splitter using exceptions and the splitter using interaction groups are both giving pretty big force errors, and the one using subtraction is giving forces that agree within a very tight tolerance with the original system. However, dynamics with the one using subtraction is unstable at 2fs...

CC'ing @c-matthews and @jchodera

Todos

Implement feature / fix bug
- Draft implementations
- Fix possible double-counting in interaction-group-based implementation
Add tests
- Add test asserting that forces are similar between reference system and split system (currently failing for 2 of 3 draft implementations)
- Make sure test is using CutoffPeriodic nonbonded method
- Test nonbonded-parameter cloning method separately
Update documentation as needed
Update changelog

test_split_nb_using_exceptions is failing with: Exception: Maximum allowable relative force error exceeded (was 0.02369251; allowed 0.00000100). test_split_nb_using_subtraction is failing with: Exception: Maximum allowable relative force error exceeded (was 0.00001687; allowed 0.00000100)

The fn split_nb_using_exceptions now works as intended, to split into protein/solvent forces.

JDC suggested to use interaction groups only here, instead of a mixture of interaction groups and exceptions. this appears cleaner and passes quick test that the cloned system is equivalent to the original system. Co-Authored-By: John Chodera <john.chodera@choderalab.org>

Thanks to JDC for spotting this. This affected `split_nb_using_interaction_groups` and `split_nb_using_exceptions` but not `split_nb_using_subtraction`... Co-Authored-By: John Chodera <john.chodera@choderalab.org>

Current test behavior: 1. test_split_nb_using_interaction_groups() Exception: Maximum allowable relative force error exceeded (was 740.98931192; allowed 0.00000100). E alchemical_force = 651130.50822817, reference_force = 878.74004125, difference = 651136.97852131 2. test_split_nb_using_exceptions(): Exception: Maximum allowable relative force error exceeded (was 746.21620950; allowed 0.00000100). E alchemical_force = 662413.13233660, reference_force = 887.70499600, difference = 662419.85727222 3. test_split_nb_using_subtraction() Exception: Maximum allowable relative force error exceeded (was 0.00001820; allowed 0.00000100). E alchemical_force = 878.52677463, reference_force = 878.52663474, difference = 0.01599347

…forces using a surrogate force Current test behavior: Exception: Maximum allowable relative force error exceeded (was 0.00001101; allowed 0.00000100). E alchemical_force = 881.84056569, reference_force = 881.84068232, difference = 0.00971232

attempted simplification Exception: Maximum allowable relative force error exceeded (was 1554.77679366; allowed 0.00000100). E alchemical_force = 1372003.46106334, reference_force = 882.44334338, difference = 1372002.43200366 Co-Authored-By: John Chodera <john.chodera@choderalab.org>

jchodera

Added a few comments!

jchodera · 2019-09-27T13:46:36Z

openmmtools/forcefactories.py

+    new_force.setCutoffDistance(nonbonded_force.getCutoffDistance())
+    #new_force.setEwaldErrorTolerance(nonbonded_force.getEwaldErrorTolerance())
+    #new_force.setForceGroup(nonbonded_force.getForceGroup())
+    new_force.setNonbondedMethod(nonbonded_force.getNonbondedMethod())


The nobonded methods may not directly translate from NonbondedForce to CustomNonbondedForce.
See: https://github.com/choderalab/openmmtools/blob/master/openmmtools/alchemy.py#L1488-L1491

Thanks! Will follow this example more closely

jchodera · 2019-09-27T13:50:11Z

openmmtools/forcefactories.py

+    for i in range(num_particles):
+        new_force.addParticle(nonbonded_force.getParticleParameters(i))
+
+    # now add all the exceptions? # TODO: check if we want to do this...


I think you definitely want to copy all exceptions -> exclusions, or else you'll end up in NaNdyland.

jchodera · 2019-09-27T13:50:52Z

openmmtools/forcefactories.py

+    # find the default nonbonded force
+    force_index, nb_force = forces.find_forces(new_system, openmm.NonbondedForce, only_one=True)
+    # create copies for each interaction. Only half in solvent/solvent and solute/solute as we double-count.
+    nb_only_solvent_solvent = clone_nonbonded_parameters(nb_force, energy_prefactor='0.5*')


Why do you want to multiple the energy by 0.5?

Mistakenly assumed double-counting -- will revert

jchodera · 2019-09-27T13:51:51Z

openmmtools/forcefactories.py

+
+    nb_only_solvent_solvent.addInteractionGroup(set1=solvent_indices, set2=solvent_indices)
+
+    nb_solute.addInteractionGroup(set1=solute_indices, set2=solute_indices)


You can just use

all_indices = solute_indices.union(solvent_indices) nb_solute.addInteractionGroup(set1=solute_indices, set2=all_indices)

Will uncomment this equivalent way (L249 and L256)

jchodera · 2019-09-27T13:52:34Z

openmmtools/forcefactories.py

+    #nb_solute.addInteractionGroup(set1=solute_indices, set2=all_indices)
+
+    # remove original force, add new forces
+    new_system.removeForce(force_index)


I think you want to leave in the original force so you don't have to copy all the exceptions (nonzero 1-4 interactions) to a CustomBondForce.

jchodera · 2019-09-27T13:52:55Z

openmmtools/forcefactories.py

+    # create copies for each interaction. Only half in solute/solute as we double count.
+    nb_only_solvent_solvent = nb_force
+    nb_only_solvent_solute = clone_nonbonded_parameters(nb_force)
+    nb_only_solute_solute = clone_nonbonded_parameters(nb_force,energy_prefactor='0.5*')


Why would you multiply the energy by 0.5?

Mistakenly assumed double-counting -- will revert

jchodera · 2019-09-27T13:53:11Z

openmmtools/forcefactories.py

+    # Remove solute-solute interactions
+    for i in solute_indices:
+        for j in solute_indices:
+            nb_only_solvent_solvent.addException(i, j, 0, 0, 0, replace=True)


This will probably be slow.

May I suggest itertools.product(solute_indices, repeat=2) instead of the nested for? Nothing to do with performance but I'd say it's prettier ✨

jchodera · 2019-09-27T13:53:27Z

openmmtools/forcefactories.py

+
+# surrogate could deviate in functional form from the NonbondedForce defaults (e.g. using soft-core)
+surrogate_energy_expression = default_energy_expression
+half_surrogate_energy_expression = '0.5*' + surrogate_energy_expression


Why use a factor of 0.5?

Mistakenly assumed double-counting -- will revert

maxentile · 2019-09-27T14:14:07Z

Thanks!

jchodera · 2019-09-27T14:20:02Z

It's probably a good idea to make sure the total potential energy matches the unmodified system pretty closely!

maxentile · 2019-09-27T14:24:55Z

It's probably a good idea to make sure the total potential energy matches the unmodified system pretty closely!

This test is currently failing, but we're aiming to match within the same tolerance as the alchemical factories

openmmtools/openmmtools/tests/test_forcefactories.py

Lines 201 to 210 in be65ec6

    
           def check_force_group_decomposition_valid(testsystem, force_group_splitter): 
        
               """force_group_splitter(system, md_topology) returns a copy of system that 
        
               splits system's forces in some convenient way. 
        
               Check that the copy describes the same system overall, by checking whether it 
        
               produces identical forces as the original on a slightly randomized configuration. 
        
               """ 
        
               new_system = force_group_splitter(testsystem.system, testsystem.mdtraj_topology) 
        
               positions = generate_new_positions(testsystem.system, testsystem.positions, nsteps=50) 
        
               compare_system_forces(testsystem.system, new_system, positions)

jchodera · 2019-09-27T18:59:11Z

It's probably much simpler just to compare the potential energy.

maxentile and others added 18 commits September 18, 2019 11:56

ignore pycharm

fee50cd

add utility for generating solvent-solute splitting string

143e067

add force factory for solvent-solute-splitting using exceptions

3fe3b6b

[WIP] add example script for g-BAOAB using solvent-solute splitting

40bc8ed

rename

00dde8d

another attempt to split nonbonded force using subtraction

1d803c3

Modified split_nb_using_exceptions

a380cf1

The fn split_nb_using_exceptions now works as intended, to split into protein/solvent forces.

Added some results for g-baoab

c0a28e1

use CutoffPeriodic in nb_split tests

7be2611

generate a slightly different conformation for testing

ac330ab

update clone_nonbonded_parameters fxn

8ec3266

remove redundant imports

fb0ce28

modified forces weren't in fact being added to the system!

458ebea

Thanks to JDC for spotting this. This affected `split_nb_using_interaction_groups` and `split_nb_using_exceptions` but not `split_nb_using_subtraction`... Co-Authored-By: John Chodera <john.chodera@choderalab.org>

jchodera reviewed Sep 27, 2019

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[WIP] Solvent-solute splitting #439

[WIP] Solvent-solute splitting #439

maxentile commented Sep 26, 2019

jchodera left a comment

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jchodera Sep 27, 2019

maxentile Sep 27, 2019

jaimergp Sep 28, 2019 •

edited

jchodera Sep 27, 2019

maxentile Sep 27, 2019

maxentile commented Sep 27, 2019

jchodera commented Sep 27, 2019

maxentile commented Sep 27, 2019

jchodera commented Sep 27, 2019


		nb_only_solvent_solvent.addInteractionGroup(set1=solvent_indices, set2=solvent_indices)

		nb_solute.addInteractionGroup(set1=solute_indices, set2=solute_indices)

[WIP] Solvent-solute splitting #439

Are you sure you want to change the base?

[WIP] Solvent-solute splitting #439

Conversation

maxentile commented Sep 26, 2019

Description

Status

Todos

jchodera 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jaimergp Sep 28, 2019 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maxentile commented Sep 27, 2019

jchodera commented Sep 27, 2019

maxentile commented Sep 27, 2019

jchodera commented Sep 27, 2019

jaimergp Sep 28, 2019 •

edited