Nothing
R/readingandfes.R
In metadynminer: Tools to Read, Analyze and Visualize Metadynamics HILLS Files from 'Plumed'
Defines functions
lines.fes
points.fes
plot.fes
prob
summary.fes
print.fes
max.fes
min.fes
`/.fes`
`*.fes`
`-.fes`
`+.fes`
fes2d21d
read.plumed
fes2.hillsfile
fes2
fes.hillsfile
fes
plotheights.hillsfile
plotheights
lines.hillsfile
points.hillsfile
plot.hillsfile
`+.hillsfile`
tail.hillsfile
head.hillsfile
summary.hillsfile
print.hillsfile
emptyhills
read.hills
Documented in
emptyhills
fes
fes2
fes2d21d
fes2.hillsfile
fes.hillsfile
head.hillsfile
lines.fes
lines.hillsfile
max.fes
min.fes
plot.fes
plotheights
plotheights.hillsfile
plot.hillsfile
points.fes
points.hillsfile
print.fes
print.hillsfile
prob
read.hills
read.plumed
summary.fes
summary.hillsfile
tail.hillsfile
#' @useDynLib metadynminer
#' @importFrom Rcpp sourceCpp
NULL
#' Read HILLS from Plumed
#'
#' `read.hills` reads a HILLS file generated by Plumed and returns a hillsfile object.
#' User can specify whether some collective variables are periodic.
#'
#' @param file HILLS file from Plumed.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.409 2.808 0.3 0.3 1.111 10"
#' l2<-"2 -2.505 2.791 0.3 0.3 1.111 10"
#' l3<-"3 -2.346 2.754 0.3 0.3 1.069 10"
#' l4<-"4 -1.198 2.872 0.3 0.3 1.074 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills(tf, per=c(TRUE,TRUE))
read.hills<-function(file="HILLS", per=c(FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), ignoretime=FALSE) {
hillsf<-read.table(file, header=F, comment.char="#")
if(ncol(hillsf)==5 || ncol(hillsf)==6) {
cat("1D HILLS file read\n")
if(ignoretime) {
cat("Warning: The time will be updated automatically from zero\n")
cat("according to the first step!\n")
hillsf[,1]<-seq(from=hillsf[1,1], by=hillsf[1,1], length.out=nrow(hillsf))
}
hills<-list(hillsfile=hillsf, time=hillsf[,1], cv1=hillsf[,2], cv2=NULL,
size=dim(hillsf), filename=file, per=per, pcv1=pcv1)
class(hills) <- "hillsfile"
return(hills)
} else if(ncol(hillsf)==7 || ncol(hillsf)==8) {
cat("2D HILLS file read\n")
if(ignoretime) {
cat("Warning: The time will be updated automatically from zero\n")
cat("according to the first step!\n")
hillsf[,1]<-seq(from=hillsf[1,1], by=hillsf[1,1], length.out=nrow(hillsf))
}
hills<-list(hillsfile=hillsf, time=hillsf[,1], cv1=hillsf[,2], cv2=hillsf[,3],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2)
class(hills) <- "hillsfile"
return(hills)
} else {
stop("Error: Number of columns in HILLS file must be 5 or 6 (1D) or 7 or 8 (2D)")
}
}
#' Generate empty HILLS
#'
#' `emptyhills` returns a hillsfile object with no hills.
#' User can specify whether some collective variables are periodic.
#'
#' @param dim dimensionality of the output.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' emptyhills(dim=2)"
emptyhills<-function(dim=2, per=c(FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi)) {
if(dim==1) {
cat("empty 1D hillsfile generated\n")
hills<-list(hillsfile=NULL, time=NULL, cv1=NULL, cv2=NULL,
size=c(0,5), filename=NULL, per=per, pcv1=pcv1)
class(hills) <- "hillsfile"
return(hills)
} else if(dim) {
cat("empty 2D hillsfile generated\n")
hills<-list(hillsfile=NULL, time=NULL, cv1=NULL, cv2=NULL,
size=c(0,7), filename=NULL, per=per, pcv1=pcv1, pcv2=pcv2)
class(hills) <- "hillsfile"
return(hills)
} else {
stop("Error: works with dim=2 or dim=3")
}
}
#' Print hillsfile
#'
#' `print.hillsfile` prints dimensionality and size of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' acealanme
print.hillsfile<-function(x,...) {
hills <- x
if(hills$size[2]==5) {
cat("1D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
}
if(hills$size[2]==7) {
cat("2D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
}
}
#' Print summary for hillsfile
#'
#' `summary.hillsfile` prints dimensionality, size and collective variable ranges of a hillsfile object.
#'
#' @param object hillsfile object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' summary(acealanme)
summary.hillsfile<-function(object,...) {
hills <- object
if(hills$size[2]==5) {
cat("1D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
cat("The CV1 ranges from ")
cat(min(hills$hillsfile[,2]))
cat(" to ")
cat(max(hills$hillsfile[,2]))
cat("\n")
}
if(hills$size[2]==7) {
cat("2D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
cat("The CV1 ranges from ")
cat(min(hills$hillsfile[,2]))
cat(" to ")
cat(max(hills$hillsfile[,2]))
cat("\nThe CV2 ranges from ")
cat(min(hills$hillsfile[,3]))
cat(" to ")
cat(max(hills$hillsfile[,3]))
cat("\n")
}
}
#' Print first n lines of hillsfile
#'
#' `head.hillsfile` prints first n lines of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' head(acealanme)
head.hillsfile<-function(x, n=10,...) {
return(head(x$hillsfile, n=n))
}
#' Print last n lines of hillsfile
#'
#' `tail.hillsfile` prints last n lines of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tail(acealanme)
tail.hillsfile<-function(x, n=10,...) {
return(tail(x$hillsfile, n=n))
}
#' @export
`+.hillsfile`<-function(hills1, hills2) {
if(ncol(hills1$hillsfile)!=ncol(hills2$hillsfile)) {
stop("Error: You can sum only hills of same dimension")
}
if(hills1$per[1]!=hills2$per[1]) {
stop("Error: You can sum only hills of same periodicity")
}
if(ncol(hills1$hillsfile)==7 || ncol(hills1$hillsfile)==8) {
if(hills1$per[2]!=hills2$per[2]) {
stop("Error: You can sum only hills of same periodicity")
}
}
hills<-list(hillsfile=rbind(hills1$hillsfile, hills2$hillsfile), size=dim(rbind(hills1$hillsfile, hills2$hillsfile)),
filename=hills1$filename, per=hills1$per, pcv1=hills1$pcv1, pcv2=hills1$pcv2)
class(hills) <- "hillsfile"
return(hills)
}
#' Plot hillsfile object
#'
#' `plot.hillsfile` plots hillsfile object. For a hillsfile with one collective variable it plots its evolution.
#' For a hillsfile with two collective variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'.
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
plot.hillsfile<-function(x, ignoretime=FALSE,
xlab=NULL, ylab=NULL,
xlim=NULL, ylim=NULL,
main=NULL, sub=NULL,
pch=1, col="black", bg="red", cex=1,
asp=NULL, lwd=1, axes=TRUE,...) {
hills <-x
xlims<-NULL
ylims<-NULL
if(!is.null(xlim)) {xlims<-xlim}
if(!is.null(ylim)) {ylims<-ylim}
if(hills$size[2]==5) {
if((hills$per[1]==T)&is.null(ylim)) {ylims<-hills$pcv1}
if(is.null(xlab)) xlab="time"
if(is.null(ylab)) ylab="CV"
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
col=col, cex=cex, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,2], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
col=col, cex=cex, lwd=lwd,
asp=asp, axes=axes)
}
}
if(hills$size[2]==7) {
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
plot(hills$hillsfile[,2], hills$hillsfile[,3], type="p",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd,
asp=asp, axes=axes)
}
}
#' Plot points for hillsfile object
#'
#' `points.hillsfile` plots points for hillsfile object. For a hillsfile with one
#' collective variable it plots its evolution. For a hillsfile with two collective
#' variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'.
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
#' points(acealanme, col="red")
points.hillsfile<-function(x, ignoretime=FALSE,
pch=1, col="black", bg="red", cex=1,
lwd=1, ...) {
hills <- x
if(hills$size[2]==5) {
if(ignoretime) {
points(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2],
col=col, cex=cex, lwd=lwd)
} else {
points(hills$hillsfile[,1], hills$hillsfile[,2],
col=col, cex=cex, lwd=lwd)
}
}
if(hills$size[2]==7) {
points(hills$hillsfile[,2], hills$hillsfile[,3],
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd)
}
}
#' Plot lines for hillsfile object
#'
#' `lines.hillsfile` plots lines for hillsfile object. For a hillsfile with one
#' collective variable it plots its evolution. For a hillsfile with two collective
#' variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
#' lines(acealanme, col="red")
lines.hillsfile<-function(x, ignoretime=FALSE,
lwd=1, col="black",...) {
hills <- x
if(hills$size[2]==5) {
if(ignoretime) {
lines(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2],
col=col, lwd=lwd)
} else {
lines(hills$hillsfile[,1], hills$hillsfile[,2],
col=col, lwd=lwd)
}
}
if(hills$size[2]==7) {
lines(hills$hillsfile[,2], hills$hillsfile[,3],
col=col, lwd=lwd)
}
}
#' Plot evolution of heights of hills (generic function for 'metadynminer' and
#' 'metadynminer3d')
#'
#' `plotheights` plots evolution of heights of hills. In well tempered metadynamics
#' hill heights decrees with flooding of the free energy surface. Evolution of heights
#' may be useful to evaluate convergence of the simulation.
#'
#' @param hills hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#'
#' @export plotheights
plotheights<-function(hills, ignoretime, xlab, ylab,
xlim, ylim, main, sub,
col, asp, lwd, axes) {
UseMethod("plotheights")
}
#' Plot evolution of heights of hills in hillsfile object
#'
#' `plotheights.hillsfile` plots evolution of heights of hills. In well tempered metadynamics
#' hill heights decrees with flooding of the free energy surface. Evolution of heights
#' may be useful to evaluate convergence of the simulation.
#'
#' @param hills hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#'
#' @export
#' @examples
#' plotheights(acealanme)
plotheights.hillsfile<-function(hills, ignoretime=FALSE,
xlab=NULL, ylab=NULL,
xlim=NULL, ylim=NULL,
main=NULL, sub=NULL,
col="black", asp=NULL, lwd=1, axes=TRUE) {
if(is.null(xlab)) xlab="time"
if(is.null(ylab)) ylab="hill height"
if(hills$size[2]==5) {
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,4], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,4], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
if(hills$size[2]==7) {
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,6], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,6], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
}
#' Calculate free energy surface by Bias Sum algorithm (generic function for
#' 'metadynminer' and 'metadynminer3d')
#'
#' `fes` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export fes
fes<-function(hills, imin, imax, xlim, ylim, zlim, npoints) {
UseMethod("fes")
}
#' Calculate free energy surface by Bias Sum algorithm
#'
#' `fes.hillsfile` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
fes.hillsfile<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imax>0 & imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
if(is.null(ylim)) {
minCV2 <- 0
maxCV2 <- 1
} else {
minCV2 <- ylim[1]
maxCV2 <- ylim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
fesm <- matrix(rep(0, npoints*npoints), nrow=npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
if(max(hills$hillsfile[,4])/min(hills$hillsfile[,4])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
if(hills$per[1]==T) {
if(min(hills$hillsfile[,2])<xlims[1]) {
stop("Error: The first collective variable outside pcv1")
}
if(max(hills$hillsfile[,2])>xlims[2]) {
stop("Error: The first collective variable outside pcv1")
}
}
if(hills$per[2]==T) {
if(min(hills$hillsfile[,3])<ylims[1]) {
stop("Error: The second collective variable outside pcv2")
}
if(max(hills$hillsfile[,3])>ylims[2]) {
stop("Error: The second collective variable outside pcv2")
}
}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
if(hills$size[2]==5) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
fesm <- rep(0, npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
if(max(hills$hillsfile[,3])/min(hills$hillsfile[,3])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
if(hills$per[1]==F) {
fesm<-hills1d1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,3])/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
if(hills$per[1]==T) {
fesm<-hills1d1p(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,3])/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
return(cfes)
}
#' Calculate free energy surface by conventional algorithm (generic function
#' for 'metadynminer' and 'metadynminer3d')
#'
#' `fes2` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export fes2
fes2<-function(hills, imin, imax, xlim, ylim, zlim, npoints) {
UseMethod("fes2")
}
#' Calculate free energy surface by conventional algorithm
#'
#' `fes2.hillsfile` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes2(acealanme, imax=1000)
fes2.hillsfile<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imax>0 & imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
if(is.null(ylim)) {
minCV2 <- 0
maxCV2 <- 1
} else {
minCV2 <- ylim[1]
maxCV2 <- ylim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
fesm <- matrix(rep(0, npoints*npoints), nrow=npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills2p1((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills2p2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills2p12((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
if(hills$size[2]==5) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
fesm <- rep(0, npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
if(hills$per[1]==F) {
fesm<-hills1d2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,3]/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
if(hills$per[1]==T) {
fesm<-hills1d2p((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,3]/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
return(cfes)
}
#' Read 1D or 2D free energy surface from PLUMED sum_hills
#'
#' `read.plumed` reads 1D or 2D free energy surface from PLUMED sum_hills.
#' The grid in the (2D) inputfile must contain the same number of points
#' for CV1 and CV2. It does not use the header of the file. Instead, user
#' must specify the dimensionality (1 or 2). Periodicity must be specified
#' as well.
#'
#' @param file input file from PLUMED sum_hills.
#' @param dim dimension (1 or 2, default 2).
#' @param per logical vector specifying periodicity of collective variables.
#' @return fes object.
#'
#' @export
#' @examples
#' l1<-"-3.142 -124.8 -44.76"
#' l2<-"-3.117 -125.9 -43.05"
#' l3<-"-3.092 -126.9 -41.22"
#' l4<-"-3.068 -127.9 -39.36"
#' l5<-"-3.043 -128.8 -37.45"
#' fourpoints<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourpoints, tf)
#' read.plumed(tf, dim=1, per=c(TRUE,TRUE))
read.plumed<-function(file="fes.dat", dim=2, per=c(F,F,F)) {
hillsf<-read.table(file, header=F, comment.char="#")
bins<-round(nrow(hillsf)^(1/dim))
if(bins^dim!=nrow(hillsf)) {
stop("Error: the number of bins cannot be determined, it must be same for all dimension, or the number of dimensions is wrong.")
}
if(dim==1) {
x <- hillsf[,1]
fesm <- hillsf[,2]
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=1, per=per, x=x, pcv1=c(min(x), max(x)))
class(cfes) <- "fes"
} else if(dim==2) {
x <- hillsf[1:bins,1]
y <- hillsf[(0:(bins-1))*bins+1,2]
fesm <- matrix(hillsf[,3], nrow=bins)
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=2, per=per, x=x, y=y,
pcv1=c(min(x), max(x)), pcv2=c(min(y), max(y)))
class(cfes) <- "fes"
} else {
stop("Error: for 3D fes use read.plumed3d from metadynminer3d, higher dimensions are not supported.")
}
return(cfes)
}
#' Calculate 1D free energy surface from hillsfile object
#'
#' `fes2d21d` calculates 2D free energy surface, converts free energies to probabilities
#' (exp(-F/kT)), sums them up along one collective variable and converts back to free
#' energy (-kT log(P)).
#'
#' @param hills hillsfile object.
#' @param remdim dimension to be removed (1 for CV1, 2 for CV2, default 2).
#' @param temp temperature in Kelvins (default 300).
#' @param eunit energy units (kJ/mol or kcal/mol, kJ/mol is default).
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes2d21d(acealanme, remdim=2, imax=5000)
fes2d21d<-function(hills, remdim=2, temp=300, eunit="kJ/mol",
imin=1, imax=NULL, xlim=NULL, ylim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(max(hills$hillsfile[,4])/min(hills$hillsfile[,4])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
binx<-(xlims[2]-xlims[1])/(npoints-1)
biny<-(ylims[2]-ylims[1])/(npoints-1)
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if(eunit=="kJ/mol") {
prob<- exp(-1000*fesm/8.314/temp)
if(remdim==1) {
fesm <- -8.314*temp*log(apply(prob, 2, sum))/1000
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[2], x=y, pcv1=hills$pcv2)
}
if(remdim==2) {
fesm <- -8.314*temp*log(apply(prob, 1, sum))/1000
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[1], x=x, pcv1=hills$pcv1)
}
}
if(eunit=="kcal/mol") {
prob<- exp(-1000*4.184*fesm/8.314/temp)
if(remdim==1) {
fesm <- -8.314*temp*log(apply(prob, 2, sum))/1000/4.184
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[2], x=y, pcv1=hills$pcv2)
}
if(remdim==2) {
fesm <- -8.314*temp*log(apply(prob, 1, sum))/1000/4.184
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[1], x=x, pcv1=hills$pcv1)
}
}
class(cfes) <- "fes"
}
if(hills$size[2]==5) {
stop("Error: Your free energy surface is already 1D")
}
return(cfes)
}
#' @export
`+.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(fes1$dimension!=fes2$dimension) {
stop("Error: Free energy surfaces have different dimension, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(fes1$dimension==2) {
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
}
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes+fes2$fes, hills=rbind(fes1$hills, fes2$hills), rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes+fes2$fes, hills=rbind(fes1$hills, fes2$hills), rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1+fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1+fes2$fes, hills=rbind(fes1$hills,fes2$hills), rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
class(cfes) <- "fes"
return(cfes)
}
#' @export
`-.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(fes1$dimension!=fes2$dimension) {
stop("Error: Free energy surfaces have different dimension, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(fes1$dimension==2) {
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
}
cat("Warning: FES obtained by subtraction of two FESes\n")
cat(" will inherit hills only from the first FES\n")
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes-fes2$fes, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes-fes2$fes, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
class(cfes) <- "fes"
return(cfes)
}
#' @export
`*.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
stop("Error: You cannot multiply fes by fes")
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
cat("Warning: multiplication of FES will multiply\n")
cat(" the FES but not hill heights\n")
class(cfes) <- "fes"
return(cfes)
}
#' @export
`/.fes`<-function(fes1, coef) {
if(inherits(fes1, "fes")&inherits(coef, "fes")) {
stop("Error: You cannot divide fes by fes")
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(coef, "fes")) {
stop("Error: You cannot divide something by fes")
}
cat("Warning: division of FES will divide\n")
cat(" the FES but not hill heights\n")
class(cfes) <- "fes"
return(cfes)
}
#' Calculate minimum of free energy surface
#'
#' `min.fes` calculates minimum of free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param na.rm a logical indicating whether missing values should be
#' removed.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' min(tfes)
min.fes<-function(inputfes, na.rm=NULL,...) {
return(min(inputfes$fes, na.rm=na.rm))
}
#' Calculate maximum of free energy surface
#'
#' `max.fes` calculates maximum of free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param na.rm a logical indicating whether missing values should be
#' removed.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' max(tfes)
max.fes<-function(inputfes, na.rm=NULL,...) {
return(max(inputfes$fes, na.rm=na.rm))
}
#' Print dimensionality, minimum and maximum of free energy surface
#'
#' `print.fes` prints dimensionality, minimum and maximum of
#' free energy in a fes object
#'
#' @param x fes object
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' tfes
print.fes<-function(x,...) {
inputfes<-x
if(inputfes$dimension==1) {
cat("1D free energy surface with ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
if(inputfes$dimension==2) {
cat("2D free energy surface with ")
cat(inputfes$rows)
cat(" x ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
}
#' Print summary of free energy surface
#'
#' `summary.fes` prints dimensionality, minimum and maximum of
#' free energy in a fes object.
#'
#' @param object fes object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' summary(tfes)
summary.fes<-function(object,...) {
inputfes <- object
if(inputfes$dimension==1) {
cat("1D free energy surface with ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
if(inputfes$dimension==2) {
cat("2D free energy surface with ")
cat(inputfes$rows)
cat(" x ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
}
#' Calculate probability of free energy surface
#'
#' `prob` calculates probability from free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param temp temperature in Kelvins (default 300).
#' @param eunit energy units (kJ/mol or kcal/mol, kJ/mol is default).
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' print(prob(tfes))
prob<-function(inputfes, temp=300, eunit="kJ/mol") {
if(inherits(inputfes, "fes")) {
if(eunit=="kJ/mol") {
if(inputfes$dimension==1) {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, pcv1=inputfes$pcv1)
}
if(inputfes$dimension==2) {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2)
}
} else if (eunit=="kJ/mol") {
if(inputfes$dimension==1) {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, pcv1=inputfes$pcv1)
}
if(inputfes$dimension==2) {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2)
}
} else {
stop("Error: Wrong eunit")
}
class(cfes) <- "fes"
return(cfes)
} else if(inherits(inputfes, "fes3d")) {
if(eunit=="kJ/mol") {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, z=inputfes$z, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2, pcv3=inputfes$pcv3)
} else if (eunit=="kJ/mol") {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, z=inputfes$z, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2, pcv3=inputfes$pcv3)
} else {
stop("Error: Wrong eunit")
}
class(cfes) <- "fes3d"
return(cfes)
} else {
stop("Error: Input must be fes or fes3d object")
}
}
#' Plot free energy surface object
#'
#' `plot.fes` plots free energy surface. For a fes with one collective variable it plots
#' a 1D profile. For a fes with two collective variables it plots 2D free energy surface
#' using image, contours or combination of both (default).
#'
#' @param x fes object.
#' @param plottype specifies whether 2D free energy surface will be plotted
#' as image, contours or both (default "both").
#' @param colscale specifies whether color scale will be plotted (default False).
#' @param colscalelab color scale label (default "free energy").
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color of the free energy surface. For 1D surface it is the color
#' of the line. For 2D it is a list of colors such as that generated by
#' 'rainbow', 'heat.colors', 'topo.colors', 'terrain.colors' or similar
#' functions (default=rainbow(135)[100:1]).
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param zlim numeric vector of length 2, giving the z coordinates range.
#' @param nlevels number of contour levels desired if 'levels' is not
#' supplied.
#' @param levels numeric vector of levels at which to draw contour lines.
#' @param labels a vector giving the labels for the contour lines. If 'NULL'
#' then the levels are used as labels, otherwise this is coerced
#' by 'as.character'.
#' @param labcex 'cex' for contour labeling. This is an absolute size, not a
#' multiple of 'par("cex")'.
#' @param drawlabels logical. Contours are labeled if 'TRUE'.
#' @param method character string specifying where the labels will be located.
#' Possible values are '"simple"', '"edge"' and '"flattest"'
#' (the default). See the 'Details' section.
#' @param lwd contour line width.
#' @param contcol contour color.
#' @param lty line type for the lines drawn.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes2d<-fes(acealanme, imax=5000)
#' plot(tfes2d)
#' tfes1d<-fes(acealanme1d)
#' plot(tfes1d)
plot.fes<-function(x, plottype="both",
colscale=F, xlim=NULL, ylim=NULL, zlim=NULL,
main=NULL, sub=NULL,
xlab=NULL, ylab=NULL,
nlevels=10, levels=NULL,
col=rainbow(135)[100:1],
labels=NULL, labcex=0.6, drawlabels=TRUE,
colscalelab="free energy",
method="flattest",
contcol=par("fg"), lty=par("lty"),
lwd=1, asp=NULL, axes=T,...) {
close.screen(all.screens=TRUE)
inputfes<-x
fes<-inputfes$fes
rows<-inputfes$rows
if(inputfes$dimension==1) {
x<-inputfes$x
if(is.null(xlab)) xlab="CV"
if(is.null(ylab)) ylab="free energy"
if(is.null(xlim)) xlim<-c(min(x),max(x))
if(is.null(ylim)) {
ylim<-range(pretty(range(fes)))
}
plot(x, fes, type="l", lwd=lwd,
col=col, xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, axes=axes,
main=main, sub=sub, asp=asp)
} else {
x<-inputfes$x
y<-inputfes$y
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
if(is.null(zlim)) {
zlim<-range(pretty(range(fes)))
}
if(is.null(levels)) {
levels<-pretty(zlim, nlevels)
}
if(is.null(xlim)) xlim<-c(min(x),max(x))
if(is.null(ylim)) ylim<-c(min(y),max(y))
if(colscale) {
#layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(4,1))
split.screen(matrix(c(0,0.75,0,1,0.75,1,0,1), byrow=T, ncol=4))
screen(2)
smat<-matrix(seq(from=zlim[1], to=zlim[2], length.out=100))
image(c(0), seq(from=zlim[1], to=zlim[2], length.out=100),
t(smat), zlim=zlim, col=col, xlab="", ylab=colscalelab, axes=F)
axis(2, lty=lty, lwd=lwd)
box(lwd=lwd)
screen(1)
}
if(plottype=="image" || plottype=="both") {
image(x, y, fes, zlim=zlim,
col=col, xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, axes=axes,
main=main, sub=sub, asp=asp)
}
if(plottype=="contour") {
contour(x, y, fes, zlim=zlim,
nlevels=nlevels, levels=levels,
labels=labels, labcex=labcex, drawlabels=drawlabels,
method=method, col=contcol, lty=lty, lwd=lwd,
main=main, sub=sub, asp=asp)
}
if(plottype=="both") {
contour(x, y, fes, zlim=zlim,
nlevels=nlevels, levels=levels,
labels=labels, labcex=labcex, drawlabels=drawlabels,
method=method, col=contcol, lty=lty, lwd=lwd, add=T)
}
}
}
#' Plots 1D free energy surface object as points
#'
#' `points.fes` plots 1D free energy surface as points.
#'
#' @param x fes object.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme1d, imax=5000)
#' plot(tfes)
#' points(tfes)
points.fes<-function(x, pch=1, col="black", bg="red", cex=1, lwd=1,...) {
fes<-x$fes
if(x$dimension==1) {
x<-x$x
points(x, fes,
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd)
} else {
stop("Error: points available only for 1D free energy surfaces\n")
}
}
#' Plots 1D free energy surface object as lines
#'
#' `lines.fes` plots 1D free energy surface as lines.
#'
#' @param x fes object.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme1d, imax=5000)
#' plot(tfes)
#' lines(tfes, lwd=4)
lines.fes<-function(x, lwd=1, col="black",...) {
fes<-x$fes
if(x$dimension==1) {
x<-x$x
lines(x, fes, lwd=lwd, col=col)
} else {
stop("Error: points available only for 1D free energy surfaces\n")
}
}
Try the metadynminer package in your browser
Any scripts or data that you put into this service are public.
metadynminer documentation built on April 14, 2022, 5:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lines.fes points.fes plot.fes prob summary.fes print.fes max.fes min.fes `/.fes` `*.fes` `-.fes` `+.fes` fes2d21d read.plumed fes2.hillsfile fes2 fes.hillsfile fes plotheights.hillsfile plotheights lines.hillsfile points.hillsfile plot.hillsfile `+.hillsfile` tail.hillsfile head.hillsfile summary.hillsfile print.hillsfile emptyhills read.hills
Documented in emptyhills fes fes2 fes2d21d fes2.hillsfile fes.hillsfile head.hillsfile lines.fes lines.hillsfile max.fes min.fes plot.fes plotheights plotheights.hillsfile plot.hillsfile points.fes points.hillsfile print.fes print.hillsfile prob read.hills read.plumed summary.fes summary.hillsfile tail.hillsfile
#' @useDynLib metadynminer
#' @importFrom Rcpp sourceCpp
NULL
#' Read HILLS from Plumed
#'
#' `read.hills` reads a HILLS file generated by Plumed and returns a hillsfile object.
#' User can specify whether some collective variables are periodic.
#'
#' @param file HILLS file from Plumed.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.409 2.808 0.3 0.3 1.111 10"
#' l2<-"2 -2.505 2.791 0.3 0.3 1.111 10"
#' l3<-"3 -2.346 2.754 0.3 0.3 1.069 10"
#' l4<-"4 -1.198 2.872 0.3 0.3 1.074 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills(tf, per=c(TRUE,TRUE))
read.hills<-function(file="HILLS", per=c(FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), ignoretime=FALSE) {
hillsf<-read.table(file, header=F, comment.char="#")
if(ncol(hillsf)==5 || ncol(hillsf)==6) {
cat("1D HILLS file read\n")
if(ignoretime) {
cat("Warning: The time will be updated automatically from zero\n")
cat("according to the first step!\n")
hillsf[,1]<-seq(from=hillsf[1,1], by=hillsf[1,1], length.out=nrow(hillsf))
}
hills<-list(hillsfile=hillsf, time=hillsf[,1], cv1=hillsf[,2], cv2=NULL,
size=dim(hillsf), filename=file, per=per, pcv1=pcv1)
class(hills) <- "hillsfile"
return(hills)
} else if(ncol(hillsf)==7 || ncol(hillsf)==8) {
cat("2D HILLS file read\n")
if(ignoretime) {
cat("Warning: The time will be updated automatically from zero\n")
cat("according to the first step!\n")
hillsf[,1]<-seq(from=hillsf[1,1], by=hillsf[1,1], length.out=nrow(hillsf))
}
hills<-list(hillsfile=hillsf, time=hillsf[,1], cv1=hillsf[,2], cv2=hillsf[,3],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2)
class(hills) <- "hillsfile"
return(hills)
} else {
stop("Error: Number of columns in HILLS file must be 5 or 6 (1D) or 7 or 8 (2D)")
}
}
#' Generate empty HILLS
#'
#' `emptyhills` returns a hillsfile object with no hills.
#' User can specify whether some collective variables are periodic.
#'
#' @param dim dimensionality of the output.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' emptyhills(dim=2)"
emptyhills<-function(dim=2, per=c(FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi)) {
if(dim==1) {
cat("empty 1D hillsfile generated\n")
hills<-list(hillsfile=NULL, time=NULL, cv1=NULL, cv2=NULL,
size=c(0,5), filename=NULL, per=per, pcv1=pcv1)
class(hills) <- "hillsfile"
return(hills)
} else if(dim) {
cat("empty 2D hillsfile generated\n")
hills<-list(hillsfile=NULL, time=NULL, cv1=NULL, cv2=NULL,
size=c(0,7), filename=NULL, per=per, pcv1=pcv1, pcv2=pcv2)
class(hills) <- "hillsfile"
return(hills)
} else {
stop("Error: works with dim=2 or dim=3")
}
}
#' Print hillsfile
#'
#' `print.hillsfile` prints dimensionality and size of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' acealanme
print.hillsfile<-function(x,...) {
hills <- x
if(hills$size[2]==5) {
cat("1D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
}
if(hills$size[2]==7) {
cat("2D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
}
}
#' Print summary for hillsfile
#'
#' `summary.hillsfile` prints dimensionality, size and collective variable ranges of a hillsfile object.
#'
#' @param object hillsfile object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' summary(acealanme)
summary.hillsfile<-function(object,...) {
hills <- object
if(hills$size[2]==5) {
cat("1D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
cat("The CV1 ranges from ")
cat(min(hills$hillsfile[,2]))
cat(" to ")
cat(max(hills$hillsfile[,2]))
cat("\n")
}
if(hills$size[2]==7) {
cat("2D hills file ")
cat(hills$filename)
cat(" with ")
cat(hills$size[1])
cat(" lines\n")
cat("The CV1 ranges from ")
cat(min(hills$hillsfile[,2]))
cat(" to ")
cat(max(hills$hillsfile[,2]))
cat("\nThe CV2 ranges from ")
cat(min(hills$hillsfile[,3]))
cat(" to ")
cat(max(hills$hillsfile[,3]))
cat("\n")
}
}
#' Print first n lines of hillsfile
#'
#' `head.hillsfile` prints first n lines of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' head(acealanme)
head.hillsfile<-function(x, n=10,...) {
return(head(x$hillsfile, n=n))
}
#' Print last n lines of hillsfile
#'
#' `tail.hillsfile` prints last n lines of a hillsfile object.
#'
#' @param x hillsfile object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tail(acealanme)
tail.hillsfile<-function(x, n=10,...) {
return(tail(x$hillsfile, n=n))
}
#' @export
`+.hillsfile`<-function(hills1, hills2) {
if(ncol(hills1$hillsfile)!=ncol(hills2$hillsfile)) {
stop("Error: You can sum only hills of same dimension")
}
if(hills1$per[1]!=hills2$per[1]) {
stop("Error: You can sum only hills of same periodicity")
}
if(ncol(hills1$hillsfile)==7 || ncol(hills1$hillsfile)==8) {
if(hills1$per[2]!=hills2$per[2]) {
stop("Error: You can sum only hills of same periodicity")
}
}
hills<-list(hillsfile=rbind(hills1$hillsfile, hills2$hillsfile), size=dim(rbind(hills1$hillsfile, hills2$hillsfile)),
filename=hills1$filename, per=hills1$per, pcv1=hills1$pcv1, pcv2=hills1$pcv2)
class(hills) <- "hillsfile"
return(hills)
}
#' Plot hillsfile object
#'
#' `plot.hillsfile` plots hillsfile object. For a hillsfile with one collective variable it plots its evolution.
#' For a hillsfile with two collective variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'.
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
plot.hillsfile<-function(x, ignoretime=FALSE,
xlab=NULL, ylab=NULL,
xlim=NULL, ylim=NULL,
main=NULL, sub=NULL,
pch=1, col="black", bg="red", cex=1,
asp=NULL, lwd=1, axes=TRUE,...) {
hills <-x
xlims<-NULL
ylims<-NULL
if(!is.null(xlim)) {xlims<-xlim}
if(!is.null(ylim)) {ylims<-ylim}
if(hills$size[2]==5) {
if((hills$per[1]==T)&is.null(ylim)) {ylims<-hills$pcv1}
if(is.null(xlab)) xlab="time"
if(is.null(ylab)) ylab="CV"
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
col=col, cex=cex, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,2], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
col=col, cex=cex, lwd=lwd,
asp=asp, axes=axes)
}
}
if(hills$size[2]==7) {
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
plot(hills$hillsfile[,2], hills$hillsfile[,3], type="p",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
xlim=xlims, ylim=ylims,
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd,
asp=asp, axes=axes)
}
}
#' Plot points for hillsfile object
#'
#' `points.hillsfile` plots points for hillsfile object. For a hillsfile with one
#' collective variable it plots its evolution. For a hillsfile with two collective
#' variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'.
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
#' points(acealanme, col="red")
points.hillsfile<-function(x, ignoretime=FALSE,
pch=1, col="black", bg="red", cex=1,
lwd=1, ...) {
hills <- x
if(hills$size[2]==5) {
if(ignoretime) {
points(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2],
col=col, cex=cex, lwd=lwd)
} else {
points(hills$hillsfile[,1], hills$hillsfile[,2],
col=col, cex=cex, lwd=lwd)
}
}
if(hills$size[2]==7) {
points(hills$hillsfile[,2], hills$hillsfile[,3],
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd)
}
}
#' Plot lines for hillsfile object
#'
#' `lines.hillsfile` plots lines for hillsfile object. For a hillsfile with one
#' collective variable it plots its evolution. For a hillsfile with two collective
#' variables it plots CV1 vs. CV2.
#'
#' @param x hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme)
#' lines(acealanme, col="red")
lines.hillsfile<-function(x, ignoretime=FALSE,
lwd=1, col="black",...) {
hills <- x
if(hills$size[2]==5) {
if(ignoretime) {
lines(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,2],
col=col, lwd=lwd)
} else {
lines(hills$hillsfile[,1], hills$hillsfile[,2],
col=col, lwd=lwd)
}
}
if(hills$size[2]==7) {
lines(hills$hillsfile[,2], hills$hillsfile[,3],
col=col, lwd=lwd)
}
}
#' Plot evolution of heights of hills (generic function for 'metadynminer' and
#' 'metadynminer3d')
#'
#' `plotheights` plots evolution of heights of hills. In well tempered metadynamics
#' hill heights decrees with flooding of the free energy surface. Evolution of heights
#' may be useful to evaluate convergence of the simulation.
#'
#' @param hills hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#'
#' @export plotheights
plotheights<-function(hills, ignoretime, xlab, ylab,
xlim, ylim, main, sub,
col, asp, lwd, axes) {
UseMethod("plotheights")
}
#' Plot evolution of heights of hills in hillsfile object
#'
#' `plotheights.hillsfile` plots evolution of heights of hills. In well tempered metadynamics
#' hill heights decrees with flooding of the free energy surface. Evolution of heights
#' may be useful to evaluate convergence of the simulation.
#'
#' @param hills hillsfile object.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#'
#' @export
#' @examples
#' plotheights(acealanme)
plotheights.hillsfile<-function(hills, ignoretime=FALSE,
xlab=NULL, ylab=NULL,
xlim=NULL, ylim=NULL,
main=NULL, sub=NULL,
col="black", asp=NULL, lwd=1, axes=TRUE) {
if(is.null(xlab)) xlab="time"
if(is.null(ylab)) ylab="hill height"
if(hills$size[2]==5) {
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,4], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,4], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
if(hills$size[2]==7) {
if(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,6], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
} else {
plot(hills$hillsfile[,1], hills$hillsfile[,6], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
}
#' Calculate free energy surface by Bias Sum algorithm (generic function for
#' 'metadynminer' and 'metadynminer3d')
#'
#' `fes` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export fes
fes<-function(hills, imin, imax, xlim, ylim, zlim, npoints) {
UseMethod("fes")
}
#' Calculate free energy surface by Bias Sum algorithm
#'
#' `fes.hillsfile` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
fes.hillsfile<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imax>0 & imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
if(is.null(ylim)) {
minCV2 <- 0
maxCV2 <- 1
} else {
minCV2 <- ylim[1]
maxCV2 <- ylim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
fesm <- matrix(rep(0, npoints*npoints), nrow=npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
if(max(hills$hillsfile[,4])/min(hills$hillsfile[,4])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
if(hills$per[1]==T) {
if(min(hills$hillsfile[,2])<xlims[1]) {
stop("Error: The first collective variable outside pcv1")
}
if(max(hills$hillsfile[,2])>xlims[2]) {
stop("Error: The first collective variable outside pcv1")
}
}
if(hills$per[2]==T) {
if(min(hills$hillsfile[,3])<ylims[1]) {
stop("Error: The second collective variable outside pcv2")
}
if(max(hills$hillsfile[,3])>ylims[2]) {
stop("Error: The second collective variable outside pcv2")
}
}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
if(hills$size[2]==5) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
fesm <- rep(0, npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
if(max(hills$hillsfile[,3])/min(hills$hillsfile[,3])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
if(hills$per[1]==F) {
fesm<-hills1d1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,3])/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
if(hills$per[1]==T) {
fesm<-hills1d1p(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,3])/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
return(cfes)
}
#' Calculate free energy surface by conventional algorithm (generic function
#' for 'metadynminer' and 'metadynminer3d')
#'
#' `fes2` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export fes2
fes2<-function(hills, imin, imax, xlim, ylim, zlim, npoints) {
UseMethod("fes2")
}
#' Calculate free energy surface by conventional algorithm
#'
#' `fes2.hillsfile` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile object.
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param zlim numeric vector of length 2, giving the CV3 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes2(acealanme, imax=1000)
fes2.hillsfile<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imax>0 & imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
if(is.null(ylim)) {
minCV2 <- 0
maxCV2 <- 1
} else {
minCV2 <- ylim[1]
maxCV2 <- ylim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
fesm <- matrix(rep(0, npoints*npoints), nrow=npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills2p1((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills2p2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills2p12((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,4]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,5]/(ylims[2]-ylims[1]),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=2, per=hills$per, x=x, y=y, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
if(hills$size[2]==5) {
if(imax==0) {
if(is.null(xlim)) {
minCV1 <- 0
maxCV1 <- 1
} else {
minCV1 <- xlim[1]
maxCV1 <- xlim[2]
}
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
fesm <- rep(0, npoints)
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
} else {
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
if(hills$per[1]==F) {
fesm<-hills1d2((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,3]/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
if(hills$per[1]==T) {
fesm<-hills1d2p((npoints-1)*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,3]/(xlims[2]-xlims[1]),
hills$hillsfile[,4],npoints,imin-1,imax-1)
}
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per, x=x, pcv1=hills$pcv1, pcv2=hills$pcv2)
class(cfes) <- "fes"
}
}
return(cfes)
}
#' Read 1D or 2D free energy surface from PLUMED sum_hills
#'
#' `read.plumed` reads 1D or 2D free energy surface from PLUMED sum_hills.
#' The grid in the (2D) inputfile must contain the same number of points
#' for CV1 and CV2. It does not use the header of the file. Instead, user
#' must specify the dimensionality (1 or 2). Periodicity must be specified
#' as well.
#'
#' @param file input file from PLUMED sum_hills.
#' @param dim dimension (1 or 2, default 2).
#' @param per logical vector specifying periodicity of collective variables.
#' @return fes object.
#'
#' @export
#' @examples
#' l1<-"-3.142 -124.8 -44.76"
#' l2<-"-3.117 -125.9 -43.05"
#' l3<-"-3.092 -126.9 -41.22"
#' l4<-"-3.068 -127.9 -39.36"
#' l5<-"-3.043 -128.8 -37.45"
#' fourpoints<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourpoints, tf)
#' read.plumed(tf, dim=1, per=c(TRUE,TRUE))
read.plumed<-function(file="fes.dat", dim=2, per=c(F,F,F)) {
hillsf<-read.table(file, header=F, comment.char="#")
bins<-round(nrow(hillsf)^(1/dim))
if(bins^dim!=nrow(hillsf)) {
stop("Error: the number of bins cannot be determined, it must be same for all dimension, or the number of dimensions is wrong.")
}
if(dim==1) {
x <- hillsf[,1]
fesm <- hillsf[,2]
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=1, per=per, x=x, pcv1=c(min(x), max(x)))
class(cfes) <- "fes"
} else if(dim==2) {
x <- hillsf[1:bins,1]
y <- hillsf[(0:(bins-1))*bins+1,2]
fesm <- matrix(hillsf[,3], nrow=bins)
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=2, per=per, x=x, y=y,
pcv1=c(min(x), max(x)), pcv2=c(min(y), max(y)))
class(cfes) <- "fes"
} else {
stop("Error: for 3D fes use read.plumed3d from metadynminer3d, higher dimensions are not supported.")
}
return(cfes)
}
#' Calculate 1D free energy surface from hillsfile object
#'
#' `fes2d21d` calculates 2D free energy surface, converts free energies to probabilities
#' (exp(-F/kT)), sums them up along one collective variable and converts back to free
#' energy (-kT log(P)).
#'
#' @param hills hillsfile object.
#' @param remdim dimension to be removed (1 for CV1, 2 for CV2, default 2).
#' @param temp temperature in Kelvins (default 300).
#' @param eunit energy units (kJ/mol or kcal/mol, kJ/mol is default).
#' @param imin index of a hill from which summation starts (default 1).
#' @param imax index of a hill from which summation stops (default the rest of hills).
#' @param xlim numeric vector of length 2, giving the CV1 coordinates range.
#' @param ylim numeric vector of length 2, giving the CV2 coordinates range.
#' @param npoints resolution of the free energy surface in number of points.
#' @return fes object.
#'
#' @export
#' @examples
#' tfes<-fes2d21d(acealanme, remdim=2, imax=5000)
fes2d21d<-function(hills, remdim=2, temp=300, eunit="kJ/mol",
imin=1, imax=NULL, xlim=NULL, ylim=NULL, npoints=256) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
cat("Warning: You requested more hills by imax than available, using all hills\n")
imax<-hills$size[1]
}
}
if(is.null(imax)) {
imax<-hills$size[1]
}
if(imin>imax) {
stop("Error: imax cannot be lower than imin")
}
if(hills$size[2]==7) {
if(max(hills$hillsfile[,4])/min(hills$hillsfile[,4])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
if(!is.null(xlim)) {xlims<-xlim}
if((hills$per[1]==T)&is.null(xlim)) {xlims<-hills$pcv1}
if(!is.null(ylim)) {ylims<-ylim}
if((hills$per[2]==T)&is.null(ylim)) {ylims<-hills$pcv2}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
binx<-(xlims[2]-xlims[1])/(npoints-1)
biny<-(ylims[2]-ylims[1])/(npoints-1)
if((hills$per[1]==F)&(hills$per[2]==F)) {
fesm<-hills1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)) {
fesm<-hills1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)) {
fesm<-hills1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)) {
fesm<-hills1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]+binx),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]+biny),
npoints*max(hills$hillsfile[,4])/(xlims[2]-xlims[1]+binx),
npoints*max(hills$hillsfile[,5])/(ylims[2]-ylims[1]+biny),
hills$hillsfile[,6],npoints,imin-1,imax-1)
}
if(eunit=="kJ/mol") {
prob<- exp(-1000*fesm/8.314/temp)
if(remdim==1) {
fesm <- -8.314*temp*log(apply(prob, 2, sum))/1000
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[2], x=y, pcv1=hills$pcv2)
}
if(remdim==2) {
fesm <- -8.314*temp*log(apply(prob, 1, sum))/1000
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[1], x=x, pcv1=hills$pcv1)
}
}
if(eunit=="kcal/mol") {
prob<- exp(-1000*4.184*fesm/8.314/temp)
if(remdim==1) {
fesm <- -8.314*temp*log(apply(prob, 2, sum))/1000/4.184
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[2], x=y, pcv1=hills$pcv2)
}
if(remdim==2) {
fesm <- -8.314*temp*log(apply(prob, 1, sum))/1000/4.184
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=1, per=hills$per[1], x=x, pcv1=hills$pcv1)
}
}
class(cfes) <- "fes"
}
if(hills$size[2]==5) {
stop("Error: Your free energy surface is already 1D")
}
return(cfes)
}
#' @export
`+.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(fes1$dimension!=fes2$dimension) {
stop("Error: Free energy surfaces have different dimension, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(fes1$dimension==2) {
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
}
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes+fes2$fes, hills=rbind(fes1$hills, fes2$hills), rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes+fes2$fes, hills=rbind(fes1$hills, fes2$hills), rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1+fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1+fes2$fes, hills=rbind(fes1$hills,fes2$hills), rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
class(cfes) <- "fes"
return(cfes)
}
#' @export
`-.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(fes1$dimension!=fes2$dimension) {
stop("Error: Free energy surfaces have different dimension, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(fes1$dimension==2) {
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
}
cat("Warning: FES obtained by subtraction of two FESes\n")
cat(" will inherit hills only from the first FES\n")
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes-fes2$fes, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes-fes2$fes, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
class(cfes) <- "fes"
return(cfes)
}
#' @export
`*.fes`<-function(fes1, fes2) {
if(inherits(fes1, "fes")&inherits(fes2, "fes")) {
stop("Error: You cannot multiply fes by fes")
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(fes2, "fes")) {
if(fes2$dimension==1) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
if(fes2$dimension==2) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, pcv1=fes2$pcv1, pcv2=fes2$pcv2)
}
}
cat("Warning: multiplication of FES will multiply\n")
cat(" the FES but not hill heights\n")
class(cfes) <- "fes"
return(cfes)
}
#' @export
`/.fes`<-function(fes1, coef) {
if(inherits(fes1, "fes")&inherits(coef, "fes")) {
stop("Error: You cannot divide fes by fes")
} else if(inherits(fes1, "fes")) {
if(fes1$dimension==1) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
if(fes1$dimension==2) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, pcv1=fes1$pcv1, pcv2=fes1$pcv2)
}
} else if(inherits(coef, "fes")) {
stop("Error: You cannot divide something by fes")
}
cat("Warning: division of FES will divide\n")
cat(" the FES but not hill heights\n")
class(cfes) <- "fes"
return(cfes)
}
#' Calculate minimum of free energy surface
#'
#' `min.fes` calculates minimum of free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param na.rm a logical indicating whether missing values should be
#' removed.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' min(tfes)
min.fes<-function(inputfes, na.rm=NULL,...) {
return(min(inputfes$fes, na.rm=na.rm))
}
#' Calculate maximum of free energy surface
#'
#' `max.fes` calculates maximum of free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param na.rm a logical indicating whether missing values should be
#' removed.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' max(tfes)
max.fes<-function(inputfes, na.rm=NULL,...) {
return(max(inputfes$fes, na.rm=na.rm))
}
#' Print dimensionality, minimum and maximum of free energy surface
#'
#' `print.fes` prints dimensionality, minimum and maximum of
#' free energy in a fes object
#'
#' @param x fes object
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' tfes
print.fes<-function(x,...) {
inputfes<-x
if(inputfes$dimension==1) {
cat("1D free energy surface with ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
if(inputfes$dimension==2) {
cat("2D free energy surface with ")
cat(inputfes$rows)
cat(" x ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
}
#' Print summary of free energy surface
#'
#' `summary.fes` prints dimensionality, minimum and maximum of
#' free energy in a fes object.
#'
#' @param object fes object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' summary(tfes)
summary.fes<-function(object,...) {
inputfes <- object
if(inputfes$dimension==1) {
cat("1D free energy surface with ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
if(inputfes$dimension==2) {
cat("2D free energy surface with ")
cat(inputfes$rows)
cat(" x ")
cat(inputfes$rows)
cat(" points, maximum ")
cat(max(inputfes$fes))
cat(" and minimum ")
cat(min(inputfes$fes))
cat("\n")
}
}
#' Calculate probability of free energy surface
#'
#' `prob` calculates probability from free energy in a fes object.
#'
#' @param inputfes fes object.
#' @param temp temperature in Kelvins (default 300).
#' @param eunit energy units (kJ/mol or kcal/mol, kJ/mol is default).
#'
#' @export
#' @examples
#' tfes<-fes(acealanme, imax=5000)
#' print(prob(tfes))
prob<-function(inputfes, temp=300, eunit="kJ/mol") {
if(inherits(inputfes, "fes")) {
if(eunit=="kJ/mol") {
if(inputfes$dimension==1) {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, pcv1=inputfes$pcv1)
}
if(inputfes$dimension==2) {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2)
}
} else if (eunit=="kJ/mol") {
if(inputfes$dimension==1) {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, pcv1=inputfes$pcv1)
}
if(inputfes$dimension==2) {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2)
}
} else {
stop("Error: Wrong eunit")
}
class(cfes) <- "fes"
return(cfes)
} else if(inherits(inputfes, "fes3d")) {
if(eunit=="kJ/mol") {
probs <- exp(-1000*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, z=inputfes$z, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2, pcv3=inputfes$pcv3)
} else if (eunit=="kJ/mol") {
probs <- exp(-1000*4.184*inputfes$fes/8.314/temp)
cfes<-list(fes=probs/sum(probs), hills=inputfes$hills, rows=inputfes$rows, dimension=inputfes$dimension, per=inputfes$per, x=inputfes$x, y=inputfes$y, z=inputfes$z, pcv1=inputfes$pcv1, pcv2=inputfes$pcv2, pcv3=inputfes$pcv3)
} else {
stop("Error: Wrong eunit")
}
class(cfes) <- "fes3d"
return(cfes)
} else {
stop("Error: Input must be fes or fes3d object")
}
}
#' Plot free energy surface object
#'
#' `plot.fes` plots free energy surface. For a fes with one collective variable it plots
#' a 1D profile. For a fes with two collective variables it plots 2D free energy surface
#' using image, contours or combination of both (default).
#'
#' @param x fes object.
#' @param plottype specifies whether 2D free energy surface will be plotted
#' as image, contours or both (default "both").
#' @param colscale specifies whether color scale will be plotted (default False).
#' @param colscalelab color scale label (default "free energy").
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param asp the y/x aspect ratio, see 'plot.window'.
#' @param col color of the free energy surface. For 1D surface it is the color
#' of the line. For 2D it is a list of colors such as that generated by
#' 'rainbow', 'heat.colors', 'topo.colors', 'terrain.colors' or similar
#' functions (default=rainbow(135)[100:1]).
#' @param xlim numeric vector of length 2, giving the x coordinates range.
#' @param ylim numeric vector of length 2, giving the y coordinates range.
#' @param zlim numeric vector of length 2, giving the z coordinates range.
#' @param nlevels number of contour levels desired if 'levels' is not
#' supplied.
#' @param levels numeric vector of levels at which to draw contour lines.
#' @param labels a vector giving the labels for the contour lines. If 'NULL'
#' then the levels are used as labels, otherwise this is coerced
#' by 'as.character'.
#' @param labcex 'cex' for contour labeling. This is an absolute size, not a
#' multiple of 'par("cex")'.
#' @param drawlabels logical. Contours are labeled if 'TRUE'.
#' @param method character string specifying where the labels will be located.
#' Possible values are '"simple"', '"edge"' and '"flattest"'
#' (the default). See the 'Details' section.
#' @param lwd contour line width.
#' @param contcol contour color.
#' @param lty line type for the lines drawn.
#' @param axes a logical value indicating whether both axes should be drawn
#' on the plot.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes2d<-fes(acealanme, imax=5000)
#' plot(tfes2d)
#' tfes1d<-fes(acealanme1d)
#' plot(tfes1d)
plot.fes<-function(x, plottype="both",
colscale=F, xlim=NULL, ylim=NULL, zlim=NULL,
main=NULL, sub=NULL,
xlab=NULL, ylab=NULL,
nlevels=10, levels=NULL,
col=rainbow(135)[100:1],
labels=NULL, labcex=0.6, drawlabels=TRUE,
colscalelab="free energy",
method="flattest",
contcol=par("fg"), lty=par("lty"),
lwd=1, asp=NULL, axes=T,...) {
close.screen(all.screens=TRUE)
inputfes<-x
fes<-inputfes$fes
rows<-inputfes$rows
if(inputfes$dimension==1) {
x<-inputfes$x
if(is.null(xlab)) xlab="CV"
if(is.null(ylab)) ylab="free energy"
if(is.null(xlim)) xlim<-c(min(x),max(x))
if(is.null(ylim)) {
ylim<-range(pretty(range(fes)))
}
plot(x, fes, type="l", lwd=lwd,
col=col, xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, axes=axes,
main=main, sub=sub, asp=asp)
} else {
x<-inputfes$x
y<-inputfes$y
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
if(is.null(zlim)) {
zlim<-range(pretty(range(fes)))
}
if(is.null(levels)) {
levels<-pretty(zlim, nlevels)
}
if(is.null(xlim)) xlim<-c(min(x),max(x))
if(is.null(ylim)) ylim<-c(min(y),max(y))
if(colscale) {
#layout(matrix(c(1,2), 1, 2, byrow = TRUE), widths=c(4,1))
split.screen(matrix(c(0,0.75,0,1,0.75,1,0,1), byrow=T, ncol=4))
screen(2)
smat<-matrix(seq(from=zlim[1], to=zlim[2], length.out=100))
image(c(0), seq(from=zlim[1], to=zlim[2], length.out=100),
t(smat), zlim=zlim, col=col, xlab="", ylab=colscalelab, axes=F)
axis(2, lty=lty, lwd=lwd)
box(lwd=lwd)
screen(1)
}
if(plottype=="image" || plottype=="both") {
image(x, y, fes, zlim=zlim,
col=col, xlim=xlim, ylim=ylim,
xlab=xlab, ylab=ylab, axes=axes,
main=main, sub=sub, asp=asp)
}
if(plottype=="contour") {
contour(x, y, fes, zlim=zlim,
nlevels=nlevels, levels=levels,
labels=labels, labcex=labcex, drawlabels=drawlabels,
method=method, col=contcol, lty=lty, lwd=lwd,
main=main, sub=sub, asp=asp)
}
if(plottype=="both") {
contour(x, y, fes, zlim=zlim,
nlevels=nlevels, levels=levels,
labels=labels, labcex=labcex, drawlabels=drawlabels,
method=method, col=contcol, lty=lty, lwd=lwd, add=T)
}
}
}
#' Plots 1D free energy surface object as points
#'
#' `points.fes` plots 1D free energy surface as points.
#'
#' @param x fes object.
#' @param pch plotting 'character', i.e., symbol to use. See 'points'
#' @param col color code or name, see 'par'.
#' @param bg background (fill) color for the open plot symbols given by
#' 'pch = 21:25'.
#' @param cex character (or symbol) expansion: a numerical vector. This
#' works as a multiple of 'par("cex")'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme1d, imax=5000)
#' plot(tfes)
#' points(tfes)
points.fes<-function(x, pch=1, col="black", bg="red", cex=1, lwd=1,...) {
fes<-x$fes
if(x$dimension==1) {
x<-x$x
points(x, fes,
pch=pch, col=col, bg=bg, cex=cex, lwd=lwd)
} else {
stop("Error: points available only for 1D free energy surfaces\n")
}
}
#' Plots 1D free energy surface object as lines
#'
#' `lines.fes` plots 1D free energy surface as lines.
#'
#' @param x fes object.
#' @param col color code or name, see 'par'.
#' @param lwd line width for drawing symbols see 'par'.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme1d, imax=5000)
#' plot(tfes)
#' lines(tfes, lwd=4)
lines.fes<-function(x, lwd=1, col="black",...) {
fes<-x$fes
if(x$dimension==1) {
x<-x$x
lines(x, fes, lwd=lwd, col=col)
} else {
stop("Error: points available only for 1D free energy surfaces\n")
}
}
Try the metadynminer package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.