--- /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'
+#ifdef RATTLE
+ 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
+#else
+ write (iout,*)
+ & "RATTLE inactive; use -DRATTLE switch at compile time."
+ stop
+#endif
+ 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'
+#ifdef RATTLE
+ 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
+#else
+ write (iout,*)
+ & "RATTLE inactive; use -DRATTLE option at compile time."
+ stop
+#endif
+ 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'
+#ifdef RATTLE
+ 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
+#else
+ write (iout,*)
+ & "RATTLE inactive; use -DRATTLE option at compile time"
+ stop
+#endif
+ end
projects / unres.git / blobdiff