simplify code

m3g · Jul 12, 2024 · 82c781d · 82c781d
1 parent 4058ca0
commit 82c781d
Show file tree

Hide file tree

Showing 7 changed files with 83 additions and 222 deletions.
diff --git a/src/computef.f90 b/src/computef.f90
@@ -21,7 +21,6 @@ subroutine computef(n,x,f)
   double precision :: v1(3), v2(3), v3(3) 
   double precision :: x(n)
   double precision :: f,fparc,fplus
-  double precision :: xtemp, ytemp, ztemp
   double precision :: xbar, ybar, zbar
   double precision :: beta, gama, teta
 
@@ -120,40 +119,21 @@ subroutine computef(n,x,f)
             ! are behind 
 
             if ( fix ) then
-              if (.not. using_pbc) then
-                call add_box_behind(iboxx-1,iboxy,iboxz)
-                call add_box_behind(iboxx,iboxy-1,iboxz)
-                call add_box_behind(iboxx,iboxy,iboxz-1)
-
-                call add_box_behind(iboxx,iboxy-1,iboxz+1)
-                call add_box_behind(iboxx,iboxy-1,iboxz-1)
-                call add_box_behind(iboxx-1,iboxy+1,iboxz)
-                call add_box_behind(iboxx-1,iboxy,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz)
-                call add_box_behind(iboxx-1,iboxy,iboxz-1)
-
-                call add_box_behind(iboxx-1,iboxy+1,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy+1,iboxz-1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz-1)
-              else
-                call add_box_behind(pbc_idx_box(iboxx-1, nboxes(1)),iboxy,iboxz)
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1, nboxes(2)),iboxz)
-                call add_box_behind(iboxx,iboxy,pbc_idx_box(iboxz-1, nboxes(3)))
-
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),iboxz)
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),iboxy,pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),iboxz)
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),iboxy,pbc_idx_box(iboxz-1,nboxes(3)))
-
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-              end if
-
+                call add_box_behind(idx_box(iboxx-1, nboxes(1)),iboxy,iboxz)
+                call add_box_behind(iboxx,idx_box(iboxy-1, nboxes(2)),iboxz)
+                call add_box_behind(iboxx,iboxy,idx_box(iboxz-1, nboxes(3)))
+
+                call add_box_behind(iboxx,idx_box(iboxy-1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(iboxx,idx_box(iboxy-1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),iboxz)
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),iboxy,idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),iboxz)
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),iboxy,idx_box(iboxz-1,nboxes(3)))
+
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
             end if
 
           end if
@@ -185,64 +165,24 @@ subroutine computef(n,x,f)
     do while( icart > 0 )
 
       if(comptype(ibtype(icart))) then
-
-        if (.not. using_pbc) then
-
-          ! Interactions inside box
-
-          f = f + fparc(icart,latomnext(icart))
-
-          ! Interactions of boxes that share faces
-
-          f = f + fparc(icart,latomfirst(i+1,j,k))
-          f = f + fparc(icart,latomfirst(i,j+1,k))
-          f = f + fparc(icart,latomfirst(i,j,k+1))
-
-          ! Interactions of boxes that share axes
-
-          f = f + fparc(icart,latomfirst(i+1,j+1,k))
-          f = f + fparc(icart,latomfirst(i+1,j,k+1))
-          f = f + fparc(icart,latomfirst(i+1,j-1,k))
-          f = f + fparc(icart,latomfirst(i+1,j,k-1))
-          f = f + fparc(icart,latomfirst(i,j+1,k+1))
-          f = f + fparc(icart,latomfirst(i,j+1,k-1))
-
-          ! Interactions of boxes that share vertices
-
-          f = f + fparc(icart,latomfirst(i+1,j+1,k+1))
-          f = f + fparc(icart,latomfirst(i+1,j+1,k-1))
-          f = f + fparc(icart,latomfirst(i+1,j-1,k+1))
-          f = f + fparc(icart,latomfirst(i+1,j-1,k-1))
-
-        else
-
           ! Interactions inside box
-
           f = f + fparc(icart,latomnext(icart))
-
           ! Interactions of boxes that share faces
-
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,k))
-          f = f + fparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),k))
-          f = f + fparc(icart,latomfirst(i,j,pbc_idx_box(k+1, nboxes(3))))
-
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,k))
+          f = f + fparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),k))
+          f = f + fparc(icart,latomfirst(i,j,idx_box(k+1, nboxes(3))))
           ! Interactions of boxes that share axes
-
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),k))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,pbc_idx_box(k+1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),k))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,pbc_idx_box(k-1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),k))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,idx_box(k+1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),k))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,idx_box(k-1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),idx_box(k-1, nboxes(3))))
           ! Interactions of boxes that share vertices
-
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-          f = f + fparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-        end if
-
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),idx_box(k-1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          f = f + fparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),idx_box(k-1, nboxes(3))))
       end if
 
       icart = latomnext(icart)

diff --git a/src/computeg.f90 b/src/computeg.f90
@@ -27,7 +27,6 @@ subroutine computeg(n,x,g)
                       dv3beta(3), dv3gama(3), dv3teta(3) 
   double precision :: v1(3), v2(3), v3(3)
   double precision :: xbar, ybar, zbar
-  double precision :: xtemp, ytemp, ztemp
   double precision :: beta, gama, teta, cb, sb, cg, sg, ct, st
 
   ! Reset gradients
@@ -117,45 +116,23 @@ subroutine computeg(n,x,g)
             lboxnext(ibox) = lboxfirst
             lboxfirst = ibox
 
-            ! Add boxes with fixed atoms which are vicinal to this box, and
-            ! are behind 
-
+            ! Add boxes with fixed atoms which are vicinal to this box, and are behind 
             if ( fix ) then
-              if (.not. using_pbc) then
-                call add_box_behind(iboxx-1,iboxy,iboxz)
-                call add_box_behind(iboxx,iboxy-1,iboxz)
-                call add_box_behind(iboxx,iboxy,iboxz-1)
-
-                call add_box_behind(iboxx,iboxy-1,iboxz+1)
-                call add_box_behind(iboxx,iboxy-1,iboxz-1)
-                call add_box_behind(iboxx-1,iboxy+1,iboxz)
-                call add_box_behind(iboxx-1,iboxy,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz)
-                call add_box_behind(iboxx-1,iboxy,iboxz-1)
-
-                call add_box_behind(iboxx-1,iboxy+1,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy+1,iboxz-1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz+1)
-                call add_box_behind(iboxx-1,iboxy-1,iboxz-1)
-              else
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),iboxy,iboxz)
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1,nboxes(2)),iboxz)
-                call add_box_behind(iboxx,iboxy,pbc_idx_box(iboxz-1, nboxes(3)))
-
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(iboxx,pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),iboxz)
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),iboxy,pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),iboxz)
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),iboxy,pbc_idx_box(iboxz-1,nboxes(3)))
-
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy+1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz+1,nboxes(3)))
-                call add_box_behind(pbc_idx_box(iboxx-1,nboxes(1)),pbc_idx_box(iboxy-1,nboxes(2)),pbc_idx_box(iboxz-1,nboxes(3)))
-              end if
-
-
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),iboxy,iboxz)
+                call add_box_behind(iboxx,idx_box(iboxy-1,nboxes(2)),iboxz)
+                call add_box_behind(iboxx,iboxy,idx_box(iboxz-1, nboxes(3)))
+
+                call add_box_behind(iboxx,idx_box(iboxy-1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(iboxx,idx_box(iboxy-1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),iboxz)
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),iboxy,idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),iboxz)
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),iboxy,idx_box(iboxz-1,nboxes(3)))
+
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy+1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),idx_box(iboxz+1,nboxes(3)))
+                call add_box_behind(idx_box(iboxx-1,nboxes(1)),idx_box(iboxy-1,nboxes(2)),idx_box(iboxz-1,nboxes(3)))
             end if                                 
 
           end if
@@ -186,63 +163,24 @@ subroutine computeg(n,x,g)
       do while ( icart .ne. 0 )
 
         if(comptype(ibtype(icart))) then
-
-          if (.not. using_pbc) then
-
-            ! Interactions inside box
-
-            call gparc(icart,latomnext(icart))
-
-            ! Interactions of boxes that share faces
-
-            call gparc(icart,latomfirst(i+1,j,k))
-            call gparc(icart,latomfirst(i,j+1,k))
-            call gparc(icart,latomfirst(i,j,k+1))
-
-            ! Interactions of boxes that share axes
-
-            call gparc(icart,latomfirst(i+1,j+1,k))
-            call gparc(icart,latomfirst(i+1,j,k+1))
-            call gparc(icart,latomfirst(i+1,j-1,k))
-            call gparc(icart,latomfirst(i+1,j,k-1))
-            call gparc(icart,latomfirst(i,j+1,k+1))
-            call gparc(icart,latomfirst(i,j+1,k-1))
-
-            ! Interactions of boxes that share vertices
-
-            call gparc(icart,latomfirst(i+1,j+1,k+1))
-            call gparc(icart,latomfirst(i+1,j+1,k-1))
-            call gparc(icart,latomfirst(i+1,j-1,k+1))
-            call gparc(icart,latomfirst(i+1,j-1,k-1))
-
-          else
           ! Interactions inside box
-
-            call gparc(icart,latomnext(icart))
-
+          call gparc(icart,latomnext(icart))
           ! Interactions of boxes that share faces
-
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,k))
-            call gparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),k))
-            call gparc(icart,latomfirst(i,j,pbc_idx_box(k+1, nboxes(3))))
-
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,k))
+          call gparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),k))
+          call gparc(icart,latomfirst(i,j,idx_box(k+1, nboxes(3))))
           ! Interactions of boxes that share axes
