Adam's unres update

[unres.git] / source / unres / src-HCD-5D / orig_frame_chain.F

diff --git a/source/unres/src-HCD-5D/orig_frame_chain.F b/source/unres/src-HCD-5D/orig_frame_chain.F
new file mode 100644 (file)
index 0000000..07053ce
--- /dev/null
+++ b/source/unres/src-HCD-5D/orig_frame_chain.F
@@ -0,0 +1,85 @@
+      subroutine orig_frame_chain(istart)
+C
+C Define the origin and orientation of the coordinate system starting
+C at residue istart and locate sites istart+1 and istart+2. The
+C coordinates of site istart and the respective dc and dc_norm must be
+C pre-defined
+C
+      implicit none
+      integer i,j,istart,ichain
+      double precision cost,sint,cosg,sing,aux
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      cost=dc_norm(1,istart)
+      aux=dsqrt(dc_norm(2,istart)**2+dc_norm(3,istart)**2)
+      cosg=dc_norm(2,istart)/aux
+      sing=dc_norm(3,istart)/aux
+      cost=dcos(theta(istart+2))
+      sint=dsin(theta(istart+2))
+      prod(1,1,istart)=cost
+      prod(1,2,istart)=-sint
+      prod(1,3,istart)=0.0d0
+      prod(2,1,istart)=sint*cosg
+      prod(2,2,istart)=cost*cosg
+      prod(2,3,istart)=-sing
+      prod(3,1,istart)=sint*sing
+      prod(3,2,istart)=cost*sing
+      prod(3,3,istart)=cosg
+      aux=prod(1,1,istart)*(prod(2,2,istart)*prod(3,3,istart)
+     & -prod(3,2,istart)*prod(2,3,istart))
+     & -prod(1,2,istart)*(prod(2,1,istart)*prod(3,3,istart)
+     & -prod(3,1,istart)*prod(2,3,istart))
+     & +prod(1,3,istart)*(prod(2,1,istart)*prod(3,2,istart)
+     & -prod(3,1,istart)*prod(2,2,istart))
+c      write (iout,*) "orig_frame_chain prod",istart
+c      do i=1,3
+c        write(iout,'(i5,3f10.5)') i,(prod(i,j,istart),j=1,3)
+c      enddo
+c      write (iout,*) "orig_frame_chain: prod",istart," determinant",aux
+      t(1,1,istart)=cost
+      t(1,2,istart)=-sint 
+      t(1,3,istart)= 0.0D0
+      t(2,1,istart)=sint
+      t(2,2,istart)=cost
+      t(2,3,istart)= 0.0D0
+      t(3,1,istart)= 0.0D0
+      t(3,2,istart)= 0.0D0
+      t(3,3,istart)= 1.0D0
+      r(1,1,istart)= 1.0D0
+      r(1,2,istart)= 0.0D0
+      r(1,3,istart)= 0.0D0
+      r(2,1,istart)= 0.0D0
+      r(2,2,istart)= 1.0D0
+      r(2,3,istart)= 0.0D0
+      r(3,1,istart)= 0.0D0
+      r(3,2,istart)= 0.0D0
+      r(3,3,istart)= 1.0D0
+      do i=1,3
+        do j=1,3
+          rt(i,j,istart)=t(i,j,istart)
+        enddo
+      enddo
+      call matmult(prod(1,1,istart),rt(1,1,istart),prod(1,1,istart+1))
+c      aux=prod(1,1,istart+1)*(prod(2,2,istart+1)*prod(3,3,istart+1)
+c     & -prod(3,2,istart+1)*prod(2,3,istart+1))
+c     & -prod(1,2,istart+1)*(prod(2,1,istart+1)*prod(3,3,istart+1)
+c     & -prod(3,1,istart+1)*prod(2,3,istart+1))
+c     & +prod(1,3,istart+1)*(prod(2,1,istart+1)*prod(3,2,istart+1)
+c     & -prod(3,1,istart+1)*prod(2,2,istart+1))
+c      write (iout,*) "orig_frame_chain prod",istart+1
+c      do i=1,3
+c        write(iout,'(i5,3f10.5)') i,(prod(i,j,istart+1),j=1,3)
+c      enddo
+c      write (iout,*)"orig_frame_chain: prod",istart+1," determinant",aux
+      do j=1,3
+        dc_norm(j,istart+1)=prod(j,1,istart+1)
+        dc(j,istart+1)=vbld(istart+2)*prod(j,1,istart+1)
+        c(j,istart+2)=c(j,istart+1)+dc(j,istart+1)
+      enddo
+      call locate_side_chain(istart+1)
+      return
+      end