[lammps-users] How to calculate the number of H bonds in polymer system?


As mentioned in the title, are there any methods to count the number of interactions between certain atoms ? For example, I have a polymer system, and I want to count the number of H bonds formed between H atoms and H atoms in this system. Do you know any method to achieve this goal ? Thank you very much.

Shuting Wang

LAMMPS doesn’t have a particular compute style suitable for this kind of calculation.
There are many (slightly) different definitions for this kind of property, so it would require some effort to have a generic and sufficiently flexible implementation. One challenge is to figure out what actually constitutes a hydrogen bond, since LAMMPS internally only “knows” atom types by their type number, and thus has no idea which atom is a hydrogen atom, which atom would be a hydrogen bond donor and which atom would be a hydrogen bond acceptor.

There is a specific hydrogen bond add-on pair style which is required for adding hydrogen bond interactions explicitly as part of the dreiding force field. If I remember correctly, some people have managed to apply these to do some hbond analysis as post-processing in LAMMPS with the “rerun” command. However, I cannot say how meaningful and applicable this is for non-Dreiding calculations.

There is a hbond tool in VMD, where the above classification is based on CHARMM force field atom type assignments (which are strings), thus you would have to construct a suitable .psf file that has the proper atom type (and atom name, and residue name) assignment(s) so that the classification works as expected. It can be done, but is non-trivial and requires sufficient knowledge of VMD Tcl scripting.

There are likely other tools in other MD analysis/post-processing packages (some are bundled with other MD codes like Gromacs, some are standalone packages). You would have to check out if you find anything suitable in the list at: https://lammps.org/prepost.html or use a search engine to look elsewhere.

Failing that, you may consider writing your own tool (possibly based on an existing implementation) to post-process your data. that gives you the flexibility to implement exactly the h-bond definition that you prefer.


You should have a look at MDAnalysis.


Simon Gravelle

Institute for Computational Physics
University of Stuttgart