How to calculate residence time of water near a specific (ionic) group?

User discussions
1

Hi,

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.

Thanks,in advance

2

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.

3

Dear Masrul,

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:

gmx hbond -f npt2.xtc -s npt2.tpr -n N-OW.ndx -ac residenceacf.xvg -a 180 -r 0.46 -b 90000`` -e 100000

in witch 0.46 is frim first rdf peak of N±Ow. But, after applying this command gromacs shows me this error:

Found 0 donors and 3520 acceptors
Making hbmap structure…done.
No Donors found

Is my procedure true?if yes, How can I fix that error?

Thanks again

4

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.

5

Dear Masrul,

Thank you very much. I got ACF from your suggested way. But, one last question. After termination of process The terminal shows line like this:

Hydrogen bond thermodynamics at T = 298.15 K
One-way 0.030 33.158 13.207
Integral 0.010 103.152 16.021
Relaxation 0.045 22.125 12.205

Does any of this numbers represent residence time or I have to fit ACF to an exponential function to get tau?

Thanks again

6

Dear Rezayani,

I never used ACF from gromacs, so I have no idea. I would fit ACF into

R(t) = Ae^{-t/\tau_1} ~+~ Be^{-t/\tau_2}

And integrate to get residence time,
\tau=\int^\infty_0 R(t) dt

You can use python for fitting and integrating. Let me know if you need help for python scripting.

Regards
Masrul

7

Dear Masrul,

Thanks again for your help

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.

Best,
Rezayani

8

Dear Masrul

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

Thanks alot

9

Can you upload your xvg/csv file here?

10

Hi,

this is my original file from Gromacs:
residenceacf.xvg (30.2 KB)

First two columns are related to ACF. I got ACF from gmx h_bond. Between N+ and water.

CSV files are not authorized. But I change my csv file format (I used for Python) to dat file here:
residence.dat (8.7 KB)