C
C $Id: setbox.F,v 1.11 1998/07/16 16:40:40 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      subroutine setbox(ierror)

      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     local:
      character dir
      dimension xyztot(3), dir(3)
      data zero /0.0/
      data dir /'X','Y','Z'/
      save dir, zero

C     in the following "local" refers to residue i
C
C     xyztot -> local vector sum of all coordinates
C     xyzmin -> global minimum
C     xyzmax -> global maximum
C     gcmin  -> global minimum coordinate geometric center residues
C     gcmax  -> global maximum coordinate geometric center residues

CC-----------------------------------------------------CC

C     calculate geometric center and gyration radius of first residue

CCC   because center of mass is needed for virial calc, we
CCC   calculate COM instead of geom. center

      i1 = irpnt(1)
      iat = iatnum(i1)
      totm = emass(iat)
      do 10 j=1,3
         xyztot(j) = emass(iat)*xyz(j,i1)
         xyzmin(j) = xyz(j,i1)
         xyzmax(j) = xyz(j,i1)
 10   continue
      i2 = irpnt(2)-1
      do 30 i=i1+1,i2
         iat = iatnum(i)
         totm = totm + emass(iat)
         do 20 j=1,3
            xyztot(j) = xyztot(j) + emass(iat)*xyz(j,i)
            xyzmin(j) = min(xyzmin(j), xyz(j,i))
            xyzmax(j) = max(xyzmax(j), xyz(j,i))
 20      continue
 30   continue
ccc      nat = irpnt(2) - i1
      do 40 j=1,3
         gc(j,1) = xyztot(j)/totm
         gcmin(j) = gc(j,1)
         gcmax(j) = gc(j,1)
 40   continue
      do 60 i=i1,i2
         do 50 j=1,3
            gcres(j,i) = gc(j,1)
 50      continue
 60   continue

C     do the same for all the other residues

      do 130 ires=2,nres
         i1 = irpnt(ires)
         iat = iatnum(i1)
         totm = emass(iat)
         do 70 j=1,3
            xyztot(j) = emass(iat)*xyz(j,i1)
            xyzmin(j) = min(xyzmin(j),xyz(j,i1))
            xyzmax(j) = max(xyzmax(j),xyz(j,i1))
 70      continue
         i2 = irpnt(ires+1)-1
         do 90 i=i1+1,i2
            iat = iatnum(i)
            totm = totm + emass(iat)
            do 80 j=1,3
               xyztot(j) = xyztot(j) + emass(iat)*xyz(j,i)
               xyzmin(j) = min(xyzmin(j),xyz(j,i))
               xyzmax(j) = max(xyzmax(j),xyz(j,i))
 80         continue
 90      continue
ccc         nat = irpnt(ires+1) - i1
         do 100 j=1,3
            gc(j,ires) = xyztot(j)/totm
            gcmin(j) = min(gcmin(j), gc(j,ires))
            gcmax(j) = max(gcmax(j), gc(j,ires))
 100     continue
         do 120 i=i1,i2
            do 110 j=1,3
               gcres(j,i) = gc(j,ires)
 110        continue
 120     continue

 130  continue

      if (resgr) then
C     .   calculate maximum gyration radius
         gyr = 0.0
         do i=1,nres
            i1 = irpnt(i)
            i2 = irpnt(i+1)-1
            do ii=i1,i2
               do k=1,3
                  gyr = max(gyr,abs(xyz(k,ii)-gc(k,i)))
               enddo
            enddo
         enddo
      endif

C     check the PBC box, calculate boxmin, boxmax
      do 140 j=1,3
c        print *,  dxyz(j) , xyzmax(j) , xyzmin(j) 
         dxyz(j) = xyzmax(j) - xyzmin(j)
ccc         if (resgr.or.mcsim.or.(stand.and.pbc)) then
c        if (resgr.or.mcsim.or.pbc) then
            dd = gcmax(j)-gcmin(j)
            d = 0.5*(dbox(j) - dd)
            boxmin(j) = gcmin(j) - d
            boxmax(j) = gcmax(j) + d
c        else
c           dd = dxyz(j)
c           d = 0.5*(dbox(j) - dd)
c           boxmin(j) = xyzmin(j) - d
c           boxmax(j) = xyzmax(j) + d
c        endif
!         if (dd.gt.dbox(j)) then
C-RDC             write(iout,'(/" BOX ERROR: BOX DIMENSION IN ",
C-RDC      &           1A," DIRECTION EXCEEDED (",F9.4,")")')
C-RDC      &           dir(j),dd
!            ierror = 1
!            return
!         endif
 140  continue

!      boxvol = dbox(1)*dbox(2)*dbox(3)

      end