newreg.f

Go to the documentation of this file.

      SUBROUTINE newreg
*
*
*       Initialization of chain regularization.
*       ---------------------------------------
*
      IMPLICIT REAL*8  (a-h,o-z)
      REAL*8  m,mij
      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),cm(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
*
*
*       Set physical momenta unless switching case (with PI in COMMON).
      IF (.NOT.switch) THEN
          DO 4 i = 1,4
              DO 2 k = 1,3
                  ik = 3*(i - 1) + k
                  pi(ik) = m(i)*xd(ik)
    2         CONTINUE
    4     CONTINUE
      END IF
*
*       Determine vector chain for regularization.
      CALL status(x,i1,i2,i3,i4)
      ip1 = 3*(i1 - 1)
      ip2 = 3*(i2 - 1)
      ip3 = 3*(i3 - 1)
      ip4 = 3*(i4 - 1)
*
*       Form appropriate momenta & relative coordinates. 
      DO 5 k = 1,3
          w(k  ) = -pi(ip1+k)
          w(k+3) = -pi(ip1+k) - pi(ip2+k)
          w(k+6) = +pi(ip4+k)
          xr(k  ) = x(ip2+k) - x(ip1+k)
          xr(k+3) = x(ip3+k) - x(ip2+k)
          xr(k+6) = x(ip4+k) - x(ip3+k)
    5 CONTINUE
*
*       Perform KS transformations.
      DO 7 l = 1,3
          l1 = 3*(l - 1) + 1
          l2 = l1 + 1
          l3 = l2 + 1
          r2 = xr(l1)**2 + xr(l2)**2 + xr(l3)**2
          r(l) = sqrt(r2)
          lq1 = 4*(l - 1) + 1
          lq2 = lq1 + 1
          lq3 = lq2 + 1
          lq4 = lq3 + 1
          q(lq4) = 0.0d0
          a = r(l) + abs(xr(l1))
          q(lq1) = sqrt(.5d0*a)
          b = 1./(2.0d0*q(lq1))
          q(lq2) = xr(l2)*b
          q(lq3) = xr(l3)*b
          IF (xr(l1).LT.0.0d0) THEN
              u1 = q(lq1)
              q(lq1) = q(lq2)
              q(lq2) = u1
              u3 = q(lq3)
              q(lq3) = q(lq4)
              q(lq4) = u3
          END IF
          p(lq1) = 2.d0*(+q(lq1)*w(l1) + q(lq2)*w(l2) + q(lq3)*w(l3))
          p(lq2) = 2.d0*(-q(lq2)*w(l1) + q(lq1)*w(l2) + q(lq4)*w(l3))
          p(lq3) = 2.d0*(-q(lq3)*w(l1) - q(lq4)*w(l2) + q(lq1)*w(l3))
          p(lq4) = 2.d0*(+q(lq4)*w(l1) - q(lq3)*w(l2) + q(lq2)*w(l3))
    7 CONTINUE
*
*       Set mass factors (note scaling by 0.25 for DERQP4).
      ta(1) = 0.25d0*(.5d0/m(i1) + .5d0/m(i2))
      ta(2) = 0.25d0*(.5d0/m(i2) + .5d0/m(i3))
      ta(3) = 0.25d0*(.5d0/m(i3) + .5d0/m(i4))
      ta(4) = -0.25d0/m(i2)
      ta(5) = -0.25d0/m(i3)
      mij(1) = m(i1)*m(i2)
      mij(2) = m(i2)*m(i3)
      mij(3) = m(i3)*m(i4)
      mij(4) = m(i1)*m(i3)
      mij(5) = m(i2)*m(i4)
      mij(6) = m(i1)*m(i4)
*
      RETURN
*
      END