Skip close atoms for avoiding division by zero as well as avoid compl…

…ex type transformations Signed-off-by: Igor S. Gerasimov <i.s.ger@ya.ru>
grimme-lab · Feb 10, 2025 · e9af142 · e9af142
1 parent 5aaab63
commit e9af142
Showing 1 changed file with 9 additions and 3 deletions.
diff --git a/src/iff/iff_energy.f90 b/src/iff/iff_energy.f90
@@ -95,6 +95,9 @@ 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
@@ -275,15 +278,17 @@ 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
-            r = dble(sqrt(real(r2)))
+            if (r2 .lt. rClose) cycle
+            r = sqrt(r2)
             esl = esl + AL1(4, i)*AL2(4, j)/(r + r0tmp(jat, iat)/r)
          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
-            r = dble(sqrt(real(r2)))
+            if (r2 .lt. rClose) cycle
+            r = sqrt(r2)
             esl = esl + AL1(4, i)*q2(j)/(r + r0tmp(j, iat)/r)
          end do
       end do
@@ -294,7 +299,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
-            r = dble(sqrt(real(r2)))
+            if (r2 .lt. rClose) cycle
+            r = sqrt(r2)
             esl = esl + AL2(4, i)*q1(j)/(r + r0tmp(iat, j)/r)
          end do
       end do