# Hydrogen bonding autocorrelations

I am currently analyzing hydrogen bonding behaviour and kinetics with molecular dynamics simulations. I've been trying to compute autocorrelation functions for hydrogen bonds in a water system. Some authors use different definitions of the HB autocorrelation function (for example, Gowers' definition has been implemented in the MDAnalysis package), but let's say I'll be using Luzar's definition:

$$c(t) = \frac{\langle h(0) h(t) \rangle}{\langle h \rangle}$$ $$h(t)= \begin{cases} 1 : \text{hydrogen bond at time } t \\ 0 : \text{otherwise} \end{cases}$$

I'm looking for some ideas on how to code this on Python. By slightly modifying the MDAnalysis library, I managed to have an array of arrays list_of_sets containing, for each frame of the simulation, the atom ids that form an hydrogen bond. Then, I do something like:

def calculate_c(list_of_sets, tau_max, window_step):
frames = len(list_of_sets)
tau_timeseries = list(range(0, tau_max))
data = [[] for _ in range(tau_max)]

# calculate correlation function
for t in range(0, frames, window_step):
print(f'Analyzing frame {t}...')
Nt = len(list_of_sets[t])

for tau in tau_timeseries:
if t + tau >= frames:
break

# calculate correlation function for each pair of frames in the window
Ntau = len(set.intersection(*list_of_sets[t:t + tau + 1]))
data[tau].append(Ntau * Nt)

# calculate the average correlation function over all frames for each tau
average_timeseries = [np.mean(x) for x in data]

average_timeseries /= average_timeseries[0]

return tau_timeseries, average_timeseries


In this case, I didn't take into account the $$\langle h \rangle$$ parameter for now.

I also tried:

def c_t_luzar(list_of_sets, tau_max=40, window_step=1):
frames = len(list_of_sets)
tau_timeseries = list(range(tau_max))  # start at 0
timeseries_data = [[] for _ in range(tau_max)]

# calculate correlation function
for t in range(0, frames, window_step):
for pair_set in list_of_sets[t]:
h_t0 = 1 if pair_set else 0  # h(0)

for tau in tau_timeseries:
if t + tau >= frames:
break

# calculate correlation function for each pair of frames in the window
h_t = 1 if pair_set in list_of_sets[t + tau] else 0  # 1 if common pair, else 0
timeseries_data[tau].append(h_t)

# calculate the average correlation function over all frames for each tau
average_timeseries = [np.mean(x) for x in timeseries_data]

return tau_timeseries, average_timeseries, timeseries_data


Do you have any opinions on these codes? I've never coded "averages" so I don't know if my train of thought is correct.

This is the methodology that I used to compute continuous hydrogen bond autocorrelation function. I had 21 nucleobases in my system, and the code is written according to it, but should be extendable for any system.

Step 1: Compute hydrogen-bond between all residues using the following code:

import sys          as sys
import numpy        as np
import pandas       as pd
import MDAnalysis   as mda
from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

try:
assert len(sys.argv) == 5, "run the code as python hbond.analysis.py <psf> <dcd> <resname = {ade, guam cyt, thy}> <simType = {add, drude}>"
simType     =   sys.argv[4]
resName     =   sys.argv[3].upper()
dcdName     =   sys.argv[2]
psfName     =   sys.argv[1]
print(resName)
except AssertionError as message:
print(message)
sys.exit(-1)

universe    =   mda.Universe(psfName,dcdName)
hbonds      =   HBA(universe=universe, d_a_cutoff=3.0, d_h_a_angle_cutoff=140.0 )
A           =   universe.select_atoms('resname {}'.format(resName))
B           =   universe.select_atoms('resname {}'.format(resName))

hbonds.hydrogens_sel = "resname TIP2 and name H1 H2"
hbonds.acceptors_sel = "resname TIP3 and name OH2"

# Run the analysis
hbonds.run()

df              =   pd.DataFrame()
df["Frame"]     =   hbonds.hbonds[:, 0].astype(int)
df["Donor"]     =   hbonds.hbonds[:, 1].astype(int)
df["Hydrogen"]  =   hbonds.hbonds[:, 2].astype(int)
df["Acceptor"]  =   hbonds.hbonds[:, 3].astype(int)
df["Distance"]  =   hbonds.hbonds[:, 4].astype(float)
df["Angle"]     =   hbonds.hbonds[:, 5].astype(float)

