Source code for metafalcon.meta
#!/usr/bin/env python
"""Contains a class for the modification of MD forces by the metadynamics scheme."""
from __future__ import print_function
import os
import numpy as np
import itertools as it
import json
from .cvs import types as cvtypes
from . import constants as c
[docs]class metadynamics:
"""
Class for performing metadynamics simulations.
Parameters
----------
symbols : list of str, optional
element symbols needed for calculation of coordination numbers (default: None)
file_extension : string, optional
string to be appended to in- and output filenames to distinguish between multiple
metadynamics objects in the same folder (default: None)
restart : bool, optional
read gaussian centers from vgcentersfile to restart previous run (default: False)
update : bool, optional
whether metadynamics potential should be updated during the simulation (default: True)
use_threshold : bool, optional
use threshold instead of fixed tau_g (default: False)
Attributes
----------
symbols : list of str
element symbols of the molecule
file_extension : string, optional
string to be appended to output filenames to distinguish between multiple metadynamics \
objects in the same folder
restart : bool
read gaussian centers from vgcentersfile to restart previous run
update : bool
whether metadynamics potential should be updated during the simulation
use_threshold : bool, optional
use threshold instead of fixed tau_g (default: False)
vg_centers : list of float
list of all gaussian centers added to vg
vg_steps : list of int
step numbers, in which gaussians have been added to vg
cfg : dict
configuration dictionary read from the input file
cvs : list
list that containing the CV objects
"""
def __init__(self, symbols=None, file_extension=None, restart=False, update=True,
use_threshold=False):
"""Construct metadynamics object."""
self.symbols = symbols
if file_extension:
self.file_extension = "_"+file_extension
else:
self.file_extension = ""
self.infile = "meta-config.json"
self.vgcentersfile = "vg_centers%s.dat" % self.file_extension
self.vgheightsfile = "vg_heights%s.dat" % self.file_extension
self.vgstepsfile = "vg_steps%s.dat" % self.file_extension
self.restart = restart
self.update = update
self.use_threshold = use_threshold
self.cvs = []
self.welltemp = False
self.vg_centers = []
self.vg_steps = [0]
[docs] def read_input(self):
"""Parse input from the input file"""
with open(self.infile) as f:
self.cfg = json.load(f)["metadynamics" + self.file_extension]
for idx, cv in enumerate(self.get_cfg("collective variables")):
self.cvs.append(cvtypes.types[cv["type"]](idx, self.symbols, **cv))
if "plot" in cv.keys():
self.cvs[-1].set_plot(cv["plot"])
# if "wall" in cv.keys():
# self.cvs[-1].set_wall(cv["wall"])
self.ncvs = len(self.cvs)
self.activecvs = [cv for cv in self.cvs if cv.active == True]
self.nactivecvs = len(self.activecvs)
if "well-tempered" in self.cfg.keys():
self.welltemp = True
self.deltat = self.get_cfg("well-tempered")["delta_T"]
self.t = self.get_cfg("well-tempered")["T"]
self.heights = []
self.tau_g = self.get_cfg("tau_g", default=0)
if "threshold" in self.cfg.keys():
self.use_threshold = True
self.height = self.get_cfg("height")*c.from_energyunit(self.get_cfg("height_unit",
default="eV"))
[docs] def get_cfg(self, key, default=None):
try:
return self.cfg[key]
except KeyError:
if default is not None:
return default
else:
print("Error: mandatory keyword %s missing in the input file!" % key)
[docs] def set_dftb_object(self, dftb):
"""For Mulliken charges, set the `DFTB.PES.PotentialEnergySurfaces` object."""
for cv in self.cvs:
if cv.type == "mullcharge":
cv.set_dftb_object(dftb)
[docs] def initialize_vg(self):
"""Initialize vg and output files"""
self.vg = 0.0
if self.restart:
with open(self.vgcentersfile) as f:
f.next()
self.vg_centers = [[float(s) for s in line.split()] for line in f]
if self.welltemp:
with open(self.vgheightsfile) as f:
f.next()
self.heights = [[float(s) for s in line.split()] for line in f]
if os.path.exists(self.vgstepsfile):
with open(self.vgstepsfile) as f:
f.next()
self.vg_steps = [int(line) for line in f]
if not self.vg_steps:
self.vg_steps.append(0)
else:
# file in which the centers of the gaussians are stored
with open(self.vgcentersfile, "w") as f:
f.write("# gaussian centers of vg, for parameters see "+str(self.infile)+"\n")
if self.welltemp:
with open(self.vgheightsfile, "w") as f:
f.write("# gaussian heights of well-tempered metadynamics vg\n")
with open(self.vgstepsfile, "w") as f:
f.write("# step numbers of added gaussians\n")
# file in which threshold values are stored, if used
if self.use_threshold:
with open("thresh_values.dat", "w") as f:
f.write("# threshold values\n")
with open("mod.dat", "w") as f:
f.write("# modication to gap\n")
with open("dmod.dat", "w") as f:
f.write("# derivative of modification to gap")
for cv in self.cvs:
cv.initialize()
[docs] def set_s_and_ds(self):
"""Calculate values of CV s and derivatives wrt coordinates ds_dr"""
for cv in self.cvs: # delete old tmp files if they exist
cv.delete_tmpfile()
self.s = np.zeros(self.ncvs)
self.ds_dr = np.zeros((self.ncvs, len(self.coordinates), 3))
for i, cv in enumerate(self.cvs):
cv.set_s_and_ds(self.coordinates)
self.s[i] = cv.get_s()
self.ds_dr[i] = cv.get_ds_dr()
[docs] def get_height(self):
"""Return the height of the gaussian in the current step."""
return self.height*np.exp(-self.old_vg/(c.kb*self.deltat))
[docs] def get_vg(self):
"""
Return vg value in the current step
Returns
-------
vg : float
value of the metadynamics potential in the current step
"""
def gaussian_exp(cv, s0, shift=0.0):
return -(cv.s-s0+shift)**2/2/cv.width**2
shifts = [p for p in it.product(*[cv.shift for cv in self.activecvs])] # periodicity of torsion angles
vg = 0.0
for i, center in enumerate(self.vg_centers):
if self.welltemp:
height = self.heights[i]
else:
height = self.height
for shift in shifts:
exps = [gaussian_exp(cv, cv_center, shift=cv_shift) for cv, cv_center, cv_shift
in zip(self.activecvs, center, shift)]
vg += height*np.exp(np.sum(exps))
## only for 1-dimensional metadynamics with periodic cv, ugly workaround
if len(self.activecvs) == 1 and self.activecvs[0].periodic:
cv = self.activecvs[0]
if cv.periodictype == "mirror":
vg += height*np.exp(gaussian_exp(cv, 2*cv.periodicstart-center[0]))
vg += height*np.exp(gaussian_exp(cv, 2*cv.periodicend-center[0]))
return vg
[docs] def get_dvg_ds(self):
"""
Return dvg_ds in the current step
Returns
-------
dvg_ds : np.array
derivative of the metadynamics potential wrt CV in the current step (shape: (Ncvs))
"""
def gaussian_exp(cv, s0, shift=0.0):
return -(cv.s-s0+shift)**2/2/cv.width**2
def gaussian_dexp(cv, s0, shift=0.0):
return -(cv.s-s0+shift)/cv.width**2
shifts = [p for p in it.product(*[cv.shift for cv in self.activecvs])] # periodicity of torsion angles
dvg_ds = np.zeros(self.nactivecvs)
for i, center in enumerate(self.vg_centers):
if self.welltemp:
height = self.heights[i]
else:
height = self.height
for shift in shifts:
exps = [gaussian_exp(cv, cv_center, shift=cv_shift) for cv, cv_center, cv_shift
in zip(self.activecvs, center, shift)]
dexps = [gaussian_dexp(cv, cv_center, shift=cv_shift) for cv, cv_center, cv_shift
in zip(self.activecvs, center, shift)]
for cv in self.activecvs:
dvg_ds[cv.idx] += height*np.exp(np.sum(exps))*dexps[cv.idx]
## only for 1-dimensional metadynamics with periodic cv, ugly workaround
if len(self.activecvs) == 1 and self.activecvs[0].periodic:
cv = self.activecvs[0]
if cv.periodictype == "mirror":
start = 2*cv.periodicstart-center[0]
end = 2*cv.periodicend-center[0]
dvg_ds[0] += height*np.exp(gaussian_exp(cv, start))*gaussian_dexp(cv, start)
dvg_ds[0] += height*np.exp(gaussian_exp(cv, end))*gaussian_dexp(cv, end)
return dvg_ds
[docs] def set_thresh_value(self, value):
"""Set value for the threshold evaluation."""
self.thresh_value = value
# print("meta.thresh_value set to ", self.thresh_value)
[docs] def update_vg(self, step=0):
"""
Add new gaussian to vg every tau_g steps.
Parameters
----------
step : int
step number of the dynamics simulation (default: 0)
"""
if self.use_threshold:
# print "meta.thresh_value is now ", self.thresh_value
# self.thresh_value += self.vg # add metadynamics potential to current gap
# print "and after addition of vg: ", self.thresh_value
with open("thresh_values.dat", "a") as f:
f.write("%12.8f \n" % self.thresh_value)
if (step - self.vg_steps[-1] == self.tau_g):
f.write("tau_g is the reason\n")
f.write("step: %d\n" % step)
f.write("last step: %d\n" % self.vg_steps[-1])
f.write("tau_g: %d\n" % self.tau_g)
if (self.use_threshold and self.thresh_value <= self.cfg["threshold"]):
f.write("threshold is the reason\n")
# if metadynamics step, update lamdvg and write s to vgcentersfile
self.updated = False
if self.tau_g:
# print "step: ", step
# print "tau_g: ", self.tau_g
# print "last step: ", self.vg_steps[-1]
# print "step - last step: ", step - self.vg_steps[-1]
do_update_tau_g = step > 0 and not((step - self.vg_steps[-1]) % self.tau_g)
# print "do_update_tau_g: ", do_update_tau_g
else:
do_update_tau_g = False
if self.use_threshold:
do_update_thresh = (self.thresh_value <= self.cfg["threshold"])
else:
do_update_thresh = False
### CV update threshold
for cv in self.activecvs:
if not cv.update:
do_update_tau_g = False
do_update_thresh = False
print("CV %s outside epsilon range in step %d" % (cv.name, step))
### end CV update threshold
if self.update and (do_update_tau_g or do_update_thresh):
if self.use_threshold:
print("threshold is reached, update metadynamics")
print("cv-values: \n", [cv.s for cv in self.activecvs])
if self.welltemp:
self.old_vg = self.get_vg()
self.heights.append(self.get_height())
with open(self.vgheightsfile, "a") as outf:
outf.write("%16.9e\n" % self.heights[-1])
self.vg_centers.append([cv.s for cv in self.activecvs])
self.vg_steps.append(step)
fmt = self.nactivecvs*"%16.9e "+"\n"
with open(self.vgcentersfile, "a") as outf:
centers_converted = [cv.convert_units(cv.s) for cv in self.activecvs]
outf.write(fmt % tuple(centers_converted))
with open(self.vgstepsfile, "a") as outf:
outf.write("%d\n" % step)
# if self.welltemp:
# with open(self.vgheightsfile, "a") as outf:
# outf.write("%12.6f\n" % (self.heights[-1]))
self.updated = True
self.vg = self.get_vg()
# get dvg_dr
self.dvg_dr = np.zeros_like(self.coordinates)
self.dvg_ds = self.get_dvg_ds()
for dvg_ds, cv in zip(self.dvg_ds, self.activecvs):
self.dvg_dr += dvg_ds*cv.ds_dr
# add wall potential
for cv in self.cvs:
if cv.use_wall:
cv.set_wall_s()
with open("wall_%d.dat" % cv.idx, "a") as f:
f.write("%16.9e\n" % cv.wallfunc(cv.s, **cv.wall_kwargs))
dwall_ds = cv.dwallfunc(cv.s, **cv.wall_kwargs)
self.dvg_dr += dwall_ds*cv.ds_dr
if self.use_threshold: #
print("save temporary files tmp_mod and tmp_dmod")
np.save("tmp_mod", np.array(self.vg))
np.save("tmp_dmod", self.dvg_dr)
with open("mod.dat", "a") as f:
f.write("%16.8f\n" % self.vg)
with open("dmod.dat", "a") as f:
f.write("%16.8f\n" % np.linalg.norm(self.dvg_dr))
[docs] def get_forces(self, coords, step):
"""
Check if vg should be updated and return metadynamics forces.
Parameters
----------
coords : np.array
coordinates (shape: (N, 3))
step : int
step number of the dynamics simulation
Returns
-------
forces : np.array
forces obtained from the current metadynamics potential (shape: (N, 3))
"""
# update step for every cv
for cv in self.cvs:
cv.set_step(step)
self.coordinates = coords
self.set_s_and_ds()
for cv in self.cvs:
cv.set_update(cv.s)
if cv.type == "energygap" and cv.do_multistate and cv.offdiag.updated:
cv.update = False
self.update_vg(step)
return -self.dvg_dr
[docs] def get_coupled_energies_grads(self, coords, en1, en2, grad1, grad2, step=0):
"""
Couple the given energies and gradients by metadynamics potential.
Parameters
----------
coords : np.array
coordinates (shape: (N, 3))
en1 : float
lower state energy to be coupled
en2 : float
upper state energy to be coupled
grad1 : np.array
lower state energy gradient to be coupled (shape: (N, 3))
grad2 : np.array
upper state energy gradient to be coupled (shape: (N, 3))
step : int, optional
step number of the dynamics simulation needed for metadynamics wih tau_g (default: 0)
Returns
-------
new_en1 : float
lower state energy after coupling
new_en2 : float
upper state energy after coupling
new_grad1 : np.array
lower state energy gradient after coupling (shape: (N, 3))
new_grad2 : np.array
upper state energy gradient after coupling (shape: (N, 3))
"""
# update step for every cv
for cv in self.cvs:
cv.set_step(step)
# get current coupling
self.coordinates = coords
self.set_s_and_ds()
# calculate new energies
sqrt = np.sqrt((en1-en2)**2+4*self.get_vg()**2) # number
new_en1 = 0.5*(en1+en2-sqrt) # number
new_en2 = 0.5*(en1+en2+sqrt) # number
# set threshold
self.set_thresh_value(abs(new_en2-new_en1))
for cv in self.cvs:
cv.set_update(abs(new_en2-new_en1))
# update metadynamics and get coupling
self.update_vg(step)
# calculate new energies if potential has been updated
if self.updated:
sqrt = np.sqrt((en1-en2)**2+4*self.get_vg()**2) # number
new_en1 = 0.5*(en1+en2-sqrt) # number
new_en2 = 0.5*(en1+en2+sqrt) # number
# calculate new gradients
dsqrt = 1/sqrt*((en2-en1)*(grad2-grad1)+4*self.get_vg()*self.dvg_dr) # coords-shapes (N, 3)
new_grad1 = 0.5*(grad1+grad2-dsqrt) # coords-shaped (N, 3)
new_grad2 = 0.5*(grad1+grad2+dsqrt) # coords-shaped (N, 3)
# with open("gap_derivatives.dat", "a") as f:
# f.write("%4d %12.8f %12.8f %12.8f %12.8f \n" % (step,
# (en2-en1)/(new_en2-new_en1),
# np.linalg.norm(grad2-grad1),
# np.linalg.norm(4*self.dvg_dr),
# np.linalg.norm(new_grad2-new_grad1)))
# with open("dnewgap_dx.dat", "a") as f:
# f.write("dnewgap_dx\n")
# for gapgrad in (new_grad2-new_grad1):
# f.write("%12.8f %12.8f %12.8f\n" % (gapgrad[0], gapgrad[1], gapgrad[2]))
# with open("dgap_dx.dat", "a") as f:
# f.write("dgap_dx\n")
# for gapgrad in (grad2-grad1):
# f.write("%12.8f %12.8f %12.8f\n" % (gapgrad[0], gapgrad[1], gapgrad[2]))
# with open("4dj_dx.dat", "a") as f:
# f.write("4dj_dx\n")
# for jgrad in (4*self.dvg_dr):
# f.write("%12.8f %12.8f %12.8f\n" % (jgrad[0], jgrad[1], jgrad[2]))
return new_en1, new_en2, new_grad1, new_grad2