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.

Sign up for GitHub

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

Jump to bottom

Add pbc kwarg to Contacts #2646

Merged
merged 8 commits into from Apr 22, 2020
Merged

Add pbc kwarg to Contacts #2646

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
921eb56
Fixes #2368: added option to calculate distances using PBC
ss62171 Mar 25, 2020
d796742
Added test to check distance calculated using PBC
ss62171 Mar 25, 2020
856344c
test added
ss62171 Mar 29, 2020
05a442c
documented changes
ss62171 Mar 29, 2020
c97cacc
added box_dim and test modified
ss62171 Mar 31, 2020
8b204d2
added lambda func to get dimension of box
ss62171 Mar 31, 2020
ee62918
edit changelog and docstring of class
ss62171 Apr 8, 2020
043819b
fixed typo
ss62171 Apr 9, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
24 changes: 19 additions & 5 deletions package/MDAnalysis/analysis/contacts.py
View file
Open in desktop
Expand Up @@ -424,18 +424,28 @@ def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
self.select = select
self.grA = u.select_atoms(select[0])
self.grB = u.select_atoms(select[1])

self.is_box = self.fraction_kwargs.get('is_box')
Copy link
Contributor

Choose a reason for hiding this comment

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

The name is_box isn't particularly intuitive; something like pbc would be better.

Also, since it is only used in __init__ and not in the fraction_contacts methods where fraction_kwargs is passed, it'd make more sense to directly include as an argument of __init__ rather than through kwargs. This would also make it easier to include in the docs here.


# contacts formed in reference
self.r0 = []
self.initial_contacts = []

if isinstance(refgroup[0], AtomGroup):
refA, refB = refgroup
self.r0.append(distance_array(refA.positions, refB.positions))
if(self.is_box):
self.r0.append(distance_array(refA.positions, refB.positions,
box=refA.universe.dimensions))
else:
self.r0.append(distance_array(refA.positions, refB.positions))
self.initial_contacts.append(contact_matrix(self.r0[-1], radius))
else:
for refA, refB in refgroup:
self.r0.append(distance_array(refA.positions, refB.positions))
if(self.is_box):
self.r0.append(distance_array(refA.positions, refB.positions,
box=refA.universe.dimensions))
else:
self.r0.append(distance_array(refA.positions, refB.positions))

self.initial_contacts.append(contact_matrix(self.r0[-1],
radius))

Expand All @@ -444,9 +454,13 @@ def _prepare(self):

def _single_frame(self):
self.timeseries[self._frame_index][0] = self._ts.frame

# compute distance array for a frame
d = distance_array(self.grA.positions, self.grB.positions)
if(self.is_box):
d = distance_array(self.grA.positions, self.grB.positions,
self._ts.dimensions)
else:
d = distance_array(self.grA.positions, self.grB.positions)

for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts,
self.r0), 1):
Expand Down
25 changes: 25 additions & 0 deletions testsuite/MDAnalysisTests/analysis/test_contacts.py
View file
Open in desktop
Expand Up @@ -37,6 +37,8 @@
from MDAnalysisTests.datafiles import (
PSF,
DCD,
TPR,
XTC,
contacts_villin_folded,
contacts_villin_unfolded,
contacts_file
Expand Down Expand Up @@ -302,6 +304,29 @@ def test_non_callable_method(self, universe):
with pytest.raises(ValueError):
self._run_Contacts(universe, method=2, stop=2)

def test_distance_box(self):
u = mda.Universe(TPR, XTC)
sel1 = "(resname ARG LYS) and (name NH* NZ)"
sel2 = "(resname ASP GLU) and (name OE* OD*)"
acidic = u.select_atoms(sel2)
basic = u.select_atoms(sel1)
Copy link
Contributor

Choose a reason for hiding this comment

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

You could focus the test on a particular pair we know is split across the periodic boundary, to make it easier to directly check the distances measured are what we expect given PBC.


q = contacts.Contacts(u,
select=(sel1, sel2),
refgroup=(acidic, basic),
method="hard_cut",
kwargs={'is_box': True})

r = contacts.Contacts(u,
select=(sel1, sel2),
refgroup=(acidic, basic),
method="hard_cut")

q.r0 = np.array(q.r0).squeeze()
r.r0 = np.array(r.r0).squeeze()
assert_array_equal(q.r0,r.r0)
Copy link
Contributor

Choose a reason for hiding this comment

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

We expect the distance to be different with or without consideration of the periodic boundaries, so this test is currently failing; maybe you could compare directly with what we expect the values to be instead?




def test_q1q2():
u = mda.Universe(PSF, DCD)
Expand Down