Wrong force calculation

[pkoutsovasilis]

Firstly i think you are calculating wrong the potential itself

float lj = 24 * epsilon (2 * pow(sigr, 12) / (dist * 1e-10) - pow(sigr, 6) / (dist \ 1e-10)); // Newton.

i cant find out why 24 * epsilon or 2 * pow(sigr,12)

Secondly and most important the final force for this kind of potential can be found here is this R(r) = 4ε [12_σ^12/r^13 - 6_σ^6/r^7]

Finally its not clear if your bbox is implemented in angstrom units and in general its not clear in what units you are working or at least the version on github is not producing the same results as the one in video (atoms stay on bboxs corners if the computation of force is not altered)...

[acappiello]

I think you're misreading the code. The line you're questioning is actually the force, not the potential. It's just rearranged from the form you've written to pull out constants and rearrange terms to decrease the number of pow calls. The potential itself is not needed at all for the computation.

You're right, units are not as well documented as the should be. If I remember correctly, other than where I use angstroms for distance measurements, everything else should be in standard SI units.

When I run the program with the default values, I see a small number of particles stuck in the corners, but otherwise behaves reasonably. Since I was interested it running with a wide variety of parameters using the command line arguments, I wasn't too concerned with finding perfect "default" parameters. You'll notice if you pause the video around the 7 second mark, it is not being run with default parameters there either. If you're seeing something vastly different than that, I'm not sure what to tell you except maybe you changed something else.

[pkoutsovasilis]

The code you implement for Force calculation is an approximate version of Lennard-Jones for computational purposes (you should have mentioned it). Also the md.cl file doesnt has a definition for constant SIZE (opencl no compilation). Furthermore if your units are already in angstroms why are you translating them in meter in the force calculation, as you can see here (http://chemwiki.ucdavis.edu/Physical_Chemistry/Physical_Properties_of_Matter/Intermolecular_Forces/Lennard-Jones_Potential) in the solutions section, example 2, (σ/r) must be in the same measure unit and in your calculation σ is angstroms for Xenon but distance is in meters... Last but not least you make the assumption that the velocity of particles is divided in half by elastic collision with box bounds but then its not an elastic collision and the wall must have accumulatively the half kinetic energy of each particle

[acappiello]

SIZE is defined at compile time, since I wanted it to be easily adjusted without hardcoding. You'll see this done on line 96-97 of md.cpp. Obviously, you're free to define it elsewhere to suit your needs.

When I precompute sigr = sigma / dist, both values are in angstroms. sigr is now a unitless quantity, so there is no discrepancy. In effect, both are converted to meters, but the conversion factor cancels out. I acknowledge that going back and forth between angstroms and meters is clumsy and ought to have just stayed in meters.

I honestly don't remember why I made the collisions inelastic, but changing that value to 1 will address that. In fairness, I can't find anywhere that I claim the collisions to be elastic.