# 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!