********>Bugfix 7:
Authors: Viktor Hornak
Date: 05/03/2004

Programs: sander

Description: The restraintmask or bellymask strings are parsed
             before the distances are read; hence the options to generate
             subsets of atoms based on distances will not work.

Fix:  copy this file to "$AMBERHOME/src/patchfile".  Then:
         cd $AMBERHOME/src
         patch -p0 < patchfile

------------------------------------------------------------------------------
*** sander/findmask.f	2004-01-28 20:34:03.000000000 -0500
--- sander.fix/findmask.f	2004-04-26 11:35:20.000000000 -0400
***************
*** 707,713 ****
        call error1("noneq", "error parsing distance cutoff value")
     end if
  
!    if (diststr(1:1).eq.'@'.or.diststr(1:1).eq.':') then   ! atom based cutoff
        do i=1,natom
           if (mask2(i) == 1) then
              ii = 3*i - 3
--- 707,716 ----
        call error1("noneq", "error parsing distance cutoff value")
     end if
  
!    ! First, find atoms closer than 'dist' (regardless of whether the
!    ! selection is atom or residue based and regardless of the distance
!    ! operator ('<' or '>') used).
!    if (diststr(1:1).eq.'@'.or.diststr(1:1).eq.':') then
        do i=1,natom
           if (mask2(i) == 1) then
              ii = 3*i - 3
***************
*** 716,733 ****
                 r2 = (crd(ii+1)-crd(jj+1))**2 +  &
                      (crd(ii+2)-crd(jj+2))**2 +  &
                      (crd(ii+3)-crd(jj+3))**2
!                if (op.eq.'<') then
!                   if (r2 < d2) mask(j) = 1
!                else if (op.eq.'>') then
!                   if (r2 > d2) mask(j) = 1
!                else
!                   call error1("noneq","operator is not < or >")
!                end if
              end do ! j
           end if
        end do ! i
     end if
!    
     if (diststr(1:1).eq.':') then    ! residue based cutoff
        ! if at least one atom in a residue is selected (from
        ! the previous step), select the whole residue
--- 719,730 ----
                 r2 = (crd(ii+1)-crd(jj+1))**2 +  &
                      (crd(ii+2)-crd(jj+2))**2 +  &
                      (crd(ii+3)-crd(jj+3))**2
!                if (r2 < d2) mask(j) = 1
              end do ! j
           end if
        end do ! i
     end if
! 
     if (diststr(1:1).eq.':') then    ! residue based cutoff
        ! if at least one atom in a residue is selected (from
        ! the previous step), select the whole residue
***************
*** 743,748 ****
--- 740,762 ----
        end do
     end if
     
+    if (op.eq.'<') then
+       continue  ! you're done
+    else if (op.eq.'>') then
+       ! Invert the selection: the best way of thinking about it is that
+       ! if you want atoms/residues further away than 'dist', select
+       ! atoms/residues that are closer than 'dist' and neg this selection.
+       do i=1,natom
+          if (mask(i) == 1) then
+             mask(i) = 0
+          else
+             mask(i) = 1
+          end if
+       end do ! i
+    else
+       call error1("noneq","operator is not < or >")
+    end if
+    
     if (diststr(1:1).ne.'@'.and.diststr(1:1).ne.':') then
        call error1("noneq", "residue (<:,>:) or atom (<@,>@) " //  &
                    "based distance cutoff must be specified")
===============================================================================
*** sander/getcor.f	2004-03-04 11:39:03.000000000 -0500
--- sander.fix/getcor.f	2004-04-26 10:04:00.000000000 -0400
***************
*** 4,12 ****
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !+ [Enter a one-line description of subroutine getcor here]
  #ifdef QMMM
! subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0,nlink)
  #else
! subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0)
  #endif
     
     !     --- reads initial coords, vel, and box for MD.
--- 4,12 ----
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  !+ [Enter a one-line description of subroutine getcor here]
  #ifdef QMMM
! subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0,writeflag,nlink)
  #else
! subroutine getcor(nr,x,v,f,ntx,box,irest,tt,temp0,writeflag)
  #endif
     
     !     --- reads initial coords, vel, and box for MD.
***************
*** 18,23 ****
--- 18,27 ----
     !     EWALD: read a modified form of the coordinate/restart file
     
     !      IMPLICIT _REAL_ (A-H,O-Z)
+ 
+    ! writeflag was introduced to avoid false failures in test.sander
+    ! after atommask() dependency on getcor() was fixed by calling getcor()
+    ! twice: once before atommask() calls and second time in its usual place.
     
     implicit none
  #  include "files.h"
