Issue with the lx for triclinic from the thermo command

Dear all,

When comparing the cell values printed with thermo command and the trajectory files in the atom format, I found different values for lx (if lx = xhi-xlo, and using the conversion in nomenclature for a triclinic system - lattice lengths and angles, when needed).

I tested different triclinic systems and compared the cell values as thermodynamic properties in the log file, with the box parameters printed in the atom format, and with the initial box from the data file. All the other cell parameters match but not the cella or lx for all the tested systems because the xmin value differs from the initial input. This doesn’t occur with a cubic system.

I am using the lammps version of (7 Aug 2019) but I got the same results with recent previous versions. I have attached a log file and the respective trajectory file in the atom format. I have also added the input files that reproduce this simulation.

Best regards,

Ingrid Padilla
North Carolina A&T State University

From log file:
Step TotEng Temp Press Volume Cella Cellb Cellc CellAlpha CellBeta CellGamma Density Lx Ly Lz Xy Xz Yz
1 -42.693332 0 796534.02 342.7945 5.48 6.76 9.2799996 90 94.330001 90 2.9491141 5.48 6.76 9.253512 4.14e-16 -0.700648 6.11e-16

From trajectory file:
ITEM: BOX BOUNDS xy xz yz pp pp pp
-7.0064800000000005e-01 5.4800000000000004e+00 4.1400000000000002e-16
0.0000000000000000e+00 6.7600000000000007e+00 -7.0064800000000005e-01
0.0000000000000000e+00 9.2535120000000006e+00 6.1099999999999997e-16

From data file
0.0000 5.48 xlo xhi
0.0000 6.76 ylo yhi
0.0000 9.253512 zlo zhi
4.14E-16 -0.700648 6.11E-16 xy xz yz

log.r_test.MIN (17.3 KB)

r_test.MIN.atom (76.7 KB)

data.r (1.45 KB)

in.r-298K1atm (1.89 KB)

reax.ffield_CaSiAlO_Pitman (20.3 KB)

lmp_control (1009 Bytes)

lx is not computed from xhi-xlo
xlo and xhi are the box boundaries, i.e. the smallest x and the largest x. their difference is only equal to the length for orthogonal boxes. with a tilted box this difference between the boundaries is larger depending on the amount of tilt and direction of tilt.

axel.

Thanks Axel for the quick reply. I still don’t understand why the origin of the box in x in the trajectory file is different from the initial input in the data file, while all of the other parameters remain the same.

I am writing a code to calculate the rdf for triclinic systems using PBC, when calculating the volume as the determinant of the parallelepiped I’m using the vectors a = (xhi-xlo,0,0); b = (xy,yhi-ylo,0); c = (xz,yz,zhi-zlo). because the differences in xlo, my calculation differs from the thermodynamic properties from lammps. Should I then assume that xlo = 0?

My results are from a hydrostatic compression test. Is LAMMPS compressing the system in the x direction assuming the origin in zero or in the printed xlo from the trajectory file?

Regards,

Ingrid.

Thanks Axel for the quick reply. I still don’t understand why the origin of the box in x in the trajectory file is different from the initial input in the data file, while all of the other parameters remain the same.

I am writing a code to calculate the rdf for triclinic systems using PBC, when calculating the volume as the determinant of the parallelepiped I’m using the vectors a = (xhi-xlo,0,0); b = (xy,yhi-ylo,0); c = (xz,yz,zhi-zlo). because the differences in xlo, my calculation differs from the thermodynamic properties from lammps. Should I then assume that xlo = 0?

please see below, for some C code that shows how you need to adjust the box bounds of a LAMMPS file for triclinic boxes. this is what VMD uses.

/* adjust box bounds */
xdelta = MIN(0.0f,xy);
xdelta = MIN(xdelta,xz);
xdelta = MIN(xdelta,xy+xz);
xlo -= xdelta;

xdelta = MAX(0.0f,xy);
xdelta = MAX(xdelta,xz);
xdelta = MAX(xdelta,xy+xz);
xhi -= xdelta;

ylo -= MIN(0.0f,yz);
yhi -= MAX(0.0f,yz);

then to get the box length and angles:

/* convert bounding box info to real box lengths /
A = xhi-xlo;
ylohi = yhi-ylo;
B = sqrt(ylohi
ylohi + xyxy);
C = sqrt((zhi-zlo)
(zhi-zlo) + xzxz + yzyz);

/* compute angles from box lengths and tilt factors /
alpha = cosangle2deg((xy
xz + ylohi*yz)/(B * C));
beta = cosangle2deg(xz/C);
gamma = cosangle2deg(xy/B);

with:

float cosangle2deg(double val)
{
if (val < -1.0) val=-1.0;
if (val > 1.0) val= 1.0;
return (float) (90.0 - asin(val)*90.0/M_PI_2);
}

FYI, here are the results of loading your data file and your atom style dump into VMD.

as you can see, the length and angle of the box vectors match:

vmd > topo readlammpsdata data.r charge
Info) parsing LAMMPS header.
Info) readlammpsdata: detected triclinic cell.
Info) parsing LAMMPS Masses section.
Info) parsing LAMMPS Atoms section with style ‘charge’.
Info) applying atoms data. sorted by atom id.
Info) Analyzing structure …
Info) Atoms: 26
Info) Bonds: 0
Info) Angles: 0 Dihedrals: 0 Impropers: 0 Cross-terms: 0
Info) Bondtypes: 0 Angletypes: 0 Dihedraltypes: 0 Impropertypes: 0
Info) Residues: 26
Info) Waters: 0
Info) Segments: 1
Info) Fragments: 26 Protein: 0 Nucleic: 0
4
vmd > pbc get
{5.480000 6.760000 9.280000 90.000000 94.330002 90.000000}
vmd > mol new r_test.MIN.atom type lammpstrj
Info) Using plugin lammpstrj for structure file r_test.MIN.atom
lammpsplugin) New style dump with 5 data fields. Coordinate data flags: 0x102
lammpsplugin) Reconstructing atomic coordinates from fractional coordinates and box vectors.
lammpsplugin) Detected triclinic box.
Info) Using plugin lammpstrj for coordinates from file r_test.MIN.atom
Info) Determining bond structure from distance search …
Info) Analyzing structure …
Info) Atoms: 26
Info) Bonds: 9
Info) Angles: 0 Dihedrals: 0 Impropers: 0 Cross-terms: 0
Info) Bondtypes: 0 Angletypes: 0 Dihedraltypes: 0 Impropertypes: 0
Info) Residues: 17
Info) Waters: 0
Info) Segments: 1
Info) Fragments: 17 Protein: 0 Nucleic: 0
5
vmd > Info) Finished with coordinate file r_test.MIN.atom.
pbc get
{5.480000 6.760000 9.280000 90.000000 94.330002 90.000000}

Thanks for the tests, the inputs and the code. This information is very helpful and clarified my doubts.

Ingrid