-1
$\begingroup$

i am trying to implement/extend an implementation of Lennard-Jones potential simulation regarding Xenon molecules (for curious ones, the code can be found here functions (force_naive->lj_force)).

First and foremost i altered the above code to make collisions full elastic with the box bounds.

For box bounds:100 Ang, particles: 16348 and dt: 1e-15 sec i am trying to calculate the total energy of the system but the numbers doesnt make sense to me for example at beginning total energy is Xe+12 (Joules?!?!) and it continues rises until float overflow to inf.

For calculation of total energy i am using the typo found here page 5.

I think the problem is at velocity calculation. Someone any help? Before misinterpreting my intentions my question is about the continuous rising of total energy and if it makes sense? cause in my understanding its definitely not. I thought that the total energy should be pretty low and in theory at equilibrium state even zero. Am I right?

Improve this question
$\endgroup$
8
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because it is about analyzing software and not physics. Perhaps Computational Science might be better suited. $\endgroup$
    – Kyle Kanos
    Jul 21, 2015 at 16:41
  • $\begingroup$ Its analyzing physics quantities in software implementation! I am not asking you how to implement this in software i am telling that its doesnt make sense to me so much Joules in the system and the effect that it keeps rising until its inf! Thats where i want help if my sense is right. So please try to be more helpful and less judgefull $\endgroup$
    – ribben
    Jul 21, 2015 at 16:51
  • $\begingroup$ You are using a very high density of about $10^{28}~m^{-3}$, corresponding to an interparticle distance of 4 angstroms. Therefore there is a high probability of each particle being deep inside the repulsive core of another particle. Why don't you try decreasing the density? If your collisions are elastic then the energy should not rise, but I expect this is due to a numerical instability of your code. $\endgroup$
    – Mark Mitchison
    Jul 21, 2015 at 17:03
  • 1
    $\begingroup$ @ribben: Yes, analysing software is not what we do here. From the most-recent computational policy: ...we are not a programming site. If your question is about implementing computational code - in particular, if it's about writing, compiling, debugging or optimising code, or about a specific language or library - then it is off topic. $\endgroup$
    – Kyle Kanos
    Jul 21, 2015 at 17:28
  • $\begingroup$ @KyleKanos: why you refuse to understand i am asking if its normal to see an continuous rising in total energy of the system from physics aspect. i posted the implementation not to point me the error just in case anyone is curious enough to see the exact procedure please stop spamming. PS: im not asking how to write it, how to compile it, how to debug it or optimize it. I am trying to understand what the physics means and where is the error in my physical quantities $\endgroup$
    – ribben
    Jul 21, 2015 at 19:24

1 Answer 1

Reset to default
0
$\begingroup$

The effect you describe is unphysical. The energy for a randomized, undriven system should never rise when the system moves into equilibrium. If your energy rises (and especially, if it diverges!) there is a problem in your code. This is most likely due to numerical instability in your method for integrating the equations of motion. This is easy to check by reducing your stepsize. You should see the rate of divergence (measured in system time, not number of steps) drop as you use smaller steps.

Otherwise, it is most likely a bug in your code.

Improve this answer
$\endgroup$
1
  • $\begingroup$ You are right when i reduce the stepsize indeed the rate of divergence drops. So i guess its something numerically which goes terribly wrong. thx $\endgroup$
    – ribben
    Jul 21, 2015 at 20:25

Not the answer you're looking for? Browse other questions tagged or ask your own question.