Comment vous faites pivoter un système de particules (une matrice 2d de taille = num_particles * 3) de telle sorte que certaines entrées deviennent nulles

Question

suppose donc que vous avez un système de particules données par

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

Je veux tourner la système de sorte que la première particule se déplace vers l'origine, soit x1 = 0, y1 = 0, z1 = 0, la deuxième particule se déplace vers l'axe z, c'est-à-dire nouvelles coordonnées x2 = 0, y2 = 0, z2 = nouveau z2, et enfin la troisième particule se déplace vers le plan yz, c'est-à-dire, x3 = 0, y3 = nouveau y3, z3 = nouveau z3. L'important est que la distance entre toutes les particules doit être préservée.

J'ai essayé d'utiliser Givens Rotation pour mettre à zéro les coordonnées que j'ai spécifiées ci-dessus mais cette méthode change la distance entre les particules. Je code dans Fortran 90.

Ajouté: Voici un sous-programme que j'appelle la contrainte. J'ai essayé en construisant des matrices de rotation pour faire tourner le système comme je l'ai décrit ci-dessus. Comme prévu je reçois les zéros que je veux. Mais quand je mesure les distances entre les particules après appel de la contrainte, elles ne sont pas les mêmes qu'avant de les appeler (en fait, je calcule l'énergie du système qui est invariante par translation et rotation car elle ne dépend que de la séparation des particules)

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

Pendant que j'étais En écrivant la réponse, vous avez inclus votre code, qui semble être basé sur la rotation du corps rigide. Parce que mon code ci-dessous est également basé sur une rotation de corps rigide, je passerai des explications détaillées; donc s'il vous plaît comparer les deux codes si nécessaire (FYI, dans mon cas, j'effectue des rotations séquentielles Rz -> Ry -> Rz, tel que défini par les angles d'Euler).

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 

---

EDIT

Réécrire les matrices de rotation en coordonnées cartésiennes donne

!> 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) 

Answer 2

ne suis pas sûr étendre au nième point, mais pour la 1ère 3 passe par votre commande:

new p1 = (0,0,0) - donné

new p2 = (0,0,dist(p1,p2))

new p3 = (0,sin(A)*dist(p1,p3),cos(A)*dist(p1,p3))

où A est l'angle entre les côtés p1-p3 et p1-p2 qui peuvent être trouvés en utilisant la loi des cosinus (puisque vous avez la longueur entre les 3 points):

A = arccos((dist(1,3)*dist(1,3)+dist(1,2)*dist(1,2)-dist(2,3)*dist(2,3))/ (2*dist(1,3)*dist(1,2)))

Answer 3

Attend ma réponse d'origine a été supprimé ... envisager de changer votre code vers:

! Translating the whole system so that the first particle at the origin 
TransX = pos(1,1) 
TransY = pos(1,2) 
TransZ = pos(1,3) 
IF(pos(1,1) .NE. 0.0d0) THEN 
    pos(:,1) = pos(:,1) - TransX 
END IF 
IF(pos(1,2) .NE. 0.0d0) THEN 
    pos(:,2) = pos(:,2) - TransY 
END IF 
IF(pos(1,3) .NE. 0.0d0) THEN 
    pos(:,3) = pos(:,3) - TransZ 
END IF