diff --git a/demos/snow_package/mpi_snow_hello.sh b/demos/snow_package/mpi_snow_hello.sh
new file mode 100644
index 0000000..e1305a9
--- /dev/null
+++ b/demos/snow_package/mpi_snow_hello.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+mpiexec -n 1 --hostfile /kube-openmpi/generated/hostfile R --vanilla -f snow-hello.R
diff --git a/demos/snow_package/slurm_submit_snow_hello.sh b/demos/snow_package/slurm_submit_snow_hello.sh
new file mode 100644
index 0000000..7be1804
--- /dev/null
+++ b/demos/snow_package/slurm_submit_snow_hello.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+#SBATCH --ntasks=5
+
+cd $SLURM_SUBMIT_DIR
+
+mpiexec -n 1 R --vanilla -f snow-hello.R
diff --git a/demos/snow_package/snow-hello.R b/demos/snow_package/snow-hello.R
new file mode 100644
index 0000000..fcda6bd
--- /dev/null
+++ b/demos/snow_package/snow-hello.R
@@ -0,0 +1,95 @@
+### Paul Johnson
+### 2017-07-20
+
+### Garrett Mills
+### 2019-05-02
+
+### Demonstration of SNOW "Simple Network of Workstations" using MPI
+### "Message Passing Interface" (OpenMPI implementation)
+
+
+library(snow)
+
+p <- rnorm(123, m = 33)
+
+## Sub script asks for 18 cores, here cluster must
+## be one smaller
+CLUSTERSIZE <- mpi.universe.size()-1
+
+cl <- makeCluster(CLUSTERSIZE, type = "MPI")
+
+
+## One way to send function to each system.
+## Could send identical arguments to nodes
+clusterCall(cl, function() {
+    Sys.info()[c("nodename","machine")]
+    date()
+}
+)
+
+
+timeFn <- function(){
+    y1 <- Sys.info()[c("nodename","machine")]
+    y2 <- date()
+    list(y1, y2)
+}
+
+## Send function to each node
+clusterExport(cl, c("timeFn"))
+
+## Evaluates a string on a node
+x1 <- clusterEvalQ(cl, timeFn())
+
+print(x1)
+
+
+clusterCall(cl, function() rnorm(1, 33, 1) )
+
+
+myNorms <- matrix(rnorm(100*CLUSTERSIZE), ncol = CLUSTERSIZE)
+
+## goes column by column
+mypapply <- parApply(cl, myNorms, 2,  print)
+
+attributes(mypapply)
+
+mypapply <- parApply(cl, myNorms, 2, mean )
+
+mypapply
+
+myNorms <- matrix(rnorm(100*CLUSTERSIZE), ncol = CLUSTERSIZE)
+
+mySum <- function( v ){
+    s <-Sys.info()[c("nodename")]
+    s[2] <- date()
+    ms <- sum(v)
+    list(s, ms)
+}
+
+mypcapply <- parApply(cl, myNorms, 2, mySum)
+
+mypcapply
+
+myNorms <- matrix(rnorm(2500*CLUSTERSIZE), ncol = CLUSTERSIZE)
+
+## Add the system date (includes time) before and after
+## calculations to vector s
+myMeans <- function(v){
+    s <- Sys.info()[c("nodename")]
+    s[2] <- date()
+    ms <- mean(v)
+    ## want to slow this down so you can study the cluster? uncomment this:
+    ## x <- rnorm(50000000); x <- log(50+x); sum(x); mean(x); quantile(x); gg <- cut(x, quantile(x)) 
+    s[3] <- date()
+    list(s, ms)
+}
+
+mypcapply <- parApply(cl, myNorms, 2, myMeans )
+
+mypcapply
+
+
+library(snow)
+stopCluster(cl)
+Rmpi::mpi.quit()
+