Learn R Programming

MetadynMiner

Web site

https://metadynamics.cz/metadynminer/

https://spiwokv.github.io/metadynminer/ (pkgdown web)

Introduction

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.

Usage

# Install from R repository
install.packages("metadynminer")

# Install from GitHub by devtools
install.packages("devtools")
devtools::install_github("spiwokv/metadynminer")

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

# Sum two hills files
hillsf+hillsf

# Summary of a hills file
summary(hillsf)

# Plot CVs
plot(hillsf)
# Plot heights
plotheights(hillsf)
# Calculate FES by bias sum (alternatively use fes2 for conventional calculation)
tfes<-fes(hillsf)

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

# Sum two FESes
tfes+tfes

# Calculate and subtract min, max or mean from a FES
tfes<-tfes-min(tfes)

# Summary of FES
summary(tfes)

# Plot FES
plot(tfes)
# Plot FES with color scale
plot(tfes, colscale=T)
# Find minima
minima<-fesminima(tfes)

# Summary of minima
summary(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
plot(minima)
# Calculate free energy profile for minima
prof<-feprof(minima)

# Plot free energy profile for minima
plot(prof)

# Make 1D free energy surface from the 2D one
tfes1<-fes2d21d(hillsf, remdim=2) # T=300K, kJ/mol
plot(tfes1)
# Calculate transition path using Nudged Elastic Band
myneb <- neb(minima, min1="A", min2="B")
myneb
summary(myneb)

# Plot transition path
plot(myneb)
# Plot transition path on FES
plot(minima)
linesonfes(myneb)

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))
tfes<-fes(hillsf)
png("filename.png", height=8, width=8, units='cm', res=600, pointsize=6)
plot(tfes)
dev.off()

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 and concatenated by:

#install.packages("magick") <- install before the first run
library(metadynminer)
library(magick)
odir<-file.path(tempdir(), "frames")
dir.create(odir, recursive=T)
hillsf <- read.hills("HILLS", per=c(T,T))
tfes<-fes(hillsf, imax=300)
png(paste(odir, "/snap%04d.png", sep=""))
plot(tfes, zlim=c(-200,0))
for(i in 1:99) {
  tfes<-tfes+fes(hillsf, imin=300*i+1, imax=300*(i+1))
  plot(tfes, zlim=c(-200,0))
}
dev.off()
figs <- list.files(odir, full.names=T)
rfigs <- lapply(figs, image_read)
allfigs <- image_join(rfigs)
anim <- image_animate(allfigs, fps=25)
image_write(image=anim, path="anim.gif")

Package magick was used. Alternatively, files can be concatenated by outside R by any movie making program.

If you instead want to see flooding, type:

#install.packages("magick") <- install before the first run
library(metadynminer)
library(magick)
odir<-file.path(tempdir(), "frames")
dir.create(odir, recursive=T)
hillsf <- read.hills("HILLS", per=c(T,T))
tfes<-fes(hillsf)
png(paste(odir, "/snap%04d.png", sep=""))
plot(tfes, zlim=c(-200,0))
for(i in 0:99) {
  tfes<-tfes + -1*fes(hillsf, imin=300*i+1, imax=300*(i+1))
  plot(tfes, zlim=c(-200,0))
}
dev.off()
figs <- list.files(odir, full.names=T)
rfigs <- lapply(figs, image_read)
allfigs <- image_join(rfigs)
anim <- image_animate(allfigs, fps=25)
image_write(image=anim, path="flooding.gif")

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)
box()

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:

acealanmec<-acealanme
acealanmec$pcv1<-c(0,2*pi)
acealanmec$hillsfile[acealanmec$hillsfile[,2]<0,2]<-
    acealanmec$hillsfile[acealanmec$hillsfile[,2]<0,2]+2*pi
tfes<-fes(acealanmec)
plot(tfes)

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.

Training

MetadynMiner Webinar video

Plumed Masterclass tutorial

Plumed Masterclass video 1

Plumed Masterclass video 2

Contact

Vojtech Spiwok - spiwokv{youknowwhat}vscht.cz

To contribute, se CONTRIBUTING.md

Copy Link

Version

Install

install.packages('metadynminer')

Monthly Downloads

246

Version

0.1.7

License

GPL-3

Maintainer

Last Published

April 14th, 2022

Functions in metadynminer (0.1.7)

plot.minima

Plot minima object
lines.nebpath

Plot lines for Nudged Elastic Band
plot.profiles

Plot free energy profile
points.hillsfile

Plot points for hillsfile object
lines.hillsfile

Plot lines for hillsfile object
points.nebpath

Plot points for Nudged Elastic Band
plotheights

Plot evolution of heights of hills (generic function for 'metadynminer' and 'metadynminer3d')
fes2d21d

Calculate 1D free energy surface from hillsfile object
neb

Find transition path on free energy surface by Nudged Elastic Band method
fespoint.hillsfile

Calculate free energy at given point in the CV space
print.profiles

Print profiles object
fespoint

Calculate free energy at given point in the CV space (generic function for 'metadynminer' and 'metadynminer3d')
min.fes

Calculate minimum of free energy surface
print.nebpath

Print Nudged Elastic Band minima
print.hillsfile

Print hillsfile
plot.nebpath

Plot Nudged Elastic Band
print.minima

Print minima object
head.hillsfile

Print first n lines of hillsfile
linesonfes

Plot lines for Nudged Elastic Band projected onto free energy surface
fesminima

Find free energy minima in the fes object (generic function for 'metadynminer' and 'metadynminer3d')
fesminima.fes

Find free energy minima in the fes object
oneminimum

Creates one ad hoc free energy minimum for a fes object (generic function for 'metadynminer' and 'metadynminer3d')
pointsonfes

Plot points for Nudged Elastic Band projected onto free energy surface
oneminimum.fes

Creates one ad hoc free energy minimum for a fes object
print.fes

Print dimensionality, minimum and maximum of free energy surface
prob

Calculate probability of free energy surface
summary.hillsfile

Print summary for hillsfile
summary.minima

Print minima object summary
plot.fes

Plot free energy surface object
max.fes

Calculate maximum of free energy surface
read.hills

Read HILLS from Plumed
lines.fes

Plots 1D free energy surface object as lines
plot.hillsfile

Plot hillsfile object
tail.hillsfile

Print last n lines of hillsfile
plotheights.hillsfile

Plot evolution of heights of hills in hillsfile object
summary.fes

Print summary of free energy surface
read.plumed

Read 1D or 2D free energy surface from PLUMED sum_hills
points.fes

Plots 1D free energy surface object as points
summary.nebpath

Print summary for Nudged Elastic Band
summary.profiles

Print summary for free energy profile
fes

Calculate free energy surface by Bias Sum algorithm (generic function for 'metadynminer' and 'metadynminer3d')
feprof.minima

Calculate free energy profile for minima object
emptyhills

Generate empty HILLS from Plumed
feprof

Calculate free energy profile for minima object (generic function for 'metadynminer' and 'metadynminer3d')
acealanme

Hills from 30 ns metadynamics of AceAlaNme in water with two collective variables
fes2.hillsfile

Calculate free energy surface by conventional algorithm
acealanme1d

Hills from 30 ns metadynamics of AceAlaNme in water with one collective variable
fes2

Calculate free energy surface by conventional algorithm (generic function for 'metadynminer' and 'metadynminer3d')
fes.hillsfile

Calculate free energy surface by Bias Sum algorithm