CHARMM Development Project Forums CHARMM Community Script Archive hbonds.inp

*FILENAME: hbonds.inp
*PURPOSE: compute static and time-dependent hydrogen bonds from trajectory
*AUTHOR: Lennart Nilsson, Karolinska Institutet, October 2003
*
!unix environment variable CHM_HOME points to CHARMM installation directory
read rtf card name $CHM_HOME/toppar/top_all22_prot.inp
read para card name $CHM_HOME/toppar/par_all22_prot.inp
! we assume there are three segments in this psf
! prot = a protein
! lig = a bound ligand
! wat = solvent waer
read psf card name my_psf.psf
read coor card name my_structure.crd
! hydrogen bonding patterns in initial structure
! the "coor hbond" depends on DONOr/ACCEptor information in PSF to
! identify possible hbonding partners; see corman.doc for details
! we use strictly geometric criteria on the H...A distance, and optionally
! on the D-H...A angle to define the hydrogen bond
! first look at the backbone-backbone hydrogen bonds
coor hbond verbose select type HN .or. type O end select type HN .or. type O end
! protein-solvent !!NB your psf MUST have DONOr/ACCEptor specifications
! for the water molecules for this to work;
! if any of the following lines are missing from RESI TIP3
! in your RTF (and hence in the PSF), add missing bits to RTF and
! regenerate the PSF
! ACCEPTOR OH2
! DONOR H1 OH2
! DONOR H2 OH2
coor hbond verbose select segid prot end select segid wat end

! now get the timedependent hydrogen bonds; the output can be voluminous, especially
! with the verbose option
open unit 51 read unform name my_traj1.cor
open unit 52 read unform name my_traj2.cor
! the SKIP 2500 (=5ps) means that for a hydrogen bond to be counted as broken
! in the calculation of average lifetime it has to be broken for > 5ps (otherwise
! we are unlikely to detect the break)

! backbone-backbone
coor hbond verbose select type HN .or. type O end select type HN .or. type O end -
firstu 51 nunit 2 begin 10500 skip 2500

! note that the trajectory is not closed between these commands, and can be reused
! protein - solvent
coor hbond verbose select segid prot end select segid wat end -
firstu 51 nunit 2 begin 10500 skip 2500

! protein - water - ligand
coor hbond verbose select segid prot end select segid lig end -
bridge tip3 firstu 51 nunit 2 begin 10500 skip 2500

Hi.
Which column in the psf specifies donor/acceptor status?
Thanks.

It's not a column, but 2 sections, a donor list and an acceptor list; the list is a sequence of integers, which is the atom number index defined in the first section of the PSF with the charges, etc.

Write out a formatted PSF file and have a look for yourself ...

Just to make sure the explanation made is understood, i 'll use an example.

so in the example below, one always read along a row with a donor index followed by an acceptor index, donor then acceptor kind of fashion.

so atom 1 and 5, 2 and 1, 3 and 1 are examples of donor/acceptor pairs?

30279 !NBOND: bonds
1 5 2 1 3 1 4 1
5 6 7 5 7 8 7 9
10 7 10 11 10 12 13 10
13 14 13 15 16 13 16 17
16 18 19 16 19 20 19 21
19 22 23 5 23 25 24 23
25 26 25 27 27 28 29 27
29 30 31 29 31 32 31 33
31 34 35 29 35 36 35 37
35 38 39 27 39 41 40 39
41 42 41 43 43 44 45 43
45 46 45 47 48 45 49 48
49 50 51 49 51 52 53 51

NO, this is not correct (and why are you looking at the list of bonds?).

Add DONOr and ACCEptor statements as needed to your RTF and regenerate the PSF - do not play directly with the PSF.

I am looking just to make sure i have the correct output.

I ...
read parm
read rtf
read psf
read coor

write new.psf
write new.coor

This is how i regenerate psf. So i want to know if the new psf is what i want. Hence i need to know where are the donor acceptor sections.

Help.

