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.