I’m implementing a fix (a biasing potential) and want to apply it in the NPT ensemble with a triclinic box. The issue I’m having is that my biasing potential does not yield a central force, which causes
r_if_j /= r_jf_i
for some i=x,y,z, j=x,y,z.
Correct me if I’m wrong, but fix nh appears to me to require the virial to be symmetric; this symmetry derived from the assumption that the energy is invariant with respect to distance-preserving transformations of the coordinates.
Is there a way to implement this fix without having to modify the pressure and integration fixes?
Aidan is the best developer to comment on this. He may
not be reading email thru the holidays. If you are
saying for the 3x3 viral tensor that Vij != Vji, you might
simply accumulate Vij += 0.5 * (vij + vii), since nothing
else in LAMMPS is going to use the values in the
asymmetric tensor, e.g. fix NPT does not use them.
All virial calculations in LAMMPS are based on the assumption of a symmetric tensor. There is a strong physical basis for this assumption. If you wish to violate it, you should ask yourself if this is really a good physical model. If the answer is yes, then I would like to see a previously published work that treats the physics and mathematics thoroughly.
The virial is symmetric provided the interatomic forces are central forces. Also, for strain energy computation the antisymmetric contributions of a matrix inner product vanish.
To make sure that I’m not misunderstanding anything:
The virial tensor is defined as
W_ij = sum_k -dU/dH_ik H_jk
where the columns of H are the box vectors.
If U is invariant under distance preserving transformations of the system, i.e.
U(r_1, r_2, …) = U(Q r_1 + c, Q r_2 + c, …)
where Q is in O(3) and c is a vector, then U can be expressed as a function of distances between positions and so the virial is symmetric,
W_ij = sum_k sum_a<b -dU/dr_ab dr_ab/dH_ik H_jk = sum_a<b -dU/dr_ab r_ab_i r_ab_j
where the sum is over all distances between atoms (and lattice vectors), and any other position that’s used in the potential.
Is there any other reason that the virial should be symmetric?
Also, in fix nh, is the box transformed by
just derivatives to the upper triangle of H (da_x/dt, db_x/dt, …, dc_z/dt) or
all derivatives, then rotated to make H upper triangular?
Is a point with a fixed (x, y, z) coordinate still at the same position relative to the dilation center?