********>Bugfix 37:
Author: Ross Walker & Mike Crowley
Date: 08/06/2007

Programs: sander

Description:  The routine that checks whether the QM region + the cutoff
              lies within the periodic box could generate an improper
              error message due to the fact that it did not carry out
              the check on the imaged coordinates. Hence very long QMMM
              runs might stop due to the QM region diffusing out of the 
              central box. This bugfix addresses this issue.

Fix: Apply the following patches to src/sander/qm_mm.f

------------------------------------------------------------------------------
--- qm_mm.f	2007-08-06 13:57:27.000000000 -0700
+++ qm_mm.f	2007-08-06 13:58:24.000000000 -0700
@@ -538,6 +538,7 @@
 subroutine qm_check_periodic(a,b,c,alpha,beta,gamma,x,natom,bxbnd,iqmatoms)
 
    use qmmm_module, only: qmmm_nml,qmmm_struct
+   use constants, only : zero, one
 
    implicit none
    integer , intent(in) :: natom
@@ -546,8 +547,10 @@
    _REAL_ , intent(out), dimension(2,3) :: bxbnd
    integer, intent(in), dimension(qmmm_struct%nquant) :: iqmatoms
 
-   integer j,nquant
+   integer j,iqm_one, jqmatom, ier
+   _REAL_, allocatable, dimension(:,:) :: tmpcrd
    _REAL_ :: xbnd0,xbnd1,ybnd0,ybnd1,zbnd0,zbnd1
+   _REAL_ :: bxinv(3), offset(3), frac(3), xx,yy,zz
 
    if( (alpha /= 90.d0) .or.  &
        (beta  /= 90.d0) .or.  &
@@ -560,33 +563,60 @@
               "QMMM periodic error, need 90 degree angles", &
               "Exiting")
    endif
-!  ---mfc 
-!     Find the bounding box limits of the QM regoin, then find all atoms
-!     inside the box + cutoff before calculating or testing distances.
-   nquant=qmmm_struct%nquant
-   xbnd0=x(1,iqmatoms(1))
-   xbnd1=x(1,iqmatoms(1))
-   ybnd0=x(2,iqmatoms(1))
-   ybnd1=x(2,iqmatoms(1))
-   zbnd0=x(3,iqmatoms(1))
-   zbnd1=x(3,iqmatoms(1))
-   
-   do j=2,nquant
-      xbnd0=min(xbnd0,x(1,iqmatoms(j)))
-      xbnd1=max(xbnd1,x(1,iqmatoms(j)))
-      ybnd0=min(ybnd0,x(2,iqmatoms(j)))
-      ybnd1=max(ybnd1,x(2,iqmatoms(j)))
-      zbnd0=min(zbnd0,x(3,iqmatoms(j)))
-      zbnd1=max(zbnd1,x(3,iqmatoms(j)))
-   enddo
 
+   allocate(tmpcrd(3,qmmm_struct%nquant),stat=ier)
+   REQUIRE(ier==0)
+   bxinv(1)=one/a
+   bxinv(2)=one/b
+   bxinv(3)=one/c
+
+   tmpcrd(1:3,1) = zero
+
+!Move the first QM atom to be at the origin.
+   iqm_one = iqmatoms(1)
+
+   frac(1)=x(1,iqm_one)*bxinv(1)-anint(x(1,iqm_one)*bxinv(1))
+   frac(2)=x(2,iqm_one)*bxinv(2)-anint(x(2,iqm_one)*bxinv(2))
+   frac(3)=x(3,iqm_one)*bxinv(3)-anint(x(3,iqm_one)*bxinv(3))
+
+   offset(1)=a*frac(1)
+   offset(2)=b*frac(2)
+   offset(3)=c*frac(3)
+
+   do j = 1, qmmm_struct%nquant
+     jqmatom = iqmatoms(j)
+     xx = x(1,jqmatom) - offset(1)
+     yy = x(2,jqmatom) - offset(2)
+     zz = x(3,jqmatom) - offset(3)
+     tmpcrd(1,j)=a*(xx*bxinv(1)-anint(xx*bxinv(1)))
+     tmpcrd(2,j)=b*(yy*bxinv(2)-anint(yy*bxinv(2)))
+     tmpcrd(3,j)=c*(zz*bxinv(3)-anint(zz*bxinv(3)))
+   end do
+
+   xbnd0=zero; xbnd1=zero; ybnd0=zero; ybnd1=zero; zbnd0=zero; zbnd1=zero
+
+   do j=1,qmmm_struct%nquant
+      xbnd0=min(xbnd0,tmpcrd(1,j))
+      xbnd1=max(xbnd1,tmpcrd(1,j))
+      ybnd0=min(ybnd0,tmpcrd(2,j))
+      ybnd1=max(ybnd1,tmpcrd(2,j))
+      zbnd0=min(zbnd0,tmpcrd(3,j))
+      zbnd1=max(zbnd1,tmpcrd(3,j))
+   enddo
    xbnd0=xbnd0-qmmm_nml%qmcut
    ybnd0=ybnd0-qmmm_nml%qmcut
    zbnd0=zbnd0-qmmm_nml%qmcut
    xbnd1=xbnd1+qmmm_nml%qmcut
    ybnd1=ybnd1+qmmm_nml%qmcut
    zbnd1=zbnd1+qmmm_nml%qmcut
-   
+
+   bxbnd(1,1)=xbnd0+offset(1)
+   bxbnd(2,1)=xbnd1+offset(1)
+   bxbnd(1,2)=ybnd0+offset(2)
+   bxbnd(2,2)=ybnd1+offset(2)
+   bxbnd(1,3)=zbnd0+offset(3)
+   bxbnd(2,3)=zbnd1+offset(3)
+
   !------ Check if QM region plus cutoff around it is too large for this box
    if( (xbnd1-xbnd0 > a) .or. (ybnd1-ybnd0 > b) .or. (zbnd1-zbnd0 > c) )then
       write(6,*) " ****************************************************"
@@ -601,20 +631,10 @@
         "QM region + cutoff larger than box", &
         "cannot continue, need larger box.")
    endif
-  !------ Check if QM region plus cutoff box hangs out of the box
-!   if( ( (xbnd1 > a   ) .or. (ybnd1 > b   ) .or. (zbnd1 > c   ) ) .or. &
-!       ( (xbnd0 < zero) .or. (ybnd0 < zero) .or. (zbnd0 < zero) ) &
-!      )then
-!      call sander_bomb("QM_CHECK_PERIODIC<qm_mm.f>", &
-!        "QM regoin + cutoff extends out of box", &
-!        "cannot continue, need to center QM region in box.")
-!   endif
-   bxbnd(1,1)=xbnd0
-   bxbnd(2,1)=xbnd1
-   bxbnd(1,2)=ybnd0
-   bxbnd(2,2)=ybnd1
-   bxbnd(1,3)=zbnd0
-   bxbnd(2,3)=zbnd1
+
+   deallocate(tmpcrd,stat=ier)
+   REQUIRE(ier==0)
+
    return
 end subroutine qm_check_periodic

------------------------------------------------------------------------------
Temporary Workarounds: Reimage the restart file before each QMMM MD restart.