Check for DONOr and ACCEptor statements in the RESIdues of interest in your RTF. What is so difficult about this?

The HBOND sections of the PSF begin with lines containing the strings "!NDON: donors" and likewise "!NACC: acceptors", and are listed in pairs as donor,donor-antecedent or acceptor,acceptor-antecedent. It is, however, mostly irrelevant as the HBOND term of the potential isn't used except for analysis, and for that it is trivial to use the DONOR and ACCEPTOR commands to define them on the fly. You don't need a new PSF for that.

Thank you all. I know better now.

Hello!!

I just have a questions about the output of the script. I have a implicit solvent trayectory of 1fs timestep, so I commented the water H-bond part, and I changed the coor hbond setup as follows:

!backbone-backbone
coor hbond verbose select type HN .or. type O end select type HN .or. type O end -
first 51 nunit 1 begin 100 skip 100 stop 200000

Part of the output is :

READING TRAJECTORY FROM UNIT 51
NUMBER OF COORDINATE SETS IN FILE: 2000
NUMBER OF PREVIOUS DYNAMICS STEPS: 100
FREQUENCY FOR SAVING COORDINATES: 100
NUMBER OF STEPS FOR CREATION RUN: 200000
...
I-atom J-atom (Bridge) Lifetime Endtime
PROT 37 LYS HN - PROT 33 ALA O 0.60 1.00
PROT 43 LYS O - PROT 46 ALA HN 0.20 1.10
PROT 51 LEU HN - PROT 51 LEU O 0.20 0.40

I am not sure of understanding the 'lifetime' and 'endtime'. I saved structures every 100. I interpret the first line hbond (37-33) as existing during the time window 0.4-1.0, and in the second case (43-46) during 0.90-1.10, Is that correct?? I also noted cases as the third line were a Hbond is considered beetween C and N of the same residue. I read that one can define CUTHA = 180 to allow only linear Hbonds, nevertheless, the angles of the first structure (showed earlier in the output) range from 170 to 90 degrees, being the average angle ~160. How can i avoid strange angles without being too restrictive?

Thanks

Your interpretation of the lifetime and endtime columns is correct. I do not think that you have to worry about the angular dependence - in practice it works very well with just a distance dependent cutoff on the H-Acceptor distance (see J. Am. Chem. Soc., 1992. 114: p. 4028-4035). It may be that more distorted hydrogen bonds do actually occur in implicit solvent simulations (I have not looked at this, so it is just a speculation).
Note that with a 100 step resolution you may underestimate the hydrogen bond lifetimes; we normally use a resolution of 1-5 ps.

In determining prot-water-ligand bridges, i've specified 'tip3 .and. resi 136' instead of 'TIP3' to make it more specific. There is no error in the output but i cant be sure it is correct. Can someone confirm? The command are as followed:
coor hbond verbose cuthb @CUTHB cutha @CUTHA select ONE end select TWO end-
bridge tip3 .and. resi 136 firstu 51 nunit @N skip @SKIP

Incorrect, but you can achieve this effect by renaming tip3 136:
rename resn t136 select atom wat 136 * end
coor hbond ... bridge t136 ....

It is typically not necessary to define a cutoff for the angle.

hi,

i would like to know how to interpret the results of lifetime and endtime. as i read the post on hbonds, it is still not clear with the output i have.

######################################################

I-atom J-atom (Bridge) Lifetime Endtime
LEC 201 ASP OD1 - MIC 13 BOG HO4 60.00 70.00
LEC 201 ASP OD2 - MIC 13 BOG HO3 140.00 150.00
LEC 201 ASP OD2 - MIC 13 BOG HO3 120.00 280.00
LEC 201 ASP OD1 - MIC 13 BOG HO4 230.00 310.00
LEC 173 LEU O - MIC 20 BOG HO3 70.00 410.00
LEC 195 ASP OD2 - MIC 8 BOG HO4 90.00 410.00
LEC 195 ASP OD2 - MIC 8 BOG HO4 100.00 520.00
LEC 145 LYS HZ2 - MIC 16 BOG O4 90.00 740.00
LEC 195 ASP OD1 - MIC 8 BOG HO4 150.00 770.00
LEC 108 GLU OE2 - MIC 16 BOG HO2 90.00 800.00
LEC 129 ASP OD1 - MIC 16 BOG HO4 150.00 800.00
LEC 145 LYS HZ2 - MIC 16 BOG O2 190.00 800.00

