Hi Petter,
Yes, you are right, thank you! Once I switched my calculations to “coulombtype = Cut-off” and turned on cut-off in my tool as well, the values started to match with high precision. Before that, I was also trying to use that scheme but got zero force but only now I realized it was due to the fact the diagonal is sqrt(3) times longer than box size so the zero force value was indeed possible:)
Also, I found the difference between ensembles: in the NVT and NPT configs, the cut-off values were different and the cut-off value is used in the PME mesh. I’m sorry for the confusion.
All in all, it looks like my calculations were correct up to the precision of the PME method in comparison to Coulomb law. I will also try to continue my production simulations with a small duration with a modified param file where electrostatics calculated in the expensive way (Coulomb) so I can make a fair analysis of the contribution of electrostatics to the total force.
Thank you so much, Petter!