# Hydrogen Bond Autocorrelation Functions **URL:** \ **Category:** User discussions\ **Created:** [October 1, 2024, 5:55am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263 "2024-10-01T05:55:39Z")\ **Posts on this page:** 10\ **Page:** 1
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [October 1, 2024, 5:55am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/1 "2024-10-01T05:55:39Z")
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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/a/adc87eb85cf8542ec622c4b535fa6c0afbe0a955.png) 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 ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 1, 2024, 7:34am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/2 "2024-10-01T07:34:45Z")
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](http://dx.doi.org/10.1063/1.1320826) 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. ---
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [October 1, 2024, 8:48am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/3 "2024-10-01T08:48:50Z")
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. ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [October 1, 2024, 8:54am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/4 "2024-10-01T08:54:32Z")
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. ---
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [January 7, 2025, 7:03am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/5 "2025-01-07T07:03:19Z")
> [@MagnusL](#): > > writing coordinates and velocities every 100 fs 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](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/2X/a/abe7f82cd86650607ab696ccb46f7468587c0665.png) 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)? ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [January 7, 2025, 8:08am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/6 "2025-01-07T08:08:08Z")
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: ```auto ; 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. ---
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [January 9, 2025, 8:55am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/7 "2025-01-09T08:55:43Z")
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 :) ---
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [January 14, 2025, 8:52am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/8 "2025-01-14T08:52:55Z")
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. ---
**Author:** ![MagnusL](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/magnusl/32/2380_2.png) [@MagnusL](https://gromacs.bioexcel.eu/u/MagnusL)\ **Post date:** [January 14, 2025, 11:02am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/9 "2025-01-14T11:02:22Z")
Please, report an issue on gitlab ([Issues · GROMACS / GROMACS · GitLab](https://gitlab.com/gromacs/gromacs/-/issues)), preferably with detailed instructions how to reproduce the problem. ---
**Author:** ![horlust](https://avatars.discourse-cdn.com/v4/letter/h/b2d939/32.png) [@horlust](https://gromacs.bioexcel.eu/u/horlust)\ **Post date:** [February 4, 2025, 11:21am UTC](https://gromacs.bioexcel.eu/t/hydrogen-bond-autocorrelation-functions/10263/10 "2025-02-04T11:21:48Z")
I reported the error :) Btw, do you know where in the GROMACS code the n(t) computation is? I wanted to take a look.