endreg.f

Go to the documentation of this file.

      SUBROUTINE endreg
*
*
*       Transformation to physical variables.
*       -------------------------------------
*
      IMPLICIT REAL*8  (a-h,o-z)
      REAL*8  m,mij,cm(3)
      LOGICAL  switch,gtype,gtype0
      common/creg/  m(4),x(12),xd(12),p(12),q(12),time4,energy,epsr2,
     &              xr(9),w(9),r(6),ta(6),mij(6),cmx(10),rmax4,tmax,
     &              ds,tstep,eps,nstep4,name4(4),kz15,kz27,nreg,nfn
      common/tpr/   switch,gtype,gtype0
      common/savep/  pi(12)
      common/iconf/  i1,i2,i3,i4
*
*
      DO 1 l = 1,3
          l1 = 3*(l - 1) + 1
          l2 = l1 + 1
          l3 = l2 + 1
*
          k1 = 4*(l - 1) + 1
          k2 = k1 + 1
          k3 = k2 + 1
          k4 = k3 + 1
*
          xr(l1) = q(k1)**2 - q(k2)**2 - q(k3)**2 + q(k4)**2
          xr(l2) = 2.d0*(q(k1)*q(k2) - q(k3)*q(k4))
          xr(l3) = 2.d0*(q(k1)*q(k3) + q(k2)*q(k4))
          r(l) = q(k1)**2 + q(k2)**2 + q(k3)**2 + q(k4)**2
          a = 0.5d0/r(l)
          w(l1) = (q(k1)*p(k1) - q(k2)*p(k2) - q(k3)*p(k3) +
     &                                         q(k4)*p(k4))*a
          w(l2) = (q(k2)*p(k1) + q(k1)*p(k2) - q(k4)*p(k3) -
     &                                         q(k3)*p(k4))*a
          w(l3) = (q(k3)*p(k1) + q(k4)*p(k2) + q(k1)*p(k3) +
     &                                         q(k2)*p(k4))*a
    1 CONTINUE
*
      ip1 = i1*3 - 3
      ip2 = i2*3 - 3
      ip3 = i3*3 - 3
      ip4 = i4*3 - 3
*
      DO 2 k = 1,3
          pi(ip1+k) = -w(k  )
          pi(ip2+k) =  w(k  ) - w(k+3)
          pi(ip3+k) =  w(k+3) - w(k+6)
          pi(ip4+k) =  w(k+6)
          x(ip1+k) = -xr(k)
          x(ip2+k) = 0.0d0
*       Note terms X(IP1+K) = 0 & X(IP2+K) = XR(K) in early formulation.
          x(ip3+k) = x(ip2+k) + xr(k+3)
          x(ip4+k) = x(ip3+k) + xr(k+6)
    2 CONTINUE
*
*       Skip reduction to centre of mass after switching.
      IF (switch) RETURN
*
*       Obtain velocities from momenta and form c.m. coordinates.
      cm(1) = 0.0d0
      cm(2) = 0.0d0
      cm(3) = 0.0d0
      smass = m(1) + m(2) + m(3) + m(4)
      DO 4 i = 1,4
          li = 3*(i - 1)
          DO 3 k = 1,3
              xd(li+k) = pi(li+k)/m(i)
              cm(k) = cm(k) + m(i)*x(li+k)/smass
    3     CONTINUE
    4 CONTINUE
*
*       Express coordinates in the c.m. frame.
      l = 0
      DO 6 i = 1,4
          DO 5 k = 1,3
              l = l + 1
              x(l) = x(l) - cm(k)
    5     CONTINUE
    6 CONTINUE
*
      RETURN
*
      END