update
[unres.git] / source / unres / src_MD-M / chainbuild.F
index 45a1a53..1261b4a 100644 (file)
@@ -12,13 +12,37 @@ 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.
+      call chainbuild_cart
+      return
+      end
+C#ifdef DEBUG
+C      if (perbox) then
+C first three CAs and SC(2).
 C
-C Define the origin and orientation of the coordinate system and locate the
-C first three CA's and SC(2).
+      subroutine chainbuild_extconf
+C 
+C Build the virtual polypeptide chain. Side-chain centroids are moveable.
+C As of 2/17/95.
 C
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      double precision e1(3),e2(3),e3(3)
+      logical lprn /.false./,perbox,fail
+
+c      write (iout,*) "Calling chainbuild_extconf"
       call orig_frame
 *
 * Build the alpha-carbon chain.
@@ -43,6 +67,12 @@ C
       if (lprn) then
 
       call cartprint
+      write (iout,*) 'dc_norm'
+      do i=1,nres
+        write (iout,'(a3,1x,i3,3f10.5,5x,3f10.5)') 
+     &   restyp(itype(i)),i,(dc_norm(j,i),j=1,3),
+     &   (dc_norm(j,i+nres),j=1,3)
+      enddo
       write (iout,'(/a)') 'Recalculated internal coordinates'
       do i=2,nres-1
        do j=1,3
@@ -58,10 +88,11 @@ C
       enddo   
  1212 format (a3,'(',i3,')',2(f10.5,2f10.2))
 
+C      endif
       endif
-
       return
       end
+C#endif
 c-------------------------------------------------------------------------
       subroutine orig_frame
 C
@@ -126,13 +157,86 @@ 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
       end
 c-----------------------------------------------------------------------------
+      subroutine locate_prev_res(i)
+C
+C Locate CA(i) and SC(i-1)
+C
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.CHAIN'
+      include 'COMMON.LOCAL'
+      include 'COMMON.GEO'
+      include 'COMMON.VAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+C
+C Define the rotation matrices corresponding to CA(i)
+C
+#ifdef OSF
+      theti=theta(i+2)
+      if (theti.ne.theti) theti=100.0     
+      phii=phi(i+3)
+      if (phii.ne.phii) phii=180.0     
+#else
+      theti=theta(i+2)      
+      phii=phi(i+3)
+#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)=-prod(j,1,i+1)
+        dc(j,i)=-vbld(i)*prod(j,1,i+1)
+        c(j,i)=c(j,i+1)+dc(j,i)
+      enddo
+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)
+C
+      call locate_side_chain(i+1)
+      return
+      end
+c-----------------------------------------------------------------------------
       subroutine locate_next_res(i)
 C
 C Locate CA(i) and SC(i-1)
@@ -272,3 +376,255 @@ cd   &   xp,yp,zp,(xx(k),k=1,3)
       enddo
       return
       end
