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