-
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),k))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,pbc_idx_box(k+1, nboxes(3))))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),k))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),j,pbc_idx_box(k-1, nboxes(3))))
-            call gparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-            call gparc(icart,latomfirst(i,pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),k))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,idx_box(k+1, nboxes(3))))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),k))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),j,idx_box(k-1, nboxes(3))))
+          call gparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          call gparc(icart,latomfirst(i,idx_box(j+1, nboxes(2)),idx_box(k-1, nboxes(3))))
           ! Interactions of boxes that share vertices
-
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j+1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),pbc_idx_box(k+1, nboxes(3))))
-            call gparc(icart,latomfirst(pbc_idx_box(i+1, nboxes(1)),pbc_idx_box(j-1, nboxes(2)),pbc_idx_box(k-1, nboxes(3))))
-          end if
-
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j+1, nboxes(2)),idx_box(k-1, nboxes(3))))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),idx_box(k+1, nboxes(3))))
+          call gparc(icart,latomfirst(idx_box(i+1, nboxes(1)),idx_box(j-1, nboxes(2)),idx_box(k-1, nboxes(3))))
         end if
 
         icart = latomnext(icart)

diff --git a/src/fparc.f90 b/src/fparc.f90
@@ -18,7 +18,7 @@ double precision function fparc(icart,firstjcart)
   integer :: icart,firstjcart
 
   ! LOCAL SCALARS
-  integer :: jcart, icoord
+  integer :: jcart
   double precision :: datom, tol, short_tol, short_tol_penalty, short_tol_scale
   double precision :: vdiff(3)
 
@@ -50,16 +50,8 @@ double precision function fparc(icart,firstjcart)
     !
     ! Otherwise, compute distance and evaluate function for this pair
     !
-    do icoord = 1, 3
-      vdiff(icoord) = xcart(icart, icoord) - xcart(jcart, icoord)
-    end do
-    if ( using_pbc ) then
-      call pbc_vector(vdiff)
-    end if
-
-    datom = ( vdiff(1) )**2 + &
-            ( vdiff(2) )**2 + &
-            ( vdiff(3) )**2
+    vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_box)
+    datom = ( vdiff(1) )**2 + ( vdiff(2) )**2 + ( vdiff(3) )**2
     tol = (radius(icart)+radius(jcart))**2
     if ( datom < tol ) then
       fparc = fparc + fscale(icart)*fscale(jcart)*(datom-tol)**2

diff --git a/src/gparc.f90 b/src/gparc.f90
@@ -17,7 +17,7 @@ subroutine gparc(icart,firstjcart)
   integer :: icart,firstjcart
 
   ! LOCAL SCALARS
-  integer :: jcart, icoord
+  integer :: jcart
   double precision :: datom, dtemp, xdiff, tol, &
                       short_tol, short_tol_scale
   double precision :: vdiff(3)
@@ -50,15 +50,8 @@ subroutine gparc(icart,firstjcart)
     ! Otherwise, compute distance and evaluate function for this pair
     !                     
     tol = (radius(icart)+radius(jcart))**2
-    do icoord = 1, 3
-      vdiff(icoord) = xcart(icart, icoord) - xcart(jcart, icoord)
-    end do
-    if ( using_pbc ) then
-      call pbc_vector(vdiff)
-    end if
-    datom = (vdiff(1))**2 + & 
-            (vdiff(2))**2 + &
-            (vdiff(3))**2
+    vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_box)
+    datom = (vdiff(1))**2 + (vdiff(2))**2 + (vdiff(3))**2
     if( datom < tol ) then
       dtemp = fscale(icart)*fscale(jcart) * 4.d0 * (datom - tol)
       xdiff = dtemp*vdiff(1)

diff --git a/src/pbc.f90 b/src/pbc.f90
@@ -7,32 +7,28 @@ module pbc
   double precision, public :: pbc_box(3) 
   logical, public :: using_pbc = .false.
 
-  public pbc_vector, pbc_idx_box
+  public delta_vector, idx_box
 
   contains
 
-  subroutine pbc_vector(v)
-
+  elemental double precision function delta_vector(v1,v2,pbc_box)
       implicit none
-
-      integer icoord
-
-      double precision, intent(inout) :: v(3)
-
-      do icoord = 1, 3
-          v(icoord) = v(icoord) - pbc_box(icoord) * nint(v(icoord) / pbc_box(icoord))
-      end do
-
-  end subroutine pbc_vector
-
-  integer function pbc_idx_box(ibox, nbox)
-
+      double precision, intent(in) :: v1, v2, pbc_box
+      delta_vector = v1 - v2
+      if (using_pbc) then
+        delta_vector = delta_vector - pbc_box * nint(delta_vector/pbc_box)
+      end if
+ end function delta_vector
+
+  integer function idx_box(ibox, nbox)
       implicit none
       integer ibox, nbox
-
-      pbc_idx_box = modulo(ibox - 1 + nbox, nbox) + 1
-
-  end function pbc_idx_box
+      if (using_pbc) then
+        idx_box = modulo(ibox - 1 + nbox, nbox) + 1
+      else
+        idx_box = ibox
+      end if
+  end function idx_box
 
 end module pbc