********>Bugfix 46:
Author: Bob Duke (bug reported and analyzed by Petr Kulhanek)
Date: 07/8/2005
 
Programs: pmemd
 
Description: The pmemd indexing algorithm assumes that fractional coordinates
             are in the range 0 <= f < 1.0.  But it is possible to obtain a
             fractional coordinate that is exactly 1.0, which can lead to
             erroneous indexing.  This should be a very rare error.

Fix:  Apply the following patches in $AMBERHOME/src/pmemd/src to
      ew_recip_cit.f90, ew_pairlist.f90 and ew_direct_cit.f90
 
------------------------------------------------------------------------------
*** ew_recip_cit.f90	Mon Feb 16 17:31:07 2004
--- ew_recip_cit.f90	Fri Jul  8 15:52:15 2005
***************
*** 1488,1493 ****
--- 1488,1499 ----
    z_bkt_lo = int(z_lo * dble(cit_tbl_z_dim) / z_box)
    z_bkt_hi = int(z_hi * dble(cit_tbl_z_dim) / z_box)
  
+   ! An extra check to insure that rounding error above does not cause us to
+   ! index off one end or the other of the cit table
+ 
+   if (z_bkt_lo .lt. 0) z_bkt_lo = 0
+   if (z_bkt_hi .ge. cit_tbl_z_dim) z_bkt_hi = cit_tbl_z_dim - 1
+ 
    if (z_bkt_lo .eq. z_bkt_hi) z_wraps = .false.
  
    if (z_wraps) then
***************
*** 1760,1765 ****
--- 1766,1777 ----
    z_bkt_lo = int(z_lo * dble(cit_tbl_z_dim))
    z_bkt_hi = int(z_hi * dble(cit_tbl_z_dim))
  
+   ! An extra check to insure that rounding error above does not cause us to
+   ! index off one end or the other of the cit table
+ 
+   if (z_bkt_lo .lt. 0) z_bkt_lo = 0
+   if (z_bkt_hi .ge. cit_tbl_z_dim) z_bkt_hi = cit_tbl_z_dim - 1
+ 
    if (z_bkt_lo .eq. z_bkt_hi) z_wraps = .false.
  
    if (z_wraps) then


*** ew_pairlist.f90	Mon Feb 16 17:31:07 2004
--- ew_pairlist.f90	Fri Jul  8 15:52:13 2005
***************
*** 30,40 ****
                         excl_img_map, img_maskdata, img_mask, tranvec, &
                         atm_maskdata, atm_mask, atm_img_map, &
                         dont_skip_belly_pairs, &
- #ifdef MPI
                         fraction, ifail)
- #else
-                        ifail)
- #endif
  
  #include "use_ew_dat.h"
  #include "use_ew_direct_cit_dat.h"
--- 30,36 ----
***************
*** 65,73 ****
    integer               :: atm_mask(*)
    integer               :: atm_img_map(atm_cnt)
    logical               :: dont_skip_belly_pairs
- #ifdef MPI
    double precision      :: fraction(3, atm_cnt)
- #endif
    integer               :: ifail
  
  ! Local record types:
--- 61,67 ----
***************
*** 82,90 ****
  
  ! Local variables:
  
!   double precision      :: cutlist_stk, cutlist_sq
    double precision      :: x_box, y_box, z_box
!   double precision      :: x_scale, y_scale, z_scale
    double precision      :: x_i, y_i, z_i
    double precision      :: x_j, y_j, z_j
    double precision      :: dx, dy, dz
--- 76,85 ----
  
  ! Local variables:
  
!   double precision      :: cutlist_sq
    double precision      :: x_box, y_box, z_box
!   double precision      :: scale_fac_x, scale_fac_y, scale_fac_z
!   double precision      :: cut_x_frac, cut_y_frac, cut_z_frac
    double precision      :: x_i, y_i, z_i
    double precision      :: x_j, y_j, z_j
    double precision      :: dx, dy, dz
***************
*** 94,99 ****
--- 89,95 ----
    integer               :: atm_mask_idx
    integer               :: nxt_img_maskptr, nummask
    integer               :: num_packed
+   integer               :: x_i_idx, y_i_idx, z_i_idx
    integer               :: x_bkts_lo, x_bkts_hi
    integer               :: y_bkts_lo, y_bkts_hi
    integer               :: z_bkts_lo, z_bkts_hi
***************
*** 125,141 ****
    ntypes_stk = ntypes
  #endif
  
    x_box = box(1)
    y_box = box(2)
    z_box = box(3)
  
