8
$\begingroup$

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.

$\endgroup$

3 Answers 3

4
$\begingroup$

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:
    fullData    =   hbondsObject.readlines()
    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.

$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$ Commented Apr 9 at 16:49
  • $\begingroup$ 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 $\endgroup$ Commented Apr 9 at 16:52
  • $\begingroup$ I already adapted the code. In my comment, I already specified the format of my hbonds.analysis.TIP3.add.txt. Is the format correct ? $\endgroup$ Commented Apr 9 at 17:03
  • $\begingroup$ Ask this as a separate question, and I can try to fix this $\endgroup$ Commented Apr 9 at 18:20
  • 1
    $\begingroup$ 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. $\endgroup$
    – horlust
    Commented Apr 10 at 7:55
2
$\begingroup$

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)

with MethodToReadTrajectory as traj:
    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)
            np.add.at(CorrSum, lag_range[nonneg_idx:], lag_counts[nonneg_idx:])
    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.

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

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

def read_bonds_data(file_name):
    """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"
    frame_numbers, bond_values = read_bonds_data(file_name)

    # 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()
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – horlust
    Commented Apr 15 at 7:08
  • $\begingroup$ 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? $\endgroup$ Commented Apr 15 at 8:29
  • $\begingroup$ 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? $\endgroup$
    – horlust
    Commented Apr 18 at 7:34
  • $\begingroup$ @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. $\endgroup$ Commented Apr 23 at 11:01

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .