newman_mod <- function(g, weights=NULL){
  A <- get.adjacency(g) # adj
  m <- ecount(g)
  n <- vcount(g)
  
  if (is.null(weights)){
    weights <- rep(1,n)
  }
  
  # Obtain the modularity matrix
  B.node.i <- function(i){degree(g)[i]*degree(g)}
  B.node.all <- sapply(1:n, B.node.i)
  B <- A - (B.node.all/(2*m))
  # NOTE: This is identical to: modularity_matrix(g) ! Can verify with:
  # modularity_matrix(g) == B
  B.eigs <- eigen(B)
  max.lam <- B.eigs$values[1]
  s <- ifelse(B.eigs$vectors[,1]>0,1,-1)
  
  weights <- B.eigs$vectors[n]/B.eigs$vectors[,1]
  # Plotting
  #V(g)$color <- ifelse(B[1,]>0,"green","yellow")
  V(g)$color <- ifelse(B.eigs$vectors[,1]>0,"green","yellow")
  V(g)$size <- 10
  plot(g, main=paste(g1.netname, " Newman Modularity"))
  clust1 = list()
  clust2 = list()
  clusters = list()
  
  # Make list of clusters
  for(i in 1:n){
      ifelse(V(g)[i]$color=="green",
              clust1 <- append(clust1, V(g)[i]$name),
              clust2 <- append(clust2, V(g)[i]$name))}

  clusters <- list(clust1, clust2)
  
  Q.node.i <- function(i){sum(
    (((B.eigs$vectors[i])*weights[i]*s)^2)*B.eigs[i]$values)}
  Q <- (1/(4*m))*sapply(1:n, Q.node.i)
 
  return(list(Q=Q,max.lam=max.lam,weights=weights,clusters=clusters))
}