Merge branch 'homology' of mmka.chem.univ.gda.pl:unres into homology
[unres.git] / source / unres / src_MD-M / 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