LAMMPS DumpReader returns incorrect triclinic box dimensions #3383
Description
Expected behavior
Triclinic box dimensions written in the lammpsdump format should be imported correctly.
Actual behavior
The DumpReader expects triclinic boxes to be written in the format:
if triclinic:
xlo, xhi, xy = map(float, f.readline().split())
ylo, yhi, xz = map(float, f.readline().split())
zlo, zhi, yz = map(float, f.readline().split())which then make up the three triclinic box vectors a = (xhi-xlo,0,0); b = (xy,yhi-ylo,0); c = (xz,yz,zhi-zlo).
But I believe this is not how recent (maybe all!?) lammps versions dump triclinic box dimensions. They give the dimensions of the orthogonal bounding box around the triclinic cell:
ITEM: BOX BOUNDS xy xz yz pp pp pp
xlo_bound xhi_bound xy
ylo_bound yhi_bound xz
zlo_bound zhi_bound yz
which are related to the dimensions we expect by:
xlo_bound = xlo + MIN(0.0,xy,xz,xy+xz)
xhi_bound = xhi + MAX(0.0,xy,xz,xy+xz)
ylo_bound = ylo + MIN(0.0,yz)
yhi_bound = yhi + MAX(0.0,yz)
zlo_bound = zlo
zhi_bound = zhi
This means that when importing lammpsdump files with triclinic boxes, the DumpReader reader interpretes the x and y dimensions of the orthogonal bounding box as the dimensions of the triclinic box, which results in elongated x and y dimensions. In my case, the bulk crystal I am studying suddenly had an interface in the x direction.
See relevant lammps doc page on the triclinic dump output style.
Code to reproduce the behavior
triclinic_lammps.tar.gz
I put my lammps .data .dump and .in in the attached zip. The box dimensions in the .data file read:
...
-0.321 16.831 xlo xhi
-0.123 25.958 ylo yhi
-0.045 12.993 zlo zhi
1.5067 -6.266 -0.421 xy xz yz
...
After simulating and importing the dump file I get:
>>>import MDAnalysis as mda
>>>u = mda.Universe("dump.atom", topology_format='LAMMPSDUMP')
>>>u.dimensions
array([ 24.925383, 26.547276, 14.473168, 93.076546, 115.656044, 86.746315])The length of the x dimension has changed from about 16 Angstrom to 25.
Current version of MDAnalysis
- Which version are you using? 2.0.0-dev0
- Which version of Python 3.7.10
- Which operating system? Ubuntu 18.04 bionic
- Lammps version 29 Oct 2020
Fix?
The fix would be easy, but I wanted to bring this problem up here before making a pull request. The lammps documentation on github where this orthogonal bounding box approach is documented reaches back to Nov 5, 2019, but I did not find any documentation yet for older versions.
MDAnalysis/coordinates/LAMMPS.py
...
if triclinic:
xlo_bound, xhi_bound, xy = map(float, f.readline().split())
ylo_bound, yhi_bound, xz = map(float, f.readline().split())
zlo, zhi, yz = map(float, f.readline().split())
xlo = xlo_bound - min(0.0, xy, xz, xy + xz)
xhi = xhi_bound - max(0.0, xy, xz, xy + xz)
ylo = ylo_bound - min(0.0, yz)
yhi = yhi_bound - max(0.0, yz)
box = np.zeros((3, 3), dtype=np.float64)
box[0] = xhi - xlo, 0.0, 0.0
box[1] = xy, yhi - ylo, 0.0
box[2] = xz, yz, zhi - zlo
xlen, ylen, zlen, alpha, beta, gamma = mdamath.triclinic_box(*box)
...