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