diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py index 3c55d24af40466c45101f2ea62132b319f49cae6..1c436689890b0e887a3ad32cc4db6caaea05ce90 100644 --- a/pycgtool/bondset.py +++ b/pycgtool/bondset.py @@ -63,22 +63,21 @@ class Bond: if not self.values: raise ValueError("No bonds were measured between atoms {0}".format(self.atoms)) - mean, var = stat_moments(self.values) - + values = np.array(self.values) if len(self.atoms) > 2: - mean_tmp = math.radians(mean) - var *= (math.pi * math.pi) / (180. * 180.) - else: - mean_tmp = mean + values = np.radians(values) - self.eqm = mean with np.errstate(divide="raise"): + self.eqm = self._func_form.eqm(values, temp) try: - self.fconst = self._func_form(mean_tmp, var, temp) + self.fconst = self._func_form.fconst(values, temp) except FloatingPointError: # Happens when variance is 0, i.e. we only have one value self.fconst = float("inf") + if len(self.atoms) > 2: + self.eqm = math.degrees(self.eqm) + def __repr__(self): try: return "<Bond containing atoms {0} with r_0 {1:.3f} and force constant {2:.3e}>".format( diff --git a/pycgtool/functionalforms.py b/pycgtool/functionalforms.py index 919e425aad7ae3c01889cf71162d2fc0ba209f3d..987d96d87bc0105a7c64c242c070734226999f9b 100644 --- a/pycgtool/functionalforms.py +++ b/pycgtool/functionalforms.py @@ -1,6 +1,8 @@ import abc import math +import numpy as np + from pycgtool.util import SimpleEnum @@ -53,15 +55,26 @@ class FunctionalForm(object, metaclass=abc.ABCMeta): New functional forms must define a static __call__ method. """ + @staticmethod + def eqm(values, temp): + """ + Calculate equilibrium value. + May be overridden by functional forms. + + :param values: Measured internal coordinate values from which to calculate equilibrium value + :param temp: Temperature of simulation + :return: Calculated equilibrium value + """ + return np.nanmean(values) + @staticmethod @abc.abstractstaticmethod - def __call__(mean, var, temp): + def fconst(values, temp): """ Calculate force constant. Abstract static method to be defined by all functional forms. - :param mean: Mean of internal coordinate distribution - :param var: Variance of internal coordinate distribution + :param values: Measured internal coordinate values from which to calculate force constant :param temp: Temperature of simulation :return: Calculated force constant """ @@ -70,31 +83,34 @@ class FunctionalForm(object, metaclass=abc.ABCMeta): class Harmonic(FunctionalForm): @staticmethod - def __call__(mean, var, temp): + def fconst(values, temp): rt = 8.314 * temp / 1000. + var = np.nanvar(values) return rt / var class CosHarmonic(FunctionalForm): @staticmethod - def __call__(mean, var, temp): + def fconst(values, temp): rt = 8.314 * temp / 1000. + mean = CosHarmonic.eqm(values, temp) + var = np.nanvar(values) return rt / (math.sin(mean)**2 * var) class MartiniDefaultLength(FunctionalForm): @staticmethod - def __call__(mean, var, temp): + def fconst(values, temp): return 1250. class MartiniDefaultAngle(FunctionalForm): @staticmethod - def __call__(mean, var, temp): + def fconst(values, temp): return 25. class MartiniDefaultDihedral(FunctionalForm): @staticmethod - def __call__(mean, var, temp): + def fconst(values, temp): return 50. diff --git a/test/test_functionalforms.py b/test/test_functionalforms.py index 56056e53ab5076b659a21c3b4d076df387f66d49..78a4464d706c19c3fea4210207b21e1ddc6b0a3c 100644 --- a/test/test_functionalforms.py +++ b/test/test_functionalforms.py @@ -14,12 +14,17 @@ class FunctionalFormTest(unittest.TestCase): def test_functional_form_new(self): class TestFunc(FunctionalForm): @staticmethod - def __call__(mean, var, temp): - return "TestResult" + def eqm(values, temp): + return "TestResultEqm" + + @staticmethod + def fconst(values, temp): + return "TestResultFconst" funcs = FunctionalForms() self.assertIn("TestFunc", funcs) - self.assertEqual("TestResult", funcs.TestFunc(None, None, None)) + self.assertEqual("TestResultEqm", funcs.TestFunc.eqm(None, None)) + self.assertEqual("TestResultFconst", funcs.TestFunc.fconst(None, None)) if __name__ == '__main__': unittest.main()