Merge branch 'lipid' of mmka.chem.univ.gda.pl:unres into lipid

[unres.git] / source / unres / src_Eshel / bond_move.f

diff --git a/source/unres/src_Eshel/bond_move.f b/source/unres/src_Eshel/bond_move.f
new file mode 100644 (file)
index 0000000..4843f60
--- /dev/null
+++ b/source/unres/src_Eshel/bond_move.f
@@ -0,0 +1,124 @@
+      subroutine bond_move(nbond,nstart,psi,lprint,error)
+C Move NBOND fragment starting from the CA(nstart) by angle PSI.
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      integer nbond,nstart
+      double precision psi
+      logical fail,error,lprint
+      include 'COMMON.GEO'
+      include 'COMMON.CHAIN'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.MCM'
+      dimension x(3),e(3,3),e1(3),e2(3),e3(3),rot(3,3),trans(3,3)
+      error=.false.
+      nend=nstart+nbond
+      if (print_mc.gt.2) then
+      write (iout,*) 'nstart=',nstart,' nend=',nend,' nbond=',nbond
+      write (iout,*) 'psi=',psi
+      write (iout,'(a)') 'Original coordinates of the fragment'
+      do i=nstart,nend
+        write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
+      enddo
+      endif
+      if (nstart.lt.1 .or. nend .gt.nres .or. nbond.lt.2 .or. 
+     & nbond.ge.nres-1) then
+        write (iout,'(a)') 'Bad data in BOND_MOVE.'
+        error=.true.
+        return
+      endif
+C Generate the reference system.
+      i2=nend
+      i3=nstart
+      i4=nstart+1
+      call refsys(i2,i3,i4,e1,e2,e3,error) 
+C Return, if couldn't define the reference system.
+      if (error) return
+C Compute the transformation matrix.
+      cospsi=dcos(psi)
+      sinpsi=dsin(psi)
+      rot(1,1)=1.0D0
+      rot(1,2)=0.0D0
+      rot(1,3)=0.0D0
+      rot(2,1)=0.0D0
+      rot(2,2)=cospsi
+      rot(2,3)=-sinpsi
+      rot(3,1)=0.0D0
+      rot(3,2)=sinpsi
+      rot(3,3)=cospsi
+      do i=1,3
+        e(1,i)=e1(i)
+        e(2,i)=e2(i)
+        e(3,i)=e3(i)
+      enddo
+
+      if (print_mc.gt.2) then
+      write (iout,'(a)') 'Reference system and matrix r:'
+      do i=1,3
+        write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
+      enddo
+      endif
+
+      call matmult(rot,e,trans)
+      do i=1,3
+        do j=1,3
+          e(i,1)=e1(i)
+          e(i,2)=e2(i)
+          e(i,3)=e3(i)
+        enddo
+      enddo
+      call matmult(e,trans,trans)
+
+      if (lprint) then
+      write (iout,'(a)') 'The trans matrix:'
+      do i=1,3
+        write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
+      enddo
+      endif
+
+      do i=nstart,nend
+        do j=1,3
+          rij=c(j,nstart)
+          do k=1,3
+            rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
+          enddo
+          x(j)=rij
+        enddo
+        do j=1,3
+          c(j,i)=x(j)
+        enddo
+      enddo
+
+      if (lprint) then
+      write (iout,'(a)') 'Rotated coordinates of the fragment'
+      do i=nstart,nend
+        write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
+      enddo
+      endif
+
+c     call int_from_cart(.false.,lprint)
+      if (nstart.gt.1) then
+        theta(nstart+1)=alpha(nstart-1,nstart,nstart+1)
+        phi(nstart+2)=beta(nstart-1,nstart,nstart+1,nstart+2)
+        if (nstart.gt.2) phi(nstart+1)=
+     &                beta(nstart-2,nstart-1,nstart,nstart+1)
+      endif
+      if (nend.lt.nres) then
+        theta(nend+1)=alpha(nend-1,nend,nend+1)
+        phi(nend+1)=beta(nend-2,nend-1,nend,nend+1)
+        if (nend.lt.nres-1) phi(nend+2)=
+     &                beta(nend-1,nend,nend+1,nend+2)
+      endif
+      if (print_mc.gt.2) then
+      write (iout,'(/a,i3,a,i3,a/)') 
+     & 'Moved internal coordinates of the ',nstart,'-',nend,
+     & ' fragment:'
+      do i=nstart+1,nstart+2
+        write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
+      enddo
+      do i=nend+1,nend+2
+        write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
+      enddo
+      endif
+      return
+      end