Help with hydrogen bond analysis

[Rahul verma]

Dear MDAnalysis Users,

I am beginner for MDAnalysis and am almost asking something naive.

So I have a xyz trajectory which I obtain from in-house  simulation code. My simulation is of simple monoatomic ions (NaCl in liquid water).

I want to calculate the ion-water and water-water hydrogen bond lifetime. Since, I dont have much experience with MDanalysis or with python, i am bit confused about the syntax.

I will discuss them below:

1.After importing MDanalysis and hbonds and defining universe,

import MDAnalysis as mda
from MDAnalysis.analysis import hbonds

u = mda.Universe('nacl.xyz')

I need to define, atomgroup hydrogens  (H = u.select_atoms(' name H') )

Now if you just have an xyz file with atom labels, how do we know what is the resname? Moreover can we define atomgroup hydrogens as mentioned above? 

Similarly (if correct) we may define atomgroup acceptor (A = u.select_atoms(' name O')) for water-water hydrogen bond and (A = u.select_atoms(' name Cl'))   for ion-water hydrogen bond.

2.However, for atomgroup donar, the manual says that this group should be of same length as hydrogen atomgroup. " For many cases, this will mean a donor appears twice in this group."

How do we define donar atomgroup twice as each oxygen is connected to two hydrogens in water molecule?  I mean what is the syntax to increase the length of atomgroup comprising of  oxygens?

3. Then we call hb_ac as in example:

hb_ac = hbonds.HydrogenBondAutoCorrel(u, acceptors = u.atoms.O,
            hydrogens = u.atoms.Hn, donors = u.atoms.N,bond_type='continuous',
            sample_time = 2, nruns = 20, nsamples = 1000)
However, the above command does not have any flag to define the timestep in xyz trajectory. Its possibly available in tpr/xtc kind files but if  you are
using just xyz file, how do we take care of timestep at which trajectory was made?

4. Further I noted there are two utilities for hydrogen bond lifetime calculation.
https://www.mdanalysis.org/docs/documentation_pages/analysis/hbond_autocorrel.html
https://www.mdanalysis.org/docs/documentation_pages/analysis/waterdynamics.html

Is there any difference in between them. Both are based on time-correlation defination of hydrogen bonds?

Regards,
Rahul Verma,
Amity University Lucknow

[orbe...@gmail.com]

Hi Rahul,

I don't have answers for all your questions but at least a few comments.

1) If you need resnames then you can either try to build a "topology" file using information that you have (eg a PDB or PSF file). Or you can add this information to your universe afterwards but we don't have good docs yet on how to do this – it requires adding Residue objects and generally a decent understanding of the "topology system". I would try to build a simple PDB file where you have one residue per water and one per ion – you can just build a PSF with VMD/CHARMM or a TPR with Gromacs and use it, as long as the ordering is the same in your XYZ.

2) I think you should be able to define hydrogens, donors, and acceptors with selections from names, such as what you say:

> I need to define, atomgroup hydrogens  (H = u.select_atoms(' name H') )

(Don't use u.atoms.Hn – this functionality will be removed, use H = u.select_atoms(' name Hn') instead.)

I don't think that resnames are needed but I could be wrong.

Building the matching donor and hydrogen groups for the hbond correlation analysis is a recurring question (and I opened issue  #2181 to improve the docs - feel free to add comments to it, e.g., what you would like to see added to the docs). I am not an expert on this question (but hopefully Richard chimes in), and I think you need to do something like

acceptors = u.select_atoms(' name Cl')

hydrogens = u.select_atoms(' name Hn')

# we need each O twice, in the same order as the Hn (assuming that the 2 H follow or precede their O)

# [O1, O1, O2, O2, ...]

# Create [O1, O2, ...]

oxygens = u.select_atoms("name O")

# use numpy to create an index list of the right shape

# (manually work through the following to understand what happens: we flatten(zip(oxygens, oxygens))

import numpy as np

o_donor_indices = np.ravel(np.transpose([oxygens.indices, oxygens.indices]))

# now create the donor group with slicing

donors = u.atoms[o_donor_indices]

I hope this works...

Then try the analysis tool with these groups.

3) The timestep in your trajectory can be set with the dt argument (eg dt=100 ps)

u = mda.Universe(XYZ, dt=100)

If this doesn't work then file a issue.

4) Regarding the implementations:

> Further I noted there are two utilities for hydrogen bond lifetime calculation.
> https://www.mdanalysis.org/docs/documentation_pages/analysis/hbond_autocorrel.html
> https://www.mdanalysis.org/docs/documentation_pages/analysis/waterdynamics.html
> 
> Is there any difference in between them. Both are based on time-correlation defination of hydrogen bonds?

You need to read the papers in detail. I don't know the answer. The waterdynamics module is primarily for water and calculates a whole bunch of properties in a consistent manner. The hbond_autocorrel module is very flexible in what input it can take and has been used for, e.g., polymers. However, I wish I had a better and more concise answer for you – if the original authors (Richard and Alejandro) read this then it would be great if they could come up with a better answer that we can include in the docs.

