Rotation in TPI codes may be buggy in anisotropic systems
The following part in TPI.cpp
```
real angleX = 2*M_PI*dist(rng);
real angleY = 2*M_PI*dist(rng);
real angleZ = 2*M_PI*dist(rng);
rotate_conf(a_tp1-a_tp0, state_global->x.rvec_array()+a_tp0, nullptr,
angleX, angleY, angleZ);
```
should be replaced with
```
real angleX = 2*M_PI*dist(rng);
real angleY = asin(2*dist(rng)-1); //This line was modified
real angleZ = 2*M_PI*dist(rng);
rotate_conf(a_tp1-a_tp0, state_global->x.rvec_array()+a_tp0, nullptr,
angleX, angleY, angleZ);
```
to make the ration of molecule is uniform in terms of density in the solid angle space.
It is well known that if we draw uniform random numbers for angles in polar coordinates (theta and phi) and plot (sin theta cos phi, sin theta sin phi, cos theta), then we do not get uniform distribution in the solid angles. The distribution is denser near poles (see, e.g., <https://mathworld.wolfram.com/SpherePointPicking.html>). One way to avoid this bias is to draw random number for cos(theta) and phi.
The same issue occurs in the original codes. To illustrate it, we consider that the molecular coordinates (X’, Y’, and Z’) are first aligned with the laboratory coordinates X, Y, Z. We then go over how the unit vector on X’ axis distributes after random rotations by the original codes. The function of `rotate_conf` first rotates the molecule by angleX along X axis. At this point X and X’ are still aligned. Then X’Y’Z’ is rotated by angleY along Y axis. Now X’ is in XZ plane with angleY above the XY plane. The the routine rotates the molecule by angleZ along Z. As the radius of circle of the unit vector on X’ made by this last rotation is cos(angleY), we get denser points for angle Y~ 90 degree. Thus the rotation in the original codes is biased.
The infinitesimal solid angle for X’ to covers against infinitesimal angles d(angleY) d(angleZ) is
given by cos(angleY) d(angleY) d(angleZ) = d(sin(angleY)) d(angleZ). Thus we need to draw uniform random number for sin(angleY) instead of angleY itself in order to achieve the uniform distribution in the solid angles space.
I feel that the impact on results is negligible for isotropic systems, actually. However, for anisotropic systems, there could be some impact due to the biased sampling of rotation.
By the way, PME’s reciprocal contribution seems explicitly considered in TPI as of gromcs2019, which I admire. Who or which groups have contributed to the recent developmenst of TPI codes?
issue