diff --git a/src/mesh_abaqus.f90 b/src/mesh_abaqus.f90
index 0404ee4ae..336fdb095 100644
--- a/src/mesh_abaqus.f90
+++ b/src/mesh_abaqus.f90
@@ -8,6 +8,7 @@
 module mesh
  use prec
  use mesh_base
+ use geometry_plastic_nonlocal
 
  implicit none
  private
@@ -1909,6 +1910,8 @@ subroutine mesh_build_ipNeighborhood
    enddo
  enddo
  
+  call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood)
+ 
  contains
  !--------------------------------------------------------------------------------------------------
 !> @brief find face-matching element of same type
diff --git a/src/mesh_marc.f90 b/src/mesh_marc.f90
index 6f3d713c3..1b2661b1f 100644
--- a/src/mesh_marc.f90
+++ b/src/mesh_marc.f90
@@ -16,6 +16,7 @@ module mesh
  use numerics
  use FEsolving
  use element
+ use geometry_plastic_nonlocal
 #if defined(DAMASK_HDF5)
  use HDF5_utilities
  use results
@@ -1209,7 +1210,7 @@ subroutine IP_neighborhood2
    endif
    f = f +1 
  enddo
-   
+ call geometry_plastic_nonlocal_set_IPneighborhood(mesh_ipNeighborhood2)
 
   do e = 1,theMesh%nElems
     do i = 1,theMesh%elem%nIPs
diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90
index 7231c5a2b..0336215ea 100644
--- a/src/plastic_nonlocal.f90
+++ b/src/plastic_nonlocal.f90
@@ -932,8 +932,8 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
     nRealNeighbors = 0.0_pReal
     neighbor_rhoTotal = 0.0_pReal
     do n = 1,theMesh%elem%nIPneighbors
-      neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
-      neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
+      neighbor_el = IPneighborhood(1,n,ip,el)
+      neighbor_ip = IPneighborhood(2,n,ip,el)
       no = phasememberAt(1,neighbor_ip,neighbor_el)
       if (neighbor_el > 0 .and. neighbor_ip > 0) then
         neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
@@ -1626,16 +1626,16 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
     
     neighbors: do n = 1,theMesh%elem%nIPneighbors
   
-      neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
-      neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
-      neighbor_n  = mesh_ipNeighborhood(3,n,ip,el)
+      neighbor_el = IPneighborhood(1,n,ip,el)
+      neighbor_ip = IPneighborhood(2,n,ip,el)
+      neighbor_n  = IPneighborhood(3,n,ip,el)
       np = phaseAt(1,neighbor_ip,neighbor_el)
       no = phasememberAt(1,neighbor_ip,neighbor_el)
   
       opposite_neighbor = n + mod(n,2) - mod(n+1,2)
-      opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el)
-      opposite_ip = mesh_ipNeighborhood(2,opposite_neighbor,ip,el)
-      opposite_n  = mesh_ipNeighborhood(3,opposite_neighbor,ip,el)
+      opposite_el = IPneighborhood(1,opposite_neighbor,ip,el)
+      opposite_ip = IPneighborhood(2,opposite_neighbor,ip,el)
+      opposite_n  = IPneighborhood(3,opposite_neighbor,ip,el)
     
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
         neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
@@ -1931,8 +1931,8 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
   !*** Loop thrugh neighbors and check whether there is any compatibility.
   
   neighbors: do n = 1,Nneighbors
-    neighbor_e = mesh_ipNeighborhood(1,n,i,e)
-    neighbor_i = mesh_ipNeighborhood(2,n,i,e)
+    neighbor_e = IPneighborhood(1,n,i,e)
+    neighbor_i = IPneighborhood(2,n,i,e)
     
     
     !* FREE SURFACE