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