MC/TFMC fixrot code question

Dear community,

I am having a look at the MC/fix_tfmc.cpp code and I have a question about how the angular momentum is “zeroed”.

I understand the process but I am surprised that the angular momentum is calculated using the already updated positions (see code below). The positions of the atoms (x array here) already have been updated by the Monte-Carlo code at line 222. In my understanding what is computed here is the angular momentum of the collection of particles using the new positions but still the old displacements (xd array, which is computed from line 222), which seems odd to me.

In a Python code, I can zero the angular momentum using x in the formula below as the current positions but not using the updated ones.

    double p[3];
    p[0] = p[1] = p[2] = 0.0;

    for (int i = 0; i < nlocal; i++) {
      if (mask[i] & groupbit) {
        domain->unmap(x[i],image[i],unwrap);
        dx = unwrap[0] - cm[0];
        dy = unwrap[1] - cm[1];
        dz = unwrap[2] - cm[2];
        if (rmass) massone = rmass[i];
        else massone = mass[type[i]];
        p[0] += massone * (dy*xd[i][2] - dz*xd[i][1]);
        p[1] += massone * (dz*xd[i][0] - dx*xd[i][2]);
        p[2] += massone * (dx*xd[i][1] - dy*xd[i][0]);
      }
    }

Of course, I might be missing something, any clarification would be highly appreciated.

I had to read up a little bit and please take into account that this is not my area of research.

Not really. Since in the tfMC method you don’t really have velocities, the displacements fill in for them. Those are position independent, so it should not matter, whether you look at positions before or after the move. What would be important for computing the angular momentum is that you are using a consistent center of mass and that is the case. The only potential issue I could see would be the sign for the displacements, but there it just seems to be the natural choice to take the “old to new” displacements. If you would do it the other way around, you would have to change the sign in the correction. Bottom line, you are trying to subtract the rotational component of the displacements and for those it seems natural to me to have the exact same flow of control that you see.

Yes, you are right, I forgot to update the positions when calculating the inertia tensor in my Python code, when it is correctly done the angular momentum is set to zero.

When computing the full pairwise distances between the angular reduced and non-angular reduced vectors, they end up being equal using the new updated positions (As done in LAMMPS). It seems that the same is not true for the other way (The one I was thinking about), hinting that the ‘dynamics’ would get impacted when doing this way. The whole thing seems counterintuitive to me, I should probably open a classical mechanics book again.

Thanks.