From 10e87edb25607e70f18f26b354e7ae9120eba3a0 Mon Sep 17 00:00:00 2001 From: mbocus Date: Wed, 24 Jun 2026 10:14:01 +0200 Subject: [PATCH 1/3] Updated deprecated numpy functions in conprob.py The deprecated `np.trpz` function has been replaced with its modern counterpart `np.trapezoid`, solving https://github.com/molmod/ThermoLIB/issues/4 --- docs/_modules/thermolib/thermodynamics/condprob.html | 6 +++--- thermolib/thermodynamics/condprob.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/_modules/thermolib/thermodynamics/condprob.html b/docs/_modules/thermolib/thermodynamics/condprob.html index 7047c02..225684c 100644 --- a/docs/_modules/thermolib/thermodynamics/condprob.html +++ b/docs/_modules/thermolib/thermodynamics/condprob.html @@ -1213,8 +1213,8 @@

Source code for thermolib.thermodynamics.condprob

# 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 @@ -1435,7 +1435,7 @@

Source code for thermolib.thermodynamics.condprob

#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 diff --git a/thermolib/thermodynamics/condprob.py b/thermolib/thermodynamics/condprob.py index e50a493..4412c6f 100644 --- a/thermolib/thermodynamics/condprob.py +++ b/thermolib/thermodynamics/condprob.py @@ -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 @@ -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 From 746d7955e5e7102a9c67ca973ee9e427c23a3f2e Mon Sep 17 00:00:00 2001 From: mbocus Date: Wed, 24 Jun 2026 15:59:20 +0200 Subject: [PATCH 2/3] Create .gitignore Added a `.gitignore` file to prevent uploading files generated when pip installing the repo for testing. --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6970d24 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +build/ +ThermoLIB.egg-info/ +thermolib/ext.c From 50eb8124fdc900bf6a8b8a4e7a797932e7192115 Mon Sep 17 00:00:00 2001 From: mbocus Date: Wed, 24 Jun 2026 17:09:37 +0200 Subject: [PATCH 3/3] Introduced a Minimum and Maximum collective variable The old `Minimum` implementation with discontinuous derivatives has been replaced with the continuous version used by PLUMED in the ALT_MIN function of the MULTICOLVAR module. Additionally, a `Maximum` is introduced that corresponds to the soft MAX of PLUMED. Both implementations have been tested and validated with `test_CV_implementations` on a model system (`tol=1E-6`). --- thermolib/thermodynamics/cv.py | 107 ++++++++++++++++++++++++++------- 1 file changed, 84 insertions(+), 23 deletions(-) diff --git a/thermolib/thermodynamics/cv.py b/thermolib/thermodynamics/cv.py index 586490b..43a60f9 100644 --- a/thermolib/thermodynamics/cv.py +++ b/thermolib/thermodynamics/cv.py @@ -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): @@ -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 ` + :param cv_list: list of instances of child classes of :py:class:`CollectiveVariable ` + :type cv1: List - :param cv2: second collective variable in the minimum - :type cv2: any child class of :py:class:`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 @@ -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` + :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):