I-atom J-atom (Bridge) Lifetime Endtime
LEC 100 TYR O - MIC 16 BOG O4 WAT1 7106 250.00 260.00
LEC 174 THR O - MIC 20 BOG O6 WAT1 6496 60.00 470.00
LEC 192 THR OG1 - MIC 8 BOG HO6 WAT1 5614 60.00 520.00

#########################################################

as seen in the second and third row the protein and lipid atoms are the same except their two different lifetime and endtimes. which one should i consider in this case.

in general:--

the endtime --- is it the actual residence time b/n those atoms?

the lifetime --- is it the duration of active hbond ?

thank you

In this section of the verbose-mode output the endtime is the timepoint at which this particular instance of the bond between these two atoms was broken. The lifetime is the duration of this particular instance of the hydrogen bond or bridge.
Thus the information in your second and third rows means that there was a hydrogen bond between these atoms in the two intervals (10,150) and (160,280).

You are probably better off looking at the more compact summary information at the end of the output.

Hello,

I also have questions of the output od hydrogen bond analysis even I've read the above posted but it's still unclear to me.

In protein-proetin output section, I have:

Analysis of hydrogen bonds using cutoff distance= 2.40 and angle= 999.00
I-atom J-atom r(H-A) [A] D-H..A [deg]
--------------------------------------------------------------------
PROT GLY 4 O - PROT LYS 56 HN 1.90 168.7
PROT GLY 4 O - PROT THR 57 HN 1.86 157.8
PROT TRP 6 HN - PROT THR 57 O 1.88 169.8
PROT TRP 6 O - PROT HSD 59 HN 1.99 163.5
PROT GLN 7 HN - PROT VAL 22 O 1.79 167.7
PROT GLN 7 O - PROT VAL 22 HN 1.89 170.1
PROT LEU 8 HN - PROT HSD 59 O 1.77 175.0
PROT THR 11 HN - PROT ILE 18 O 1.86 165.5
PROT THR 11 O - PROT ILE 18 HN 1.91 159.1

and the lifetime section is:

I-atom J-atom (Bridge) Lifetime Endtime
PROT 2 SER HN - PROT 2 SER O 5.00 10.00
PROT 2 SER O - PROT 2 SER HN 5.00 10.00
PROT 5 ILE O - PROT 24 VAL HN 5.00 10.00
PROT 17 VAL HN - PROT 17 VAL O 5.00 10.00
PROT 17 VAL O - PROT 17 VAL HN 5.00 10.00
PROT 24 VAL HN - PROT 5 ILE O 5.00 10.00
PROT 41 GLU O - PROT 45 PHE HN 5.00 10.00
PROT 45 PHE HN - PROT 41 GLU O 5.00 10.00

My questions are:
1. For lifetime section, Does the first line mean that the bridge between HN of SER2 and O of the same residue occurs within 5ps and it will be broken after 10 ps? I just want to make sure whether I'm correctly understand.

2. Why the lines of the first section (protein-protein output section in this case) are not match with the lines of the lifetime section? For example,

If the protein-protein section starts with

PROT GLY 4 O - PROT LYS 56 HN 1.90 168.7
PROT GLY 4 O - PROT THR 57 HN 1.86 157.8
PROT TRP 6 HN - PROT THR 57 O 1.88 169.8

The lifetime section should also start with:

I-atom J-atom (Bridge) Lifetime Endtime

PROT GLY 4 O - PROT LYS 56 HN ... ....
PROT GLY 4 O - PROT THR 57 HN ... ....
PROT TRP 6 HN - PROT THR 57 O ... ....

