algorithm - 如何旋转粒子系统(大小为 =num_particles*3 的二维矩阵)，以便某些条目变为零

Question

所以假设你有一个由

给出的粒子系统

pos = [x1, y1, z1,
       x2, y2, z2,
            .
            .
       xn, yn , zn]

我想旋转系统，使第一个粒子移动到原点，即 x1 = 0 , y1 =0 , z1=0，第二个粒子移动到 z 轴，即新坐标 x2 = 0, y2 = 0, z2 = new z2, 最后第三个粒子移动到yz平面，即x3=0, y3 = new y3, z3 = new z3。重要的是所有粒子之间的距离必须保持不变。

我尝试使用 Givens Rotation用于将我上面指定的坐标归零，但此方法会改变粒子之间的距离。我正在用 Fortran 90 编码。

添加:这是一个我称之为约束的子程序。我尝试通过构建一些旋转矩阵来旋转系统，如上所述。正如预期的那样，我得到了我想要的零。但是当我在调用约束后测量粒子之间的距离时，它们与调用它之前不一样(实际上我所做的是我计算系统的能量，它在平移和旋转下是不变的，因为它只取决于粒子分离)

SUBROUTINE constraint(pos)
REAL(KIND=dp), DIMENSION(np,3), INTENT(INOUT) :: pos
REAL(KIND=dp) :: r1, r2
REAL(KIND=dp), DIMENSION(3,3) :: rotMatrix
!------------------
! Translating the whole system so that the first particle at the origin
IF(pos(1,1) .NE. 0.0d0) THEN
    pos(:,1) = pos(:,1) - pos(1,1)
END IF  
IF(pos(1,2) .NE. 0.0d0) THEN
    pos(:,2) = pos(:,2) - pos(1,2)
END IF
IF(pos(1,3) .NE. 0.0d0) THEN
    pos(:,3) = pos(:,3) - pos(1,3)
END IF

! First rotation: Roates the whole system so that the second particle is on
! the z-axis
IF(pos(2,1) .NE. 0.0d0 .OR. pos(2,2) .NE. 0.0d0) THEN
    r1 = NORM2(pos(2,:))
    r2 = NORM2(pos(2,1:2))
    r2 = r2*r2
    rotMatrix(1,1) = ( pos(2,2)*pos(2,2) + ( pos(2,1) * pos(2,1) * pos(2,3) ) /r1 ) / r2
    rotMatrix(1,2) = pos(2,1)* pos(2,2) * (-1.0d0 + pos(2,3)/r1) / r2
    rotMatrix(1,3) = - pos(2,1) / r1
    rotMatrix(2,1) = rotMatrix(1,2)
    rotMatrix(2,2) = ( pos(2,1)*pos(2,1) + ( pos(2,2) * pos(2,2) * pos(2,3) ) /r1 ) / r2
    rotMatrix(2,3) = - pos(2,2) / r1
    rotMatrix(3,1) =  pos(2,1) / r1
    rotMatrix(3,2) =  pos(2,2) / r1
    rotMatrix(3,3) =  pos(2,3) / r1

    pos = MATMUL( pos, TRANSPOSE(rotMatrix) )
END IF

! Second rotation: Roates the whole system around the z-axis so that the 
! third particle is on the zy-plane
! the z-axis
IF( pos(3,1) .NE. 0.0d0 ) THEN
    r1 = NORM2(pos(3,1:2))

    rotMatrix(1,1) = pos(3,2) / r1
    rotMatrix(1,2) = - pos(3,1) / r1
    rotMatrix(1,3) = 0.0d0
    rotMatrix(2,1) = pos(3,1) / r1
    rotMatrix(2,2) = - pos(3,2) / r1
    rotMatrix(2,3) = 0.0d0
    rotMatrix(3,1) = 0.0d0
    rotMatrix(3,2) = 0.0d0
    rotMatrix(3,3) = 1.0d0

    pos = MATMUL( pos, TRANSPOSE(rotMatrix) )
END IF
END SUBROUTINE constraint

