+++ /dev/null
- subroutine rattle1
-c RATTLE algorithm for velocity Verlet - step 1, UNRES
-c AL 9/24/04
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision gginv(maxres2,maxres2),
- & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
- & Cmat(MAXRES2,MAXRES2),x(MAXRES2),xcorr(3,MAXRES2)
- common /przechowalnia/ GGinv,gdc,Cmat,nbond
- integer max_rattle /5/
- logical lprn /.false./, lprn1 /.false./,not_done
- double precision tol_rattle /1.0d-5/
- if (lprn) write (iout,*) "RATTLE1"
- nbond=nct-nnt
- do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
- enddo
-c Make a folded form of the Ginv-matrix
- ind=0
- ii=0
- do i=nnt,nct-1
- ii=ii+1
- do j=1,3
- ind=ind+1
- ind1=0
- jj=0
- do k=nnt,nct-1
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
- enddo
- endif
- enddo
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ii=ii+1
- do j=1,3
- ind=ind+1
- ind1=0
- jj=0
- do k=nnt,nct-1
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=Ginv(ind,ind1)
- enddo
- endif
- enddo
- enddo
- endif
- enddo
- if (lprn1) then
- write (iout,*) "Matrix GGinv"
- call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
- endif
- not_done=.true.
- iter=0
- do while (not_done)
- iter=iter+1
- if (iter.gt.max_rattle) then
- write (iout,*) "Error - too many iterations in RATTLE."
- stop
- endif
-c Calculate the matrix C = GG**(-1) dC_old o dC
- ind1=0
- do i=nnt,nct-1
- ind1=ind1+1
- do j=1,3
- dC_uncor(j,ind1)=dC(j,i)
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind1=ind1+1
- do j=1,3
- dC_uncor(j,ind1)=dC(j,i+nres)
- enddo
- endif
- enddo
- do i=1,nbond
- ind=0
- do k=nnt,nct-1
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
- enddo
- endif
- enddo
- enddo
-c Calculate deviations from standard virtual-bond lengths
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- x(ind)=vbld(i+1)**2-vbl**2
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
- endif
- enddo
- if (lprn) then
- write (iout,*) "Coordinates and violations"
- do i=1,nbond
- write(iout,'(i5,3f10.5,5x,e15.5)')
- & i,(dC_uncor(j,i),j=1,3),x(i)
- enddo
- write (iout,*) "Velocities and violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i+nres,ind,(d_t_new(j,i+nres),j=1,3),
- & scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
- endif
- enddo
-c write (iout,*) "gdc"
-c do i=1,nbond
-c write (iout,*) "i",i
-c do j=1,nbond
-c write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
-c enddo
-c enddo
- endif
- xmax=dabs(x(1))
- do i=2,nbond
- if (dabs(x(i)).gt.xmax) then
- xmax=dabs(x(i))
- endif
- enddo
- if (xmax.lt.tol_rattle) then
- not_done=.false.
- goto 100
- endif
-c Calculate the matrix of the system of equations
- do i=1,nbond
- do j=1,nbond
- Cmat(i,j)=0.0d0
- do k=1,3
- Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
- enddo
- enddo
- enddo
- if (lprn1) then
- write (iout,*) "Matrix Cmat"
- call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
- endif
- call gauss(Cmat,X,MAXRES2,nbond,1,*10)
-c Add constraint term to positions
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- xx=0.5d0*xx
- dC(j,i)=dC(j,i)-xx
- d_t_new(j,i)=d_t_new(j,i)-xx/d_time
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- xx=0.5d0*xx
- dC(j,i+nres)=dC(j,i+nres)-xx
- d_t_new(j,i+nres)=d_t_new(j,i+nres)-xx/d_time
- enddo
- endif
- enddo
-c Rebuild the chain using the new coordinates
- call chainbuild_cart
- if (lprn) then
- write (iout,*) "New coordinates, Lagrange multipliers,",
- & " and differences between actual and standard bond lengths"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- xx=vbld(i+1)**2-vbl**2
- write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
- & i,(dC(j,i),j=1,3),x(ind),xx
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
- write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
- & i,(dC(j,i+nres),j=1,3),x(ind),xx
- endif
- enddo
- write (iout,*) "Velocities and violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i+nres,ind,(d_t_new(j,i+nres),j=1,3),
- & scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
- endif
- enddo
- endif
- enddo
- 100 continue
- return
- 10 write (iout,*) "Error - singularity in solving the system",
- & " of equations for Lagrange multipliers."
- stop
- end
-c------------------------------------------------------------------------------
- subroutine rattle2
-c RATTLE algorithm for velocity Verlet - step 2, UNRES
-c AL 9/24/04
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision gginv(maxres2,maxres2),
- & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
- & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
- common /przechowalnia/ GGinv,gdc,Cmat,nbond
- integer max_rattle /5/
- logical lprn /.false./, lprn1 /.false./,not_done
- double precision tol_rattle /1.0d-5/
- if (lprn) write (iout,*) "RATTLE2"
- if (lprn) write (iout,*) "Velocity correction"
-c Calculate the matrix G dC
- do i=1,nbond
- ind=0
- do k=nnt,nct-1
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC(j,k)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC(j,k+nres)
- enddo
- endif
- enddo
- enddo
-c if (lprn) then
-c write (iout,*) "gdc"
-c do i=1,nbond
-c write (iout,*) "i",i
-c do j=1,nbond
-c write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
-c enddo
-c enddo
-c endif
-c Calculate the matrix of the system of equations
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- do j=1,nbond
- Cmat(ind,j)=0.0d0
- do k=1,3
- Cmat(ind,j)=Cmat(ind,j)+dC(k,i)*gdc(k,ind,j)
- enddo
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- do j=1,nbond
- Cmat(ind,j)=0.0d0
- do k=1,3
- Cmat(ind,j)=Cmat(ind,j)+dC(k,i+nres)*gdc(k,ind,j)
- enddo
- enddo
- endif
- enddo
-c Calculate the scalar product dC o d_t_new
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- x(ind)=scalar(d_t(1,i),dC(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- x(ind)=scalar(d_t(1,i+nres),dC(1,i+nres))
- endif
- enddo
- if (lprn) then
- write (iout,*) "Velocities and violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i,ind,(d_t(j,i),j=1,3),x(ind)
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind)
- endif
- enddo
- endif
- xmax=dabs(x(1))
- do i=2,nbond
- if (dabs(x(i)).gt.xmax) then
- xmax=dabs(x(i))
- endif
- enddo
- if (xmax.lt.tol_rattle) then
- not_done=.false.
- goto 100
- endif
- if (lprn1) then
- write (iout,*) "Matrix Cmat"
- call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
- endif
- call gauss(Cmat,X,MAXRES2,nbond,1,*10)
-c Add constraint term to velocities
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- d_t(j,i)=d_t(j,i)-xx
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- d_t(j,i+nres)=d_t(j,i+nres)-xx
- enddo
- endif
- enddo
- if (lprn) then
- write (iout,*)
- & "New velocities, Lagrange multipliers violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- if (lprn) write (iout,'(2i5,3f10.5,5x,2e15.5)')
- & i,ind,(d_t(j,i),j=1,3),x(ind),scalar(d_t(1,i),dC(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,2e15.5)')
- & i+nres,ind,(d_t(j,i+nres),j=1,3),x(ind),
- & scalar(d_t(1,i+nres),dC(1,i+nres))
- endif
- enddo
- endif
- 100 continue
- return
- 10 write (iout,*) "Error - singularity in solving the system",
- & " of equations for Lagrange multipliers."
- stop
- end
-c------------------------------------------------------------------------------
- subroutine rattle_brown
-c RATTLE/LINCS algorithm for Brownian dynamics, UNRES
-c AL 9/24/04
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision gginv(maxres2,maxres2),
- & gdc(3,MAXRES2,MAXRES2),dC_uncor(3,MAXRES2),
- & Cmat(MAXRES2,MAXRES2),x(MAXRES2)
- common /przechowalnia/ GGinv,gdc,Cmat,nbond
- integer max_rattle /5/
- logical lprn /.true./, lprn1 /.true./,not_done
- double precision tol_rattle /1.0d-5/
- if (lprn) write (iout,*) "RATTLE_BROWN"
- nbond=nct-nnt
- do i=nnt,nct
- if (itype(i).ne.10) nbond=nbond+1
- enddo
-c Make a folded form of the Ginv-matrix
- ind=0
- ii=0
- do i=nnt,nct-1
- ii=ii+1
- do j=1,3
- ind=ind+1
- ind1=0
- jj=0
- do k=nnt,nct-1
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
- enddo
- endif
- enddo
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ii=ii+1
- do j=1,3
- ind=ind+1
- ind1=0
- jj=0
- do k=nnt,nct-1
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1) GGinv(ii,jj)=fricmat(ind,ind1)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- jj=jj+1
- do l=1,3
- ind1=ind1+1
- if (j.eq.1 .and. l.eq.1)GGinv(ii,jj)=fricmat(ind,ind1)
- enddo
- endif
- enddo
- enddo
- endif
- enddo
- if (lprn1) then
- write (iout,*) "Matrix GGinv"
- call MATOUT(nbond,nbond,MAXRES2,MAXRES2,GGinv)
- endif
- not_done=.true.
- iter=0
- do while (not_done)
- iter=iter+1
- if (iter.gt.max_rattle) then
- write (iout,*) "Error - too many iterations in RATTLE."
- stop
- endif
-c Calculate the matrix C = GG**(-1) dC_old o dC
- ind1=0
- do i=nnt,nct-1
- ind1=ind1+1
- do j=1,3
- dC_uncor(j,ind1)=dC(j,i)
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind1=ind1+1
- do j=1,3
- dC_uncor(j,ind1)=dC(j,i+nres)
- enddo
- endif
- enddo
- do i=1,nbond
- ind=0
- do k=nnt,nct-1
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k)
- enddo
- enddo
- do k=nnt,nct
- if (itype(k).ne.10) then
- ind=ind+1
- do j=1,3
- gdc(j,i,ind)=GGinv(i,ind)*dC_old(j,k+nres)
- enddo
- endif
- enddo
- enddo
-c Calculate deviations from standard virtual-bond lengths
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- x(ind)=vbld(i+1)**2-vbl**2
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- x(ind)=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
- endif
- enddo
- if (lprn) then
- write (iout,*) "Coordinates and violations"
- do i=1,nbond
- write(iout,'(i5,3f10.5,5x,e15.5)')
- & i,(dC_uncor(j,i),j=1,3),x(i)
- enddo
- write (iout,*) "Velocities and violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i,ind,(d_t(j,i),j=1,3),scalar(d_t(1,i),dC_old(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i+nres,ind,(d_t(j,i+nres),j=1,3),
- & scalar(d_t(1,i+nres),dC_old(1,i+nres))
- endif
- enddo
- write (iout,*) "gdc"
- do i=1,nbond
- write (iout,*) "i",i
- do j=1,nbond
- write (iout,'(i5,3f10.5)') j,(gdc(k,j,i),k=1,3)
- enddo
- enddo
- endif
- xmax=dabs(x(1))
- do i=2,nbond
- if (dabs(x(i)).gt.xmax) then
- xmax=dabs(x(i))
- endif
- enddo
- if (xmax.lt.tol_rattle) then
- not_done=.false.
- goto 100
- endif
-c Calculate the matrix of the system of equations
- do i=1,nbond
- do j=1,nbond
- Cmat(i,j)=0.0d0
- do k=1,3
- Cmat(i,j)=Cmat(i,j)+dC_uncor(k,i)*gdc(k,i,j)
- enddo
- enddo
- enddo
- if (lprn1) then
- write (iout,*) "Matrix Cmat"
- call MATOUT(nbond,nbond,MAXRES2,MAXRES2,Cmat)
- endif
- call gauss(Cmat,X,MAXRES2,nbond,1,*10)
-c Add constraint term to positions
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- xx=-0.5d0*xx
- d_t(j,i)=d_t(j,i)+xx/d_time
- dC(j,i)=dC(j,i)+xx
- enddo
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- do j=1,3
- xx=0.0d0
- do ii=1,nbond
- xx = xx+x(ii)*gdc(j,ind,ii)
- enddo
- xx=-0.5d0*xx
- d_t(j,i+nres)=d_t(j,i+nres)+xx/d_time
- dC(j,i+nres)=dC(j,i+nres)+xx
- enddo
- endif
- enddo
-c Rebuild the chain using the new coordinates
- call chainbuild_cart
- if (lprn) then
- write (iout,*) "New coordinates, Lagrange multipliers,",
- & " and differences between actual and standard bond lengths"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- xx=vbld(i+1)**2-vbl**2
- write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
- & i,(dC(j,i),j=1,3),x(ind),xx
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- xx=vbld(i+nres)**2-vbldsc0(1,itype(i))**2
- write (iout,'(i5,3f10.5,5x,f10.5,e15.5)')
- & i,(dC(j,i+nres),j=1,3),x(ind),xx
- endif
- enddo
- write (iout,*) "Velocities and violations"
- ind=0
- do i=nnt,nct-1
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i,ind,(d_t_new(j,i),j=1,3),scalar(d_t_new(1,i),dC_old(1,i))
- enddo
- do i=nnt,nct
- if (itype(i).ne.10) then
- ind=ind+1
- write (iout,'(2i5,3f10.5,5x,e15.5)')
- & i+nres,ind,(d_t_new(j,i+nres),j=1,3),
- & scalar(d_t_new(1,i+nres),dC_old(1,i+nres))
- endif
- enddo
- endif
- enddo
- 100 continue
- return
- 10 write (iout,*) "Error - singularity in solving the system",
- & " of equations for Lagrange multipliers."
- stop
- end
projects / unres.git / blobdiff