!   x_scale = dble(cit_tbl_x_dim) / x_box
!   y_scale = dble(cit_tbl_y_dim) / y_box
!   z_scale = dble(cit_tbl_z_dim) / z_box
! 
!   cutlist_stk = cutlist
! 
!   cutlist_sq = cutlist_stk * cutlist_stk
  
    ! Set up the bucket and translation index mapping arrays.  We really only need
    ! to do this each time if the cit table dimensions are changing, but there
--- 121,139 ----
    ntypes_stk = ntypes
  #endif
  
+   cutlist_sq = cutlist * cutlist
+ 
    x_box = box(1)
    y_box = box(2)
    z_box = box(3)
  
!   scale_fac_x = dble(cit_tbl_x_dim)
!   scale_fac_y = dble(cit_tbl_y_dim)
!   scale_fac_z = dble(cit_tbl_z_dim)
! 
!   cut_x_frac = cutlist / x_box + 0.0001d0 ! fudge to thwart rounding errs
!   cut_y_frac = cutlist / y_box + 0.0001d0 ! fudge to thwart rounding errs
!   cut_z_frac = cutlist / z_box + 0.0001d0 ! fudge to thwart rounding errs
  
    ! Set up the bucket and translation index mapping arrays.  We really only need
    ! to do this each time if the cit table dimensions are changing, but there
***************
*** 179,192 ****
      iaci = ntypes_stk * (img(img_i)%iac - 1)
  #endif
  
!     ! Determine the bucket ranges that need to be searched:
  
      x_i = img(img_i)%imgcrd(1)
      y_i = img(img_i)%imgcrd(2)
      z_i = img(img_i)%imgcrd(3)
  
-     atm_i = img_atm_map(img_i)
- 
      do i = 0, 17
        i_tranvec(1, i) = tranvec(1, i) - x_i 
        i_tranvec(2, i) = tranvec(2, i) - y_i 
--- 177,190 ----
      iaci = ntypes_stk * (img(img_i)%iac - 1)
  #endif
  
!     atm_i = img_atm_map(img_i)
! 
!     ! These are imaged coordinates:
  
      x_i = img(img_i)%imgcrd(1)
      y_i = img(img_i)%imgcrd(2)
      z_i = img(img_i)%imgcrd(3)
  
      do i = 0, 17
        i_tranvec(1, i) = tranvec(1, i) - x_i 
        i_tranvec(2, i) = tranvec(2, i) - y_i 
***************
*** 209,227 ****
        excl_img_map(img_j) = 1                   ! .ne. 0 = excluded
      end do
  
!     x_bkts_lo = int((x_i - cutlist_stk + x_box) * x_scale)
!     x_bkts_hi = int((x_i + cutlist_stk + x_box) * x_scale)
  
!     y_bkts_lo = int((y_i - cutlist_stk + y_box) * y_scale)
!     y_bkts_hi = int((y_i + cutlist_stk + y_box) * y_scale)
  
!     z_bkts_lo = int((z_i - cutlist + z_box) * z_scale)
!     z_bkts_hi = int((z_i + z_box) * z_scale)
  
-     i_bkt = x_bkts(int((x_i + x_box) * x_scale)) + &
-             y_bkts(int((y_i + y_box) * y_scale)) + &
-             z_bkts(z_bkts_hi)
-                                         
      saved_i_bkt_img_j_hi = flat_cit(i_bkt)%img_hi
      flat_cit(i_bkt)%img_hi = img_i - 1
  
--- 207,231 ----
        excl_img_map(img_j) = 1                   ! .ne. 0 = excluded
      end do
  
!     ! These are indexes into the 3D cit:
  
!     x_i_idx = int(fraction(1, atm_i) * scale_fac_x)
!     y_i_idx = int(fraction(2, atm_i) * scale_fac_y)
!     z_i_idx = int(fraction(3, atm_i) * scale_fac_z)
  
!     ! This is the flat cit bucket index:
! 
!     i_bkt = x_i_idx + cit_tbl_x_dim * (y_i_idx + z_i_idx * cit_tbl_y_dim)
! 
!     ! Determine the bucket ranges that need to be searched:
! 
!     x_bkts_lo = int((1.d0 + fraction(1, atm_i) - cut_x_frac) * scale_fac_x)
!     x_bkts_hi = int((1.d0 + fraction(1, atm_i) + cut_x_frac) * scale_fac_x)
!     y_bkts_lo = int((1.d0 + fraction(2, atm_i) - cut_y_frac) * scale_fac_y)
!     y_bkts_hi = int((1.d0 + fraction(2, atm_i) + cut_y_frac) * scale_fac_y)
!     z_bkts_lo = int((1.d0 + fraction(3, atm_i) - cut_z_frac) * scale_fac_z)
!     z_bkts_hi = z_i_idx + cit_tbl_z_dim
  
      saved_i_bkt_img_j_hi = flat_cit(i_bkt)%img_hi
      flat_cit(i_bkt)%img_hi = img_i - 1
  
