Is there way in LAMMPS prevent rescaling of frozen atoms at constant
pressure (NPT) without effecting the pressure regulation of the
remaining atoms?
Running npt simulations on a partial subset of the atoms in a
simulation effectively freezes the motion of the remaining atoms. IF
"dilate partial" is used, these remaining atoms are not rescaled
during pressure equilibration. In this specific case, the forces
acting on them (including long-range coulombic forces) should not
contribute to the virial. (However the forces they exert on the mobile
atoms should remain in the virial. The attached PDF file includes a
more detailed derivation. Feel free to ignore it. Please let me know
if it is not readable.)
Alternately, one can think of it another way: The forces that would
have been necessary to immobilize the atoms (to negate their movement)
should be included in the virial.
How does one do this in LAMMPS?
Rephrasing my question, what is the proper way to run LAMMPS (with the
correct virial) when using "dilate partial"?
Right now, the "setforce" fix does not seem to effect the virial.
(Under normal circumstances, this is actually what most people would
want. Unless "dilate partial" is used, setting the forces to zero
should not effect the virial. See equation 8 of the attached PDF file
for a detailed discussion.)
Thank you in advance for your feedback.
Cheers
Andrew
P.S.
The "rigid" fix might be another alternative, but vaguely remember it
has issues with NPT also, and is not appropriate for large periodic
objects like walls and slabs.
"fix wall" is useful, but it is not general enough to handle the
systems that we usually run into (immobilized atoms may be charged, or
interact with the other atoms using more complicated potentials).
"neigh_modify exclude" does effect the virial, but I don't think it
solves this problem. (We need to exclude the forces which act on the
immobilized atoms from the virial, but not the forces in the other
direction.)
I suppose I could edit the LAMMPS code to manually subtract the
contribution that these forces make on the immobilized atoms from the
virial as it's currently calculated (...in order to remove
contributions which should not have been added). Please let me know
if you think this is a good/bad strategy.
I know another research group which uses LAMMPS which ran into the
same problem. I think they solved the problem by gluing the atoms
together with strong harmonic pairwise bonds, and ran the simulation
in the ordinary way (without "dilate partial" or "setforce", and using
"fix 1 all npt ...").
virial_immobile_atoms_sm.pdf (101 KB)