diff --git a/.RData b/.RData
new file mode 100644
index 0000000..e9d8f16
Binary files /dev/null and b/.RData differ
diff --git a/.Rbuildignore b/.Rbuildignore
index 664987a..7ff65b7 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,4 +1,15 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
-.github
-BIGf90_0.3.0.pdf
\ No newline at end of file
+^\.github$
+^\.git$
+^\.gitignore$
+^\.Rhistory$
+^\.RData$
+^\.DS_Store$
+^.*\.DS_Store$
+^BIGf90.*\.pdf$
+^inst/software
+^inst/exec_mac
+^inst/exec_windows
+^inst/data
+^inst/test_data_salmon
diff --git a/.gitignore b/.gitignore
index 3b06184..ff1696a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,9 +1,21 @@
-.DS_store
-.Rproj.user
-.Rbuildignore
-BIGf90.Rproj
-inst/exec_mac
-inst/exec_windows/
-inst/data
-tests/
-.Rhistory
\ No newline at end of file
+# ---- OS / editor cruft ----
+.DS_Store
+**/.DS_Store
+.Rhistory
+.RData
+.Rproj.user
+
+# ---- RStudio project file ----
+BIGf90.Rproj
+
+# ---- Bundled BLUPf90 executables: keep locally, never push (platform-specific, ~165 MB) ----
+inst/software/
+inst/exec_mac/
+inst/exec_windows/
+inst/data/
+
+# ---- CVA/validation run outputs & local test data: keep locally, never push (~427 MB) ----
+inst/test_data_salmon/
+
+# NOTE: .Rbuildignore and tests/ are intentionally NOT ignored so they stay tracked
+#       (.Rbuildignore controls the built package; tests/ will hold testthat unit tests).
diff --git a/DESCRIPTION b/DESCRIPTION
index 70a540a..1687265 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,22 +1,25 @@
 Package: BIGf90
-Title: R face front for running K-fold crossvalidation, estimating ebvs and variance component estimation with Blupf90 modules
-Version: 0.3.2
+Title: R Face Front for Running K-Fold Crossvalidation, Estimating EBVs and Variance Component Estimation with Blupf90 Modules
+Version: 0.4.0
 Authors@R: c(person("Josue", "Chinchilla-Vargas", role = c("aut", "cre"), email = "jc3635@cornell.edu"),
              person("Cristiane", "Taniguti", role = "aut"),
              person("Alexander", "Sandercock", role = "aut"),
              person("Breeding Insight Team", role = "ctb"))
 Maintainer: Josue Chinchilla-Vargas <jc3635@cornell.edu>
-Description: This package provides R functions to run several BLUPf90 modules. Along with an R function to run K-fold crossvalidation for univariate analyses through blupf90+. 
-    The output table from the K-fold crossvalidation function calculates accuracy as cor(y,ebv^) and cor(y*,ebv^) along with bias of ebvs calculated as reg(y,ebv.)
-    You will need to create the .par file to feed into Renumf90 manually. Once this is done, there are functions to run Renumf90, Blupf90+, Predictf90, Gibbsf90+ and Postgibbsf90.
+Description: Provides R functions to run several BLUPf90 modules. Along with an R function to run K-fold crossvalidation for univariate analyses through blupf90+. 
+    The output table from the K-fold crossvalidation function reports predictive ability as cor(y*, ebv^) and bias as the regression of y* on ebv^, with accuracy as cor(y*, ebv^)/sqrt(h2), where y* is the phenotype corrected for all fixed effects.
+    You can create the .par file to feed into Renumf90 manually or with write_par(). Once this is done, there are functions to run Renumf90, Blupf90+, Predictf90, Gibbsf90+ and Postgibbsf90.
     If you need to learn how to use the blupf90 suite of programs refer to nce.ads.uga.edu/wiki/doku.php?id=start. Please remember to cite Blupf90 appropriately along with this package  when used for publications.
     To run the functions in this package you will need to have a directory with all the blupf90 executables so that you can indicate the path for R to find the execs.
     As of version 0.3.0, all functions have been tested in PC and Unix environments, we have noticed that PC environments tend to have issues about permissions, so it is advised to run RStudio as administrator when planning to use this package.
 License: Apache License (== 2.0)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
-Imports: 
-    base (>= 4.3.1),
+Imports:
     dplyr (>= 1.1.4),
+    stats,
     utils (>= 4.3.1)
+Suggests:
+    testthat (>= 3.0.0)
+Config/testthat/edition: 3
+Config/roxygen2/version: 8.0.0
diff --git a/NAMESPACE b/NAMESPACE
index 757a28f..e79aff4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,4 +7,5 @@ export(run_gibbs)
 export(run_postgibbs)
 export(run_predict)
 export(run_renum)
+export(write_par)
 import(dplyr)
diff --git a/NEWS.md b/NEWS.md
index 7b14e6b..7dc1c30 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,14 @@
+# BIGf90 0.4.0
+
+* New function write_par() to build the raw RENUMF90 parameter (.par) file from R instead of writing it by hand
+* Effects are given as a named list keyed by a label; each effect gives its data-file column (col, or pos as an alias), type (cross/cov) and class (numer/alpha); an optional comment overrides the EFFECT #... label
+* random marks the random effect by its column (or label); the random effect is written last automatically so RANDOM follows it, with the pedigree FILE / FILE_POS block after
+* Supports genomic models (snp_file), permanent-environment / repeatability terms (optional = "pe") and (CO)VARIANCES starting values
+* Multi-trait models: give a trait vector with matrix RESIDUAL_VARIANCE / (CO)VARIANCES, and effect columns repeat per trait
+* Reproduces the sealice test .par files (pedigree and genomic)
+* Added a testthat suite covering write_par() and the create_folds() helper
+
+
 # BIGf90 0.3.1
 
 * Users now can have executable files, input files, and output files in different locations
diff --git a/R/bf90_cv.R b/R/bf90_cv.R
index 623e17b..dd32ac0 100644
--- a/R/bf90_cv.R
+++ b/R/bf90_cv.R
@@ -4,7 +4,7 @@
 #'
 #' This function sets up and runs a K-fold cross-validation analysis (CVA) using blupf90+ and predictf90.
 #' The function run_renumf90 needs to be used beforehand to process a .par file created by the user.
-#' This function calculates 2 accuracy estimates: correlation between raw phenotypes and ebvs divided by the square-root of narrow sense heritabilty and correlation between corrected phenotypes and ebvs along with bias estimations calculated as the regression of the phenotypes on the ebvs.
+#' Using the phenotype corrected for all fixed effects (y*), the function reports predictive ability as the correlation between y* and the ebvs, accuracy as that correlation divided by the square-root of the narrow-sense heritability, and bias as the regression of y* on the ebvs.
 #'
 #' @param missing_value_code code used in the .par file after OPTION MISSING to indicate missing phenotype, if this option is no use, this value must be 0.
 #' @param random_effect_col Column where random effect is located, found under RANDOM_GROUP in the renf90.par file.
@@ -12,8 +12,6 @@
 #' @param num_runs Number of independent cross-validation runs to be performed.
 #' @param num_folds Number of folds to be generated within each independent run.
 #' @param output_table_name Name of the final tab-separated out-up file. This field should be in quotes "".
-#' @param renf90_ped_name Name of pedigree file generated by renumf90. This field should be in quotes "".
-#' @param snp_file_name the name of the genotype file to be used. If use, this field should be in quotes"". Default for the function is no genotype file.
 #' @param path_2_execs path to a folder that holds all blupf90 executables that ill be used (blupf90+,predictf90). This field should be in quotes "".
 #' @param input_files_dir directory containing files renf90.par renf90.fields renf90.inb renf90.tables renf90.dat
 #' @param output_files_dir path to directory to store the output files
@@ -33,8 +31,8 @@
 #' #      num_runs = 5,
 #' #      num_folds = 10,
 #' #      output_table_name = "example_run",
-#' #     renf90_ped_name = "renadd03.ped",
-#' #     snp_file_name = "my_genos.geno" )
+#' #      input_files_dir = ".",
+#' #      output_files_dir = ".")
 #'
 #'
 #' @export
