********>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