       for (unsigned int j=0; j < n_cols; ++j) {
         smem[ref_id*n_cols+j] = coords[gid*n_cols+j];
       }
-      // load current best mindists into private mem
+      // load current best mindists into registers
       float nh_mindist = nh_dist_ndx[gid];
       float nhhd_mindist = nh_dist_ndx[gid];
       float nh_minndx = nh_dist_ndx[n_rows+gid];
       float nhhd_minndx = nhhd_dist_ndx[n_rows+gid];
-
-      //TODO run over comp frames and compare dist -> result
+      // compare squared distances of reference
+      // to (other) frames in shared mem
       for (unsigned int i=0; i < bsize/2; ++i) {
         float dist2=0.0f;
         for (unsigned int j=0; j < n_cols; ++j) {
           float c = smem[ref_id*n_cols+j] - smem[(i+bsize/2)*n_cols+j];
           dist2 = fma(c, c, dist2);
         }
-
-        //TODO
-        //float mindist = min(dist2, nh_mindist);
-        //if mindist == 0 take max else min
-
-
-      }
+        // frame with min distance (i.e. nearest neighbor)
+        if ((nh_mindist == 0)
+         || (dist2 != 0 && dist2 < nh_mindist)) {
+          nh_mindist = dist2;
+          nh_minndx = i+offset;
+        }
+        // frame with min distance and lower energy
+        if ((nhhd_mindist == 0)
+         || (dist2 != 0 && dist2 < nhhd_mindist)) {
+          nhhd_mindist = dist2;
+          nhhd_minndx = i+offset;
+        }
+      }
+      // write results (dist & ndx) to global buffers
+      nh_dist_ndx[gid] = nh_mindist;
+      nh_dist_ndx[n_rows+gid] = nh_minndx;
+      nhhd_dist_ndx[gid] = nhhd_mindist;
+      nhhd_dist_ndx[n_rows+gid] = nhhd_minndx;
     }
   }
 

      private(i)\
      firstprivate(n_gpus,n_rows,n_cols,gpu_range)\
      shared(partial_pops,radii,coords,sorted_coords,blimits)\
      num_threads(n_gpus)
    for (i=0; i < n_gpus; ++i) {
      // compute partial populations in parallel
      // on all available GPUs

    cudaMalloc((void**) &d_fe
             , sizeof(float) * n_rows);
    for (unsigned int i=0; i < N_STREAMS_NH; ++i) {
      cudaMalloc((void**) &d_nh[i]
               , sizeof(float) * n_rows * 2);
      cudaMalloc((void**) &d_nhhd[i]
               , sizeof(float) * n_rows * 2);
      cudaMemset(&d_nh[i]
               , 0
               , sizeof(float) * n_rows * 2);
      cudaMemset(&d_nhhd[i]
               , 0
               , sizeof(float) * n_rows * 2);
    }