@@ -47,21 +45,19 @@ bf90_cv <- function(missing_value_code = NULL,
                     input_files_dir = ".",
                     output_files_dir = ".",
                     output_table_name = NULL,
-                    renf90_ped_name= NULL,
-                    snp_file_name = NULL,
                     seed = 101919,
                     verbose = TRUE) {
   
+  input_files_dir <- normalizePath(input_files_dir)
   
-  #check id pedigree and genotype file are present
-  if(is.null(renf90_ped_name)) stop("Specify renf90_ped_name")
-  
-  if (!file.exists(renf90_ped_name)) {
-    stop("Parameter file not found at: ", renf90_ped_name)
+  # Checks
+  output_files_dir <- normalizePath(output_files_dir)
+  if(file.exists(output_files_dir)){
+    check_files <- list.files(output_files_dir)
+    if(length(check_files) > 0) warning(paste("Directory", output_files_dir, "is not empty. Some files may be replaced."))
+  } else {
+    stop(paste("Directory", output_files_dir, "does not exist. Create it before running the function."))
   }
-  renf90_ped_name <- normalizePath(renf90_ped_name)
-  
-  input_files_dir <- normalizePath(input_files_dir)
   
   if(is.null(missing_value_code)) stop("Specify missing_value_code")
   if(is.null(random_effect_col)) stop("Specify random_effect_col")
@@ -70,36 +66,42 @@ bf90_cv <- function(missing_value_code = NULL,
   if(is.null(num_folds)) stop("Specify num_folds")
   
   if (!file.exists(file.path(input_files_dir,"renf90.par"))) stop("File 'renf90.par' not found at: ", input_files_dir)
-  if (!file.exists(file.path(input_files_dir,"renf90.inb"))) stop("File 'renf90.inb' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.fields"))) stop("File 'renf90.fields' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.tables"))) stop("File 'renf90.tables' not found at: ", input_files_dir)
   if (!file.exists(file.path(input_files_dir,"renf90.dat"))) stop("File 'renf90.dat' not found at: ", input_files_dir)
   
