Chimie Théorique » OpenPath

subroutine Std_ori(natom,rx,ry,rz,a,b,c,rmass)

  !c

  !c      rotational constants for all replicas by inline diagonalization

  !c      of inertial tensor

  !c

!!!!!!! Input

! natom    : number of atoms

! rx,ry,rz : arrays natom with cartesian coordinates of natom particles

! rmass    : vector of length natom with particle masses

!

!!!!!!! Output

! a,b,c    : rotational constants a, b, and c

!

!!!!!! Dependencies

! tensor: calculate the inertial tensor

! Jacobi, trie : diaganalization subroutines

  implicit none

  INTEGER, PARAMETER :: KINT=KIND(1)

  INTEGER, PARAMETER :: KREAL=KIND(1.0D0)

  INTEGER(KINT), INTENT(IN) :: Natom

  REAL(KREAL), INTENT(INOUT) :: rx(Natom),ry(Natom),rz(Natom)

  REAL(KREAL), INTENT(IN) :: rmass(Natom)

  REAL(KREAL), INTENT(OUT) :: a,b,c

  real(KREAL) :: Tinert(3,3),D(3),V(3,3),CoordG(3)

  real(KREAL) :: ca,sa,na,cb,sb,nb

  INTEGER(KINT) :: I

  REAL(KREAL) :: Vt,Rxt,Ryt

  call tensor(natom,rx,ry,rz,Tinert,CoordG,rmass)

  call Jacobi(Tinert,3,D,V,3)

  call Trie(3,D,V,3)

!1000 FORMAT(3F15.3)

!1010 FORMAT(4F15.3)

  a=D(1)

!     WRITE(*,1010) a,(V(i,1),i=1,3)

  b=D(2)

!     WRITE(*,1010) b,(V(i,2),i=1,3)

  c=D(3)

!     WRITE(*,1010) c,(V(i,3),i=1,3)

!C On commence par mettre le centre de masse a l'origine

  Do I=1,natom

!     WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I)

     rx(i)=rx(i)-CoordG(1)

     ry(i)=ry(i)-CoordG(2)

     rz(i)=rz(i)-CoordG(3)

!     WRITE(*,*) 'DBG Tensor', I, rx(I),rY(I),rz(I)

  END DO

  ! C Ensuite, on va faire coincider les axes d'inertie avec x,y,z

  ! C Pour etre coherent, et bien que ce ne soit pas le plus logique,

  ! C on classe les vecteurs dans le sens croissant, et on fait coincider

  ! C les axes dans l'ordre x,y,z pour l'instant.

  ! C Coincidation du premier axe avec x

  ! C Rotation autour de z pour l'amener dans le plan xOz

  ! C puis Rotation autour de y pour l'amener suivant x

  na=sqrt(V(1,1)*V(1,1)+V(2,1)*V(2,1))

  nb=1

  if (na.gt.1D-10) then

     ca=V(1,1)/na

     sa=-V(2,1)/na

  else

     ca=1

     sa=0

  END IF

  cb=na/nb

  sb=-V(3,1)/nb

  DO I=1,natom

     rxt=rx(i)

     rx(i)=rx(i)*ca-ry(i)*sa

     ry(i)=rxt*sa+ry(i)*ca

     rxt=rx(i)

     rx(i)=rx(i)*cb-rz(i)*sb

     rz(i)=rxt*sb+rz(i)*cb

  END DO

  !C on applique aussi la rotation aux deux autres axes...

  Vt=V(1,2)

  V(1,2)=V(1,2)*ca-V(2,2)*sa

  V(2,2)=Vt*sa+V(2,2)*ca

  Vt=V(1,2)

  V(1,2)=V(1,2)*cb-V(3,2)*sb

  V(3,2)=Vt*sb+V(3,2)*cb

  !c     WRITE(IOOUT,1010) b,(V(i,2),i=1,3)

  !c      Vt=V(1,3)

  !c      V(1,3)=V(1,3)*ca-V(2,3)*sa

  !c      V(2,3)=Vt*sa+V(2,3)*ca

  !c      Vt=V(1,3)

  !c      V(1,3)=V(1,3)*cb-V(3,3)*sb

  !c      V(3,3)=Vt*sb+V(3,3)*cb

  !c     WRITE(IOOUT,1010) c,(V(i,3),i=1,3)

  !C et on repart pour le 2e axe sur y

  !C commes les vecteurs sont orthonormes, il est

  !C deja dans le plan yOz d'ou une seule rotation

  !c      na=sqrt(V(2,2)*V(2,2)+V(3,2)*V(3,2))

  na=1

  nb=1

  ca=V(2,2)/na

  sa=-V(3,2)/na

  !c      write(IOOUT,*) 'ca:',ca,sa,cb,sb,na

  DO I=1,natom

     ryt=ry(i)

     ry(i)=ry(i)*ca-rz(i)*sa

     rz(i)=ryt*sa+rz(i)*ca

  END DO

  !C on applique aussi la rotation au dernier axe...

  !c      Vt=V(2,3)

  !c      V(2,3)=V(2,3)*ca-V(3,3)*sa

  !c      V(3,3)=Vt*sa+V(3,3)*ca

  !c      WRITE(IOOUT,1010) c,(V(i,3),i=1,3)

  !C et le 3e axe est deja sur z car  ils sont orthogonaux...

  return

end subroutine Std_ori