Regarding writing a code to calculate hydrogen bond autocorrelation function

[cy16f01]

Hello, I am trying to write a code using mdanalysis library to write a code to calculate HB autocorrelation function.... I tried the following in the jupyter notebook (with python 2.7) as:-

from MDAnalysis.analysis import hbonds
import matplotlib.pyplot as plt
H = u.select_atoms('name H1')
O = u.select_atoms('name OW')
N = u.select_atoms('name N1')
hb_ac = hbonds.HydrogenBondAutoCorrel(u, acceptors = u.atoms.OW,hydrogens = u.atoms.H1, donors = u.atoms.N1,bond_type = 'continuous',sample_time = 0,nruns = 1,  nsamples = 2)
hb_ac.run()
hb_ac.solve() 
tau = hb_ac.solution['tau']
time = hb_ac.solution['time']
results = hb_ac.solution['results']
estimate = hb_ac.solution['estimate']
plt.plot(time, results, 'ro')
plt.plot(time, estimate)
plt.show()

but i am getting error as"- File "", line 10 hb_ac.solve() ^ SyntaxError: invalid syntax

How can i solve this.. i am a newbie here....

Is this the full code to calculate continuous HB autocorrelation function.....

Any suggestions are helpful...

Thank you...

[kain88-de]

Can you upload the notebook? Best unmodified.

[orbeckst]

@cy16f01 without the actual input files and more information it is difficult to help. I tried to reproduce your example with some of the test trajectories and I can run the analysis (although the test data are not long enough to do some real calculations). I am also not an expert on this code ( @richardjgowers is the one who wrote it). Perhaps the following is useful:

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import TRZ, TRZ_psf
from MDAnalysis.analysis import hbonds
import matplotlib.pyplot as plt

u = mda.Universe(TRZ_psf, TRZ)

H = u.select_atoms('name Hn')
O = u.select_atoms('name O')
N = u.select_atoms('name N')

# match donors and hydrogens
# (not really necessary here because in the test trajectory H and N are 
# already matched but this is a useful trick in more general cases)
NH = N.bonds.atomgroup_intersection(H)
# now NH.atom1 are the N that are matched to the H
# and NH.atom2 are the H that are matched to the N

# use atom groups for acceptors, donors, hydrogens
hb_ac = hbonds.HydrogenBondAutoCorrel(u, acceptors=O, hydrogens=NH.atom2, donors=NH.atom1, 
             bond_type='continuous', sample_time=0.06)
hb_ac.run()
print(hb_ac.solution['results'])
# [ 1.          0.92668623  0.83137828  0.74486804  0.67741936  0.60263932]

# the next line does not work properly with the limited input data
hb_ac.solve()

After all this, we get

print(hb_ac.solution)

in the results dict

{'results': array([ 1.        ,  0.92668623,  0.83137828,  0.74486804,  0.67741936,
        0.60263932], dtype=float32), 'time': array([ 0.  ,  0.01,  0.02,  0.03,  0.04,  0.05], dtype=float32), 'fit': array([ 0.02759426,  0.1023183 ,  0.10231667]), 'tau': 0.10231671451853262, 'estimate': array([ 1.        ,  0.90688848,  0.82244676,  0.74586761,  0.67641872,
        0.61343646], dtype=float32), 'infodict': {'fvec': array([ 0.        ,  0.0197977 ,  0.00893148, -0.00099953,  0.00100063,
       -0.01079709]), 'nfev': 65, 'fjac': array([[ -4.66872576e+00,  -1.29349952e-01,   3.48206897e-05,
         -4.45182582e-01,  -5.38308029e-01,  -6.10231751e-01],
       [ -1.32501715e-01,   1.01019597e-06,  -4.36898835e-06,
          3.49412235e-01,  -6.53572132e-02,  -6.18111982e-01],
       [  3.56691364e-05,  -1.04085755e-07,  -3.71252212e-07,
          4.11899991e-02,  -7.91960400e-01,   4.49856122e-01]]), 'ipvt': array([3, 2, 1], dtype=int32), 'qtf': array([  4.46978842e-07,  -2.00122487e-02,   2.23453496e-03])}, 'mesg': 'Both actual and predicted relative reductions in the sum of squares\n  are at most 0.000000', 'ier': 1}

[cy16f01]

I am new to this notebook,,can u kindly help me how to upload the notebook unmodified..?? The above code is as the same...