+c------------------------------------------
+      subroutine returnbox
+      include 'DIMENSIONS'
+#ifdef MPI
+      include 'mpif.h'
+#endif
+      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'
+      include 'COMMON.REMD'
+      include 'COMMON.SETUP'
+      include 'COMMON.MUCA'
+      include 'COMMON.HAIRPIN'
+C change suggested by Ana - begin
+      integer allareout
+C change suggested by Ana - end
+        j=1
+        chain_beg=1
+C        do i=1,nres
+C       write(*,*) 'initial', i,j,c(j,i)
+C        enddo
+C change suggested by Ana - begin
+        allareout=1
+C change suggested by Ana -end
+        do i=1,nres-1
+         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+          chain_end=i
+          if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxxsize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,chain_end
+              c(j,k)=c(j,k)-ireturnval*boxxsize
+              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
+            enddo
+C Suggested by Ana
+            if (chain_beg.eq.1) 
+     &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
+C Suggested by Ana -end
+           endif
+           chain_beg=i+1
+           allareout=1
+         else
+          if (int(c(j,i)/boxxsize).eq.0) allareout=0
+         endif
+        enddo
+         if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxxsize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,nres
+              c(j,k)=c(j,k)-ireturnval*boxxsize
+              c(j,k+nres)=c(j,k+nres)-ireturnval*boxxsize
+            enddo
+          endif
+C NO JUMP 
+C        do i=1,nres
+C        write(*,*) 'befor no jump', i,j,c(j,i)
+C        enddo
+        nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+C             print *,'diff', difference
+             if (difference.gt.boxxsize/2.0) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+              else
+              nojumpval=0
+              endif
+              endif
+              c(j,i)=c(j,i)+nojumpval*boxxsize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
+         enddo
+       nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+             if (difference.gt.boxxsize/2.0) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+              else
+              nojumpval=0
+              endif
+             endif
+              c(j,i)=c(j,i)+nojumpval*boxxsize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxxsize
+         enddo
+
+C        do i=1,nres
+C        write(*,*) 'after no jump', i,j,c(j,i)
+C        enddo
+
+C NOW Y dimension
+C suggesed by Ana begins
+        allareout=1
+C suggested by Ana ends
+        j=2
+        chain_beg=1
+        do i=1,nres-1
+         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+          chain_end=i
+          if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxysize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,chain_end
+              c(j,k)=c(j,k)-ireturnval*boxysize
+             c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
+            enddo
+C Suggested by Ana
+            if (chain_beg.eq.1)
+     &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
+C Suggested by Ana -end
+           endif
+           chain_beg=i+1
+           allareout=1
+         else
+          if (int(c(j,i)/boxysize).eq.0) allareout=0
+         endif
+        enddo
+         if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxysize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,nres
+              c(j,k)=c(j,k)-ireturnval*boxysize
+              c(j,k+nres)=c(j,k+nres)-ireturnval*boxysize
+            enddo
+          endif
+        nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+             if (difference.gt.boxysize/2.0) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+             else
+              nojumpval=0
+              endif
+           endif
+              c(j,i)=c(j,i)+nojumpval*boxysize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
+         enddo
+      nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+             if (difference.gt.boxysize/2.0) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+              else
+              nojumpval=0
+              endif
+            endif
+              c(j,i)=c(j,i)+nojumpval*boxysize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxysize
+         enddo
+C Now Z dimension
+C Suggested by Ana -begins
+        allareout=1
+C Suggested by Ana -ends
+       j=3
+        chain_beg=1
+        do i=1,nres-1
+         if ((itype(i).eq.ntyp1).and.(itype(i+1).eq.ntyp1)) then
+          chain_end=i
+          if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxysize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,chain_end
+              c(j,k)=c(j,k)-ireturnval*boxzsize
+              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
+            enddo
+C Suggested by Ana
+            if (chain_beg.eq.1)
+     &      dc_old(1,0)=dc_old(1,0)-ireturnval*boxxsize
+C Suggested by Ana -end
+           endif
+           chain_beg=i+1
+           allareout=1
+         else
+          if (int(c(j,i)/boxzsize).eq.0) allareout=0
+         endif
+        enddo
+         if (allareout.eq.1) then
+            ireturnval=int(c(j,i)/boxzsize)
+            if (c(j,i).le.0) ireturnval=ireturnval-1
+            do k=chain_beg,nres
+              c(j,k)=c(j,k)-ireturnval*boxzsize
+              c(j,k+nres)=c(j,k+nres)-ireturnval*boxzsize
+            enddo
+          endif
+        nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+             if (difference.gt.(boxzsize/2.0)) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+              else
+              nojumpval=0
+              endif
+            endif
+              c(j,i)=c(j,i)+nojumpval*boxzsize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
+         enddo
+       nojumpval=0
+        do i=2,nres
+           if (itype(i).eq.ntyp1 .and. itype(i-1).eq.ntyp1) then
+             difference=abs(c(j,i-1)-c(j,i))
+             if (difference.gt.boxzsize/2.0) then
+                if (c(j,i-1).gt.c(j,i)) then
+                  nojumpval=1
+                 else
+                   nojumpval=-1
+                 endif
+              else
+              nojumpval=0
+              endif
+            endif
+             c(j,i)=c(j,i)+nojumpval*boxzsize
+              c(j,i+nres)=c(j,i+nres)+nojumpval*boxzsize
+         enddo
+
+        return
+        end
+