I have two questions. I just want to check whether I am working correctly or not. I have done protein-ligand simulation in Gromacs. For generating the hbonds i am creating the index file for donor-hydrogen-acceptor. Now, the atom number of donor in em.gro file is 2121, the atom number of hydrogen in 2122, and the atom number of acceptor is 1158. So to create the index file I am giving the command gmx make_ndx -f em.gro -n index.ndx
Then
a 2121 | a 2122 creates group 25
a 1158 creates group 26
25 | 26
q
Is this the correct way of making index for analyzing h bond between an acceptor and donor? Can someone please confirm.
Secondly, when i created the hbnum.xvg file it shows three columns: first column is for time (ns), second column - Hydrogen bonds (S0) is for hydrogen bonds that satisfy the distance (<0.35 nm) and the angle (<30 degrees) criteria, and third column (S1) - Pairs within 0.35 nm is for donor-acceptor pairs that satisfy the distance criteria ignoring the angle criteria. Now for plotting the h bonds with time which column should I be selecting? I am thinking about going with S0 column (Hydrogen bonds). Is it the correct way? What is the purpose of column S1. Where can I use the data from S1 column and what information can be inferred from that.
Your approach for defining the group is fine, though you can do it either with all three atoms in one group (and selecting it twice when prompted) or to have the donor-H and acceptor groups separate. Either works.
The output is what it says. The first column is H-bonds (satisfy the distance and angle criteria) and the second column is just atom pairs that satisfy the distance criterion. This distinction can be useful for studying the directionality of the interactions. Plot what is useful to you. If you want actual H-bonds, it’s the first data column.
Dr Lemkul thanks a lot for earlier responses. Just a quick follow-up question. If the number of interactions are more in the third column as compared to the second column then will it be appropriate to consider that the interaction between donor-hydrogen-acceptor is not in any specific direction (due to different angles formed between X-H-Y) and hence the resulting h bond might be weaker.
The number of contacts will always be larger than the number of actual hydrogen bonds (or at least equivalent to it). Strength of an interaction is an energetic/thermodynamic concept, not something you can get from these kinds of data, but a hydrogen bond is always stronger the more linear it is.
Hi @jalemkul, I’m a bit confused about this. Do you mean that the S1 legend “Pairs within 0.35 nm” has a higher value than the S0 legend “Hydrogen bonds” in hbnum.xvg?
Yes. The pairs will be greater than, or equal to, the number of H-bonds in a given frame. The pairs column indicates all the interatomic pairs that satisfy the distance criterion, but not necessarily the angle criterion.
selecting the PEG and PSS groups for the calculations. In the resulting hbnum.xvg file, I noticed that at certain picoseconds, the value for “Pairs within 0.35 nm” is lower than that for “Hydrogen bonds.” Is this possible, or does it indicate an error in my calculations?
@jalemkul@Ema2022 I’ve always noticed the same in my results also. The “Pairs within 0.35 nm” column values have been lesser than the “Hydrogen bonds” column. Is this some error?
The pairs column would refer to what we would get using mindist and cut-off of 0.35 nm right?