Calculating the number of hydrogen bonds in the system using hbond/dreiding pair_style

LAMMPS LAMMPS General Discussion
,
1

Dear LAMMPS users,

I want to calculate the number of hydrogen bonds in my system. As I understand, there is no direct way to do this via any computes (is it correct?). From an old discussion Hydrogen bonding distance and angle cutoff I conclude that I can try to use hbond/dreiding/lj pair_style to define the Hbonds in the system (by defining the A-H-D cutoff angle and A-D cutoff distance) during postprocessing. As I understand, I can use

pair_style hybrid/overlay [potential that I used during the simulation] hbond/dreiding/lj 4 2 3.5 120 (assuming 3.5 and 120 are the distance and the angle in the Hbond definition)

and get the number of hydrogen bonds in the system by

compute hb all pair hbond/dreiding/lj
variable n_hbond equal c_hb[1] #number hbonds
thermo_style custom step temp v_n_hbond

Is it okay to use this pair_style with other potentials, say lj/cut/tip4p/long, with hybrid/overlay assuming that I do not run any kind of simulation but just do this Hbond counting for a single frame as postprocessing? Does it only compute the number of interactions that satisfy the distance/angle criteria?

Thank you in advance.

2

The pair style does what its documentation says it does. If in doubt, set up some test cases (e.g. with just a couple water molecules in different relative positions and orientations) and check the results against manual computation.

Please note that this will give you the hydrogen bond counts according to the specific definition of hydrogen bonds of the Dreiding force field interactions. Those may or may not coincide with the specific definition of how hydrogen bond should be counted that you prefer (there are several in use).

There are plenty of alternatives, though. VMD, for example, has hydrogen bond computation tools (you only need to manufacture a proper .psf file or find some other way to convert the numeric LAMMPS atom types into atom types/names that VMD will recognize as hydrogen bond donors and acceptors. VMD usually applies conventions from the CHARMM force field. I know the MDAnalysis package also has tools specifically for hydrogen bond analysis and there are likely others, too.

3

Thanks a lot for your reply.