c bd_nohydro1.f

c 1st algorithm
c gc(n,3) is the forces
c rran(n,3) is the random variation
	subroutine bd_nohydro1(n,ntime)
	parameter(maxa=500,maxn=maxa*3)
	implicit double precision(a-h, o-z)
	integer n,n3
	double precision gc(maxn),rran(maxn)
      common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
      common/matrix/x(maxa,3),bangle(maxa)
      common/myrandom/amean,var(maxn),fvar(maxn)
       common/ischeme/icircular,iele,ihydro,iclose,iknot
     +,neuler
 
       common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
       common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)

        call force4(n,gc,ntime)
        n3=n*3
      if ((icircular.eq.1).and.(iclose.eq.1)) then
          iout = 0
	  call torsion_bd(n,ftor,phiold,iout)
      endif

	call rannv2(n3,rran,amean,var)

	ii=0
	do i=1,n
	   do j=1,3
	      ii=ii+1
           x(i,j)=x(i,j)+coef*gc(ii)+rran(ii)
c	  x(i,j)=x(i,j)+coef*gc(ii)
	 enddo
	enddo


	call update_rij(n)
        call update_vector(n)
        if ((icircular .eq. 1).and.(iclose.eq.1)) then
         call update_euler(n)
        endif

       	return
	end


c-----------------------------------------------------------
c 2nd algorithm
c gc(n,3) is the forces
c rran(n,3) is the random variation
	subroutine bd_nohydro2(n,ntime)
	parameter(maxa=500,maxn=maxa*3)
	implicit double precision(a-h, o-z)
	integer n,nout,ntime,n3
	double precision gc(maxn),rran(maxn),xold(maxa,3)
	double precision gc2(maxn),phiold(maxa),ftor(maxa,3)
	double precision fold(maxa,3),vold(maxa,3),uold(maxa,3)
      common/param/temp,bendc,r0,beadr,ss1,interval,dtime,gamma,
     +amass,am(maxn),coef,dtime2
      common/matrix/x(maxa,3),bangle(maxa)
      common/myrandom/amean,var(maxn),fvar(maxn)
      common/ischeme/icircular,iele,ihydro,iclose,iknot
     +,neuler
      common/bf_coords/f(maxa,3),v(maxa,3),u(maxa,3)
      common/euler/alpha(maxa),beta(maxa),gama(maxa),phi(maxa)
 
	n3=n*3
      do i = 1, n
         xold(i,1) = x(i,1)
         xold(i,2) = x(i,2)
         xold(i,3) = x(i,3)
         if (icircular .eq. 1) then
            do j = 1, 3
               fold(i,j) = f(i,j) 
               uold(i,j) = u(i,j) 
               vold(i,j) = v(i,j)
            enddo
            phiold(i) = phi(i)
         endif
       enddo
        call force4(n,gc,ntime)

      if ((icircular.eq.1).and.(iclose.eq.1)) then
          iout = 0
	  call torsion_bd(n,ftor,phiold,iout)
      endif


	call rannv2(n3,rran,amean,var)

	ii=0
	do i=1,n
	   do j=1,3
	      ii=ii+1
           x(i,j)=x(i,j)+coef*gc(ii)+rran(ii)
	  enddo
	enddo
	call update_rij(n)
        call update_vector(n)

      if ((icircular .eq. 1).and.(iclose.eq.1)) then
         call update_euler(n)
      endif


c xold is the auxiliary position vector now
c x=xnew
c gc2 and rran2 save the force and ran for coordinates xold
	
 
        call force4(n,gc2,ntime)
       if ((icircular .eq. 1).and.(iclose.eq.1)) then
         do i = 1, n
            do j = 1, 3
               f(i,j) = fold(i,j)
               v(i,j) = vold(i,j)
               u(i,j) = uold(i,j)
            enddo
         enddo
         iout = 1
         call torsion_bd(n,ftor,phiold,iout)
        endif

  	coef05=coef/2
	ii=0
	do i=1,n
	   do j=1,3
	      ii=ii+1
	      x(i,j)=xold(i,j)+coef05*(gc(ii)+gc2(ii))+rran(ii)
	   enddo
	enddo

	call update_rij(n)	
        call update_vector(n)

        if ((icircular .eq. 1).and.(iclose.eq.1)) then
         call update_euler(n)
        endif
    
	
c debug
c	write(*,*) 'config:',ntime
cc	do i=1,n
c	   write(*,*) x(i,1),x(i,2),x(i,3)
c	enddo

	return
	end