gromacs_funs.bash

run_gromacs_standard ()
{
# run_gromacs_standard {current run} {last run}
    if [ -f $1.gro ] ; then
        echo 'found $1.gro'
    else
        gmx grompp -f $1.mdp \
            -c $2 \
            -t $2 \
            -p topol_solvated.top \
            -o $1.tpr
        gmx mdrun -nt 10 -deffnm $1 -v
        gmx energy -f $1.edr -o $1_energy.xvg <<< $'Potential\nKinetic\nTotal\nTemperature\nPressure\nDensity\nVolume\n\n'
    fi
}

continue_gromacs_standard () 
{
# run_gromacs_standard {current run} {last run}
# skips run if {current_run.gro} is newer than {last_run.gro}
    if [ -f $1.gro ] ; then
        echo 'found $1.gro'
        if [ $1.gro -ot $2.mdp ] ; then
            echo 'found $1 (current set) to be older than $2 (past set), re-running.....'
            run_gromacs_standard $1 $2
        else
            if [ $1.gro -ot $1.mdp ] ; then
                echo 'found $1.gro to be older than $1.mdp, re-running.....'
                run_gromacs_standard $1 $2
            else
                echo 'found $1 to be newer than $2 and $1.mdp, skipping step'
            fi
        fi
    else
        echo 'could not find $1, running'
        run_gromacs_standard $1 $2
    fi
}

make_movie_standard ()
{
    SKIPS=${2:-1}
    SUBSYSTEM=${3:-'0'}
    gmx trjconv -f $1.xtc -o $1.pdb -s $1.tpr -skip $SKIPS -pbc mol <<< $SUBSYSTEM
}

make_standard_energy ()
{
    gmx make_ndx -f $1.gro -o $1_index.ndx <<< $'a OW\na HW*\na C\na O\na CH*\n q\n'
    gmx rdf -n $1_index.ndx -o RDF_C_OW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'C\nOW\n'
    gmx rdf -n $1_index.ndx -o RDF_C_HW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'C\nHW\n'
    gmx rdf -n $1_index.ndx -o RDF_O_OW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'O\nOW\n'
    gmx rdf -n $1_index.ndx -o RDF_O_HW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'O\nHW\n'
    gmx rdf -n $1_index.ndx -o RDF_CH_OW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'CH\nOW\n'
    gmx rdf -n $1_index.ndx -o RDF_CH_HW$1 -f $1.xtc -s $1.tpr -xvg none <<< $'CH\nHW\n'
    gmx rdf -n $1_index.ndx -o RDF_CH_CH$1 -f $1.xtc -s $1.tpr -xvg none <<< $'CH\nCH\n'
}

# run_gromacs_standard {current run} {last run}

xvg_tools.R

#' @export read_xvg
#' @title Read output data from a XVG format file.
#' 
#' @param file A .XVG output file of the GROMACS molecular dynamics package
#' 
#' @description XVG is the default format file of the GROMACS molecular dynamics package, contains data formatted to be imported into the Grace 2-D plotting program.
#' @references  Pronk, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., ... & Lindahl, E. (2013). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, 29 (7), 845-854.
#' @author Daniel Osorio <daniel.osorio@correo.uis.edu.co>
#' @details GROMACS (GROningen MAchine for Chemical Simulations) is a molecular dynamics package designed for simulations of proteins, lipids and nucleic acids. It is free, open source software released under the GNU General Public License. 
#' The file format used by GROMACS is XVG. This format can be displayed in graphical form through the GRACE program on UNIX/LINUX systems and the GNUPlot program on Windows. XVG files are plain text files containing tabular data separated by tabulators and two types of comments which contain data labels. Although manual editing is possible, this is not a viable option when working with multiple files of this type. 
#' For ease of reading, information management and data plotting, the functions \code{read.xvg} and \code{plot.xvg} were incorporated.
#' @examples # READING FILE
#' XVGfile <- system.file("xvg-files/epot.xvg",package="Peptides")
#' read_xvg(XVGfile)
#' 
#' #    Time (ps)  Potential
#' #  1         1 6672471040
#' #  2         2 6516461568
#' #  3         3 6351947264
#' #  4         4 6183133184
#' #  5         5 6015310336
#' #  6         6 5854271488

read.xvg <- function(file) {
  ## Helper functions
  # Remove quotes, starting and ending whitespaces
  unquote <- function(x, ...) { gsub("\\\"|^\\s|\\s$", "", x, ...) }
  
  # Subsets a list for elments that match a regex and then removes that regex
  perlgsub <- function(pattern, x, replacement = "", ...) {
    gsub(pattern = pattern,
         replacement = replacement, 
         x = grep(pattern = pattern, x, value = TRUE, perl = TRUE, ... ),
         perl = TRUE, ...)
  }
  
  # Read flat file
  content <- readLines(file)

  # Read colnames and title
  header <- grep('^[@#]', content, value = TRUE)
  variables <- unquote(perlgsub("^@ s[0-9]+ legend ", header))
  x_axis_label <- unquote(perlgsub("^@ \\s+xaxis\\s+label ", header))
  x_axis_label <- perlgsub("(?!=\\w+)\\W+\\(\\w*\\)$", x_axis_label)
  xvg_labels <- perlgsub("^@ \\s+[a-z]axis\\s+label ", header)
  xvg_labels <- unquote(unlist(strsplit(xvg_labels, ',')))
  title <- unquote(perlgsub("^@ \\s+title ", header))

  # Extracting the data
  content <- perlgsub('^\\s+', content)
  content <- as.data.frame(
    t(sapply(
      content, 
      (function(x) {unlist(strsplit(x, "\\s+"))}), 
      USE.NAMES = FALSE)
      ), stringsAsFactors = FALSE)

  # Asign colnames
  colnames(content) <- c(x_axis_label, variables)

  # Add units and title as attribute
  attr(content, 'xvg_labels') <- xvg_labels
  attr(content, 'title') <- title

  # Return a matrix
  return(content)
}


