Regarding the atom/residue naming and selection in HBond autocorrelation function code in water dynamics analysis

[cy16f01]

Hello, The following is the .gro file of glycine-water mixture

Glycine-Water
XXX
1GLYN N 1 1.554 0.192 1.237 0.2337 0.7678 -0.4403
1GLYN HT1 2 1.595 0.113 1.190 1.2957 1.3150 -0.4663
1GLYN HT2 3 1.567 0.275 1.179 0.2555 1.4287 0.4977
1GLYN CA 4 1.596 0.220 1.375 -0.1659 0.2593 -0.2142
1GLYN HA1 5 1.690 0.279 1.386 0.3912 -0.2223 -2.2182
1GLYN HA2 6 1.527 0.296 1.419 -0.2036 -0.1122 0.3771
1GLYN C 7 1.590 0.109 1.479 -0.4882 0.5503 0.0814
1GLYN OT1 8 1.489 0.061 1.528 -0.1025 -0.2671 0.0940
1GLYN OT2 9 1.711 0.048 1.513 -0.4918 0.6753 0.3137
1GLYN HO2 10 1.698 -0.035 1.560 -0.8295 -0.2734 -1.4207
2SOL OW 11 1.974 1.991 1.614 0.2492 -0.3487 0.3687
2SOL HW1 12 2.040 1.918 1.594 -0.6780 -1.3661 0.9449
2SOL HW2 13 1.884 1.951 1.630 -0.5310 1.0328 -0.4632
3SOL OW 14 0.331 0.052 1.858 -0.2124 0.4141 -0.8687
3SOL HW1 15 0.310 -0.007 1.936 0.1912 -0.3242 -1.3072
3SOL HW2 16 0.293 0.011 1.775 0.6684 0.4113 -1.2785
4SOL OW 17 0.482 1.719 1.590 -0.5622 -0.3605 0.2368
4SOL HW1 18 0.478 1.632 1.541 0.8012 -1.4225 1.9691
4SOL HW2 19 0.444 1.707 1.682 0.8493 0.5546 0.9465
5SOL OW 20 0.168 1.334 1.467 0.1289 -0.3488 0.2064
5SOL HW1 21 0.149 1.375 1.556 3.3113 1.1012 0.2810
5SOL HW2 22 0.167 1.234 1.475 1.0811 -0.2295 2.0747
.
.
.

And i am trying to analysie the HBL(continuous and intermittent) for only water molecules, from the water dynamics code..as follows:-

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import HydrogenBondLifetimes as HBL
import matplotlib.pyplot as plt
u = MDAnalysis.Universe('md1.gro','md1.xtc')
selection1 = "resname SOL"
selection2 = "resname SOL"
HBL_analysis = HBL(u, selection1, selection2, 0, 2000, 30)
HBL_analysis.run()
time = 0
#now we print the data ready to plot. The first two columns are the HBLc vs t
#plot and the second two columns are the HBLi vs t graph
for HBLc, HBLi in HBL_analysis.timeseries:
    print("{time} {HBLc} {time} {HBLi}".format(time=time, HBLc=HBLc, HBLi=HBLi))
    time += 1
#we can also plot our data
plt.figure(1,figsize=(18, 6))
#HBL continuos
plt.subplot(121)
plt.xlabel('time')
plt.ylabel('HBLc')
plt.title('HBL Continuos')
plt.plot(range(0,time),[column[0] for column in HBL_analysis.timeseries])
#HBL intermitent
plt.subplot(122)
plt.xlabel('time')
plt.ylabel('HBLi')
plt.title('HBL Intermitent')
plt.plot(range(0,time),[column[1] for column in HBL_analysis.timeseries])
plt.show()

I am getting the values,

But suppose if i give different naming as in selection1 = "name OW " , selection2 = "name HW1 and name HW2" since this will calculate between intermolecular HBL..i am getting the values as:-

0 1.0 0 0.0
1 0.0 1 0.0
2 0.0 2 0.0
3 0.0 3 0.0
4 0.0 4 0.0
5 0.0 5 0.0
6 0.0 6 0.0
7 0.0 7 0.0
.
.

which is not correct and has some error, due to wrong selection/naming of atoms/residues......which i am not able to find it out... So i needed help regarding this...

Thank you...

[orbeckst]

I am not sure how the waterdynamics.HydrogenBondLifetimes class identifies donors and acceptors. Have a look at the code itself.

Perhaps @alejob or @cyanezstange can comment.

[alejob]

We are using the hbond_analysis module of MDAnalysis to identify donors and acceptors, and in this way to count the number of hydrogen bonds in a determined timestep.

edit: You should remember that an h-bond is defined as an electrostatic force and not confounding with a covalent bond.

I don't know if hbond_analysis is able to calculate h-bond inside a molecule, does it @orbeckst ?

[alejob]

By the way, we are not using a custom definition of donors and acceptors as in HydrogenBondAutoCorrel, we are using the default values given by hbond_analysis, you can see it here. In other words, you can't set yours own aceptors and donors.

[orbeckst]

I don't think there is a problem calculating h bonds inside a molecule, only the definitions of donors and acceptors matter.

[orbeckst]

@cy16f01 did these comments answer your question?

[orbeckst]

No more comments. Closing.