Question about NPT integration

Greetings LAMMPS list/developers.

I have been attempting to diagnose some issues in a simulation related to instability in the NPT ensemble. While I do not think the integrator is at fault, I have been examining the code in fix_nh.cpp to ascertain what behavior leads to the instability. The following code struck me as odd:

void FixNH::compute_deviatoric()
{
  // generate upper-triangular part of h*sigma*h^t
  // units of fdev are are PV, e.g. atm*A^3 
  // [ 0 5 4 ]   [ 0 5 4 ] [ 0 5 4 ] [ 0 - - ]
  // [ 5 1 3 ] = [ - 1 3 ] [ 5 1 3 ] [ 5 1 - ]
  // [ 4 3 2 ]   [ - - 2 ] [ 4 3 2 ] [ 4 3 2 ]
  
  double* h = domain->h;
   
  fdev[0] = 
    h[0]*(sigma[0]*h[0]+sigma[5]*h[5]+sigma[4]*h[4]) +
    h[5]*(sigma[5]*h[0]+sigma[1]*h[5]+sigma[3]*h[4]) +
    h[4]*(sigma[4]*h[0]+sigma[3]*h[5]+sigma[2]*h[4]);
  fdev[1] = 
    h[1]*(              sigma[1]*h[1]+sigma[3]*h[3]) +
    h[3]*(              sigma[3]*h[1]+sigma[2]*h[3]);
  fdev[2] = 
    h[2]*(                            sigma[2]*h[2]);
  fdev[3] = 
    h[1]*(                            sigma[3]*h[2]) +
    h[3]*(                            sigma[2]*h[2]);
  fdev[4] = 
    h[0]*(                            sigma[4]*h[2]) +
    h[5]*(                            sigma[3]*h[2]) +
    h[4]*(                            sigma[2]*h[2]);
  fdev[5] = 
    h[0]*(              sigma[5]*h[1]+sigma[4]*h[3]) +
    h[5]*(              sigma[1]*h[1]+sigma[3]*h[3]) +
    h[4]*(              sigma[3]*h[1]+sigma[2]*h[3]);
}

I cannot determine why the computations of matrix elements fdev[0,1,2] are in any way different. Symmetry would appear to dictate that the matrix fdev[i][j] = h[i][l] * sigma[l][k] * h[j][k] (using Einstein summation) so that the same number of terms appears in each expression—some repeated terms should occur, but this will lead to multiplicative factors in front of said terms. Moreover, exchanging (e.g.) x for z results in a different expression for fdev, and thus a different equation of motion for the box. I note that Voigt notation is used, but in this notation h is a symmetric tensor and not an upper-triangular one, so the point still stands.

Surely there is a simple reason for this I am missing, but I would appreciate an explanation. Thank you.

Aidan can comment on this.

Steve

Here is a simple explanation: The h matrix is not symmetric. In the
special case of LAMMPS it is upper triangular, so we can represent it
using Voigt notation.