#' @export ggplot_xvg
#' @title Plot gromacs exported data files
#'
#' @param XVGfile A .XVG output file of the GROMACS molecular dynamics package
#' @inherit read.xvg
#'
#' @description Read and plot output data from a XVG format file.
#' @references  Pronk, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., ... & Lindahl, E. (2013). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, 29 (7), 845-854.
#' @details GROMACS (GROningen MAchine for Chemical Simulations) is a molecular dynamics package designed for simulations of proteins, lipids and nucleic acids. It is free, open source software released under the GNU General Public License. 
#' The file format used by GROMACS is XVG. This format can be displayed in graphical form through the GRACE program on UNIX/LINUX systems and the GNUPlot program on Windows. XVG files are plain text files containing tabular data separated by tabulators and two types of comments which contain data labels. Although manual editing is possible, this is not a viable option when working with multiple files of this type. 
#' For ease of reading, information management and data plotting, the functions \code{read.xvg} and \code{plot.xvg} were incorporated.
#' @examples XVGfile <- system.file("xvg-files/epot.xvg",package="Peptides")
#' ggplot_xvg(XVGfile)

ggplot_xvg <- function(XVGfile, filter = TRUE,...) {
  # Read flat file
  content <- read.xvg(XVGfile)
  xlabel <- colnames(content)[[1]]
  
  try({
    # overwrites the x axis name if it exists
    xlabel <- attr(content, 'xvg_labels')[1]
    })
  
  data <- purrr::map_df(content, as.character)
  data <- reshape2::melt(content, 1)
  
  data[[1]] <- as.numeric(data[[1]])
  data[['value']] <- as.numeric(data[['value']])
  
  aesthetics <- ggplot2::aes_string(
    x = colnames(data)[[1]], 
    y = 'value')
  
  g <- ggplot2::ggplot(
    data, 
    mapping = aesthetics ) 
  g <- g + 
    ggplot2::coord_cartesian(expand = FALSE) +
    ggplot2::facet_wrap(
      ~ variable, scales = "free_y", 
      strip.position = 'left') +
    ggplot2::theme_bw() +
    ggplot2::ggtitle(attr(content, 'title'))
  
  if (filter) {
    filter <- min(max(filter + 1, ceiling(nrow(content)/100)), nrow(content))
    smooth_data <- Reduce(
      rbind,
      lapply(split(data, data$variable), 
             function(x){
               x$value <- zoo::rollmean(
                 x$value,
                 k = filter,
                 align = 'center',
                 na.pad = TRUE )
               return(x)
               }))
    g <- g + ggplot2::geom_line(alpha = 0.35, inherit.aes = TRUE) +
      ggplot2::geom_line(data = smooth_data, mapping = aesthetics)
  } else {
    g <- g + ggplot2::geom_line(inherit.aes = TRUE)
  }
  return(g)
}


#' @export plot.xvg
#' @title Plot time series from GROMACS XVG files
#'
#' @param XVGfile A .XVG output file of the GROMACS molecular dynamics package
#' @param ... Arguments to be passed to methods, such as graphical parameters.
#'
#' @description Read and plot output data from a XVG format file.
#' @references  Pronk, S., Pall, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., ... & Lindahl, E. (2013). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, 29 (7), 845-854.
#' @details GROMACS (GROningen MAchine for Chemical Simulations) is a molecular dynamics package designed for simulations of proteins, lipids and nucleic acids. It is free, open source software released under the GNU General Public License. 
#' The file format used by GROMACS is XVG. This format can be displayed in graphical form through the GRACE program on UNIX/LINUX systems and the GNUPlot program on Windows. XVG files are plain text files containing tabular data separated by tabulators and two types of comments which contain data labels. Although manual editing is possible, this is not a viable option when working with multiple files of this type. 
#' For ease of reading, information management and data plotting, the functions \code{read.xvg} and \code{plot.xvg} were incorporated.
#' @author Daniel Osorio <daniel.osorio@correo.uis.edu.co>
#' @examples XVGfile <- system.file("xvg-files/epot.xvg",package="Peptides")
#' plot.xvg(XVGfile)

plot.xvg <- function(XVGfile, ...) {
  # Read flat file
  content <- read.xvg(XVGfile)
  xlabel <- colnames(content)[[1]]
  
  try({
    # overwrites the x axis name if it exists
    xlabel <- attr(content, 'xvg_labels')[1]
    })
  
  # Get graphical parameters to revert to them later
  original_par <- par(no.readonly = TRUE)
  
  # Plot plot maximum 4 facets per row
  graphics::par(mfcol = c((ncol(content) - 1.5) %/% 4 + 1,
                          min((ncol(content) - 1), 4)), 
                oma = c(0, 0, 2.5, 0))
  for (i in seq_len(ncol(content) - 1)) {
    
    graphics::plot(
      x = content[, 1],
      y = content[, i + 1],
      type = "l",
      ylab = colnames(content)[[i + 1]],
      xlab = xlabel,
      ...
    )
  }
  title(attr(content, 'title'), outer = TRUE)

  # Revert to the previous graphical parameters
  par(original_par)
}