https://metadynamics.cz/metadynminer3d/
https://spiwokv.github.io/metadynminer3d/ (pkgdown web)
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.
MetadynMiner3d is its addendum for plotting 3D free energy surfaces. It uses RGL package. MetadynMiner3d installs and loads metadynminer automatically.
# Install from R repository
install.packages("metadynminer3d")
# Install from GitHub by devtools
install.packages("devtools")
devtools::install_github("spiwokv/metadynminer3d")
# Load library
library(metadynminer3d)
# Read hills file
hillsf<-read.hills3d("HILLS", per=c(TRUE, TRUE, TRUE)) # HILLS with periodicity on CV1, CV2 and CV3
# 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, level=20)
# Find minima
minima<-fesminima(tfes)
# Summary of minima
summary(minima)
# Plot free energy minima
plot(minima)
# Calculate free energy profile for minima
prof<-feprof(minima)
# Plot free energy profile for minima
plot(prof)
Following script can be used to generate a publication quality figure:
hillsf<-read.hills3d("HILLS", per=c(TRUE, TRUE, TRUE))
tfes<-fes(hillsf)
plot(tfes)
Change window size, zoom in and rotate if necessary, then (without closing the window) type:
rgl.snapshot(filename="plot.png")
You can save free energy surface in WebGL and present it on a web site by typing:
writeWebGL(filename="fes3d.html")
It is possible to make movie of rotation of the plot by:
movie3d(spin3d(axis=c(0,0,1)), duration=3, dir=".", type="gif")
This will create a three-second animated GIF movie.gif
file of rotation around z-axis.
The aspect ratio of the box can be modified by:
aspect3d(2, 1, 1)
after plot
command. This makes the x-axis twice longer than other axes.
Vojtech Spiwok - spiwokv{youknowwhat}vscht.cz
To contribute, se CONTRIBUTING.md