with open("hbonds.analysis.{}.{}.txt".format(resName, simType),"w") as log:
for data in zip(df["Frame"], df["Donor"], df["Hydrogen"], df["Acceptor"], df["Distance"], df["Angle"]):
log.write("{:10.0f}{:10.0f}{:10.0f}{:10.0f}{:10.3f}{:10.3f}\n".format(data[0], data[1]+1, data[2]+1, data[3]+1, data[4], data[5]))


Step 2: Now compute the autocorrelation function using the following code:

import numpy    as  np
import sys      as  sys
import math     as  mt

nFrames     =   int(sys.argv[1]) # Total number of frames
frameStride =   int(sys.argv[2]) # How many frames do you want to consider in one go
nMolecules  =   int(sys.argv[3]) # Total number of molecules

class Frame:
def __init__(self, frames) -> None:
self.frameData      =   [values for data in frames for values in data]
self.frameIndex     =   set(self.frameData[0::6])
self.donorIndex     =   self.frameData[1::6]
self.acceptorIndex  =   self.frameData[3::6]
self.hbondMatrix    =   np.zeros((nMolecules, nMolecules))

for donor, acceptor in zip(self.donorIndex, self.acceptorIndex):
self.hbondMatrix[int(donor)][int(acceptor)] = 1
self.hbondMatrix[int(acceptor)][int(donor)] = 1

with open("hbonds_analysis_water_298K.txt") as hbondsObject:
firstLine   =   fullData[0].split()
currFrame   =   int(firstLine[0])
prevFrame   =   currFrame
currData    =   np.zeros(6)
frameData   =   []
trajData    =   []
for index, data in enumerate(firstLine):
currData[index] =   data
frameData.append(currData)

# Now read the rest of the lines
for line in fullData[1:]:
currFrame   =   int(line.split()[0])
if prevFrame == currFrame:
currData    =   np.zeros(6)
for index, data in enumerate(line.split()):
currData[index] =   data
frameData.append(currData)
else:
trajData.append(Frame(frameData))
currData    =   np.zeros(6)
frameData   =   []
for index, data in enumerate(line.split()):
currData[index] =   data
frameData.append(currData)
prevFrame = currFrame
trajData.append(Frame(frameData))

analysisSections    =   [trajData[i:i+frameStride] for i in range(nFrames - frameStride + 1)]
hbondFunction       =   np.zeros(frameStride)

for analysisSection in analysisSections:
hbondMatrix =   analysisSection[0].hbondMatrix
for index, frame in enumerate(analysisSection):
hbondMatrix             =   np.multiply(hbondMatrix, frame.hbondMatrix)
hbondFunction[index]    =   sum(sum(hbondMatrix)) + hbondFunction[index]

hbondFunction   =   hbondFunction / (hbondFunction[0])
print(hbondFunction)


You can use scipy.optimize.curvefit to fit the curve to a biexponential function to obtain the slow and fast autocorrelation times.

• I attempted to use this code to calculate an autocorrelation function between two residues, but unfortunately, it didn't work. The file I'm working with has the frame number in the first column, donor index in the second column, hydrogen bond value in the third column, and receptor index in the fourth column. Apr 9 at 16:49
• I gave this as an example code. It was no way meant to work out of the box for your system, because I do not know how your file is organized. You should be able to adapt the code for your system, provided you know python Apr 9 at 16:52
• I already adapted the code. In my comment, I already specified the format of my hbonds.analysis.TIP3.add.txt. Is the format correct ? Apr 9 at 17:03
• Ask this as a separate question, and I can try to fix this Apr 9 at 18:20
• Hi, thanks for your code. I also tried to implement it to my systems. After having the .txt file containing the information on the HBs, I tried to run the second part with some modifications. At the end I had an array of NaNs because the hbondMatrix is empty. What type of information is fullTraj supposed to carry? To me it seems like the hbondMatrix is not being updated so it is always zero. Apr 10 at 7:55

Do you have any opinions on these codes?

