[lammps-users] NPH: cell matrix h & cell momenta p_g

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

lammps uses in its nph/npt/nvt schemes the algorithm proposed in Shinoda, Shiga, and Mikami, Phys Rev B, 69, 134103 (2004).
Basically, the cell matrix h = [a b c], with the three vectors a, b and c denoting the three edges of the cell box, is transformed into a p_g, the “modularly invariant form of the cell momenta”.

The quantity p_g is symmetric and lammps uses the governing equations given in Eq. (1) of mentioned paper.

I have following questions since in FixNH::nh_omega_dot() a dilation[6] is calculated

  1. How is the cell momentum p_g linked to the cell matrix h? Can they be transformed to each other?

  2. In FixNH::remap() the box is adjusted depending to new size/shape, which is governed by dilation[6].
    Why ( in FixNH::nh_omega_dot() ) are dilation[0]-[2] calculated like dilation[i] = exp(dtoomega_dot[i])
    and dilation[3]-[5] like dilation[i] = dto    omega_dot[i]?

Thanks for help and thank you for the answers to my last questions,
Manfred

2

I'll let Aidan answer these Qs as he did the implementation
of the latest NPT version.

Steve

3

Manfred,

Here are responses to your two questions.

1) p_g = W*omega_dot, as in the code comment below:

// barostat energy is equivalent to Eq. (8) in
  // Martyna, Tuckerman, Tobias, Klein, Mol Phys, 87, 1117
  // Sum(0.5*p_omega^2/W + P*V),
  // where N = natoms
  // p_omega = W*omega_dot
  // W = N*k*T/p_freq^2
  // sum is over barostatted dimensions

2) The reason for the difference is that the dilations correspond to
increments of the h-matrix according to the following ODEs:

h_dot[0] = omega_dot[0]*h[0]
h_dot[1] = omega_dot[1]*h[1]
h_dot[2] = omega_dot[2]*h[2]

h_dot[3] = omega_dot[3]*h[2]
h_dot[4] = omega_dot[4]*h[2]
h_dot[5] = omega_dot[5]*h[1]

The indices are those used for omega_dot[6] and dilation[6] in the code:
(psuedo-Voigt notation): 0 = xx, 1 = yy, 2 = zz, 3 = yz, 4 = xz, 5 = zy

Aidan

4

Hello,

thank you for the answers.
That brings me to another question:

The implementation of h is very smart, vector a is always pointing in x direction and obtaining h with
a = (xhi-xlo,0,0), b = (xy,yhi-ylo,0) and c = (xz,yz,zhi-zlo) leads always to six non-zero components.

When h_dot is updated it follows with p_g = W*omega_dot:
h_dot = p_g/W * h
[ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ]
[ - 1 3 ] = [ - 1 3 ] * [ - 1 3 ]
[ - - 2 ] [ - - 2 ] [ - - 2 ]
( “-” indicates zero).
Shouldn’t it be something like:
[ 0 5 4 ] [ 0 5 4 ] [ 0 5 4 ]
[ 8 1 3 ] = [ 5 1 3 ] * [ - 1 3 ]
[ 7 6 2 ] [ 4 3 2 ] [ - - 2 ]
p_g is symmetric and hence all positions of h_dot change (h_dot is non-symmetric).

I know that the dilation[6] is taken to update the cell shape, but h_dot is even different in positions 0-5 in the preceding two equations.

Thanks, Manfred

2010/6/7 Thompson, Aidan <[email protected]>