Compute once those that does not depend on k

Signed-off-by: Igor S. Gerasimov <i.s.ger@ya.ru>
grimme-lab · Nov 27, 2024 · 955ff88 · 955ff88
1 parent 294a774
commit 955ff88
Showing 1 changed file with 16 additions and 17 deletions.
diff --git a/src/constrain_pot.f90 b/src/constrain_pot.f90
@@ -106,14 +106,14 @@ subroutine qpothess2(fix,n,at,xyz,h)
       ii = (fix%atoms(i)-1)*3
       do j = 1, fix%n !inner loop for sum over all atoms
          if (i.ne.j) then
+            rij = xyz(:,fix%atoms(j))-xyz(:,fix%atoms(i))
+            r  = norm2(rij)
+            r2 = r*r
+            r3 = r2*r
+            ij = lin(i-1,j-1)+1
+            r0 = fix%val(ij)
+            dr = r-r0
             do k = 1, 3  !loop diagonal elements (xyz components)
-               rij = xyz(:,fix%atoms(j))-xyz(:,fix%atoms(i))
-               r  = norm2(rij)
-               r2 = r*r
-               r3 = r2*r
-               ij = lin(i-1,j-1)+1
-               r0 = fix%val(ij)
-               dr = r-r0
                dx = xyz(k,fix%atoms(i))-xyz(k,fix%atoms(j))
                dx2 = dx*dx
                iik = lin(ii+k,ii+k)
@@ -127,23 +127,22 @@ subroutine qpothess2(fix,n,at,xyz,h)
          endif
       enddo !end inner loop for sum over all atoms
       do j = i+1, fix%n !loop over the rest (mixed atoms)
+         rij = xyz(:,fix%atoms(j))-xyz(:,fix%atoms(i))
+         r  = norm2(rij)
+         r2 = r*r
+         r3 = r2*r
+         ij = lin(i-1,j-1)+1
+         r0 = fix%val(ij)
+         dr = r-r0
+         jj = (fix%atoms(j)-1)*3
          do k = 1, 3  !loop diagonal elements (xyz components)
-            rij = xyz(:,fix%atoms(j))-xyz(:,fix%atoms(i))
-            r  = norm2(rij)
-            r2 = r*r
-            r3 = r2*r
-            ij = lin(i-1,j-1)+1
-            r0 = fix%val(ij)
-            dr = r-r0
-            jj = (fix%atoms(j)-1)*3
+            dx = xyz(k,fix%atoms(i))-xyz(k,fix%atoms(j))
             do m = 1, 3
                if (k.eq.m) then !same component case
-                  dx = xyz(k,fix%atoms(i))-xyz(k,fix%atoms(j))
                   dx2 = dx*dx
                   ijk = lin(ii+k,jj+k)
                   h(ijk) = h(ijk) - 2.0_wp*fix%fc*(1.0_wp+dx2/r2-dx2*dr/r3-r0/r)
                else !different component case
-                  dx = xyz(k,fix%atoms(i))-xyz(k,fix%atoms(j))
                   dy = xyz(m,fix%atoms(i))-xyz(m,fix%atoms(j))
                   ikjm = lin(ii+k,jj+m)
                   h(ikjm) = h(ikjm) - 2.0_wp*fix%fc*r0*dx*dy/r3