Commit 8a41f6df authored by Charlles Abreu's avatar Charlles Abreu
Browse files

Better data encapsulation

parent 0b5d58d7
...@@ -123,25 +123,29 @@ class Metadynamics(object): ...@@ -123,25 +123,29 @@ class Metadynamics(object):
self.saveFrequency = saveFrequency self.saveFrequency = saveFrequency
self._id = np.random.randint(0x7FFFFFFF) self._id = np.random.randint(0x7FFFFFFF)
self._saveIndex = 0 self._saveIndex = 0
self._expansion_variables = []
for v in variables: for v in variables:
v._expanded = v.periodic and len(variables) > 1 expanded = v.periodic and len(variables) > 1
v._extraWidth = min(gridExpansion, v.gridWidth - 1) if v._expanded else 0 extraWidth = min(gridExpansion, v.gridWidth - 1) if expanded else 0
extraRange = v._extraWidth*(v.maxValue - v.minValue)/(v.gridWidth - 1) extraRange = extraWidth*(v.maxValue - v.minValue)/(v.gridWidth - 1)
v._actualWidth = v.gridWidth + 2*v._extraWidth actualWidth = v.gridWidth + 2*extraWidth
v._actualMin = v.minValue - extraRange actualMin = v.minValue - extraRange
v._actualMax = v.maxValue + extraRange actualMax = v.maxValue + extraRange
v._slice = slice(v._extraWidth, v.gridWidth + v._extraWidth) arraySlice = slice(extraWidth, v.gridWidth + extraWidth)
self._selfBias = np.zeros(tuple(v._actualWidth for v in reversed(variables))) self._expansion_variables.append(
self._totalBias = np.zeros(tuple(v._actualWidth for v in reversed(variables))) _ExpansionData(expanded, extraWidth, actualWidth, actualMin, actualMax, arraySlice),
)
self._selfBias = np.zeros(tuple(v.actualWidth for v in reversed(self._expansion_variables)))
self._totalBias = np.zeros(tuple(v.actualWidth for v in reversed(self._expansion_variables)))
self._loadedBiases = {} self._loadedBiases = {}
self._deltaT = temperature*(biasFactor-1) self._deltaT = temperature*(biasFactor-1)
varNames = ['cv%d' % i for i in range(len(variables))] varNames = ['cv%d' % i for i in range(len(variables))]
self._force = mm.CustomCVForce('table(%s)' % ', '.join(varNames)) self._force = mm.CustomCVForce('table(%s)' % ', '.join(varNames))
for name, var in zip(varNames, variables): for name, var in zip(varNames, variables):
self._force.addCollectiveVariable(name, var.force) self._force.addCollectiveVariable(name, var.force)
widths = [v._actualWidth for v in variables] widths = [ev.actualWidth for ev in self._expansion_variables]
mins = [v._actualMin for v in variables] mins = [ev.actualMin for ev in self._expansion_variables]
maxs = [v._actualMax for v in variables] maxs = [ev.actualMax for ev in self._expansion_variables]
if len(variables) == 1: if len(variables) == 1:
self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0], variables[0].periodic) self._table = mm.Continuous1DFunction(self._totalBias.flatten(), mins[0], maxs[0], variables[0].periodic)
elif len(variables) == 2: elif len(variables) == 2:
...@@ -194,7 +198,7 @@ class Metadynamics(object): ...@@ -194,7 +198,7 @@ class Metadynamics(object):
if len(self.variables) == 1: if len(self.variables) == 1:
return f return f
else: else:
s = [v._slice for v in self.variables] s = [ev.slice for ev in self._expansion_variables]
return f[s[1], s[0]] if len(self.variables) == 2 else f[s[2], s[1], s[0]] return f[s[1], s[0]] if len(self.variables) == 2 else f[s[2], s[1], s[0]]
def getCollectiveVariables(self, simulation): def getCollectiveVariables(self, simulation):
...@@ -206,7 +210,7 @@ class Metadynamics(object): ...@@ -206,7 +210,7 @@ class Metadynamics(object):
# Compute a Gaussian along each axis. # Compute a Gaussian along each axis.
axisGaussians = [] axisGaussians = []
for i,v in enumerate(self.variables): for i, (v, ev) in enumerate(zip(self.variables, self._expansion_variables)):
x = (position[i]-v.minValue) / (v.maxValue-v.minValue) x = (position[i]-v.minValue) / (v.maxValue-v.minValue)
if v.periodic: if v.periodic:
x = x % 1.0 x = x % 1.0
...@@ -214,8 +218,8 @@ class Metadynamics(object): ...@@ -214,8 +218,8 @@ class Metadynamics(object):
if v.periodic: if v.periodic:
dist = np.min(np.array([dist, np.abs(dist-1)]), axis=0) dist = np.min(np.array([dist, np.abs(dist-1)]), axis=0)
values = np.exp(-0.5*dist*dist/v._scaledVariance) values = np.exp(-0.5*dist*dist/v._scaledVariance)
if v._expanded: if ev.expanded:
n = v._extraWidth + 1 n = ev.extraWidth + 1
values = np.hstack((values[-n:-1], values, values[1:n])) values = np.hstack((values[-n:-1], values, values[1:n]))
axisGaussians.append(values) axisGaussians.append(values)
...@@ -231,9 +235,9 @@ class Metadynamics(object): ...@@ -231,9 +235,9 @@ class Metadynamics(object):
height = height.value_in_unit(unit.kilojoules_per_mole) height = height.value_in_unit(unit.kilojoules_per_mole)
self._selfBias += height*gaussian self._selfBias += height*gaussian
self._totalBias += height*gaussian self._totalBias += height*gaussian
widths = [v._actualWidth for v in self.variables] widths = [ev.actualWidth for ev in self._expansion_variables]
mins = [v._actualMin for v in self.variables] mins = [ev.actualMin for ev in self._expansion_variables]
maxs = [v._actualMax for v in self.variables] maxs = [ev.actualMax for ev in self._expansion_variables]
if len(self.variables) == 1: if len(self.variables) == 1:
self._totalBias[-1] = self._totalBias[0] self._totalBias[-1] = self._totalBias[0]
self._table.setFunctionParameters(self._totalBias.flatten(), mins[0], maxs[0]) self._table.setFunctionParameters(self._totalBias.flatten(), mins[0], maxs[0])
...@@ -329,3 +333,7 @@ class BiasVariable(object): ...@@ -329,3 +333,7 @@ class BiasVariable(object):
return quantity return quantity
_LoadedBias = namedtuple('LoadedBias', ['id', 'index', 'bias']) _LoadedBias = namedtuple('LoadedBias', ['id', 'index', 'bias'])
_ExpansionData = namedtuple(
'ExpansionData',
['expanded', 'extraWidth', 'actualWidth', 'actualMin', 'actualMax', 'slice'],
)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment