X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Frattle.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Frattle.F;h=5a8ed0c8871cefadad1c889d77d482539ab6422c;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-M-newcorr/rattle.F b/source/unres/src_MD-M-newcorr/rattle.F
new file mode 100644
index 0000000..5a8ed0c
--- /dev/null
+++ b/source/unres/src_MD-M-newcorr/rattle.F
@@ -0,0 +1,724 @@
+      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'
+#ifdef RATTLE
+      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
+#else
+      write (iout,*) 
+     & "RATTLE inactive; use -DRATTLE switch at compile time."
+      stop
+#endif
+      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'
+#ifdef RATTLE
+      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
+#else
+      write (iout,*) 
+     & "RATTLE inactive; use -DRATTLE option at compile time."
+      stop
+#endif
+      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'
+#ifdef RATTLE
+      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
+#else
+      write (iout,*) 
+     & "RATTLE inactive; use -DRATTLE option at compile time"
+      stop
+#endif
+      end