MCstuff

find_lj_maximum.R

require(tidyverse)

energies <- (function(){
    tempdata <- data.table::fread('./energy.dat')
    colnames(tempdata) <- c('Step', 'Energy')
    return(tempdata)
})()

max_step <- max(
    energies$Energy[2:length(energies$Energy)] - energies$Energy[1:(length(energies$Energy) - 1)]
                )
print(max_step)

makeinputfiles.R

# 2.5           ! cut-off radius
# 1000          ! Number of particles
# 0.40          ! Density of system (will determine volume)
# 2.00          ! Temperature
# 0             ! start (0 start, 1 resume)
# 0.50          ! maximum displacement (step size)
# 0  100000     ! number of steps to equilibrate system -- number of steps in the production run
# 100           ! frequency of sampling for energy and pressure
# 500           ! frequency of adjustment of max dissplacement
# 

base <- c(cutoff = '2.5', 
          npart = '1000', 
          density = '0.40',
          temperature = '2.00',
          start = '0',
          step_size = '0.5',
          eq_steps_prod_steps = '1000000   100000',
          sample_freq = '100',
          adj_freo = '500')

temps <- c(
  '0.7', '0.9', '1.1', '1.3', '1.5', '1.7', '1.9',                    
  '2.1', '2.3', '2.5', '2.7', '2.9', '3.1', '3.3',
  '3.5', '3.7', '3.9', '4.1', '4.3', '4.5', '4.7',      
  '4.9', '5.1', '5.3', '5.5', '5.7', '5.9')    

densities <-  c(
  '0.1', '0.2', '0.3', 
  '0.4', '0.5', '0.6', 
  '0.7', '0.8')

repetitions <- 1

for (t in temps) {
  for (d in densities) {
    for (r in repetitions) {
      name <- paste0('T_', t, 
                     '_D_', d,
                     '_R_',r,
                     '.input')
      base[['temperature']] <- t
      base[['density']] <- d
      print(base)
      print(name)
      write.table(base, 
                  file = name,
                  row.names = FALSE,
                  col.names = FALSE,
                  quote = FALSE)
    }
  }
}

plot_my_xvg.Rscript

#!/usr/bin/env Rscript
# Usage example
# >for i in $( find | grep -E '.xvg$') ; do Rscript --vanilla plot_my_xvg.Rscript $i ; done
args = commandArgs(trailingOnly = TRUE)

require(Peptides)

data <- read.xvg(args[[1]])
colnames(data) <- make.names(colnames(data))
timecollumn <- colnames(data)[[1]]
data <- reshape2::melt(data, timecollumn)

require(ggplot2)

g <- ggplot(data, aes_string(x = timecollumn, y = 'value')) + 
    geom_line() + 
    coord_cartesian(expand = FALSE) + 
    facet_wrap(~ variable, scales = "free_y") + 
    theme_bw()

ggsave(plot = g, filename = paste0(args[[1]],'.png'),
        width = 10, height = 10, 
        units = 'cm')

read_folder.R

require(tidyverse)

input_names <- c("cutoff", "npart" ,
                 "density" , "temperature" , 
                 "start" , "step_size" ,
                 "eq_steps_prod_steps" , "sample_freq" ,
                 "adj_freo")

parse_log <- function(directory) {
  metadata_file <- paste0(directory, '/log.log')
  metadata <- read_lines(metadata_file) %>%
    grep(':[^$]', .,value = TRUE) %>%
    strsplit(':')
  metadata_names <- purrr::map_chr(metadata,
                                   function(X){
                                     make.names(X[[1]])
                                     }) %>%
    gsub('\\.+', '.', . ) %>%
    gsub('^X\\.|\\.$', '', . ) %>%
    gsub('\\.', '_', . ) 
  metadata_values <- purrr::map_dbl(metadata, function(X){as.numeric(X[[2]])})
  names(metadata_values) <- metadata_names
  return(metadata_values)
}

read_montecarlo_folder <- function(directory) {
  ret_list <- list()
  input_file <- paste0(directory,'/input')
  energy_file <- paste0(directory, '/energy.dat')
  pressure_file <- paste0(directory, '/pressure.dat')
  
  print(input_file)
  ret_list[['input']] <- unlist(
    read.table(input_file,
               as.is = TRUE, 
               colClasses = 'character',
               sep = '\n'),
    use.names = FALSE)
  names(ret_list[['input']]) <- input_names
  pressures <- data.table::fread(input = pressure_file)
  colnames(pressures) <- c('Step', 'Pressure')
  energies  <- data.table::fread(input = energy_file)
  colnames(energies) <- c('Step', 'Energy')
  ret_list[['df']] <- dplyr::full_join(energies, pressures)
  ret_list[['df']][['Temperature']] <- ret_list[['input']][['temperature']]
  ret_list[['df']][['Density']] <- ret_list[['input']][['density']]
  ret_list[['df']][['Iteration']] <- make.names(gsub('.*_R_', "", directory))
  return(ret_list[['df']])
}


dirs <- dir('./multi_inputs/',
            pattern = 'T_',
            full.names = TRUE, 
            include.dirs = TRUE)

my_data <- purrr::map(dirs, read_montecarlo_folder)
my_data <- purrr::reduce(my_data, dplyr::full_join)

my_logs <- purrr::map(dirs, parse_log) %>% purrr::reduce(rbind) %>% tbl_df()


full_data <- my_data %>% 
  group_by(Temperature, Density, Iteration) %>%
  summarise(meanEnergy = mean(Energy),
            meanPressure = mean(Pressure),
            Pressure_sd = sd(Pressure),
            Energy_sd = sd(Energy)) %>%
  ungroup() %>%
  full_join(my_logs %>% mutate(Density = as.character(Density),
                               Temperature = as.character(Temperature))) %>%
  mutate(Pressure = meanPressure + Pressure_correction,
         Energy = meanEnergy + Energy_correction)

g <- ggplot(full_data, 
            aes(x = Density, 
                y = Pressure, 
                colour = Temperature,
                group = Temperature)) + 
  geom_line(alpha = 0.6) + theme_bw() + 
  #geom_errorbar(aes(ymin = meanPressure - Pressure_sd, ymax = meanPressure + Pressure_sd)) +
  coord_cartesian(expand = FALSE)
g

g <- ggplot(full_data, 
            aes(x = Density, 
                y = Pressure, 
                colour = Temperature,
                group = Temperature)) + 
  geom_line(alpha = 0.6) + theme_bw() + 
  #geom_errorbar(aes(ymin = meanPressure - Pressure_sd, ymax = meanPressure + Pressure_sd)) +
  coord_cartesian(expand = FALSE)
g