CustomHBondForce with repeated atom in different donor groups.

[Carlos Bueno]

I'm working on a forcefield that uses the CustomHBondForce for multiple body interactions. I obtain unexpected results when I include the same atom in multiple donor groups. Some particles that are outside the CustomHBondForce cutoff experience high forces.

Is there a reason why I shouldn't include the same atom in multiple donor groups?

It seems that the problem can be solved if I create multiple instances of the CustomHBond force, repeating the acceptor groups, but separating the donors. That way I don't repeat the atom index in the donor groups on the same instance.

[Peter Eastman]

I don't know of any reason that would cause problems. Which platform are you using? Try telling it to use the Reference platform. Does that also give the same result?

Can you create an example that demonstrates the problem? If I can reproduce it, then it's much easier to figure out what's going on.

[Carlos Bueno]

I was using the OpenCL platform. I tried to use the CPU platform and the Reference platform but I could not reproduce the unexpected results.

I attach an example of the problem. Some CAM particles outside the non-bonded cutoff seem to experience a high force.

[Peter Eastman]

What exactly should I be looking for in your example? It just prints out the energy and then runs a simulation. The initial energy is almost identical for OpenCL and CPU platforms. I also tried computing the initial force on every atom with both platforms and comparing them. No atom has an absolute difference in force between the two platforms of more 0.1 kJ/mole/nm.

[Carlos Bueno]

The temperature should raise faster on OpenCl than on the CPU platform.

Temperature_increase.png (29.66 KiB) Viewed 364 times

For example on the CPU platform at 3 000 000 steps:
#"Step","Time (ps)","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)"
3000000,300000.0000019568,-355.8944037050504,1451.904116005339,1096.0097123002884,454.7492718491549
while on the OpenCL platform at 3 000 000 steps:
#"Step","Time (ps)","Potential Energy (kJ/mole)","Kinetic Energy (kJ/mole)","Total Energy (kJ/mole)","Temperature (K)"
3000000,300000.0000019568,-460.6856689453125,4271.379672840238,3810.694003894925,1337.834072238581

I think the error is more perceptible on the simulation trajectory and not at the beginning.
Some particles will accelerate with a preferred direction, which is what is causing the raise of temperature.

[Peter Eastman]

I was hoping for something more direct. You said the CustomHBondForce was applying forces between atom pairs that should be beyond the cutoff. What exactly led you to that conclusion?

[Carlos Bueno]

When I traced back the particles that move fast during the simulation, they accelerate when they are still far from the Hbond aceptor. They don't seem to accelerate as a result of a collision, and they don't seem to be attracted to any other particle.
The cutoff for the CustomHBondFoce is 12 nm, while the distance to the ACT when the particles accelerate may be much bigger than that (from 36 nm to more than 100 nm).
When I turn off the CustomHbond or when I don't repeat the donor index I don't have that problem.
My conclusion was that the repeated index on the CustomHBond donor creates the problem.

On the example I sent, each Cc atom from the CAM appears 6 times as part of a Hbond donor.

[Peter Eastman]

I've had it running for a while now. It's over 1 million steps, and the temperature is still just fluctuating between about 350 and 500 K.

Can you create a test that more clearly shows incorrect behavior? For example, if you can identify a specific time step where you think the forces are incorrect, call simulation.saveState() to write out the current state to a file and then post it.

[Carlos Bueno]

Yes, I just noticed that it may be easier to see the error with the sum of force on the system.

Mean_Force.png (21.76 KiB) Viewed 328 times

On this simulation the error becomes noticeable around frames 75-76.
The particle with indices 1235, 1236, 1237 starts drifting away.
I attached the states.

[Peter Eastman]

Thanks! I can confirm that the sum of all forces is nonzero for state 76. But that's actually coming from one of the CustomNonbondedForces, not from the CustomHBondForce. I'll investigate why.