Build Status

Build Status

CRAN downloads



MetadynMiner is R packages for reading, analysis and visualization of metadynamics HILLS files produced by Plumed. It reads HILLS files from Plumed, calculates free energy surface by fast Bias Sum algorithm, finds minima and analyses transition paths by Nudged Elastic Band method.


# Install from R repository

# Install from GitHub by devtools

# Load library
library(metadynminer) # in future will be added to R repository

# Read hills file
hillsf<-read.hills("HILLS", per=c(TRUE, TRUE)) # HILLS with periodicity on CV1 and CV2

# Sum two hills files

# Summary of a hills file

# Plot CVs

# Plot heights

# Calculate FES by bias sum (alternatively use fes2 for conventional calculation)

# Calculate FES for given range (indexes of hills)
tfes<-fes(hillsf, imin=5000, imax=10000)

# Sum two FESes

# Calculate and subtract min, max or mean from a FES

# Summary of FES

# Plot FES

# Plot FES with color scale
plot(tfes, colscale=T)

# Find minima

# Summary of minima

# Create empty minima list, create ad hoc minimum, add minima
minima<-oneminimum(tfes, cv1=0, cv2=0)
fesminima(fes) + oneminimum(tfes, cv1=0, cv2=0)

# Plot free energy minima

# Calculate free energy profile for minima

# Plot free energy profile for minima

# Make 1D free energy surface from the 2D one
tfes1<-fes2d21d(hillsf, remdim=2) # T=300K, kJ/mol

# Calculate transition path using Nudged Elastic Band
myneb <- neb(minima, min1="A", min2="B")

# Plot transition path

# Plot transition path on FES

Tips and Tricks

Publication quality figures

Following script can be used to generate a publication quality figure (8x8 cm, 600 dpi):

hillsf <- read.hills("HILLS", per=c(T,T))
png("filename.png", height=8, width=8, units='cm', res=600, pointsize=6)

Making FES relative to the global minimum

You can set free energy minimum to zero by typing:

tfes <- tfes - min(tfes)

Hills from restarted simulations

There are two ways to cope with hills file from restarted simulations, i.e. simulations where the time column starts from zero at every restart. It is possible to set ignoretime=T in the read.hills function. It will take the time of the first hill and use it as the uniform step.

Alternatively, it is possible to use ignoretime=T for other ploting functions.

MetadynMiner does not support hills files with variable PACE.

Making movie

Individual snapshots of a movie can be generated by:

hillsf <- read.hills("HILLS", per=c(T,T))
tfes<-fes(hillsf, tmax=100)
plot(tfes, zlim=c(-200,0))
for(i in 1:299) {
 tfes<-tfes+fes(hillsf, imin=100*i+1, imax=100*(i+1))
 plot(tfes, zlim=c(-200,0))

These files can be concatenated by a movie making program such as mencoder.

If you instead want to see flooding, type:

hillsf <- read.hills("HILLS", per=c(T,T))
plot(tfes, zlim=c(-200,0))
for(i in 0:299) {
  tfes<-tfes + -1*fes(hillsf, imin=100*i+1, imax=100*(i+1))
  plot(tfes, zlim=c(-200,0))

Evaluation of convergence of one CV

You can use function fes2d21d to convert a 2D surface to 1D and to evaluate the evolution:

hillsf <- read.hills("HILLS", per=c(T,T))
tfes1<-fes2d21d(hillsf, remdim=2)
plot(tfes1-min(tfes1), ylim=c(0,80), lwd=4, col="black")
for(i in 1:10) {
 tfes1<-fes2d21d(hillsf, imax=3000*i)
 lines(tfes1-min(tfes1), col=rainbow(13)[i])

Transforming CVs

If you want to use degrees instead of radians on axes, set axes=F in the plot function and then plot (without closing the plot window!) both axes separately.

plot(tfes, axes=F)
axis(2, at=-3:3*pi/3, labels=-3:3*60)
axis(1, at=-3:3*pi/3, labels=-3:3*60)

The expression -3:3 will generate a vector {-3,-2,-1,0,1,2,3}, which can be multiplied by pi/3 (tick positions in radians) or by 60 (tick positions in degrees). If you want to transform just one axis, e.g. the horizontal one while keeping the vertical unchanged, simply type axis(2) for the vertical one. box() redraws a box.

kcal vs kJ

MetadynMiner works in kJ/mol by defauls. If your MD engine uses kcal/mol instead, you can either multiply your free energy surface by 4.184 to get kJ/mol. If you prefer to keep kcal/mol, you can set eunit="kcal/mol" for functions fes2d21d or summary of minima object. Other units are not supported.

Shifting a periodic CV

It may happen that some simulations with a torsion CV it may be difficult to analyze and visualize it in the range -pi - +pi. However, this problem is not very common so we did not make any user friendly way how to solve this and it can be solved in a user unfriendly way. Let us consider we want to shift the first collective variable to be in the range 0 - 2pi. First we will make copy of acealanme. We will change its pcv1 to c(0,2*pi). Finally we can add 2pi to the first collective variable:


The hills file object has several instances including hillsfile, which contains the HILLS file, and pcv1 with collective variable periodicity. They can be printed by $ operator. The expression acealanmec$hillsfile[,2] prints all values of the first collective variable. The expression acealanmec$hillsfile[,2]<0 prints the same number of TRUE or FALSE values depending whether the first collective variable is positive or negative. The expression acealanmec$hillsfile[acealanmec$hillsfile[,2]<0,2] prints only negative values of the first collective variable. They can be replaced by the same value + 2pi.