Nothing
R/readingandfes.R
In metadynminer3d: Tools to Read, Analyze and Visualize Metadynamics 3D HILLS Files from 'Plumed'
Defines functions
plot.fes3d
summary.fes3d
print.fes3d
max.fes3d
min.fes3d
`/.fes3d`
`*.fes3d`
`-.fes3d`
`+.fes3d`
read.plumed3d
fes2.hillsfile3d
fes.hillsfile3d
plotheights.hillsfile3d
plot.hillsfile3d
`+.hillsfile3d`
tail.hillsfile3d
head.hillsfile3d
summary.hillsfile3d
print.hillsfile3d
read.hills3d2
read.hills3d
Documented in
fes2.hillsfile3d
fes.hillsfile3d
head.hillsfile3d
max.fes3d
min.fes3d
plot.fes3d
plotheights.hillsfile3d
plot.hillsfile3d
print.fes3d
print.hillsfile3d
read.hills3d
read.plumed3d
summary.fes3d
summary.hillsfile3d
tail.hillsfile3d
#' @useDynLib metadynminer3d
#' @import metadynminer
#' @importFrom Rcpp sourceCpp
#' @importFrom misc3d contour3d
#' @importFrom rgl axes3d title3d box3d text3d points3d
NULL
#' Read 3D HILLS from Plumed
#'
#' `read.hills3d` reads a HILLS file generated by Plumed and returns a hillsfile3d 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 pcv3 periodicity of CV3.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.587 -2.969 3.013 0.3 0.3 0.3 1.111 10"
#' l2<-"2 -1.067 2.745 2.944 0.3 0.3 0.3 1.109 10"
#' l3<-"3 -1.376 2.697 3.049 0.3 0.3 0.3 1.080 10"
#' l4<-"4 -1.663 2.922 -3.065 0.3 0.3 0.3 1.072 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills3d(tf, per=c(TRUE,TRUE))
read.hills3d<-function(file="HILLS", per=c(FALSE, FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), pcv3=c(-pi,pi), ignoretime=FALSE) {
hillsf<-read.table(file, header=F, comment.char="#")
if(ncol(hillsf)==5 || ncol(hillsf)==6) {
warning("looks like you trying to load 1D HILLS, use read.hills from metadynminer\n")
} else if(ncol(hillsf)==7 || ncol(hillsf)==8) {
warning("looks like you trying to load 2D HILLS, use read.hills from metadynminer\n")
} else if(ncol(hillsf)==9 || ncol(hillsf)==10) {
cat("3D HILLS file read\n")
if(ignoretime) {
warning("Warning: The time will be updated automatically from zero 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], cv3=hillsf[,4],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2, pcv3=pcv3)
class(hills) <- "hillsfile3d"
return(hills)
} else {
stop("Error: Number of columns in HILLS file must be 9 or 10 (3D)")
}
}
#' Read 3D HILLS from any program
#'
#' `read.hills3d2` reads a HILLS file generated by any program and returns a hillsfile3d object.
#' User can specify whether some collective variables are periodic.
#'
#' @param file HILLS file from any program.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @param pcv3 periodicity of CV3.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param times the index of a column containing time values.
#' @param centres vector of indexes of columns containing centres of hills.
#' @param widths vector of indexes of columns containing widths of hills.
#' @param heights the index of a column containing heights of hills.
#' @param sep the field separator character.
#' @param dec the character used in the file for decimal points.
#' @param na.strings a character vector of strings which are to be interpreted as NA values.
#' @param comment.char a character vector of length one containing a single character or an empty string.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.587 -2.969 3.013 0.3 0.3 0.3 1.111 10"
#' l2<-"2 -1.067 2.745 2.944 0.3 0.3 0.3 1.109 10"
#' l3<-"3 -1.376 2.697 3.049 0.3 0.3 0.3 1.080 10"
#' l4<-"4 -1.663 2.922 -3.065 0.3 0.3 0.3 1.072 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills3d2(tf, per=c(TRUE,TRUE), times=1, centres=2:4, widths=5:7, heights=8)
read.hills3d2<-function(file="HILLS", per=c(FALSE, FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), pcv3=c(-pi,pi),
ignoretime=FALSE, times=1, centres=2:4, widths=5:7, heights=8,
sep="", dec=".", na.strings="NA", comment.char="#") {
hillsf<-read.table(file, header=F, sep=sep, dec=dec, na.strings=na.strings, comment.char=comment.char)
if(length(times)!=1) {
stop("Error: Parameter values must be set to one value")
}
if(length(centres)!=3) {
stop("Error: Parameter centres must be set to three values")
}
if(length(widths)!=3) {
stop("Error: Parameter widths must be set to three values")
}
if(length(heights)!=1) {
stop("Error: Parameter heights must be set to one value")
}
cat("3D HILLS file read\n")
newhillsf <- data.frame(rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)))
newhillsf[,1]<-hillsf[,times]
newhillsf[,2]<-hillsf[,centres[1]]
newhillsf[,3]<-hillsf[,centres[2]]
newhillsf[,4]<-hillsf[,centres[3]]
newhillsf[,5]<-hillsf[,widths[1]]
newhillsf[,6]<-hillsf[,widths[2]]
newhillsf[,7]<-hillsf[,widths[3]]
newhillsf[,8]<-hillsf[,heights]
if(ignoretime) {
warning("Warning: The time will be updated automatically from zero according to the first step!\n")
newhillsf[,1]<-seq(from=hillsf[1,times], by=hillsf[1,times], length.out=nrow(hillsf))
}
hills<-list(hillsfile=newhillsf, time=newhillsf[,1], cv1=newhillsf[,2], cv2=newhillsf[,3], cv3=newhillsf[,4],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2, pcv3=pcv3)
class(hills) <- "hillsfile3d"
return(hills)
}
#' Print hillsfile3d
#'
#' `print.hillsfile3d` prints dimensionality and size of a hillsfile object.
#'
#' @param x hillsfile3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' acealanme3d
print.hillsfile3d<-function(x,...) {
hills <- x
toprint <- "3D hills file "
toprint <- paste(toprint, toString(hills$filename), sep="")
toprint <- paste(toprint, " with ", sep="")
toprint <- paste(toprint, toString(hills$size[1]), sep="")
toprint <- paste(toprint, " lines", sep="")
message(toprint)
}
#' Print summary for hillsfile3d
#'
#' `summary.hillsfile3d` prints dimensionality, size and collective variable ranges of a hillsfile3d object.
#'
#' @param object hillsfile3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' summary(acealanme3d)
summary.hillsfile3d<-function(object,...) {
hills <- object
toprint <- "3D hills file "
toprint <- paste(toprint, toString(hills$filename), sep="")
toprint <- paste(toprint, " with ", sep="")
toprint <- paste(toprint, toString(hills$size[1]), sep="")
toprint <- paste(toprint, " lines", sep="")
message(toprint)
toprint <- "The CV1 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,2])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,2])), sep="")
message(toprint)
toprint <- "The CV2 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,3])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,3])), sep="")
message(toprint)
toprint <- "The CV3 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,4])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,4])), sep="")
message(toprint)
}
#' Print first n lines of hillsfile3d
#'
#' `head.hillsfile3d` prints first n lines of a hillsfile3d object.
#'
#' @param x hillsfile3d object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' head(acealanme3d)
head.hillsfile3d<-function(x, n=10,...) {
return(head(x$hillsfile, n=n))
}
#' Print last n lines of hillsfile3d
#'
#' `tail.hillsfile3d` prints last n lines of a hillsfile3d object.
#'
#' @param x hillsfile3d object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tail(acealanme3d)
tail.hillsfile3d<-function(x, n=10,...) {
return(tail(x$hillsfile, n=n))
}
#' @export
`+.hillsfile3d`<-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")
}
}
if(ncol(hills1$hillsfile)==9 || ncol(hills1$hillsfile)==10) {
if(hills1$per[2]!=hills2$per[2]) {
stop("Error: You can sum only hills of same periodicity")
}
if(hills1$per[3]!=hills2$per[3]) {
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, pcv3=hills1$pcv3)
class(hills) <- "hillsfile3d"
return(hills)
}
#' Plot hillsfile3d object
#'
#' `plot.hillsfile3d` plots hillsfile object. It plots CV1 vs CV2 vs CV3.
#'
#' @param x hillsfile object.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param zlab a title for the z axis: see 'title'.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param col color code or name, 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 zlim numeric vector of length 2, giving the z coordinates range.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme3d)
plot.hillsfile3d<-function(x,
xlab="CV1", ylab="CV2", zlab="CV3",
main=NULL, sub=NULL,
col="orange",...) {
hills <-x
points3d(hills$hillsfile[,2:4], col=col)
axes3d()
title3d(xlab=xlab, ylab=ylab, zlab=zlab,
main=main, sub=sub)
box3d()
}
#' Plot evolution of heights of hills in hillsfile3d object
#'
#' `plotheights.hillsfile3d` 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(acealanme3d)
plotheights.hillsfile3d<-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(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,8], 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[,8], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
#' Calculate 3D free energy surface by Bias Sum algorithm
#'
#' `fes.hillsfile3d` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile3d 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 x lim 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(acealanme3d, imax=5000)
fes.hillsfile3d<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=NULL) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
warning("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(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,6])/min(hills$hillsfile[,6])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,7])/min(hills$hillsfile[,7])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(is.null(npoints)) {
npoints <- 64
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
minCV3 <- min(hills$hillsfile[,4])
maxCV3 <- max(hills$hillsfile[,4])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
zlims<-c(minCV3-0.05*(maxCV3-minCV3), maxCV3+0.05*(maxCV3-minCV3))
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(!is.null(zlim)) {zlims<-zlim}
if((hills$per[3]==T)&is.null(zlim)) {zlims<-hills$pcv3}
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")
}
}
if(hills$per[3]==T) {
if(min(hills$hillsfile[,4])<zlims[1]) {
stop("Error: The third collective variable outside pcv3")
}
if(max(hills$hillsfile[,4])>zlims[2]) {
stop("Error: The third collective variable outside pcv3")
}
}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
z<-0:(npoints-1)*(zlims[2]-zlims[1])/(npoints-1)+zlims[1]
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d1p3(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d1p13(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d1p23(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d1p123(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
fesm<-aperm(array(fesm, c(npoints, npoints, npoints)))
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=3, per=hills$per,
x=x, y=y, z=z, pcv1=hills$pcv1, pcv2=hills$pcv2, pcv3=hills$pcv3)
class(cfes) <- "fes3d"
return(cfes)
}
#' Calculate 3D free energy surface by conventional algorithm
#'
#' `fes2.hills3d` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile3d 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(acealanme3d, imax=100)
fes2.hillsfile3d<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=NULL) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
warning("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(is.null(npoints)) {
npoints <- 64
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
minCV3 <- min(hills$hillsfile[,4])
maxCV3 <- max(hills$hillsfile[,4])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
zlims<-c(minCV3-0.05*(maxCV3-minCV3), maxCV3+0.05*(maxCV3-minCV3))
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(!is.null(zlim)) {zlims<-zlim}
if((hills$per[3]==T)&is.null(zlim)) {zlims<-hills$pcv3}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
z<-0:(npoints-1)*(zlims[2]-zlims[1])/(npoints-1)+zlims[1]
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d2((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d2p1((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d2p2((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d2p3((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d2p12((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d2p13((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d2p23((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d2p123((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
fesm<-aperm(array(fesm, c(npoints, npoints, npoints)))
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=3, per=hills$per,
x=x, y=y, z=z, pcv1=hills$pcv1, pcv2=hills$pcv2, pcv3=hills$pcv3)
class(cfes) <- "fes3d"
return(cfes)
}
#' Read 3D free energy surface from PLUMED sum_hills
#'
#' `read.plumed3d` reads 3D free energy surface from PLUMED sum_hills.
#' The grid in the input file must contain the same number of points
#' for CV1, CV2 and CV3. It does not use the header of the file.
#' Periodicity must be specified.
#'
#' @param file input file from PLUMED sum_hills.
#' @param per logical vector specifying periodicity of collective variables.
#' @return fes object.
#'
#' @export
#' @examples
#' l1<-" -3.14 -3.14 -3.14 -61.13 -47.43 19.00 2.04"
#' l2<-" -1.05 -3.14 -3.14 -70.72 25.95 25.78 2.43"
#' l3<-" 1.05 -3.14 -3.14 -65.58 8.34 2.82 -3.09"
#' l4<-" -3.14 -1.05 -3.14 -51.31 -43.88 -19.91 1.51"
#' l5<-" -1.05 -1.05 -3.14 -66.43 7.67 -22.45 -0.39"
#' l6<-" 1.05 -1.05 -3.14 -61.08 -7.50 -7.36 -0.83"
#' l7<-" -3.14 1.05 -3.14 -53.07 -55.12 0.19 -0.28"
#' l8<-" -1.05 1.05 -3.14 -62.81 36.19 1.65 0.45"
#' l9<-" 1.05 1.05 -3.14 -65.28 22.84 11.47 0.59"
#' l10<-" -3.14 -3.14 -1.05 -13.03 -32.17 8.24 -35.25"
#' l11<-" -1.05 -3.14 -1.05 -21.88 17.89 21.91 -51.20"
#' l12<-" 1.05 -3.14 -1.05 -14.49 3.60 6.04 -44.05"
#' l13<-" -3.14 -1.05 -1.05 -2.26 -7.00 -7.01 -10.65"
#' l14<-" -1.05 -1.05 -1.05 -8.21 3.69 -22.89 -28.48"
#' l15<-" 1.05 -1.05 -1.05 -1.10 0.52 3.59 -1.99"
#' l16<-" -3.14 1.05 -1.05 -3.75 -11.70 -5.65 -15.36"
#' l17<-" -1.05 1.05 -1.05 -1.15 5.75 1.05 -2.42"
#' l18<-" 1.05 1.05 -1.05 -10.67 8.23 -10.42 -36.77"
#' l19<-" -3.14 -3.14 1.05 -4.64 -13.79 10.51 14.96"
#' l20<-" -1.05 -3.14 1.05 -7.80 12.24 20.59 23.03"
#' l21<-" 1.05 -3.14 1.05 -5.32 3.46 3.17 21.99"
#' l22<-" -3.14 -1.05 1.05 -2.06 -6.59 0.17 10.04"
#' l23<-" -1.05 -1.05 1.05 -9.69 8.43 -0.97 36.97"
#' l24<-" 1.05 -1.05 1.05 -0.19 -0.44 -0.26 0.91"
#' l25<-" -3.14 1.05 1.05 -7.98 -23.02 3.97 26.98"
#' l26<-" -1.05 1.05 1.05 -4.64 13.66 -9.74 10.15"
#' l27<-" 1.05 1.05 1.05 -13.42 15.78 16.36 41.60"
#' twentysevenpoints<-c(l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,
#' l11,l12,l13,l14,l15,l16,l17,l18,l19,l20,
#' l21,l22,l23,l24,l25,l26,l27)
#' tf <- tempfile()
#' writeLines(twentysevenpoints, tf)
#' read.plumed3d(tf, per=c(TRUE,TRUE,TRUE))
read.plumed3d<-function(file="fes.dat", per=c(FALSE,FALSE,FALSE)) {
hillsf<-read.table(file, header=F, comment.char="#")
bins<-round(nrow(hillsf)^(1/3))
if(bins^3!=nrow(hillsf)) {
stop("Error: the number of bins cannot be determined.")
}
x <- hillsf[1:bins,1]
y <- hillsf[(0:(bins-1))*bins+1,2]
z <- hillsf[(0:(bins-1))*bins*bins+1,3]
fesm <- array(hillsf[,4], c(bins, bins, bins))
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=3, per=per, x=x, y=y, z=z,
pcv1=c(min(x), max(x)), pcv2=c(min(y), max(y)), pcv3=c(min(z), max(z)))
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`+.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
if(sum(fes1$z!=fes2$z)>0) {
stop("Error: Free energy surfaces have different CV3 axes, exiting")
}
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, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
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, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
} else if(inherits(fes1, "fes")) {
stop("Error: cannot sum 3D free energy surfaces with 1D or 2D, exiting")
} else if(inherits(fes2, "fes")) {
stop("Error: cannot sum 3D free energy surfaces with 1D or 2D, exiting")
}
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`-.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
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(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
if(sum(fes1$z!=fes2$z)>0) {
stop("Error: Free energy surfaces have different CV3 axes, exiting")
}
warning("Warning: FES obtained by subtraction of two FESes will inherit hills only from the first FES\n")
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, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
} else if(inherits(fes1, "fes")) {
stop("Error: cannot subtract 3D free energy surfaces with 1D or 2D, exiting")
} else if(inherits(fes2, "fes")) {
stop("Error: cannot subtract 3D free energy surfaces with 1D or 2D, exiting")
}
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`*.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
stop("Error: You cannot multiply fes by fes")
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
}
warning("Warning: multiplication of FES will multiply the FES but not hill heights\n")
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`/.fes3d`<-function(fes1, coef) {
if(inherits(fes1, "fes3d")&inherits(coef, "fes3d")) {
stop("Error: You cannot divide fes by fes")
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(coef, "fes3d")) {
stop("Error: You cannot divide something by fes")
}
warning("Warning: division of FES will divide the FES but not hill heights\n")
class(cfes) <- "fes3d"
return(cfes)
}
#' Calculate minimum of 3D free energy surface
#'
#' `min.fes3d` calculates minimum of free energy in a fes3d object.
#'
#' @param inputfes fes3d 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(acealanme3d, imax=5000)
#' min(tfes)
min.fes3d<-function(inputfes, na.rm=NULL,...) {
return(min(inputfes$fes, na.rm=na.rm))
}
#' Calculate maximum of 3D free energy surface
#'
#' `max.fes3d` calculates maximum of free energy in a fes3d object.
#'
#' @param inputfes fes3d 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(acealanme3d, imax=5000)
#' max(tfes)
max.fes3d<-function(inputfes, na.rm=NULL,...) {
return(max(inputfes$fes, na.rm=na.rm))
}
#' Print minimum and maximum of 3D free energy surface
#'
#' `print.fes3d` prints dimensionality, minimum and maximum of
#' free energy in a fes object
#'
#' @param x fes3d object
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme3d, imax=5000)
#' tfes
print.fes3d<-function(x,...) {
inputfes<-x
toprint <- "3D free energy surface with "
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " points, maximum ", sep="")
toprint <- paste(toprint, toString(max(inputfes$fes)), sep="")
toprint <- paste(toprint, " and minimum ", sep="")
toprint <- paste(toprint, toString(min(inputfes$fes)), sep="")
message(toprint)
}
#' Print summary of 3D free energy surface
#'
#' `summary.fes3d` prints minimum and maximum of free energy
#' in a fes object.
#'
#' @param object fes3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme3d, imax=5000)
#' summary(tfes)
summary.fes3d<-function(object,...) {
inputfes <- object
toprint <- "3D free energy surface with "
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " points, maximum ", sep="")
toprint <- paste(toprint, toString(max(inputfes$fes)), sep="")
toprint <- paste(toprint, " and minimum ", sep="")
toprint <- paste(toprint, toString(min(inputfes$fes)), sep="")
message(toprint)
}
#' Plot 3D free energy surface object
#'
#' `plot.fes3d` plots 3D free energy surface using .
#'
#' @param x fes3d object.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param zlab a title for the z axis: see 'title'.
#' @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 col color of the free energy surface. It can be a single color
#' or a vector with multiple colors for multiple 3D isosurfaces.
#' @param alpha number or numeric vector of alpha levels (transparency) of
#' 3D isosurfaces.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param level number or numeric vector of levels at which to draw 3D isosurface.
#' @param fill a logical value indicating whether 3D isosurface is ploted as
#' solid surface (True) or wireframe (False).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes3d<-fes(acealanme3d, imax=5000)
#' plot(tfes3d)
plot.fes3d<-function(x, xlab=NULL, ylab=NULL, zlab=NULL,
xlim=NULL, ylim=NULL, zlim=NULL,
level=NULL, col=NULL, alpha=NULL,
main=NULL, sub=NULL,
fill=TRUE,...) {
close.screen(all.screens=TRUE)
inputfes<-x
fes<-inputfes$fes
rows<-inputfes$rows
x<-inputfes$x
y<-inputfes$y
z<-inputfes$z
if(is.null(xlim)) {
xlim=range(x)
} else {
fes <- fes[x >= xlim[1] & x <= xlim[2],,]
x <- x[x >= xlim[1] & x <= xlim[2]]
}
if(is.null(ylim)) {
ylim=range(y)
} else {
fes <- fes[,y >= ylim[1] & y <= ylim[2],]
y <- y[y >= ylim[1] & y <= ylim[2]]
}
if(is.null(zlim)) {
zlim=range(z)
} else {
fes <- fes[,,z >= zlim[1] & z <= zlim[2]]
z <- z[z >= zlim[1] & z <= zlim[2]]
}
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
if(is.null(zlab)) zlab="CV3"
if(is.null(level)) level=(max(fes)+min(fes))/2
if(length(level)>1) {
if(is.null(col)) col<-rainbow(1.35*length(level))[length(level):1]
if(is.null(alpha)) {
alpha<-length(level):1/length(level)
level<-sort(level)
}
} else {
if(is.null(col)) col<-"orange"
if(is.null(alpha)) alpha<-1
}
rgl::plot3d(c(), c(), c(), xlim=xlim, ylim=ylim, zlim=zlim,
xlab=xlab, ylab=ylab, zlab=zlab, main=main, sub=sub)
contour3d(f=fes, level=level, x=x, y=y, z=z,
color=col, alpha=alpha, fill=fill, add=T)
rgl::aspect3d(1,1,1)
axes3d()
box3d()
}
Try the metadynminer3d package in your browser
Any scripts or data that you put into this service are public.
metadynminer3d documentation built on April 14, 2022, 5:08 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 plot.fes3d summary.fes3d print.fes3d max.fes3d min.fes3d `/.fes3d` `*.fes3d` `-.fes3d` `+.fes3d` read.plumed3d fes2.hillsfile3d fes.hillsfile3d plotheights.hillsfile3d plot.hillsfile3d `+.hillsfile3d` tail.hillsfile3d head.hillsfile3d summary.hillsfile3d print.hillsfile3d read.hills3d2 read.hills3d
Documented in fes2.hillsfile3d fes.hillsfile3d head.hillsfile3d max.fes3d min.fes3d plot.fes3d plotheights.hillsfile3d plot.hillsfile3d print.fes3d print.hillsfile3d read.hills3d read.plumed3d summary.fes3d summary.hillsfile3d tail.hillsfile3d
#' @useDynLib metadynminer3d
#' @import metadynminer
#' @importFrom Rcpp sourceCpp
#' @importFrom misc3d contour3d
#' @importFrom rgl axes3d title3d box3d text3d points3d
NULL
#' Read 3D HILLS from Plumed
#'
#' `read.hills3d` reads a HILLS file generated by Plumed and returns a hillsfile3d 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 pcv3 periodicity of CV3.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.587 -2.969 3.013 0.3 0.3 0.3 1.111 10"
#' l2<-"2 -1.067 2.745 2.944 0.3 0.3 0.3 1.109 10"
#' l3<-"3 -1.376 2.697 3.049 0.3 0.3 0.3 1.080 10"
#' l4<-"4 -1.663 2.922 -3.065 0.3 0.3 0.3 1.072 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills3d(tf, per=c(TRUE,TRUE))
read.hills3d<-function(file="HILLS", per=c(FALSE, FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), pcv3=c(-pi,pi), ignoretime=FALSE) {
hillsf<-read.table(file, header=F, comment.char="#")
if(ncol(hillsf)==5 || ncol(hillsf)==6) {
warning("looks like you trying to load 1D HILLS, use read.hills from metadynminer\n")
} else if(ncol(hillsf)==7 || ncol(hillsf)==8) {
warning("looks like you trying to load 2D HILLS, use read.hills from metadynminer\n")
} else if(ncol(hillsf)==9 || ncol(hillsf)==10) {
cat("3D HILLS file read\n")
if(ignoretime) {
warning("Warning: The time will be updated automatically from zero 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], cv3=hillsf[,4],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2, pcv3=pcv3)
class(hills) <- "hillsfile3d"
return(hills)
} else {
stop("Error: Number of columns in HILLS file must be 9 or 10 (3D)")
}
}
#' Read 3D HILLS from any program
#'
#' `read.hills3d2` reads a HILLS file generated by any program and returns a hillsfile3d object.
#' User can specify whether some collective variables are periodic.
#'
#' @param file HILLS file from any program.
#' @param per logical vector specifying periodicity of collective variables.
#' @param pcv1 periodicity of CV1.
#' @param pcv2 periodicity of CV2.
#' @param pcv3 periodicity of CV3.
#' @param ignoretime time in the first column of the HILLS file will be ignored.
#' @param times the index of a column containing time values.
#' @param centres vector of indexes of columns containing centres of hills.
#' @param widths vector of indexes of columns containing widths of hills.
#' @param heights the index of a column containing heights of hills.
#' @param sep the field separator character.
#' @param dec the character used in the file for decimal points.
#' @param na.strings a character vector of strings which are to be interpreted as NA values.
#' @param comment.char a character vector of length one containing a single character or an empty string.
#' @return hillsfile object.
#'
#' @export
#' @examples
#' l1<-"1 -1.587 -2.969 3.013 0.3 0.3 0.3 1.111 10"
#' l2<-"2 -1.067 2.745 2.944 0.3 0.3 0.3 1.109 10"
#' l3<-"3 -1.376 2.697 3.049 0.3 0.3 0.3 1.080 10"
#' l4<-"4 -1.663 2.922 -3.065 0.3 0.3 0.3 1.072 10"
#' fourhills<-c(l1,l2,l3,l4)
#' tf <- tempfile()
#' writeLines(fourhills, tf)
#' read.hills3d2(tf, per=c(TRUE,TRUE), times=1, centres=2:4, widths=5:7, heights=8)
read.hills3d2<-function(file="HILLS", per=c(FALSE, FALSE, FALSE), pcv1=c(-pi,pi), pcv2=c(-pi,pi), pcv3=c(-pi,pi),
ignoretime=FALSE, times=1, centres=2:4, widths=5:7, heights=8,
sep="", dec=".", na.strings="NA", comment.char="#") {
hillsf<-read.table(file, header=F, sep=sep, dec=dec, na.strings=na.strings, comment.char=comment.char)
if(length(times)!=1) {
stop("Error: Parameter values must be set to one value")
}
if(length(centres)!=3) {
stop("Error: Parameter centres must be set to three values")
}
if(length(widths)!=3) {
stop("Error: Parameter widths must be set to three values")
}
if(length(heights)!=1) {
stop("Error: Parameter heights must be set to one value")
}
cat("3D HILLS file read\n")
newhillsf <- data.frame(rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)),
rep(0, times=nrow(hillsf)), rep(0, times=nrow(hillsf)))
newhillsf[,1]<-hillsf[,times]
newhillsf[,2]<-hillsf[,centres[1]]
newhillsf[,3]<-hillsf[,centres[2]]
newhillsf[,4]<-hillsf[,centres[3]]
newhillsf[,5]<-hillsf[,widths[1]]
newhillsf[,6]<-hillsf[,widths[2]]
newhillsf[,7]<-hillsf[,widths[3]]
newhillsf[,8]<-hillsf[,heights]
if(ignoretime) {
warning("Warning: The time will be updated automatically from zero according to the first step!\n")
newhillsf[,1]<-seq(from=hillsf[1,times], by=hillsf[1,times], length.out=nrow(hillsf))
}
hills<-list(hillsfile=newhillsf, time=newhillsf[,1], cv1=newhillsf[,2], cv2=newhillsf[,3], cv3=newhillsf[,4],
size=dim(hillsf), filename=file, per=per, pcv1=pcv1, pcv2=pcv2, pcv3=pcv3)
class(hills) <- "hillsfile3d"
return(hills)
}
#' Print hillsfile3d
#'
#' `print.hillsfile3d` prints dimensionality and size of a hillsfile object.
#'
#' @param x hillsfile3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' acealanme3d
print.hillsfile3d<-function(x,...) {
hills <- x
toprint <- "3D hills file "
toprint <- paste(toprint, toString(hills$filename), sep="")
toprint <- paste(toprint, " with ", sep="")
toprint <- paste(toprint, toString(hills$size[1]), sep="")
toprint <- paste(toprint, " lines", sep="")
message(toprint)
}
#' Print summary for hillsfile3d
#'
#' `summary.hillsfile3d` prints dimensionality, size and collective variable ranges of a hillsfile3d object.
#'
#' @param object hillsfile3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' summary(acealanme3d)
summary.hillsfile3d<-function(object,...) {
hills <- object
toprint <- "3D hills file "
toprint <- paste(toprint, toString(hills$filename), sep="")
toprint <- paste(toprint, " with ", sep="")
toprint <- paste(toprint, toString(hills$size[1]), sep="")
toprint <- paste(toprint, " lines", sep="")
message(toprint)
toprint <- "The CV1 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,2])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,2])), sep="")
message(toprint)
toprint <- "The CV2 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,3])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,3])), sep="")
message(toprint)
toprint <- "The CV3 ranges from "
toprint <- paste(toprint, toString(min(hills$hillsfile[,4])), sep="")
toprint <- paste(toprint, " to ", sep="")
toprint <- paste(toprint, toString(max(hills$hillsfile[,4])), sep="")
message(toprint)
}
#' Print first n lines of hillsfile3d
#'
#' `head.hillsfile3d` prints first n lines of a hillsfile3d object.
#'
#' @param x hillsfile3d object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' head(acealanme3d)
head.hillsfile3d<-function(x, n=10,...) {
return(head(x$hillsfile, n=n))
}
#' Print last n lines of hillsfile3d
#'
#' `tail.hillsfile3d` prints last n lines of a hillsfile3d object.
#'
#' @param x hillsfile3d object.
#' @param n number of lines (default 10).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tail(acealanme3d)
tail.hillsfile3d<-function(x, n=10,...) {
return(tail(x$hillsfile, n=n))
}
#' @export
`+.hillsfile3d`<-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")
}
}
if(ncol(hills1$hillsfile)==9 || ncol(hills1$hillsfile)==10) {
if(hills1$per[2]!=hills2$per[2]) {
stop("Error: You can sum only hills of same periodicity")
}
if(hills1$per[3]!=hills2$per[3]) {
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, pcv3=hills1$pcv3)
class(hills) <- "hillsfile3d"
return(hills)
}
#' Plot hillsfile3d object
#'
#' `plot.hillsfile3d` plots hillsfile object. It plots CV1 vs CV2 vs CV3.
#'
#' @param x hillsfile object.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param zlab a title for the z axis: see 'title'.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param col color code or name, 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 zlim numeric vector of length 2, giving the z coordinates range.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' plot(acealanme3d)
plot.hillsfile3d<-function(x,
xlab="CV1", ylab="CV2", zlab="CV3",
main=NULL, sub=NULL,
col="orange",...) {
hills <-x
points3d(hills$hillsfile[,2:4], col=col)
axes3d()
title3d(xlab=xlab, ylab=ylab, zlab=zlab,
main=main, sub=sub)
box3d()
}
#' Plot evolution of heights of hills in hillsfile3d object
#'
#' `plotheights.hillsfile3d` 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(acealanme3d)
plotheights.hillsfile3d<-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(ignoretime) {
plot(seq(from=hills$hillsfile[1,1],by=hills$hillsfile[1,1],length.out=nrow(hills$hillsfile)),
hills$hillsfile[,8], 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[,8], type="l",
xlab=xlab, ylab=ylab,
main=main, sub=sub,
col=col, lwd=lwd,
asp=asp, axes=axes)
}
}
#' Calculate 3D free energy surface by Bias Sum algorithm
#'
#' `fes.hillsfile3d` sums up hills using fast Bias Sum algorithm.
#'
#' @param hills hillsfile3d 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 x lim 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(acealanme3d, imax=5000)
fes.hillsfile3d<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=NULL) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
warning("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(max(hills$hillsfile[,5])/min(hills$hillsfile[,5])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,6])/min(hills$hillsfile[,6])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(max(hills$hillsfile[,7])/min(hills$hillsfile[,7])>1.00000000001) {
stop("Error: Bias Sum algorithm works only with hills of the same sizes")
}
if(is.null(npoints)) {
npoints <- 64
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
minCV3 <- min(hills$hillsfile[,4])
maxCV3 <- max(hills$hillsfile[,4])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
zlims<-c(minCV3-0.05*(maxCV3-minCV3), maxCV3+0.05*(maxCV3-minCV3))
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(!is.null(zlim)) {zlims<-zlim}
if((hills$per[3]==T)&is.null(zlim)) {zlims<-hills$pcv3}
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")
}
}
if(hills$per[3]==T) {
if(min(hills$hillsfile[,4])<zlims[1]) {
stop("Error: The third collective variable outside pcv3")
}
if(max(hills$hillsfile[,4])>zlims[2]) {
stop("Error: The third collective variable outside pcv3")
}
}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
z<-0:(npoints-1)*(zlims[2]-zlims[1])/(npoints-1)+zlims[1]
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d1p1(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d1p2(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d1p3(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d1p12(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d1p13(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d1p23(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d1p123(npoints*(hills$hillsfile[,2]-xlims[1])/(xlims[2]-xlims[1]),
npoints*(hills$hillsfile[,3]-ylims[1])/(ylims[2]-ylims[1]),
npoints*(hills$hillsfile[,4]-zlims[1])/(zlims[2]-zlims[1]),
npoints*max(hills$hillsfile[,5])/(xlims[2]-xlims[1]),
npoints*max(hills$hillsfile[,6])/(ylims[2]-ylims[1]),
npoints*max(hills$hillsfile[,7])/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
fesm<-aperm(array(fesm, c(npoints, npoints, npoints)))
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=3, per=hills$per,
x=x, y=y, z=z, pcv1=hills$pcv1, pcv2=hills$pcv2, pcv3=hills$pcv3)
class(cfes) <- "fes3d"
return(cfes)
}
#' Calculate 3D free energy surface by conventional algorithm
#'
#' `fes2.hills3d` sums up hills using slow conventional algorithm. It can be used
#' as a reference or when hill widths are variable.
#'
#' @param hills hillsfile3d 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(acealanme3d, imax=100)
fes2.hillsfile3d<-function(hills, imin=1, imax=NULL, xlim=NULL, ylim=NULL, zlim=NULL, npoints=NULL) {
if(!is.null(imax)) {
if(hills$size[1]<imax) {
warning("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(is.null(npoints)) {
npoints <- 64
}
minCV1 <- min(hills$hillsfile[,2])
maxCV1 <- max(hills$hillsfile[,2])
minCV2 <- min(hills$hillsfile[,3])
maxCV2 <- max(hills$hillsfile[,3])
minCV3 <- min(hills$hillsfile[,4])
maxCV3 <- max(hills$hillsfile[,4])
xlims<-c(minCV1-0.05*(maxCV1-minCV1), maxCV1+0.05*(maxCV1-minCV1))
ylims<-c(minCV2-0.05*(maxCV2-minCV2), maxCV2+0.05*(maxCV2-minCV2))
zlims<-c(minCV3-0.05*(maxCV3-minCV3), maxCV3+0.05*(maxCV3-minCV3))
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(!is.null(zlim)) {zlims<-zlim}
if((hills$per[3]==T)&is.null(zlim)) {zlims<-hills$pcv3}
x<-0:(npoints-1)*(xlims[2]-xlims[1])/(npoints-1)+xlims[1]
y<-0:(npoints-1)*(ylims[2]-ylims[1])/(npoints-1)+ylims[1]
z<-0:(npoints-1)*(zlims[2]-zlims[1])/(npoints-1)+zlims[1]
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d2((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==F)) {
fesm<-hills3d2p1((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d2p2((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d2p3((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==F)) {
fesm<-hills3d2p12((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==F)&(hills$per[3]==T)) {
fesm<-hills3d2p13((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==F)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d2p23((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
if((hills$per[1]==T)&(hills$per[2]==T)&(hills$per[3]==T)) {
fesm<-hills3d2p123((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]-zlims[1])/(zlims[2]-zlims[1]),
(npoints-1)*hills$hillsfile[,5]/(xlims[2]-xlims[1]),
(npoints-1)*hills$hillsfile[,6]/(ylims[2]-ylims[1]),
(npoints-1)*hills$hillsfile[,7]/(zlims[2]-zlims[1]),
hills$hillsfile[,8],npoints,imin-1,imax-1)
}
fesm<-aperm(array(fesm, c(npoints, npoints, npoints)))
cfes<-list(fes=fesm, hills=hills$hillsfile, rows=npoints, dimension=3, per=hills$per,
x=x, y=y, z=z, pcv1=hills$pcv1, pcv2=hills$pcv2, pcv3=hills$pcv3)
class(cfes) <- "fes3d"
return(cfes)
}
#' Read 3D free energy surface from PLUMED sum_hills
#'
#' `read.plumed3d` reads 3D free energy surface from PLUMED sum_hills.
#' The grid in the input file must contain the same number of points
#' for CV1, CV2 and CV3. It does not use the header of the file.
#' Periodicity must be specified.
#'
#' @param file input file from PLUMED sum_hills.
#' @param per logical vector specifying periodicity of collective variables.
#' @return fes object.
#'
#' @export
#' @examples
#' l1<-" -3.14 -3.14 -3.14 -61.13 -47.43 19.00 2.04"
#' l2<-" -1.05 -3.14 -3.14 -70.72 25.95 25.78 2.43"
#' l3<-" 1.05 -3.14 -3.14 -65.58 8.34 2.82 -3.09"
#' l4<-" -3.14 -1.05 -3.14 -51.31 -43.88 -19.91 1.51"
#' l5<-" -1.05 -1.05 -3.14 -66.43 7.67 -22.45 -0.39"
#' l6<-" 1.05 -1.05 -3.14 -61.08 -7.50 -7.36 -0.83"
#' l7<-" -3.14 1.05 -3.14 -53.07 -55.12 0.19 -0.28"
#' l8<-" -1.05 1.05 -3.14 -62.81 36.19 1.65 0.45"
#' l9<-" 1.05 1.05 -3.14 -65.28 22.84 11.47 0.59"
#' l10<-" -3.14 -3.14 -1.05 -13.03 -32.17 8.24 -35.25"
#' l11<-" -1.05 -3.14 -1.05 -21.88 17.89 21.91 -51.20"
#' l12<-" 1.05 -3.14 -1.05 -14.49 3.60 6.04 -44.05"
#' l13<-" -3.14 -1.05 -1.05 -2.26 -7.00 -7.01 -10.65"
#' l14<-" -1.05 -1.05 -1.05 -8.21 3.69 -22.89 -28.48"
#' l15<-" 1.05 -1.05 -1.05 -1.10 0.52 3.59 -1.99"
#' l16<-" -3.14 1.05 -1.05 -3.75 -11.70 -5.65 -15.36"
#' l17<-" -1.05 1.05 -1.05 -1.15 5.75 1.05 -2.42"
#' l18<-" 1.05 1.05 -1.05 -10.67 8.23 -10.42 -36.77"
#' l19<-" -3.14 -3.14 1.05 -4.64 -13.79 10.51 14.96"
#' l20<-" -1.05 -3.14 1.05 -7.80 12.24 20.59 23.03"
#' l21<-" 1.05 -3.14 1.05 -5.32 3.46 3.17 21.99"
#' l22<-" -3.14 -1.05 1.05 -2.06 -6.59 0.17 10.04"
#' l23<-" -1.05 -1.05 1.05 -9.69 8.43 -0.97 36.97"
#' l24<-" 1.05 -1.05 1.05 -0.19 -0.44 -0.26 0.91"
#' l25<-" -3.14 1.05 1.05 -7.98 -23.02 3.97 26.98"
#' l26<-" -1.05 1.05 1.05 -4.64 13.66 -9.74 10.15"
#' l27<-" 1.05 1.05 1.05 -13.42 15.78 16.36 41.60"
#' twentysevenpoints<-c(l1,l2,l3,l4,l5,l6,l7,l8,l9,l10,
#' l11,l12,l13,l14,l15,l16,l17,l18,l19,l20,
#' l21,l22,l23,l24,l25,l26,l27)
#' tf <- tempfile()
#' writeLines(twentysevenpoints, tf)
#' read.plumed3d(tf, per=c(TRUE,TRUE,TRUE))
read.plumed3d<-function(file="fes.dat", per=c(FALSE,FALSE,FALSE)) {
hillsf<-read.table(file, header=F, comment.char="#")
bins<-round(nrow(hillsf)^(1/3))
if(bins^3!=nrow(hillsf)) {
stop("Error: the number of bins cannot be determined.")
}
x <- hillsf[1:bins,1]
y <- hillsf[(0:(bins-1))*bins+1,2]
z <- hillsf[(0:(bins-1))*bins*bins+1,3]
fesm <- array(hillsf[,4], c(bins, bins, bins))
cfes <- list(fes=fesm, hills=NULL, rows=bins, dimension=3, per=per, x=x, y=y, z=z,
pcv1=c(min(x), max(x)), pcv2=c(min(y), max(y)), pcv3=c(min(z), max(z)))
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`+.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
if(fes1$rows!=fes2$rows) {
stop("Error: Free energy surfaces have different numbers of points, exiting")
}
if(sum(fes1$x!=fes2$x)>0) {
stop("Error: Free energy surfaces have different CV1 axes, exiting")
}
if(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
if(sum(fes1$z!=fes2$z)>0) {
stop("Error: Free energy surfaces have different CV3 axes, exiting")
}
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, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes+fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
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, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
} else if(inherits(fes1, "fes")) {
stop("Error: cannot sum 3D free energy surfaces with 1D or 2D, exiting")
} else if(inherits(fes2, "fes")) {
stop("Error: cannot sum 3D free energy surfaces with 1D or 2D, exiting")
}
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`-.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
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(sum(fes1$y!=fes2$y)>0) {
stop("Error: Free energy surfaces have different CV2 axes, exiting")
}
if(sum(fes1$z!=fes2$z)>0) {
stop("Error: Free energy surfaces have different CV3 axes, exiting")
}
warning("Warning: FES obtained by subtraction of two FESes will inherit hills only from the first FES\n")
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, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes-fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
cfes<-list(fes=fes1-fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
} else if(inherits(fes1, "fes")) {
stop("Error: cannot subtract 3D free energy surfaces with 1D or 2D, exiting")
} else if(inherits(fes2, "fes")) {
stop("Error: cannot subtract 3D free energy surfaces with 1D or 2D, exiting")
}
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`*.fes3d`<-function(fes1, fes2) {
if(inherits(fes1, "fes3d")&inherits(fes2, "fes3d")) {
stop("Error: You cannot multiply fes by fes")
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes*fes2, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(fes2, "fes3d")) {
cfes<-list(fes=fes1*fes2$fes, hills=fes2$hills, rows=fes2$rows, dimension=fes2$dimension, per=fes2$per, x=fes2$x, y=fes2$y, z=fes2$z, pcv1=fes2$pcv1, pcv2=fes2$pcv2, pcv3=fes2$pcv3)
}
warning("Warning: multiplication of FES will multiply the FES but not hill heights\n")
class(cfes) <- "fes3d"
return(cfes)
}
#' @export
`/.fes3d`<-function(fes1, coef) {
if(inherits(fes1, "fes3d")&inherits(coef, "fes3d")) {
stop("Error: You cannot divide fes by fes")
} else if(inherits(fes1, "fes3d")) {
cfes<-list(fes=fes1$fes/coef, hills=fes1$hills, rows=fes1$rows, dimension=fes1$dimension, per=fes1$per, x=fes1$x, y=fes1$y, z=fes1$z, pcv1=fes1$pcv1, pcv2=fes1$pcv2, pcv3=fes1$pcv3)
} else if(inherits(coef, "fes3d")) {
stop("Error: You cannot divide something by fes")
}
warning("Warning: division of FES will divide the FES but not hill heights\n")
class(cfes) <- "fes3d"
return(cfes)
}
#' Calculate minimum of 3D free energy surface
#'
#' `min.fes3d` calculates minimum of free energy in a fes3d object.
#'
#' @param inputfes fes3d 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(acealanme3d, imax=5000)
#' min(tfes)
min.fes3d<-function(inputfes, na.rm=NULL,...) {
return(min(inputfes$fes, na.rm=na.rm))
}
#' Calculate maximum of 3D free energy surface
#'
#' `max.fes3d` calculates maximum of free energy in a fes3d object.
#'
#' @param inputfes fes3d 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(acealanme3d, imax=5000)
#' max(tfes)
max.fes3d<-function(inputfes, na.rm=NULL,...) {
return(max(inputfes$fes, na.rm=na.rm))
}
#' Print minimum and maximum of 3D free energy surface
#'
#' `print.fes3d` prints dimensionality, minimum and maximum of
#' free energy in a fes object
#'
#' @param x fes3d object
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme3d, imax=5000)
#' tfes
print.fes3d<-function(x,...) {
inputfes<-x
toprint <- "3D free energy surface with "
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " points, maximum ", sep="")
toprint <- paste(toprint, toString(max(inputfes$fes)), sep="")
toprint <- paste(toprint, " and minimum ", sep="")
toprint <- paste(toprint, toString(min(inputfes$fes)), sep="")
message(toprint)
}
#' Print summary of 3D free energy surface
#'
#' `summary.fes3d` prints minimum and maximum of free energy
#' in a fes object.
#'
#' @param object fes3d object.
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes<-fes(acealanme3d, imax=5000)
#' summary(tfes)
summary.fes3d<-function(object,...) {
inputfes <- object
toprint <- "3D free energy surface with "
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " x ", sep="")
toprint <- paste(toprint, toString(inputfes$rows), sep="")
toprint <- paste(toprint, " points, maximum ", sep="")
toprint <- paste(toprint, toString(max(inputfes$fes)), sep="")
toprint <- paste(toprint, " and minimum ", sep="")
toprint <- paste(toprint, toString(min(inputfes$fes)), sep="")
message(toprint)
}
#' Plot 3D free energy surface object
#'
#' `plot.fes3d` plots 3D free energy surface using .
#'
#' @param x fes3d object.
#' @param xlab a title for the x axis: see 'title'.
#' @param ylab a title for the y axis: see 'title'.
#' @param zlab a title for the z axis: see 'title'.
#' @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 col color of the free energy surface. It can be a single color
#' or a vector with multiple colors for multiple 3D isosurfaces.
#' @param alpha number or numeric vector of alpha levels (transparency) of
#' 3D isosurfaces.
#' @param main an overall title for the plot: see 'title'.
#' @param sub a sub title for the plot: see 'title'.
#' @param level number or numeric vector of levels at which to draw 3D isosurface.
#' @param fill a logical value indicating whether 3D isosurface is ploted as
#' solid surface (True) or wireframe (False).
#' @param ... further arguments passed to or from other methods.
#'
#' @export
#' @examples
#' tfes3d<-fes(acealanme3d, imax=5000)
#' plot(tfes3d)
plot.fes3d<-function(x, xlab=NULL, ylab=NULL, zlab=NULL,
xlim=NULL, ylim=NULL, zlim=NULL,
level=NULL, col=NULL, alpha=NULL,
main=NULL, sub=NULL,
fill=TRUE,...) {
close.screen(all.screens=TRUE)
inputfes<-x
fes<-inputfes$fes
rows<-inputfes$rows
x<-inputfes$x
y<-inputfes$y
z<-inputfes$z
if(is.null(xlim)) {
xlim=range(x)
} else {
fes <- fes[x >= xlim[1] & x <= xlim[2],,]
x <- x[x >= xlim[1] & x <= xlim[2]]
}
if(is.null(ylim)) {
ylim=range(y)
} else {
fes <- fes[,y >= ylim[1] & y <= ylim[2],]
y <- y[y >= ylim[1] & y <= ylim[2]]
}
if(is.null(zlim)) {
zlim=range(z)
} else {
fes <- fes[,,z >= zlim[1] & z <= zlim[2]]
z <- z[z >= zlim[1] & z <= zlim[2]]
}
if(is.null(xlab)) xlab="CV1"
if(is.null(ylab)) ylab="CV2"
if(is.null(zlab)) zlab="CV3"
if(is.null(level)) level=(max(fes)+min(fes))/2
if(length(level)>1) {
if(is.null(col)) col<-rainbow(1.35*length(level))[length(level):1]
if(is.null(alpha)) {
alpha<-length(level):1/length(level)
level<-sort(level)
}
} else {
if(is.null(col)) col<-"orange"
if(is.null(alpha)) alpha<-1
}
rgl::plot3d(c(), c(), c(), xlim=xlim, ylim=ylim, zlim=zlim,
xlab=xlab, ylab=ylab, zlab=zlab, main=main, sub=sub)
contour3d(f=fes, level=level, x=x, y=y, z=z,
color=col, alpha=alpha, fill=fill, add=T)
rgl::aspect3d(1,1,1)
axes3d()
box3d()
}
Try the metadynminer3d 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.