PME: Fix instability in get_center on i386

cp2k · Dec 29, 2021 · c6a8d5e · c6a8d5e
1 parent ebb2aed
commit c6a8d5e
Show file tree

Hide file tree

Showing 6 changed files with 91 additions and 127 deletions.
diff --git a/src/common/structure_factor_types.F b/src/common/structure_factor_types.F
@@ -25,6 +25,7 @@ MODULE structure_factor_types
       COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: shell_ex, shell_ey, shell_ez
       COMPLEX(KIND=dp), DIMENSION(:, :), POINTER :: core_ex, core_ey, core_ez
       INTEGER, DIMENSION(:, :), POINTER :: centre, core_centre, shell_centre
+      REAL(KIND=dp), DIMENSION(:, :), POINTER :: delta, core_delta, shell_delta
       INTEGER :: lb(3)
    END TYPE structure_factor_type
 

diff --git a/src/common/structure_factors.F b/src/common/structure_factors.F
@@ -38,6 +38,7 @@ SUBROUTINE structure_factor_init(exp_igr)
       NULLIFY (exp_igr%shell_ex, exp_igr%shell_ey, exp_igr%shell_ez)
       NULLIFY (exp_igr%core_ex, exp_igr%core_ey, exp_igr%core_ez)
       NULLIFY (exp_igr%centre, exp_igr%shell_centre, exp_igr%core_centre)
+      NULLIFY (exp_igr%delta, exp_igr%shell_delta, exp_igr%core_delta)
 
    END SUBROUTINE structure_factor_init
 
@@ -63,19 +64,17 @@ SUBROUTINE structure_factor_deallocate(exp_igr)
          DEALLOCATE (exp_igr%core_ez)
       END IF
       IF (ASSOCIATED(exp_igr%centre)) THEN
-         DEALLOCATE (exp_igr%centre)
+         DEALLOCATE (exp_igr%centre, exp_igr%delta)
       END IF
       IF (ASSOCIATED(exp_igr%shell_centre)) THEN
-         DEALLOCATE (exp_igr%shell_centre)
+         DEALLOCATE (exp_igr%shell_centre, exp_igr%shell_delta)
       END IF
       IF (ASSOCIATED(exp_igr%core_centre)) THEN
-         DEALLOCATE (exp_igr%core_centre)
+         DEALLOCATE (exp_igr%core_centre, exp_igr%core_delta)
       END IF
 
    END SUBROUTINE structure_factor_deallocate
 
-! **************************************************************************************************
-
 ! **************************************************************************************************
 !> \brief ...
 !> \param bds ...
@@ -100,15 +99,15 @@ SUBROUTINE structure_factor_allocate(bds, nparts, exp_igr, &
       ALLOCATE (exp_igr%ex(bds(1, 1):bds(2, 1) + 1, nparts))
       ALLOCATE (exp_igr%ey(bds(1, 2):bds(2, 2) + 1, nparts))
       ALLOCATE (exp_igr%ez(bds(1, 3):bds(2, 3) + 1, nparts))
-      NULLIFY (exp_igr%centre)
+      NULLIFY (exp_igr%centre, exp_igr%delta)
 
       exp_igr%lb(1) = LBOUND(exp_igr%ex, 1)
       exp_igr%lb(2) = LBOUND(exp_igr%ey, 1)
       exp_igr%lb(3) = LBOUND(exp_igr%ez, 1)
 
       IF (PRESENT(allocate_centre)) THEN
          IF (allocate_centre) THEN
-            ALLOCATE (exp_igr%centre(3, nparts))
+            ALLOCATE (exp_igr%centre(3, nparts), exp_igr%delta(3, nparts))
          END IF
       END IF
 
@@ -117,24 +116,25 @@ SUBROUTINE structure_factor_allocate(bds, nparts, exp_igr, &
             ALLOCATE (exp_igr%shell_ex(bds(1, 1):bds(2, 1) + 1, nshell))
             ALLOCATE (exp_igr%shell_ey(bds(1, 2):bds(2, 2) + 1, nshell))
             ALLOCATE (exp_igr%shell_ez(bds(1, 3):bds(2, 3) + 1, nshell))
-            NULLIFY (exp_igr%shell_centre)
+            NULLIFY (exp_igr%shell_centre, exp_igr%shell_delta)
 
             ALLOCATE (exp_igr%core_ex(bds(1, 1):bds(2, 1) + 1, nshell))
             ALLOCATE (exp_igr%core_ey(bds(1, 2):bds(2, 2) + 1, nshell))
             ALLOCATE (exp_igr%core_ez(bds(1, 3):bds(2, 3) + 1, nshell))
-            NULLIFY (exp_igr%core_centre)
+            NULLIFY (exp_igr%core_centre, exp_igr%core_delta)
 
             IF (PRESENT(allocate_shell_centre)) THEN
                IF (allocate_shell_centre) THEN
-                  ALLOCATE (exp_igr%shell_centre(3, nshell))
-                  ALLOCATE (exp_igr%core_centre(3, nshell))
+                  ALLOCATE (exp_igr%shell_centre(3, nshell), exp_igr%shell_delta(3, nshell))
+                  ALLOCATE (exp_igr%core_centre(3, nshell), exp_igr%core_delta(3, nshell))
                END IF
             END IF
          END IF
       ELSE
          NULLIFY (exp_igr%shell_ex, exp_igr%shell_ey, exp_igr%shell_ez)
          NULLIFY (exp_igr%core_ex, exp_igr%core_ey, exp_igr%core_ez)
          NULLIFY (exp_igr%shell_centre, exp_igr%core_centre)
+         NULLIFY (exp_igr%shell_delta, exp_igr%core_delta)
       END IF
 
    END SUBROUTINE structure_factor_allocate

diff --git a/src/ewald_methods_tb.F b/src/ewald_methods_tb.F
@@ -103,6 +103,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: center
       INTEGER, DIMENSION(3)                              :: nd, npts
       REAL(KIND=dp)                                      :: alpha, dvols, fat(3), ffa, ffb, fint, vgc
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: delta
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: rhos
       REAL(KIND=dp), DIMENSION(3, 3)                     :: f_stress, h_stress
       TYPE(cp_para_env_type), POINTER                    :: para_env
@@ -138,8 +139,8 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
       CALL rs_grid_set_box(grid_spme, rs=rden)
       CALL rs_grid_zero(rden)
 
-      ALLOCATE (center(3, npart))
-      CALL get_center(particle_set, box, center, npts, n)
+      ALLOCATE (center(3, npart), delta(3, npart))
+      CALL get_center(particle_set, box, center, delta, npts, n)
 
       !-------------- DENSITY CALCULATION ----------------
       ipart = 0
@@ -148,7 +149,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
          IF (p1 == 0) EXIT
 
          ! calculate function on small boxes
-         CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., &
+         CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
                         is_shell=.FALSE., unit_charge=.TRUE.)
          rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
 
@@ -200,7 +201,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
                CALL set_list(particle_set, npart, center, p1, rden, ipart)
                IF (p1 == 0) EXIT
                ! calculate function on small boxes
-               CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., &
+               CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
                               is_shell=.FALSE., unit_charge=.TRUE.)
                CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
                atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*mcharge(p1)*fint*dvols
@@ -232,7 +233,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
                      CALL set_list(particle_set, npart, center, p1, rden, ipart)
                      IF (p1 == 0) EXIT
                      ! calculate function on small boxes
-                     CALL get_patch(particle_set, box, green, npts, p1, rhos, &
+                     CALL get_patch(particle_set, delta, green, p1, rhos, &
                                     is_core=.FALSE., is_shell=.FALSE., unit_charge=.TRUE.)
                      ! integrate box and potential
                      CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
@@ -300,7 +301,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
          IF (p1 == 0) EXIT
 
          ! calculate function on small boxes
-         CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., &
+         CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
                         is_shell=.FALSE., unit_charge=.TRUE.)
 
          CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
@@ -330,7 +331,7 @@ SUBROUTINE tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, &
          END IF
       END IF
       DEALLOCATE (rhos)
-      DEALLOCATE (center)
+      DEALLOCATE (center, delta)
 
       CALL timestop(handle)
 
@@ -433,6 +434,7 @@ SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcha
       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: center
       INTEGER, DIMENSION(3)                              :: npts
       REAL(KIND=dp)                                      :: alpha, dvols, fat(3), fint, vgc
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: delta
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: rhos
       TYPE(cp_para_env_type), POINTER                    :: para_env
       TYPE(greens_fn_type), POINTER                      :: green
@@ -467,8 +469,8 @@ SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcha
       CALL rs_grid_set_box(grid_spme, rs=rden)
       CALL rs_grid_zero(rden)
 
-      ALLOCATE (center(3, npart))
-      CALL get_center(particle_set, box, center, npts, n)
+      ALLOCATE (center(3, npart), delta(3, npart))
+      CALL get_center(particle_set, box, center, delta, npts, n)
 
       !-------------- DENSITY CALCULATION ----------------
       ipart = 0
@@ -477,7 +479,7 @@ SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcha
          IF (p1 == 0) EXIT
 
          ! calculate function on small boxes