Thanks in advance.

1/ Yes, your understanding is correct.
2/ You are apparently looking at two different COOR HBOND outputs. The one you call protein-protein is from a static analysis of a single coordinate set; a similar output, but with hydrogen bond lifetime data, is printed from the analysis of a trajectory. Your other output, which you call lifetime section, is simply a chronological record of all hydrogen bond forming/breaking events in the trajectory, as such it is time-ordered; it is usually not necessary to look at this output. The summary information, giving the average occupancies and lifetimes is more informative and easier to grasp.

Hello,

Thank you Prof. Lennart. However I still have questions.

I've looked at the summary section and I've found something like this:

I-atom J-atom <occupancy> <time> # events
---------------------------------------------------------------------------
PROT SER 2 HN - PROT SER 2 O 0.44 18.3 6
PROT SER 2 O - PROT SER 2 HN 0.44 18.3 6
PROT SER 2 O - PROT GLY 4 HN 0.06 7.5 2
PROT SER 2 O - PROT ILE 5 HN 0.28 7.8 9
PROT SER 2 O Summary: 0.78 11.2 17
PROT GLY 4 HN - PROT SER 2 O 0.06 7.5 2
PROT GLY 4 O - PROT LYS 56 HN 1.00 250.0 1
PROT GLY 4 O - PROT THR 57 HN 0.98 122.5 2
PROT GLY 4 O Summary: 1.98 186.2 3
PROT ILE 5 HN - PROT SER 2 O 0.28 7.8 9
PROT ILE 5 O - PROT VAL 24 HN 0.08 5.0 4
PROT TRP 6 HN - PROT TRP 6 O 0.10 5.0 5
PROT TRP 6 HN - PROT THR 57 O 1.00 250.0 1
PROT TRP 6 HN Summary: 1.10 127.5 6


My questions are:

1. what's the meaning of time and events? Are the time and event refer to the duration of h-bond (ps) and a numbers of time that we can find this h-bond, respectively? And what's the maximum value of "summary of occupancy" per residue?

2. In the last section of this summary, I've found:

OVERALL
Average number (over ISEL!)= 0.6
<lifetime>= 83.52ps
Total number/frame= 175.3

Is the "lifetime" refer to the total lifetime or what?

Thanks in advance.

<...> means average. The symbol # means "number". So "<time>" means "average [life]time" (in ps) and "# events" means "number [of] events", ie how many times this hydrogen bond was formed. I do not know exactly what the maximum number of the occupancy would be since it depends on how many simultaneous hydrogen bonds that particular atom can form; the summary is not per residue. The final numbers are the averages over all selected donors/acceptors over the whole trajectory. Check corman.doc.

Thank you to Prof. Lennart for your clarification.

For my question below, I'm not sure whether I should post it here or should be in another section.

I just have problem on h-bond analysis when reading more than 1 trajectory. For example, I used charmm 29b2 to simulate the first trajectory and used c32b1 for the other trajectories (because we just have c32b1).

I can not use neither c29b2 nor c32b1 to analyse all trajectories produced by both c29b2 and c32b1.The header produced by c29b2 and c32b1 does not match each other. However, if I use c29b2 to analyse only trajectory produced by c29b2 or vice versa for c32b1, it works well. I'd like to know that are there any possiblities to read trajectories produced by different versions of charmm wihtout any configtion.

Thanks in advance.

Your CHARMM versions probably have been compiled with different settings for little_endian vs big_endian representation of binary data in unformatted files; for c32b1 on gnu/linux systems the default is big_endian with the intel and pgi compilers. It is easy to modify install.com to get the native format /little_endian/ to be used instead.
If your system uses the intel compiler this behavior can be changed at runtime without reinstallation; you can also convert the files you already have. Search the forums, or look at some recent posts on format of dcd files.

Thank you very much Prof. Lennart.
I've checked the format dcd file, as you mentioned.
In my case c29b2 writes little_endian and c32b1 writes big_endian. After I changed dcd file (from c32) in big_endian format to little_endian (or little, c29, to big) but still can not analyse both trajectories together.

