Thread: Please help for a numerical calculation

Hello everybody! A pleasure to be on this board.

My problem is as follows: As part of a molecular dynamics simulation on a simple system, I already have a basic open source MD-program. Thing is though, I need to edit it to be able to do one more thing. The experiment itself is carried out in a cubic box, where atoms float around happily. For the experiment, I need the density function in one single coordinate, which is akin to placing two sheets of paper inside the box, and counting how many atoms are in between.

Now, I am aware that the problem will probably go into the differential region, and that I will probably have to integrate something (that much mathematics even a chemist like me knows ), but I lack the knowledge of how I could do this, or how I could implement it.

My idea so far was to simply let the computer literally "count" how many of the atoms have an x-coordinate between two numbers (say: 3 and 3.1), but this is a computational disaster and absolutely inadequate to do statistics.

(For those who are interested: The experiment I am trying to reproduce is "Surface structure of a liquid film", M. Rao and D. Levesque, 1976. The paper is very nice and fluid to read, but they unfortunately don't mention the formula by which they calculated the surface density profile, or I have overlooked it.)

I hope somebody can help me out.

Thanks very much in advance.

Is your density uniform? You say that the atoms are "floating around happily". That would seem to indicate a uniform density.

Hi there!

Thank you for your quick response! Yes I indeed forgot to mention this.

In this case, the distribution should be uniform on average at first. The main trouble is that the goal of this project is to simulate the liquid film. Some atoms are expected to leave the liquid state and become gasious. So what you could expect if you had a density profile would be in fact sort of a sigma-curve, with the density jumping up radically, since you have more atoms in a liquid than in a gas.

If you would like me to explain fully, just let me know, but I did not want to go into too much detail in the interest of not wasting anyone's time.

Thanks again

Well, the more data we have, the better the model. Can you draw a diagram of the box, with gravity noted, and the level of the liquid noted, and the direction in which you want the density function noted, please? Also, is the liquid uniform density? And is the gas uniform density? Thanks!

Wow, ok... Then here I go.

The way the simulation is done, is this: You take a cubic box, all sides equal length. You then will up the box full of argon atoms (around 5000), and you let those move a bit (this is called equilibration). Now the density is not perfectly uniform of course. A few areas have more atoms than others, if you take a snapshot. On a mean, the atoms should be uniformly distributed however.

What you then do is to double the box along one of the axes in both directions. So now you have 5000 argon atoms sitting in a box, where you have more space. Now, some of the atoms will stick in the "liquid", by internal forces. A few of them however will break out and go into the empty region.

The density profile thus is not uniform once you take the snapshot. You first have a low density (a few atoms in the otherwise empty space. Those are the gaseous ones), then, a sharp jump, where you have a high density in the tightly packed pseudo-cube, and then again a sharp jump down again where you have gas again, if you can picture this.

Gravity in this case plays no role, since atoms are so light that it's hard to even show evidence for gravity on this level. The atoms are just controlled via electromagnetic forces, which you can approximate by a hard ball model, which has the atoms shooting around space like billiard balls, and on collision, the atoms move in new directions again.


I am rather convinced that this is calculated via an integral (counting, as I noted, the number of atoms in a small area), but I am not certain.

I hope this helps. If you really need a drawing, I will try to make one (though I am horrible at drawing...)

Thanks again.

Thank you for the clarifications. I still have more questions, though.

What you then do is to double the box along one of the axes in both directions.

So you double the lengths of only two perpendicular sides? Or all three? Do you end up with a cube again? I guess what I'm asking is this: does the volume quadruple, or octuple (not sure this is a word, but I'll coin it for now)? And you do this without increasing the total number of argon atoms, correct? So now the same number of argon atoms are free to move around in a much larger box. Correct?

Now, some of the atoms will stick in the "liquid", by internal forces.

Not sure I understand this. Do you mean they condense? So you have mixed phase?

The density profile thus is not uniform once you take the snapshot.

Not uniform in time? Or not uniform in space? If the latter, what is the function that describes the density as a function of x, y, and z? It sounds more like you're talking about a non-uniform distribution in time, which I can certainly imagine to be the case.

Gravity in this case plays no role,...

Even on the liquid?

Ok. I think there is something we can say, even now. Let's say we have a density function that depends on the three cartesian coordinates We'll denote this density function and it has units of say. We wish to have instead a single-dimensional density function, call it that does not depend on It needs to have units of Is that all correct so far? In that case, I believe you would define your density function as follows:



Here, is the area of the cross-section perpendicular to the direction .

Make sense? This expression has the correct units, because you're multiplying the integrand's units by two spatial dimensions (the differentials).

Now, if the area is constant along the entire cube/rectangular prism, and the density is uniform, then you get the following simplification:



Make sense?

Finally, if you wanted to know how many atoms are between two sheets (in this special case), positioned at and , then you would perform the following computation:



Otherwise, if you're not in the special case, you would compute



Any questions about any of this?