Hydrogen Bond Autocorrelation Functions

User discussions
1

GROMACS version: 2024
GROMACS modification: No

Hello everyone!

I recently found out that GROMACS can compute the forward and backward constants of the theory of Luzar and Chandler. I have a simulation box containing 1414 molecules of water (TIP4P) and I did gmx hbond-legacy -f file.xtc -s file.tpr -ac ac.xvg. The aucotorrelation function is ok, I compared this with a Python code that I have for this. However, the constant values are weird:

image

This is for 298 K. According to Luzar and Chandler, k = 0.6 and k’ = 0.9. Why don’t I have these values? The simulation length is of 50 ps, maybe it should be longer?

Also, where can I find the code for the calculation of these constants?

Edit: I am also calculating it for simulations of 500 ps, but still get weird results (negative values for example).

Thank you very much!
Best,
Maria

2

There are a few factors that could be relevant.

  1. Do you know the expected values for the water model you are using?
  2. How often are you writing the simulation trajectory output? In Resolving the hydrogen bond dynamics conundrum | The Journal of Chemical Physics | AIP Publishing Luzar writes “This is why it is important to use sampling intervals well below 10 fs.” I haven’t checked if his analyses are the same, but it’s quite possible that you should have a very high trajectory output frequency. He also writes “Therefore, to obtain good statistics for long time tails, one must collect data over hundreds of picoseconds (hundreds of nanoseconds at low T).” So, 500 ps sounds like a reasonable lower limit.
3

My .mdp file for production is:

; Run control
integrator = md
tinit = 0
dt = 0.002
nsteps = 250000 ; 500 ps
nstcomm = 100

; Output control
nstxout = 0
nstvout = 50
nstfout = 0
nstlog = 50
nstenergy = 50
nstxout-compressed = 50

I will try to run a 5 ns simulation and improve the trajectory output frequency.

4

You are writing coordinates and velocities every 100 fs. Luzar said that sampling intervals well below 10 fs should be used.

I’m usually suggesting to write coordinates (and velocities, if any) at low frequencies (every 10-200 ps depending on use case). But in this case, it seems like it’s important to write frequently.

Hopefully it will help.

5

Here’s a concise and clear version of your message:

Hi,

I saved the coordinates every 10 fs, but I’m still encountering some issues.

image
image408×118 28.7 KB

According to Luzar, the expected values are approximately k=0.6 and k’ = 1.0. The .mdp file is something like:

; Run control
integrator = md
tinit = 0
dt = 0.001
nsteps = 25000 ; 25 ps
nstcomm = 100

; Output control
nstxout = 10
nstvout = 10
nstfout = 0
nstlog = 500
nstenergy = 500
nstxout-compressed = 10

Maybe I should run longer simulations? Also, do you know where I can find the code for the calculation of c(t) and n(t)?

6

The recommendations by Luzar was to sample well below every 10 fs and to run for several hundreds of ps. Does it make any difference with:

; Run control
integrator = md
tinit = 0
dt = 0.001
nsteps = 500000 ; 500 ps
nstcomm = 100

; Output control
nstxout = 5
; nstvout = 10
nstfout = 0
nstlog = 500
nstenergy = 500
; nstxout-compressed = 10

?

I also suspect that it might be worth using an uncompressed trr file in this case (as per the settings above). However, keep in mind that the output files will be large.

As far as I can see, the analyses can be found in the file src/gromacs/gmxana/gmx_hbond.cpp, see the analyse_corr function.

7

Still incorrect (k’ = 0.0 from what I’ve calculated). Maybe it’s because I’m using gmx hbond-legacy? Anyways, I’m taking a look at the code to see what I can do :)

8

I computed k and k’ on Python and got the same results as Luzar’s. I think there is something wrong with this code on GROMACS.

9

Please, report an issue on gitlab (Issues · GROMACS / GROMACS · GitLab), preferably with detailed instructions how to reproduce the problem.

10

I reported the error :) Btw, do you know where in the GROMACS code the n(t) computation is? I wanted to take a look.