PROT 110 VAL HN - PROT 106 ILE O 5.00 250.00
PROT 111 ARG O - PROT 114 ALA HN 5.00 250.00
PROT 114 ALA HN - PROT 111 ARG O 5.00 250.00
***** ERROR ***** HEADERS DO NOT MATCH ô^A^@^@ CORD VELD
*** LEVEL -1 WARNING *** BOMLEV IS -3

I have no idea how to deal with the header. Please help. Thank you.

Dear CHARMMERS,

From the hbond analysis output, I have seen some part of the output which looks like this:
Analysis of hydrogen bonds using cutoff distance= 2.80 and angle= 999.00
Total time= 600.0ps. Resolution= 5.00ps
Occupancy cutoff= 0.00 Lifetime cutoff= 0.0ps

I-atom J-atom <occupancy> <time> # events
---------------------------------------------------------------------------
PROT SER 2 O - PROT GLY 4 HN 0.62 11.9 31
PROT SER 2 O - PROT ILE 5 HN 0.83 29.1 17
PROT SER 2 O Summary: 1.44 20.5 48
PROT PRO 3 O - PROT ILE 5 HN 0.02 5.0 2
PROT GLY 4 HN - PROT SER 2 O 0.62 11.9 31
PROT GLY 4 O - PROT LYS 56 HN 1.00 600.0 1
PROT GLY 4 O - PROT THR 57 HN 1.00 600.0 1
PROT GLY 4 O Summary: 2.00 600.0 2

If I would like to the find % occupation of hbond , for example, % hbond occupation of the O-Ser2 and other residues, how can I obtain it?

I am not sure if I can use "1.44*100" because it's higher than 100%. If I think that 100% means this kind of hbond occurs entire the simulation time, how can I calculate the % hbond occupation? Any advice would be appreciated. Thank you.

The carbonyl of Ser 2 apparently hydrogen bonds to the amides of Gly4 and Ile5, partly at the same time. 1.44*100 is correct if you want to express this total occpancy as a percentage.

Thank you Prof. Lennart. From your reply, it means that % hbond occupation can higher than 100%, is my understanding correct? First I though that the highest value could be 100% but I am not sure.

The maximum fractional occupancy is 1.0 when there is a 1:1 correspondence between donor and acceptor, but a good H-bond acceptor can have 2 donors, and the fractional occupancy can exceed 1.0 in that case.

HI,

I will do this by hand.

Hi,

1) When i was doing analysis for a trajectory for a 20ns simulation with 2000 frames with:
"coor hbond verbose sele monitor end sele .not. monitor end firstu 20 nunits 1 begin 10000 skip 10000"

the following was observed in the output:
Analysis of hydrogen bonds using cutoff distance= 2.40 and angle= 999.00
Total time= 977.8ps. Resolution= 0.49ps
Occupancy cutoff= 0.00 Lifetime cutoff= 0.0ps

>> why only 977.8 ps?

2) the spec "skip" and "nskip" for the traj-spec is same?

regards,
Bryan

Dear Prof. Nilsson,
I would like to analyze the hydrogen bond between lipid and peptide to compare with the literature. The criteria that they use is the hydrogen bond said to exist if donor -acceptor distance is < 2.5 angstrom and the donor-hydrogen..acceptor hydrogen bond form angle > 145 degree. to set that, I would say cuthb 2.5, but to set the angle I could not simply say cutha 145, correct? how could I reach that angle? thank you

Why not? From corman.doc:
"[hydrogen bond exists if ] the angle formed by D-H..A is greater than this CUTAngle (in degrees, 180 is a linear H-bond)".

Note that adding the angular criterion makes almost no difference in practice.

perhaps it is a naive question.

Quote:

---

It is typically not necessary to define a cutoff for the angle.

---

Since there are many difference criteria have been used in the lipid simulation, is this criteria is also reasonable for the lipid system. I am not so sure becuse we are here mainly focused on the peptide/protein systems.
thank you very much

