I want to calculate water residence time near N+ as a fixed ionic group in my polymer system. Is there any direct way to calculate average residence time of water near N+? If not, could you please guide me how to do it? I found some clues but they were a bit confusing.

Gromacs does not have such facility to calculate residence time directly. Typically, one needs to calculate autocorrelation function and integrate it to get residence time. Gromacs calculates autocorrelation function for hydrogen bonds, which is almost similar to residence autocorrelation. You might want to tweak gmx hbond distance criteria to get autocorrelation. I would say, use first peak of rdf as rcut and ignore angle criteria by using angle cutoff 180.

Thank you very much for your reply. As I understand I can use gmx hbond with some modifications to get ACF between N+ and OW. Then integrate ACF to get residence time. I used:

Gromacs first finds donor-acceptor-hydrogen triplet for hydrogen bond analysis. So mere OW-N will not give you hydrogen bond. I would say use entire water molecule instead of OW.

You use double exponential for more accuracy, is it true? And, to be honest, I use python but I am beginner and I try to apply R(t) to my ACF and integrate it to calculate tau. That would be very kind of you if you help me with this.

I tried to write a script for fitting and integrating. But, I could not define double exponential function and its bound. I attached file here. Could you please correct it?

import numpy as np
import csv
import matplotlib.pylab as plt
from scipy.optimize import curve_fit
from scipy import integrate as intg
with open("residence.csv",'r') as i: #open a file in directory of this script for reading
rawdata = list(csv.reader(i,delimiter=",")) #make a list of data in file
exampledata = np.array(rawdata[0:],dtype=np.float) #convert to data array
xdata = exampledata[:,1]
ydata = exampledata[:,0]
plt.figure(1,dpi=120)
plt.yscale('linear')# plt.xscale('linear')
plt.xlim(0,1) # x axis start from 0 to 1( C(t))
plt.ylim(0,5000) # y axis start from 0 to 5000 (time ps)
plt.title("Example Data")
plt.xlabel(rawdata[0][0])
plt.ylabel(rawdata[0][1])
plt.plot(xdata,ydata,label="Experimental Data")
def func(x,b):
return a0 * np.exp(-(1/b) * x) + a1 #how to define second exp?
a0 = 3500
a1 = 0
funcdata = func(xdata,0.1)
plt.plot(xdata,funcdata,label="Model")
plt.legend()
popt, pcov = curve_fit(func,xdata,ydata,bounds=(0,1))
perr = np.sqrt(np.diag(pcov))
TotalInt = intg.trapz(ydata,xdata) #Compute numerical integral
TotalInt_func = intg.quad(func,0,1, args=(1))[0] #Compute integral of function