Chimie Théorique » OpenPath

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

  ! c

  ! c      elements of inertial tensor for a set of atomic systems.

  ! c

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

  ! c      txx,tyy,tzz,txy,txz,tyz: elements of inertial tensors ; output

  ! c      a,b,c    :  returning center of mass coordinates.

  ! c      rmass    : vector of length natom with particle masses; input

  ! c      natom    : number of particles per system; input

  ! c      IOOUT     : unit number for messages, silent for IOOUT.le.0; input

  ! c      subroutines called: cmshft                   m. lewerenz jul/92

  ! c

  use Io_module

  IMPLICIT NONE

  REAL(KREAL), parameter :: zero=0.d0,one=1.d0

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

  real(KREAL) :: rx(Natom),ry(Natom),rz(Natom),rmass(Natom)

  real(KREAL) :: a,b,c,Tinert(3,3),CoordG(3)

  INTEGER(KINT) :: jover, I, J

  REAL(KREAL) :: Txx, Tyy, Tzz, Txy, Txz, Tyz, Txxi

  REAL(KREAL) :: Ryb, Rzc, Rxa, Rmj2, ryb2, rzc2, rxa2, Rmj

  if(natom.le.0) then

     if(IOOUT.gt.0) write(IOOUT,'(/2a/)')   &

          '  *** error, illegal array dimension(s) in tensor,', &

          ' return without action ***'

  else if(natom.eq.1) then

     a=rx(1)

     b=ry(1)

     c=rz(1)

     txx=zero

     tyy=zero

     tzz=zero

     txy=zero

     txz=zero

     tyz=zero

  else

     !c

     !c      first center of mass coordinates, then tensor elements

     !c

     call cmshft(natom,rx,ry,rz,rmass,a,b,c,0)

     CoordG(1)=a

     CoordG(2)=b

     CoordG(3)=c

1000 FORMAT(3F15.3)

     !c       WRITE(IOOUT,1000) (CoordG(i),i=1,3)

     jover=mod(natom,2)

     if(jover.eq.0) then

        txx=zero

        tyy=zero

        tzz=zero

        txy=zero

        txz=zero

        tyz=zero

     else

        rmj=rmass(1)

        ryb=b-ry(1)

        rzc=c-rz(1)

        tyz=-rmj*ryb*rzc

        rxa=a-rx(1)

        txz=(-rmj*rxa)*rzc

        txy=(-rmj*rxa)*ryb

        txx=rmj*((ryb*ryb)+(rzc*rzc))

        tyy=rmj*((rxa*rxa)+(rzc*rzc))

        tzz=rmj*((rxa*rxa)+(ryb*ryb))

     end if

!c

!c      force vectorization of inner loop if necessary (ibm-3090)

!c

     do j=jover+1,natom,2

        rmj=rmass(j)

        rmj2=rmass(j+1)

        ryb=b-ry(j)

        rzc=c-rz(j)

        txxi=txx+rmj*((ryb*ryb)+(rzc*rzc))

        ryb2=b-ry(j+1)

        rzc2=c-rz(j+1)

        txx=txxi+rmj2*((ryb2*ryb2)+(rzc2*rzc2))

        tyz=tyz-rmj*ryb*rzc-rmj2*ryb2*rzc2

        rxa=a-rx(j)

        rxa2=a-rx(j+1)

        txy=txy-(rmj*rxa)*ryb-(rmj2*rxa2)*ryb2

        txz=txz-(rmj*rxa)*rzc-(rmj2*rxa2)*rzc2

        tyy=tyy+rmj*((rxa*rxa)+(rzc*rzc))  &

             +rmj2*((rxa2*rxa2)+(rzc2*rzc2))

        tzz=tzz+rmj*((rxa*rxa)+(ryb*ryb))   &

             +rmj2*((rxa2*rxa2)+(ryb2*ryb2)) 

     END DO

  end if

  Tinert(1,1)=txx

  Tinert(1,2)=txy

  Tinert(1,3)=txz

  Tinert(2,1)=txy

  Tinert(2,2)=tyy

  Tinert(2,3)=tyz

  Tinert(3,1)=txz

  Tinert(3,2)=tyz

  Tinert(3,3)=tzz

  return

end subroutine tensor