***************
*** 31,37 ****
     integer lun,nr,ntx,irest
     _REAL_ x(*),v(*),f(*),box(3),tt,temp0,ltemp0
     
!    logical form
     _REAL_ fvar(7),a,b,c,alpha,beta,gamma
     integer i,ivar,ifld(7),ihol(1)
     integer natom,nr3,ier
--- 35,41 ----
     integer lun,nr,ntx,irest
     _REAL_ x(*),v(*),f(*),box(3),tt,temp0,ltemp0
     
!    logical form, writeflag
     _REAL_ fvar(7),a,b,c,alpha,beta,gamma
     integer i,ivar,ifld(7),ihol(1)
     integer natom,nr3,ier
***************
*** 41,48 ****
     ltemp0 = 0.0d0
     nr3 = 3*nr
     lun = 9 !FIXME: use lun=-1 for new amopen()
!    
!    write(6,9108)
     
     !     ----- OPEN THE COORDINATE FILE -----
     
--- 45,52 ----
     ltemp0 = 0.0d0
     nr3 = 3*nr
     lun = 9 !FIXME: use lun=-1 for new amopen()
! 
!    if (writeflag) write(6,9108)
     
     !     ----- OPEN THE COORDINATE FILE -----
     
***************
*** 95,102 ****
           do i = 1,nr3
              v(i) = 0.d0
           end do
!          write(6,9008) title1
!          write(6,9009) tt
           close(lun, iostat=ier)
           return
        end if
--- 99,108 ----
           do i = 1,nr3
              v(i) = 0.d0
           end do
!          if (writeflag) then
!             write(6,9008) title1
!             write(6,9009) tt
!          end if
           close(lun, iostat=ier)
           return
        end if
***************
*** 108,115 ****
  #else
        read(lun,9028,end=1010) (v(i),i=1,nr3)
  #endif
!       write(6,9008) title1
!       write(6,9009) tt
        close(lun, iostat=ier)
        return
        
--- 114,123 ----
  #else
        read(lun,9028,end=1010) (v(i),i=1,nr3)
  #endif
!       if (writeflag) then
!          write(6,9008) title1
!          write(6,9009) tt
!       end if
        close(lun, iostat=ier)
        return
        
===============================================================================
*** sander/mdread.f	2004-03-04 11:39:04.000000000 -0500
--- sander.fix/mdread.f	2004-04-26 10:29:26.000000000 -0400
***************
*** 1581,1586 ****
--- 1581,1599 ----
  
        call rdrest(natom,ntx,ntrx,x(lcrdr))
        close(10)
+ 
+       ! inserted here to fix the bug that coords are not available
+       ! yet when distance based selection (<,>) is requested
+ #ifdef QMMM
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0, &
+                   .FALSE.,0)
+ #else
+ #ifdef LES
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.FALSE.)
+ #else
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.)
+ #endif
+ #endif
        
        ! VH - tgtmd change: preferably call atommask() instead of rgroup()
        if (konst) then
***************
*** 1673,1678 ****
--- 1686,1703 ----
     belly = ibelly > 0
     ngrp = 0
     if(belly) then
+       ! inserted here to fix the bug that coords are not available
+       ! yet when distance based selection (<,>) is requested
+ #ifdef QMMM
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0, &
+                   .FALSE.,0)
+ #else
+ #ifdef LES
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.FALSE.)
+ #else
+       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.FALSE.)
+ #endif
+ #endif
        write(6,9418)
        if( len_trim(bellymask) <= 0 ) then
           call rgroup(natom,natbel,nres,ngrp,ix(i02),ih(m02), &
===============================================================================
*** sander/sander.f	2004-03-12 17:03:31.000000000 -0500
--- sander.fix/sander.f	2004-04-26 10:15:45.000000000 -0400
***************
*** 780,791 ****
  #ifdef QMMM
        nlink = 0
        call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0, &
!                  nlink)
  #else
  #ifdef LES
!       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les)
  #else
!       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0)
  #endif
  #endif
  #ifdef MPI
--- 780,791 ----
  #ifdef QMMM
        nlink = 0
        call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0, &
!                   .TRUE.,nlink)
  #else
  #ifdef LES
!       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0les,.TRUE.)
  #else
!       call getcor(nr,x(lcrd),x(lvel),x(lforce),ntx,box,irest,t,temp0,.TRUE.)
  #endif
  #endif
  #ifdef MPI
------------------------------------------------------------------------------

Temporary workarounds: don't use the distance features of the mask routines