subroutine friction_0
     x  (idnode,mxnode,natms,imcon,tstep,engke,weight,cell,
     x  xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz,buffer,stress,fric)
      
c     
c***********************************************************************
c     
c     dl_poly subroutine for integrating newtonian equations of
c     motion in molecular dynamics - verlet leapfrog.
c     
c     parallel replicated data version
c     
c     copyright - daresbury laboratory 1992
c     author    - w. smith march 1992.
c     
c     $Author: thiels $
c     $Date: 2004-08-19 13:51:33 +0200 (Thu, 19. Aug 2004) $
c     $Revision: 2490 $
c     $State$
c     
c     FRICTION DYNAMICS
c     Propagate with a constant friction proportional to the velocity.
c     The friction parameter fric is defined as
c        fric = friction coefficient * Dt / 2
c
c     Hans Martin Senn (hms), MPI fuer Kohlenforschung, Nov. 2002
c***********************************************************************
c     
      
      implicit real*8 (a-h,o-z)
      include 'dl_params.inc'
      
      dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
      dimension vxx(mxatms),vyy(mxatms),vzz(mxatms)
      dimension fxx(mxatms),fyy(mxatms),fzz(mxatms)
      dimension weight(mxatms),buffer(mxbuff),cell(9)
      dimension stress(9),stres1(9)
c     
c     set kinetic energy accumulator
      
      engke=0.d0
#ifdef STRESS
c     
c     temporary stress tensor accumulators
      do i = 1,9
        stres1(i) = 0.d0
      enddo
#endif
c     
c     block indices
      
      iatm0 = (idnode*natms)/mxnode + 1
      iatm1 = ((idnode+1)*natms)/mxnode
      
c     
c     move atoms by leapfrog algorithm
      svar = 1.d0/(1.d0 + fric)

      do i=iatm0,iatm1
c     
c     store old velocities
        
        uxx=vxx(i)
        uyy=vyy(i)
        uzz=vzz(i)
c     
c     calculate new velocities
        svar2 = tstep/weight(i)
        vxx(i)=svar * ((1.d0 - fric)*vxx(i) + svar2*fxx(i))
        vyy(i)=svar * ((1.d0 - fric)*vyy(i) + svar2*fyy(i))
        vzz(i)=svar * ((1.d0 - fric)*vzz(i) + svar2*fzz(i))
c     
c     update positions
        
        xxx(i)=xxx(i)+tstep*vxx(i)
        yyy(i)=yyy(i)+tstep*vyy(i)
        zzz(i)=zzz(i)+tstep*vzz(i)
c
c     current velocity
        vxt = 0.5d0*(vxx(i)+uxx)
        vyt = 0.5d0*(vyy(i)+uyy)
        vzt = 0.5d0*(vzz(i)+uzz)
c     
c     calculate kinetic energy
        engke=engke+0.5d0*weight(i)*(vxt*vxt+vyt*vyt+vzt*vzt)
        
#ifdef STRESS
c     
c     kinetic contribution to stress tensor
        
        stres1(1) = stres1(1) + weight(i)*vxt*vxt
        stres1(2) = stres1(2) + weight(i)*vxt*vyt
        stres1(3) = stres1(3) + weight(i)*vxt*vzt
        stres1(5) = stres1(5) + weight(i)*vyt*vyt
        stres1(6) = stres1(6) + weight(i)*vyt*vzt
        stres1(9) = stres1(9) + weight(i)*vzt*vzt
#endif
      enddo
      
c     
c     periodic boundary condition
      
      call images(imcon,idnode,mxnode,natms,cell,xxx,yyy,zzz)
      
c     
c     global exchange of configuration data
      
      if(mxnode.gt.1)then
        
        nbuff=mxbuff
        call gdsum(engke,1,buffer)
        call merge(idnode,mxnode,natms,nbuff,xxx,yyy,zzz,buffer)
        call merge(idnode,mxnode,natms,nbuff,vxx,vyy,vzz,buffer)
        
      endif
#ifdef STRESS
c     
c     complete stress tensor
      
      if(mxnode.gt.1) call gdsum(stres1,9,buffer)
      
      stress(1) = stress(1) + stres1(1)
      stress(2) = stress(2) + stres1(2)
      stress(3) = stress(3) + stres1(3)
      stress(4) = stress(4) + stres1(2)
      stress(5) = stress(5) + stres1(5)
      stress(6) = stress(6) + stres1(6)
      stress(7) = stress(7) + stres1(3)
      stress(8) = stress(8) + stres1(6)
      stress(9) = stress(9) + stres1(9)
#endif
      return
      end