Chemical Reactions plugin update - discussion

[Piotr Janik]

Hi,

I don't want to add noise to existing PT story, so I'm starting a thread here.

I've been looking at latest changes, discussions and the general idea looks very good to me. Previously, when the bond was created, we were adding a bit less KE to atoms to compensate change in LJ potential (which is not calculated between bonded atoms). So the bond formation was causing “random” increase in KE, depending on how close atoms were to each other. As chemical reactions plugin is ran quite rarely, it could be really random. Now the KE increase should be constant, so it’s perhaps much easier to setup model according to our needs.

It will definitely affect our models, e.g. this test doesn't work as it used to work before:

http://lab.concord.org/branch/chemical-reactions-energy-conservation/interactives.html#interactives/conversion-tests/chemical-reactions/simplest-free-radical-reaction.json

...but it makes sense - now we add a bit more KE, so the bond breaks. I can tweak the model so it works like before. If Dan reviewed other models and they worked good enough, that seems totally fine.

I would say that perhaps the cleanest solution would be to just calculate LJ all the time, even between bonded atoms. But that's bigger change, I don't think it's worth it, and I'm not 100% sure if it's totally correct as well.

When it comes to divergence issue (e.g. in water reaction model): 

1. I think I traced the problem and actually it's caused by the new addKEToAtomPairAndConserveMomentum function. The quadratic equation that we try to solve there has no solution sometimes (disc < 0). I get the general idea, but I don't understand all the steps we perform there. We can either fix this function or use old one that doesn't conserve momentum.

What's interesting, this model diverges when timestep is equal to 0.15, but works totally fine when it's 0.17, 0.20 or 0.30. ;) So it doesn't seem to be caused by too big accelerations etc.

2. In a few places we have:

  if (creationPotential > 0) { return; }

Out of curiosity - did it fix some real issues or it's just supposed to prevent possible issues?

3. I remember the idea that a bond can break while atoms are very close to each other, then huge LJ interaction is applied and model can easily diverge (perhaps Sam mentioned that on Slack).

Similar question - did someone observe that or was it speculation?

If this problem really exists in practice (it does sound likely, but perhaps mostly depends on parameters we use - bond strength, LJ interaction etc.), we could ensure that a bond can break only when it's stretched, not when it's squeezed. That should ensure that enormous LJ interaction can't happen.

Regards,

Piotr

[Piotr Janik]

If this problem really exists in practice (it does sound likely, but perhaps mostly depends on parameters we use - bond strength, LJ interaction etc.), we could ensure that a bond can break only when it's stretched, not when it's squeezed. That should ensure that enormous LJ interaction can't happen.

Correction to my previous email - actually we already do it for breaking bonds (current bond length has to be > bond length). However this still could happen during bond transfer.