***************
*** 436,446 ****
  ! Local variables:
  
    double precision      :: cutlist_sq
!   double precision      :: dble_x_tbl_dim, dble_y_tbl_dim, dble_z_tbl_dim
    double precision      :: x_i, y_i, z_i
    double precision      :: x_j, y_j, z_j
    double precision      :: x_fract_i, y_fract_i, z_fract_i
-   double precision      :: x_fract_cut, y_fract_cut, z_fract_cut
    double precision      :: f1, f2, f3
    double precision      :: ucell_stk(3, 3)
    double precision      :: dx, dy, dz
--- 440,450 ----
  ! Local variables:
  
    double precision      :: cutlist_sq
!   double precision      :: scale_fac_x, scale_fac_y, scale_fac_z
!   double precision      :: cut_x_frac, cut_y_frac, cut_z_frac
    double precision      :: x_i, y_i, z_i
    double precision      :: x_j, y_j, z_j
    double precision      :: x_fract_i, y_fract_i, z_fract_i
    double precision      :: f1, f2, f3
    double precision      :: ucell_stk(3, 3)
    double precision      :: dx, dy, dz
***************
*** 451,456 ****
--- 455,461 ----
    integer               :: atm_mask_idx
    integer               :: nxt_img_maskptr, nummask
    integer               :: num_packed
+   integer               :: x_i_idx, y_i_idx, z_i_idx
    integer               :: x_bkts_lo, x_bkts_hi
    integer               :: y_bkts_lo, y_bkts_hi
    integer               :: z_bkts_lo, z_bkts_hi
***************
*** 482,498 ****
    ntypes_stk = ntypes
  #endif
  
!   dble_x_tbl_dim = dble(cit_tbl_x_dim)
!   dble_y_tbl_dim = dble(cit_tbl_y_dim)
!   dble_z_tbl_dim = dble(cit_tbl_z_dim)
  
    ucell_stk(:,:) = ucell(:,:)
  
!   x_fract_cut = cut_factor(1) * (cutlist / ucell_stk(1,1))
!   y_fract_cut = cut_factor(2) * (cutlist / ucell_stk(2,2))
!   z_fract_cut = cut_factor(3) * (cutlist / ucell_stk(3,3))
! 
!   cutlist_sq = cutlist * cutlist
  
    ! Set up the bucket and translation index mapping arrays.  We really only need
    ! to do this each time if the cit table dimensions are changing, but there
--- 487,505 ----
    ntypes_stk = ntypes
  #endif
  
!   cutlist_sq = cutlist * cutlist
  
    ucell_stk(:,:) = ucell(:,:)
  
!   scale_fac_x = dble(cit_tbl_x_dim)
!   scale_fac_y = dble(cit_tbl_y_dim)
!   scale_fac_z = dble(cit_tbl_z_dim)
! 
!   ! Here we again fudge cutoff fractions by 0.0001d0 to thwart rounding errs.
! 
!   cut_x_frac = cut_factor(1) * (cutlist / ucell_stk(1,1)) + 0.0001d0
!   cut_y_frac = cut_factor(2) * (cutlist / ucell_stk(2,2)) + 0.0001d0
!   cut_z_frac = cut_factor(3) * (cutlist / ucell_stk(3,3)) + 0.0001d0
  
    ! Set up the bucket and translation index mapping arrays.  We really only need
    ! to do this each time if the cit table dimensions are changing, but there
***************
*** 536,542 ****
      iaci = ntypes_stk * (img(img_i)%iac - 1)
  #endif
          
!     ! Determine the bucket ranges that need to be searched:
  
      x_i = img(img_i)%imgcrd(1)
      y_i = img(img_i)%imgcrd(2)
--- 543,551 ----
      iaci = ntypes_stk * (img(img_i)%iac - 1)
  #endif
          
!     atm_i = img_atm_map(img_i)
! 
!     ! These are imaged coordinates:
  
      x_i = img(img_i)%imgcrd(1)
      y_i = img(img_i)%imgcrd(2)
***************
*** 548,559 ****
        i_tranvec(3, i) = tranvec(3, i) - z_i 
      end do
  
-     atm_i = img_atm_map(img_i)
- 
-     x_fract_i = fraction(1, atm_i)
-     y_fract_i = fraction(2, atm_i)
-     z_fract_i = fraction(3, atm_i)
- 
      ! Convert the excluded atom mask info into excluded image mask info.  At
      ! this point, the image mask info will include all excluded images.  After
      ! final processing, it will only include those processed by this process
