density profile

Hello all,

Maybe someone here will have some insight into my problem. I am calculating the density of a confined liquid using the following command in LAMMPS :

fix 1 liquid ave/spatial 10 20000 200000 z lower ${DZ} density/mass norm all units box file production/liquid.profile

I get a very good density curve and am satisfied with the way it looks.

I then calculate the integral of the density curve with Matlab to check that it corresponds to the liquid mass present in the simulation.

X*Y ∫ρ(z)dz = mass

XYdz ∑ρ(i) = mass

where X, Y and dz are the dimensions in X, Y and the slab thickness in Z.

XYdz = 40.92 Angstrom *41.877 Angstrom *0.3191 Angstrom *1e-24 (cm/Angstrom)^3

40.92*41.877 0.31911e-24 ∑ρ(i) = mass (in grams)

When I take only calculate the density profile for 1 frame (effectively NOT averaging over many frames), I obtain the appropriate mass of the liquid system ( 2.0966e-020 grams). When I average over 200000 femtoseconds however, I get a much lower than expected value (1.2630e-020 grams)

Is this to be suspected ? Possibly round-off errors due to the huge number of averages I am completing ? I am surprised though, that I do not obtain the same mass.

Thank you for any help,

Tom

also, as a sidenote, I realize I did not specify really what the system looks like. I have 2 slabs of metal (parallel to the XY plane) confining 2.0966e-020 grams of liquid. There are periodic boundary conditions in x and y, and an aperiodic boundary condition in z (the slabs prevent the escape of the liquid).

Thank you

If the integral of the curve is different by nearly 2x, then the
curves must look radically different. Which one is right?
Where did the atoms go if they were not included in the
spatial profile? I doubt round-off could be an issue
since these values are accumulated in double precision.
You can also play with the Nfreq param of fix ave/spatial
and average over smaller windows (or a running average)
to see when things go bad.

Steve

These are indeed good questions. The problem arises after I compress the metal plates. But I do not know how that could change the amount of liquid seen by the density profile.

I have calculated the total mass of the liquid after compression using

compute 1 liquid property/atom mass
compute 2 liquid reduce sum c_m

and I receive the correct value, so no atoms have disappeared from the simulation box, just from the density profile…I am at a loss at where they could have gone.

In any case, thank you for your help Steve.

Tom.

I suggest also outputting number density and doing it
in “ave one” mode, not running, and verify that
you count all the atoms, both at the start and
end of the run.

Steve

Interesting.

Summing the count does give that all of the atoms in the liquid are accounted for in the density profile…I have also dumped the atomic masses, and they are correct.

I’m thinking I may have to dump all the liquid atomic positions, and calculate myself the density profile in order to figure out what exactly is going wrong.

So now do mass density (instead of num density)
the same way (ave one), and verify it is correct.

Then move on to ave running for both densities and
see what goes wrong. Maybe there is some
bug or round-off issue I am not aware of, since
the mass density values are much than counts.

Steve

Thanks a lot Steve, I’ll do this right now.

So I did the same thing (ave one), but this time with mass density, and the calculated mass is incorrect.

That being said, this may well be an already resolved bug, as I just realized my institution is using the June 2010 version of LAMMPS. I am not an administrator, so I cannot update the version being used on our cluster.

In any case, I know that the density number works, and I can most definitely use this output.

Thanks for your input Steve.

-Tom

So I did the same thing (ave one), but this time with mass density, and the calculated mass is incorrect.

That being said, this may well be an already resolved bug, as I just just realized my institution is using the June 2010 version of LAMMPS. I am not an administrator, so I cannot update the version being used on our cluster.

In any case, I know that the density number works, and I can most definitely use this output.

Thanks for your input Steve.

-Tom

So I did the same thing (ave one), but this time with mass density, and the
calculated mass is incorrect.

That being said, this may well be an already resolved bug, as I just
realized my institution is using the June 2010 version of LAMMPS. I am not
an administrator, so I cannot update the version being used on our cluster.

sure you can compile your own version on that cluster and use it.
there is nothing in LAMMPS that requires any administrator privileges
to work. people do this all the time.

axel.

This was indeed an bug associated with the version of LAMMPS I was using. I have updated LAMMPS and it works like a charm. I did not think to look at the version until a bit late in the game :wink:

Thank you for all your input !