The forces between atoms are computed in CHARMM w/o regard to the molecule type; for H-bonds, the partial charges and VDW radii play a major role, esp. in the close approach distance.

If you are comparing to another simulation result, you should probably try to use the same distance and angle criteria.

Otherwise, I suggest trying the COOR HBOND defaults; we've done that in a publication to evaluate lipid:water and lipid:sugar H-bonds.

Making sure you've properly accounted for image interactions may be a more important issue, esp. if the peptide reaches the edge of the simulation cell during the simulation.

Hi,
I used the code listed here for estimating the number of hydrogen bonds and lifetime. I am wondering is there any way to get the average distance between the donors and acceptors as well ?

Thanks a lot
Diwakar

As of CHARMM version 32 the COOR HBONd command can also give the distance distribution, from which you should easily be able to compute the average. From corman.doc:
"Distance and lifetime histograms can be computed for all (putative) hydrogen
bonds encountered in the analysis"

Thanks a lot for the help.
I have a couple more questions about calculating the number and lifetime of hydrogen bonds.
1) Is it possible to exclude the intra-molecular hbonds?
I calculated the hbonds by selecting each molecule separately and then counting the hbonds of this molecule with other molecule. This process excludes the intra-molecular bonds but I was wondering if there is a better way of doing it.

2) How can I get the number of hbonds as a function of time?

Regards,
Diwakar

1/ You did the right thing, there is no other way to do it. Using the two selections is not all that difficult, is it?
2/ Use a COOR HBOND (without keyword VERBose) in a CHARMM loop and write the number of hydrogen bonds (available as a CHARMM substitution variable, see subst.doc) to a file.

I am trying to get water water hydrogen bonds. I have added the donor list to TIP3 water model as described in the hbonds.inp. I get the same psf without any extra donors on regenerating the psf using this code:

read rtf card name top_all27.rtf
read para card name par_all27.prm

read psf card name @pdbfile@qual.psf
read coor card name @pdbfile@qual.crd

write psf card name @pdbfile@qual-water-donor.psf
write coor card name @pdbfile@qual-water-donor.crd

I want to know if anything more needs to be done to regenrate the psf?

Thanks
Diwakar

Yes, you must use the GENERATE command instead of READ PSF

An alternative is the DONOR and ACCEPTOR commands; they seem to work as long as there only one of each in a given CHARMM input.

DONOr/ACCEptor ADD have some issues, and do not work properly in older versions of CHARMM. It is really very easy to regenerate the PSF so I strongly recommend that route; if your water segment is the last one all it takes is:
read rtf card name ... ! the new RTF w/ DONO, ACCE for TIP3
read param card name ..
read psf card name old.psf
delet atom sele segid wat
read seque tip3 345
gene wat noangle nodihe
write psf card name new.psf
* PSF for ... with DONO/ACCE for TIP3
*

If water residue numbering in the original PSF was not sequential one may change it to sequential with the JOIN command (struct.doc).

A good point; you don't have to start over from the beginning of your model building. What's most important is that the number and ordering of atom types in the PSF does not change. Reading the old PSF and deleting the water segment alone won't work if additional segments were added after the water, such as ions. In that case, you'd need to also delete the segments that follow the water segment in the PSF, and then add them back via GENERATE, in the same order as was done initially.

Dear all,

Sorry if this is a naive question, but I feel a little bit confused about ?NHBOND and ?AVNOHB.

I want to calculate the total number of hydrogen bond between lipid and water and average number of hydrogen bonds per lipid.

Could I just simply say

coor hbond sele resn lip end -
sele resn tip3 end firstu 18 nunit 1 skip 5000
echo ?NHBOND
to get the total number of hydrogen bond?

However, I am not sure if I can apply just ?avnohb to get the correct value for number of hydrogen bond per residue. In addition, is it possible with the charmm command to get the standard deviation for that?

thank you

