source/unres/src_MD/rattle.F

   1       subroutine rattle1
   2 c RATTLE algorithm for velocity Verlet - step 1, UNRES
   3 c AL 9/24/04
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.VAR'
   8       include 'COMMON.MD'
   9 #ifndef LANG0
  10       include 'COMMON.LANGEVIN'
  11 #else
  12       include 'COMMON.LANGEVIN.lang0'
  13 #endif
  14       include 'COMMON.CHAIN'
  15       include 'COMMON.DERIV'
  16       include 'COMMON.GEO'
  17       include 'COMMON.LOCAL'
  18       include 'COMMON.INTERACT'
  19       include 'COMMON.IOUNITS'
  20       include 'COMMON.NAMES'
  21       include 'COMMON.TIME1'
  22       double precision gginv(maxres2,maxres2),
  23      & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
  24      & Cmat(MAXRES2,MAXRES2),x(MAXRES2),xcorr(3,MAXRES2)
  25       common /przechowalnia/ GGinv,gdc,Cmat,nbond
  26       integer max_rattle /5/
  27       logical lprn /.false./, lprn1 /.false./,not_done
  28       double precision tol_rattle /1.0d-5/
  29       if (lprn) write (iout,*) "RATTLE1"
  30       nbond=nct-nnt
  31       do i=nnt,nct
  32         if (itype(i).ne.10) nbond=nbond+1
  33       enddo
  34 c Make a folded form of the Ginv-matrix
  35       ind=0
  36       ii=0
  37       do i=nnt,nct-1
  38         ii=ii+1
  39         do j=1,3
  40           ind=ind+1
  41           ind1=0
  42           jj=0
  43           do k=nnt,nct-1
  44             jj=jj+1
  45             do l=1,3
  46               ind1=ind1+1
  47               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
  48             enddo
  49           enddo
  50           do k=nnt,nct
  51             if (itype(k).ne.10) then
  52               jj=jj+1
  53               do l=1,3
  54                 ind1=ind1+1
  55                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
  56               enddo
  57             endif
  58           enddo
  59         enddo
  60       enddo
  61       do i=nnt,nct
  62         if (itype(i).ne.10) then
  63           ii=ii+1
  64           do j=1,3
  65             ind=ind+1
  66             ind1=0
  67             jj=0
  68             do k=nnt,nct-1
  69               jj=jj+1
  70               do l=1,3
  71                 ind1=ind1+1
  72                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
  73               enddo
  74             enddo
  75             do k=nnt,nct
  76               if (itype(k).ne.10) then
  77                 jj=jj+1
  78                 do l=1,3
  79                   ind1=ind1+1
  80                   if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
  81                 enddo
  82               endif
  83             enddo
  84           enddo
  85         endif
  86       enddo
  87       if (lprn1) then
  88         write (iout,*) "Matrix GGinv"
  89         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
  90       endif
  91       not_done=.true.
  92       iter=0
  93       do while (not_done)
  94       iter=iter+1
  95       if (iter.gt.max_rattle) then
  96         write (iout,*) "Error - too many iterations in RATTLE."
  97         stop
  98       endif
  99 c Calculate the matrix C = GG**(-1) dC_old o dC
 100       ind1=0
 101       do i=nnt,nct-1
 102         ind1=ind1+1
 103         do j=1,3
 104           dC_uncor(j,ind1)=dC(j,i)
 105         enddo
 106       enddo
 107       do i=nnt,nct
 108         if (itype(i).ne.10) then
 109           ind1=ind1+1
 110           do j=1,3
 111             dC_uncor(j,ind1)=dC(j,i+nres)
 112           enddo
 113         endif
 114       enddo
 115       do i=1,nbond
 116         ind=0
 117         do k=nnt,nct-1
 118           ind=ind+1
 119           do j=1,3
 120             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
 121           enddo
 122         enddo
 123         do k=nnt,nct
 124           if (itype(k).ne.10) then
 125             ind=ind+1
 126             do j=1,3
 127               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
 128             enddo
 129           endif
 130         enddo
 131       enddo
 132 c Calculate deviations from standard virtual-bond lengths
 133       ind=0
 134       do i=nnt,nct-1
 135         ind=ind+1
 136         x(ind)=vbld(i+1)**2-vbl**2
 137       enddo
 138       do i=nnt,nct
 139         if (itype(i).ne.10) then
 140           ind=ind+1
 141           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
 142         endif
 143       enddo
 144       if (lprn) then
 145         write (iout,*) "Coordinates and violations"
 146         do i=1,nbond
 147           write(iout,'(i5,3f10.5,5x,e15.5)')
 148      &     i,(dC_uncor(j,i),j=1,3),x(i)
 149         enddo
 150         write (iout,*) "Velocities and violations"
 151         ind=0
 152         do i=nnt,nct-1
 153           ind=ind+1
 154           write (iout,'(2i5,3f10.5,5x,e15.5)')
 155      &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
 156         enddo
 157         do i=nnt,nct
 158           if (itype(i).ne.10) then
 159             ind=ind+1
 160             write (iout,'(2i5,3f10.5,5x,e15.5)')
 161      &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
 162      &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
 163           endif
 164         enddo
 165 c        write (iout,*) "gdc"
 166 c        do i=1,nbond
 167 c          write (iout,*) "i",i
 168 c          do j=1,nbond
 169 c            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
 170 c          enddo
 171 c        enddo
 172       endif
 173       xmax=dabs(x(1))
 174       do i=2,nbond
 175         if (dabs(x(i)).gt.xmax) then
 176           xmax=dabs(x(i))
 177         endif
 178       enddo
 179       if (xmax.lt.tol_rattle) then
 180         not_done=.false.
 181         goto 100
 182       endif
 183 c Calculate the matrix of the system of equations
 184       do i=1,nbond
 185         do j=1,nbond
 186           Cmat(i,j)=0.0d0
 187           do k=1,3
 188             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
 189           enddo
 190         enddo
 191       enddo
 192       if (lprn1) then
 193         write (iout,*) "Matrix Cmat"
 194         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
 195       endif
 196       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
 197 c Add constraint term to positions
 198       ind=0
 199       do i=nnt,nct-1
 200         ind=ind+1
 201         do j=1,3
 202           xx=0.0d0
 203           do ii=1,nbond
 204             xx = xx+x(ii)*gdc(j,ind,ii)
 205           enddo
 206           xx=0.5d0*xx
 207           dC(j,i)=dC(j,i)-xx
 208           d_t_new(j,i)=d_t_new(j,i)-xx/d_time
 209         enddo
 210       enddo
 211       do i=nnt,nct
 212         if (itype(i).ne.10) then
 213           ind=ind+1
 214           do j=1,3
 215             xx=0.0d0
 216             do ii=1,nbond
 217               xx = xx+x(ii)*gdc(j,ind,ii)
 218             enddo
 219             xx=0.5d0*xx
 220             dC(j,i+nres)=dC(j,i+nres)-xx
 221             d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
 222           enddo
 223         endif
 224       enddo
 225 c Rebuild the chain using the new coordinates
 226       call chainbuild_cart
 227       if (lprn) then
 228         write (iout,*) "New coordinates, Lagrange multipliers,",
 229      &  " and differences between actual and standard bond lengths"
 230         ind=0
 231         do i=nnt,nct-1
 232           ind=ind+1
 233           xx=vbld(i+1)**2-vbl**2
 234           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 235      &        i,(dC(j,i),j=1,3),x(ind),xx
 236         enddo
 237         do i=nnt,nct
 238           if (itype(i).ne.10) then
 239             ind=ind+1
 240             xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
 241             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 242      &       i,(dC(j,i+nres),j=1,3),x(ind),xx
 243           endif
 244         enddo
 245         write (iout,*) "Velocities and violations"
 246         ind=0
 247         do i=nnt,nct-1
 248           ind=ind+1
 249           write (iout,'(2i5,3f10.5,5x,e15.5)')
 250      &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
 251         enddo
 252         do i=nnt,nct
 253           if (itype(i).ne.10) then
 254             ind=ind+1
 255             write (iout,'(2i5,3f10.5,5x,e15.5)')
 256      &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
 257      &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
 258           endif
 259         enddo
 260       endif
 261       enddo
 262   100 continue
 263       return
 264    10 write (iout,*) "Error - singularity in solving the system",
 265      & " of equations for Lagrange multipliers."
 266       stop
 267       end
 268 c------------------------------------------------------------------------------
 269       subroutine rattle2
 270 c RATTLE algorithm for velocity Verlet - step 2, UNRES
 271 c AL 9/24/04
 272       implicit real*8 (a-h,o-z)
 273       include 'DIMENSIONS'
 274       include 'COMMON.CONTROL'
 275       include 'COMMON.VAR'
 276       include 'COMMON.MD'
 277 #ifndef LANG0
 278       include 'COMMON.LANGEVIN'
 279 #else
 280       include 'COMMON.LANGEVIN.lang0'
 281 #endif
 282       include 'COMMON.CHAIN'
 283       include 'COMMON.DERIV'
 284       include 'COMMON.GEO'
 285       include 'COMMON.LOCAL'
 286       include 'COMMON.INTERACT'
 287       include 'COMMON.IOUNITS'
 288       include 'COMMON.NAMES'
 289       include 'COMMON.TIME1'
 290       double precision gginv(maxres2,maxres2),
 291      & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
 292      & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
 293       common /przechowalnia/ GGinv,gdc,Cmat,nbond
 294       integer max_rattle /5/
 295       logical lprn /.false./, lprn1 /.false./,not_done
 296       double precision tol_rattle /1.0d-5/
 297       if (lprn) write (iout,*) "RATTLE2"
 298       if (lprn) write (iout,*) "Velocity correction"
 299 c Calculate the matrix G dC
 300       do i=1,nbond
 301         ind=0
 302         do k=nnt,nct-1
 303           ind=ind+1
 304           do j=1,3
 305             gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
 306           enddo
 307         enddo
 308         do k=nnt,nct
 309           if (itype(k).ne.10) then
 310             ind=ind+1
 311             do j=1,3
 312               gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
 313             enddo
 314           endif
 315         enddo
 316       enddo
 317 c      if (lprn) then
 318 c        write (iout,*) "gdc"
 319 c        do i=1,nbond
 320 c          write (iout,*) "i",i
 321 c          do j=1,nbond
 322 c            write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
 323 c          enddo
 324 c        enddo
 325 c      endif
 326 c Calculate the matrix of the system of equations
 327       ind=0
 328       do i=nnt,nct-1
 329         ind=ind+1
 330         do j=1,nbond
 331           Cmat(ind,j)=0.0d0
 332           do k=1,3
 333             Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
 334           enddo
 335         enddo
 336       enddo
 337       do i=nnt,nct
 338         if (itype(i).ne.10) then
 339           ind=ind+1
 340           do j=1,nbond
 341             Cmat(ind,j)=0.0d0
 342             do k=1,3
 343               Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
 344             enddo
 345           enddo
 346         endif
 347       enddo
 348 c Calculate the scalar product dC o d_t_new
 349       ind=0
 350       do i=nnt,nct-1
 351         ind=ind+1
 352         x(ind)=scalar(d_t(1,i),dC(1,i))
 353       enddo
 354       do i=nnt,nct
 355         if (itype(i).ne.10) then
 356           ind=ind+1
 357           x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
 358         endif
 359       enddo
 360       if (lprn) then
 361         write (iout,*) "Velocities and violations"
 362         ind=0
 363         do i=nnt,nct-1
 364           ind=ind+1
 365           write (iout,'(2i5,3f10.5,5x,e15.5)')
 366      &     i,ind,(d_t(j,i),j=1,3),x(ind)
 367         enddo
 368         do i=nnt,nct
 369           if (itype(i).ne.10) then
 370             ind=ind+1
 371             write (iout,'(2i5,3f10.5,5x,e15.5)')
 372      &       i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
 373           endif
 374         enddo
 375       endif
 376       xmax=dabs(x(1))
 377       do i=2,nbond
 378         if (dabs(x(i)).gt.xmax) then
 379           xmax=dabs(x(i))
 380         endif
 381       enddo
 382       if (xmax.lt.tol_rattle) then
 383         not_done=.false.
 384         goto 100
 385       endif
 386       if (lprn1) then
 387         write (iout,*) "Matrix Cmat"
 388         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
 389       endif
 390       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
 391 c Add constraint term to velocities
 392       ind=0
 393       do i=nnt,nct-1
 394         ind=ind+1
 395         do j=1,3
 396           xx=0.0d0
 397           do ii=1,nbond
 398             xx = xx+x(ii)*gdc(j,ind,ii)
 399           enddo
 400           d_t(j,i)=d_t(j,i)-xx
 401         enddo
 402       enddo
 403       do i=nnt,nct
 404         if (itype(i).ne.10) then
 405           ind=ind+1
 406           do j=1,3
 407             xx=0.0d0
 408             do ii=1,nbond
 409               xx = xx+x(ii)*gdc(j,ind,ii)
 410             enddo
 411             d_t(j,i+nres)=d_t(j,i+nres)-xx
 412           enddo
 413         endif
 414       enddo
 415       if (lprn) then
 416         write (iout,*)
 417      &    "New velocities, Lagrange multipliers violations"
 418         ind=0
 419         do i=nnt,nct-1
 420           ind=ind+1
 421           if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)')
 422      &       i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
 423         enddo
 424         do i=nnt,nct
 425           if (itype(i).ne.10) then
 426             ind=ind+1
 427             write (iout,'(2i5,3f10.5,5x,2e15.5)')
 428      &        i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),
 429      &        scalar(d_t(1,i+nres),dC(1,i+nres))
 430           endif
 431         enddo
 432       endif
 433   100 continue
 434       return
 435    10 write (iout,*) "Error - singularity in solving the system",
 436      & " of equations for Lagrange multipliers."
 437       stop
 438       end
 439 c------------------------------------------------------------------------------
 440       subroutine rattle_brown
 441 c RATTLE/LINCS algorithm for Brownian dynamics, UNRES
 442 c AL 9/24/04
 443       implicit real*8 (a-h,o-z)
 444       include 'DIMENSIONS'
 445       include 'COMMON.CONTROL'
 446       include 'COMMON.VAR'
 447       include 'COMMON.MD'
 448 #ifndef LANG0
 449       include 'COMMON.LANGEVIN'
 450 #else
 451       include 'COMMON.LANGEVIN.lang0'
 452 #endif
 453       include 'COMMON.CHAIN'
 454       include 'COMMON.DERIV'
 455       include 'COMMON.GEO'
 456       include 'COMMON.LOCAL'
 457       include 'COMMON.INTERACT'
 458       include 'COMMON.IOUNITS'
 459       include 'COMMON.NAMES'
 460       include 'COMMON.TIME1'
 461       double precision gginv(maxres2,maxres2),
 462      & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
 463      & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
 464       common /przechowalnia/ GGinv,gdc,Cmat,nbond
 465       integer max_rattle /5/
 466       logical lprn /.true./, lprn1 /.true./,not_done
 467       double precision tol_rattle /1.0d-5/
 468       if (lprn) write (iout,*) "RATTLE_BROWN"
 469       nbond=nct-nnt
 470       do i=nnt,nct
 471         if (itype(i).ne.10) nbond=nbond+1
 472       enddo
 473 c Make a folded form of the Ginv-matrix
 474       ind=0
 475       ii=0
 476       do i=nnt,nct-1
 477         ii=ii+1
 478         do j=1,3
 479           ind=ind+1
 480           ind1=0
 481           jj=0
 482           do k=nnt,nct-1
 483             jj=jj+1
 484             do l=1,3
 485               ind1=ind1+1
 486               if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
 487             enddo
 488           enddo
 489           do k=nnt,nct
 490             if (itype(k).ne.10) then
 491               jj=jj+1
 492               do l=1,3
 493                 ind1=ind1+1
 494                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
 495               enddo
 496             endif
 497           enddo
 498         enddo
 499       enddo
 500       do i=nnt,nct
 501         if (itype(i).ne.10) then
 502           ii=ii+1
 503           do j=1,3
 504             ind=ind+1
 505             ind1=0
 506             jj=0
 507             do k=nnt,nct-1
 508               jj=jj+1
 509               do l=1,3
 510                 ind1=ind1+1
 511                 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
 512               enddo
 513             enddo
 514             do k=nnt,nct
 515               if (itype(k).ne.10) then
 516                 jj=jj+1
 517                 do l=1,3
 518                   ind1=ind1+1
 519                   if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
 520                 enddo
 521               endif
 522             enddo
 523           enddo
 524         endif
 525       enddo
 526       if (lprn1) then
 527         write (iout,*) "Matrix GGinv"
 528         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
 529       endif
 530       not_done=.true.
 531       iter=0
 532       do while (not_done)
 533       iter=iter+1
 534       if (iter.gt.max_rattle) then
 535         write (iout,*) "Error - too many iterations in RATTLE."
 536         stop
 537       endif
 538 c Calculate the matrix C = GG**(-1) dC_old o dC
 539       ind1=0
 540       do i=nnt,nct-1
 541         ind1=ind1+1
 542         do j=1,3
 543           dC_uncor(j,ind1)=dC(j,i)
 544         enddo
 545       enddo
 546       do i=nnt,nct
 547         if (itype(i).ne.10) then
 548           ind1=ind1+1
 549           do j=1,3
 550             dC_uncor(j,ind1)=dC(j,i+nres)
 551           enddo
 552         endif
 553       enddo
 554       do i=1,nbond
 555         ind=0
 556         do k=nnt,nct-1
 557           ind=ind+1
 558           do j=1,3
 559             gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
 560           enddo
 561         enddo
 562         do k=nnt,nct
 563           if (itype(k).ne.10) then
 564             ind=ind+1
 565             do j=1,3
 566               gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
 567             enddo
 568           endif
 569         enddo
 570       enddo
 571 c Calculate deviations from standard virtual-bond lengths
 572       ind=0
 573       do i=nnt,nct-1
 574         ind=ind+1
 575         x(ind)=vbld(i+1)**2-vbl**2
 576       enddo
 577       do i=nnt,nct
 578         if (itype(i).ne.10) then
 579           ind=ind+1
 580           x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
 581         endif
 582       enddo
 583       if (lprn) then
 584         write (iout,*) "Coordinates and violations"
 585         do i=1,nbond
 586           write(iout,'(i5,3f10.5,5x,e15.5)')
 587      &     i,(dC_uncor(j,i),j=1,3),x(i)
 588         enddo
 589         write (iout,*) "Velocities and violations"
 590         ind=0
 591         do i=nnt,nct-1
 592           ind=ind+1
 593           write (iout,'(2i5,3f10.5,5x,e15.5)')
 594      &     i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
 595         enddo
 596         do i=nnt,nct
 597           if (itype(i).ne.10) then
 598             ind=ind+1
 599             write (iout,'(2i5,3f10.5,5x,e15.5)')
 600      &       i+nres,ind,(d_t(j,i+nres),j=1,3),
 601      &       scalar(d_t(1,i+nres),dC_old(1,i+nres))
 602           endif
 603         enddo
 604         write (iout,*) "gdc"
 605         do i=1,nbond
 606           write (iout,*) "i",i
 607           do j=1,nbond
 608             write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
 609           enddo
 610         enddo
 611       endif
 612       xmax=dabs(x(1))
 613       do i=2,nbond
 614         if (dabs(x(i)).gt.xmax) then
 615           xmax=dabs(x(i))
 616         endif
 617       enddo
 618       if (xmax.lt.tol_rattle) then
 619         not_done=.false.
 620         goto 100
 621       endif
 622 c Calculate the matrix of the system of equations
 623       do i=1,nbond
 624         do j=1,nbond
 625           Cmat(i,j)=0.0d0
 626           do k=1,3
 627             Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
 628           enddo
 629         enddo
 630       enddo
 631       if (lprn1) then
 632         write (iout,*) "Matrix Cmat"
 633         call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
 634       endif
 635       call gauss(Cmat,X,MAXRES2,nbond,1,*10)
 636 c Add constraint term to positions
 637       ind=0
 638       do i=nnt,nct-1
 639         ind=ind+1
 640         do j=1,3
 641           xx=0.0d0
 642           do ii=1,nbond
 643             xx = xx+x(ii)*gdc(j,ind,ii)
 644           enddo
 645           xx=-0.5d0*xx
 646           d_t(j,i)=d_t(j,i)+xx/d_time
 647           dC(j,i)=dC(j,i)+xx
 648         enddo
 649       enddo
 650       do i=nnt,nct
 651         if (itype(i).ne.10) then
 652           ind=ind+1
 653           do j=1,3
 654             xx=0.0d0
 655             do ii=1,nbond
 656               xx = xx+x(ii)*gdc(j,ind,ii)
 657             enddo
 658             xx=-0.5d0*xx
 659             d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
 660             dC(j,i+nres)=dC(j,i+nres)+xx
 661           enddo
 662         endif
 663       enddo
 664 c Rebuild the chain using the new coordinates
 665       call chainbuild_cart
 666       if (lprn) then
 667         write (iout,*) "New coordinates, Lagrange multipliers,",
 668      &  " and differences between actual and standard bond lengths"
 669         ind=0
 670         do i=nnt,nct-1
 671           ind=ind+1
 672           xx=vbld(i+1)**2-vbl**2
 673           write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 674      &        i,(dC(j,i),j=1,3),x(ind),xx
 675         enddo
 676         do i=nnt,nct
 677           if (itype(i).ne.10) then
 678             ind=ind+1
 679             xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
 680             write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
 681      &       i,(dC(j,i+nres),j=1,3),x(ind),xx
 682           endif
 683         enddo
 684         write (iout,*) "Velocities and violations"
 685         ind=0
 686         do i=nnt,nct-1
 687           ind=ind+1
 688           write (iout,'(2i5,3f10.5,5x,e15.5)')
 689      &     i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
 690         enddo
 691         do i=nnt,nct
 692           if (itype(i).ne.10) then
 693             ind=ind+1
 694             write (iout,'(2i5,3f10.5,5x,e15.5)')
 695      &       i+nres,ind,(d_t_new(j,i+nres),j=1,3),
 696      &       scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
 697           endif
 698         enddo
 699       endif
 700       enddo
 701   100 continue
 702       return
 703    10 write (iout,*) "Error - singularity in solving the system",
 704      & " of equations for Lagrange multipliers."
 705       stop
 706       end