source/unres/src-HCD-5D/rattle.F

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