My doubt is that.... I have done a MD simulation of glycine and water mixture as an example in gromacs, and now i need to calculate the hydrogen bond auto correlation function (HBond lifetime-intermittent and continuous). I searched the program/code for HBond autocorrelation function of the following links, but i am not able to find out which is the right code to calculate the HBond auto correlation function - intermittent and continuous. 1] http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/waterdynamics.html#id2 2] http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/hbond_autocorrel.html 3] http://www.mdanalysis.org/mdanalysis/_modules/MDAnalysis/analysis/hbonds/hbond_autocorrel.html which code in the following pages are the correct code in MDAnalysis library to write the code to calculate HBond auto correlation function - intermittent and continuous. ... Can anybody help me regarding this....??

Thank you...

[cy16f01]

Thank you @orbeckst.. how can i contact/put a query to @richardjgowers thn..? if i need to write the code..?? Regarding limited input data:- what more data is needed for this (to calculate HBond autocorrelation function)..?? , can u kindly explain/elaborate..?? Thank you..

[orbeckst]

If @richardjgowers finds time to reply he will likely do it here; whenever you add a GitHub username to an issue report then that person is notified.

You could also ask on the user mailing list, which is read by the developers, too. http://users.mdanalysis.org

Am Aug 2, 2017 um 23:02 schrieb cy16f01 notifications@github.com:

Thank you @orbeckst.. how can i contact/put a query to @richardjgowers thn..? if i need to write the code..?? Thank you..

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[orbeckst]

I searched the program/code for HBond autocorrelation function of the following links, but i am not able to find out which is the right code to calculate the HBond auto correlation function - intermittent and continuous.

From your question it is not clear exactly what you want to calculate. What mathematical quantity do you want to calculate?

Then read the descriptions and references listed in

• http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/hbond_autocorrel.html
• http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/waterdynamics.html

and compare to what you need.

EDIT: Once you know which code you want to use, we can help with what you need as input etc.

[orbeckst]

Regarding limited input data:- what more data is needed for this (to calculate HBond autocorrelation function)..?? , can u kindly explain/elaborate..??

The test data that I used in the example above does not contain a sufficient number of trajectory frames to produce enough data points for a good fit. Hopefully, when you do this on your data you will have more frames.

[cy16f01]

Regarding From your question it is not clear exactly what you want to calculate. What mathematical quantity do you want to calculate? @orbeckst 1] I have done a MD simulation of glycine and water mixture as an example in gromacs, and now i need to calculate the hydrogen bond auto correlation function (intermittent and continuous) . For HBond autocorrelation function, which of the following links is the correct one..?? which is the right code to calculate the HBond auto correlation function (intermittent and continuous). What are the differences between these codes/what these two codes are going to calculate exactly, are there any differences in results which i get using these two codes ..?? http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/hbond_autocorrel.html http://www.mdanalysis.org/mdanalysis/documentation_pages/analysis/waterdynamics.html since both the codes are very different, and i am finding difficult to understand the code...

2] what are the input files that are needed to be given, (a brief elaboration with basics would be highly beneficial since i am a novice/newbie here.. ) like, between which donors and acceptors i need to calculate, what does the meaning and significance of this statement:- "hydrogens and donors selections must be aligned, that is hydrogens[0] and donors[0] must represent a bonded pair."..??

[cy16f01]

and in this code:-

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import HydrogenBondLifetimes as HBL
import matplotlib.pyplot as plt

u = MDAnalysis.Universe('md.gro','md.trr')  # .gro and .trr file
selection1 ="byres name OH2 and sphzone 6.0 protein and resid 42"
selection2 = "resid 38"
HBL_analysis = HBL(u, selection1, selection2, 0, 200, 30)
HBL_analysis.run()
time = 0
#now we print the data ready to plot. The first two columns are the HBLc vs t
#plot and the second two columns are the HBLi vs t graph
for HBLc, HBLi in HBL_analysis.timeseries:
    print("{time} {HBLc} {time} {HBLi}".format(time=time, HBLc=HBLc, HBLi=HBLi))
    time += 1

#we can also plot our data
plt.figure(1,figsize=(18, 6))

#HBL continuos
plt.subplot(121)
plt.xlabel('time')
plt.ylabel('HBLc')
plt.title('HBL Continuos')
plt.plot(range(0,time),[column[0] for column in HBL_analysis.timeseries])

the code is calculating for only between two water molecules for the whole frame, but i need to calculate for the whole set of water molecules for all the time frame.....how can i do tht.??

[orbeckst]

1] I have done a MD simulation of glycine and water mixture as an example in gromacs, and now i need to calculate the hydrogen bond auto correlation function (intermittent and continuous) . For HBond autocorrelation function, which of the following links is the correct one..?? which is the right code to calculate the HBond auto correlation function (intermittent and continuous). What are the differences between these codes/what these two codes are going to calculate exactly, are there any differences in results which i get using these two codes ..??

