Source code for metafalcon.analysis.reconstruct
#!/usr/bin/env python
"""Reconstruct and plot the metadynamics potential."""
from __future__ import print_function
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from .. import meta
from .. import constants as c
[docs]class reconstruct_metadynamics(meta.metadynamics):
[docs] def reconstruct_vg(self, x0=None, ngauss=None, plot=True, animate=False, stride=1, unit="eV",
fes=True):
"""
Reconstruct and plot potential from data in vg_centers.dat and meta-config.py. Grid and
VG are saved as numpy arrays to vg_grid.npy and vg_values.npy
Parameters
----------
x0 : list of numpy arrays
manually specified x-values for each collective variable (default: None)
ngauss : number
first ngauss gaussian functions of the simulation will be evaluated (default: all)
plot : boolean
plot single image of vg with the specified number of gaussians (default: True)
animate : boolean
animate evolution of vg with number of gaussians (default: False)
stride : number
number of gaussians to be added in each frame of the animation (default: 1)
fes : boolean
plot the free energy surface instead of the metadynamics potential (default: True)
"""
self.unit = unit
self.fes = fes
self.centers = np.loadtxt(self.vgcentersfile, ndmin=1)
if not x0: # if x0 not set manually, use standard ranges attributed to the CVs
x0 = [cv.x for cv in self.activecvs]
if not ngauss: # if ngauss not set manually, use all gaussians present in vg_centers.dat
ngauss = len(self.centers)
if self.welltemp:
h = list(np.loadtxt(self.vgheightsfile)*c.to_energyunit(self.unit))
else:
h = self.get_cfg("height")*c.from_energyunit(self.get_cfg("height_unit", default="eV"))
h *= c.to_energyunit(self.unit)
# w = [cv.convert_units(cv.width) for cv in self.activecvs]
w = [cv.width for cv in self.activecvs]
if self.nactivecvs == 1:
x, w, cv = x0[0], w[0], self.activecvs[0]
vg = self.get_1dvg(x, h, w, ngauss)
if plot:
self.plot_1dvg(vg, x, cv, ngauss)
if animate:
self.animate_1dvg(vg, x, h, w, cv, ngauss, stride)
else:
x = np.meshgrid(*x0)
vg = self.get_ndvg(x, h, w, ngauss)
if self.nactivecvs == 2:
if plot:
self.plot_2dvg(vg, x, self.activecvs, ngauss)
if animate:
self.animate_2dvg(vg, x, h, w, self.activecvs, ngauss, stride)
else:
print("The metadynamics potential can only be plotted for up to two CVs.")
# fmt = "%i: %s"
# for cv in self.cvs:
# print fmt % (cv.idx, cv.name)
# plotted = raw_input("Please choose two CVs to be plotted (space-separated):\n").split()
# plotted = [int(r) for r in plotted]
# not_plotted = [i for i in range(self.nactivecvs) if not i in plotted]
# if len(not_plotted) == 1:
# not_plotted = not_plotted[0]
# fmt = "%i: %6.2f"
# for i, x_i in enumerate(x0[not_plotted]):
# print fmt % (i, x_i)
# not_plotted_value = int(raw_input("Please choose value for the CV %s:\n" % \
# self.cvs[not_plotted].name))
# # choose right part of vg for plotting from vg, x and not_plotted_value
# if not_plotted == 0:
# x = [x_i[:, not_plotted_value] for i, x_i in enumerate(x) if not i == 0]
# vg = vg[:, not_plotted_value]
# vg = vg.T
# x = [x_i.T for x_i in x]
# elif not_plotted == 1:
# x = [x_i[not_plotted_value] for i, x_i in enumerate(x) if not i == 1]
# vg = vg[not_plotted_value]
# vg = vg.T
# x = [x_i.T for x_i in x]
# elif not_plotted == 2:
# x = [x_i[:, :, not_plotted_value] for i, x_i in enumerate(x) if not i == 2]
# vg = vg[:, :, not_plotted_value]
# if plot:
# cvs = [cv for i, cv in enumerate(self.activecvs) if i in plotted]
# self.plot_2dvg(vg, x, cvs, ngauss)
# if animate:
## self.animate_2dvg(vg, x, h, w, self.activecvs, ngauss, stride)
# print "animation is currently not supported for more than 2 CVs."
# else:
# print "Plotting the metadynamics potential is only supported for up to 2 CVs"
np.save("vg_grid.npy", x)
print("Grid written to vg_grid.npy")
np.save("vg_values.npy", vg)
print("Metadynamics potential written to vg_values.npy")
[docs] def get_1dvg(self, x, heights, w, ngauss):
"""
Reconstruct vg on the 1D-grid x
Parameters
----------
x : numpy array
x-values for the collective variable
h : number or list
height of the added gaussians (single number for original, list for well-tempered)
w : number
width of the added gaussians
ngauss : number
first ngauss gaussian functions of the simulation will be evaluated
Returns
-------
vg : numpy array
metadynamics potential on the grid x
"""
if not type(heights) == type([]):
heights = [heights]*ngauss
cv = self.activecvs[0]
vg = np.zeros(x.size)
for h, g in zip(heights, self.centers[:ngauss]):
vg += h*self.gauss(x, g, w)
if cv.periodic and cv.periodictype == "continue":
shift = cv.periodicend - cv.periodicstart
vg += h*self.gauss(x + shift, g, w)
vg += h*self.gauss(x - shift, g, w)
if cv.periodic:
if cv.periodictype == "mirror":
vg += h*self.gauss(2*cv.periodicstart-x, g, w)
vg += h*self.gauss(2*cv.periodicend-x, g, w)
return vg
[docs] def get_ndvg(self, x, heights, w, ngauss):
"""
Reconstruct vg on the nD-grid x
Parameters
----------
x : list of numpy arrays
x-values for each collective variable
h : number or list
height of the added gaussians (single number for original, list for well-tempered)
w : list of numbers
widths of the added gaussians (one per CV)
ngauss : number
first ngauss gaussian functions of the simulation will be evaluated
Returns
-------
vg : numpy array
metadynamics potential on the grid x
"""
if not type(heights) == type([]):
heights = [heights]*ngauss
vg = np.zeros(x[0].shape)
for h, g in zip(heights, self.centers[:ngauss]): # add one gaussian in each step
glist = np.array([self.gauss(x[i], g[i], w[i]) for i in range(self.nactivecvs)])
vg += h*np.prod(glist, axis=0) # combine single gaussians
for j, cv in enumerate(self.activecvs):
if cv.periodic and cv.periodictype == "continue":
shifts = [cv.periodicstart-cv.periodicend, cv.periodicend-cv.periodicstart]
for shift in shifts:
xs = [x_i.copy() for x_i in x]
xs[j] += shift
glist = np.array([self.gauss(xs[i], g[i], w[i]) for i in
range(self.nactivecvs)])
vg += h*np.prod(glist, axis=0) # combine single gaussians
return vg
[docs] def plot_1dvg(self, vg, x, cv, ngauss):
"""
Plot 1D metadynamics potential and save to file vg_ngauss-g.png
Parameters
----------
vg : 1d numpy array
metadynamics potential on the grid x
x : 1d numpy array
grid for the metadynamics potential
cv : cv object
parameters of the collective variable read from input file
ngauss: number
number of gaussians contained in vg
"""
if cv.ticks:
plt.xticks(np.arange(x.min(), x.max()+cv.ticks, cv.ticks))
plt.xlabel("CV: "+cv.name)
if self.fes:
plt.ylabel("free energy / " + self.unit)
plt.title("Free Energy Surface ("+str(ngauss)+" gaussians)")
if self.welltemp:
plt.plot(x, -(self.t + self.deltat)/self.deltat*vg, color="blue", lw=2.5)
plt.fill_between(x, -(self.t + self.deltat)/self.deltat*vg, color="blue", alpha=0.1)
else:
plt.plot(x, -vg, color="blue", lw=2.5)
plt.fill_between(x, -vg, color="blue", alpha=0.1)
else:
plt.ylabel("metadynamics potential / " + self.unit)
plt.title("Metadynamics Potential ("+str(ngauss)+" gaussians)")
plt.plot(x, vg, color="blue", lw=2.5)
plt.fill_between(x, vg, color="blue", alpha=0.1)
plt.tight_layout(pad=0.1)
plt.savefig("vg_"+str(ngauss)+"-g.png", dpi=600)
plt.show()
[docs] def plot_2dvg(self, vg, x, cvs, ngauss):
"""
Plot 2D metadynamics potential and save to file vg_ngauss-g.png
Parameters
----------
vg : 2d numpy array
metadynamics potential on the grid x
x : 2d numpy meshgrid
grid for the metadynamics potential
cvs : list of cv objects
list of parameters of the collective variables read from input file
ngauss : number
number of gaussians contained in vg
"""
xmin, xmax, ymin, ymax = x[0].min(), x[0].max(), x[1].min(), x[1].max()
asp = (xmax-xmin)/(ymax-ymin)
if self.fes:
if self.welltemp:
vg = (self.t + self.deltat)/self.deltat*vg
vg -= np.min(vg)
plt.imshow(-vg, origin='lower', aspect=asp,
extent=(xmin, xmax, ymin, ymax), cmap="jet")
else:
vg -= np.min(vg)
plt.imshow(-vg, origin='lower', aspect=asp, extent=(xmin, xmax, ymin, ymax),
cmap="jet")
plt.title("Free Energy Surface ("+str(ngauss)+" gaussians)")
else:
vg -= np.min(vg)
plt.imshow(vg, origin='lower', aspect=asp, extent=(xmin, xmax, ymin, ymax), cmap="jet")
plt.title("Metadynamics Potential ("+str(ngauss)+" gaussians)")
plt.colorbar()
type1, type2 = [cv.type for cv in cvs]
if cvs[0].ticks:
plt.xticks(np.arange(xmin, xmax+cvs[0].ticks, cvs[0].ticks)) # CV 1
if cvs[1].ticks:
plt.yticks(np.arange(ymin, ymax+cvs[1].ticks, cvs[1].ticks)) # CV 2
plt.xlabel("CV 1: "+cvs[0].name)
plt.ylabel("CV 2: "+cvs[1].name)
plt.savefig("vg_"+str(ngauss)+"-g.png", dpi=600)
plt.show()
[docs] def animate_1dvg(self, vg, x, h, w, cv, ngauss, stride=1):
"""
Animate 1D metadynamics potential and save frames to directory animation
Parameters
----------
vg : 1d numpy array
metadynamics potential on the grid x
x : 1d numpy array
grid for the metadynamics potential
h : number or list
height of the added gaussians (single number for original, list for well-tempered)
w : number
width of the added gaussians
cv : cv object
parameters of the collective variable read from input file
ngauss: number
number of gaussians contained in vg
stride : number
number of gaussians to be added in each frame of the animation (default: 1)
"""
# animation function
def animate(ng):
vg_new = self.get_1dvg(x, h, w, ng)
line.set_ydata(vg_new)
for coll in ax.collections:
ax.collections.remove(coll)
ax.fill_between(x, vg_new, color="blue", alpha=0.1)
plt.title("Metadynamics Potential (%i gaussians)" % (ng))
plt.savefig("animation/%04i.png" % (ng))
return line,
# configure plot
fig, ax = plt.subplots()
if cv.ticks:
plt.xticks(np.arange(x.min(), x.max()+cv.ticks, cv.ticks))
plt.xlabel("CV: "+cv.name)
plt.ylabel("metadynamics potential / " + self.unit)
# plot and animate
line, = ax.plot(x, vg, color="blue", lw=2.5)
ax.fill_between(x, vg, color="blue", alpha=0.1)
if not os.path.isdir("animation"):
os.mkdir("animation")
ani = animation.FuncAnimation(fig, animate, np.arange(0, ngauss+1, stride), interval=25,
repeat=False)
plt.show()
[docs] def animate_2dvg(self, vg, x, h, w, cvs, ngauss, stride=1):
"""
Animate 2D metadynamics potential and save frames to directory animation
Parameters
----------
vg : 2d numpy array
metadynamics potential on the grid x
x : 2d numpy array
grid for the metadynamics potential
h : number or list
height of the added gaussians (single number for original, list for well-tempered)
w : list of numbers
widths of the added gaussians (1 per CV)
cvs : list of cv objects
parameters of the collective variables read from input file
ngauss: number
number of gaussians contained in vg
stride : number
number of gaussians to be added in each frame of the animation (default: 1)
"""
# animation function
def animate(ng):
im.set_array(self.get_ndvg(x, h, w, ng))
plt.title("Metadynamics Potential (%i gaussians)" % (ng))
plt.savefig("animation/%04i.png" % (ng))
return im,
# configure plot
fig, ax = plt.subplots()
xmin, xmax, ymin, ymax = x[0].min(), x[0].max(), x[1].min(), x[1].max()
asp = (xmax-xmin)/(ymax-ymin)
type1, type2 = [cv.type for cv in cvs]
if cvs[0].ticks:
plt.xticks(np.arange(xmin, xmax+cvs[0].ticks, cvs[0].ticks)) # CV 1
if cvs[1].ticks:
plt.yticks(np.arange(ymin, ymax+cvs[1].ticks, cvs[1].ticks)) # CV 2
plt.xlabel("CV 1: "+cvs[0].name)
plt.ylabel("CV 2: "+cvs[1].name)
# plot and animate
im = plt.imshow(vg, origin='lower', aspect=asp, extent=(xmin, xmax, ymin, ymax))
plt.colorbar()
if not os.path.isdir("animation"):
os.mkdir("animation")
ani = animation.FuncAnimation(fig, animate, np.arange(0, ngauss+1, stride), interval=25,
repeat=False)
plt.show()