[lammps-users] Density of ellipsoidal particles

Dear LAMMPS users,

I have a question on how to set the density for ellipsoidal particles that results in correct masses.

Normally in the data file, for ellipsoidal particles, you set the ellipsoid mass by its density, which is mass/(4/3piabc) with a, b, c are the semi-axis lengths. For unit type “real”, the mass has a unit of g/mol, distance has a unit of Angstrom, while density has a unit of (gram/cm^3). So, it means that density = mass/(4/3piabcN_A(10e-8)^3) with N_A is the Avogadro number and the (10e-8)^3 for the conversion of distance unit from Angstrom to cm.

However, I could not find a way to do this correctly. For example, for an ellipsoid CG model of benzene, mass is 78.11g/mol and a, b, c are ~2.25,2.25,0.75 Angstroms respectively. If you do the math correctly, it gives a density of 8.155316 g/(cm^3). However, when running a simulation, the mass in the dump file is not 78.11g/mol. Instead, when setting the density to ~4.911g/(cm^3) the mass is 78.11 g/mol in the dump file. More curiously, 8.155316/4.911 is 1.6055, which reminds me of the atomic mass unit (1.6055*10^(-27) kg)

Of course, I can set the mass manually in the input script, but I was wondering if I had done something incorrectly, or if there is a mistake in the code?

I had a look at the source code, and I can see that for 3D ellipsoids the mass is 4.0*MY_PI/3.0 * shape[0]*shape[1]*shape[2].

What LAMMPS version is this with?

Can you provide a minimal input deck demonstrating the problem and post this as a bug report issue on GitHub?




Dear Axel,

I just figured out it was a misunderstanding. The density set for ellipsoids has a different unit compared to the system density. In particular, for unit real, system density has a unit of g/cm^3, but the ellipsoid density has a unit of mass/(distance^3) which is g/mol/Angstrom^3 (this is specified in the manual but not very easy to be realised)

So, for calculation of ellipsoid density, the correct math is density = mass/((4/3)piabc) with a, b, c being the ellipsoidal semiaxis lengths. I tried this and it works perfectly!