Read the associated papers (1) Gowers 2015 and (2) Rapaport 1983 together with Araya-Secchi 2014 and Milischuk 2011 to understand what the code is trying achieve. Then pick the one that best fits your purpose. There is no substitute for reading the primary literature.

[orbeckst]

what are the input files that are needed to be given, (a brief elaboration with basics would be highly beneficial since i am a novice/newbie here.. )

You need topology and trajectory files to create a universe and you need to create AtomGroups for hydrogens, donors, and acceptors by doing selections. All of this is explained in the MDAnalysis Tutorial.

like, between which donors and acceptors i need to calculate,

This is up to you – it depends on your scientific question. Read the original papers to understand how this was used originally.

We can help with technical questions such as "How do I select the water oxygen atoms". You have to do the science yourself.

what does the meaning and significance of this statement:- "hydrogens and donors selections must be aligned, that is hydrogens[0] and donors[0] must represent a bonded pair."..??

hydrogens and donors are two AtomGroups. You can think of them as lists. The first atom in hydrogens, namely hydrogens[0] must be connected to the first atom donors[0] through a bond, i.e., it must be the hydrogen that comes from that donor. It is the user's responsibility to construct these AtomGroups in such a way that this is true. In my https://github.com/MDAnalysis/mdanalysis/issues/1563#issuecomment-319540775 I showed you a trick how one can do this.

[cy16f01]

@orbeckst, @richardjgowers, @dotsdl, @arose ... I have ran the simulation in GROMACS for 10ns {integrator= md; leap-frog integrator,nsteps= 10000000 ; 0.001 * 10000000 = 10000 ps (10 ns),dt = 0.001; 1e+07 fs} and saving coordinates for every 0.5 ps (nstxout = 500; save coordinates every 0.5 ps) This is my .gro (final .gro file after the md run) file which contains glycine and water molecules... Glycine-Water XXX 1GLYN N 1 1.554 0.192 1.237 0.2337 0.7678 -0.4403 1GLYN HT1 2 1.595 0.113 1.190 1.2957 1.3150 -0.4663 1GLYN HT2 3 1.567 0.275 1.179 0.2555 1.4287 0.4977 1GLYN CA 4 1.596 0.220 1.375 -0.1659 0.2593 -0.2142 1GLYN HA1 5 1.690 0.279 1.386 0.3912 -0.2223 -2.2182 1GLYN HA2 6 1.527 0.296 1.419 -0.2036 -0.1122 0.3771 1GLYN C 7 1.590 0.109 1.479 -0.4882 0.5503 0.0814 1GLYN OT1 8 1.489 0.061 1.528 -0.1025 -0.2671 0.0940 1GLYN OT2 9 1.711 0.048 1.513 -0.4918 0.6753 0.3137 1GLYN HO2 10 1.698 -0.035 1.560 -0.8295 -0.2734 -1.4207 2SOL OW 11 1.974 1.991 1.614 0.2492 -0.3487 0.3687 2SOL HW1 12 2.040 1.918 1.594 -0.6780 -1.3661 0.9449 2SOL HW2 13 1.884 1.951 1.630 -0.5310 1.0328 -0.4632 3SOL OW 14 0.331 0.052 1.858 -0.2124 0.4141 -0.8687 3SOL HW1 15 0.310 -0.007 1.936 0.1912 -0.3242 -1.3072 3SOL HW2 16 0.293 0.011 1.775 0.6684 0.4113 -1.2785 4SOL OW 17 0.482 1.719 1.590 -0.5622 -0.3605 0.2368 4SOL HW1 18 0.478 1.632 1.541 0.8012 -1.4225 1.9691 4SOL HW2 19 0.444 1.707 1.682 0.8493 0.5546 0.9465 5SOL OW 20 0.168 1.334 1.467 0.1289 -0.3488 0.2064 5SOL HW1 21 0.149 1.375 1.556 3.3113 1.1012 0.2810 5SOL HW2 22 0.167 1.234 1.475 1.0811 -0.2295 2.0747 . . . 1] So which code should i follow in order to calculate the HBond auto correlation function(intermittent and continuous). 2] I want to calculate HBond (say, between two water molecules) or (say, between glycine atom and other water molecule) for the whole set/full time frame (30 ns), what are the modifications tht i need to do..?? can u kindly give/share a full sample code in order to exactly understand wht input i need to give and output will be getting,by running the code. etc.,

Thank you...

[orbeckst]

For (2) hydrogen bonds look at http://www.mdanalysis.org/docs/documentation_pages/analysis/hbond_analysis.html

For (1) I would try http://www.mdanalysis.org/docs/documentation_pages/analysis/hbond_autocorrel.html

[orbeckst]

It seems that your questions have moved to other issues. I'm closing this one but if you have more questions or data to show feel free to add to this issue and we can reopen it.