does not seem like a useful question to me. My opinion might be that those codes are the most beautiful thing ever written, but I might also be someone with zero knowledge of Python.

Having said that, I'd optimize this calculation using what I know about the data structure.

Suppose I analyzed your trajectory and wanted to tell you when each relevant hydrogen was bonded and when it wasn't. I might start with a Numpy ndarray with nHyd rows, one for each hydrogen; each row would have nFrame entries, and each entry would be either 1 if the hydrogen were bonded or 0 if the hydrogen weren't.

But if you know for sure that each entry takes only 0 or 1, wouldn't it be far more efficient to encode the array as a series of flip frames? After all, suppose I give you this list for hydrogen number 312:

[0, 10, 42, 145, 162, 221]

with the key that the hydrogen was H-bonded from frames 0 to 9, and then from frames 42 to 144, and then ... I would have collapsed 300+ 1s and 0s into only six integers!

Hence:

flip_frames = [[] for _ in range(nHyd)]
frame_bonded = np.zeros(nHyd)

for nFrame, frame in enumerate(traj):
sorted_list_of_bonded_H = GetHBonds(frame)
frame_change = -1 * frame_bonded
# frame_change[indH] = -1 if bonded last frame and 0 if not
for indH in sorted_list_of_bonded_H:
frame_change[indH] += 1
# frame_change [indH] = 0 if bonded last frame and this frame, etc.
for changed_indH in np.nonzero(frame_change):
flip_frames[changed_indH].append(nFrame)
frame_bonded += frame_change

# now how to reconstruct:

from itertools import batched
def FrameHBondHistory(indH, nFrames):
HBond_History = np.zeros(nFrames)
indH_flip_frames = flip_frames[indH][:] # copy value, not ref
if (len(indH_flip_frames) % 2):
indH_flip_frames.append(nFrames) # if odd number of flips, add last frame
for start_frame, end_frame in batched(indH_flip_frames, 2):
HBond_History[start_frame:end_frame] = 1.
return HBond_History


Now you can process the trajectory efficiently (i.e. frame-by-frame), with a compact representation of the information that you can reprocess into the hydrogen bond history of any one hydrogen.

This lets you then calculate the autocorrelation for a single hydrogen's history any way you like -- my recommendation is to use the Fourier transform approach enabled by the Wiener-Khinchin theorem, coded (for example) as the function autocorr_func_1d from Dan Foreman-Mackey's blog.

EDIT:

There's even a very clean way to calculate the correlation integral directly from a flip-frames representation:

from itertools import batched
def HBondCorrSum(indH_flip_frames, nFrames):
CorrSum = np.zeros(nFrames)
if (len(indH_flip_frames) % 2):
indH_flip_frames.append(nFrames) # if odd number of flips, add last frame
HBonded_Intervals = batched(indH_flip_frames, 2)
for nInt, (Int1_start, Int1_end) in enumerate(HBonded_Intervals):
for Int0_start, Int0_end in HBonded_Intervals[:nInt]:
min_lag = Int1_start - Int0_end
max_lag = Int1_end - Int0_start
base_range = np.array(range(max_lag - min_lag + 1))
lag_range = base_range + min_lag
lag_counts = np.minimum(base_range, base_range[::-1]) + 1
nonneg_idx = 0 - min(min_lag, 0)
return CorrSum


The array CorrSum then contains the values:

$$\mathrm{CorrSum}_i[\tau] = \sum_{t = 0}^{t_{max}-\tau} c_i(t) c_i(t + \tau)$$

and if you look carefully, this function can even be implemented inline during a run -- simply checking if a "H-bonded interval" has finished on the current timestep and then updating the correlation sum based on all previous intervals recorded, fairly cheaply.

• Hi, thanks for your answer. Your code has been running for days now, is it normal? Apr 12 at 6:30
• I have no way of knowing. It is your responsibility to understand, test, and debug the code that you use to analyse your data. Apr 12 at 9:08

I noticed you've selected Hemanth's answer. However, I've devised my code differently and would like to verify if it generates the same results as yours. To do so, I need a reference. Could you please try to run my code with the input file, where the first column represents frame numbers and the second column represents the hydrogen bond values between two residues for example? Let me know the outcome once you've tested it.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

