waterdynamics does NOT respect PBC

[JosipL]

(EDIT): The original issue report (see below) indicates that analysis.waterdynamics required PBC-remapped trajectories, which is apparently not clear from the docs. More generally, it should be able to deal with trajectories under PBC directly.

---

Expected behavior

Hi,
First, I am not sure if this is bug, I just think it is :).
I am using MDAnalysis AngularDistribution library to analyze some water/organic systems. Code produce figure without error but I get a strange peak around costheta=0 for all try types of vectors. I do not expect this peak.
Is due to some 0 dividing? Any ideas?

Code to reproduce the behavior

import MDAnalysis
import matplotlib
from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD
import matplotlib.pyplot as plt

u = MDAnalysis.Universe('../water_oct10mMoctacid10mM_surface.tpr','../trr1.trr')
selection = "byres resname SOL and (prop z < 42)"
bins = 30
AD_analysis = AD(u,selection,bins,axis='z')
AD_analysis.run()
#now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector ,
#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns
#are P(cos(theta)) vs cos(theta) for dipole vector
#for bin in range(bins):
#      print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

#and if we want to graph our results
plt.figure(1,figsize=(18, 6))

#AD OH
plt.subplot(131)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for OH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]])

#AD HH
plt.subplot(132)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for HH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]])

#AD dip
plt.subplot(133)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for dipole')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]])

plt.savefig('angular.png')

Current version of MDAnalysis

• Which version are you using? --->2.0.0
• Which version of Python? ---->Python 3.8.6
• Which operating system? ---->Linu
  
  x

[p-j-smith]

Hi @JosipL I'm not very familiar with the code for AngularDistribution, but I've seen artefacts like this when periodic boundaries are not taken into account. For example, if you are finding the orientation of the dipole moment of water with respect to the z-axis but you're not taking periodic boundaries into account, then all water molecules split across the x or y boundaries will have dipole vectors almost size of the x or y dimension. This will mean the dipole moment of all water molecules broken across periodic boundaries will be oriented at 90 degrees to the z axis.. Is there a pbc argument for AngularDistribution?

[orbeckst]

@alejob could you please look into this question regarding waterdynamics?

[JosipL]

Hi,
thank you for your answer p j smith.
There could be something related to PBC, but as shown in the figure there is also a peak around 0 for HH and OH vectors. Also, I have done the same analysis on a similar system and there is no peak around 0 for any vector.

[p-j-smith]

Sorry my mentioning the dipole moment was just an example. The same applies to the OH and HH vectors - if they're split across a periodic boundary in a given dimension, then their length will be approximately equal to the length of that dimension, and the orientation will be in the plane of that dimension.

For your similar system, had you already made the molecules whole, using e.g. gromacs trjconv? If so then there would be no molecules split across periodic boundaries and you wouldn't see this artefact.

Looking at the code for AngularDistribution I think this is indeed the issue:

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]

So periodic boundaries are not being considered here, and there is no option to do so.

Until there is a fix for this, the easiest solution might be to make your molecules whole, either using gromacs or an on-the-fly transformation in MDAnalsis (see here for details of this and Example 1 for how to make your molecules whole: https://www.mdanalysis.org/2020/03/09/on-the-fly-transformations/ although currently this will probably be quite slow to make all of your water molecules whole at each frame).

[JosipL]

Thanks for your time. You were right. PBC was the issue.
What still confuses me is, another case where I do not see this peak even do PBC splitting is not removed. I do not have time to deal with this now in details. If somebody is interested I can share code and trajectory files.
Best,
Josip

[orbeckst]

Hi Paul, Thanks for the detailed analysis (and solving the problem!). Would you mind raising an issue for the lack of PBC awareness in the AngularDistribution code? Just pointing someone to the problematic code might motivate someone to fix the code. On that note: do we have a library function to make vectors conform to minimum-image (similar to what calc_bonds() does internally)? Oliver

…

On Nov 12, 2021, at 2:41 AM, Paul Smith ***@***.***> wrote: Looking at the code for AngularDistribution <https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/analysis/waterdynamics.py#L652-L661> I think this is indeed the issue: 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] So periodic boundaries are not being considered here, and there is no option to do so.

-- Oliver Beckstein (he/his/him) email: ***@***.*** twitter: @orbeckst MDAnalysis – a NumFOCUS fiscally sponsored project https://www.mdanalysis.org/

[orbeckst]

I renamed the issue and update the report. If someone gets to adding PBC to waterdynamics please reference this issue.

[p-j-smith]

Sorry for not getting back to you - things have been quite hectic. If no one gets around to fixing it by next week then I should have some time then to do it.

Would it be okay to call calc_bonds() directly, or would there be a better way to do it?

[richardjgowers]

@p-j-smith I've implemented the function I suspect you need over in MDAnalysis/mdanalysis#3472

[p-j-smith]

Thanks @richardjgowers ! minimise_vectors is going to be useful for so many things, I can't wait to use it