add outputs to waterdynamics.AngularDistribution ?

[yannclaveau]

Hello,
The angularDistribution (AD) class gives too "high level" information according to me. Here is what I mean:

I would like to plot a 2D histogram representing $\theta(\mu, \vec{R})$ where $\mu$ is the dipolar moment of water molecules and $\vec{R}$ the radial position with respect to z-axis.
My point is that, I would to see how water molecules are oriented when they are close to the interface (water molecules in a meso-porous material).

I do not see how I can achieve that with angularDistribution (AD). Even if it wasn't exactly what I would like to do, I tried to select atoms near the interface, however the cylayer keywords doesn't select the whole molecules (using it, I do not have 2 H for 1 O). Maybe is there a way to say "if atoms with molecule ID ("resid XXX" in my case) is in cylayer, then select all the atoms with the same molecule ID" ? But still, it would'nt be exactly what I am looking for (but a possible workaround).

I think, it would be more general if AD output could be similar to hbonds class:
<O index, O x-position, O y-position, O z-position, $\theta$&gt;

Then it would be our responsability to plot an histogram, as it is with hbonds, and why not propose the actual AD class as a "higher level" .

For instance, with hbonds, I did this:

for frame, acceptor_ix, *_ in hbonds.results.hbonds:
    u.trajectory[frame.astype(int)]
    acceptor = u.atoms[acceptor_ix.astype(int)]
    R     = np.linalg.norm(acceptor.position[0:2])
    hist, *_ = np.histogram(R, bins=bin_edges)
    counts  += hist * 2  # multiply by two as each hydrogen bond involves two water molecules

vol = 2*np.pi*step*(zhi-zlo)*bin_centers
rhoHB = counts / vol / hbonds.n_frames

I would like to do something similar but with the dipolar-moment orientation.

Thanks,
Yann

[yannclaveau]

For now, I have create a function based on _getCosTheta. It does the job. I can plot a 2D histogram depending on radial position and $\cos(\theta)$

import numpy as np
from MDAnalysis.lib.log import ProgressBar

#### from _getCosTheta of MDAnalysis.analysis.waterdynamics.AngularDistribution
def getCosTheta(selection, axis):

    line = selection.positions

    Ot0 = line[::3]
    H1t0 = line[1::3]
    H2t0 = line[2::3]

    RH2O = np.linalg.norm(Ot0[:,0:2], axis=1)

    OHVector0 = H1t0 - Ot0
    HHVector0 = H1t0 - H2t0
    dipVector0 = (H1t0 + H2t0) * 0.5 - Ot0

    unitOHVector0  = OHVector0  / np.linalg.norm(OHVector0, axis=1)[:, None]
    unitHHVector0  = HHVector0  / np.linalg.norm(HHVector0, axis=1)[:, None]
    unitdipVector0 = dipVector0 / np.linalg.norm(dipVector0, axis=1)[:, None]

    if axis == "z":
        valOH  = (unitOHVector0[:,2])
        valHH  = (unitHHVector0[:,2])
        valdip = (unitdipVector0[:,2])

    elif axis == "x":
        valOH  = (unitOHVector0[:,0])
        valHH  = (unitHHVector0[:,0])
        valdip = (unitdipVector0[:,0])

    elif axis == "y":
        valOH.append(unitOHVector0[:,1])
        valHH.append(unitHHVector0[:,1])
        valdip.append(unitdipVector0[:,1])

    return (valOH, valHH, valdip, RH2O)

  
    
cosOH, cosHH, cosdip, RH2O = [], [] , [], []     

for ts in ProgressBar(u.trajectory, verbose=True, total=u.trajectory.n_frames):
    valOH, valHH, valdip, rH2O = getCosTheta(H2O.ts, 'z')
    
    cosOH  = np.append(cosOH, valOH)
    cosHH  = np.append(cosHH, valHH)
    cosdip = np.append(cosdip, valdip)
    RH2O   = np.append(RH2O, rH2O)
 

[hmacdope]

Hi Yann, thanks for your suggestions! I'm glad you found a solution that works for you.

We always welcome suggestions for API changes, however as I am not the author or a heavy user of waterdynamics I am poorly placed to comment on them. Perhaps other @MDAnalysis/coredevs will know more.

[orbeckst]

You're very welcome to propose changes. We definitely like code being cleaned up and modularized. We (or rather: our users) don't generally like breaking existing code so we tend to minimize any API breaks. But if you find a way to maintain the existing code from the user perspective (even though it might change under the hood) while at the same time adding new functionality then that would be interesting.

Typically we can really only give more detailed feedback when we see code so the best way forward is to put up a pull request. At this point we see what your idea is and we can point out issues in the broader context of MDAnalysis, if we see any. We'll then also see what new tests would be needed. There's no guarantee that we add code although from past experience, most PRs will eventually get merged if the author remains committed to seeing them through until the end. But it can be a long-ish process because we want to maintain high code quality and minimize technical debt.

[yannclaveau]

I think, in order to keep the existing code and add new functionnality would be to modify the _getCosTheta so that it can be accessed "outside" the function.
using waterdynamics.AngularDistribution.getCosTheta
I tried to use it that way, but because of how the selection works in this class I didn't succed. But, as I am a noob in pyhton, I guess it is probably because I did something wrong:

H2O = u.select_atoms('type 1 or type 2')

bins = 100

AD_analysis = AD(universe=u, select='type 1 or type 2', bins=bins,axis='z')

valOH, valHH, valdip = AD_analysis._getCosTheta(u,H2O,'z')

return Atom has no attribute positions. H2O.positions returns the atom positions of the first frame, but the method uses H2O[i].positions. I understand that it comes from the _selection_serial method, but I do not understand how (and if it is possible) to access the internal method _getCosTheta().

I will create a PR as you suggested.