Hydrogen bonding distance and angle cutoff

Thanks Axel.

Dear Rezaei Mohammadhasan,

   I'm sorry that I did not read to the end of the doc page, so I was
wrong. If this is what you were doing, I'm sorry. If it helps,
perhaps I can explain Axel's "rerun" suggestion to you with more
detail:

   Axel said that you can use "pair_style hbond/dreiding/lj" and
"compute hbond/dreiding/lj" for analysis. When you do this, you can
use the cutoffs you desire (for example a 3.5 angstrom distance
cutoff, and an A-H-D angle > 120 degrees). See example below.

  But you must do this after you run your simulation.
(postprocessing) You would not want to use these cutoffs during your
simulation because they would generate unrealistic hydrogen-bonding
forces. (The Dreiding authors recommend using a 9.0 Angstrom cutoff
during the simulation, for example, not 3.5 Angstroms.) You don't
have to use "hbond/dreiding/lj" during your simulation at all.
Perhaps you realize this. (If so, my apologies, again.)

  To use "compute hbond/dreiding/lj" in this way, you can use the
LAMMPS "rerun" command. The "rerun" command allows you to
post-process a DUMP file trajectory that you have created during a
previous simulation:
http://lammps.sandia.gov/doc/rerun.html

To use these cutoffs, try putting this at the beginning of your LAMMPS
input script file:

    pair_style hbond/dreiding/lj 4 3.4999 3.5 120

and then, as Axel suggested, add these commands somewhere in your
input script before the "rerun" command:

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

Cheers
Andrew

Thanks Andrew for your reply,
I have successfully used the commands pair_style hbond/dreiding/lj and compute hbond/dreiding/lj and computed the number of hydrogen bonds in my system containing an aromatic polymer and some water molecules.
in my input file, i used the command pair_style hbond/dreiding/lj 4 2.5 3.5 120, that is, i set the global inner spline cuttoff ditstance to 2.5 angstrom, the global outer cutoff distance to 3.5 angstrom and the global angle cutoff to 120 degree. it should be noted that the values of 3.5 and 120 are those used in geometric hbond deffinition. therefore, according to this setting and hydrogen bond deffinition, whenever the distance A-D is less than 3.5 and the angle A-H-D is more than 120, a hydrogen bond between acceptor and donor should be formed (is this true?). however, i do not know how the value of global inner spline cutoff is determined and also how it affects the number of hbonds, because in original DREIDING paper this paramer has not been discussed.
for calculating the nonbonded interactions (lj and coulomb), a cutoff value of 12 angstrom was used and 3.5 angstrom is only used for hbond interaction (using rerun command).
thanks again for your detailed descriptions
best regards

Thanks Andrew for your reply,
I have successfully used the commands pair_style hbond/dreiding/lj and
compute hbond/dreiding/lj and computed the number of hydrogen bonds in my
system containing an aromatic polymer and some water molecules.
in my input file, i used the command pair_style hbond/dreiding/lj 4 2.5 3.5
120, that is, i set the global inner spline cuttoff ditstance to 2.5
angstrom, the global outer cutoff distance to 3.5 angstrom and the global
angle cutoff to 120 degree. it should be noted that the values of 3.5 and
120 are those used in geometric hbond deffinition. therefore, according to
this setting and hydrogen bond deffinition, whenever the distance A-D is
less than 3.5 and the angle A-H-D is more than 120, a hydrogen bond between
acceptor and donor should be formed (is this true?).

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

It looks like these commands just count the number of hydrogen bonding
interactions computed by LAMMPS.

A hydrogen-bonding interaction is calculated by LAMMPS for every
triplet of atoms i,j,k which satisfy the atom type requirements, and
distance and angle cutoffs you have specified (3.5 and 120). These
are the correct cutoffs, then yes, this is the number of hydrogen
bonds as you have defined them.

After looking at the source code, I suspect this number of hydrogen
bond interactions is stored in the "hbcount" variable. (See
"pair_hbond_dreiding_lj.cpp", line 230) That variable counts the
total number of times inner-loop of PairHbondDreidingLJ::compute()
function is visited at a particular moment in time. (It is later
stored in "pvector[0]")

Hopefully, this number is calculated correctly when you use these
commands above with the "rerun" command. You will have to test it and
see. (I am not sure if LAMMPS invokes any of the Pair::compute()
functions if you are using the "rerun" command. If not it will be
obvious, because the hydrogen bond count will always be zero.)
"Rerun" is a new feature and it may not behave the way we expect it to
for every pair style.

You are using features of LAMMPS (rerun and compute hb all pair
hbond/dreiding/lj) which may not have been tested when used together.
Bugs are possible.

Hopefully this works for you.

however, i do not know
how the value of global inner spline cutoff is determined and also how it
affects the number of hbonds, because in original DREIDING paper this
paramer has not been discussed.

The inner cutoff is is only used to make the energy function more
smooth near the outer cutoff. And this only matters when you are
using pair hbond/dreiding/lj to add forces to your atoms during an
ordinary simulation. It does not effect the number of hydrogen bonds
calculated during post-processing (using the "rerun" command).

This parameter does not effect the number of times the inner-loop of
the dreiding force calculator is visited.
It will definitely not effect your hydrogen bond count. I just set it
to 3.499, but you could set it to 3.141592653589793. It does not
matter.

Cheers!

I will be mostly away from email for the next week.
Good luck

Andrew

Many thanks Andrew,
best regards