I would like to calculate the number of hydrogen bonds between water molecules and a membrane in my simulations. In several papers, hydrogen bonds are defined using the following criteria:
The distance between the donor and acceptor atoms is less than 3.5 Å.
The hydrogen–donor–acceptor angle is smaller than 30°.
I would appreciate guidance on the following points:
What is the simplest and most reliable way to compute hydrogen bonds using these criteria?
Is there a built-in method in LAMMPS to perform this calculation?
If not, what output data (e.g., trajectories, atom IDs, bond information) should be saved from LAMMPS for post-processing?
Are there recommended post-processing tools or scripts for this purpose?
Any comments or suggestions would be greatly appreciated.
From your description I think you could make use of fix bond/create/angle that will essentially “create” the hydrogen bonds. You can use a bond style that has no interactions so that it will not affect the simulation, and then count the new bonds using the relevant computes provided by lammps
That is not as straightforward as you might think, since adding new bonded interactions to the bond topology will change the list of excluded neighbors and could thus change the force field quite significantly in unexpected ways and would only work reasonably for a force field with implicit bonds. See Pulling on a carbon nanotube - for a discussion.
Thank you for the clarification — that makes sense.
My goal is only to analyze hydrogen bonds between water and the membrane (with LJ potential) using geometric criteria (distance < 3.5 Å and angle < 30°), not to introduce real bonded interactions into the force field.
Given that, would you recommend doing this entirely in post-processing rather than using fix bond/create? If so, what would you consider the most reliable workflow with LAMMPS — for example, dumping trajectories and using tools like MDAnalysis, VMD, or a custom script to apply the geometric criteria?
I’m mainly looking for a robust and standard approach that does not perturb the simulation physics. Any guidance on best practice would be greatly appreciated.
Since you asked this is not something I know how to do easily in LAMMPS. However there is an MDA module and a VMD module from which the former is inspired to do so as a post-process.
However note that MDA LAMMPS parsers are very limited at the moment but a data file can be used as a topology file. So this should work. The documentation allows you to tune your distance and angle criterion but you would have to work a bit if you want something more specific than the initial donor-selector guess from mass and charge.[1]
By default hydrogens are guessed as the atoms with a mass between 1.1 a.u. and 0.9 a.u. and a minimum charge of 0.3 ↩︎
I used to say to some colleagues of mine “the only important thing that LAMMPS can’t do yet is coffee[1], but you can simulate hydrated caffeine molecules.”
Now that I am about to ditch out OVITO with your new features, I’m waiting for the day I can say that “LAMMPS coffee may not be the best, but it makes coffee.”
Which is partially incorrect as there are now many scripts and building plans for coffee machines connected to a network through ssh with a simple computer such as a raspberry pi. So you could execute the scripts in a LAMMPS session through the shell. ↩︎
This is amusing, since I myself have to avoid caffeine like the plague.
That said, I could see that - given a suitable API - you could operate a french press coffee maker with fix deform for example.
I better stop here, we’re getting off-topic very fast.
P.S.: don’t throw out OVITO entirely. There are many things it can do well, that LAMMPS will not be able to.
P.P.S: I feel like we just wasted a great opportunity for an April’s fools joke.
): 