c ********************************************************************
      subroutine torsion_bd(n,ftor,phiold,iout) 
c
c  Torsional BD
c  update f,v by calculation dphi(z,i)
c  iout to determine if 1st or 2nd step in BD

      implicit double precision(a-h, o-z)
      parameter(maxa=500,maxn=maxa*3)
      integer n,iout
      double precision ftor(maxa),torque(maxa),rz(maxa)
     +	    , phiold(maxa)
      common/matrix/x(maxa,3),bangle(maxa),r(maxa,maxa)
      common/tor_param/dr,cg,phi0,dlk,phii,subseg
      common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2      
      common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
      common/bf_old/fold(maxa,3),vold(maxa,3),uold(maxa,3)
      common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)
       common/myrandom/amean,var(maxn),fvar(maxa)
 
      tmp = dr*dtime*cg

c calculate dphi(z,i)=DT*tmp+Q
c iout=0: first step in 2nd BD  set random ftor
      if (iout .eq. 0) then
         call rannv2(n,ftor,amean,fvar)  

         rz(1) = tmp*(phi(1)-phi(n)) + ftor(1)
         do i = 2, n
             rz(i)=tmp*(phi(i)-phi(i-1))+ftor(i)
        enddo 

c iout=1: second step in 2nd BD  using same ftor as the first step
c ftor is passed as variable
      else
         rz(1)=tmp*(phi(1)-phi(n)+phiold(1)-phiold(n))*0.5d0+ftor(1)
         do i = 2, n
         rz(i)=tmp*(phi(i)-phi(i-1)+phiold(i)-phi(i-1))*0.5d0+ftor(i)
         enddo
      endif

c  update vectors f and v (normalized)

      do i = 1, n
         ftotal = 0.d0
         vtotal = 0.d0
         do j = 1, 3
            tmp = f(i,j)
            f(i,j) = tmp + rz(i)*v(i,j)
            v(i,j) = v(i,j) - rz(i)*tmp
            ftotal = ftotal + f(i,j)*f(i,j)
            vtotal = vtotal + v(i,j)*v(i,j)
         enddo
         ftotal = dsqrt(ftotal)
         vtotal = dsqrt(vtotal)
         do j = 1, 3
            f(i,j) = f(i,j)/ftotal
            v(i,j) = v(i,j)/ftotal
         enddo 
      enddo

      return
      end