Skip to content

Commit

Permalink
bug fix: double counting water in Water bridge analysis (#3120)
Browse files Browse the repository at this point in the history
Fixes #3119

* Fixes bug linked with double counting of waters in WaterBridgeAnalysis
  • Loading branch information
xiki-tempula committed May 7, 2021
1 parent da46e84 commit d7b3c24
Show file tree
Hide file tree
Showing 3 changed files with 26 additions and 1 deletion.
1 change: 1 addition & 0 deletions package/CHANGELOG
Expand Up @@ -23,6 +23,7 @@ The rules for this file:
* 2.0.0

Fixes
* WaterBridgeAnalysis double counts the water (Issue #3119, PR #3120)
* NCDFReader now defaults to a dt value of 1.0 ps when it cannot estimate
it from the first two frames of the file (Issue #3166)
* Get chi1_selection to match canonical gamma atoms (Issue #2044)
Expand Down
Expand Up @@ -1349,7 +1349,8 @@ def traverse_water_network(graph, node, end, route, maxdepth, result):
else:
if node in end:
# check if any duplication happens
if len(route) == len(set(route)):
heavy_atom = [line[3] or line[2] for line in route]
if len(heavy_atom) == len(set(heavy_atom)):
add_route(result, route)
else:
for new_node in graph[node]:
Expand Down
23 changes: 23 additions & 0 deletions testsuite/MDAnalysisTests/analysis/test_wbridge.py
Expand Up @@ -723,3 +723,26 @@ def test_timesteps_by_type(self, wb_multiframe):

timesteps = sorted(wb_multiframe.timesteps_by_type())
assert_array_equal(timesteps[3], [1, 12, 'ALA', 1, 'H', 'ALA', 6, 'O', 0, 2])

def test_duplicate_water(self):
'''A case #3119 where
Acceptor···H−O···H-Donor
|
H···O-H
will be recognised as 3rd order water bridge.
'''
grofile = '''Test gro file
7
1LEU O 1 1.876 0.810 1.354
117SOL HW1 2 1.853 0.831 1.162
117SOL OW 3 1.877 0.890 1.081
117SOL HW2 4 1.908 0.828 1.007
135SOL OW 5 1.924 0.713 0.845
1LEU H 6 1.997 0.991 1.194
1LEU N 7 2.041 1.030 1.274
2.22092 2.22092 2.22092'''
u = MDAnalysis.Universe(StringIO(grofile), format='gro')
wb = WaterBridgeAnalysis(u, 'resname LEU and name O',
'resname LEU and name N H', order=4)
wb.run()
assert len(wb.timeseries[0]) == 2

0 comments on commit d7b3c24

Please sign in to comment.