# Hydrogen bonding distance and angle cutoff

Hello professor Steve, Axel and all other lammps users
I have a question regarding calculation of hydrogen bonding in lammps using : pair_style hbond/dreiding/lj command
i have read the documentation carefully for many times. using this command for hbond in lammps, we need two cutoff distances (global inner spline and global cutoff for donor-acceptor interactions) and one global angle cutoff for acceptor-hydrogen-donor interactions.
As all we know, hydrogen bond is usually defined via geometric criteria for example: the A-D distance <3.5 angstrom, and A-H-D angle >120.
I would like to know how we can apply this definition in above lammps command, are these two geometric values the same as global cutoff for A-D interactions and global angle cutoff for A-H-D interactions? if so, what is the value of global inner spline cutoff?

Any help is highly appreciated,

Hello professor Steve, Axel and all other lammps users
I have a question regarding calculation of hydrogen bonding in lammps using
: pair_style hbond/dreiding/lj command
i have read the documentation carefully for many times. using this command
for hbond in lammps, we need two cutoff distances (global inner spline and
global cutoff for donor-acceptor interactions)

(Interesting. I did not notice this earlier. Why do the docs call it
a "spline cutoff"?)

and one global angle cutoff
for acceptor-hydrogen-donor interactions.
As all we know, hydrogen bond is usually defined via geometric criteria for
example: the A-D distance <3.5 angstrom, and A-H-D angle >120.
I would like to know how we can apply this definition in above lammps
command, are these two geometric values the same as global cutoff for A-D
interactions and global angle cutoff for A-H-D interactions? if so, what is
the value of global inner spline cutoff?
Any help is highly appreciated,

Are you asking why we recommend using a 9.0 - 11.0 cutoff distance
(when hydrogen bonds are < 3.5 Angstroms)?

Actually, in the the original Mayo JPC 1990 paper they did not use any
cutoff at all. (This is equivalent to an infinite cutoff.) However
on page 8903 of that paper (footnote, left-hand column), they
recommend using a 9.0 Angstrom cutoff. That's why I use 9.0 Angstroms
for the inner cutoff. They recommended it in their paper.

Journal of Physical Chemistry, Vol. 94, No. 26, 1990, 8897

The outer cutoff (11.0 in the example online) does not matter as long
as it is at least 0.5 Angstroms greater than the inner cutoff. The
inner and outer cutoff distances are usually very similar to each
other.

--- discussion ---

You are correct that hydrogen-bond distances are usually much smaller:
on the order of 2.5-4.0 Angstroms. Beyond this, the interaction
between the donor and acceptor atoms weakens rapidly. But you should
not confuse this with the cutoff distance(s).

Keep in mind that the "Morse(r)" and "LJ(r)" functions never decay to
http://lammps.sandia.gov/doc/pair_hbond_dreiding.html
If we did not use cutoffs, then LAMMPS would have to spend time
evaluating the Morse(r) or LJ(r) functions for every possible pair of
donor and accpetor atoms in the entire system.
As you know, there is no reason to wait for LAMMPS to spending time
evaluating the Morse(r) function when the distance between the donor
and acceptor atoms is extremely large. So the cutoff distance(s) you
outer cutoff distance, the interaction energy between the donor and
acceptor atom is negligible. This does not occur at 3.5 Angstroms (at
least, not if you are using pair_style hbond/dreiding).

You have probably seen this before if you have run other simulations
in the past. For example, the cutoff distance(s) used in the
van-der-Walls (Lennard-Jones) interaction is usually much larger than
"sigma". (Recall that "sigma" is the sum of the atom radii, which is
also near the separation distance where they are most stable).

I hope my explanation was clear.
Cheers!

Andrew

my question is not about the cutoff values for non-bonded LJ and electrostatic interaction
but i would like to calculate the number of hydrogen bonds using pair_style hbond/dreiding/lj command.
this lammps command for determination of hydrogen bonds needs cutoff values(spline cutoff, global acceptor-donor cutoff, global acceptor-hydrogen-donor angle cutoff). on the other, it is known that a hydrogen bond has geometric conditions such as acceptor-donor distance (which is usually less than 3.5) and acceptor-hydrogen-donor angle. i do not know how i can determine the cutoff values in this lammps command, especially the spline cutoff, and if the global acceptor-hydrogen-donor angle cutoff is exactly the same as the acceptor-hyrogen-donor angle used in hydrogen bond definition.
best regards

my question is not about the cutoff values for non-bonded LJ and
electrostatic interaction
but i would like to calculate the number of hydrogen bonds using pair_style
hbond/dreiding/lj command.

