From 7cddfda4afa2b096d592e28632c56cfbb8ad11d1 Mon Sep 17 00:00:00 2001
From: James Graham <J.A.Graham@soton.ac.uk>
Date: Fri, 24 Feb 2017 15:16:07 +0000
Subject: [PATCH] Give all measured values to FunctionalForm - no longer
limited to using mean, var
---
pycgtool/bondset.py | 15 +++++++--------
pycgtool/functionalforms.py | 32 ++++++++++++++++++++++++--------
test/test_functionalforms.py | 11 ++++++++---
3 files changed, 39 insertions(+), 19 deletions(-)
diff --git a/pycgtool/bondset.py b/pycgtool/bondset.py
index 3c55d24..1c43668 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 919e425..987d96d 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 56056e5..78a4464 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()
--
GitLab