Avoid division by zero by updating expression

Signed-off-by: Igor S. Gerasimov <[email protected]>
grimme-lab · Feb 10, 2025 · 48b195e · 48b195e
1 parent 6d86ebb
commit 48b195e
Showing 1 changed file with 3 additions and 9 deletions.
diff --git a/src/iff/iff_energy.f90 b/src/iff/iff_energy.f90
@@ -95,9 +95,6 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
       logical :: linear, ATM, ijh, ijx
       character(len=4) :: pgroup
 
-      !> threshold
-      real(wp), parameter :: rClose = 1e-12_wp
-
       !> Initial stuff
       metal = 1
       metal(1:2) = 0
@@ -278,18 +275,16 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
             r2 = (AL1(1, i) - AL2(1, j))**2&
               &+ (AL1(2, i) - AL2(2, j))**2&
               &+ (AL1(3, i) - AL2(3, j))**2
-            if (r2 .lt. rClose) cycle
             r = sqrt(r2)
-            esl = esl + AL1(4, i)*AL2(4, j)/(r + r0tmp(jat, iat)/r)
+            esl = esl + AL1(4, i)*AL2(4, j)*r/(r*r + r0tmp(jat, iat))
          end do
 !        LP atom corrections
          do j = 1, n2
             r2 = (A2(1, j) - AL1(1, i))**2&
               &+ (A2(2, j) - AL1(2, i))**2&
               &+ (A2(3, j) - AL1(3, i))**2
-            if (r2 .lt. rClose) cycle
             r = sqrt(r2)
-            esl = esl + AL1(4, i)*q2(j)/(r + r0tmp(j, iat)/r)
+            esl = esl + AL1(4, i)*q2(j)*r/(r*r + r0tmp(j, iat))
          end do
       end do
 !     LP atom corrections
@@ -299,9 +294,8 @@ subroutine iff_e(env, n, n1, n2, at1, at2, neigh, A1, c02, q01, q02, c6ab,&
             r2 = (A1(1, j) - AL2(1, i))**2&
               &+ (A1(2, j) - AL2(2, i))**2&
               &+ (A1(3, j) - AL2(3, i))**2
-            if (r2 .lt. rClose) cycle
             r = sqrt(r2)
-            esl = esl + AL2(4, i)*q1(j)/(r + r0tmp(iat, j)/r)
+            esl = esl + AL2(4, i)*q1(j)*r/(r*r + r0tmp(iat, j))
          end do
       end do