From 469c0899dea79dbe54665e1fdb25b544e615ceb5 Mon Sep 17 00:00:00 2001 From: Max Linke Date: Mon, 28 Mar 2016 20:02:44 +0200 Subject: [PATCH] Add q1q2 analysis wrapper --- package/MDAnalysis/analysis/contacts.py | 27 ++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/analysis/contacts.py b/package/MDAnalysis/analysis/contacts.py index a1aecb276a0..7135838939f 100644 --- a/package/MDAnalysis/analysis/contacts.py +++ b/package/MDAnalysis/analysis/contacts.py @@ -33,6 +33,8 @@ Classes are available for two somewhat different ways to perform a contact analysis: + + 1. Contacts between two groups of atoms are defined with :class:`ContactAnalysis1`), which allows one to calculate *q(t)* over time. This is especially useful in order to look at native contacts during @@ -359,17 +361,18 @@ def _single_frame(self): # compute distance array for a frame d = distance_array(self.grA.positions, self.grB.positions) - y = np.empty(len(self.r0)) + y = np.empty(len(self.r0) + 1) + y[0] = self._ts.frame for i, (initial_contacts, r0) in enumerate(zip(self.initial_contacts, self.r0)): # select only the contacts that were formed in the reference state r = d[initial_contacts] r0 = r0[initial_contacts] - y[i] = self.fraction_contacts(r, r0, **self.fraction_kwargs) + y[i + 1] = self.fraction_contacts(r, r0, **self.fraction_kwargs) if len(y) == 1: y = y[0] - self.timeseries.append((self._ts.frame, y)) + self.timeseries.append(y) def save(self, outfile): """save contacts timeseries @@ -387,6 +390,24 @@ def save(self, outfile): f.write("{frame:4d} {q1:8.6f}\n".format(**vars())) +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) + + +def q1q2(u, selection=('all', 'all'), radius=4.5, method='cutoff', + start=None, stop=None, step=None, **kwargs): + """Helper function to create q1q2 contact analysis""" + first_frame_refs = _new_selections(u, selection, 0) + last_frame_refs = _new_selections(u, selection, -1) + return Contacts(u, selection, + (first_frame_refs, last_frame_refs), + radius=radius, method=method, + start=start, stop=stop, step=step, + **kwargs) + ################################################################################ ################################################################################ ################################################################################