npt integration scheme

Hello lammps developers,

I have a question about the integrate equations for npt ensemble, for the
fix npt is use the Trotter factorization of the Liouville operator. If I
understand the scheme of integration can summarized as following:

nhc_press_integrate(); (thermostat and propagated the barostat velocities)

nhc_temp_integrate(); (thermostat and propagated the particles velocities)

nh_omega_dot(); (barostat propagation the barostat velocities)

nh_v_press(); (barostat scaled and propagation the the particles velocities)

nve_v(); (verlet propagation of particles velocities)

remap(); (barostat propagation of box)

nve_x(); (verlet propagation the particles positions)

remap(); (barostat propagation of box)

nve_v(); (verlet propagation of particles velocities)

nh_v_press(); (barostat propagation the the particles velocities)

nh_omega_dot(); (barostat propagation the barostat velocities)

nhc_temp_integrate(); (thermostat and propagated the particles velocities)

nhc_press_integrate(); (thermostat and propagated the barostat velocities)

I read the another list archives related to this issue, you mention the
Tuckerman paper (Tuckerman, Alejandre, Lopez-Rendon, Jochim, and Martyna,
J Phys A: Math Gen, 39, 5629 (2006)) the factorization is similar to, but
not identical to, the factorization described by Tuckerman et al., in the
motions equations to sampling the npt ensemble the barostat scaled the
particles positions similar to the particles velocities, in which part of
code take account this. I have another doubt the lammps use upper
triangular matrix to describe the box and symmetric virial stress tensor,
is possible use both together?, because with only one of those
constriction the cell rotation is eliminated.

Thanks
Samuel

This is a Q for Aidan.

Steve

  1. On the doc page, there is a section on the Liouville factorization that carefully defines the order of operations. Your listing (presumably obtained by reading the source code) seems to be consistent with that, except for one thing. The doc page is slightly incorrect in that the L_{\epsilon, 1} operator should symmetrically straddle the L_1 operator. The code does this by making two calls to remap() before and after the call to nve_x(). The biggest difference between LAMMPS and Tuckerman(2006) is in the combination of the non-hamiltonian barostat and the rRESPA propagators. If you are not doing rRESPA, this is not an issue.

  2. The instanteous stress tensor in MD is always symmetric. I think it follows from the assumption of point masses.

  3. The restriction of the box to an upper-triangular h-matrix effectively fixes the orientation of the cell in the reference frame. Generally, MD systems are rotationally invariant, and so this is not an issue.

Aidan