Momentum contribution to stress/mop & stress/mop/profile

Hello all,

I am studying the source code of the compute stress/mop and stress/mop/profile commands and trying to understand the implementation for the momentum contribution to the stress tensor. It seems to me there is a disparity between the two commands when dealing with finite-size particles. My question is how stress/mop handles the case of finite size particles which have rmass mass per particle instead of mass mass per type?
In compute_stress_mop.cpp the momentum contribution is computed in lines computed in lines 460-472 :

        if (rmass) {
          vcross[0] = vi[0]-fi[0]/rmass[i]/2.0*ftm2v*dt;
          vcross[1] = vi[1]-fi[1]/rmass[i]/2.0*ftm2v*dt;
          vcross[2] = vi[2]-fi[2]/rmass[i]/2.0*ftm2v*dt;
        } else {
          vcross[0] = vi[0]-fi[0]/mass[itype]/2.0*ftm2v*dt;
          vcross[1] = vi[1]-fi[1]/mass[itype]/2.0*ftm2v*dt;
          vcross[2] = vi[2]-fi[2]/mass[itype]/2.0*ftm2v*dt;
        }
        values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
        values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
        values_local[m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;

In compute_stress_mop.cpp the momentum contribution is computed in lines computed in lines 458-476 :

          if (rmass) {
            double vcross[3];
            vcross[0] = vi[0]-fi[0]/rmass[i]/2*ftm2v*dt;
            vcross[1] = vi[1]-fi[1]/rmass[i]/2*ftm2v*dt;
            vcross[2] = vi[2]-fi[2]/rmass[i]/2*ftm2v*dt;

            values_local[ibin][m] += rmass[i]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
            values_local[ibin][m+1] += rmass[i]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
            values_local[ibin][m+2] += rmass[i]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
          } else {
            double vcross[3];
            vcross[0] = vi[0]-fi[0]/mass[itype]/2*ftm2v*dt;
            vcross[1] = vi[1]-fi[1]/mass[itype]/2*ftm2v*dt;
            vcross[2] = vi[2]-fi[2]/mass[itype]/2*ftm2v*dt;

            values_local[ibin][m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
            values_local[ibin][m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
            values_local[ibin][m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
          }

In the stress/mop command, the mass is retrieved both by mass and rmass vectors while in the stress/mop/profile, the mass is retrieved only in the rmass vector. I would expect something like the following block of code for the stress/mop case:

        if (rmass) {
          vcross[0] = vi[0]-fi[0]/rmass[i]/2.0*ftm2v*dt;
          vcross[1] = vi[1]-fi[1]/rmass[i]/2.0*ftm2v*dt;
          vcross[2] = vi[2]-fi[2]/rmass[i]/2.0*ftm2v*dt;

          values_local[m] += rmass[i]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
          values_local[m+1] += rmass[i]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
          values_local[m+2] += rmass[i]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;

        } else {
          vcross[0] = vi[0]-fi[0]/mass[itype]/2.0*ftm2v*dt;
          vcross[1] = vi[1]-fi[1]/mass[itype]/2.0*ftm2v*dt;
          vcross[2] = vi[2]-fi[2]/mass[itype]/2.0*ftm2v*dt;

          values_local[m] += mass[itype]*vcross[0]*sgn/dt/area*nktv2p/ftm2v;
          values_local[m+1] += mass[itype]*vcross[1]*sgn/dt/area*nktv2p/ftm2v;
          values_local[m+2] += mass[itype]*vcross[2]*sgn/dt/area*nktv2p/ftm2v;
        }

The initial commit that made these commands compatible with per-atom masses dates back to 2018:

Kind regards
E. Voyiatzis

This is obviously a bug. Should be fixed in commit make handling of masses consistent and simplify code · akohlmey/lammps@d813519 · GitHub

There is also potential to simplify the code (so the compiler may do a better job optimizing it).

P.S.: a simple way to test this is to add fix property/atom rmass (or not) to enable (or distable) per-atom masses with a per-type system and then experiment with changing masses to see if the right and consistent changes are observed in the computed output.

Thanks a lot for the prompt response and your quick action to fix this bug!
I will follow your advice and experiment with fix property/atom rmass as my next step!