--- 557,562 ----
***************
*** 570,588 ****
        excl_img_map(img_j) = 1                   ! .ne. 0 = excluded
      end do
  
!     x_bkts_lo = int((x_fract_i - x_fract_cut + 1.d0) * dble_x_tbl_dim)
!     x_bkts_hi = int((x_fract_i + x_fract_cut + 1.d0) * dble_x_tbl_dim)
  
!     y_bkts_lo = int((y_fract_i - y_fract_cut + 1.d0) * dble_y_tbl_dim)
!     y_bkts_hi = int((y_fract_i + y_fract_cut + 1.d0) * dble_y_tbl_dim)
  
!     z_bkts_lo = int((z_fract_i - z_fract_cut + 1.d0) * dble_z_tbl_dim)
!     z_bkts_hi = int((z_fract_i + 1.d0) * dble_z_tbl_dim)
  
-     i_bkt = x_bkts(int((x_fract_i + 1.d0) * dble_x_tbl_dim)) + &
-             y_bkts(int((y_fract_i + 1.d0) * dble_y_tbl_dim)) + &
-             z_bkts(z_bkts_hi)
-                                         
      z_loop: &
      do z_bkts_idx = z_bkts_lo, z_bkts_hi
  
--- 573,599 ----
        excl_img_map(img_j) = 1                   ! .ne. 0 = excluded
      end do
  
!     ! These are indexes into the 3D cit:
! 
!     x_i_idx = int(fraction(1, atm_i) * scale_fac_x)
!     y_i_idx = int(fraction(2, atm_i) * scale_fac_y)
!     z_i_idx = int(fraction(3, atm_i) * scale_fac_z)
! 
!     ! This is the flat cit bucket index:
! 
!     i_bkt = x_i_idx + cit_tbl_x_dim * (y_i_idx + z_i_idx * cit_tbl_y_dim)
! 
!     ! Determine the bucket ranges that need to be searched:
  
!     x_bkts_lo = int((1.d0 + fraction(1, atm_i) - cut_x_frac) * scale_fac_x)
!     x_bkts_hi = int((1.d0 + fraction(1, atm_i) + cut_x_frac) * scale_fac_x)
!     y_bkts_lo = int((1.d0 + fraction(2, atm_i) - cut_y_frac) * scale_fac_y)
!     y_bkts_hi = int((1.d0 + fraction(2, atm_i) + cut_y_frac) * scale_fac_y)
!     z_bkts_lo = int((1.d0 + fraction(3, atm_i) - cut_z_frac) * scale_fac_z)
!     z_bkts_hi = z_i_idx + cit_tbl_z_dim
  
!     ! We have not yet implemented img_hi fiddling per the orthog version...
  
      z_loop: &
      do z_bkts_idx = z_bkts_lo, z_bkts_hi
  

*** ew_direct_cit.f90	Mon Feb 16 17:31:07 2004
--- ew_direct_cit.f90	Fri Jul  8 15:52:12 2005
***************
*** 130,140 ****
                         cit_tranvec, atm_nb_maskdata, atm_nb_mask, &
                         cit_atm_img_map, &
                         dont_skip_belly_pairs, &
- #ifdef MPI
                         fraction, ifail)
- #else
-                        ifail)
- #endif
      else
        call get_nb_list_nonorthog(atm_cnt, crd_idx_tbl, cit_img, cit_img_atm_map, &
  #ifdef DIRFRC_BIGCACHE_OPT
--- 130,136 ----
***************
*** 790,795 ****
--- 786,806 ----
  
    end if
  
+   ! We must have fractional coordinates in the range of 0.0 - 0.999...
+   ! The algorithm used above will produce fractionals in the range of 0.0 -
+   ! 1.0, with some possibility of slight underflow (neg value) due to
+   ! rounding error (confirmed). SO we force fractionals to be nonredundant
+   ! and within the anticipated range here.
+ 
+   do i = 1, atm_cnt
+     if (fraction(1, i) .lt. 0.d0) fraction(1, i) = fraction(1, i) + 1.d0
+     if (fraction(1, i) .ge. 1.d0) fraction(1, i) = fraction(1, i) - 1.d0
+     if (fraction(2, i) .lt. 0.d0) fraction(2, i) = fraction(2, i) + 1.d0
+     if (fraction(2, i) .ge. 1.d0) fraction(2, i) = fraction(2, i) - 1.d0
+     if (fraction(3, i) .lt. 0.d0) fraction(3, i) = fraction(3, i) + 1.d0
+     if (fraction(3, i) .ge. 1.d0) fraction(3, i) = fraction(3, i) - 1.d0
+   end do
+ 
    return
  
  end subroutine get_fract_crds

------------------------------------------------------------------------------
 
Temporary workarounds: none