c-------------------------------------------------------------------------------
      subroutine brown_step(itime)
c------------------------------------------------
c  Perform a single Euler integration step of Brownian dynamics
c------------------------------------------------
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
#ifdef MPI
      include 'mpif.h'
#endif
      include 'COMMON.CONTROL'
      include 'COMMON.VAR'
      include 'COMMON.MD'
#ifndef LANG0
      include 'COMMON.LANGEVIN'
#else
      include 'COMMON.LANGEVIN.lang0'
#endif
      include 'COMMON.CHAIN'
      include 'COMMON.DERIV'
      include 'COMMON.GEO'
      include 'COMMON.LOCAL'
      include 'COMMON.INTERACT'
      include 'COMMON.IOUNITS'
      include 'COMMON.NAMES'
      include 'COMMON.TIME1'
      double precision zapas(MAXRES6)
      integer ilen,rstcount
      external ilen
      double precision stochforcvec(MAXRES6)
      double precision Bmat(MAXRES6,MAXRES2),Cmat(maxres2,maxres2),
     & Cinv(maxres2,maxres2),GBmat(MAXRES6,MAXRES2),
     & Tmat(MAXRES6,MAXRES2),Pmat(maxres6,maxres6),Td(maxres6),
     & ppvec(maxres2)
      common /stochcalc/ stochforcvec
      common /gucio/ cm
      integer itime
      logical lprn /.false./,lprn1 /.false./
      integer maxiter /5/
      double precision difftol /1.0d-5/
      nbond=nct-nnt
      do i=nnt,nct
        if (itype(i).ne.10) nbond=nbond+1
      enddo
