Re: [ESPResSo-devel] Pressure virials changes

[Rudolf Weeber]

Hi Josh,
On Mon, Jun 24, 2013 at 11:46:19AM +0200, Josh Berryman wrote:
> .....I've come across a change in the pressure calculation on the main
> branch, fairly recent:
> 
> 
> Rudolf Weeber (Hi Rudolf!)
> commit 464f2bf67f82a79b56134fe39c4b3a85467df592
> 
> - #ifdef ROTATION
> -  virials.data.e[0] /= (6.0*volume*time_step*time_step);
> - #else
>    virials.data.e[0] /= (3.0*volume*time_step*time_step);
> - #endif
> 
> ....as you might expect, this means that my merged code (if ROTATION is
> defined) gives half the ideal contribution to the pressure that the old
> code did, which is a fairly big behaviour change.
> 
> Have I messed up my merging, or is there something meaningful happening
> here which I need to understand?
So, this change was already made in May 2011.
It removes the contribution of the rotational kinetic energy from the pressure.
See commit a46afe for the original.
The physical reasoning is that the amount of momentum transferred, when a 
particle runs into the wall does not depend on how fast the particle is 
rotating.
So, to my understanding, the current version should be correct.

The change does two things:
* Don't add contribution like viriaul += 1/2 J omega^2 (in 
add_kinetic_virials())
* Divide only by number of translational degrees of freedom (hence 3 instead of 
6).
In most cases, the simulation result will not change either way, because every 
degree of freedom carries 1/2kT. So if I have twice as many degrees of freedom 
and divide by 6 instead of 3, the result is the same.

It's strange, that the changes appear again in the commits merging the 
"rotation per particle"  feature, though.


Hope that helps,
Regards, Rudolf