Adam's unres update
[unres.git] / source / unres / src-HCD-5D / chainbuild.F
index 51419ef..a60e2bd 100644 (file)
@@ -108,11 +108,11 @@ C
       include 'COMMON.VAR'
       cost=dcos(theta(3))
       sint=dsin(theta(3))
-      t(1,1,1)=-cost
+      t(1,1,1)=cost
       t(1,2,1)=-sint 
       t(1,3,1)= 0.0D0
-      t(2,1,1)=-sint
-      t(2,2,1)= cost
+      t(2,1,1)=sint
+      t(2,2,1)=cost
       t(2,3,1)= 0.0D0
       t(3,1,1)= 0.0D0
       t(3,2,1)= 0.0D0
@@ -170,7 +170,7 @@ C
 C Locate CA(i) and SC(i-1)
 C
       implicit none
-      integer i,j
+      integer i,j,k
       double precision theti,phii,cost,sint,cosphi,sinphi
       include 'DIMENSIONS'
       include 'COMMON.CHAIN'
@@ -196,14 +196,15 @@ C
       sint=dsin(theti)
       cosphi=dcos(phii)
       sinphi=dsin(phii)
+c      write (iout,*) "locate_next_res i",i
 * Define the matrices of the rotation about the virtual-bond valence angles
 * theta, T(i,j,k), virtual-bond dihedral angles gamma (miscalled PHI in this
 * program), R(i,j,k), and, the cumulative matrices of rotation RT
       t(1,1,i-2)=-cost
-      t(1,2,i-2)=-sint 
+      t(1,2,i-2)= sint 
       t(1,3,i-2)= 0.0D0
       t(2,1,i-2)=-sint
-      t(2,2,i-2)= cost
+      t(2,2,i-2)=-cost
       t(2,3,i-2)= 0.0D0
       t(3,1,i-2)= 0.0D0
       t(3,2,i-2)= 0.0D0
@@ -212,26 +213,36 @@ C
       r(1,2,i-2)= 0.0D0
       r(1,3,i-2)= 0.0D0
       r(2,1,i-2)= 0.0D0
-      r(2,2,i-2)=-cosphi
-      r(2,3,i-2)= sinphi
+      r(2,2,i-2)= cosphi
+      r(2,3,i-2)=-sinphi
       r(3,1,i-2)= 0.0D0
       r(3,2,i-2)= sinphi
       r(3,3,i-2)= cosphi
       rt(1,1,i-2)=-cost
-      rt(1,2,i-2)=-sint
+      rt(1,2,i-2)= sint
       rt(1,3,i-2)=0.0D0
-      rt(2,1,i-2)=sint*cosphi
+      rt(2,1,i-2)=-sint*cosphi
       rt(2,2,i-2)=-cost*cosphi
-      rt(2,3,i-2)=sinphi
+      rt(2,3,i-2)=-sinphi
       rt(3,1,i-2)=-sint*sinphi
-      rt(3,2,i-2)=cost*sinphi
-      rt(3,3,i-2)=cosphi
+      rt(3,2,i-2)=-cost*sinphi
+      rt(3,3,i-2)= cosphi
       call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+c      write (iout,*) "prod",i-2
+c      do j=1,3
+c        write (iout,*) (prod(j,k,i-2),k=1,3)
+c      enddo
+c      write (iout,*) "prod",i-1
+c      do j=1,3
+c        write (iout,*) (prod(j,k,i-1),k=1,3)
+c      enddo
       do j=1,3
         dc_norm(j,i-1)=prod(j,1,i-1)
         dc(j,i-1)=vbld(i)*prod(j,1,i-1)
         c(j,i)=c(j,i-1)+dc(j,i-1)
       enddo
+c      write (iout,*) "dc",i-1,(dc(j,i-1),j=1,3)
+c      write (iout,*) "c",i,(dc(j,i),j=1,3)
 cd    print '(2i3,2(3f10.5,5x))', i-1,i,(dc(j,i-1),j=1,3),(c(j,i),j=1,3)
 C 
 C Now calculate the coordinates of SC(i-1)
@@ -240,6 +251,86 @@ C
       return
       end
 c-----------------------------------------------------------------------------
