Hi Peter,
I'd like to be able to count Hydrogen bonds between water molecules in a pure water simulation. I've seen the CustomHbondForce class which allows me to model Hydrogen bonding, but what I am interested in is not the modelling itself (I already have that done differently and have to keep it that way for the moment), instead, I am interested in counting the active Hydrogen bonds at every time step, as the simulation evolves.
As the CustomHbondForce allows for custom interactions between triplets of particles, I thought I would use it to do that in the following way:
1) define all water molecules (H,O,H) as donors (d1,d2,d3), define all Oxygens also as acceptors (a1), and define exclusions for bonded particles, so that the Oxygen in the same water molecule does not serve as a donor and acceptor at the same time;
2) define CustomHbondForce expression as "(step(theta - angle(d1,d2,a2)) + step(theta - angle(d3,d2,a2)))*step(maxDistance - distance(d2,a2))", where
- theta is the minimum angle needed between the O1-H1 in water molecule 1 and O1-O2 (between two interacting water molecules) for the hydrogen bond to exist
- maxDistance is the maximum separation distance between two Oxygens for the hydrogen bond to exist.
I would then run the original water simulation and also define a separate system/context/integrator triplet, where the only force would be the CustomHbondForce defined above, and effectively pass the particle positions between the two simulations to get the CustomHbondForce to count the number of Hydrogen bonds per water molecule after the first simulation completed a time step.
It sounds a little long-winded and I was wondering if you had any comments on whether this seems reasonable or suggestions on how this could be simplified. I do not see how I could do it within the same system/context/integrator, as all the force contributions get summed together.
Thank you for your comments.
Spela
Hydrogen Bonding for water
-
Peter Eastman
- Posts: 2583
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Hydrogen Bonding for water
Hi Spela,
That's an interesting idea. It never occurred to me to use a Force to do analysis like that. The more conventional approach would be to just write a routine that takes the positions and identifies the hydrogen bonds. But if you're going to do that every time step, you really want to do that on the GPU or else it will slow things down a lot.
So first: the approach you describe sounds like it would work. The main issue again will be performance, since swapping data back and forth between Contexts at every time step will be expensive.
Because your custom expression consists entirely of step functions, its derivative is zero everywhere. That means it won't actually apply any forces. So one option would be to add it to the same System as everything else. Put it in a separate force group, and to query the number of hydrogen bonds, just call getState(), specifying that you only want that one force group. A problem with this approach is if you are using a MonteCarloBarostat. If so, this wouldn't work correctly. The barostat cares about the energy, not the forces, and your CustomHbondForce will modify the energy.
Another thing to note is that CustomHbondForce is implemented with a simple loop comparing every donor to every acceptor. It doesn't use a neighbor list or anything like that. Usually that's not a problem, because the number of donor/acceptor groups is small compared to the total number of atoms. But since your number of donors and acceptors is close to the total number of atoms, it might be slow. This will of course depend on the size of your system.
Something to consider with all of these approaches is whether you really need to compute this at every time step. What are you going to use the information for?
Peter
That's an interesting idea. It never occurred to me to use a Force to do analysis like that. The more conventional approach would be to just write a routine that takes the positions and identifies the hydrogen bonds. But if you're going to do that every time step, you really want to do that on the GPU or else it will slow things down a lot.
So first: the approach you describe sounds like it would work. The main issue again will be performance, since swapping data back and forth between Contexts at every time step will be expensive.
Because your custom expression consists entirely of step functions, its derivative is zero everywhere. That means it won't actually apply any forces. So one option would be to add it to the same System as everything else. Put it in a separate force group, and to query the number of hydrogen bonds, just call getState(), specifying that you only want that one force group. A problem with this approach is if you are using a MonteCarloBarostat. If so, this wouldn't work correctly. The barostat cares about the energy, not the forces, and your CustomHbondForce will modify the energy.
Another thing to note is that CustomHbondForce is implemented with a simple loop comparing every donor to every acceptor. It doesn't use a neighbor list or anything like that. Usually that's not a problem, because the number of donor/acceptor groups is small compared to the total number of atoms. But since your number of donors and acceptors is close to the total number of atoms, it might be slow. This will of course depend on the size of your system.
Something to consider with all of these approaches is whether you really need to compute this at every time step. What are you going to use the information for?
Peter
-
Spela Ivekovic
- Posts: 26
- Joined: Thu Mar 17, 2011 4:27 am
Re: Hydrogen Bonding for water
Hi Peter,
thank you for your reply.
I'm not using the barostat, so it looks like this could just do the job.
> Another thing to note is that CustomHbondForce is implemented with a simple loop comparing every donor to every acceptor...
Yes, the way I've set it up, the number of donors and the number of acceptors are equal to the number of water molecules in the system. I noticed a mention of an optional cutoff in the implementation of the CustomHbondForce which made me think that using it would be as expensive as doing, say, the nonbonded force calculation on the entire system.
Our system is a water droplet on a substrate which will have pretty much as many water molecules as we can fit into the GPU memory (bar a few additional atoms used to model the substrate). In our case, hydrogen bond counting is very expensive on the CPU. It will also be expensive on the GPU, but I am hoping that, in comparison, it will be faster.
> Something to consider with all of these approaches is whether you really need to compute this at every time step.
I can relax the frequency to every N timesteps, but it will still be quite frequent as I am interested in the transient behaviour during the non-equilibrium stage of the simulation, e.g., droplet spreading on the substrate.
As I am calculating the forces on the GPU, it seems to make a lot of sense to do the hydrogen bond counting on GPU as well. Otherwise I am duplicating the work on the CPU and would prefer to just do the counting inside the force calculation loop on the CPU and eliminate the GPU altogether.
> What are you going to use the information for?
We need it for validation (both of the water model and the code) that our simulated liquid is behaving like water and that the code is giving results that match with benchmarks we have available.
Spela
P.S.: I made a mistake in the first post - the a2 acceptor in the energy expression should have been called a1 to match the problem description.
thank you for your reply.
I'm not using the barostat, so it looks like this could just do the job.
> Another thing to note is that CustomHbondForce is implemented with a simple loop comparing every donor to every acceptor...
Yes, the way I've set it up, the number of donors and the number of acceptors are equal to the number of water molecules in the system. I noticed a mention of an optional cutoff in the implementation of the CustomHbondForce which made me think that using it would be as expensive as doing, say, the nonbonded force calculation on the entire system.
Our system is a water droplet on a substrate which will have pretty much as many water molecules as we can fit into the GPU memory (bar a few additional atoms used to model the substrate). In our case, hydrogen bond counting is very expensive on the CPU. It will also be expensive on the GPU, but I am hoping that, in comparison, it will be faster.
> Something to consider with all of these approaches is whether you really need to compute this at every time step.
I can relax the frequency to every N timesteps, but it will still be quite frequent as I am interested in the transient behaviour during the non-equilibrium stage of the simulation, e.g., droplet spreading on the substrate.
As I am calculating the forces on the GPU, it seems to make a lot of sense to do the hydrogen bond counting on GPU as well. Otherwise I am duplicating the work on the CPU and would prefer to just do the counting inside the force calculation loop on the CPU and eliminate the GPU altogether.
> What are you going to use the information for?
We need it for validation (both of the water model and the code) that our simulated liquid is behaving like water and that the code is giving results that match with benchmarks we have available.
Spela
P.S.: I made a mistake in the first post - the a2 acceptor in the energy expression should have been called a1 to match the problem description.
-
Peter Eastman
- Posts: 2583
- Joined: Thu Aug 09, 2007 1:25 pm
Re: Hydrogen Bonding for water
Even if you need to do this, say, every 10 time steps, that will be a lot faster than doing it every time step. For example, if counting bonds is 5 times slower than taking a single time step, then your simulation will run at about 2/3 the speed it would have otherwise. Whereas if you did it every time step, it would run at 1/6 the speed it would have otherwise!
Peter
Peter