"""Reads bond data from the file."""
frame_numbers = []
bond_values = []
with open(file_name, 'r') as file:
for line in file:
frame, bond = map(float, line.strip().split())
frame_numbers.append(frame)
bond_values.append(bond)
return np.array(frame_numbers), np.array(bond_values)

def autocorrelation_function(frame_numbers, bond_values, time_interval):
"""
Calculates the autocorrelation function.

Autocorrelation Equation:
Autocorrelation = (1 / N) * âˆ‘[x(t) * x(t + time_interval)]

Parameters:
frame_numbers (array-like): Array of frame numbers.
bond_values (array-like): Array of bond values.
time_interval (int): Lag between observations.

Returns:
float: Autocorrelation value.
"""
autocorrelation = np.sum(bond_values[:-time_interval] * bond_values[time_interval:]) / len(frame_numbers)
return autocorrelation

def third_order_exponential(x, a, b, c, d):
"""Third order exponential function."""
return a * np.exp(-b * x) + c * np.exp(-d * x)

def main():
file_name = "hbonds2.dat"

# Determine max_value for max_lag
max_value = int(frame_numbers[-1]) + 125

# Choose a suitable time interval
time_interval = 1

best_autocorrelation = None
best_max_lag = None

# Initialize previous autocorrelation value
prev_autocorrelation = None

# Iterate over different max_lag values
tolerance = 0.001  # Adjust as needed
for max_lag in range(100, max_value, 100):  # Vary max_lag from 100 to maximum frame number + 100
lags = np.arange(1, max_lag)
autocorrelation_values = np.array([autocorrelation_function(frame_numbers, bond_values, lag) for lag in lags])

# Check if autocorrelation starts decreasing consistently
if prev_autocorrelation is not None:
min_length = min(len(prev_autocorrelation), len(autocorrelation_values))
if np.all(autocorrelation_values[:min_length] <= prev_autocorrelation[:min_length] + tolerance):
# If autocorrelation values decrease consistently, select the previous max_lag
best_max_lag = max_lag - 100
best_autocorrelation = autocorrelation_values[:min_length]
else:
# If autocorrelation values start increasing or stay constant, break the loop
break

# Update previous autocorrelation value
prev_autocorrelation = autocorrelation_values

if best_max_lag is None:
print("No convergence found within the specified range of max_lag.")
return

# Fit the autocorrelation function with a third-order exponential using the best max_lag
lags = np.arange(1, best_max_lag)
popt, pcov = curve_fit(third_order_exponential, lags, best_autocorrelation)

# Plot the results
plt.plot(lags, best_autocorrelation / best_autocorrelation[1], 'b-', label='Autocorrelation')  # Normalize by the first value
plt.plot(lags, third_order_exponential(lags, *popt) / best_autocorrelation[1], 'r--', label='Fitted Curve')  # Normalize by the first value
plt.xlabel('Time Lag')
plt.ylabel('Normalized Autocorrelation')
plt.legend()
plt.title('Hydrogen Bond Autocorrelation with Fitted Curve')

# Show the plot
plt.show()

if __name__ == "__main__":
main()

• Before running your code, what do you mean by bond values? I have a file that contains: frame number, donor id, hydrogen id, acceptor id, distance and angle. Apr 15 at 7:08
• Then, It is different from yours. Bond values mean if the bond exist, it will be 1, if the bond doesn't exist, it will be 0. You used the script in the answer to extract the distance and the angle? Apr 15 at 8:29
• I used the script in the answer to create a .txt file with this information. I could make another script adapted to your code, but I don't understand it very well. I mean, if the bond exists the bond value is 1 and 0 otherwise, but how to know which atoms are hydrogen bonded if in your definition there is no atom id? Or do you want to compute the autocorrelation function for a certain pair of atoms only? Apr 18 at 7:34
• @horlust. Yes, as you said I want to compute the autocorrelation function for a certain pair of atoms only. In my code, It is not important which atoms are forming the hydrogen bond, but I am interested in the presence of this hydrogen bond between two residues or not, as well you can specify the atoms. For example, by extracting only the hydrogen bond between two certain atoms, and therefore you have 0 for no bond and 1 for hydrogen bond between these two atoms. Apr 23 at 11:01