+      subroutine sc_coord_rebuild(i)
+C 
+C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
+C
+      implicit none
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      integer i,j,k
+      double precision xx(3)
+      double precision dsci,dsci_inv,alphi,omegi,cosalphi,sinalphi,
+     &  cosomegi,sinomegi,rp(3),theta2,cost2,sint2,rj
+      double precision scalar
+      double precision ref(3,3),scalp,sscalp,refnorm
+      double precision alpha,beta
+c      dsci=dsc(itype(i))
+c      dsci_inv=dsc_inv(itype(i))
+      dsci=vbld(i+nres)
+      dsci_inv=vbld_inv(i+nres)
+#ifdef OSF
+      alphi=alph(i)
+      omegi=omeg(i)
+      if (alphi.ne.alphi) alphi=100.0
+      if (omegi.ne.omegi) omegi=-100.0
+#else
+      alphi=alph(i)
+      omegi=omeg(i)
+#endif
+      cosalphi=dcos(alphi)
+      sinalphi=dsin(alphi)
+      cosomegi=dcos(omegi)
+      sinomegi=dsin(omegi) 
+      rp(1)= cosalphi
+      rp(2)= sinalphi*cosomegi
+      rp(3)=-sinalphi*sinomegi
+c Build the reference system
+      do j=1,3
+        ref(j,1)=-dc_norm(j,i-1)+dc_norm(j,i)
+      enddo
+      refnorm=dsqrt(scalar(ref(1,1),ref(1,1)))
+      do j=1,3
+        ref(j,1)=ref(j,1)/refnorm
+      enddo
+      scalp=scalar(ref(1,1),dc_norm(1,i))
+      sscalp=1.0d0/dsqrt(1.0d0-scalp*scalp)
+      do j=1,3
+        ref(j,2)=(dc_norm(j,i)-scalp*ref(j,1))*sscalp
+      enddo
+      ref(1,3)= ref(2,1)*ref(3,2)-ref(3,1)*ref(2,2)
+      ref(2,3)=-ref(1,1)*ref(3,2)+ref(3,1)*ref(1,2)
+      ref(3,3)= ref(1,1)*ref(2,2)-ref(2,1)*ref(1,2)
+c      do j=1,3
+c        write (iout,*) j,scalar(ref(1,j),ref(1,1)),
+c     &  scalar(ref(1,j),ref(1,2)),scalar(ref(1,j),ref(1,3))
+c      enddo
+c Bring the coordinates to the global reference system
+      do j=1,3
+        dc_norm(j,nres+i)=0.0d0
+        do k=1,3
+          dc_norm(j,nres+i)=dc_norm(j,nres+i)+ref(j,k)*rp(k)
+        enddo
+        dc(j,nres+i)=dc_norm(j,nres+i)*dsci
+        c(j,nres+i)=c(j,i)+dc(j,nres+i)
+      enddo 
+c      write (iout,*) scalar(dc_norm(1,i+nres),dc_norm(1,i+nres)),
+c     &  dsqrt(scalar(dc(1,i+nres),dc(1,i+nres)))
+c Check the internal coordinates
+c      c(:,2*nres+1)=ref(:,1)+c(:,i)
+c      write (iout,*) "alpha",rad2deg*alphi,
+c     & rad2deg*alpha(nres+i,i,2*nres+1)
+c      write (iout,*) "omega",rad2deg*omegi,
+c     & rad2deg*beta(nres+i,i,2*nres+1,i+1)
+      return
+      end
+c-----------------------------------------------------------------------------
       subroutine locate_side_chain(i)
 C 
 C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
@@ -296,7 +387,7 @@ cd   &   xp,yp,zp,(xx(k),k=1,3)
       xx(2)=xloc(2,i)*r(2,2,i-1)+xloc(3,i)*r(2,3,i-1)
       xx(3)=xloc(2,i)*r(3,2,i-1)+xloc(3,i)*r(3,3,i-1)
       do j=1,3
-       xrot(j,i)=xx(j)
+        xrot(j,i)=xx(j)
       enddo
       do j=1,3
         rj=0.0D0