newsys.f

Go to the documentation of this file.

      SUBROUTINE newsys(X,XD,M,NP,ENERGY,GAM)
*
*
*       Total energy of subsystem.
*       --------------------------
*
      IMPLICIT REAL*8  (a-h,o-z)
      REAL*8  m(*),x(*),xd(*)
*
*
*       Calculate total energy in physical units.
      t = 0.0d0
      v = 0.0d0
      DO 10 i = 1,np
          k1 = 3*(i - 1) + 1
          k2 = k1 + 1
          k3 = k2 + 1
          t = t + 0.5d0*m(i)*(xd(k1)**2 + xd(k2)**2 + xd(k3)**2)
          IF (i.EQ.np) go to 10
          j1 = i + 1
          DO 9 j = j1,np
              ki = 3*(i - 1)
              kj = 3*(j - 1)
              r2 = 0.0d0
              DO 8 k = 1,3
                  ki = ki + 1
                  kj = kj + 1
                  r2 = r2 + (x(ki) - x(kj))**2
    8         CONTINUE
              v = v - m(i)*m(j)/sqrt(r2)
    9     CONTINUE
   10 CONTINUE
*
      energy = t + v
      gam = t - v
*
      RETURN
*
      END