Yes, ?NHBOND gives the total number The average for ?avnohb is however over selected atoms, not residues (but the conversion should be easy. No, there is no additional statistics available; it is straightforward to output the numbers of interest for relevant timeframes using a CHARMM loop and post-process that with your favorite graphics/spreadsheet/databases/whatever program. Plese not that in this case the VERBose keyword for CORR HBON should not be used.

Sorry for the late response.
I tried to generate the new psf file based on the suggestions by lennart and rmv. The system consists of the solute molecules in a water box. I started from a water box and then I add the solute molecules randomly. After the addition of solute molecules, I remove overlapping water molecules. If I setup the system again using the same procedure, the solute molecules are placed randomly again. Therefore, it is not possible to generate the same system again.
However, I found a quick way to add water donors in the psf.
I edited the donor section of the psf using awk and it works. I need to change the Ndonors and the list of the donor atoms only.

Thanks a lot for all your suggestions.

Diwakar

Hand editing of the PSF should not be necessary. Note that for regenerating the PSF you do not need to perform any coordinate manipulations, all you need to know is how many molecules you have of a certain type, not where they are. Let's say you start with 100 water molecules, then add 30 solute molecules and delete 50 water molecules. To regenerate this PSF you can either
read seque tip3 50
gene wat noangle nodihe
read seque myml 30
gene solu

Sequences could also be read from coordinate files (one for each segment) prepared from the final state of your inital setup.

could you please give me a hint how to get the hydrogen bond autocorrelation function of i.e. N–H---O=C atoms ?

Thank you in advance.

Not directly implemented, but there are a couple of thing you can try.

1. If you are interested in specific donor-acceptor pairs you can use CORRel to extract the timeseries of the distance between the atoms, and the process that with your own program.

2. The MRD option of coor anal (corman.doc) does calculate the intermittent correlation function (ie, no requirement that the "donor"/"acceptor" are within cutoff distance at all times). Something like this (not tested):
COOR ANAL SELE first-atom(s) END SITE SELE second-atom END -
FIRSTU 51 NUNIT 2 NCORS 100 IMRD 21 RRES 3.0

This should output the correlation function in the fourth column for atom(s) in the first selection around the atom in the second selection.

Note that correlation times extracted from these correlation functions are not the same thing (and possiblby more approximate) as the (exact) averages calculated by COOR HBONd. You can also get lifetime distributions from COOR HBONd, which should be related.

1. If you are interested in specific donor-acceptor pairs you can use CORRel to extract the timeseries of the distance between the atoms, and the process that with your own program.

It should be noted that CORREL can also easily calculate an autocorrelation function from a distance time series; there's no need to write a program for this.

Well, in this case it is usually not <r(t)r(t+T)> that is sought (not even <dr(t)dr(t+T)>), but rather something like <h(t)h(t+T)>, where h=1 if r<rcut, and 0 otherwise.

It appears that one of the pieces needed is present in c33 and later, but was never documented. What's missing from correl.doc is a description of the following MANTime command--

STATE realmin realmax

Returns a state function in Q(t), with 1.0 for values between realmin and realmax, and 0.0 for values outside the range.

Hello!

! ACCEPTOR OH2
! DONOR H1 OH2
! DONOR H2 OH2

Regarding with adding the donor and acceptor list for water, I'm wondering whether this also applies to SWM4 drude model. There is only ACCEPTOR OH2 in original SWM4.

RESI SWM4 0.000 ! SWM4-NDP water
GROUP
ATOM OH2 ODW 0.00000 TYPE DOH2 ALPHA -0.97825258 THOLE 1.3
ATOM OM LPDW -1.11466
ATOM H1 HDW 0.55733
ATOM H2 HDW 0.55733
BOND OH2 H1
BOND OH2 H2
BOND OH2 OM
BOND H1 H2 ! for SHAKE
ANGLE H1 OH2 H2
ACCEPTOR OH2

Yes. coor hbond will only recognize hbonds between atoms that have been designated as DONOrs and ACCEptors.