c
      if (lprn1) then
        write (iout,*) "Generalized inverse of fricmat"
        call matout(dimen,dimen,MAXRES6,MAXRES6,fricmat)
      endif 
      do i=1,dimen
        do j=1,nbond
          Bmat(i,j)=0.0d0
        enddo
      enddo
      ind=3
      ind1=0
      do i=nnt,nct-1
        ind1=ind1+1
        do j=1,3
          Bmat(ind+j,ind1)=dC_norm(j,i)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind1=ind1+1
          do j=1,3
            Bmat(ind+j,ind1)=dC_norm(j,i+nres)
          enddo
          ind=ind+3
        endif
      enddo
      if (lprn1) then 
        write (iout,*) "Matrix Bmat"
        call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Bmat)
      endif
      do i=1,dimen
        do j=1,nbond
          GBmat(i,j)=0.0d0
          do k=1,dimen
            GBmat(i,j)=GBmat(i,j)+fricmat(i,k)*Bmat(k,j)
          enddo
        enddo
      enddo   
      if (lprn1) then
        write (iout,*) "Matrix GBmat"
        call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Gbmat)
      endif
      do i=1,nbond
        do j=1,nbond
          Cmat(i,j)=0.0d0
          do k=1,dimen
            Cmat(i,j)=Cmat(i,j)+Bmat(k,i)*GBmat(k,j)
          enddo
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Matrix Cmat"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
      endif
      call matinvert(nbond,MAXRES2,Cmat,Cinv) 
      if (lprn1) then
        write (iout,*) "Matrix Cinv"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cinv)
      endif
      do i=1,dimen
        do j=1,nbond
          Tmat(i,j)=0.0d0
          do k=1,nbond
            Tmat(i,j)=Tmat(i,j)+GBmat(i,k)*Cinv(k,j)
          enddo
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Matrix Tmat"
        call MATOUT(nbond,dimen,MAXRES6,MAXRES2,Tmat)
      endif
      do i=1,dimen
        do j=1,dimen
          if (i.eq.j) then
            Pmat(i,j)=1.0d0
          else
            Pmat(i,j)=0.0d0
          endif
          do k=1,nbond
            Pmat(i,j)=Pmat(i,j)-Tmat(i,k)*Bmat(j,k)
          enddo
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Matrix Pmat"
        call MATOUT(dimen,dimen,MAXRES6,MAXRES6,Pmat)
      endif
      do i=1,dimen
        Td(i)=0.0d0
        ind=0
        do k=nnt,nct-1
          ind=ind+1
          Td(i)=Td(i)+vbl*Tmat(i,ind)
        enddo
        do k=nnt,nct
          if (itype(k).ne.10) then
            ind=ind+1
            Td(i)=Td(i)+vbldsc0(1,itype(k))*Tmat(i,ind)
          endif
        enddo
      enddo 
      if (lprn1) then
        write (iout,*) "Vector Td"
        do i=1,dimen
          write (iout,'(i5,f10.5)') i,Td(i)
        enddo
      endif
      call stochastic_force(stochforcvec)
      if (lprn) then
        write (iout,*) "stochforcvec"
        do i=1,dimen
          write (iout,*) i,stochforcvec(i)
        enddo
      endif
      do j=1,3
        zapas(j)=-gcart(j,0)+stochforcvec(j)
        d_t_work(j)=d_t(j,0)
        dC_work(j)=dC_old(j,0)
      enddo
      ind=3      
      do i=nnt,nct-1
        do j=1,3
          ind=ind+1
          zapas(ind)=-gcart(j,i)+stochforcvec(ind)
          dC_work(ind)=dC_old(j,i)
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            ind=ind+1
            zapas(ind)=-gxcart(j,i)+stochforcvec(ind)
            dC_work(ind)=dC_old(j,i+nres)
          enddo
        endif
      enddo

      if (lprn) then
        write (iout,*) "Initial d_t_work"
        do i=1,dimen
          write (iout,*) i,d_t_work(i)
        enddo
      endif

      do i=1,dimen
        d_t_work(i)=0.0d0
        do j=1,dimen
          d_t_work(i)=d_t_work(i)+fricmat(i,j)*zapas(j)
        enddo
      enddo

      do i=1,dimen
        zapas(i)=Td(i)
        do j=1,dimen
          zapas(i)=zapas(i)+Pmat(i,j)*(dC_work(j)+d_t_work(j)*d_time)
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Final d_t_work and zapas"
        do i=1,dimen
          write (iout,*) i,d_t_work(i),zapas(i)
        enddo
      endif

      do j=1,3
        d_t(j,0)=d_t_work(j)
        dc(j,0)=zapas(j)
        dc_work(j)=dc(j,0)
      enddo
      ind=3
      do i=nnt,nct-1
        do j=1,3
          d_t(j,i)=d_t_work(i)
          dc(j,i)=zapas(ind+j)
          dc_work(ind+j)=dc(j,i)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        do j=1,3
          d_t(j,i+nres)=d_t_work(ind+j)
          dc(j,i+nres)=zapas(ind+j)
          dc_work(ind+j)=dc(j,i+nres)
        enddo
        ind=ind+3
      enddo
      if (lprn) then
        call chainbuild_cart
        write (iout,*) "Before correction for rotational lengthening"
        write (iout,*) "New coordinates",
     &  " and differences between actual and standard bond lengths"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          xx=vbld(i+1)-vbl
          write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &        i,(dC(j,i),j=1,3),xx
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            xx=vbld(i+nres)-vbldsc0(1,itype(i))
            write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &       i,(dC(j,i+nres),j=1,3),xx
          endif
        enddo
      endif
