X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fchainbuild.F;h=c4970d37dc1370f0977da443354c7190502afb77;hb=0bb81c1c3180a2079d70af7dd534295e1e0b1e4c;hp=45a1a53837e28da1c8198f61d4d9973d13076818;hpb=2ec03757b68a8fb1a604a7f84f83086291fc75db;p=unres.git

diff --git a/source/unres/src_MD-M/chainbuild.F b/source/unres/src_MD-M/chainbuild.F
index 45a1a53..c4970d3 100644
--- a/source/unres/src_MD-M/chainbuild.F
+++ b/source/unres/src_MD-M/chainbuild.F
@@ -12,9 +12,121 @@ C
       include 'COMMON.IOUNITS'
       include 'COMMON.NAMES'
       include 'COMMON.INTERACT'
-      logical lprn
+      double precision e1(3),e2(3),e3(3)
+      logical lprn,perbox,fail
 C Set lprn=.true. for debugging
       lprn = .false.
+      perbox=.false.
+      fail=.false.
+       if (perbox) then
+      cost=dcos(theta(3))
+      sint=dsin(theta(3))
+      print *,'before refsys'
+      call refsys(2,3,4,e1,e2,e3,fail)
+      print *,'after refsys'
+          if (fail) then
+            e2(1)=0.0d0
+            e2(2)=1.0d0
+            e2(3)=0.0d0
+          endif
+      dc(1,0)=c(1,1)
+      dc(2,0)=c(2,1)
+      dc(3,0)=c(3,1)
+      print *,'dc',dc(1,0),dc(2,0),dc(3,0)
+      dc(1,1)=c(1,2)-c(1,1)
+      dc(2,1)=c(2,2)-c(2,1)
+      dc(3,1)=c(3,2)-c(3,1)
+      dc(1,2)=c(1,3)-c(1,2)
+      dc(2,2)=c(2,3)-c(2,2)
+      dc(3,2)=c(3,3)-c(3,2)
+      t(1,1,1)=e1(1)
+      t(1,2,1)=e1(2)
+      t(1,3,1)=e1(3)
+      t(2,1,1)=e2(1)
+      t(2,2,1)=e2(2)
+      t(2,3,1)=e2(3)
+      t(3,1,1)=e3(1)
+      t(3,2,1)=e3(2)
+      t(3,3,1)=e3(3)
+      veclen=0.0d0
+      do i=1,3
+       veclen=veclen+(c(i,2)-c(i,1))**2
+      enddo
+      veclen=sqrt(veclen)
+      r(1,1,1)= 1.0D0
+      r(1,2,1)= 0.0D0
+      r(1,3,1)= 0.0D0
+      r(2,1,1)= 0.0D0
+      r(2,2,1)= 1.0D0
+      r(2,3,1)= 0.0D0
+      r(3,1,1)= 0.0D0
+      r(3,2,1)= 0.0D0
+      r(3,3,1)= 1.0D0
+      do i=1,3
+        do j=1,3
+          rt(i,j,1)=t(i,j,1)
+        enddo
+      enddo
+      do i=1,3
+        do j=1,3
+          prod(i,j,1)=0.0D0
+          prod(i,j,2)=t(i,j,1)
+        enddo
+        prod(i,i,1)=1.0D0
+      enddo
+        call locate_side_chain(2)
+      do i=4,nres
+#ifdef OSF
+      theti=theta(i)
+      if (theti.ne.theti) theti=100.0
+      phii=phi(i)
+      if (phii.ne.phii) phii=180.0
+#else
+      theti=theta(i)
+      phii=phi(i)
+#endif
+      cost=dcos(theti)
+      sint=dsin(theti)
+      cosphi=dcos(phii)
+      sinphi=dsin(phii)
+* 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,3,i-2)= 0.0D0
+      t(2,1,i-2)=-sint
+      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
+      t(3,3,i-2)= 1.0D0
+      r(1,1,i-2)= 1.0D0
+      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(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,3,i-2)=0.0D0
+      rt(2,1,i-2)=sint*cosphi
+      rt(2,2,i-2)=-cost*cosphi
+      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
+        call matmult(prod(1,1,i-2),rt(1,1,i-2),prod(1,1,i-1))
+      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)
+      enddo
+        call locate_side_chain(i-1)
+       enddo
+      else
 C
 C Define the origin and orientation of the coordinate system and locate the
 C first three CA's and SC(2).
@@ -59,7 +171,7 @@ C
  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 
       endif
-
+      endif
       return
       end
 c-------------------------------------------------------------------------
@@ -126,8 +238,8 @@ C
       dc_norm(3,1)=0.0D0
       do j=1,3
         dc_norm(j,2)=prod(j,1,2)
-	dc(j,2)=vbld(3)*prod(j,1,2)
-	c(j,3)=c(j,2)+dc(j,2)
+        dc(j,2)=vbld(3)*prod(j,1,2)
+        c(j,3)=c(j,2)+dc(j,2)
       enddo
       call locate_side_chain(2)
       return