Hi @noam.bernstein,
Please have a look at other threads on the topic, for example this one. The paper I link is previous work where we convergence issues in the NPT ensemble to compute the full elastic constant tensor.
I believe there were some corrections in the algorithm since (my PhD advisor told me so) but I never tested the new version and I still think some of our results hold: i.e. there is some impact from the coupling constant of the NH integrator on the sampling. This is something that also stemmed out from discussion with other people in the field. So I don’t think your analyses is incorrect but this requires careful processing to be sure about convergence anyway.
If you have issues with press/langevin please tell or open a PR on Github, I did the implementation and would be eager to correct any problem encountered.