CHARMM Development Project Forums CHARMM Community Script Archive hbonds.inp

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