c Second correction (rotational lengthening)
c      do iter=1,maxiter
      diffmax=0.0d0
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        blen2 = scalar(dc(1,i),dc(1,i))
        ppvec(ind)=2*vbl**2-blen2
        diffbond=dabs(vbl-dsqrt(blen2))
        if (diffbond.gt.diffmax) diffmax=diffbond
        if (ppvec(ind).gt.0.0d0) then
          ppvec(ind)=dsqrt(ppvec(ind))
        else
          ppvec(ind)=0.0d0
        endif
        if (lprn) then
          write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
        endif
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          blen2 = scalar(dc(1,i+nres),dc(1,i+nres))
          ppvec(ind)=2*vbldsc0(1,itype(i))**2-blen2
          diffbond=dabs(vbldsc0(1,itype(i))-dsqrt(blen2))
          if (diffbond.gt.diffmax) diffmax=diffbond
          if (ppvec(ind).gt.0.0d0) then
            ppvec(ind)=dsqrt(ppvec(ind))
          else
            ppvec(ind)=0.0d0
          endif
          if (lprn) then
            write (iout,'(i5,3f10.5)') ind,diffbond,ppvec(ind)
          endif
        endif
      enddo
      if (lprn) write (iout,*) "iter",iter," diffmax",diffmax
      if (diffmax.lt.difftol) goto 10
      do i=1,dimen
        Td(i)=0.0d0
        do j=1,nbond
          Td(i)=Td(i)+ppvec(j)*Tmat(i,j)
        enddo
      enddo 
      do i=1,dimen
        zapas(i)=Td(i)
        do j=1,dimen
          zapas(i)=zapas(i)+Pmat(i,j)*dc_work(j)
        enddo
      enddo
      do j=1,3
        dc(j,0)=zapas(j)
        dc_work(j)=zapas(j)
      enddo
      ind=3
      do i=nnt,nct-1
        do j=1,3
          dc(j,i)=zapas(ind+j)
          dc_work(ind+j)=zapas(ind+j)
        enddo
        ind=ind+3
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          do j=1,3
            dc(j,i+nres)=zapas(ind+j)
            dc_work(ind+j)=zapas(ind+j)
          enddo
          ind=ind+3
        endif
      enddo 
c   Building the chain from the newly calculated coordinates	
      call chainbuild_cart
      if(ntwe.ne.0) then
      if (large.and. mod(itime,ntwe).eq.0) then
        write (iout,*) "Cartesian and internal coordinates: step 1"
        call cartprint
        call intout
        write (iout,'(a)') "Potential forces"
        do i=0,nres
      	  write (iout,'(i3,3f10.5,3x,3f10.5)') i,(-gcart(j,i),j=1,3),
     &    (-gxcart(j,i),j=1,3)
        enddo
        write (iout,'(a)') "Stochastic forces"
        do i=0,nres
      	  write (iout,'(i3,3f10.5,3x,3f10.5)') i,(stochforc(j,i),j=1,3),
     &    (stochforc(j,i+nres),j=1,3)
        enddo
        write (iout,'(a)') "Velocities"
        do i=0,nres
      	  write (iout,'(i3,3f10.5,3x,3f10.5)') i,(d_t(j,i),j=1,3),
     &    (d_t(j,i+nres),j=1,3)
        enddo
      endif
      endif
      if (lprn) then
        write (iout,*) "After correction for rotational lengthening"
        write (iout,*) "New coordinates",
     &  " and differences between actual and standard bond lengths"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          xx=vbld(i+1)-vbl
          write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &        i,(dC(j,i),j=1,3),xx
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            xx=vbld(i+nres)-vbldsc0(1,itype(i))
            write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &       i,(dC(j,i+nres),j=1,3),xx
          endif
        enddo
      endif
c      ENDDO
c      write (iout,*) "Too many attempts at correcting the bonds"
c      stop
   10 continue
#ifdef MPI
      tt0 =MPI_Wtime()
#else
      tt0 = tcpu()
#endif
c Calculate energy and forces
      call zerograd
      call etotal(potEcomp)
      potE=potEcomp(0)-potEcomp(20)
      call cartgrad
      totT=totT+d_time
c  Calculate the kinetic and total energy and the kinetic temperature
      call kinetic(EK)
#ifdef MPI
      t_enegrad=t_enegrad+MPI_Wtime()-tt0
#else
      t_enegrad=t_enegrad+tcpu()-tt0
#endif
      totE=EK+potE
      kinetic_T=2.0d0/(dimen*Rb)*EK
      return
      end