code: def dist(x, y, z, L, periodicX, periodicY, periodicZ):
      " " "Computes distances between all particles and places the result in a matrix such that the ij th matrix entry corresponds to the distance between particle i and j" " "
      N = len(x)
      xtemp = np.tile(x, (N, 1))
      dx = xtemp - xtemp.T
      ytemp = np.tile(y, (N, 1))
      dy = ytemp - ytemp.T
      ztemp = np.tile(z, (N, 1))
      dz = ztemp - ztemp.T
  
      # Particles 'feel' each other across the periodic boundaries
      if periodicX:
          dx[dx > L / 2] = dx[dx > L / 2] - L
          dx[dx < -L / 2] = dx[dx < -L / 2] + L
  
      if periodicY:
          dy[dy > L / 2] = dy[dy > L / 2] - L
          dy[dy < -L / 2] = dy[dy < -L / 2] + L
  
      if periodicZ:
          dz[dz > L / 2] = dz[dz > L / 2] - L
          dz[dz < -L / 2] = dz[dz < -L / 2] + L
  
      # Total Distances
      d = np.sqrt(dx**2 + dy**2 + dz**2)
  
      # Mark zero entries with negative 1 to avoid divergences
      d[d == 0] = -1
  
      return d, dx, dy, dz
  
label_desc: Use image node_id 1 for calling a global function with example usage: dist(x, y, z, L, periodicX, periodicY, periodicZ) and returns: d,  dx,  dy,  dz