subroutine rattle1
c RATTLE algorithm for velocity Verlet - step 1, UNRES
c AL 9/24/04
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      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 gginv(maxres2,maxres2),
     & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
     & Cmat(MAXRES2,MAXRES2),x(MAXRES2),xcorr(3,MAXRES2)
      common /przechowalnia/ GGinv,gdc,Cmat,nbond
      integer max_rattle /5/
      logical lprn /.false./, lprn1 /.false./,not_done
      double precision tol_rattle /1.0d-5/
      if (lprn) write (iout,*) "RATTLE1"
      nbond=nct-nnt
      do i=nnt,nct
        if (itype(i).ne.10) nbond=nbond+1
      enddo
c Make a folded form of the Ginv-matrix
      ind=0
      ii=0
      do i=nnt,nct-1
        ii=ii+1
        do j=1,3
          ind=ind+1
          ind1=0
          jj=0
          do k=nnt,nct-1
            jj=jj+1
            do l=1,3 
              ind1=ind1+1
              if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
            enddo
          enddo
          do k=nnt,nct
            if (itype(k).ne.10) then
              jj=jj+1
              do l=1,3
                ind1=ind1+1
                if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
              enddo
            endif 
          enddo
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ii=ii+1
          do j=1,3
            ind=ind+1
            ind1=0
            jj=0
            do k=nnt,nct-1
              jj=jj+1
              do l=1,3 
                ind1=ind1+1
                if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
              enddo
            enddo
            do k=nnt,nct
              if (itype(k).ne.10) then
                jj=jj+1
                do l=1,3
                  ind1=ind1+1
                  if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
                enddo
              endif 
            enddo
          enddo
        endif
      enddo
      if (lprn1) then
        write (iout,*) "Matrix GGinv"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
      endif
      not_done=.true.
      iter=0
      do while (not_done)
      iter=iter+1
      if (iter.gt.max_rattle) then
        write (iout,*) "Error - too many iterations in RATTLE."
        stop
      endif
c Calculate the matrix C = GG**(-1) dC_old o dC
      ind1=0
      do i=nnt,nct-1
        ind1=ind1+1
        do j=1,3
          dC_uncor(j,ind1)=dC(j,i)
        enddo
      enddo 
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind1=ind1+1
          do j=1,3
            dC_uncor(j,ind1)=dC(j,i+nres)
          enddo
        endif
      enddo 
      do i=1,nbond
        ind=0
        do k=nnt,nct-1
          ind=ind+1
          do j=1,3
            gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
          enddo
        enddo
        do k=nnt,nct
          if (itype(k).ne.10) then
            ind=ind+1
            do j=1,3
              gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
            enddo
          endif
        enddo
      enddo
c Calculate deviations from standard virtual-bond lengths
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        x(ind)=vbld(i+1)**2-vbl**2
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
        endif
      enddo
      if (lprn) then
        write (iout,*) "Coordinates and violations"
        do i=1,nbond
          write(iout,'(i5,3f10.5,5x,e15.5)') 
     &     i,(dC_uncor(j,i),j=1,3),x(i)
        enddo
        write (iout,*) "Velocities and violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
     &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
          endif
        enddo
c        write (iout,*) "gdc"
c        do i=1,nbond
c          write (iout,*) "i",i
c          do j=1,nbond
c            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
c          enddo
c        enddo
      endif
      xmax=dabs(x(1))
      do i=2,nbond
        if (dabs(x(i)).gt.xmax) then
          xmax=dabs(x(i))
        endif
      enddo
      if (xmax.lt.tol_rattle) then
        not_done=.false.
        goto 100
      endif
c Calculate the matrix of the system of equations
      do i=1,nbond
        do j=1,nbond
          Cmat(i,j)=0.0d0
          do k=1,3
            Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
          enddo
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Matrix Cmat"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
      endif
      call gauss(Cmat,X,MAXRES2,nbond,1,*10) 
c Add constraint term to positions
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        do j=1,3
          xx=0.0d0
          do ii=1,nbond
            xx = xx+x(ii)*gdc(j,ind,ii)
          enddo
          xx=0.5d0*xx
          dC(j,i)=dC(j,i)-xx
          d_t_new(j,i)=d_t_new(j,i)-xx/d_time
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          do j=1,3
            xx=0.0d0
            do ii=1,nbond
              xx = xx+x(ii)*gdc(j,ind,ii)
            enddo
            xx=0.5d0*xx
            dC(j,i+nres)=dC(j,i+nres)-xx
            d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time 
          enddo
        endif
      enddo
