Dear LAMMPS Users and Developers,
Has anyone tried to implement a general feature for swapping molecule COM velocities, like would be needed to work with SHAKE? This would be good for small rigid solvents, e.g., SPC/E has a good thermal conductivity (0.7-0.8 W/mK) and 2.0 fs timestep.
My understanding is that as long as the relative velocities of each atom are maintained, the system will still conserve energy with SHAKE. I modified fix thermal_conductivity to switch the COM velocity for atoms defined in an angle (e.g., 3-site rigid water). It seemed to work well (give linear temperature profile) in a few serial tests.
Tim
No one has done it so far as I know except you. It
would be nice to release in LAMMPS if what you
did is general. How do you know what atoms are
in the rigid body, and how to select them and how
to reset all their velocities together?
Steve
Steve, I think my way is not very general.
I only needed to switch velocity for 3-site water using SHAKE for bonds/angles. The strategy was to loop over local atoms (in a fix) looking for a particular type (oxygen), which I know are constrained. When oxygen is found, I retrieve the global ids of the other two atoms in the angle from angle_atom1, angle_atom2, angle_atom3 then use map() to find their local id. Then check the molecule’s KE (instead of atom’s KE).
When its time to switch velocity, the COM velocity of the original molecule is subtracted off all three atoms, and the new COM velocity is added on (per the Zhang reference).
Not sure how this would work in general. One idea might be to loop through local (and ghost?) atoms and sum the KE for each molecule, store atom IDs for lowest and highest molecular KE, then switch COM velocities.
Tim
Rather than use bonds/angles to find the other atoms, it
would probably be more general to just use the molecule ID.
Do you still need to identify a specific atom in the molecule (e.g. O),
or is the idea just to swap COM velocities of entire molecules?
Steve
Yes- using molecule ID would be much better.
The idea is to swap COM velocity of the entire molecule, but there has to be a check to see where the molecule is at. I used the O position in SPC/E for this, but the COM position would be better (although a bit more expensive). I’d guess that its not important for performance unless the swaps are happening very frequently or the solvent has alot of atoms per molecule. About one swap every 50-150 timesteps seems pretty reasonable for water.
Tim