Add pbc kwarg to Contacts #2646

ss62171 · 2020-03-25T20:36:40Z

Changes made in this Pull Request:

Added a kwargs boolean option 'is_box' whether to take PBC into account or not.

PR Checklist

Tests?
Docs?
CHANGELOG updated?
Issue raised/referenced?

fiona-naughton · 2020-03-29T12:39:51Z

package/MDAnalysis/analysis/contacts.py

@@ -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')


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.

fiona-naughton · 2020-03-29T13:04:35Z

testsuite/MDAnalysisTests/analysis/test_contacts.py

+        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)


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.

fiona-naughton · 2020-03-29T13:07:22Z

testsuite/MDAnalysisTests/analysis/test_contacts.py

+
+        q.r0 = np.array(q.r0).squeeze()
+        r.r0 = np.array(r.r0).squeeze()
+        assert_array_equal(q.r0,r.r0)


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?

ss62171 · 2020-03-29T17:51:47Z

@fiona-naughton @orbeckst can you please review these changes?

lilyminium · 2020-03-29T21:51:11Z

Just FYI, distance_array will take a box=None argument.

ss62171 · 2020-03-29T22:09:17Z

Just FYI, distance_array will take a box=None argument.

@lilyminium I didn't get this. Can you please clarify :)

codecov · 2020-03-29T22:09:30Z

Codecov Report

Merging #2646 into develop will decrease coverage by 0.04%.
The diff coverage is 100.00%.

@@             Coverage Diff             @@
##           develop    #2646      +/-   ##
===========================================
- Coverage    91.01%   90.96%   -0.05%     
===========================================
  Files          174      174              
  Lines        23595    23314     -281     
  Branches      3082     3083       +1     
===========================================
- Hits         21474    21207     -267     
+ Misses        1498     1487      -11     
+ Partials       623      620       -3

Impacted Files	Coverage Δ
package/MDAnalysis/analysis/contacts.py	`100.00% <100.00%> (ø)`
package/MDAnalysis/auxiliary/base.py	`88.19% <0.00%> (-0.74%)`	⬇️
package/MDAnalysis/coordinates/chain.py	`87.37% <0.00%> (-0.59%)`	⬇️
package/MDAnalysis/coordinates/TRZ.py	`83.46% <0.00%> (-0.38%)`	⬇️
package/MDAnalysis/coordinates/base.py	`93.46% <0.00%> (-0.36%)`	⬇️
package/MDAnalysis/coordinates/DMS.py	`88.23% <0.00%> (-0.34%)`	⬇️
package/MDAnalysis/coordinates/GMS.py	`85.41% <0.00%> (-0.30%)`	⬇️
package/MDAnalysis/coordinates/INPCRD.py	`90.62% <0.00%> (-0.29%)`	⬇️
package/MDAnalysis/coordinates/GSD.py	`86.95% <0.00%> (-0.28%)`	⬇️
package/MDAnalysis/coordinates/TXYZ.py	`90.10% <0.00%> (-0.22%)`	⬇️
... and 37 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 5f51383...043819b. Read the comment docs.

fiona-naughton

A couple comments.

Following @lilyminium's comment, you could set a value box based on pbc (ie. None or dimensions, as appropriate) and run distance_array with that, rather than needing an if/else every time.

fiona-naughton · 2020-03-31T03:41:17Z

package/MDAnalysis/analysis/contacts.py

@@ -392,6 +392,9 @@ def __init__(self, u, select, refgroup, method="hard_cut", radius=4.5,
        method : string | callable (optional)
            Can either be one of ``['hard_cut' , 'soft_cut', 'radius_cut']`` or a callable
            with call signature ``func(r, r0, **kwargs)`` (the "Contacts API").
+        pbc : bool (optional)
+            Uses periodic boundary conditions to calculate distances if set to ``True``; the
+             default is ``False``.


I would set the default to True instead, since this is probably what most users would expect the behaviour to be

fiona-naughton · 2020-03-31T03:46:41Z

package/MDAnalysis/analysis/contacts.py

@@ -461,7 +479,8 @@ def _new_selections(u_orig, selections, frame):
    """create stand alone AGs from selections at frame"""
    u = MDAnalysis.Universe(u_orig.filename, u_orig.trajectory.filename)
    u.trajectory[frame]
-    return [u.select_atoms(s) for s in selections]
+    temp_u = u.copy()


The idea of using copy() here would be to avoid having to reload a universe from the filenames - you can copy u_orig directly.

can I fix this in another pull request?

fiona-naughton · 2020-03-31T03:55:58Z

testsuite/MDAnalysisTests/analysis/test_contacts.py

+        average_contacts_r = np.mean(r.timeseries[:, 1])
+        average_contacts_q = np.mean(q.timeseries[:, 1])
+
+        assert not average_contacts_q == average_contacts_r


A test like this isn't ideal as inequality could be due to any number of things going wrong. I'd find a specific value(s) you expect to be returned and check against that - say by measuring the expected distances with and without PBC in another program or directly using distance_array

lilyminium · 2020-03-31T19:23:27Z

you could set a value box based on pbc (ie. None or dimensions, as appropriate) and run distance_array with that, rather than needing an if/else every time.

Yep, thanks for clarifying @fiona-naughton. However the box dimensions will change for every frame in an NPT ensemble so it might be better to set a box getting function to use in _single_frame instead e.g.

if pbc:
    self._get_box = lambda ts: ts.dimensions
else:
    self._get_box = lambda ts: None

and call it with box=self._get_box(self._ts). This is less readable than if/else-ing though.

ss62171 · 2020-03-31T19:48:51Z

@lilyminium thanks for clarifying and providing better way to implement this :)