c Rebuild the chain using the new coordinates
      call chainbuild_cart
      if (lprn) then
        write (iout,*) "New coordinates, Lagrange multipliers,",
     &  " and differences between actual and standard bond lengths"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          xx=vbld(i+1)**2-vbl**2
          write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &        i,(dC(j,i),j=1,3),x(ind),xx
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
            write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &       i,(dC(j,i+nres),j=1,3),x(ind),xx
          endif
        enddo
        write (iout,*) "Velocities and violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
     &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
          endif
        enddo
      endif
      enddo
  100 continue
      return
   10 write (iout,*) "Error - singularity in solving the system",
     & " of equations for Lagrange multipliers."
      stop
      end
c------------------------------------------------------------------------------
      subroutine rattle2
c RATTLE algorithm for velocity Verlet - step 2, UNRES
c AL 9/24/04
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      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 gginv(maxres2,maxres2),
     & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
     & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
      common /przechowalnia/ GGinv,gdc,Cmat,nbond
      integer max_rattle /5/
      logical lprn /.false./, lprn1 /.false./,not_done
      double precision tol_rattle /1.0d-5/
      if (lprn) write (iout,*) "RATTLE2"
      if (lprn) write (iout,*) "Velocity correction"
c Calculate the matrix G dC
      do i=1,nbond
        ind=0
        do k=nnt,nct-1
          ind=ind+1
          do j=1,3
            gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
          enddo
        enddo
        do k=nnt,nct
          if (itype(k).ne.10) then
            ind=ind+1
            do j=1,3
              gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
            enddo
          endif
        enddo
      enddo
c      if (lprn) then
c        write (iout,*) "gdc"
c        do i=1,nbond
c          write (iout,*) "i",i
c          do j=1,nbond
c            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
c          enddo
c        enddo
c      endif
c Calculate the matrix of the system of equations
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        do j=1,nbond
          Cmat(ind,j)=0.0d0
          do k=1,3
            Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
          enddo
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          do j=1,nbond
            Cmat(ind,j)=0.0d0
            do k=1,3
              Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
            enddo
          enddo
        endif
      enddo
c Calculate the scalar product dC o d_t_new
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        x(ind)=scalar(d_t(1,i),dC(1,i))
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
        endif
      enddo
      if (lprn) then
        write (iout,*) "Velocities and violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &     i,ind,(d_t(j,i),j=1,3),x(ind)
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &       i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
          endif
        enddo
      endif
      xmax=dabs(x(1))
      do i=2,nbond
        if (dabs(x(i)).gt.xmax) then
          xmax=dabs(x(i))
        endif
      enddo
      if (xmax.lt.tol_rattle) then
        not_done=.false.
        goto 100
      endif
      if (lprn1) then
        write (iout,*) "Matrix Cmat"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
      endif
      call gauss(Cmat,X,MAXRES2,nbond,1,*10) 
c Add constraint term to velocities
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        do j=1,3
          xx=0.0d0
          do ii=1,nbond
            xx = xx+x(ii)*gdc(j,ind,ii)
          enddo
          d_t(j,i)=d_t(j,i)-xx
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          do j=1,3
            xx=0.0d0
            do ii=1,nbond
              xx = xx+x(ii)*gdc(j,ind,ii)
            enddo
            d_t(j,i+nres)=d_t(j,i+nres)-xx
          enddo
        endif
      enddo
      if (lprn) then
        write (iout,*) 
     &    "New velocities, Lagrange multipliers violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)') 
     &       i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,2e15.5)') 
     &        i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),
     &        scalar(d_t(1,i+nres),dC(1,i+nres))
          endif
        enddo
      endif
  100 continue
      return
   10 write (iout,*) "Error - singularity in solving the system",
     & " of equations for Lagrange multipliers."
      stop
      end
c------------------------------------------------------------------------------
      subroutine rattle_brown
c RATTLE/LINCS algorithm for Brownian dynamics, UNRES
c AL 9/24/04
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      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 gginv(maxres2,maxres2),
     & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
     & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
      common /przechowalnia/ GGinv,gdc,Cmat,nbond
      integer max_rattle /5/
      logical lprn /.true./, lprn1 /.true./,not_done
      double precision tol_rattle /1.0d-5/
      if (lprn) write (iout,*) "RATTLE_BROWN"
      nbond=nct-nnt
      do i=nnt,nct
        if (itype(i).ne.10) nbond=nbond+1
      enddo
