The format of timerseries for hydrogen bond and water bridge analysis

[xiki-tempula]

When one is interested in understanding the hydrogen bond from a hydrogen bond donor to a hydrogen bond acceptor, one might want to know the name of the heavy atom corresponding to the hydrogen bond donor. The current practice implemented in hydrogen bond analysis is that when a hydrogen bond is detected only the donor and acceptor were recorded and the heavy atom corresponding the donor is added later, which cause problems such as #1687.

@orbeckst has suggested that a workaround will be incorporating the heavy atom information directly into the timeseries object so as to save the trouble afterword. I'm also in favour of this option. I have a PR open which adds multi-order water bridge support to the existing water bridge analysis which will, rather unfortunately, break the existing time series attribute of the water bridge analysis.

So the question is how should we reformulate the time series attribute. Currently, the format has this form.

[<atom1 index>, <atom2 index>, (atom1 resname, atom1 resid, atom1 name), (atom2 resname, atom2 resid, atom2 name), <distance>, <angle>]

In hydrogen bond analysis, atom 1 is the donor and atom 2 is the receptor. In water bridge analysis, since two donors can be joined by a water molecule. The atom1 is the atoms more closely linked to selection 1 and atom2 is the atom more closely linked to selection 2. So if hydrogen bond acceptor α from selection 1 and hydrogen bond donor atom β from selection B is joined by atom γ and δ from a water molecule. The time series will have this form.

[<α index>, <γ index>, (α resname, α resid, α name), (γ resname, γ resid, γ name), <distance>, <angle>],
[<δ index>, <β index>, (δ resname, δ resid, δ name), (β resname, β resid, β name), <distance>, <angle>],

My suggestion of how we incorporate heavy atom information will be adding the heavy atom name to the last column of the tuple. So that the time series will have the form:

[<α index>, <γ index>, (α resname, α resid, α name), (γ resname, γ resid, γ name, **γ heavy atom name**), <distance>, <angle>],
[<δ index>, <β index>, (δ resname, δ resid, δ name), (β resname, β resid, β name, **β heavy atom name**), <distance>, <angle>],

For hydrogen-bond analysis, instead of taking the form of where we have donor precede acceptor.

[<donor index>, <acceptor index>, (donor resname, donor resid, donor name), (acceptor resname, acceptor resid, acceptor name), <distance>, <angle>]

We have atoms from selection 1 precede the selection 2.

[<sele1 index>, <sele2 index>, (sele1 resname, sele1 resid, sele1 name), (sele2 resname, sele2 resid, sele2 name), <distance>, <angle>]

The heavy atom name is added to the atom which is hydrogen bond donor. If sele1 is hydrogen bond donor, we have:

[<sele1 index>, <sele2 index>, (sele1 resname, sele1 resid, sele1 name, **sele1 heavy atom name**), (sele2 resname, sele2 resid, sele2 name), <distance>, <angle>]

If sele2 is hydrogen bond donor, we have:

[<sele1 index>, <sele2 index>, (sele1 resname, sele1 resid, sele1 name), (sele2 resname, sele2 resid, sele2 name, **sele2 heavy atom name**), <distance>, <angle>]

Thus, we have a consistent format between hydrogen bond analysis and water bridge analysis and the two functionality can be merged together. The donor and acceptor can be relatively easy distinguished by the length of the third and fourth column as the donor atom has a length of four (atom resname, atom resid, atom name, atom heavy atom name) and the acceptor has a length of three (atom resname, atom resid, atom name).

I'm welcoming suggestions.

[orbeckst]

I am in favor of making the output format more uniform. A couple of questions/comments

