search

openmm / openmm

OpenMM is a toolkit for molecular simulation using high performance GPU code.
1.42k stars 505 forks source link

Massless atoms destroying/altering motion of system #4527

Open chiara-corsini7 opened 2 months ago

chiara-corsini7 commented 2 months ago

I am encountering problems in the dynamics of my system due to massless atoms. As found before by #2113 the system energy explodes when two atoms overlap even if one is massless and has no LJ site, due to a bug in the NonbondedForce and in my case also in a CustomNonbondedForce. This was also reported here https://simtk.org/tracker/index.php?func=detail&aid=1591&group_id=161&atid=435 in 2012. And if does not explode the system dynamic of my system is unrealistic as all molecules try to "avoid" the massless sites. One solution would be to add Exclusions or Exception lists but due to size and nature of my system it would be too onerous as explained by #2571.

Could this bug be fixed without having to resort to exclusion lists?

peastman commented 2 months ago

This isn't something that can be "fixed" in any simple way. We could add a check to the nonbonded kernel for atoms at exactly the same position, but that would add overhead to every interaction. It would make every single simulation run in OpenMM slower to deal with an obscure corner case, which obviously isn't acceptable.

Can you describe more about your simulation? There may be a reasonable workaround in your case. What are you using massless particles to represent? Why do the particles have no LJ term? Why are there particles at exactly the same position?

chiara-corsini7 commented 2 months ago

I am using the massless particles to restrain the x and y coordinates of a few atoms of two solid slabs which I am approaching along z to perform PMF simulations through umbrella sampling.

As I am placing the particles on top of the atoms to restrain the system explodes. By displacing them a bit it I can still run the simulation but as the rest of the solid moves along z it deforms to avoid the massless particles. I can shift in z the massless particles to not interfere with the solid slabs when approaching, but still I have to create an exclusion list between the massless particles and the atoms of the solvents surrounding my solid slabs. Mind that all these simulations are polarizable via drude particles.

Without solvent I have around 17000 particles in my simulation which can easily increase to ~150000 (~130000 particles describing the solvent) . So the exclusion list would be enormous even if I try to restrain just 30 atoms in each slab of the solid.

Why do the particles have no LJ term? I did not express myself well. They do have LJ terms as they cannot not have one, but it is set to sigma 1 and epsilon 0.

peastman commented 2 months ago

I think the problem in your simulation is different from what you think. The issue you linked only applies when two particles are at exactly the same position, down to the last bit in all three coordinates. If they are just very close, everything works as expected and the forces are zero.

You said the other particles are deforming to avoid the massless ones. That indicates there is a force between them. What is that force? It sounds like you have some kind of interaction in your simulation that you didn't intend.

It also isn't clear to me if you really need the massless particles at all. If your goal is to apply a restraint to a few atoms, could you use a CustomExternalForce with per-particle parameters to specify the location they're restrained to? See this cookbook entry for how to do that.