-  if (!is.null(snp_file_name)) {
-    if (!file.exists(snp_file_name)) {
-      stop("Parameter file not found at:", snp_file_name)
-      snp_file_name <- normalizePath(snp_file_name)
-    }
+  # Check files in the parameter files
+  parfile <- readLines(file.path(input_files_dir,"renf90.par"))
+  
+  snpfile <- grep("SNP_file", parfile)
+  if(length(snpfile) != 0) {
+    snp_file_name <- sapply(strsplit(parfile[snpfile], " "), function(x) x[length(x)])
+    if(!file.exists(file.path(input_files_dir, snp_file_name))) stop(paste("File", file.path(input_files_dir, snp_file_name),
+                                                                           "especified in parameter file",file.path(input_files_dir,"renf90.par"),
+                                                                           "line", snpfile, "does not exist."))
+    snp_file_name <- normalizePath(file.path(input_files_dir, snp_file_name))
+  } else snp_file_name <- NULL
+
+  # cleaned-SNP and saved G-inverse filenames (built once on the full data, reused per fold)
+  if(!is.null(snp_file_name)){
+    snp_base  <- base::basename(snp_file_name)      # e.g. sealice_match.geno
+    clean_snp <- base::paste0(snp_base, "_clean")   # produced by saveCleanSNPs
   }
   
-  if(is.null(output_table_name)) stop("Define output table name.")
+  pedfile <- grep(" FILE", parfile) + 1
+  if(length(pedfile) != 0) {
+    renf90_ped_name <- gsub(" ", "", parfile[pedfile])
+    if(!file.exists(file.path(input_files_dir, renf90_ped_name))) stop(paste("File", file.path(input_files_dir, renf90_ped_name),
+                                                                             "especified in parameter file",file.path(input_files_dir,"renf90.par"),
+                                                                             "line", pedfile, "does not exist."))
+    renf90_ped_name <- normalizePath(file.path(input_files_dir, renf90_ped_name))
+  } else renf90_ped_name <- NULL
   
-  if (!file.exists(renf90_ped_name)) {
-    stop("Parameter file not found at: ", renf90_ped_name)
-  }
+  if(is.null(output_table_name)) stop("Define output table name.")
   
   # set seed for reproducibility
   base::set.seed(seed)
   
-  # Checks
-  output_files_dir <- normalizePath(output_files_dir)
-  if(file.exists(output_files_dir)){
-    check_files <- list.files(output_files_dir)
-    if(length(check_files) > 0) warning(paste("Directory", output_files_dir, "is not empty. Some files may be replaced."))
-  } else {
-    stop(paste("Directory", output_files_dir, "does not exist. Create it before running the function."))
-  }
-  
   path_2_execs <- normalizePath(path_2_execs)
   
   # Inform parameters and directories set
@@ -130,19 +132,52 @@ bf90_cv <- function(missing_value_code = NULL,
   
   # Run BF90 programs for the whole dataset
   setwd(input_files_dir)
-  output <- execute_command(command = paste0(file.path(path_2_execs, blup)," ", "renf90.par"), logfile = "run_blup.log")
-  command_predict <- paste0(file.path(path_2_execs, predict), " ", paste0("renf90.par"))
+  # For genomic analyses, make sure the G-inverse (Gi) and cleaned SNPs are saved on
+  # this full-data run so every fold can reuse them (built once here).
+  blup_parfile <- "renf90.par"
+  if(!is.null(snp_file_name)){
+    blup_par <- base::readLines("renf90.par")
+    if(!any(base::grepl("saveGInverse",  blup_par))) blup_par <- c(blup_par, "OPTION saveGInverse")
+    if(!any(base::grepl("saveCleanSNPs", blup_par))) blup_par <- c(blup_par, "OPTION saveCleanSNPs")
+    base::writeLines(blup_par, "renf90_blup.par")
+    blup_parfile <- "renf90_blup.par"
+  }
+  output <- execute_command(command = paste0(file.path(path_2_execs, blup)," ", blup_parfile), logfile = "run_blup.log")
+
+  # predictf90 on a copy of renf90.par that adjusts the phenotype for all effects
+  # except the random (animal) effect -> yhat = corrected phenotype (y*)
+  predict_par <- c(base::readLines("renf90.par"), base::paste("OPTION include_effects", random_effect_col))
+  base::writeLines(predict_par, "renf90_predict.par")
+  command_predict <- paste0(file.path(path_2_execs, predict), " ", "renf90_predict.par")
   output <- execute_command(command = command_predict, logfile = "run_predict.log")
   
-  files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions")
-  for(i in 1:length(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
+  if(is.null(snp_file_name)){
+    files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions")
+  } else {
+    files_res <- c("run_blup.log", "bvs.dat", "bvs2.dat", "yhat_residual", "solutions",
+                   "freqdata.count", "freqdata.count.after.clean", "Gen_call_rate", "Gen_conflicts",
+                   "run_predict.log", "sum2pq",
+                   "Gi", clean_snp, base::paste0(clean_snp, "_XrefID"))   # reused by every fold
+  }
+
+  for(i in 1:length(files_res)) if(file.exists(files_res[i])) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
   
   # Prepare files for each BLUP run
   renf90 <- base::readLines(paste0("renf90.par"))
   #ped_file <- dirname(renf90_ped_name)
   fields_file <- base::readLines(paste0("renf90.fields"))
-  inb_file <- base::readLines(paste0("renf90.inb"))
   tables_file <- base::readLines(paste0("renf90.tables"))
+
+  # Genomic: point each fold at the G-inverse + cleaned SNPs built on the full data,
+  # so folds skip genotype QC and G construction (identical results, much faster).
+  if(!is.null(snp_file_name)){
+    renf90 <- base::gsub("saveGInverse",  "readGInverse",       renf90, fixed = TRUE)
+    renf90 <- base::gsub("saveCleanSNPs", "no_quality_control", renf90, fixed = TRUE)
+    if(!any(base::grepl("readGInverse",       renf90))) renf90 <- c(renf90, "OPTION readGInverse")
+    if(!any(base::grepl("no_quality_control", renf90))) renf90 <- c(renf90, "OPTION no_quality_control")
+    renf90 <- base::gsub(snp_base, clean_snp, renf90, fixed = TRUE)   # SNP_file -> cleaned set
+    reuse_files <- base::file.path(output_files_dir, c("Gi", clean_snp, base::paste0(clean_snp, "_XrefID")))
+  } else reuse_files <- NULL
   
   # Read and preprocess the phenotype data
   bf90_phenos <- utils::read.table(paste0("renf90.dat"), sep = " ", header = FALSE) %>%
@@ -153,15 +188,20 @@ bf90_cv <- function(missing_value_code = NULL,
                                   dplyr::select((random_effect_col + 1)))
   
   folds <- base::lapply(data_shuffled, function(x) create_folds(x, num_folds))
-  mutated_data <- base::lapply(folds, function(f) mutate_folds(bf90_phenos, f, num_folds,missing_value_code))
+  mutated_data <- base::lapply(folds, function(f) mutate_folds(bf90_phenos, f, num_folds, missing_value_code, random_effect_col + 1))  # id col = random_effect_col+1
   
   setwd(output_files_dir)
   for (run in 1:num_runs) {
     for (fold in 1:num_folds) {
       dir_path <- base::sprintf("run%d/fold%d", run, fold)
-      create_cv_datasets(run, fold, mutated_data[[run]][[fold]],
-                         dir_path, renf90, renf90_ped_name,
-                         input_files_dir, snp_file_name)
+      create_cv_datasets(run, 
+                         fold, 
+                         data_frame =  mutated_data[[run]][[fold]],
+                         dir_path, 
+                         renf90, 
+                         renf90_ped_name,
+                         input_files_dir,
+                         reuse_files)
     }
   }
   
@@ -178,7 +218,7 @@ bf90_cv <- function(missing_value_code = NULL,
       data_file <- base::sprintf("renf90_run%d_fold%d.dat", run, fold)
       masked_ids <- utils::read.table(data_file) %>%
         dplyr::filter(V1 == missing_value_code) %>%
-        dplyr::select(4) %>%
+        dplyr::select(random_effect_col + 1) %>%   # id column
         base::unlist()
       
       ebvs_for_cv <- utils::read.table("solutions", header = FALSE, skip = 1) %>%
@@ -195,7 +235,6 @@ bf90_cv <- function(missing_value_code = NULL,
   
   # Calculate correlations and bias
   corrected_phenos <- utils::read.table(paste0(output_files_dir, "/yhat_residual"), header = FALSE) %>% dplyr::select(1, 2)
-  raw_phenos <- utils::read.table(paste0(input_files_dir, "/renf90.dat")) %>% dplyr::select(4, 1)
   
   ystar_correlations <- base::numeric(num_runs)
   bias_list <- base::numeric(num_runs)
@@ -205,9 +244,7 @@ bf90_cv <- function(missing_value_code = NULL,
     ystar <- dplyr::inner_join(ebvs_for_cv_runs[[base::sprintf("ebvs_for_cv_run%d", i)]], corrected_phenos, by = c("V3" = "V1"))
     ystar_correlations[i] <- base::round(stats::cor(ystar$V4, ystar$V2), 3)
     
-    yraw <- dplyr::inner_join(ebvs_for_cv_runs[[base::sprintf("ebvs_for_cv_run%d", i)]], raw_phenos, by = c("V3" = "V4"))
-    
-    model <- stats::lm(V1 ~ V4, data = yraw)
+    model <- stats::lm(V2 ~ V4, data = ystar)   # bias (LR b1): corrected phenotype (V2) regressed on EBV (V4)
     bias_list[i] <- base::round(stats::coefficients(model)["V4"], 3)
     
     accuracy_list[i] <- base::round(ystar_correlations[i] / base::sqrt(h2), 3)
diff --git a/R/run_blup.R b/R/run_blup.R
index f28cc0a..0447ed1 100644
--- a/R/run_blup.R
+++ b/R/run_blup.R
@@ -33,15 +33,15 @@ run_blup <- function(path_2_execs) {
   }
 
   # Check if executable exists
-  if (!file.exists(paste0(path_2_execs, blup))) {
-    stop("Executable not found at: ", paste0(path_2_execs, blup))
+  if (!file.exists(file.path(path_2_execs, blup))) {
+    stop("Executable not found at: ", file.path(path_2_execs, blup))
   }
   if (!file.exists("renf90.par")) {
     stop("Parameter file not found: renf90.par")
   }
 
   # Run the command
-  output <- execute_command(command = paste0(path_2_execs, blup, " renf90.par"), logfile = "run_blup.log")
+  output <- execute_command(command = paste0(file.path(path_2_execs, blup), " renf90.par"), logfile = "run_blup.log")
 
   # Capture and print the log file content
   if (file.exists("run_blup.log")) {
diff --git a/R/run_postgibbs.R b/R/run_postgibbs.R
index a620ae7..9bbb0ce 100644
--- a/R/run_postgibbs.R
+++ b/R/run_postgibbs.R
@@ -77,14 +77,20 @@ run_postgibbs <- function(path_2_execs,
   # Remove the temporary file
   unlink(temp_input_file)
 
+  # Move result files to output directory if different from input
+  if(output_files_dir != input_files_dir){
+    files_res <- c(base::list.files(pattern = "^post"), "run_postgibbs.log")
+    for(i in seq_along(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir, files_res[i]))
+  }
+
   # Debugging: Print the output
   if(verbose) cat("Command output:", output, "\n")
 
   # Capture and print the log file content
   if(verbose){
-    if (file.exists("run_postgibbs.log")) {
+    if (file.exists(file.path(output_files_dir,"run_postgibbs.log"))) {
       cat("Log file content:\n")
-      cat(readLines("run_postgibbs.log"), sep = "\n")
+      cat(readLines(file.path(output_files_dir,"run_postgibbs.log")), sep = "\n")
     } else {
       cat("Log file not created.\n")
     }
diff --git a/R/run_predict.R b/R/run_predict.R
index 9ff5c81..56f2d2e 100644
--- a/R/run_predict.R
+++ b/R/run_predict.R
@@ -35,15 +35,15 @@ run_predict <- function(path_2_execs) {
   }
 
   # Check if executable exists
-  if (!file.exists(paste0(path_2_execs, predict))) {
-    stop("Executable not found at: ", paste0(path_2_execs, predict))
+  if (!file.exists(file.path(path_2_execs, predict))) {
+    stop("Executable not found at: ", file.path(path_2_execs, predict))
   }
   if (!file.exists("renf90.par")) {
     stop("Parameter file not found: renf90.par")
   }
 
   # Construct the command
-  command_predict <- paste0(path_2_execs, predict, " renf90.par")
+  command_predict <- paste0(file.path(path_2_execs, predict), " renf90.par")
 
   # Run the command and log the output
   output <- execute_command(command = command_predict, logfile = "run_predict.log")
diff --git a/R/run_renum.R b/R/run_renum.R
index faba8f7..7076b44 100644
--- a/R/run_renum.R
+++ b/R/run_renum.R
@@ -13,7 +13,8 @@
 #'
 #' \donttest{
 #'  #run_renum(path_2_execs = "path/bf90_execs/",
-#'  #input_files_dir = "weight_2022_no_cov_cv.par")
+#'  #          raw_par_file = "my_analysis.par",
+#'  #          output_files_dir = "results")
 #' }
 #'
 #'
@@ -22,10 +23,10 @@ run_renum <- function(path_2_execs = ".",
                       raw_par_file = NULL,
                       output_files_dir = NULL,
                       verbose = TRUE) {
-
+  
   # Checks
   if(is.null(output_files_dir)) stop("Define a output directory in output_files_dir")
-
+  
   if(file.exists(output_files_dir)){
     check_files <- list.files(output_files_dir)
     if(length(check_files) > 0) warning(paste("Directory", output_files_dir, "is not empty. Some files may be replaced."))
@@ -33,52 +34,109 @@ run_renum <- function(path_2_execs = ".",
   } else {
     stop(paste("Directory '", output_files_dir, "' does not exist. Create it before running the function."))
   }
-
+  
   path_2_execs <- normalizePath(path_2_execs)
-
+  
   if(is.null(raw_par_file)) stop("Define raw_par_file.")
   raw_file <- normalizePath(raw_par_file)
   input_files_dir <- dirname(raw_file)
-
+  
   if (!file.exists(raw_file)) {
     stop("Parameter file not found at: ", raw_file)
   }
-
+  
+  #### Check par files ####
+  parfile <- base::readLines(raw_file)
+  
+  ## Check DATAFILE
+  datafile <- parfile[grep("DATAFILE", parfile) + 1]
+  if(length(datafile) == 0) stop(paste("DATAFILE was not defined in the parameter file", raw_file))
+  
+  # Check if file exist
+  if(!file.exists(file.path(input_files_dir, datafile))) stop(paste("File:", datafile, "defined in the file", raw_file,"line", grep("DATAFILE", parfile) + 1,"does not exist."))
+  
+  ## Check pedfile
+  pedfile <- parfile[which(grepl("^FILE$", parfile) | grepl("^FILE ", parfile) | grepl("^FILE#", parfile)) + 1]
+  
+  if(length(pedfile) != 0){
+    pedfile <- pedfile[1]   # first match only -> scalar (avoids length>1 in conditions)
+    # Check if file exist
+    if(!file.exists(file.path(input_files_dir, pedfile))) stop(paste("File:", pedfile, "defined in the file", raw_file,"line", which(grepl("^FILE$", parfile) | grepl("^FILE ", parfile) | grepl("^FILE#", parfile)) + 1,"does not exist."))
+  } else pedfile <- "not defined"
+  
+  ## Check snpfile
+  snpfile <- parfile[grep("SNP_FILE", parfile) + 1]
+  
+  if(length(snpfile) != 0){
+    snpfile <- snpfile[1]   # first match only -> scalar (avoids length>1 in conditions)
+    # Check if file exist
+    if(!file.exists(file.path(input_files_dir, snpfile))) stop(paste("File:", snpfile, "defined in the file", raw_file,"line", grep("SNP_FILE", parfile) + 1, "does not exist."))
+  } else snpfile <- "not defined"
+  
+  if(verbose){
+    cat(paste("DATAFILE:", datafile, "\n"))
+    cat(paste("FILE (pedigree file):", pedfile, "\n"))
+    cat(paste("SNP_FILE:", snpfile, "\n"))
+  }
+  
+  #####
+  
   #Assign .exe or not based on OS
   if (.Platform$OS.type == "unix") {
     renum = "renumf90"
   } else if (.Platform$OS.type == "windows") {
     renum = "renumf90.exe"
   }
-
+  
   cur_dir <- getwd() # save working directory location
-
+  
   # Construct the command
   command_renum <- paste0(file.path(path_2_execs, renum)," ", paste0("'",raw_file, "'"))
-
+  
   # Check if executable and parameter files exist
   if (!file.exists(file.path(path_2_execs, renum))) {
-    stop("Executable not found at: ", paste0(path_2_execs, renum))
+    stop("Executable not found at: ", paste0("'",path_2_execs, renum,"'"))
   }
-
+  
   # Run the command and log the output
   setwd(input_files_dir)
   output <- execute_command(command = command_renum, logfile = "run_renum.log")
-
+  
   # Capture and print the log file content
+  result <- if(file.exists("run_renum.log")) readLines("run_renum.log") else character()
+  
   if(verbose){
     if (file.exists("run_renum.log")) {
       cat("Log file content:\n")
-      cat(readLines("run_renum.log"), sep = "\n")
+      cat(result, sep = "\n")
     } else {
       cat("Log file not created.\n")
     }
   }
-
-  # Move generated files to working directory
-  files_res <- c("renadd03.ped", "renf90.dat", "renf90.fields", "renf90.inb", "renf90.par" ,"renf90.tables", "run_renum.log")
-  for(i in 1:length(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
-
+  
+  if(any(grepl("EFFECT: not found", result))) {
+    files_res <- c("run_renum.log")
+    file.rename(from = files_res, to = file.path(output_files_dir,files_res))
+    stop(paste("Effect not found. Check ",file.path(output_files_dir,files_res)))
+  } else {
+    # Move generated files to working directory
+    renadd <- base::list.files(pattern = "^renadd[0-9]+\\.ped$")   # renumf90 names ped by effect position
+    if (file.exists("renf90.inb")) {   # renumf90 just wrote it in the current dir; files not moved yet
+      files_res <- c(renadd, "renf90.dat", "renf90.fields", "renf90.inb", "renf90.par", "renf90.tables", "run_renum.log")
+    } else {
+      files_res <- c(renadd, "renf90.dat", "renf90.fields", "renf90.par" ,"renf90.tables", "run_renum.log")
+    }
+    
+    for(i in 1:length(files_res)) file.rename(from = files_res[i], to = file.path(output_files_dir,files_res[i]))
+  }
+  
+  # copy snp file
+  if(snpfile != "not defined" && !file.exists(file.path(output_files_dir, basename(snpfile)))){
+    snpfile <- normalizePath(snpfile)
+    file.symlink(base::file.path(snpfile), base::file.path(output_files_dir, basename(snpfile)))
+    file.symlink(paste0(base::file.path(snpfile),"_XrefID"), base::file.path(output_files_dir, paste0(basename(snpfile), "_XrefID")))
+  }
+  
   # Return to past working directory
   setwd(cur_dir)
 }
diff --git a/R/utils.R b/R/utils.R
index f55150f..a7dcdf5 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -2,7 +2,8 @@
 #' 
 #' @param command comment line used to run executable file
 #' @param logfile logfile name
-#' 
+#' @noRd
+#'
 execute_command <- function(command, logfile) {
   if (.Platform$OS.type == "unix") {
     output <- system(paste(command, "2>&1 | tee -a", logfile), intern = TRUE)
@@ -15,6 +16,7 @@ execute_command <- function(command, logfile) {
 }
 
 #' 
+#' @noRd
 create_folds <- function(data, num_folds) {
   n <- base::nrow(data)
   fold_size <- n %/% num_folds
@@ -28,15 +30,15 @@ create_folds <- function(data, num_folds) {
 }
 
 # Function to mutate phenotypes to missing value for given folds
-mutate_folds <- function(phenos, folds, num_folds, missing_value_code) {
+mutate_folds <- function(phenos, folds, num_folds, missing_value_code, id_col) {
   base::lapply(1:num_folds, function(i) {
     phenos %>%
-      dplyr::mutate(V2 = base::ifelse(V5 %in% folds[[i]], missing_value_code, V2))
+      dplyr::mutate(V2 = base::ifelse(phenos[[id_col]] %in% base::unlist(folds[[i]]), missing_value_code, V2))  # id_col = position of id column; unlist -> compare IDs, not a data.frame
   })
 }
 
 # Create cross-validation datasets
-create_cv_datasets <- function(run, fold, data_frame, dir_path, renf90, renf90_ped_name, input_files_dir, snp_file_name) {
+create_cv_datasets <- function(run, fold, data_frame, dir_path, renf90, renf90_ped_name, input_files_dir, reuse_files) {
   file_name <- base::sprintf("renf90_run%d_fold%d.dat", run, fold)
   utils::write.table(data_frame, file = file_name, sep = " ", row.names = FALSE, col.names = FALSE, quote = FALSE)
   
@@ -51,12 +53,12 @@ create_cv_datasets <- function(run, fold, data_frame, dir_path, renf90, renf90_p
   # Create symbolic links instead of copying files
   file.symlink(base::file.path(renf90_ped_name), base::file.path(dir_path, basename(renf90_ped_name)))
   file.symlink(base::file.path(input_files_dir, "renf90.fields"), base::file.path(dir_path, "renf90.fields"))
-  file.symlink(base::file.path(input_files_dir, "renf90.inb"), base::file.path(dir_path, "renf90.inb"))
+  if (base::file.exists(base::file.path(input_files_dir, "renf90.inb")))   # renf90.inb is optional
+    file.symlink(base::file.path(input_files_dir, "renf90.inb"), base::file.path(dir_path, "renf90.inb"))
   file.symlink(base::file.path(input_files_dir, "renf90.tables"), base::file.path(dir_path, "renf90.tables"))
   
-  if (!is.null(snp_file_name)) {
-    Xref_file <- paste0(snp_file_name, "_XrefID")
-    file.symlink(base::file.path(snp_file_name), base::file.path(dir_path, basename(snp_file_name)))
-    file.symlink(base::file.path(Xref_file), base::file.path(dir_path, basename(Xref_file)))
+  # Genomic: link the pre-built G-inverse (Gi) + cleaned SNPs so the fold reuses them
+  if (!is.null(reuse_files)) {
+    for (f in reuse_files) file.symlink(f, base::file.path(dir_path, base::basename(f)))
   }
 }
diff --git a/R/write_par.R b/R/write_par.R
new file mode 100644
index 0000000..d829caf
--- /dev/null
+++ b/R/write_par.R
@@ -0,0 +1,201 @@
+#' Write a RENUMF90 parameter (.par) file
+#'
+#' This function writes the raw RENUMF90 parameter (.par) file.
+#'
+#' This function assembles the raw parameter file that you feed to \code{\link{run_renum}}
+#' (i.e. the file you would otherwise write by hand), NOT the renf90.par file that
+#' RENUMF90 produces. It is built from the standard RENUMF90 blocks in order:
+#' DATAFILE, TRAITS, FIELDS_PASSED TO OUTPUT, WEIGHT(S), RESIDUAL_VARIANCE, the EFFECT
+#' blocks, the random-effect block, an optional SNP_FILE, an optional (CO)VARIANCES block
+#' and any OPTION lines.
+#'
+#' Effects are supplied as a named list. Each list name is the effect's label, written as the
+#' "#..." comment on its EFFECT line (a per-effect \code{comment} overrides it). Inside each
+#' element give the data-file column with
+#' \code{col} (\code{pos} is accepted as an alias), the \code{type} ("cross" or "cov") and, for
+#' cross-classified effects, the \code{class} ("numer" or "alpha"). For multi-trait models \code{col}
+#' may be a vector with one column per trait; a single value is repeated across traits (so two
+#' traits give, e.g., \code{2 2 cross alpha}).
+#'
+#' Mark the random effect with \code{random}, giving its column number (e.g. \code{random = 1})
+#' or its label. Because RENUMF90 requires RANDOM to directly follow its effect, the random effect
+#' is always written last (it is moved to the end automatically, with a message, if you list it
+#' earlier). RANDOM (with \code{random_type}) and an optional OPTIONAL line (e.g.
+#' \code{optional = "pe"} for a permanent-environment / repeatability term) belong to that effect
+#' block; they are followed by the separate pedigree block FILE / FILE_POS, then an optional
+#' SNP_FILE and an optional (CO)VARIANCES block. The order of OPTION lines is irrelevant to RENUMF90.
+#'
+#' @param file path/name of the .par file to create. This field should be in quotes "".
+#' @param datafile name of the data file written under DATAFILE. This field should be in quotes "".
+#' @param traits column position(s) of the trait(s). For multi-trait give a vector, e.g. c(6, 7).
+#' @param effects named list of effects, in model order. Each list name is the effect's label
+#'   (written as "EFFECT #name"). Inside each
+#'   element supply \code{col} (the effect's data-file column; \code{pos} is an alias, and a vector
+#'   gives one column per trait), \code{type} ("cross"/"cov") and, for cross effects, \code{class}
+#'   ("numer"/"alpha"). An optional \code{comment} overrides the "#..." label.
+#' @param random which effect is the random effect, given as its column number (e.g.
+#'   \code{random = 1}) or its label. Matching is by \code{col}, so the effect's "#..." label can
+#'   be anything you like. This effect is written last (RENUMF90 requires RANDOM to follow it),
+#'   with RANDOM immediately after it. Leave NULL for a fixed-effects-only model.
+#' @param random_effect deprecated alias for \code{random} (name-based); use \code{random}.
+#' @param pedigree_file pedigree file written under FILE. Required when a random effect is set.
+#' @param pedigree_pos integer vector written under FILE_POS (animal, sire, dam, then two 0s).
+#'   Defaults to c(1, 2, 3, 0, 0).
+#' @param random_type keyword written under RANDOM. Defaults to "animal".
+#' @param optional keyword(s) written under an OPTIONAL line after RANDOM (e.g. "pe" for a
+#'   permanent-environment / repeatability term). NULL to omit.
+#' @param residual_variance starting residual (co)variance written under RESIDUAL_VARIANCE. A
+#'   scalar for one trait, or a matrix for multi-trait (one row per line).
+#' @param covariances starting (co)variance for the random effect, written under (CO)VARIANCES.
+#'   Scalar or matrix. NULL to omit.
+#' @param snp_file genotype file written under SNP_FILE (inside the random-effect block). Leave
+#'   NULL for a pedigree-only model.
+#' @param missing_value value written as "OPTION missing <value>". Leave NULL to omit.
+#' @param options character vector of extra OPTION lines given WITHOUT the leading "OPTION "
+#'   (e.g. c("callrate 0.9", "cat 0 0")).
+#' @param fields_passed column position(s) written under FIELDS_PASSED TO OUTPUT (NULL = blank).
+#' @param weights column position written under WEIGHT(S) (NULL = blank).
+#' @param overwrite logical; if FALSE (default) the function stops when \code{file} already exists.
+#'
+#' @return (invisibly) the character vector of lines written to \code{file}.
+#' @examples
+#'
+#' \donttest{
+#'  # two-trait repeatability model: animal + permanent environment, genomic
+#'  # write_par(file     = "act.par",
+#'  #           datafile = "act_act_auc.csv",
+#'  #           traits   = c(6, 7),
+#'  #           residual_variance = diag(2),
+#'  #           effects  = list(lact_group         = list(col = 2, class = "alpha"),
+#'  #                           breed_group        = list(col = 3, class = "alpha"),
+#'  #                           DIM_group          = list(col = 4, class = "alpha"),
+#'  #                           contemporary_group = list(col = 5, class = "alpha"),
+#'  #                           animal             = list(col = 1, class = "alpha")),
+#'  #           random        = 1,          # the animal effect is column 1 (label can be anything)
+#'  #           optional      = "pe",
+#'  #           pedigree_file = "ponderosa.ped",
+#'  #           snp_file      = "ponderosa_bulls.geno",
+#'  #           covariances   = matrix(c(1, 0.1, 0.1, 1), 2, 2),
+#'  #           options       = "cat 0 0")
+#' }
+#'
+#' @export
+write_par <- function(file = NULL,
+                      datafile = NULL,
+                      traits = NULL,
+                      effects = NULL,
+                      random = NULL,
+                      random_effect = NULL,
+                      pedigree_file = NULL,
+                      pedigree_pos = c(1, 2, 3, 0, 0),
+                      random_type = "animal",
+                      optional = NULL,
+                      residual_variance = NULL,
+                      covariances = NULL,
+                      snp_file = NULL,
+                      missing_value = NULL,
+                      options = NULL,
+                      fields_passed = NULL,
+                      weights = NULL,
+                      overwrite = FALSE) {
+
+  # Checks
+  if(is.null(file))              stop("Define the output .par file name in 'file'.")
+  if(is.null(datafile))          stop("Define the data file name in 'datafile'.")
+  if(is.null(traits))            stop("Define the trait column position(s) in 'traits'.")
+  if(is.null(residual_variance)) stop("Define a starting value in 'residual_variance'.")
+  if(is.null(effects) || !is.list(effects) || length(effects) == 0)
+    stop("'effects' must be a non-empty named list (see ?write_par).")
+  if(is.null(names(effects)) || any(names(effects) == ""))
+    stop("Every element of 'effects' must be named; the names are used as the effect labels.")
+  rnd <- if(!is.null(random)) random else random_effect      # random effect: its col number or its label
+  if(!is.null(rnd) && is.null(pedigree_file))
+    stop("'pedigree_file' is required when a random effect is set (see 'random').")
+  if(!overwrite && file.exists(file))
+    stop(paste("File", file, "already exists. Set overwrite = TRUE to replace it."))
+
+  n_traits <- length(traits)
+
+  ## Format a scalar/vector/matrix (co)variance as one line per matrix row
+  var_lines <- function(x) {
+    if(is.matrix(x)) apply(x, 1, function(r) paste(r, collapse = " ")) else paste(x, collapse = " ")
+  }
+
+  ## Format one EFFECT definition line: one column per trait, then cross/cov [class]
+  effect_line <- function(e, cols, nm) {
+    type <- if(is.null(e[["type"]])) "cross" else e[["type"]]
+    if(!type %in% c("cross", "cov")) stop(paste0("Effect '", nm, "': 'type' must be 'cross' or 'cov'."))
+    posdef <- paste(cols, collapse = " ")
+    if(type == "cov") return(paste(posdef, "cov"))
+    cls <- if(is.null(e[["class"]])) "numer" else e[["class"]]
+    if(!cls %in% c("numer", "alpha")) stop(paste0("Effect '", nm, "': 'class' must be 'numer' or 'alpha'."))
+    paste(posdef, "cross", cls)
+  }
+
+  # Assemble the file, block by block, in RENUMF90 order
+  lines <- c("DATAFILE", datafile,
+             "TRAITS",   paste(traits, collapse = " "),
+             "FIELDS_PASSED TO OUTPUT", if(is.null(fields_passed)) "" else paste(fields_passed, collapse = " "),
+             "WEIGHT(S)",              if(is.null(weights))       "" else paste(weights, collapse = " "),
+             "RESIDUAL_VARIANCE",      var_lines(residual_variance))
+
+  # Resolve each effect's data-file column(s) and find the random effect. 'col' (or 'pos' alias)
+  # may be a vector, one column per trait; a single value is repeated across traits.
+  cols   <- list()
+  rnd_nm <- NULL
+  for(nm in names(effects)) {
+    e   <- base::as.list(effects[[nm]])
+    col <- if(!is.null(e[["col"]])) e[["col"]] else e[["pos"]]                  # data-file column ('pos' = alias)
+    if(is.null(col))
+      stop(paste0("Effect '", nm, "': supply its data-file column with col = (e.g. col = 1)."))
+    col <- suppressWarnings(as.integer(col))
+    if(any(is.na(col)) || any(col < 1))
+      stop(paste0("Effect '", nm, "': 'col' must be a positive integer (one per trait)."))
+    if(length(col) == 1 && n_traits > 1) col <- rep(col, n_traits)             # same column across traits
+    if(length(col) != n_traits)
+      stop(paste0("Effect '", nm, "': give one column (recycled) or one per trait (", n_traits, ")."))
+    cols[[nm]] <- col
+    # match by col number (label-independent) or by name (alias); keep the first match
+    if(!is.null(rnd) && is.null(rnd_nm) &&
+       ((is.numeric(rnd) && rnd %in% col) || identical(as.character(rnd), nm))) rnd_nm <- nm
+  }
+  if(!is.null(rnd) && is.null(rnd_nm))
+    stop(paste0("No effect matches 'random' = ", rnd,
+                ". Give the col number (or label) of the random effect."))
+
+  # RENUMF90 requires RANDOM to directly follow its effect, so the random effect is written last.
+  ord <- names(effects)
+  if(!is.null(rnd_nm)) {
+    if(ord[length(ord)] != rnd_nm)
+      message("write_par: moving random effect '", rnd_nm, "' to the end (RANDOM must follow it).")
+    ord <- c(setdiff(ord, rnd_nm), rnd_nm)
+  }
+
+  # EFFECT blocks in order; RANDOM (+ OPTIONAL) belong to the random effect, then the separate
+  # pedigree block FILE / FILE_POS, then SNP_FILE and (CO)VARIANCES.
+  snp_done <- FALSE
+  for(nm in ord) {
+    e      <- base::as.list(effects[[nm]])
+    cmt    <- if(is.null(e[["comment"]])) nm else e[["comment"]]                # "#..." label (defaults to the name)
+    header <- if(is.na(cmt) || !nzchar(cmt)) "EFFECT" else paste0("EFFECT #", cmt)
+    lines  <- c(lines, header, effect_line(e, cols[[nm]], nm))
+
+    if(identical(nm, rnd_nm)) {
+      lines <- c(lines, "RANDOM", random_type)                                 # RANDOM (+ OPTIONAL): part of the effect block
+      if(!is.null(optional)) lines <- c(lines, "OPTIONAL", paste(optional, collapse = " "))
+      lines <- c(lines, "FILE", pedigree_file, "FILE_POS", paste(pedigree_pos, collapse = " "))  # separate pedigree block
+      if(!is.null(snp_file))    { lines <- c(lines, "SNP_FILE", snp_file); snp_done <- TRUE }
+      if(!is.null(covariances))   lines <- c(lines, "(CO)VARIANCES", var_lines(covariances))
+    }
+  }
+
+  # SNP_FILE for a genomic model with no random pedigree effect (rare)
+  if(!is.null(snp_file) && !snp_done) lines <- c(lines, "SNP_FILE", snp_file)
+
+  # OPTION lines (order is irrelevant to RENUMF90)
+  if(!is.null(missing_value)) lines <- c(lines, paste("OPTION missing", missing_value))
+  if(!is.null(options))       lines <- c(lines, paste("OPTION", options))
+
+  base::writeLines(lines, con = file)
+  invisible(lines)
+}
diff --git a/README.md b/README.md
index c8caa83..0bf748f 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,20 @@
 # BIGf90  
 ## (B)reeding (I)nsight (G)enomics f90
+
+<!-- badges: start -->
+[![Development Status](https://img.shields.io/badge/status-active%20development-yellow)](https://github.com/josuechinchilla/BIGf90)
+[![R-CMD-check](https://github.com/josuechinchilla/BIGf90/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/josuechinchilla/BIGf90/actions/workflows/R-CMD-check.yaml)
+[![R](https://img.shields.io/badge/R-%3E%3D%204.3-blue)](https://www.r-project.org/)
+[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://www.apache.org/licenses/LICENSE-2.0)
+[![GitHub issues](https://img.shields.io/github/issues/josuechinchilla/BIGf90)](https://github.com/josuechinchilla/BIGf90/issues)
+[![GitHub pull requests](https://img.shields.io/github/issues-pr/josuechinchilla/BIGf90)](https://github.com/josuechinchilla/BIGf90/pulls)
+[![GitHub Release](https://img.shields.io/github/v/release/josuechinchilla/BIGf90?include_prereleases)](https://github.com/josuechinchilla/BIGf90/releases/latest)
+<!-- badges: end -->
   
 BIGf90 is a wrapper package for BLUPf90 family of programs, please use the following citations:
 
 #### BIGf90 Reference:
-  Chinchilla-Vargas J, Taniguti C, Sandercock A, Breeding Insight Team (2024). BIGf90: Breeding Insight Genomics R front face to blupf90 modules. R package version 0.3.2, https://github.com/josuechinchilla/BIGf90
+  Chinchilla-Vargas J, Taniguti C, Sandercock A, Breeding Insight Team (2024). BIGf90: Breeding Insight Genomics R front face to blupf90 modules. R package version 0.4.0, https://github.com/josuechinchilla/BIGf90
   
 #### BLUPF90 Reference:
   Misztal, I., S. Tsuruta, D.A.L. Lourenco, I. Aguilar, A. Legarra, and Z. Vitezica. 2014. Manual for BLUPF90 family of programs: http://nce.ads.uga.edu/wiki/lib/exe/fetch.php?media=blupf90_all2.pdf")
@@ -14,6 +24,7 @@ please contact Ignacy Misztal (ignacy@uga.edu) or Daniela Lourenco (danilino@uga
   
 ### Overview
 This package has functions to:
+* Write RENUMF90 parameter (.par) files with write_par().
 * Process parameter files with Renumf90.
 * Calculate EBVs through BLUP, GBLUP,ssGBLUP with Blupf90+
 * Calculate adjusted phenotypes using Predictf90
diff --git a/inst/no_cov_script.R b/inst/no_cov_script.R
deleted file mode 100644
index 1c97d98..0000000
--- a/inst/no_cov_script.R
+++ /dev/null
@@ -1,47 +0,0 @@
-#Set seed, working directory and load packages needed
-set.seed(101919)
-library("BIGf90")
-setwd("~/Desktop/weight_length_2024/weight_2022phenos_2022ped/no_cov")
-#set path to executables
-#path_2_execs = system.file("exec_windows",package = "BIGf90")
-path_2_execs = system.file("exec_mac",package = "BIGf90")
-
-####    1) Process raw parameter file with renumf90    ####
-run_renum(path_2_execs,
-          output_files_dir = "results",
-          raw_par_file = "weight22_ped22.par")
-
-####    2) Estimate  variance components with gibbsf90+ and postgibbsf90    ####
-
-run_gibbs(path_2_execs,
-          input_files_dir = "results",
-          gibbs_iter = 250000,
-          gibbs_burn = 50000,
-          gibbs_keep = 1)
-
-run_postgibbs(path_2_execs, input_files_dir = "results",
-              postgibbs_burn = 1,
-              postgibbs_keep = 100)
-
-####    3) Run a cross-validation analysis (CVA)  to estimate the accuracy of your EBV predictions  ####
-
-bf90_cv(
-  path_2_execs,
-  missing_value_code = -999,
-  random_effect_col = 3,
-  h2 = 0.95,
-  num_runs = 10,
-  num_folds = 5,
-  renf90_ped_name = "results/renadd03.ped",
-  input_files_dir = "results",
-  output_files_dir = "bf90_cv_results/",
-  output_table_name = "CVA_results_weight_2022phenos_2022ped_no_cov.txt"
-  )
-
-# Check in your version if the files generated have the same results below (with some tolerance)
-yhat <- read.table("bf90_cv_results/yhat_residual")
-dim(yhat)
-sum(yhat$V2) # 89.84
-
-cv_results <- read.table("bf90_cv_results/CVA_results_weight_2022phenos_2022ped_no_cov.txt", skip = 1, nrows = 30)
-sum(cv_results$V4) # 25.746
diff --git a/man/bf90_cv.Rd b/man/bf90_cv.Rd
index cf70c6c..6dd628f 100644
--- a/man/bf90_cv.Rd
+++ b/man/bf90_cv.Rd
@@ -14,8 +14,6 @@ bf90_cv(
   input_files_dir = ".",
   output_files_dir = ".",
   output_table_name = NULL,
-  renf90_ped_name = NULL,
-  snp_file_name = NULL,
   seed = 101919,
   verbose = TRUE
 )
@@ -39,10 +37,6 @@ bf90_cv(
 
 \item{output_table_name}{Name of the final tab-separated out-up file. This field should be in quotes "".}
 
-\item{renf90_ped_name}{Name of pedigree file generated by renumf90. This field should be in quotes "".}
-
-\item{snp_file_name}{the name of the genotype file to be used. If use, this field should be in quotes"". Default for the function is no genotype file.}
-
 \item{seed}{set seed for the stochastic process}
 
 \item{verbose}{logical defining if information will be printed on the console}
@@ -56,7 +50,7 @@ This function runs a K-fold cross-validation analysis (CVA) using blupf90 module
 \details{
 This function sets up and runs a K-fold cross-validation analysis (CVA) using blupf90+ and predictf90.
 The function run_renumf90 needs to be used beforehand to process a .par file created by the user.
-This function calculates 2 accuracy estimates: correlation between raw phenotypes and ebvs divided by the square-root of narrow sense heritabilty and correlation between corrected phenotypes and ebvs along with bias estimations calculated as the regression of the phenotypes on the ebvs.
+Using the phenotype corrected for all fixed effects (y*), the function reports predictive ability as the correlation between y* and the ebvs, accuracy as that correlation divided by the square-root of the narrow-sense heritability, and bias as the regression of y* on the ebvs.
 }
 \examples{
 ## Example for a CVA with 5 independent runs dividing the data in 10 folds.
@@ -69,8 +63,8 @@ This function calculates 2 accuracy estimates: correlation between raw phenotype
 #      num_runs = 5,
 #      num_folds = 10,
 #      output_table_name = "example_run",
-#     renf90_ped_name = "renadd03.ped",
-#     snp_file_name = "my_genos.geno" )
+#      input_files_dir = ".",
+#      output_files_dir = ".")
 
 
 }
diff --git a/man/execute_command.Rd b/man/execute_command.Rd
deleted file mode 100644
index b79f891..0000000
--- a/man/execute_command.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{execute_command}
-\alias{execute_command}
-\title{Function to run commands on the terminal and log output}
-\usage{
-execute_command(command, logfile)
-}
-\arguments{
-\item{command}{comment line used to run executable file}
-
-\item{logfile}{logfile name}
-}
-\description{
-Function to run commands on the terminal and log output
-}
diff --git a/man/run_renum.Rd b/man/run_renum.Rd
index 15af7f1..eaec6c5 100644
--- a/man/run_renum.Rd
+++ b/man/run_renum.Rd
@@ -31,7 +31,8 @@ The outputs will be the standard output files produced by renumf90. A log file c
 
 \donttest{
  #run_renum(path_2_execs = "path/bf90_execs/",
- #input_files_dir = "weight_2022_no_cov_cv.par")
+ #          raw_par_file = "my_analysis.par",
+ #          output_files_dir = "results")
 }
 
 
diff --git a/man/write_par.Rd b/man/write_par.Rd
new file mode 100644
index 0000000..8c6c01b
--- /dev/null
+++ b/man/write_par.Rd
@@ -0,0 +1,129 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/write_par.R
+\name{write_par}
+\alias{write_par}
+\title{Write a RENUMF90 parameter (.par) file}
+\usage{
+write_par(
+  file = NULL,
+  datafile = NULL,
+  traits = NULL,
+  effects = NULL,
+  random = NULL,
+  random_effect = NULL,
+  pedigree_file = NULL,
+  pedigree_pos = c(1, 2, 3, 0, 0),
+  random_type = "animal",
+  optional = NULL,
+  residual_variance = NULL,
+  covariances = NULL,
+  snp_file = NULL,
+  missing_value = NULL,
+  options = NULL,
+  fields_passed = NULL,
+  weights = NULL,
+  overwrite = FALSE
+)
+}
+\arguments{
+\item{file}{path/name of the .par file to create. This field should be in quotes "".}
+
+\item{datafile}{name of the data file written under DATAFILE. This field should be in quotes "".}
+
+\item{traits}{column position(s) of the trait(s). For multi-trait give a vector, e.g. c(6, 7).}
+
+\item{effects}{named list of effects, in model order. Each list name is the effect's label
+(written as "EFFECT #name"). Inside each
+element supply \code{col} (the effect's data-file column; \code{pos} is an alias, and a vector
+gives one column per trait), \code{type} ("cross"/"cov") and, for cross effects, \code{class}
+("numer"/"alpha"). An optional \code{comment} overrides the "#..." label.}
+
+\item{random}{which effect is the random effect, given as its column number (e.g.
+\code{random = 1}) or its label. Matching is by \code{col}, so the effect's "#..." label can
+be anything you like. This effect is written last (RENUMF90 requires RANDOM to follow it),
+with RANDOM immediately after it. Leave NULL for a fixed-effects-only model.}
+
+\item{random_effect}{deprecated alias for \code{random} (name-based); use \code{random}.}
+
+\item{pedigree_file}{pedigree file written under FILE. Required when a random effect is set.}
+
+\item{pedigree_pos}{integer vector written under FILE_POS (animal, sire, dam, then two 0s).
+Defaults to c(1, 2, 3, 0, 0).}
+
+\item{random_type}{keyword written under RANDOM. Defaults to "animal".}
+
+\item{optional}{keyword(s) written under an OPTIONAL line after RANDOM (e.g. "pe" for a
+permanent-environment / repeatability term). NULL to omit.}
+
+\item{residual_variance}{starting residual (co)variance written under RESIDUAL_VARIANCE. A
+scalar for one trait, or a matrix for multi-trait (one row per line).}
+
+\item{covariances}{starting (co)variance for the random effect, written under (CO)VARIANCES.
+Scalar or matrix. NULL to omit.}
+
+\item{snp_file}{genotype file written under SNP_FILE (inside the random-effect block). Leave
+NULL for a pedigree-only model.}
+
+\item{missing_value}{value written as "OPTION missing \if{html}{\out{<value>}}". Leave NULL to omit.}
+
+\item{options}{character vector of extra OPTION lines given WITHOUT the leading "OPTION "
+(e.g. c("callrate 0.9", "cat 0 0")).}
+
+\item{fields_passed}{column position(s) written under FIELDS_PASSED TO OUTPUT (NULL = blank).}
+
+\item{weights}{column position written under WEIGHT(S) (NULL = blank).}
+
+\item{overwrite}{logical; if FALSE (default) the function stops when \code{file} already exists.}
+}
+\value{
+(invisibly) the character vector of lines written to \code{file}.
+}
+\description{
+This function writes the raw RENUMF90 parameter (.par) file.
+}
+\details{
+This function assembles the raw parameter file that you feed to \code{\link{run_renum}}
+(i.e. the file you would otherwise write by hand), NOT the renf90.par file that
+RENUMF90 produces. It is built from the standard RENUMF90 blocks in order:
+DATAFILE, TRAITS, FIELDS_PASSED TO OUTPUT, WEIGHT(S), RESIDUAL_VARIANCE, the EFFECT
+blocks, the random-effect block, an optional SNP_FILE, an optional (CO)VARIANCES block
+and any OPTION lines.
+
+Effects are supplied as a named list. Each list name is the effect's label, written as the
+"#..." comment on its EFFECT line (a per-effect \code{comment} overrides it). Inside each
+element give the data-file column with
+\code{col} (\code{pos} is accepted as an alias), the \code{type} ("cross" or "cov") and, for
+cross-classified effects, the \code{class} ("numer" or "alpha"). For multi-trait models \code{col}
+may be a vector with one column per trait; a single value is repeated across traits (so two
+traits give, e.g., \code{2 2 cross alpha}).
+
+Mark the random effect with \code{random}, giving its column number (e.g. \code{random = 1})
+or its label. Because RENUMF90 requires RANDOM to directly follow its effect, the random effect
+is always written last (it is moved to the end automatically, with a message, if you list it
+earlier). RANDOM (with \code{random_type}) and an optional OPTIONAL line (e.g.
+\code{optional = "pe"} for a permanent-environment / repeatability term) belong to that effect
+block; they are followed by the separate pedigree block FILE / FILE_POS, then an optional
+SNP_FILE and an optional (CO)VARIANCES block. The order of OPTION lines is irrelevant to RENUMF90.
+}
+\examples{
+
+\donttest{
+ # two-trait repeatability model: animal + permanent environment, genomic
+ # write_par(file     = "act.par",
+ #           datafile = "act_act_auc.csv",
+ #           traits   = c(6, 7),
+ #           residual_variance = diag(2),
+ #           effects  = list(lact_group         = list(col = 2, class = "alpha"),
+ #                           breed_group        = list(col = 3, class = "alpha"),
+ #                           DIM_group          = list(col = 4, class = "alpha"),
+ #                           contemporary_group = list(col = 5, class = "alpha"),
+ #                           animal             = list(col = 1, class = "alpha")),
+ #           random        = 1,          # the animal effect is column 1 (label can be anything)
+ #           optional      = "pe",
+ #           pedigree_file = "ponderosa.ped",
+ #           snp_file      = "ponderosa_bulls.geno",
+ #           covariances   = matrix(c(1, 0.1, 0.1, 1), 2, 2),
+ #           options       = "cat 0 0")
+}
+
+}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..4897341
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+library(BIGf90)
+
+test_check("BIGf90")
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
new file mode 100644
index 0000000..039d9db
--- /dev/null
+++ b/tests/testthat/test-utils.R
@@ -0,0 +1,27 @@
+# Tests for internal helpers in utils.R.
+# create_folds() splits a data frame into contiguous folds; the last fold
+# absorbs the remainder when num_folds does not divide the row count evenly.
+# (Two-column frames are used so row subsetting keeps its data.frame shape.)
+
+test_that("create_folds splits rows into equal contiguous folds", {
+  d <- data.frame(id = 1:10, x = 11:20)
+  f <- create_folds(d, 5)
+  expect_length(f, 5)
+  expect_equal(vapply(f, nrow, integer(1)), rep(2L, 5))
+  expect_equal(sort(unlist(lapply(f, `[[`, "id"))), 1:10)   # every row kept exactly once (no dups/missing)
+})
+
+test_that("create_folds puts the remainder in the last fold", {
+  d <- data.frame(id = 1:11, x = 21:31)
+  f <- create_folds(d, 5)
+  expect_equal(vapply(f, nrow, integer(1)), c(2L, 2L, 2L, 2L, 3L))
+  expect_equal(sort(unlist(lapply(f, `[[`, "id"))), 1:11)
+})
+
+test_that("create_folds with a single fold returns all rows", {
+  d <- data.frame(id = 1:7, x = 8:14)
+  f <- create_folds(d, 1)
+  expect_length(f, 1)
+  expect_equal(nrow(f[[1]]), 7)
+  expect_equal(sort(f[[1]]$id), 1:7)
+})
diff --git a/tests/testthat/test-write_par.R b/tests/testthat/test-write_par.R
new file mode 100644
index 0000000..47fdc7d
--- /dev/null
+++ b/tests/testthat/test-write_par.R
@@ -0,0 +1,163 @@
+# Tests for write_par(): assembling RENUMF90 parameter (.par) files.
+# write_par() returns the written lines invisibly, so we test the return value
+# (and, in one case, the round-tripped file) instead of eyeballing text.
+
+test_that("single-trait genomic model produces the expected .par lines", {
+  tmp <- tempfile(fileext = ".par")
+  res <- write_par(
+    file     = tmp,
+    datafile = "sealice_match.txt",
+    traits   = 4,
+    residual_variance = 0.157,
+    effects  = list(mean   = list(col = 2, type = "cross", class = "numer"),
+                    tank   = list(col = 3, type = "cross", class = "numer"),
+                    animal = list(col = 1, type = "cross", class = "alpha")),
+    random        = 1,
+    pedigree_file = "sealice.ped",
+    snp_file      = "sealice_match.geno",
+    missing_value = -999,
+    options = c("callrate 0.9", "saveGInverse")
+  )
+
+  expected <- c(
+    "DATAFILE", "sealice_match.txt",
+    "TRAITS", "4",
+    "FIELDS_PASSED TO OUTPUT", "",
+    "WEIGHT(S)", "",
+    "RESIDUAL_VARIANCE", "0.157",
+    "EFFECT #mean", "2 cross numer",
+    "EFFECT #tank", "3 cross numer",
+    "EFFECT #animal", "1 cross alpha",
+    "RANDOM", "animal",
+    "FILE", "sealice.ped",
+    "FILE_POS", "1 2 3 0 0",
+    "SNP_FILE", "sealice_match.geno",
+    "OPTION missing -999",
+    "OPTION callrate 0.9",
+    "OPTION saveGInverse"
+  )
+  expect_identical(res, expected)             # returned lines
+  expect_identical(readLines(tmp), expected)  # round-trips to disk
+})
+
+test_that("random matches by column, independent of the effect label", {
+  res <- write_par(
+    file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+    residual_variance = 1,
+    effects = list(mu = list(col = 2, class = "numer"),
+                   id = list(col = 1, class = "alpha")),   # label 'id', not 'animal'
+    random = 1, pedigree_file = "ped.txt"
+  )
+  i <- match("EFFECT #id", res)               # the col-1 effect, whatever its label
+  expect_identical(res[i + 1], "1 cross alpha")
+  expect_identical(res[i + 2], "RANDOM")
+  expect_identical(res[i + 3], "animal")
+  expect_identical(res[match("FILE", res) + 1], "ped.txt")
+})
+
+test_that("a random effect not listed last is moved to the end (with a message)", {
+  expect_message(
+    res <- write_par(
+      file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+      residual_variance = 1,
+      effects = list(plant = list(col = 2, class = "numer"),   # random, listed FIRST
+                     sky   = list(col = 3, class = "numer"),
+                     fish  = list(col = 1, class = "alpha")),
+      random = 2, pedigree_file = "ped.txt"
+    ),
+    "moving random effect"
+  )
+  expect_identical(res[grepl("^EFFECT", res)],
+                   c("EFFECT #sky", "EFFECT #fish", "EFFECT #plant"))  # plant now last
+  i <- match("EFFECT #plant", res)
+  expect_identical(res[i + 1], "2 cross numer")
+  expect_identical(res[i + 2], "RANDOM")      # RANDOM directly follows its effect
+})
+
+test_that("'pos' is an alias for 'col'", {
+  a <- write_par(file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+                 residual_variance = 1,
+                 effects = list(mu = list(pos = 2, class = "numer"),
+                                id = list(pos = 1, class = "alpha")),
+                 random = 1, pedigree_file = "ped.txt")
+  b <- write_par(file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+                 residual_variance = 1,
+                 effects = list(mu = list(col = 2, class = "numer"),
+                                id = list(col = 1, class = "alpha")),
+                 random = 1, pedigree_file = "ped.txt")
+  expect_identical(a, b)
+})
+
+test_that("'comment' overrides the label on the EFFECT line", {
+  res <- write_par(file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+                   residual_variance = 1,
+                   effects = list(a = list(col = 1, class = "numer", comment = "my label")))
+  expect_true("EFFECT #my label" %in% res)
+})
+
+test_that("multi-trait: matrix variances, per-trait columns, OPTIONAL and (CO)VARIANCES", {
+  res <- write_par(
+    file = tempfile(fileext = ".par"), datafile = "act.csv", traits = c(6, 7),
+    residual_variance = matrix(c(1, 0, 0, 1), 2, 2),
+    effects = list(cg     = list(col = 5, class = "alpha"),
+                   animal = list(col = 1, class = "alpha")),
+    random = 1, optional = "pe", pedigree_file = "p.ped",
+    snp_file = "g.geno", covariances = matrix(c(1, 0.1, 0.1, 1), 2, 2),
+    options = "cat 0 0"
+  )
+  expect_identical(res[match("TRAITS", res) + 1], "6 7")
+  rv <- match("RESIDUAL_VARIANCE", res)
+  expect_identical(res[c(rv + 1, rv + 2)], c("1 0", "0 1"))   # one matrix row per line
+  expect_true("5 5 cross alpha" %in% res)                     # column repeated per trait
+  expect_true("1 1 cross alpha" %in% res)
+  expect_identical(res[match("RANDOM", res) + 1], "animal")
+  expect_identical(res[match("OPTIONAL", res) + 1], "pe")
+  cv <- match("(CO)VARIANCES", res)
+  expect_identical(res[c(cv + 1, cv + 2)], c("1 0.1", "0.1 1"))
+  expect_true("OPTION cat 0 0" %in% res)
+})
+
+test_that("fixed-effects-only model writes no RANDOM/FILE block", {
+  res <- write_par(file = tempfile(fileext = ".par"), datafile = "d.txt", traits = 1,
+                   residual_variance = 1,
+                   effects = list(mu = list(col = 2, class = "numer")))
+  expect_false(any(c("RANDOM", "FILE", "FILE_POS") %in% res))
+})
+
+test_that("informative errors on bad input", {
+  expect_error(                                              # missing column
+    write_par(file = tempfile(), datafile = "d", traits = 1, residual_variance = 1,
+              effects = list(a = list(class = "numer"))),
+    "data-file column"
+  )
+  expect_error(                                              # random matches nothing
+    write_par(file = tempfile(), datafile = "d", traits = 1, residual_variance = 1,
+              effects = list(a = list(col = 1, class = "numer")),
+              random = 9, pedigree_file = "p"),
+    "No effect matches"
+  )
+  expect_error(                                              # pedigree required with random
+    write_par(file = tempfile(), datafile = "d", traits = 1, residual_variance = 1,
+              effects = list(a = list(col = 1, class = "alpha")),
+              random = 1),
+    "pedigree_file"
+  )
+  expect_error(                                              # bad class
+    write_par(file = tempfile(), datafile = "d", traits = 1, residual_variance = 1,
+              effects = list(a = list(col = 1, class = "bogus"))),
+    "'class'"
+  )
+})
+
+test_that("refuses to overwrite an existing file unless overwrite = TRUE", {
+  tmp <- tempfile(fileext = ".par")
+  writeLines("x", tmp)                                       # file now exists
+  expect_error(
+    write_par(file = tmp, datafile = "d", traits = 1, residual_variance = 1,
+              effects = list(a = list(col = 1, class = "numer"))),
+    "already exists"
+  )
+  res <- write_par(file = tmp, datafile = "d", traits = 1, residual_variance = 1,
+                   effects = list(a = list(col = 1, class = "numer")), overwrite = TRUE)
+  expect_true("DATAFILE" %in% res)
+})