• Do we really need to store information that could be deduced via the atom index, i.e., would it be sufficient to only store indices? For example, when I have <α index>, I can get (α resname, α resid, α name) from the topology.

  We could then just store

  <α index>, <γ index>, <distance>, <angle>]

  (in the old scheme without the heavy atom information. This has the advantages that it saves memory and is DRY.

  If we need data that can be analyzed without the topology then we should write an export() method that builds a dataframe with all the data and dumps it to a CSV or other file (using pandas).

• I really dislike the list-of-list data structure. It would be much better if we store in something that is closer to tidy data, i.e., "one line – one observation". It then becomes trivial to do pandas-style analysis. In this case, each line should just be tagged with the trajectory frame number (or perhaps the time, but frame numbers are less error prone). One can then groupby to get per-time-step information. So something like the "table format" that we also generate:

  (<frame>, <α index>, <γ index>, <distance>, <angle>)

  In principle, our data structure should already look like a tidy DataFrame.

• For the heavy atom information: As @xiki-tempula pointed out, this is difficult to ascertain after the fact (not the least because some atoms can act as either donors or acceptors). It seems much more convenient to add this information to the data structure. Again thinking in terms of a tidy DataFrame in which either α or γ can be the donor, this means that we need two additional columns for connected heavy atoms. For a discussion I offer

  (<frame>, <α index>, <α heavy index>, <γ index>, <γ heavy index>, <distance>, <angle>)

  where <α heavy index> and <γ heavy index> contain either

  • index of the connected heavy atom, if <index> is a hydrogen, or (i.e., a donor)
  • None, if <index> is already a heavy atom (i.e., a acceptor).

  We still retain a tidy DataFrame (even though it is not completely DRY anymore) so in this way one can easily (i.e., with operations that work on columns of a dataframe or array) detect donors and acceptors (test if the heavy index value is None) . This approach is easier than the proposed one where one would have to test the length of a tuple (length 3: acceptor, length 4: donor), which requires unpacking of the nested data structure.

I think that the above format can be used to unify hydrogenbond and hydrogenbridge analysis, please comment.

Final thought: if we are changing the format then we might as well do it in a radical fashion that includes all the lessons that we have learned over time.

[kain88-de]

Before you make any changes to the water bridge output please consider it is a very widely used feature of MDAnalysis. It is the most frequently cited analysis in papers! I don’t have exact numbers but I looked at almost every paper this year citing MDAnalysis.

Maybe we keep the old time series and add the improved version with a new name.

I really appreciate the work on this module! Such Breaking changes unfortunately keep users from upgrading as that means they have to rewrite their own code (everyone likes to avoid this because the next paper deadline is close by)

[xiki-tempula]

I agree that the [<α index>, <γ index>, <distance>, <angle>] format is a very good idea since it saves memory and adheres to a good software principle. However, considering that timeseries is a public attribute and people who are not very familiar with python and MDAnalysis might also want to learn something from it, a more informative representation might be more helpful. I, personally, had people emailing me asking the details of using water bridge analysis, so a compact display might cause unnecessary confusion for the user. Furthermore, to be honest, this type of analysis is not very efficient and it takes months to generate even several hundred megabytes of data, I believe the memory usage probably won't be such an issue.

I like the idea of incorporating the frame number and treat each entry as a single line, it might come in handy if we eventually will implement the multithreading support for the analysis module.

[orbeckst]

Thanks for this information.

Perhaps we change the format internally and then do what we are already doing: convert the internal format to the legacy format if necessary. It means that the classes accrue some crud but that's the (or a) difference between a widely used product and a new project where one has to pay no attention to others.

An alternative approach is to leave the old code in place, flag as deprecated, and rewrite the module from scratch.

[xiki-tempula]

I think changing the format internally will be a good practice.

Personally, I favour deprecating the hydrogen bond analysis. As the new multi-order water bridge analysis implicitly calculates the direct hydrogen bond and will give the same output when the order is set to zero.

The only two difference is hydrogen bond analysis supports the heuristic method of finding hydrogens which shouldn't be used anyway. The more concerning difference is hydrogen bond analysis puts donor before acceptor while water bridge analysis puts selection 1 before selection 2. However, if we were to add heavy atom support then the two formats can be switched with a single option.