Skip to content

Added Minimum and Maximum CV + bug fixes #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mbocus wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Added Minimum and Maximum CV + bug fixes #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
10e87ed
Updated deprecated numpy functions in conprob.py
mbocus Jun 24, 2026
746d795
Create .gitignore
mbocus Jun 24, 2026
50eb812
Introduced a Minimum and Maximum collective variable
mbocus Jun 24, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
build/
ThermoLIB.egg-info/
thermolib/ext.c
6 changes: 3 additions & 3 deletions docs/_modules/thermolib/thermodynamics/condprob.html
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1213,8 +1213,8 @@ <h1>Source code for thermolib.thermodynamics.condprob</h1><div class="highlight"
<span class="c1"># Construct 1D FEP</span>
<span class="k">def</span><span class="w"> </span><span class="nf">transform</span><span class="p">(</span><span class="n">fs</span><span class="p">,</span> <span class="n">pconds</span><span class="p">):</span>
<span class="n">mask</span> <span class="o">=</span> <span class="o">~</span><span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">fs</span><span class="p">)</span>
<span class="n">ps</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapz</span><span class="p">(</span><span class="n">pconds</span><span class="p">[:,</span><span class="n">mask</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span><span class="o">*</span><span class="n">fs</span><span class="p">[</span><span class="n">mask</span><span class="p">]),</span> <span class="n">x</span><span class="o">=</span><span class="n">cvs</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
<span class="n">ps</span> <span class="o">/=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapz</span><span class="p">(</span><span class="n">ps</span><span class="p">,</span> <span class="n">x</span><span class="o">=</span><span class="n">qs</span><span class="p">)</span>
<span class="n">ps</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapezoid</span><span class="p">(</span><span class="n">pconds</span><span class="p">[:,</span><span class="n">mask</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span><span class="o">*</span><span class="n">fs</span><span class="p">[</span><span class="n">mask</span><span class="p">]),</span> <span class="n">x</span><span class="o">=</span><span class="n">cvs</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
<span class="n">ps</span> <span class="o">/=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapezoid</span><span class="p">(</span><span class="n">ps</span><span class="p">,</span> <span class="n">x</span><span class="o">=</span><span class="n">qs</span><span class="p">)</span>
<span class="n">fs_new</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">([</span><span class="nb">len</span><span class="p">(</span><span class="n">qs</span><span class="p">)],</span> <span class="nb">float</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">nan</span>
<span class="n">fs_new</span><span class="p">[</span><span class="n">ps</span><span class="o">&gt;</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">ps</span><span class="p">[</span><span class="n">ps</span><span class="o">&gt;</span><span class="mi">0</span><span class="p">])</span><span class="o">/</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span>
<span class="k">return</span> <span class="n">fs_new</span>
Expand Down Expand Up @@ -1435,7 +1435,7 @@ <h1>Source code for thermolib.thermodynamics.condprob</h1><div class="highlight"
<span class="c1">#construct 2D FES</span>
<span class="k">def</span><span class="w"> </span><span class="nf">transform</span><span class="p">(</span><span class="n">fs</span><span class="p">,</span> <span class="n">pconds</span><span class="p">):</span>
<span class="n">mask</span> <span class="o">=</span> <span class="o">~</span><span class="n">np</span><span class="o">.</span><span class="n">isnan</span><span class="p">(</span><span class="n">fs</span><span class="p">)</span>
<span class="n">ps</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapz</span><span class="p">(</span><span class="n">pconds</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">mask</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span><span class="o">*</span><span class="n">fs</span><span class="p">[</span><span class="n">mask</span><span class="p">]),</span> <span class="n">x</span><span class="o">=</span><span class="n">fep</span><span class="o">.</span><span class="n">cvs</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
<span class="n">ps</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">trapezoid</span><span class="p">(</span><span class="n">pconds</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">mask</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span><span class="o">*</span><span class="n">fs</span><span class="p">[</span><span class="n">mask</span><span class="p">]),</span> <span class="n">x</span><span class="o">=</span><span class="n">fep</span><span class="o">.</span><span class="n">cvs</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span>
<span class="n">ps</span> <span class="o">/=</span> <span class="n">integrate2d</span><span class="p">(</span><span class="n">ps</span><span class="p">,</span> <span class="n">x</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">qs</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">y</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">qs</span><span class="p">[</span><span class="mi">1</span><span class="p">])</span>
<span class="n">fs_new</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">ps</span><span class="o">.</span><span class="n">shape</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">nan</span>
<span class="n">fs_new</span><span class="p">[</span><span class="n">ps</span><span class="o">&gt;</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">ps</span><span class="p">[</span><span class="n">ps</span><span class="o">&gt;</span><span class="mi">0</span><span class="p">])</span><span class="o">/</span><span class="n">fep</span><span class="o">.</span><span class="n">beta</span>
Expand Down
6 changes: 3 additions & 3 deletions thermolib/thermodynamics/condprob.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -784,8 +784,8 @@ def transform(self, fep, f_output_unit=None, f_label=None, f_output_class=BaseFr
# Construct 1D FEP
def transform(fs, pconds):
mask = ~np.isnan(fs)
ps = np.trapz(pconds[:,mask]*np.exp(-fep.beta*fs[mask]), x=cvs[mask])
ps /= np.trapz(ps, x=qs)
ps = np.trapezoid(pconds[:,mask]*np.exp(-fep.beta*fs[mask]), x=cvs[mask])
ps /= np.trapezoid(ps, x=qs)
fs_new = np.zeros([len(qs)], float)*np.nan
fs_new[ps>0] = -np.log(ps[ps>0])/fep.beta
return fs_new
Expand Down Expand Up @@ -991,7 +991,7 @@ def deproject(self, fep, f_output_unit=None, f_label=None, f_output_class=FreeEn
#construct 2D FES
def transform(fs, pconds):
mask = ~np.isnan(fs)
ps = np.trapz(pconds[...,mask]*np.exp(-fep.beta*fs[mask]), x=fep.cvs[mask])
ps = np.trapezoid(pconds[...,mask]*np.exp(-fep.beta*fs[mask]), x=fep.cvs[mask])
ps /= integrate2d(ps, x=self.qs[0], y=self.qs[1])
fs_new = np.zeros(ps.shape, float)*np.nan
fs_new[ps>0] = -np.log(ps[ps>0])/fep.beta
Expand Down
107 changes: 84 additions & 23 deletions thermolib/thermodynamics/cv.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@

__all__ = [
'CenterOfMass', 'CenterOfPosition', 'NormalizedAxis', 'NormalToPlane', 'Distance', 'DistanceCOP', 'CoordinationNumber', 'OrthogonalDistanceToPore',
'Average', 'Difference', 'Minimum', 'LinearCombination', 'DotProduct', 'DistOrthProjOrig'
'Average', 'Difference', 'Minimum', 'Maximum', 'LinearCombination', 'DotProduct', 'DistOrthProjOrig'
]

class CollectiveVariable(object):
Expand Down Expand Up @@ -759,34 +759,38 @@ def compute(self, atoms, deriv=True):

class Minimum(CollectiveVariable):
'''
Class to implement a collective variable representing the minimum of two other collective variables:
Class to implement a collective variable representing the minimum of a set of collective variables. To make its definition continuous, it is implemented as:

.. math:: CV &= \\min\\left(CV_1,CV_2\\right)
.. math:: \\min(\\mathbf{s}) = -\\frac{1}{\\beta} \\log \\left( \\sum_i \\exp(-\\beta s_i) \\right)

This definition is consistent with the ALT_MIN function of the MULTICOLVAR Plumed module.
'''

type = 'scalar'

def __init__(self, cv1, cv2, name=None):
def __init__(self, cvs, beta=50.0, name=None):
'''
:param cv1: first collective variable in the minimum
:type cv1: any child class of :py:class:`CollectiveVariable <thermolib.thermodynamics.cv.CollectiveVariable>`
:param cv_list: list of instances of child classes of :py:class:`CollectiveVariable <thermolib.thermodynamics.cv.CollectiveVariable>`
:type cv1: List

:param cv2: second collective variable in the minimum
:type cv2: any child class of :py:class:`CollectiveVariable <thermolib.thermodynamics.cv.CollectiveVariable>`
:param beta: parameter to regulate the softness of the minimum.
:type beta: float | optional, default=50.0

:param name: Name of CV for printing/logging purposes. If None, the default implemented in the ``_default_name`` routine will be used.
:type name: str | None, optional, default=None
'''
self.cv1 = cv1
self.cv2 = cv2
self.cvs = cvs
for cv in cvs:
assert cv.type == 'scalar', "The Minimum function is only implemented for scalar collective variables."
self.beta = beta
CollectiveVariable.__init__(self, name=name)

def _default_name(self):
return 'Min(%s,%s)' %(self.cv1.name, self.cv2.name)
return 'Min(' + ",".join([cv.name for cv in self.cvs]) + ")"

def compute(self, atoms, deriv=True):
'''
Compute the minimum (and optionally gradient) of the two CVs for the given atomic coordinates
Compute the minimum (and optionally gradient) of the CV list for the given atomic coordinates

:param atoms: ASE Atoms object on which the CV needs to be computed
:type atoms: ase.Atoms
Expand All @@ -797,19 +801,76 @@ def compute(self, atoms, deriv=True):
:return: CV value and potentially the gradient
:rtype: np.ndarray(3) or float,np.ndarray([3,Natoms,3])
'''

if not deriv:
cv1 = self.cv1.compute(atoms, deriv=False)
cv2 = self.cv2.compute(atoms, deriv=False)
return min(cv1,cv2)
values = [cv.compute(atoms, deriv=False) for cv in self.cvs]
S = sum(np.exp(-self.beta * v) for v in values)
return - (1.0 / self.beta) * np.log(S)
else:
cv1, grad1 = self.cv1.compute(atoms, deriv=True)
cv2, grad2 = self.cv2.compute(atoms, deriv=True)
value = min(cv1,cv2)
if cv1<cv2: # note that right not gradient is discontinuous, can be fixed using plumed MIN CV
grad = grad1
else:
grad = grad2
return value, grad
vs, gs = zip(*(cv.compute(atoms, deriv=True) for cv in self.cvs))
weights = [np.exp(-self.beta * v) for v in vs]
S = sum(weights)
grad = sum(w * g for w, g in zip(weights, gs)) / S
value = - (1.0 / self.beta) * np.log(S)
return value, grad


class Maximum(CollectiveVariable):
'''
Class to implement a collective variable representing the maximum of a set of collective variables. To make its definition continuous, it is implemented as:

.. math:: \\max(\\mathbf{s}) = \\beta \\log \\left( \\sum_i \\exp(s_i / \\beta) \\right)

This definition is consistent with the MAX function of the MULTICOLVAR Plumed module.
'''

type = 'scalar'

def __init__(self, cvs, beta=50.0, name=None):
'''
:param cv_list: list of instances of child classes of :py:class:`CollectiveVariable <thermolib.thermodynamics.cv.CollectiveVariable>`
:type cv1: List

:param beta: parameter to regulate the softness of the maximum.
:type beta: float | optional, default=50.0

:param name: Name of CV for printing/logging purposes. If None, the default implemented in the ``_default_name`` routine will be used.
:type name: str | None, optional, default=None
'''
self.cvs = cvs
for cv in cvs:
assert cv.type == 'scalar', "The Maximum function is only implemented for scalar collective variables."
self.beta = beta
CollectiveVariable.__init__(self, name=name)

def _default_name(self):
return 'Max(' + ",".join([cv.name for cv in self.cvs]) + ")"

def compute(self, atoms, deriv=True):
'''
Compute the maximum (and optionally gradient) of the CV list for the given atomic coordinates

:param atoms: ASE Atoms object on which the CV needs to be computed
:type atoms: ase.Atoms

:param deriv: if True, also compute and return the gradient of the CV towards all atomic coordinates
:type deriv: bool, optional, default=True

:return: CV value and potentially the gradient
:rtype: np.ndarray(3) or float,np.ndarray([3,Natoms,3])
'''

if not deriv:
values = [cv.compute(atoms, deriv=False) for cv in self.cvs]
S = sum(np.exp(v / self.beta) for v in values)
return self.beta * np.log(S)
else:
vs, gs = zip(*(cv.compute(atoms, deriv=True) for cv in self.cvs))
weights = [np.exp(v / self.beta) for v in vs]
S = sum(weights)
grad = sum(w * g for w, g in zip(weights, gs)) / S
value = self.beta * np.log(S)
return value, grad


class LinearCombination(CollectiveVariable):
Expand Down