2 c RATTLE algorithm for velocity Verlet - step 1, UNRES
6 include 'COMMON.IOUNITS'
8 include 'COMMON.CONTROL'
12 include 'COMMON.LAGRANGE.5diag'
14 include 'COMMON.LAGRANGE'
17 include 'COMMON.LANGEVIN'
20 include 'COMMON.LANGEVIN.lang0.5diag'
22 include 'COMMON.LANGEVIN.lang0'
25 include 'COMMON.CHAIN'
26 include 'COMMON.DERIV'
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"
42 if (itype(i).ne.10) nbond=nbond+1
44 c Make a folded form of the Ginv-matrix
57 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
61 if (itype(k).ne.10) then
65 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
72 if (itype(i).ne.10) then
82 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
86 if (itype(k).ne.10) then
90 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
98 write (iout,*) "Matrix GGinv"
99 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
105 if (iter.gt.max_rattle) then
106 write (iout,*) "Error - too many iterations in RATTLE."
109 c Calculate the matrix C = GG**(-1) dC_old o dC
114 dC_uncor(j,ind1)=dC(j,i)
118 if (itype(i).ne.10) then
121 dC_uncor(j,ind1)=dC(j,i+nres)
130 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
134 if (itype(k).ne.10) then
137 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
142 c Calculate deviations from standard virtual-bond lengths
146 x(ind)=vbld(i+1)**2-vbl**2
149 if (itype(i).ne.10) then
151 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
155 write (iout,*) "Coordinates and violations"
157 write(iout,'(i5,3f10.5,5x,e15.5)')
158 & i,(dC_uncor(j,i),j=1,3),x(i)
160 write (iout,*) "Velocities and violations"
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))
168 if (itype(i).ne.10) then
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))
175 c write (iout,*) "gdc"
177 c write (iout,*) "i",i
179 c write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
185 if (dabs(x(i)).gt.xmax) then
189 if (xmax.lt.tol_rattle) then
193 c Calculate the matrix of the system of equations
198 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
203 write (iout,*) "Matrix Cmat"
204 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
206 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
207 c Add constraint term to positions
214 xx = xx+x(ii)*gdc(j,ind,ii)
218 d_t_new(j,i)=d_t_new(j,i)-xx/d_time
222 if (itype(i).ne.10) then
227 xx = xx+x(ii)*gdc(j,ind,ii)
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
235 c Rebuild the chain using the new coordinates
238 write (iout,*) "New coordinates, Lagrange multipliers,",
239 & " and differences between actual and standard bond lengths"
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
248 if (itype(i).ne.10) then
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
255 write (iout,*) "Velocities and violations"
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))
263 if (itype(i).ne.10) then
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))
274 10 write (iout,*) "Error - singularity in solving the system",
275 & " of equations for Lagrange multipliers."
279 & "RATTLE inactive; use -DRATTLE switch at compile time."
283 c------------------------------------------------------------------------------
285 c RATTLE algorithm for velocity Verlet - step 2, UNRES
289 include 'COMMON.IOUNITS'
291 include 'COMMON.CONTROL'
295 include 'COMMON.LAGRANGE.5diag'
297 include 'COMMON.LAGRANGE'
300 include 'COMMON.LANGEVIN'
303 include 'COMMON.LANGEVIN.lang0.5diag'
305 include 'COMMON.LANGEVIN.lang0'
308 include 'COMMON.CHAIN'
309 include 'COMMON.DERIV'
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
331 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
335 if (itype(k).ne.10) then
338 gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
344 c write (iout,*) "gdc"
346 c write (iout,*) "i",i
348 c write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
352 c Calculate the matrix of the system of equations
359 Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
364 if (itype(i).ne.10) then
369 Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
374 c Calculate the scalar product dC o d_t_new
378 x(ind)=scalar(d_t(1,i),dC(1,i))
381 if (itype(i).ne.10) then
383 x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
387 write (iout,*) "Velocities and violations"
391 write (iout,'(2i5,3f10.5,5x,e15.5)')
392 & i,ind,(d_t(j,i),j=1,3),x(ind)
395 if (itype(i).ne.10) then
397 write (iout,'(2i5,3f10.5,5x,e15.5)')
398 & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
404 if (dabs(x(i)).gt.xmax) then
408 if (xmax.lt.tol_rattle) then
413 write (iout,*) "Matrix Cmat"
414 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
416 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
417 c Add constraint term to velocities
424 xx = xx+x(ii)*gdc(j,ind,ii)
430 if (itype(i).ne.10) then
435 xx = xx+x(ii)*gdc(j,ind,ii)
437 d_t(j,i+nres)=d_t(j,i+nres)-xx
443 & "New velocities, Lagrange multipliers violations"
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))
451 if (itype(i).ne.10) then
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))
461 10 write (iout,*) "Error - singularity in solving the system",
462 & " of equations for Lagrange multipliers."
466 & "RATTLE inactive; use -DRATTLE option at compile time."
470 c------------------------------------------------------------------------------
471 subroutine rattle_brown
472 c RATTLE/LINCS algorithm for Brownian dynamics, UNRES
474 implicit real*8 (a-h,o-z)
477 include 'COMMON.CONTROL'
481 include 'COMMON.LAGRANGE.5diag'
483 include 'COMMON.LAGRANGE'
486 include 'COMMON.LANGEVIN'
489 include 'COMMON.LANGEVIN.lang0.5diag'
491 include 'COMMON.LANGEVIN.lang0'
494 include 'COMMON.CHAIN'
495 include 'COMMON.DERIV'
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"
512 if (itype(i).ne.10) nbond=nbond+1
514 c Make a folded form of the Ginv-matrix
527 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
531 if (itype(k).ne.10) then
535 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
542 if (itype(i).ne.10) then
552 if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
556 if (itype(k).ne.10) then
560 if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
568 write (iout,*) "Matrix GGinv"
569 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
575 if (iter.gt.max_rattle) then
576 write (iout,*) "Error - too many iterations in RATTLE."
579 c Calculate the matrix C = GG**(-1) dC_old o dC
584 dC_uncor(j,ind1)=dC(j,i)
588 if (itype(i).ne.10) then
591 dC_uncor(j,ind1)=dC(j,i+nres)
600 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
604 if (itype(k).ne.10) then
607 gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
612 c Calculate deviations from standard virtual-bond lengths
616 x(ind)=vbld(i+1)**2-vbl**2
619 if (itype(i).ne.10) then
621 x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
625 write (iout,*) "Coordinates and violations"
627 write(iout,'(i5,3f10.5,5x,e15.5)')
628 & i,(dC_uncor(j,i),j=1,3),x(i)
630 write (iout,*) "Velocities and violations"
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))
638 if (itype(i).ne.10) then
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))
649 write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
655 if (dabs(x(i)).gt.xmax) then
659 if (xmax.lt.tol_rattle) then
663 c Calculate the matrix of the system of equations
668 Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
673 write (iout,*) "Matrix Cmat"
674 call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
676 call gauss(Cmat,X,MAXRES2,nbond,1,*10)
677 c Add constraint term to positions
684 xx = xx+x(ii)*gdc(j,ind,ii)
687 d_t(j,i)=d_t(j,i)+xx/d_time
692 if (itype(i).ne.10) then
697 xx = xx+x(ii)*gdc(j,ind,ii)
700 d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
701 dC(j,i+nres)=dC(j,i+nres)+xx
705 c Rebuild the chain using the new coordinates
708 write (iout,*) "New coordinates, Lagrange multipliers,",
709 & " and differences between actual and standard bond lengths"
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
718 if (itype(i).ne.10) then
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
725 write (iout,*) "Velocities and violations"
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))
733 if (itype(i).ne.10) then
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))
744 10 write (iout,*) "Error - singularity in solving the system",
745 & " of equations for Lagrange multipliers."
749 & "RATTLE inactive; use -DRATTLE option at compile time"