-         CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., &
+         CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
                         is_shell=.FALSE., unit_charge=.TRUE.)
          rhos(:, :, :) = rhos(:, :, :)*mcharge(p1)
 
@@ -544,7 +546,7 @@ SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcha
          IF (p1 == 0) EXIT
 
          ! calculate function on small boxes
-         CALL get_patch(particle_set, box, green, npts, p1, rhos, is_core=.FALSE., &
+         CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.FALSE., &
                         is_shell=.FALSE., unit_charge=.TRUE.)
 
          CALL dg_sum_patch_force_1d(rpot, rhos, center(:, p1), fint)
@@ -570,7 +572,7 @@ SUBROUTINE tb_spme_zforce(ewald_env, ewald_pw, particle_set, box, gmcharge, mcha
          DEALLOCATE (drpot)
       END IF
       DEALLOCATE (rhos)
-      DEALLOCATE (center)
+      DEALLOCATE (center, delta)
 
       CALL timestop(handle)
 

diff --git a/src/pme.F b/src/pme.F
@@ -171,11 +171,11 @@ SUBROUTINE pme_evaluate(ewald_env, ewald_pw, box, particle_set, vg_coulomb, &
       CALL rs_grid_zero(rden)
 
       IF (rden%desc%parallel .AND. rden%desc%distributed) THEN
-         CALL get_center(particle_set, box, exp_igr%centre, grid_b%npts, 1)
+         CALL get_center(particle_set, box, exp_igr%centre, exp_igr%delta, grid_b%npts, 1)
       END IF
       IF (PRESENT(shell_particle_set) .AND. rden%desc%parallel .AND. rden%desc%distributed) THEN
-         CALL get_center(shell_particle_set, box, exp_igr%shell_centre, grid_b%npts, 1)
-         CALL get_center(core_particle_set, box, exp_igr%core_centre, grid_b%npts, 1)
+         CALL get_center(shell_particle_set, box, exp_igr%shell_centre, exp_igr%shell_delta, grid_b%npts, 1)
+         CALL get_center(core_particle_set, box, exp_igr%core_centre, exp_igr%core_delta, grid_b%npts, 1)
       END IF
 
       !-------------- DENSITY CALCULATION ----------------

diff --git a/src/pme_tools.F b/src/pme_tools.F
@@ -123,37 +123,46 @@ END FUNCTION in_slice
 !> \param part ...
 !> \param box ...
 !> \param center ...
+!> \param delta ...
 !> \param npts ...
 !> \param n ...
 ! **************************************************************************************************
-   SUBROUTINE get_center(part, box, center, npts, n)
+   SUBROUTINE get_center(part, box, center, delta, npts, n)
 
       TYPE(particle_type), DIMENSION(:), INTENT(IN)      :: part
       TYPE(cell_type), POINTER                           :: box
       INTEGER, DIMENSION(:, :), INTENT(OUT)              :: center
+      REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: delta
       INTEGER, DIMENSION(:), INTENT(IN)                  :: npts
       INTEGER, INTENT(IN)                                :: n
 
       INTEGER                                            :: ipart, mp
       REAL(KIND=dp)                                      :: rmp
-      REAL(KIND=dp), DIMENSION(3)                        :: gp, s
+      REAL(KIND=dp), DIMENSION(3)                        :: ca, gp, s
 
+      ! The pbc algorithm is sensitive to numeric noise and compiler optimization because of ANINT.
+      ! Therefore center and delta have to be computed simultaneously to ensure they are consistent.
       mp = MAXVAL(npts(:))
       rmp = REAL(mp, KIND=dp)
       DO ipart = 1, SIZE(part)
          ! compute the scaled coordinate of atom ipart
          CALL real_to_scaled(s, part(ipart)%r, box)
-         s = s - NINT(s)
+         s = s - ANINT(s)
+         ! find the continuous ``grid'' point
          gp = REAL(npts, KIND=dp)*s
          ! find the closest grid point (on big grid)
          IF (MOD(n, 2) == 0) THEN
             center(:, ipart) = INT(gp + rmp) - mp
+            ca(:) = REAL(center(:, ipart), KIND=dp) + 0.5_dp
          ELSE
             center(:, ipart) = NINT(gp)
+            ca(:) = REAL(center(:, ipart), KIND=dp)
          END IF
          center(:, ipart) = center(:, ipart) + npts(:)/2
          center(:, ipart) = MODULO(center(:, ipart), npts(:))
          center(:, ipart) = center(:, ipart) - npts(:)/2
+         ! find the distance vector
+         delta(:, ipart) = gp - ca(:)
       END DO
 
    END SUBROUTINE get_center