c Make a folded form of the Ginv-matrix
      ind=0
      ii=0
      do i=nnt,nct-1
        ii=ii+1
        do j=1,3
          ind=ind+1
          ind1=0
          jj=0
          do k=nnt,nct-1
            jj=jj+1
            do l=1,3 
              ind1=ind1+1
              if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
            enddo
          enddo
          do k=nnt,nct
            if (itype(k).ne.10) then
              jj=jj+1
              do l=1,3
                ind1=ind1+1
                if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
              enddo
            endif 
          enddo
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ii=ii+1
          do j=1,3
            ind=ind+1
            ind1=0
            jj=0
            do k=nnt,nct-1
              jj=jj+1
              do l=1,3 
                ind1=ind1+1
                if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
              enddo
            enddo
            do k=nnt,nct
              if (itype(k).ne.10) then
                jj=jj+1
                do l=1,3
                  ind1=ind1+1
                  if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
                enddo
              endif 
            enddo
          enddo
        endif
      enddo
      if (lprn1) then
        write (iout,*) "Matrix GGinv"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
      endif
      not_done=.true.
      iter=0
      do while (not_done)
      iter=iter+1
      if (iter.gt.max_rattle) then
        write (iout,*) "Error - too many iterations in RATTLE."
        stop
      endif
c Calculate the matrix C = GG**(-1) dC_old o dC
      ind1=0
      do i=nnt,nct-1
        ind1=ind1+1
        do j=1,3
          dC_uncor(j,ind1)=dC(j,i)
        enddo
      enddo 
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind1=ind1+1
          do j=1,3
            dC_uncor(j,ind1)=dC(j,i+nres)
          enddo
        endif
      enddo 
      do i=1,nbond
        ind=0
        do k=nnt,nct-1
          ind=ind+1
          do j=1,3
            gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
          enddo
        enddo
        do k=nnt,nct
          if (itype(k).ne.10) then
            ind=ind+1
            do j=1,3
              gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
            enddo
          endif
        enddo
      enddo
c Calculate deviations from standard virtual-bond lengths
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        x(ind)=vbld(i+1)**2-vbl**2
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
        endif
      enddo
      if (lprn) then
        write (iout,*) "Coordinates and violations"
        do i=1,nbond
          write(iout,'(i5,3f10.5,5x,e15.5)') 
     &     i,(dC_uncor(j,i),j=1,3),x(i)
        enddo
        write (iout,*) "Velocities and violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &     i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &       i+nres,ind,(d_t(j,i+nres),j=1,3),
     &       scalar(d_t(1,i+nres),dC_old(1,i+nres))
          endif
        enddo
        write (iout,*) "gdc"
        do i=1,nbond
          write (iout,*) "i",i
          do j=1,nbond
            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
          enddo
        enddo
      endif
      xmax=dabs(x(1))
      do i=2,nbond
        if (dabs(x(i)).gt.xmax) then
          xmax=dabs(x(i))
        endif
      enddo
      if (xmax.lt.tol_rattle) then
        not_done=.false.
        goto 100
      endif
c Calculate the matrix of the system of equations
      do i=1,nbond
        do j=1,nbond
          Cmat(i,j)=0.0d0
          do k=1,3
            Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
          enddo
        enddo
      enddo
      if (lprn1) then
        write (iout,*) "Matrix Cmat"
        call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
      endif
      call gauss(Cmat,X,MAXRES2,nbond,1,*10) 
c Add constraint term to positions
      ind=0
      do i=nnt,nct-1
        ind=ind+1
        do j=1,3
          xx=0.0d0
          do ii=1,nbond
            xx = xx+x(ii)*gdc(j,ind,ii)
          enddo
          xx=-0.5d0*xx
          d_t(j,i)=d_t(j,i)+xx/d_time
          dC(j,i)=dC(j,i)+xx
        enddo
      enddo
      do i=nnt,nct
        if (itype(i).ne.10) then
          ind=ind+1
          do j=1,3
            xx=0.0d0
            do ii=1,nbond
              xx = xx+x(ii)*gdc(j,ind,ii)
            enddo
            xx=-0.5d0*xx
            d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time 
            dC(j,i+nres)=dC(j,i+nres)+xx
          enddo
        endif
      enddo
c Rebuild the chain using the new coordinates
      call chainbuild_cart
      if (lprn) then
        write (iout,*) "New coordinates, Lagrange multipliers,",
     &  " and differences between actual and standard bond lengths"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          xx=vbld(i+1)**2-vbl**2
          write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &        i,(dC(j,i),j=1,3),x(ind),xx
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
            write (iout,'(i5,3f10.5,5x,f10.5,e15.5)') 
     &       i,(dC(j,i+nres),j=1,3),x(ind),xx
          endif
        enddo
        write (iout,*) "Velocities and violations"
        ind=0
        do i=nnt,nct-1
          ind=ind+1
          write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
        enddo
        do i=nnt,nct
          if (itype(i).ne.10) then
            ind=ind+1
            write (iout,'(2i5,3f10.5,5x,e15.5)') 
     &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
     &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
          endif
        enddo
      endif
      enddo
  100 continue
      return
   10 write (iout,*) "Error - singularity in solving the system",
     & " of equations for Lagrange multipliers."
      stop
      end