Oliver

[Rahul verma]

Dear Oliver,

Thanks for the detailed reply. I found that we don't need to know resnames to perform any hydrogen bond lifetime calculations with MDanalysis.

 By using np.ravel, I could calculate the water-water hydrogen bond lifetime. I did not get the error for donor/hydrogen size.

However, my simulation comprises of 10 sodium and 10 chloride ions with 88 water molecules. For calculating ion-water HB lifetime,

I had to keep the number of donors and hydrogens same. So instead of considering all 176 hydrogens, I chose 160 hydrogens which could match with 10*16 (np.repeat (10 chloride ions)*16) donor atoms length.

So I modified the script a bit as shown below:

import MDAnalysis as mda
from MDAnalysis.analysis import hbonds

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
u = mda.Universe('nacl.xyz',dt=0.5)
H = u.select_atoms('bynum 20:260 and (name H)')
Cl = u.select_atoms('name Cl')
acceptors = u.select_atoms('bynum 19:260 and (name H)')
c_donor_indices = np.ravel(np.repeat(Cl.indices,16,axis=0))
print ([Cl.indices])
print ([H.indices])
donar = u.atoms[c_donor_indices]
print ([donar.indices])
hb_ac = hbonds.HydrogenBondAutoCorrel(u, acceptors = u.atoms.H,
            hydrogens = u.atoms.H, donors = donar, bond_type='continuous',angle_crit=125,dist_crit=5.0,
            sample_time = 5, nruns = 20, nsamples = 10000)
hb_ac.run()

However, for ion-water HB lifetime, I continue to get the same donor/hydrogen number size error.

I even checked the size of donor/hydrogen by printing them out which is also shown below. But I can't remove the error. Any suggestion would be highly appreciated.

[array([10, 11, 12, 13, 14, 15, 16, 17, 18, 19])]  # number of chloride ions
[array([ 21,  22,  24,  25,  27,  28,  30,  31,  33,  34,  36,  37,  39,
        40,  42,  43,  45,  46,  48,  49,  51,  52,  54,  55,  57,  58,
        60,  61,  63,  64,  66,  67,  69,  70,  72,  73,  75,  76,  78,
        79,  81,  82,  84,  85,  87,  88,  90,  91,  93,  94,  96,  97,
        99, 100, 102, 103, 105, 106, 108, 109, 111, 112, 114, 115, 117,
       118, 120, 121, 123, 124, 126, 127, 129, 130, 132, 133, 135, 136,
       138, 139, 141, 142, 144, 145, 147, 148, 150, 151, 153, 154, 156,
       157, 159, 160, 162, 163, 165, 166, 168, 169, 171, 172, 174, 175,
       177, 178, 180, 181, 183, 184, 186, 187, 189, 190, 192, 193, 195,
       196, 198, 199, 201, 202, 204, 205, 207, 208, 210, 211, 213, 214,
       216, 217, 219, 220, 222, 223, 225, 226, 228, 229, 231, 232, 234,
       235, 237, 238, 240, 241, 243, 244, 246, 247, 249, 250, 252, 253,

       255, 256, 258, 259])]  # number of hydrogens considers  upto 259 although the actual number would go upto 283

[array([10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11,
       11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12,
       12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13,
       13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14,
       14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15,
       15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16,
       16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 17,
       17, 17, 17, 17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 18, 18,
       18, 18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 19,

       19, 19, 19, 19, 19, 19, 19])]  # number of chloride ions repeated 16 times to match the size of hydrogen array

/home/rverma/.local/lib/python2.7/site-packages/MDAnalysis/core/topologyattrs.py:507: DeprecationWarning: Instant selector AtomGroup['<name>'] or AtomGroup.<name> is deprecated and will be removed in 1.0. Use AtomGroup.select_atoms('name <name>') instead.
  DeprecationWarning)
Traceback (most recent call last):
  File "shb_ion.py", line 21, in <module>
    sample_time = 5, nruns = 20, nsamples = 1000)
  File "/home/rverma/.local/lib/python2.7/site-packages/MDAnalysis/analysis/hbonds/hbond_autocorrel.py", line 228, in __init__
    raise ValueError("Donors and Hydrogen groups must be matched")
ValueError: Donors and Hydrogen groups must be matched

Thanks,

Rahul

--
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

[Oliver Beckstein]

Hi Rahul,

I am not 100% sure but it looks to me that you are pairing Cl- ions with H atoms, which implies that you are considering the Cl- ions as bonded to the H atoms. I don’t think that this makes sense. If you think that there could be H-bonds between Cl- and water molecules then you need to set up

donors: water oxygens

hydrogens: water hydrogens

acceptors: Cl- (and water oxygens)

The development version of MDAnalysis also has a new function to help with finding donors (and updated docs): https://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/hbond_autocorrel.html#input

Best,

Oliver

--

Oliver Beckstein

orbe...@gmail.com * orbe...@gmx.net