C
C $Id: boxmove.F,v 1.7 1998/07/16 16:39:25 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
      subroutine boxmove(ires,movedres,subpicked1)

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

C     does the MC box move for a NPT ensemble

      dimension ires(maxres)
      logical subpicked1(maxsub)

CCC-------------------------------------------------------------CCC

C     move the box (same move in all directions)

      sdrbox = drboxmc*(2.0*xrandom(iseed)-1.0)
      sdrbox2 = 0.5*sdrbox

      do 210 i=1,3
         oldbox = dbox(i)
         oldmin = boxmin(i)

         dbox(i) = oldbox + sdrbox
         dhalf(i) = 0.5*dbox(i)
         boxmax(i) = boxmax(i) + sdrbox2
         boxmin(i) = oldmin - sdrbox2

         xyzmin(i) = boxmax(i)
         xyzmax(i) = boxmin(i)
         gcmin(i) = boxmax(i)
         gcmax(i) = boxmin(i)
C     .  scale the old geometric center
         do 150 ir=1,nres
            oldgc = gcmc(i,ir)
            gc(i,ir) = ((oldgc-oldmin)*dbox(i))/oldbox+boxmin(i)
            gcmax(i) = max(gcmax(i), gc(i,ir))
            gcmin(i) = min(gcmin(i), gc(i,ir))
            dgc = gc(i,ir) - oldgc
C     .     move the atoms
            jbeg = irpnt(ir)
            jend = irpnt(ir+1)-1
            do 140 j=jbeg,jend
               xyz(i,j) = xyzmc(i,j) + dgc
               gcres(i,j) = gc(i,ir)
               xyzmin(i) = min(xyz(i,j),xyzmin(i))
               xyzmax(i) = max(xyz(i,j),xyzmax(i))
 140        continue
 150     continue
 210  continue

      do i=1,3
         dxyz(i) = xyzmax(i)-xyzmin(i)
      enddo
      boxvol = dbox(1)*dbox(2)*dbox(3)

      movedres = nres
      do 1000 i=1,nres
         ires(i) = i
 1000 continue

C     all subsystems moved, so always do a full SCF iteration
      if (frozenmc.and.(.not.stand)) then
         do i=1,nsub
            notfroz(i) = .true.
            subpicked1(i) = .true.
         enddo
      endif

      recip1(1) = 1.0/dbox(1)
      recip2(2) = 1.0/dbox(2)
      recip3(3) = 1.0/dbox(3)

      dens = dble(nres)/boxvol
      denstemp = dens*tempmc

      if (setbpme) then
         coulco = min(dhalf(1),dhalf(2),dhalf(3))
         coulco = coulco-1.0E-3
         betapme = 3.5/coulco
         betapme2 = betapme*betapme
      endif	

      end