diff --git a/R/AAA.R b/R/AAA.R
index 3b1309c..b12fd32 100644
--- a/R/AAA.R
+++ b/R/AAA.R
@@ -1,3 +1,10 @@
+.onLoad <- function(lib, pkg) {
+  n <- suppressWarnings(as.integer(Sys.getenv("OMP_NUM_THREADS", unset = NA_character_)))
+  if (!is.na(n) && n >= 1L) {
+    resemble_set_omp_threads(n)
+  }
+}
+
 .onAttach <- function(lib, pkg) {
   pkg_v <- pkg_info()
   
diff --git a/R/RcppExports.R b/R/RcppExports.R
index b9be50a..4c2e310 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1150,3 +1150,9 @@ ith_pred_cpp <- function(plslib, xscale, Xu, dxrxu = NULL) {
     .Call(`_resemble_ith_pred_cpp`, plslib, xscale, Xu, dxrxu)
 }
 
+#' @keywords internal
+#' @noRd
+resemble_set_omp_threads <- function(n) {
+    invisible(.Call(`_resemble_resemble_set_omp_threads`, n))
+}
+
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index d096d2b..53f5b41 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -566,6 +566,16 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// resemble_set_omp_threads
+void resemble_set_omp_threads(int n);
+RcppExport SEXP _resemble_resemble_set_omp_threads(SEXP nSEXP) {
+BEGIN_RCPP
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type n(nSEXP);
+    resemble_set_omp_threads(n);
+    return R_NilValue;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_resemble_fast_diss", (DL_FUNC) &_resemble_fast_diss, 3},
@@ -605,6 +615,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_resemble_ith_local_fit", (DL_FUNC) &_resemble_ith_local_fit, 10},
     {"_resemble_final_fits_cpp", (DL_FUNC) &_resemble_final_fits_cpp, 9},
     {"_resemble_ith_pred_cpp", (DL_FUNC) &_resemble_ith_pred_cpp, 4},
+    {"_resemble_resemble_set_omp_threads", (DL_FUNC) &_resemble_resemble_set_omp_threads, 1},
     {NULL, NULL, 0}
 };
 
diff --git a/src/diss_helpers.cpp b/src/diss_helpers.cpp
index 7df1e7a..aa09985 100644
--- a/src/diss_helpers.cpp
+++ b/src/diss_helpers.cpp
@@ -1253,6 +1253,21 @@ List top_k_neighbors(
   for (int j = 0; j < m; j++) {
     out[j] = wrap(results[j]);
   }
-  
+
   return out;
 }
+
+//' Set the number of OpenMP threads
+//'
+//' Calls \code{omp_set_num_threads()} so that the limit takes effect for all
+//' subsequent parallel regions, regardless of when \code{OMP_NUM_THREADS} was
+//' set in the environment.
+//' @param n Number of threads.
+//' @keywords internal
+//' @noRd
+// [[Rcpp::export(visible = false)]]
+void resemble_set_omp_threads(int n) {
+#ifdef _OPENMP
+  omp_set_num_threads(n);
+#endif
+}
diff --git a/vignettes/ag-evolutionary-subset-search.qmd b/vignettes/ag-evolutionary-subset-search.qmd
index 14314ce..aa14827 100644
--- a/vignettes/ag-evolutionary-subset-search.qmd
+++ b/vignettes/ag-evolutionary-subset-search.qmd
@@ -34,10 +34,10 @@ vignette: >
 ```{r}
 #| message: false
 #| echo: false
-library(resemble)
 Sys.setenv(OMP_NUM_THREADS = 2)
 Sys.setenv(OMP_THREAD_LIMIT = 2)
 options(mc.cores = 2)
+library(resemble)
 ```
 
 
@@ -332,18 +332,18 @@ my_control <- gesearch_control(
 
 ```{r}
 #| eval: true
-#| echo: false 
+#| echo: false
 gs <- gesearch(
   Xr = train_x[!is.na(train_y), ], # the spectra of the reference "set/library"
   Yr = train_y[!is.na(train_y)],   # the response of the reference "set/library"
   Xu = test_x_ref, # the available target domain spectra
   Yu = test_y_ref, # the available target domain response (accepts NAs)
-  k = 50, b = 100, retain = 0.95,
+  k = 50, b = 20, retain = 0.95,
   target_size = 120,
   fit_method = fit_pls(ncomp = 20, method = "simpls", scale = TRUE),
   optimization = c("reconstruction", "similarity", "response"),
   control = my_control,
-  seed = 1124, 
+  seed = 1124,
   verbose = FALSE
 )
 ```
@@ -466,18 +466,18 @@ The search process can be conducted when no information on the response variable
 
 ```{r}
 #| eval: true
-#| echo: false 
+#| echo: false
 gs_nr <- gesearch(
   Xr = train_x[!is.na(train_y), ], # the spectra of the reference "set/library"
   Yr = train_y[!is.na(train_y)],   # the response of the reference "set/library"
   Xu = test_x_ref, # the available target domain spectra
   Yu = NULL, # NO available target domain response (accepts NAs)
-  k = 50, b = 80, retain = 0.98,
+  k = 50, b = 15, retain = 0.98,
   target_size = 160,
   fit_method = fit_pls(ncomp = 20, method = "simpls", scale = TRUE),
   optimization = c("reconstruction", "similarity"),
   control = my_control,
-  seed = 1124, 
+  seed = 1124,
   verbose = FALSE
 )
 ```
@@ -515,19 +515,19 @@ In some cases, the response values for the target set may not be available, but
 
 ```{r}
 #| eval: true
-#| echo: false 
+#| echo: false
 gs_nr2 <- gesearch(
   Xr = train_x[!is.na(train_y), ], # the spectra of the reference "set/library"
   Yr = train_y[!is.na(train_y)],   # the response of the reference "set/library"
   Xu = test_x_ref, # the available target domain spectra
   Yu = NULL, # NO available target domain response (accepts NAs)
   Yu_lims = c(0, 5), # the response range in the target domain
-  k = 50, b = 80, retain = 0.98,
+  k = 50, b = 15, retain = 0.98,
   target_size = 160,
   fit_method = fit_pls(ncomp = 20, method = "simpls", scale = TRUE),
   optimization = c("reconstruction", "similarity", "range"),
   control = my_control,
-  seed = 1124, 
+  seed = 1124,
   verbose = FALSE
 )
 ```