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