ss62171 · 2020-04-02T19:57:03Z

@fiona-naughton any other changes needed?

ss62171 · 2020-04-06T22:16:46Z

@lilyminium can you please review this.

lilyminium

Please also update CHANGELOG and add a .. versionchanged to the docstring of Contacts.

lilyminium · 2020-04-08T07:04:21Z

testsuite/MDAnalysisTests/analysis/test_contacts.py

+                        refgroup=(acidic, basic), radius=6.0, pbc=pbc)
+        r.run()
+        if pbc:
+            expected = [1., 0.43138152, 0.3989021, 0.43824337, 0.41948765,


Looks good @ss62171 -- personally I would change the test so that the expected values are also part of parametrize, e.g.

@pytest.mark.parametrize("pbc,expected", [ (True, [1., 0.43138152, 0.3989021, 0.43824337, 0.41948765, 0.42223239, 0.41354071, 0.43641354, 0.41216834, 0.38334858]), (False, [1., 0.42327791, 0.39192399, 0.40950119, 0.40902613, 0.42470309, 0.41140143, 0.42897862, 0.41472684, 0.38574822]) ])

A new push would also restart the Appveyor tests, which as far as I can tell have failed due to running out of time. I can't seem to restart these jobs individually without a commit -- @MDAnalysis/coredevs have any insight?

lilyminium · 2020-04-09T08:18:56Z

package/MDAnalysis/analysis/contacts.py

@@ -373,6 +373,8 @@ class Contacts(AnalysisBase):
    .. versionchanged:: 1.0.0
       ``save()`` method has been removed. Use ``np.savetxt()`` on
       :attr:`Contacts.timeseries` instead.
+    .. versionchanged:: 2.0.0


Thanks @ss62171, just one final thing -- this should be for version 1.0.0 which hasn't been released yet.

lilyminium · 2020-04-10T11:25:08Z

Thanks @ss62171, looks good to me! Are you happy with this @fiona-naughton?

ss62171 · 2020-04-18T15:37:33Z

@lilyminium any updates on this?

lilyminium · 2020-04-20T06:23:46Z

Sorry for the delay @ss62171, I've been busy. I've requested another review from @fiona-naughton and if she doesn't do so within 24 hours I'll presume that she is too and merge the PR.

Fiona's comments re: pbc=True kwarg name and default, and the test equality have been addressed; the copy() comment refers to now-removed code.

lilyminium · 2020-04-22T06:07:13Z

Thanks for your patience @ss62171, this will make Contacts analysis much nicer to use.

ss62171 added 2 commits March 26, 2020 01:46

Fixes MDAnalysis#2368: added option to calculate distances using PBC

921eb56

Added test to check distance calculated using PBC

d796742

fiona-naughton requested changes Mar 29, 2020

View reviewed changes

test added

856344c

documented changes

05a442c

fiona-naughton previously requested changes Mar 31, 2020

View reviewed changes

added box_dim and test modified

c97cacc

added lambda func to get dimension of box

8b204d2

lilyminium requested changes Apr 8, 2020

View reviewed changes

edit changelog and docstring of class

ee62918

lilyminium requested changes Apr 9, 2020

View reviewed changes

fixed typo

043819b

lilyminium approved these changes Apr 10, 2020

View reviewed changes

lilyminium requested a review from fiona-naughton April 20, 2020 06:22

lilyminium changed the title ~~Contact~~ Add pbc kwarg to Contacts Apr 22, 2020

lilyminium merged commit 400f93c into MDAnalysis:develop Apr 22, 2020

fiona-naughton added the defect label Sep 26, 2023

fiona-naughton added the Component-Analysis label Sep 26, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add pbc kwarg to Contacts #2646

Add pbc kwarg to Contacts #2646

ss62171 commented Mar 25, 2020 •

edited

fiona-naughton Mar 29, 2020

fiona-naughton Mar 29, 2020

fiona-naughton Mar 29, 2020

ss62171 commented Mar 29, 2020

lilyminium commented Mar 29, 2020

ss62171 commented Mar 29, 2020

codecov bot commented Mar 29, 2020 •

edited

fiona-naughton left a comment

fiona-naughton Mar 31, 2020

fiona-naughton Mar 31, 2020

ss62171 Mar 31, 2020

fiona-naughton Mar 31, 2020

lilyminium commented Mar 31, 2020

ss62171 commented Mar 31, 2020

ss62171 commented Apr 2, 2020

ss62171 commented Apr 6, 2020

lilyminium left a comment

lilyminium Apr 8, 2020

lilyminium Apr 9, 2020

lilyminium commented Apr 10, 2020

ss62171 commented Apr 18, 2020

lilyminium commented Apr 20, 2020

lilyminium commented Apr 22, 2020

Add pbc kwarg to Contacts #2646

Add pbc kwarg to Contacts #2646

Conversation

ss62171 commented Mar 25, 2020 • edited

PR Checklist

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ss62171 commented Mar 29, 2020

lilyminium commented Mar 29, 2020

ss62171 commented Mar 29, 2020

codecov bot commented Mar 29, 2020 • edited

Codecov Report

fiona-naughton 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

lilyminium commented Mar 31, 2020

ss62171 commented Mar 31, 2020

ss62171 commented Apr 2, 2020

ss62171 commented Apr 6, 2020

lilyminium left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lilyminium commented Apr 10, 2020

ss62171 commented Apr 18, 2020

lilyminium commented Apr 20, 2020

lilyminium commented Apr 22, 2020

ss62171 commented Mar 25, 2020 •

edited

codecov bot commented Mar 29, 2020 •

edited