I don't think that LAMMPS can count the number of hydrogen bonds
during the simulation (at least not using the criteria you described).
You will need to calculate this using post-processing tools.

To clarify the misunderstanding (I think), the "pair_style" command is
never used to calculate the number of hydrogen bonds. It is used to
control the forces between pairs of atoms during the simulation (which
are not directly bonded together covalently). The "pair_style
hbond/dreiding/lj" command command tells LAMMPS that you want to use
the "hbond/dreiding/lj" force-field to influence the movement of
certain pairs of donor/acceptor atoms.
http://lammps.sandia.gov/doc/pair_style.html
The energy corresponding to this interaction is not directly
proportional to the number of hydrogen bonds. You can't count the
number of hydrogen bonds this way.

Unless somebody writes a "compute hbonds" compute for LAMMPS, it
sounds like you will have to:

1) Use a 3rd-party script to read a LAMMPS DUMP trajectory and
calculate the list of hydrogen bonds (at every dump interval in your
simulation). VMD can read LAMMPS trajectories, and can calculate
hydrogen bonds. Take a look at these web pages:
http://www.ks.uiuc.edu/Research/vmd/plugins/hbonds/
I attached a file containing instructions how to load/view LAMMPS

2) Save your trajectory in DCD format and use ambertools' ptraj to
load the file and count the hydrogen bonds. Take a look at these web
pages:
http://ambermd.org/tutorials/basic/tutorial3/section6.htm
http://archive.ambermd.org/201203/0173.html

3) write your own script to read a LAMMPS dump file and calculate the
hydrogen bonds. (Or, even better yet, if you are familiar with C/C++
then write a "compute hbonds" for LAMMPS. Many features get added to
LAMMPS this way.)

The cutoffs for the Deiding force-fields are not relevant to what it
sounds like you are trying to do.
In a hurry. Must go.

Good luck.
Andrew

this looks like a copy and modify oversight.
further down in the docs it is just called an
inner cutoff.

axel.

I don't think that LAMMPS can count the number of hydrogen bonds
during the simulation (at least not using the criteria you described).

it can. but only those that are added by the hbond/dreiding pair styles
and that follow the dreiding specs and assumptions. for starters, i don't
think that the dreiding definition of the hbond angle is a good one. they
use the A-H-D angle, but using the A-D-H angle is considered the better
choice for a number of reasons.

from the docs:

These pair styles tally a count of how many hydrogen bonding
interactions they calculate each timestep and the hbond energy. These
quantities can be accessed via the compute pair command as a vector of
values of length 2.

To print these quantities to the log file (with a descriptive column
heading) the following commands could be included in an input script:

compute hb all pair hbond/dreiding/lj
variable n_hbond equal c_hb #number hbonds
variable E_hbond equal c_hb #hbond energy
thermo_style custom step temp epair v_E_hbond

You will need to calculate this using post-processing tools.

that being said, using post-process is very *much* recommended
for this kind of analysis, exactly because you can use a simpler
simulation (no need to compute interactions that you cannot add).
however, this analysis could probably be done with LAMMPS
using the rerun feature and the compute pair thing.

again, this exchange proves how important to formulate
questions well and make an effort to explain the big picture

axel.

as i said in my previous e-mails, i have read carefully the doc pages for this command, and i was able to compute the number of hbonds in my system,
however, i am not still sure that whether the values i used for inner spline cutoff and global cutoff for A-D, and global angle cutoff for A-H-D are correct. for inner spline cutoff, i used the default value of lammps doc page, for global distance cutoff i used 3.5 Angstrom and 120 for global angle cutoff (the distance value of 3.5 and angle value of 120 are the most widely used ones for hbond deffinitions in literature).

best regards

as i said in my previous e-mails, i have read carefully the doc pages for

huh? when? where? do you have a split personality?

this command, and i was able to compute the number of hbonds in my system,

you obviously didn't read it "carefully".
you may have computed some number,
but i doubt that it was a correct number.

however, i am not still sure that whether the values i used for inner spline
cutoff and global cutoff for A-D, and global angle cutoff for A-H-D are
correct. for inner spline cutoff, i used the default value of lammps doc
page, for global distance cutoff i used 3.5 Angstrom and 120 for global
angle cutoff (the distance value of 3.5 and angle value of 120 are the most
widely used ones for hbond deffinitions in literature).

as was explained several times already. pair_style hbond/dreiding/*
is a method to *add* a hydrogen bonding energy term (and the resulting
forces!) to a system according to the definitions of the dreiding force field.
whether this complies with what you want to do is of no apparent
concern to anybody but you.