X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD_DFA%2Frattle.F;fp=source%2Funres%2Fsrc_MD_DFA%2Frattle.F;h=0000000000000000000000000000000000000000;hb=664d495e70d14eed4e97f7b8efd2e107dee2fd4e;hp=a2e503426790bcbdc9e3ca6dbac5b2d7107f36f1;hpb=77276d5043cefaf512451c3ad9f669ed22b90d04;p=unres.git

diff --git a/source/unres/src_MD_DFA/rattle.F b/source/unres/src_MD_DFA/rattle.F
deleted file mode 100644
index a2e5034..0000000
--- a/source/unres/src_MD_DFA/rattle.F
+++ /dev/null
@@ -1,706 +0,0 @@
-      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