subroutine chainbuild
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,perbox,fail
C Set lprn=.true. for debugging
      lprn = .false.
      perbox=.false.
      fail=.false.
      print *, 'enter chainbuild' 
      call chainbuild_cart
      return
      end
#ifdef DEBUG
      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).
C
      call orig_frame
*
* Build the alpha-carbon chain.
*
      do i=4,nres
	call locate_next_res(i)
      enddo     
C
C First and last SC must coincide with the corresponding CA.
C
      do j=1,3
	dc(j,nres+1)=0.0D0
        dc_norm(j,nres+1)=0.0D0
	dc(j,nres+nres)=0.0D0
        dc_norm(j,nres+nres)=0.0D0
        c(j,nres+1)=c(j,1)
        c(j,nres+nres)=c(j,nres)
      enddo
*
* Temporary diagnosis
*
      if (lprn) then

      call cartprint
      write (iout,'(/a)') 'Recalculated internal coordinates'
      do i=2,nres-1
	do j=1,3
	  c(j,maxres2)=0.5D0*(c(j,i-1)+c(j,i+1))
        enddo
        be=0.0D0
        if (i.gt.3) be=rad2deg*beta(i-3,i-2,i-1,i)
        be1=rad2deg*beta(nres+i,i,maxres2,i+1)
        alfai=0.0D0
        if (i.gt.2) alfai=rad2deg*alpha(i-2,i-1,i)
        write (iout,1212) restyp(itype(i)),i,dist(i-1,i),
     &  alfai,be,dist(nres+i,i),rad2deg*alpha(nres+i,i,maxres2),be1
      enddo   
 1212 format (a3,'(',i3,')',2(f10.5,2f10.2))

      endif
      endif
      return
      end
#endif
c-------------------------------------------------------------------------
      subroutine orig_frame
C
C Define the origin and orientation of the coordinate system and locate 
C the first three atoms.
C
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.LOCAL'
      include 'COMMON.GEO'
      include 'COMMON.VAR'
      cost=dcos(theta(3))
      sint=dsin(theta(3))
      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,3,1)= 0.0D0
      t(3,1,1)= 0.0D0
      t(3,2,1)= 0.0D0
      t(3,3,1)= 1.0D0
      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   
      c(1,1)=0.0D0
      c(2,1)=0.0D0
      c(3,1)=0.0D0
      c(1,2)=vbld(2)
      c(2,2)=0.0D0
      c(3,2)=0.0D0
      dc(1,0)=0.0d0
      dc(2,0)=0.0D0
      dc(3,0)=0.0D0
      dc(1,1)=vbld(2)
      dc(2,1)=0.0D0
      dc(3,1)=0.0D0
      dc_norm(1,0)=0.0D0
      dc_norm(2,0)=0.0D0
      dc_norm(3,0)=0.0D0
      dc_norm(1,1)=1.0D0
      dc_norm(2,1)=0.0D0
      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)
      enddo
      call locate_side_chain(2)
      return
      end
c-----------------------------------------------------------------------------
      subroutine locate_next_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)
      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)
        c(j,i)=c(j,i-1)+dc(j,i-1)
      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_side_chain(i)
C 
C Locate the side-chain centroid i, 1 < i < NRES. Put in C(*,NRES+i).
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'
      dimension xx(3)

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) 
      xp= dsci*cosalphi
      yp= dsci*sinalphi*cosomegi
      zp=-dsci*sinalphi*sinomegi
* Now we have to rotate the coordinate system by 180-theta(i)/2 so as to get its
* X-axis aligned with the vector DC(*,i)
      theta2=pi-0.5D0*theta(i+1)
      cost2=dcos(theta2)
      sint2=dsin(theta2)
      xx(1)= xp*cost2+yp*sint2
      xx(2)=-xp*sint2+yp*cost2
      xx(3)= zp
cd    print '(a3,i3,3f10.5,5x,3f10.5)',restyp(itype(i)),i,
cd   &   xp,yp,zp,(xx(k),k=1,3)
      do j=1,3
        xloc(j,i)=xx(j)
      enddo
* Bring the SC vectors to the common coordinate system.
      xx(1)=xloc(1,i)
      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)
      enddo
      do j=1,3
        rj=0.0D0
        do k=1,3
          rj=rj+prod(j,k,i-1)*xx(k)
        enddo
        dc(j,nres+i)=rj
        dc_norm(j,nres+i)=rj*dsci_inv
        c(j,nres+i)=c(j,i)+rj
      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'
        j=1
        chain_beg=1
C        do i=1,nres
C       write(*,*) 'initial', i,j,c(j,i)
C        enddo
        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
           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
        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
           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
       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
           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