[lammps-users] confusion about scaled coordinates for triclinic cells

LAMMPS LAMMPS Mailing List Mirror
1

Hi! I am trying to calculate interparticle distances from a configuration-dump file with header

ITEM: TIMESTEP

14000000

ITEM: NUMBER OF ATOMS

100000

ITEM: BOX BOUNDS xy xz yz pp pp pp

-2.7226138904763260e+01 2.7158265908439908e+01 -4.3446974625327561e-01

-2.6695721711161067e+01 2.6681029457996679e+01 3.6659674992992375e-01

-2.6524840642222021e+01 2.6524840642222021e+01 -1.4692253164388701e-02

Using the definitions on the doc page (https://docs.lammps.org/Howto_triclinic.html), I get

xlo_bound = -27.2261
xhi_bound = 27.1583
xy = -0.43447
ylo_bound = -26.6957
yhi_bound = 26.681
xz = 0.366597
zlo_bound = -26.5248
zhi_bound = 26.5248
yz = 0.0146923

Converting these to the cell shape parameters gives

xlo = xlo_bound - min(0, xy, xz, xy + xz) = -26.7917
xhi = xhi_bound - max(0, xy, xz, xy + xz) = 26.7917
ylo = ylo_bound - min(0,yz) = -26.681
yhi = yhi_bound - max(0,yz) = 26.681

zlo = zlo_bound = -26.5248
zhi = zhi_bound = 26.5248

Using these definitions, the cell-edge (a, b, and c) vectors are

a = {xhi - xlo, 0, 0} = {53.5833, 0, 0}
b = {xy, yhi - ylo, 0} = {-0.43447, 53.3621, 0}
c = {xz, yz, zhi - zlo} = {0.366597, -0.0146923, 53.0497}

Then the h-matrix (the shape matrix whose column vectors are a, b, c) is

53.5833 -0.43447 0.366597
0 53.3621 -0.0146923
0 0 53.0497

and its inverse h_inv is

0.0186625 0.000151949 -0.000128924
0. 0.0187399 5.19007*10^-6
0 0 0.0188503

Here’s where I ran into trouble. I would expect the scaled coordinates of the cell corners {xlo, ylo, zlo} and {xhi, yhi, zhi} to be respectively

hinv {xlo, ylo, zlo} = {-.5, -.5, -.5}

and

hinv {xhi, yhi, zhi} = {.5, .5, .5}.

However, for the above hinv, they are

hinv {xlo, ylo, zlo} = {-.500634, -.500138, -.5}

and

hinv {xhi, yhi, zhi} = {.500634, .500138, .5}.

Thus the scaled cell isn’t a unit cube, so I can’t calculate interparticle distances using the usual textbook method of multiplying differences in unit-cube-scaled coordinates by h :frowning:

Any idea what’s going on here? I’ve checked the math several times so I assume I’ve misunderstood something, but can’t figure out what that might be…

Thanks in advance,
Rob

2

You are working with numbers that have about 5 digits precision. because you take differences of values of different magnitude, for smaller magnitude values that precision drops further.
So the precision on your results is disappointingly low, but I would suspect that if you output the data with full 15 digits precision that is available to floating point numbers, that error should drop.

3

Hi Axel. Thanks for the quick reply. And good suggestion, but that doesn’t seem to be the issue here. The C++ code I used to perform this calculation (attached in case you’re interested) reads in the data in the LAMMPS header to double precision. When I edited the code to output the h-matrix and scaled box bounds to high precision I got

rshoy@hoyoffice ONanalyze-triclinic_cell % ./ONanalyze-triclinic 100000 1 10 10 10 jammedconfig.fastcompression.purelyrepulsive

Double-precision h matrix:

53.58333831702 -0.434469746253276 0.366596749929924

0 53.3620589159934 -0.0146922531643887

0 0 53.049681284444

Scaled box bounds: xlo -0.50063445673268 xhi 0.50063445673268 ylo -0.500137665726013 yhi 0.500137665726013 zlo -0.5 zhi 0.5

Execution took 2 seconds.

So, ???

Thanks,
Rob

ONanalyze-triclinic.cpp (13.9 KB)

4

I don’t have time to look into the code, but could it be that you are missing a transpose somewhere?