I’m attempting to compute the velocity autocorrelation function (VACF) for the center of mass (CoM) of molecules using the latest LAMMPS version. I found the command compute vac all vacf com yes molecule yes, but it’s undocumented, so I’m unsure if it’s outdated or unsupported.
Here’s what I’ve tried for CoM velocities:
compute cc1 all chunk/atom molecule
compute vchunk all vcm/chunk cc1
But I’m stuck with getting fix ave/correlate to work with chunks. I know exporting velocities for postprocessing is an option, but I’d rather handle everything within LAMMPS. Could someone point me in the right direction for computing VACF with the latest features?
but it takes no arguments. So adding “com yes” or “molecule yes” will have no effect.
Since the VACF is a rather noisy function, why don’t you just define a group with the atom in each molecule that is closest to the center of mass and compute the VACF for that? Given the margin of noise and statistical fluctuations, the result should be indistinguishable. Or define the group with two atoms where the average position would be close to the center of mass?
You have to explain that and ideally provide some evidence of how this would produce incorrect results. Many times people claim that things would not work without proof.
Compute vacf makes a copy of the original velocities of all atoms in the selected group.
During the run, it computes the dot product between that original velocity vector and the current velocity vector and then averages this over all atoms in the group. That is how auto-correlation functions are defined. See, e.g., Autocorrelation page on SklogWiki - a wiki for statistical mechanics and thermodynamics
I’m not claiming that this would produce an incorrect result; I’m merely following what you just said about finding an atom close to the center of mass.
I understand how the VACF is computed. What I’m confused about is this: if I have two atoms in a group, will the VACF be computed using just these two atoms? In other words, is a group equivalent to my ensemble? Am I missing something?
But that makes no sense then. How would symmetry make this impossible? And what kind of symmetry?
I have just explained this. I don’t know how to do this differently. So if my way of explaining is not sufficient, you can always look at the source code and see what it does. That is the ultimate documentation. The source for compute vacf is rather short and easy to understand.
The issue is with a diatomic molecular fluid. We want to calculate the VACF of the centre of mass. If one uses the atomic VACF, one gets a big oscillation into negative VACF from the rotation. This is correct, of course, but its not what we want.
i.e. for setting a group of two atoms, we want the dot product of average velocity of all atoms in the group with its original value. Whereas (I think) the code provides the average of the dot products of each of the atoms with its original value.
[sum(v_i(t)).sum(v_i(0)] and not sum[v_i(t).v_i(0)] for i in the group.
Note that compute msd/chunk will calculate the mean squared displacement of a chunk’s centre of mass over time. Since displacement is the time integral of velocity, the MSD can be Fourier-transformed back into the velocity autocorrelation function relatively straightforwardly.