BUG: Floating point overflow lead to wrong neighbor list, close #30

Atomistica · May 27, 2020 · 84ae53d · 84ae53d
1 parent f541da6
commit 84ae53d
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 6 deletions.
diff --git a/src/python/f90/python_neighbors.f90 b/src/python/f90/python_neighbors.f90
@@ -638,7 +638,7 @@ recursive subroutine fill_neighbor_list(this, p, error)
     !$omp do
     i_loop: do i = 1, p%nat
        ! Compute the 3-index of the cell of atom i
-       celli = floor(matmul(this%rec_cell_size, PNC3(p, i) - p%lower_with_border) + 1)
+       celli = floor(matmul(this%rec_cell_size, PNC3(p, i) - p%lower_with_border)) + 1
        shift = 0
 
        ! Map current cell back to box and keep track of cell shift
@@ -685,7 +685,7 @@ recursive subroutine fill_neighbor_list(this, p, error)
 
                 do while (j /= -1)
                    ! Compute the 3-index of the cell of atom j
-                   cellj = floor(matmul(this%rec_cell_size, PNC3(p, j) - p%lower_with_border) + 1)
+                   cellj = floor(matmul(this%rec_cell_size, PNC3(p, j) - p%lower_with_border)) + 1
 
                    ! Map current cell back to box and keep track of cell shift
                    shift2 = shift1
@@ -843,8 +843,8 @@ subroutine neighbors_binning_init(this, p, error)
     ny  = cv * ny / dot_product(ny, ny)
     nz  = cv * nz / dot_product(nz, nz)
 
-    ! Now dx, dy, dz needs to be adjusted such that a sphere of size cutoff
-    ! fits into the box spanned by (dx,dy,dz)*(nx,ny,nz)
+    ! Now dx, dy, dz needs to be adjusted such dx*|nx| > cutoff, dy*|ny| > cutoff
+    ! and dz*|nz| > cutoff.
     dx  = int(this%cutoff/sqrt(dot_product(nx, nx))) + 1
     dy  = int(this%cutoff/sqrt(dot_product(ny, ny))) + 1
     dz  = int(this%cutoff/sqrt(dot_product(nz, nz))) + 1

diff --git a/src/standalone/neighbors.f90 b/src/standalone/neighbors.f90
@@ -676,7 +676,7 @@ recursive subroutine fill_neighbor_list(this, p, error)
     !$omp do
     i_loop: do i = 1, p%nat
        ! Compute the 3-index of the current cell
-       cell = floor( matmul(this%rec_cell_size, PNC3(p, i) - p%lower_with_border) + 1 )
+       cell = floor(matmul(this%rec_cell_size, PNC3(p, i) - p%lower_with_border)) + 1
 
        this%seed(i) = cur
 
@@ -707,7 +707,7 @@ recursive subroutine fill_neighbor_list(this, p, error)
           ! re-evaluate the cell index considering the boundary conditions
           if (lebc .and. c(3) /= 0) then
              ! Compute 3-index of the cell for the (Lees-Edwards) shifted grid:
-             cell2     = floor( matmul(this%rec_cell_size, PNC3(p, i) - c(3)*shear_dx - p%lower_with_border) + 1 )
+             cell2     = floor(matmul(this%rec_cell_size, PNC3(p, i) - c(3)*shear_dx - p%lower_with_border)) + 1
              cur_cell  = cell2 + this%d(1:3, x)
 
              c         = 0