Answer 1

在我写答案时，您已经包含了您的代码，它似乎是基于刚体旋转的。因为我下面的代码也是基于刚体旋转的，所以就不详细解释了；所以请在必要时比较这两个代码(仅供引用，在我的例子中，我执行顺序 Rz -> Ry -> Rz 旋转，由欧拉角定义)。

program rotation
    implicit none
    integer, parameter :: N = 10, x=1, y=2, z=3
    real, parameter :: pi = acos(-1.0)
    real :: pos( 3, N ), alpha, beta, gamma, phi, ref( 3 ), rot(3,3)
    integer i

!> Initial coordinates.
    do i = 1, N
        phi = 2.0 * pi / N * (i - 1)
        pos( :, i ) = [ cos( phi ), sin( phi ), 0. ] &
                      * ( 3.0 + 2.0 * mod(i,2) ) * 0.55
    enddo
    pos(2,:) = pos(2,:) + 5.0

!> Translate the system such that pos(:,1) = 0.
    ref(:) = pos( :, 1 )
    do i = 1, N
        pos( :, i ) = pos( :, i ) - ref(:)
    enddo

!> Get the polar coordinates of pos(:, 2).
    beta  = acos( pos( z, 2 ) / norm2( pos(:, 2) ) )  !! in [0,pi]
    alpha = atan2( pos( y, 2 ), pos( x, 2 ) )         !! in [-pi,pi]

!> Apply Rz( -alpha ).
    pos = matmul( get_Rz( -alpha ), pos )

!> Apply Ry( -beta ).
    pos = matmul( get_Ry( -beta ), pos )

!> Get the azimuthal angle of pos(:, 3).
    gamma = atan2( pos( y, 3 ), pos( x, 3 ) )

!> Apply Rz( -gamma + pi/2 ).
    pos = matmul( get_Rz( -gamma + pi/2 ), pos )

!> Result.
    print *, "new coord:"
    do i = 1, N
        print "(3f10.5)", pos( :, i )
    enddo

    rot = matmul( get_Rz( -gamma + pi/2 ), &
          matmul( get_Ry( -beta ), get_Rz( -alpha ) ) )

    print *, "full rotational matrix (to be applied after translation):"
    do i = 1, 3
        print "(3f10.5)", rot( i, : )
    enddo

contains

    function get_Rz( ang ) result( R )
        real :: ang, R(3,3)
        R( 1, : ) = [ cos( ang ), -sin( ang ), 0. ]
        R( 2, : ) = [ sin( ang ),  cos( ang ), 0. ]
        R( 3, : ) = [    0.,          0.,      1. ]
    endfunction

    function get_Ry( ang ) result( R )
        real :: ang, R(3,3)
        R( 1, : ) = [  cos( ang ), 0., sin( ang ) ]
        R( 2, : ) = [    0.,       1.,    0.      ]
        R( 3, : ) = [ -sin( ang ), 0., cos( ang ) ]
    endfunction
endprogram

---

编辑

根据笛卡尔坐标重写旋转矩阵给出

!> Apply Ry( -beta ) * Rz( -alpha ).
    p(:) = pos( :, 2 )
    r1 = norm2( p(:) )
    L  = norm2( p( 1:2 ) )
    Lr = L * r1

    rot( 1, : ) = [ p(z)*p(x) / Lr, p(z)*p(y) / Lr, - L / r1   ]
    rot( 2, : ) = [    - p(y) / L,       p(x) / L,     0.      ]
    rot( 3, : ) = [      p(x) / r1,      p(y) / r1, p(z) / r1  ]

    pos = matmul( rot, pos )

!> Apply Rz( -gamma + pi/2 ).
    p(:) = pos( :, 3 )
    L = norm2( p( 1:2 ) )
    rot( 1, : ) = [  p(y) / L, p(x) / L, 0. ]
    rot( 2, : ) = [ -p(x) / L, p(y) / L, 0. ]
    rot( 3, : ) = [   0.,        0.,     1. ]

    pos = matmul( rot, pos )